Source code for dials.algorithms.indexing.compare_orientation_matrices

import math

from scitbx import matrix
from scitbx.math import r3_rotation_axis_and_angle_from_matrix


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


def rotation_matrix_differences(
    crystal_models, miller_indices=None, comparison="pairwise"
):
    assert comparison in ("pairwise", "sequential"), comparison
    output = []
    for i, cm_i in enumerate(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(
                cm_i, crystal_models[j]
            )
            output.append(f"Change of basis op: {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(
                f"Rotation of {angle:.3f} degrees about axis ({axis[0]:.3f}, {axis[1]:.3f}, {axis[2]:.3f})"
            )
            if miller_indices is not None:
                for hkl in miller_indices:
                    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(
                        f"({hkl[0]},{hkl[1]},{hkl[2]}): %.2f deg"
                        % v_i.angle(v_j, deg=True)
                    )
            output.append("")
    return "\n".join(output)