Source code for dials.algorithms.indexing.basis_vector_search.combinations

import logging
import math

from dxtbx.model import Crystal
from scitbx.array_family import flex

from dials.algorithms.indexing.compare_orientation_matrices import (
    difference_rotation_matrix_axis_angle,
)
from dials.algorithms.indexing.symmetry import find_matching_symmetry

logger = logging.getLogger(__name__)


[docs]def candidate_orientation_matrices(basis_vectors, max_combinations=None): # select unique combinations of input vectors to test # the order of combinations is such that combinations comprising vectors # nearer the beginning of the input list will appear before combinations # comprising vectors towards the end of the list n = len(basis_vectors) # hardcoded limit on number of vectors, fixes issue #72 # https://github.com/dials/dials/issues/72 n = min(n, 100) basis_vectors = basis_vectors[:n] combinations = flex.vec3_int(flex.nested_loop((n, n, n))) combinations = combinations.select( flex.sort_permutation(combinations.as_vec3_double().norms()) ) # select only those combinations where j > i and k > j i, j, k = combinations.as_vec3_double().parts() sel = flex.bool(len(combinations), True) sel &= j > i sel &= k > j combinations = combinations.select(sel) if max_combinations is not None and max_combinations < len(combinations): combinations = combinations[:max_combinations] half_pi = 0.5 * math.pi min_angle = 20 / 180 * math.pi # 20 degrees, arbitrary cutoff for i, j, k in combinations: a = basis_vectors[i] b = basis_vectors[j] angle = a.angle(b) if angle < min_angle or (math.pi - angle) < min_angle: continue a_cross_b = a.cross(b) gamma = a.angle(b) if gamma < half_pi: # all angles obtuse if possible please b = -b a_cross_b = -a_cross_b c = basis_vectors[k] if abs(half_pi - a_cross_b.angle(c)) < min_angle: continue alpha = b.angle(c, deg=True) if alpha < half_pi: c = -c if a_cross_b.dot(c) < 0: # we want right-handed basis set, therefore invert all vectors a = -a b = -b c = -c model = Crystal(a, b, c, space_group_symbol="P 1") uc = model.get_unit_cell() cb_op_to_niggli = uc.change_of_basis_op_to_niggli_cell() model = model.change_basis(cb_op_to_niggli) uc = model.get_unit_cell() params = uc.parameters() if uc.volume() > (params[0] * params[1] * params[2] / 100): # unit cell volume cutoff from labelit 2004 paper yield model
[docs]def filter_known_symmetry( crystal_models, target_symmetry, relative_length_tolerance=0.1, absolute_angle_tolerance=5, max_delta=5, ): """Filter crystal models for known symmetry. Args: crystal_models (list): A list of :class:`dxtbx.model.Crystal` objects. target_symmetry (cctbx.crystal.symmetry): The target symmetry for filtering. relative_length_tolerance (float): Relative tolerance for unit cell lengths in unit cell comparison (default value is 0.1). absolute_angle_tolerance (float): Angular tolerance (in degrees) in unit cell comparison (default value is 5). max_delta (float): Maximum allowed Le Page delta used in searching for basis vector combinations that are consistent with the given symmetry (default value is 5). """ n_matched = 0 cb_op_ref_to_primitive = target_symmetry.change_of_basis_op_to_primitive_setting() if target_symmetry.unit_cell() is not None: target_symmetry_primitive = target_symmetry.change_basis(cb_op_ref_to_primitive) else: target_symmetry_primitive = target_symmetry.customized_copy( space_group_info=target_symmetry.space_group_info().change_basis( cb_op_ref_to_primitive ) ) for model in crystal_models: uc = model.get_unit_cell() best_subgroup = find_matching_symmetry( uc, target_symmetry_primitive.space_group(), max_delta=max_delta ) if best_subgroup is not None: if target_symmetry.unit_cell() is not None and not ( best_subgroup["best_subsym"] .unit_cell() .is_similar_to( target_symmetry.as_reference_setting().best_cell().unit_cell(), relative_length_tolerance=relative_length_tolerance, absolute_angle_tolerance=absolute_angle_tolerance, ) ): logger.debug( "Rejecting crystal model inconsistent with input symmetry:\n" f" Unit cell: {str(model.get_unit_cell())}" ) continue n_matched += 1 yield model if not n_matched: logger.warning( "No crystal models remaining after comparing with known symmetry" )
[docs]def filter_similar_orientations( crystal_models, other_crystal_models, minimum_angular_separation=5 ): for cryst in crystal_models: orientation_too_similar = False for cryst_a in other_crystal_models: R_ab, axis, angle, cb_op_ab = difference_rotation_matrix_axis_angle( cryst_a, cryst ) if abs(angle) < minimum_angular_separation: # degrees orientation_too_similar = True break if orientation_too_similar: logger.debug("skipping crystal: too similar to other crystals") continue yield cryst