Allow dials.find_rotation_axis to work with multi-axis goniometers (d…

…ials#2658)

* For a multi-axis goniometer we set the overall rotation axis by altering
the orientation of the base axis. That way we are not changing the
relationship between axes of the goniometer (so we still trust the header
angles), we are just rotating the whole goniometer to fit the observations

newsfragments/2658.bugfix

@@ -0,0 +1 @@
+``dials.find_rotation_axis``: Correctly set the orientation of the rotation axis for a multi-axis goniometer.

src/dials/command_line/find_rotation_axis.py

@@ -26,6 +26,7 @@
 import dxtbx.flumpy as flumpy
 import libtbx
 import libtbx.phil
+from scitbx import matrix
 
 import dials.util
 import dials.util.log
@@ -522,7 +523,27 @@ def run(args=None, phil=phil_scope):
     logger.info(
         f"\nRotation axis found: {azimuth_deg:.2f} deg. / {azimuth_rad:.3f} rad."
     )
-    expt.goniometer.set_rotation_axis(rotation_axis_to_xyz(azimuth_rad))
+
+    current_rotation_axis = matrix.col(expt.goniometer.get_rotation_axis())
+    new_rotation_axis = matrix.col(rotation_axis_to_xyz(azimuth_rad))
+
+    # Rotation matrix to take the current_rotation_axis to the new_rotation_axis
+    angle = current_rotation_axis.angle(new_rotation_axis)
+    if angle:
+        axis = matrix.col(current_rotation_axis.cross(new_rotation_axis)).normalize()
+        R = axis.axis_and_angle_as_r3_rotation_matrix(angle)
+    else:
+        R = matrix.sqr(matrix.identity(3))
+
+    try:
+        # For a multi-axis goniometer, rotate the base axis by R
+        axes = expt.goniometer.get_axes()
+        base_axis = matrix.col(axes[-1])
+        axes[-1] = R * base_axis
+        expt.goniometer.set_axes(axes)
+    except AttributeError:
+        # For a single-axis goniometer, just set the new_rotation_axis
+        expt.goniometer.set_rotation_axis(new_rotation_axis)
     logger.info(str(expt.goniometer))
 
     logger.info(f"Saving optimised experiments to {params.output.experiments}")

tests/command_line/test_find_rotation_axis.py

@@ -1,7 +1,9 @@
 from __future__ import annotations
 
+from dxtbx.model.goniometer import GoniometerFactory
 from dxtbx.serialize import load
 from scitbx import matrix
+from scitbx.array_family import flex
 
 import dials.command_line.find_rotation_axis as dials_find_rotation_axis
 
@@ -22,3 +24,32 @@ def test_find_rotation_axis(dials_data, run_in_tmp_path):
     expected = matrix.col((-0.627963, -0.778243, 0))
 
     assert axis.angle(expected) < 1e-6
+
+
+def test_find_rotation_axis_multi_axis_goniometer(dials_data, run_in_tmp_path):
+    myd88 = dials_data("MyD88_processed", pathlib=True)
+    experiments = load.experiment_list(str(myd88 / "imported.expt"), check_format=False)
+
+    # Replace single-axis goniometer with multi-axis goniometer
+    axes = flex.vec3_double([(0, -1, 0), (0, -0.642788, 0.766044), (0, -1, 0)])
+    angles = flex.double((0, 0, -62))
+    names = flex.std_string(("PHI", "KAPPA=CHI", "OMEGA"))
+    scan_axis = 2
+    experiments[0].goniometer = GoniometerFactory.multi_axis(
+        axes, angles, names, scan_axis
+    )
+    experiments.as_file("modified.expt")
+
+    # Now run find_rotation_axis, and expect the same result.
+    dials_find_rotation_axis.run(args=["modified.expt", str(myd88 / "strong.refl")])
+    assert run_in_tmp_path.joinpath("optimised.expt").is_file()
+    assert run_in_tmp_path.joinpath("find_rotation_axis-histogram.png").is_file()
+    assert run_in_tmp_path.joinpath("find_rotation_axis-projection.png").is_file()
+    experiments = load.experiment_list(
+        run_in_tmp_path / "optimised.expt", check_format=False
+    )
+
+    axis = matrix.col(experiments[0].goniometer.get_rotation_axis())
+    expected = matrix.col((-0.627963, -0.778243, 0))
+
+    assert axis.angle(expected) < 1e-6