This documentation page refers to a previous release of DIALS (1.14).
Click here to go to the corresponding page for the latest version of DIALS

Source code for dials.algorithms.indexing.compare_orientation_matrices

from __future__ import absolute_import, division, print_function

import math

from scitbx import matrix
from scitbx.math import r3_rotation_axis_and_angle_from_matrix


[docs]def difference_rotation_matrix_axis_angle(crystal_a, crystal_b, target_angle=0): from cctbx import sgtbx # assert crystal_a.get_space_group() == crystal_b.get_space_group() space_group = crystal_b.get_space_group() best_R_ab = None best_cb_op = None best_axis = None best_angle = 1e8 # iterate over space group ops to find smallest differences for i_op, op in enumerate(space_group.build_derived_laue_group().all_ops()): if op.r().determinant() < 0: continue elif not op.t().is_zero(): continue cb_op = sgtbx.change_of_basis_op(op.inverse()) crystal_b_sym = crystal_b.change_basis(cb_op) U_a = matrix.sqr(crystal_a.get_U()) U_b = matrix.sqr(crystal_b_sym.get_U()) assert U_a.is_r3_rotation_matrix() assert U_b.is_r3_rotation_matrix() # the rotation matrix to transform from U_a to U_b R_ab = U_b * U_a.transpose() axis_angle = r3_rotation_axis_and_angle_from_matrix(R_ab) axis = axis_angle.axis angle = axis_angle.angle() * 180.0 / math.pi for sign in (+1, -1): if abs(sign * angle - target_angle) < abs(best_angle - target_angle): best_angle = sign * angle best_axis = tuple(sign * a for a in axis) best_R_ab = R_ab if sign > 0 else R_ab.inverse() best_cb_op = cb_op if sign > 0 else cb_op.inverse() return best_R_ab, best_axis, best_angle, best_cb_op
[docs]def show_rotation_matrix_differences( crystal_models, out=None, miller_indices=None, comparison="pairwise" ): if out is None: import sys out = sys.stdout print( "Warning: Use of the .show_rotation_matrix_differences() method is deprecated.\n" "Use print(rotation_matrix_differences()) instead.", file=out, ) print( rotation_matrix_differences( crystal_models, miller_indices=miller_indices, comparison=comparison ), file=out, )
[docs]def rotation_matrix_differences( crystal_models, miller_indices=None, comparison="pairwise" ): assert comparison in ("pairwise", "sequential"), comparison output = [] for i in range(len(crystal_models)): for j in range(i + 1, len(crystal_models)): if comparison == "sequential" and j > i + 1: break R_ij, axis, angle, cb_op = difference_rotation_matrix_axis_angle( crystal_models[i], crystal_models[j] ) output.append("Change of basis op: %s" % cb_op) output.append( "Rotation matrix to transform crystal %i to crystal %i:" % (i + 1, j + 1) ) output.append(R_ij.mathematica_form(format="%.3f", one_row_per_line=True)) output.append( ("Rotation of %.3f degrees " % angle) + ("about axis (%.3f, %.3f, %.3f)" % axis) ) if miller_indices is not None: for hkl in miller_indices: cm_i = crystal_models[i] cm_j = crystal_models[j].change_basis(cb_op) A_i = cm_i.get_A() A_j = cm_j.get_A() a_star_i = matrix.col(A_i[0:3]) b_star_i = matrix.col(A_i[3:6]) c_star_i = matrix.col(A_i[6:9]) a_star_j = matrix.col(A_j[0:3]) b_star_j = matrix.col(A_j[3:6]) c_star_j = matrix.col(A_j[6:9]) v_i = hkl[0] * a_star_i + hkl[1] * b_star_i + hkl[2] * c_star_i v_j = hkl[0] * a_star_j + hkl[1] * b_star_j + hkl[2] * c_star_j output.append( ("(%i,%i,%i): " % hkl) + ("%.2f deg" % v_i.angle(v_j, deg=True)) ) output.append("") return "\n".join(output)