Click here to go to the corresponding page for the latest version of DIALS
Source code for dials.algorithms.symmetry
"""Methods for symmetry determination.
This module provides a base class for symmetry determination algorithms.
"""
from __future__ import division, absolute_import, print_function
import logging
logger = logging.getLogger(__name__)
from cStringIO import StringIO
import libtbx
from scitbx.array_family import flex
from cctbx import adptbx
from cctbx import sgtbx
from cctbx import uctbx
from cctbx.sgtbx.lattice_symmetry import metric_subgroups
from mmtbx import scaling
from mmtbx.scaling import absolute_scaling
from mmtbx.scaling import matthews
[docs]class symmetry_base(object):
"""Base class for symmetry analysis."""
def __init__(
self,
intensities,
normalisation="ml_aniso",
lattice_symmetry_max_delta=2.0,
d_min=libtbx.Auto,
min_i_mean_over_sigma_mean=4,
min_cc_half=0.6,
relative_length_tolerance=None,
absolute_angle_tolerance=None,
):
"""Initialise a symmetry_base object.
Args:
intensities (cctbx.miller.array): The intensities on which to perform
symmetry anaylsis.
normalisation (str): The normalisation method to use. Possible choices are
'kernel', 'quasi', 'ml_iso' and 'ml_aniso'. Set to None to switch off
normalisation altogether.
lattice_symmetry_max_delta (float): The maximum value of delta for
determining the lattice symmetry using the algorithm of Le Page (1982).
d_min (float): Optional resolution cutoff to be applied to the input
intensities. If set to :data:`libtbx.Auto` then d_min will be
automatically determined according to the parameters
``min_i_mean_over_sigma_mean`` and ``min_cc_half``.
min_i_mean_over_sigma_mean (float): minimum value of |I|/|sigma(I)| for
automatic determination of resolution cutoff.
min_cc_half (float): minimum value of CC1/2 for automatic determination of
resolution cutoff.
relative_length_tolerance (float): Relative length tolerance in checking
consistency of input unit cells against the median unit cell.
absolute_angle_tolerance (float): Absolute angle tolerance in checking
consistency of input unit cells against the median unit cell.
"""
self.input_intensities = intensities
uc_params = [flex.double() for i in range(6)]
for d in self.input_intensities:
for i, p in enumerate(d.unit_cell().parameters()):
uc_params[i].append(p)
self.median_unit_cell = uctbx.unit_cell(
parameters=[flex.median(p) for p in uc_params]
)
self._check_unit_cell_consistency(
relative_length_tolerance, absolute_angle_tolerance
)
self.intensities = self.input_intensities[0]
self.dataset_ids = flex.double(self.intensities.size(), 0)
for i, d in enumerate(self.input_intensities[1:]):
self.intensities = self.intensities.concatenate(
d, assert_is_similar_symmetry=False
)
self.dataset_ids.extend(flex.double(d.size(), i + 1))
self.intensities = self.intensities.customized_copy(
unit_cell=self.median_unit_cell
)
self.intensities.set_observation_type_xray_intensity()
sys_absent_flags = self.intensities.sys_absent_flags(integral_only=True).data()
self.intensities = self.intensities.select(~sys_absent_flags)
self.dataset_ids = self.dataset_ids.select(~sys_absent_flags)
self.cb_op_inp_min = self.intensities.change_of_basis_op_to_niggli_cell()
self.intensities = (
self.intensities.change_basis(self.cb_op_inp_min)
.customized_copy(space_group_info=sgtbx.space_group_info("P1"))
.map_to_asu()
.set_info(self.intensities.info())
)
self.lattice_symmetry_max_delta = lattice_symmetry_max_delta
self.subgroups = metric_subgroups(
self.intensities.crystal_symmetry(),
max_delta=self.lattice_symmetry_max_delta,
bravais_types_only=False,
)
self.cb_op_min_best = self.subgroups.result_groups[0]["cb_op_inp_best"]
self.lattice_group = self.subgroups.result_groups[0][
"best_subsym"
].space_group()
self.lattice_group = self.lattice_group.change_basis(
self.cb_op_min_best.inverse()
).make_tidy()
self.patterson_group = (
self.lattice_group.build_derived_patterson_group().make_tidy()
)
logger.info("Patterson group: %s" % self.patterson_group.info())
sel = self.patterson_group.epsilon(self.intensities.indices()) == 1
self.intensities = self.intensities.select(sel)
self.dataset_ids = self.dataset_ids.select(sel)
# Correct SDs by "typical" SD factors
self._correct_sigmas(sd_fac=2.0, sd_b=0.0, sd_add=0.03)
self._resolution_filter(d_min, min_i_mean_over_sigma_mean, min_cc_half)
def _check_unit_cell_consistency(
self, relative_length_tolerance, absolute_angle_tolerance
):
for d in self.input_intensities:
if (
relative_length_tolerance is not None
and absolute_angle_tolerance is not None
):
assert d.unit_cell().is_similar_to(
self.median_unit_cell,
relative_length_tolerance,
absolute_angle_tolerance,
), (str(d.unit_cell()), str(self.median_unit_cell))
def _normalise(self):
if normalisation is None:
return
elif normalisation == "kernel":
normalise = self.kernel_normalisation
elif normalisation == "quasi":
normalise = self.quasi_normalisation
elif normalisation == "ml_iso":
normalise = self.ml_iso_normalisation
elif normalisation == "ml_aniso":
normalise = self.ml_aniso_normalisation
for i in range(int(flex.max(self.dataset_ids) + 1)):
logger.info("Normalising intensities for dataset %i" % (i + 1))
sel = self.dataset_ids == i
intensities = self.intensities.select(self.dataset_ids == i)
if i == 0:
normalised_intensities = normalise(intensities)
else:
normalised_intensities = normalised_intensities.concatenate(
normalise(intensities)
)
self.intensities = normalised_intensities.set_info(
self.intensities.info()
).set_observation_type_xray_intensity()
def _correct_sigmas(self, sd_fac, sd_b, sd_add):
# sd' = SDfac * Sqrt(sd^2 + SdB * I + (SDadd * I)^2)
sigmas = sd_fac * flex.sqrt(
flex.pow2(
self.intensities.sigmas()
+ (sd_b * self.intensities.data())
+ flex.pow2(sd_add * self.intensities.data())
)
)
variance = flex.pow2(self.intensities.sigmas())
si2 = flex.pow2(sd_add * self.intensities.data())
ssc = variance + sd_b * self.intensities.data() + si2
MINVARINFRAC = 0.1
ssc.set_selected(ssc < MINVARINFRAC * variance, MINVARINFRAC * variance)
sd = sd_fac * flex.sqrt(ssc)
self.intensities = self.intensities.customized_copy(sigmas=sd).set_info(
self.intensities.info()
)
[docs] @staticmethod
def kernel_normalisation(intensities):
"""Kernel normalisation of the input intensities.
Args:
intensities (cctbx.miller.array): The intensities to be normalised.
Returns:
cctbx.miller.array: The normalised intensities.
"""
normalisation = absolute_scaling.kernel_normalisation(
intensities, auto_kernel=True
)
return normalisation.normalised_miller.deep_copy().set_info(intensities.info())
[docs] @staticmethod
def quasi_normalisation(intensities):
"""Quasi-normalisation of the input intensities.
Args:
intensities (cctbx.miller.array): The intensities to be normalised.
Returns:
cctbx.miller.array: The normalised intensities.
"""
# handle negative reflections to minimise effect on mean I values.
intensities.data().set_selected(intensities.data() < 0.0, 0.0)
# set up binning objects
if intensities.size() > 20000:
n_refl_shells = 20
elif intensities.size() > 15000:
n_refl_shells = 15
else:
n_refl_shells = 10
d_star_sq = intensities.d_star_sq().data()
step = (flex.max(d_star_sq) - flex.min(d_star_sq) + 1e-8) / n_refl_shells
binner = intensities.setup_binner_d_star_sq_step(d_star_sq_step=step)
normalisations = intensities.intensity_quasi_normalisations()
return intensities.customized_copy(
data=(intensities.data() / normalisations.data()),
sigmas=(intensities.sigmas() / normalisations.data()),
)
[docs] @staticmethod
def ml_aniso_normalisation(intensities):
"""Anisotropic maximum-likelihood normalisation of the input intensities.
Args:
intensities (cctbx.miller.array): The intensities to be normalised.
Returns:
cctbx.miller.array: The normalised intensities.
"""
return symmetry_base._ml_normalisation(intensities, aniso=True)
[docs] @staticmethod
def ml_iso_normalisation(intensities):
"""Isotropic maximum-likelihood normalisation of the input intensities.
Args:
intensities (cctbx.miller.array): The intensities to be normalised.
Returns:
cctbx.miller.array: The normalised intensities.
"""
return symmetry_base._ml_normalisation(intensities, aniso=False)
@staticmethod
def _ml_normalisation(intensities, aniso):
# estimate number of residues per unit cell
mr = matthews.matthews_rupp(intensities.crystal_symmetry())
n_residues = mr.n_residues
# estimate B-factor and scale factors for normalisation
if aniso:
normalisation = absolute_scaling.ml_aniso_absolute_scaling(
intensities, n_residues=n_residues
)
u_star = normalisation.u_star
else:
normalisation = absolute_scaling.ml_iso_absolute_scaling(
intensities, n_residues=n_residues
)
u_star = adptbx.b_as_u(
adptbx.u_iso_as_u_star(intensities.unit_cell(), normalisation.b_wilson)
)
# record output in log file
if aniso:
b_cart = normalisation.b_cart
logger.info("ML estimate of overall B_cart value:")
logger.info(
"""\
%5.2f, %5.2f, %5.2f
%12.2f, %5.2f
%19.2f"""
% (b_cart[0], b_cart[3], b_cart[4], b_cart[1], b_cart[5], b_cart[2])
)
else:
logger.info("ML estimate of overall B value:")
logger.info(" %5.2f A**2" % normalisation.b_wilson)
logger.info("ML estimate of -log of scale factor:")
logger.info(" %5.2f" % (normalisation.p_scale))
s = StringIO()
mr.show(out=s)
normalisation.show(out=s)
logger.debug(s.getvalue())
# apply scales
return intensities.customized_copy(
data=scaling.ml_normalise_aniso(
intensities.indices(),
intensities.data(),
normalisation.p_scale,
intensities.unit_cell(),
u_star,
),
sigmas=scaling.ml_normalise_aniso(
intensities.indices(),
intensities.sigmas(),
normalisation.p_scale,
intensities.unit_cell(),
u_star,
),
)
def _resolution_filter(self, d_min, min_i_mean_over_sigma_mean, min_cc_half):
if d_min is libtbx.Auto and (
min_i_mean_over_sigma_mean is not None or min_cc_half is not None
):
from dials.util import Resolutionizer
rparams = Resolutionizer.phil_defaults.extract().resolutionizer
rparams.nbins = 20
rparams.plot = False
resolutionizer = Resolutionizer.resolutionizer(self.intensities, rparams)
d_min_isigi = 0
d_min_cc_half = 0
if min_i_mean_over_sigma_mean is not None:
try:
d_min_isigi = resolutionizer.resolution_i_mean_over_sigma_mean(
min_i_mean_over_sigma_mean
)
except RuntimeError as e:
logger.info(
"I/sigI resolution filter failed with the following error:"
)
logger.error(e)
else:
logger.info(
"Resolution estimate from <I>/<sigI> > %.1f : %.2f"
% (min_i_mean_over_sigma_mean, d_min_isigi)
)
if min_cc_half is not None:
try:
d_min_cc_half = resolutionizer.resolution_cc_half(min_cc_half)
except RuntimeError as e:
logger.info(
"CChalf resolution filter failed with the following error:"
)
logger.error(e)
else:
logger.info(
"Resolution estimate from CC1/2 > %.2f: %.2f"
% (min_cc_half, d_min_cc_half)
)
valid_d_mins = list(set([d_min_cc_half, d_min_isigi]).difference(set([0])))
if valid_d_mins:
d_min = min(valid_d_mins)
logger.info("High resolution limit set to: %.2f" % d_min)
if d_min is not None:
sel = self.intensities.resolution_filter_selection(d_min=d_min)
self.intensities = self.intensities.select(sel).set_info(
self.intensities.info()
)
self.dataset_ids = self.dataset_ids.select(sel)
logger.info(
"Selecting %i reflections with d > %.2f"
% (self.intensities.size(), d_min)
)
