"""Methods for symmetry determination.
This module provides a base class for symmetry determination algorithms.
"""
import logging
logger = logging.getLogger(__name__)
from io import StringIO
import libtbx
from cctbx import adptbx, sgtbx, uctbx
from cctbx.sgtbx.lattice_symmetry import metric_subgroups
from mmtbx import scaling
from mmtbx.scaling import absolute_scaling, matthews
from scitbx.array_family import flex
from dials.util import resolution_analysis
from dials.util.normalisation import quasi_normalisation
[docs]class symmetry_base:
"""Base class for symmetry analysis."""
[docs] 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,
best_monoclinic_beta=True,
):
"""Initialise a symmetry_base object.
Args:
intensities (cctbx.miller.array): The intensities on which to perform
symmetry analysis.
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 :math:`|I|/|sigma(I)|` for
automatic determination of resolution cutoff.
min_cc_half (float): minimum value of CC½ 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.
best_monoclinic_beta (bool): If True, then for monoclinic centered cells, I2
will be preferred over C2 if it gives a less oblique cell (i.e. smaller
beta angle).
"""
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.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,
best_monoclinic_beta=best_monoclinic_beta,
)
self.cb_op_inp_min = self.subgroups.cb_op_inp_minimum
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_group = (
self.subgroups.result_groups[0]["subsym"].space_group().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._normalise(normalisation)
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
):
if not d.unit_cell().is_similar_to(
self.median_unit_cell,
relative_length_tolerance,
absolute_angle_tolerance,
):
raise ValueError(
f"Incompatible unit cell: {d.unit_cell()}\n"
+ f" median unit cell: {self.median_unit_cell}"
)
def _normalise(self, method):
if method is None:
return
elif method == "kernel":
normalise = self.kernel_normalisation
elif method == "quasi":
normalise = quasi_normalisation
elif method == "ml_iso":
normalise = self.ml_iso_normalisation
elif method == "ml_aniso":
normalise = self.ml_aniso_normalisation
for i in range(int(flex.max(self.dataset_ids) + 1)):
logger.info("\n" + "-" * 80 + "\n")
logger.info("Normalising intensities for dataset %i\n", i + 1)
intensities = self.intensities.select(self.dataset_ids == i)
try:
intensities = normalise(intensities)
# Catch any of the several errors that can occur when there are too few
# reflections for the selected normalisation routine.
except (AttributeError, IndexError, RuntimeError):
logger.warning(
"A problem occurred when trying to normalise the intensities by "
"the %s method.\n"
"There may be too few unique reflections. Resorting to using "
"un-normalised intensities instead.",
method,
exc_info=True,
)
return
if i == 0:
normalised_intensities = intensities
else:
normalised_intensities = normalised_intensities.concatenate(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)
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 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):
logger.info("\n" + "-" * 80 + "\n")
logger.info("Estimation of resolution for Laue group analysis\n")
if d_min is libtbx.Auto and (
min_i_mean_over_sigma_mean is not None or min_cc_half is not None
):
d_min = resolution_filter_from_array(
self.intensities, min_i_mean_over_sigma_mean, min_cc_half
)
if d_min is not None:
logger.info("High resolution limit set to: %.2f", d_min)
else:
logger.info("High resolution limit set to: None")
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
)
[docs]def resolution_filter_from_array(intensities, min_i_mean_over_sigma_mean, min_cc_half):
"""Run the resolution filter using miller array data format."""
rparams = resolution_analysis.phil_defaults.extract().resolution
resolutionizer = resolution_analysis.Resolutionizer(intensities, rparams)
return _resolution_filter(resolutionizer, min_i_mean_over_sigma_mean, min_cc_half)
[docs]def resolution_filter_from_reflections_experiments(
reflections, experiments, min_i_mean_over_sigma_mean, min_cc_half
):
"""Run the resolution filter using native dials data formats."""
rparams = resolution_analysis.phil_defaults.extract().resolution
resolutionizer = (
resolution_analysis.Resolutionizer.from_reflections_and_experiments(
reflections, experiments, rparams
)
)
return _resolution_filter(resolutionizer, min_i_mean_over_sigma_mean, min_cc_half)
def _resolution_filter(resolutionizer, min_i_mean_over_sigma_mean, min_cc_half):
"""Use a resolutionizer to perform filtering for symmetry analysis."""
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(
resolution_analysis.metrics.I_MEAN_OVER_SIGMA_MEAN,
limit=min_i_mean_over_sigma_mean,
).d_min
except RuntimeError as e:
logger.info("I/σ(I) resolution filter failed with the following error:")
logger.error(e)
else:
if d_min_isigi:
logger.info(
"Resolution estimate from <I>/<σ(I)> > %.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(
resolution_analysis.metrics.CC_HALF, limit=min_cc_half
).d_min
except RuntimeError as e:
logger.info("CC½ resolution filter failed with the following error:")
logger.error(e)
else:
if d_min_cc_half:
logger.info(
"Resolution estimate from CC½ > %.2f: %.2f",
min_cc_half,
d_min_cc_half,
)
valid = [d for d in (d_min_cc_half, d_min_isigi) if d]
if valid:
return min(valid)
