From df37962dffbac6807e1ad0e0d853ab217c537a55 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 2 Jun 2026 10:54:06 -0400
Subject: [PATCH 01/37] Double quaternion/versor code

---
 hoomd-vector/src/double_quaternion.rs | 31 +++++++++++++++++++++++++++
 hoomd-vector/src/lib.rs               |  2 ++
 2 files changed, 33 insertions(+)
 create mode 100644 hoomd-vector/src/double_quaternion.rs

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
new file mode 100644
index 000000000..24188a542
--- /dev/null
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -0,0 +1,31 @@
+use crate::{Quaternion, Versor};
+
+/// A pair of [`Versor`]s that represent a 4D rotation.
+///
+/// Each [`Versor`] represents an independent rotation about a plane in R^4. We assume
+/// that the left part of the rotation is about the `XY` plane and the right part is
+/// about the `ZW plane`.
+pub struct DoubleVersor {
+    /// The left-isoclinic part of the rotation.
+    l: Versor,
+    /// The right-isoclinic part of the rotation.
+    r: Versor,
+}
+
+impl Default for DoubleVersor {
+    /// Create an identity rotation.
+    ///
+    /// # Example
+    /// ```
+    /// use hoomd_vector::DoubleVersor;
+    ///
+    /// let v = DoubleVersor::default();
+    /// ```
+    #[inline]
+    fn default() -> Self {
+        Self {
+            l: Versor::default(),
+            r: Versor::default(),
+        }
+    }
+}
diff --git a/hoomd-vector/src/lib.rs b/hoomd-vector/src/lib.rs
index 0ffcc3996..be9b3733d 100644
--- a/hoomd-vector/src/lib.rs
+++ b/hoomd-vector/src/lib.rs
@@ -207,6 +207,8 @@ mod cartesian;
 pub mod distribution;
 mod quaternion;
 
+pub mod double_quaternion;
+
 pub use angle::Angle;
 pub use cartesian::{Cartesian, RotationMatrix};
 pub use quaternion::{Quaternion, Versor};

From 6d918a58c6a628f27cbd774755725549428c85f4 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 2 Jun 2026 14:41:37 -0400
Subject: [PATCH 02/37] properly export

---
 hoomd-vector/src/lib.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd-vector/src/lib.rs b/hoomd-vector/src/lib.rs
index be9b3733d..8e33fad58 100644
--- a/hoomd-vector/src/lib.rs
+++ b/hoomd-vector/src/lib.rs
@@ -208,6 +208,7 @@ pub mod distribution;
 mod quaternion;
 
 pub mod double_quaternion;
+pub use double_quaternion::DoubleVersor;
 
 pub use angle::Angle;
 pub use cartesian::{Cartesian, RotationMatrix};

From 3d3fe33294ace7dd7a280a8f1c2b364f51c15816 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 2 Jun 2026 14:41:54 -0400
Subject: [PATCH 03/37] from, l/r, etc

---
 hoomd-vector/src/double_quaternion.rs | 59 ++++++++++++++++++++++++++-
 1 file changed, 58 insertions(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 24188a542..08aced5fe 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,4 +1,4 @@
-use crate::{Quaternion, Versor};
+use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///
@@ -29,3 +29,60 @@ impl Default for DoubleVersor {
         }
     }
 }
+
+impl From<(Versor, Versor)> for DoubleVersor {
+    #[inline]
+    fn from(value: (Versor, Versor)) -> Self {
+        Self {
+            l: value.0,
+            r: value.1,
+        }
+    }
+}
+
+impl DoubleVersor {
+    #[inline]
+    fn left_isoclinic(&self) -> Versor {
+        self.l
+    }
+    #[inline]
+    fn right_isoclinic(&self) -> Versor {
+        self.r
+    }
+}
+
+impl Rotate<Cartesian<4>> for DoubleVersor {
+    type Matrix = RotationMatrix<4>;
+
+    /// Rotate a [`Cartesian<4>`] by a [`DoubleVersor`]
+    ///
+    /// ```math
+    /// \mathbf{q} \vec{a} \mathbf{q}^*
+    /// ```
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use approxim::assert_relative_eq;
+    /// use hoomd_vector::{Cartesian, Rotate, Rotation, DoubleVersor, Versor};
+    /// use std::f64::consts::PI;
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// let a = Cartesian::from([-1.0, 0.0, 0.0, 0.0]);
+    /// let v = DoubleVersor::from_left_isoclinic(
+    ///     Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0)
+    /// );
+    ///
+    /// // Initializing from left isoclinic implies the right isoclinic is unit
+    /// assert_eq!(v.right_isoclinic(), Versor::default());
+    ///
+    /// let b = v.rotate(&a);
+    /// assert_relative_eq!(b, [0.0, -1.0, 0.0, 0.0].into());
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    fn rotate(&self, vector: &Cartesian<4>) -> Cartesian<4> {
+        todo!()
+    }
+}

From 53e18cf08b9ccfa2d85d30cf3890148d459c630b Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 2 Jun 2026 15:57:18 -0400
Subject: [PATCH 04/37] add l/r initialization

---
 hoomd-vector/src/double_quaternion.rs | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 08aced5fe..9f2c48f04 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -49,6 +49,20 @@ impl DoubleVersor {
     fn right_isoclinic(&self) -> Versor {
         self.r
     }
+    #[inline]
+    fn from_left_isoclinic(&self, l: Versor) -> Self {
+        Self {
+            l,
+            r: Versor::default(),
+        }
+    }
+    #[inline]
+    fn from_right_isoclinic(&self, r: Versor) -> Self {
+        Self {
+            l: Versor::default(),
+            r,
+        }
+    }
 }
 
 impl Rotate<Cartesian<4>> for DoubleVersor {

From 791ebbeac84ed4095947a30cbd6d9339989c4ae3 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 9 Jun 2026 16:07:06 -0400
Subject: [PATCH 05/37] pub and docs

---
 hoomd-vector/src/double_quaternion.rs | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 9f2c48f04..4536764f8 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -42,22 +42,26 @@ impl From<(Versor, Versor)> for DoubleVersor {
 
 impl DoubleVersor {
     #[inline]
-    fn left_isoclinic(&self) -> Versor {
+    /// Get the left-isoclinic part of the rotation.
+    pub fn left_isoclinic(&self) -> Versor {
         self.l
     }
     #[inline]
-    fn right_isoclinic(&self) -> Versor {
+    /// Get the right-isoclinic part of the rotation.
+    pub fn right_isoclinic(&self) -> Versor {
         self.r
     }
     #[inline]
-    fn from_left_isoclinic(&self, l: Versor) -> Self {
+    /// Create a purely left-isoclinic double versor.
+    pub fn from_left_isoclinic(&self, l: Versor) -> Self {
         Self {
             l,
             r: Versor::default(),
         }
     }
     #[inline]
-    fn from_right_isoclinic(&self, r: Versor) -> Self {
+    /// Create a purely right-isoclinic double versor.
+    pub fn from_right_isoclinic(&self, r: Versor) -> Self {
         Self {
             l: Versor::default(),
             r,

From 7ba6df81ac409a19fe83e03abb156c101d1e0dc1 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 9 Jun 2026 16:07:24 -0400
Subject: [PATCH 06/37] must_use

---
 hoomd-vector/src/double_quaternion.rs | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 4536764f8..b23813c40 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -42,16 +42,19 @@ impl From<(Versor, Versor)> for DoubleVersor {
 
 impl DoubleVersor {
     #[inline]
+    #[must_use]
     /// Get the left-isoclinic part of the rotation.
     pub fn left_isoclinic(&self) -> Versor {
         self.l
     }
     #[inline]
+    #[must_use]
     /// Get the right-isoclinic part of the rotation.
     pub fn right_isoclinic(&self) -> Versor {
         self.r
     }
     #[inline]
+    #[must_use]
     /// Create a purely left-isoclinic double versor.
     pub fn from_left_isoclinic(&self, l: Versor) -> Self {
         Self {
@@ -60,6 +63,7 @@ impl DoubleVersor {
         }
     }
     #[inline]
+    #[must_use]
     /// Create a purely right-isoclinic double versor.
     pub fn from_right_isoclinic(&self, r: Versor) -> Self {
         Self {

From fb80995320354c26f85890572a2ee4cf7937171d Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 12:13:02 -0400
Subject: [PATCH 07/37] Correct text

---
 hoomd-vector/src/double_quaternion.rs | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index b23813c40..8430873a7 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -2,9 +2,7 @@ use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///
-/// Each [`Versor`] represents an independent rotation about a plane in R^4. We assume
-/// that the left part of the rotation is about the `XY` plane and the right part is
-/// about the `ZW plane`.
+/// Each [`Versor`] represents an independent rotation about a plane in R^4.
 pub struct DoubleVersor {
     /// The left-isoclinic part of the rotation.
     l: Versor,

From e8462c26ac2f089c3ba6d3819b7056246b5a53a6 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 13:05:23 -0400
Subject: [PATCH 08/37] Implement conversions from Matrix44 to
 RotationMatrix<4>

---
 hoomd-vector/src/cartesian.rs         | 9 +++++++++
 hoomd-vector/src/double_quaternion.rs | 2 ++
 2 files changed, 11 insertions(+)

diff --git a/hoomd-vector/src/cartesian.rs b/hoomd-vector/src/cartesian.rs
index 09c2d1c8f..2ba8ba106 100644
--- a/hoomd-vector/src/cartesian.rs
+++ b/hoomd-vector/src/cartesian.rs
@@ -634,6 +634,15 @@ impl<const N: usize> From<RotationMatrix<N>> for Matrix<N, N> {
     }
 }
 
+impl<const N: usize> From<Matrix<N, N>> for RotationMatrix<N> {
+    #[inline]
+    fn from(value: Matrix<N, N>) -> Self {
+        Self {
+            rows: value.rows.map(Cartesian::from),
+        }
+    }
+}
+
 impl<const N: usize> RotationMatrix<N> {
     /// Get the rows of the rotation matrix.
     ///
diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 8430873a7..35fcbf1f5 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,3 +1,5 @@
+use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
+
 use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.

From 04618a7412dcedd1eb4546ceacbc650ef1de5626 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 13:06:03 -0400
Subject: [PATCH 09/37] Implement the direct rotation form for double
 quaternion rotations

---
 hoomd-vector/src/double_quaternion.rs | 48 ++++++++++++++++++++++++---
 1 file changed, 44 insertions(+), 4 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 35fcbf1f5..4e4a10cbc 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -73,13 +73,51 @@ impl DoubleVersor {
     }
 }
 
+impl From<DoubleVersor> for RotationMatrix<4> {
+    /// Construct a rotation matrix equivalent to this double versor's rotation.
+    ///
+    /// When rotating many vectors by the same [`DoubleVersor`], improve performance
+    /// by converting to a matrix first and applying that matrix to the vectors.
+    ///
+    /// # Example
+    /// ```
+    /// use approxim::assert_relative_eq;
+    /// use hoomd_vector::{Cartesian, Rotate, RotationMatrix, Versor};
+    /// use std::f64::consts::PI;
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    fn from(versor: DoubleVersor) -> RotationMatrix<4> {
+        let (&q_l, &q_r) = (
+            versor.left_isoclinic().get(),
+            versor.right_isoclinic().get(),
+        );
+        let (a, [b, c, d]) = (q_l.scalar, q_l.vector.coordinates);
+        let (p, [q, r, s]) = (q_r.scalar, q_l.vector.coordinates);
+
+        // Construct the left-isoclinic matrix L(Q_L)
+        let l_mat = [[a, -b, -c, -d], [b, a, -d, c], [c, d, a, -b], [d, -c, b, a]];
+
+        // Construct the right-isoclinic matrix R(Q_R)
+        let r_mat = [[p, -q, -r, -s], [q, p, s, -r], [r, -s, p, q], [s, r, -q, p]];
+
+        // Combine the left and right isoclinic parts as L@R
+        Matrix44 { rows: l_mat }
+            .matmul(&Matrix44 { rows: r_mat })
+            .into()
+    }
+}
+
 impl Rotate<Cartesian<4>> for DoubleVersor {
     type Matrix = RotationMatrix<4>;
 
-    /// Rotate a [`Cartesian<4>`] by a [`DoubleVersor`]
+    /// Rotate a [`Cartesian<4>`] by a [`DoubleVersor`].
     ///
     /// ```math
-    /// \mathbf{q} \vec{a} \mathbf{q}^*
+    /// \mathbf{q_l} \vec{a} \mathbf{q_r}
     /// ```
     ///
     /// # Example
@@ -95,7 +133,7 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     ///     Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0)
     /// );
     ///
-    /// // Initializing from left isoclinic implies the right isoclinic is unit
+    /// // Initializing from left isoclinic implies the right isoclinic is [1 0 0 0]
     /// assert_eq!(v.right_isoclinic(), Versor::default());
     ///
     /// let b = v.rotate(&a);
@@ -105,6 +143,8 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     /// ```
     #[inline]
     fn rotate(&self, vector: &Cartesian<4>) -> Cartesian<4> {
-        todo!()
+        let q = *self.l.get() * Quaternion::from(vector.coordinates) * *self.r.get();
+        let [x, y, z] = q.vector.coordinates;
+        [q.scalar, x, y, z].into()
     }
 }

From f5307ced4a17e1a09c77b6c8bd4799ae0e23cad4 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 13:09:18 -0400
Subject: [PATCH 10/37] Bugfix

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 4e4a10cbc..0a655bc60 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -96,7 +96,7 @@ impl From<DoubleVersor> for RotationMatrix<4> {
             versor.right_isoclinic().get(),
         );
         let (a, [b, c, d]) = (q_l.scalar, q_l.vector.coordinates);
-        let (p, [q, r, s]) = (q_r.scalar, q_l.vector.coordinates);
+        let (p, [q, r, s]) = (q_r.scalar, q_r.vector.coordinates);
 
         // Construct the left-isoclinic matrix L(Q_L)
         let l_mat = [[a, -b, -c, -d], [b, a, -d, c], [c, d, a, -b], [d, -c, b, a]];

From 4c2cd4d9805c557ae142d585b660f952bb39733f Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 13:12:23 -0400
Subject: [PATCH 11/37] Very basic test for roations

---
 hoomd-vector/src/double_quaternion.rs | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 0a655bc60..e5d88bfe5 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -148,3 +148,24 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
         [q.scalar, x, y, z].into()
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use approxim::assert_relative_eq;
+    use std::f64::consts::PI;
+
+    use crate::{Cartesian, DoubleVersor, Rotate, RotationMatrix, Versor};
+
+    #[test]
+    fn rotation_matrix_matches_direct_rotation() {
+        let l = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into().unwrap(), PI / 3.0);
+        let r = Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), PI / 5.0);
+        let dv = DoubleVersor::from((l, r));
+        let v = Cartesian::from([1.0, 2.0, -3.0, 4.0]);
+
+        let direct = dv.rotate(&v);
+        let via_matrix = RotationMatrix::from(dv).rotate(&v);
+
+        assert_relative_eq!(direct, via_matrix);
+    }
+}

From b086d1a82b3b0f9ddfe415ca8091c8b585f2f2cd Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 14:12:12 -0400
Subject: [PATCH 12/37] Fix from_l/r

---
 hoomd-vector/src/double_quaternion.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index e5d88bfe5..078a8c0cc 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -56,7 +56,7 @@ impl DoubleVersor {
     #[inline]
     #[must_use]
     /// Create a purely left-isoclinic double versor.
-    pub fn from_left_isoclinic(&self, l: Versor) -> Self {
+    pub fn from_left_isoclinic(l: Versor) -> Self {
         Self {
             l,
             r: Versor::default(),
@@ -65,7 +65,7 @@ impl DoubleVersor {
     #[inline]
     #[must_use]
     /// Create a purely right-isoclinic double versor.
-    pub fn from_right_isoclinic(&self, r: Versor) -> Self {
+    pub fn from_right_isoclinic(r: Versor) -> Self {
         Self {
             l: Versor::default(),
             r,

From e36084b4e9d440272a5058c5c5456dc6aa8f8d51 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 14:22:04 -0400
Subject: [PATCH 13/37] WIP on distributions

---
 hoomd-vector/src/double_quaternion.rs | 33 ++++++++++++++++++++++++---
 1 file changed, 30 insertions(+), 3 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 078a8c0cc..9918cc24e 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,4 +1,6 @@
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
+use rand::Rng;
+use rand_distr::{Distribution, StandardUniform};
 
 use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
 
@@ -128,16 +130,18 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     /// use std::f64::consts::PI;
     ///
     /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
-    /// let a = Cartesian::from([-1.0, 0.0, 0.0, 0.0]);
+    /// let a = Cartesian::from([1.0, 2.0, 0.0, 0.0]);
+    ///
+    /// // A rotation of PI/2 radians about the `xy` and `zw` planes
     /// let v = DoubleVersor::from_left_isoclinic(
-    ///     Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0)
+    ///     Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI)
     /// );
     ///
     /// // Initializing from left isoclinic implies the right isoclinic is [1 0 0 0]
     /// assert_eq!(v.right_isoclinic(), Versor::default());
     ///
     /// let b = v.rotate(&a);
-    /// assert_relative_eq!(b, [0.0, -1.0, 0.0, 0.0].into());
+    /// assert_relative_eq!(b, [0.0, 0.0, 2.0, 1.0].into());
     /// # Ok(())
     /// # }
     /// ```
@@ -149,6 +153,29 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     }
 }
 
+impl Distribution<DoubleVersor> for StandardUniform {
+    /// Sample a random [`DoubleVersor`] from the uniform distribution over all rotations in SO(4).
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use hoomd_vector::DoubleVersor;
+    /// use rand::{RngExt, SeedableRng, rngs::StdRng};
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// let mut rng = StdRng::seed_from_u64(1);
+    /// let v: DoubleVersor = rng.random();
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> DoubleVersor {
+        todo!()
+    }
+}
+
+// TODO: implement Rotate and/or Rotation?
+
 #[cfg(test)]
 mod tests {
     use approxim::assert_relative_eq;

From 4d1c5f469e656dde68970f95807a5f89c4e9e0a9 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 15:54:28 -0400
Subject: [PATCH 14/37] Rng

---
 hoomd-vector/src/double_quaternion.rs | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 9918cc24e..057f36e3c 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,5 +1,5 @@
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
-use rand::Rng;
+use rand::{Rng, RngExt};
 use rand_distr::{Distribution, StandardUniform};
 
 use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
@@ -134,14 +134,14 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     ///
     /// // A rotation of PI/2 radians about the `xy` and `zw` planes
     /// let v = DoubleVersor::from_left_isoclinic(
-    ///     Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI)
+    ///     Versor::from_axis_angle([1.0, 0.0, 0.0].try_into()?, PI)
     /// );
     ///
     /// // Initializing from left isoclinic implies the right isoclinic is [1 0 0 0]
     /// assert_eq!(v.right_isoclinic(), Versor::default());
     ///
     /// let b = v.rotate(&a);
-    /// assert_relative_eq!(b, [0.0, 0.0, 2.0, 1.0].into());
+    /// assert_relative_eq!(b, [-2.0, 1.0, 0.0, 0.0].into());
     /// # Ok(())
     /// # }
     /// ```
@@ -156,6 +156,9 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
 impl Distribution<DoubleVersor> for StandardUniform {
     /// Sample a random [`DoubleVersor`] from the uniform distribution over all rotations in SO(4).
     ///
+    /// This is implemented as a random sampling of *pairs* of [`Versor`]'s, which is
+    /// equivalent to a uniform sampling of the full manifold.
+    ///
     /// # Example
     ///
     /// ```
@@ -170,7 +173,7 @@ impl Distribution<DoubleVersor> for StandardUniform {
     /// ```
     #[inline]
     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> DoubleVersor {
-        todo!()
+        (rng.random(), rng.random()).into()
     }
 }
 

From 1fa8f5b32b6cd642771b97411270106fad631a7d Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 16:41:05 -0400
Subject: [PATCH 15/37] Add missing expect

---
 hoomd-vector/src/double_quaternion.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 057f36e3c..952330ad0 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -92,6 +92,7 @@ impl From<DoubleVersor> for RotationMatrix<4> {
     /// # }
     /// ```
     #[inline]
+    #[expect(clippy::many_single_char_names, reason = "Clarity.")]
     fn from(versor: DoubleVersor) -> RotationMatrix<4> {
         let (&q_l, &q_r) = (
             versor.left_isoclinic().get(),

From b4dcce42b68924881ddb8609f4ed4967dd99dee6 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 20:15:41 -0400
Subject: [PATCH 16/37] Try implementing DoubleQuaternion trial moves

---
 hoomd-mc/src/rotate.rs              |  1 +
 hoomd-mc/src/rotate/doubleversor.rs | 77 +++++++++++++++++++++++++++++
 2 files changed, 78 insertions(+)
 create mode 100644 hoomd-mc/src/rotate/doubleversor.rs

diff --git a/hoomd-mc/src/rotate.rs b/hoomd-mc/src/rotate.rs
index 63d4eaf59..a8b948ce0 100644
--- a/hoomd-mc/src/rotate.rs
+++ b/hoomd-mc/src/rotate.rs
@@ -9,6 +9,7 @@ use std::{fmt, marker::PhantomData};
 use hoomd_utility::valid::PositiveReal;
 
 mod angle;
+mod doubleversor;
 mod versor;
 
 /// Change the orientation of a body by a small amount.
diff --git a/hoomd-mc/src/rotate/doubleversor.rs b/hoomd-mc/src/rotate/doubleversor.rs
new file mode 100644
index 000000000..3c1e672d2
--- /dev/null
+++ b/hoomd-mc/src/rotate/doubleversor.rs
@@ -0,0 +1,77 @@
+use std::f64::consts::PI;
+
+use hoomd_microstate::property::Orientation;
+use hoomd_utility::valid::PositiveReal;
+use hoomd_vector::DoubleVersor;
+use rand::{Rng, RngExt};
+use rand_distr::Distribution;
+
+use crate::{Adjust, LocalTrial, Rotate, rotate::versor::VersorDisplacement};
+
+/// A normal distribution of random [`DoubleVersor`]s, centered on a mean with
+/// some standard deviation.
+#[derive(Debug, Clone, Copy)]
+pub(crate) struct DoubleVersorDisplacement {
+    /// The standard deviation of the normal distribution of quaternions around the mean.
+    std_dev: f64,
+}
+
+impl From<f64> for DoubleVersorDisplacement {
+    #[inline]
+    fn from(value: f64) -> Self {
+        Self { std_dev: value }
+    }
+}
+impl Distribution<DoubleVersor> for DoubleVersorDisplacement {
+    /// Sample a random [`DoubleVersor`] displacement from a provided mean.
+    ///
+    /// Mathematically, we sample from a 6-dimensional Normal distribution
+    /// in the tangent space of SO(4), lift to the manifold, then rotate to center on
+    /// the mean of the [`DoubleVersorDisplacement`]. The result is a small displacement
+    /// from a rotational input, with fast decay in the tails that make large
+    /// displacements unlikely. This is desirable for Monte Carlo, as large moves are
+    /// very likely to be rejected.
+    #[inline]
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> DoubleVersor {
+        let single_displacement = VersorDisplacement::from(self.std_dev);
+        // Sample two independent small rotations
+        let q_l = single_displacement.sample(rng);
+        let q_r = single_displacement.sample(rng);
+
+        // Combine them into a DoubleVersor
+        (q_l, q_r).into()
+    }
+}
+
+impl<B> LocalTrial<B> for Rotate<DoubleVersor>
+where
+    B: Orientation<Rotation = DoubleVersor>,
+{
+    /// Perturb a body's orientation by a random amount.
+    #[inline]
+    fn propose<R: Rng>(&self, rng: &mut R, body_properties: B) -> B {
+        let mut trial = body_properties;
+        let displacement = DoubleVersorDisplacement {
+            std_dev: self.maximum_rotation.get(),
+        };
+
+        let delta_quat = displacement.sample(rng);
+        *trial.orientation_mut() = delta_quat.combine(trial.orientation());
+
+        trial
+    }
+}
+
+impl Adjust for Rotate<DoubleVersor> {
+    /// Change the maximum trial move size by the given scale factor.
+    #[inline]
+    fn adjust(&mut self, factor: PositiveReal) {
+        self.maximum_rotation *= factor;
+
+        if self.maximum_rotation.get() > PI / 2.0 {
+            self.maximum_rotation = (PI / 2.0)
+                .try_into()
+                .expect("PI/2.0 should be a positive real");
+        }
+    }
+}

From 9f6bdf7dc5e9136979b220db0931c8354d755305 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 20:27:22 -0400
Subject: [PATCH 17/37] Crate docs

---
 hoomd-vector/src/double_quaternion.rs | 86 +++++++++++++++++++++++++--
 1 file changed, 80 insertions(+), 6 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 952330ad0..d08853837 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,12 +1,18 @@
+//! Implement [`DoubleVersor`], a representation of rotations in four dimensions.
+//! Similar to [`Versor`] in 3D, this approach is more numerically stable and space
+//! efficient than the equivalent matrix representation, but slower when applying
+//! rotations.
+
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
 use rand::{Rng, RngExt};
 use rand_distr::{Distribution, StandardUniform};
 
-use crate::{Cartesian, Quaternion, Rotate, RotationMatrix, Versor};
+use crate::{Cartesian, Quaternion, Rotate, Rotation, RotationMatrix, Versor};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///
 /// Each [`Versor`] represents an independent rotation about a plane in R^4.
+#[derive(Clone, Copy, Debug)]
 pub struct DoubleVersor {
     /// The left-isoclinic part of the rotation.
     l: Versor,
@@ -127,16 +133,17 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     ///
     /// ```
     /// use approxim::assert_relative_eq;
-    /// use hoomd_vector::{Cartesian, Rotate, Rotation, DoubleVersor, Versor};
+    /// use hoomd_vector::{Cartesian, DoubleVersor, Rotate, Rotation, Versor};
     /// use std::f64::consts::PI;
     ///
     /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
     /// let a = Cartesian::from([1.0, 2.0, 0.0, 0.0]);
     ///
     /// // A rotation of PI/2 radians about the `xy` and `zw` planes
-    /// let v = DoubleVersor::from_left_isoclinic(
-    ///     Versor::from_axis_angle([1.0, 0.0, 0.0].try_into()?, PI)
-    /// );
+    /// let v = DoubleVersor::from_left_isoclinic(Versor::from_axis_angle(
+    ///     [1.0, 0.0, 0.0].try_into()?,
+    ///     PI,
+    /// ));
     ///
     /// // Initializing from left isoclinic implies the right isoclinic is [1 0 0 0]
     /// assert_eq!(v.right_isoclinic(), Versor::default());
@@ -178,7 +185,74 @@ impl Distribution<DoubleVersor> for StandardUniform {
     }
 }
 
-// TODO: implement Rotate and/or Rotation?
+impl Rotation for DoubleVersor {
+    /// Combine two rotations.
+    ///
+    /// The resulting versor is obtained by left and right quaternion multiplications.
+    /// ```math
+    /// \mathbf{q}_{l_{ab}} = \mathbf{q}_{l_a} \mathbf{q}_{l_b}
+    /// \mathbf{q}_{r_{ab}} = \mathbf{q}_{r_a} \mathbf{q}_{r_b}
+    /// ```
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use hoomd_vector::{Rotation, Versor};
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    fn combine(&self, other: &Self) -> Self {
+        Self {
+            l: self.l.combine(&other.l),
+            r: self.r.combine(&other.r),
+        }
+    }
+
+    /// Create the identity [`DoubleVersor`]: ([1, [0, 0, 0]], [1, [0, 0, 0]])
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use hoomd_vector::{DoubleVersor, Rotation};
+    ///
+    /// let identity = DoubleVersor::identity();
+    /// ```
+    #[inline]
+    fn identity() -> Self {
+        Self::default()
+    }
+
+    /// Create a [`DoubleVersor`] that performs the inverse rotation of the given double versor.
+    ///
+    /// ```math
+    /// \mathbf{q}^*
+    /// ```
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use hoomd_vector::{DoubleVersor, Rotation};
+    ///
+    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
+    /// let v = DoubleVersor::from_left_isoclinic(Versor::from_axis_angle(
+    ///     [0.0, 1.0, 0.0].try_into()?,
+    ///     1.5,
+    /// ));
+    /// let v_star = v.inverted();
+    /// # Ok(())
+    /// # }
+    /// ```
+    #[inline]
+    fn inverted(self) -> Self {
+        Self {
+            l: self.l.inverted(),
+            r: self.r.inverted(),
+        }
+    }
+}
 
 #[cfg(test)]
 mod tests {

From f137e9dda027d9b6cebdd19f1f00e264c6969b83 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 20:29:17 -0400
Subject: [PATCH 18/37] lint

---
 hoomd-mc/src/rotate/doubleversor.rs   | 5 ++++-
 hoomd-vector/src/double_quaternion.rs | 3 +++
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/hoomd-mc/src/rotate/doubleversor.rs b/hoomd-mc/src/rotate/doubleversor.rs
index 3c1e672d2..4008c4d20 100644
--- a/hoomd-mc/src/rotate/doubleversor.rs
+++ b/hoomd-mc/src/rotate/doubleversor.rs
@@ -1,8 +1,11 @@
+// Copyright (c) 2024-2026 The Regents of the University of Michigan.
+// Part of hoomd-rs, released under the BSD 3-Clause License.
+
 use std::f64::consts::PI;
 
 use hoomd_microstate::property::Orientation;
 use hoomd_utility::valid::PositiveReal;
-use hoomd_vector::DoubleVersor;
+use hoomd_vector::{DoubleVersor, Rotation};
 use rand::{Rng, RngExt};
 use rand_distr::Distribution;
 
diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index d08853837..a92c658bc 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -1,3 +1,6 @@
+// Copyright (c) 2024-2026 The Regents of the University of Michigan.
+// Part of hoomd-rs, released under the BSD 3-Clause License.
+
 //! Implement [`DoubleVersor`], a representation of rotations in four dimensions.
 //! Similar to [`Versor`] in 3D, this approach is more numerically stable and space
 //! efficient than the equivalent matrix representation, but slower when applying

From 3d60d1e6694cbbda02e26d810bc8c6b2453f2cf5 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 20:51:54 -0400
Subject: [PATCH 19/37] Add missing derives

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index a92c658bc..56592e41f 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -15,7 +15,7 @@ use crate::{Cartesian, Quaternion, Rotate, Rotation, RotationMatrix, Versor};
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///
 /// Each [`Versor`] represents an independent rotation about a plane in R^4.
-#[derive(Clone, Copy, Debug)]
+#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
 pub struct DoubleVersor {
     /// The left-isoclinic part of the rotation.
     l: Versor,

From 76dcd7725a0a9a5c58161478af8256ddcfd28629 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 21:00:22 -0400
Subject: [PATCH 20/37] Better tests

---
 hoomd-vector/src/double_quaternion.rs | 55 ++++++++++++++++++++++++++-
 1 file changed, 54 insertions(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 56592e41f..c8ef5959d 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -6,9 +6,11 @@
 //! efficient than the equivalent matrix representation, but slower when applying
 //! rotations.
 
+use approxim::RelativeEq;
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
 use rand::{Rng, RngExt};
 use rand_distr::{Distribution, StandardUniform};
+use serde::{Deserialize, Serialize};
 
 use crate::{Cartesian, Quaternion, Rotate, Rotation, RotationMatrix, Versor};
 
@@ -260,9 +262,11 @@ impl Rotation for DoubleVersor {
 #[cfg(test)]
 mod tests {
     use approxim::assert_relative_eq;
+    use rand::{RngExt, SeedableRng, rngs::StdRng};
+    use rstest::rstest;
     use std::f64::consts::PI;
 
-    use crate::{Cartesian, DoubleVersor, Rotate, RotationMatrix, Versor};
+    use crate::{Cartesian, DoubleVersor, Rotate, Rotation, RotationMatrix, Versor};
 
     #[test]
     fn rotation_matrix_matches_direct_rotation() {
@@ -276,4 +280,53 @@ mod tests {
 
         assert_relative_eq!(direct, via_matrix);
     }
+
+    #[rstest]
+    fn random_rotations_match_matrix(
+        #[values([0.0, 0.0, 0.0, 0.0], [1.0,2.0,3.0,4.0], [-3.0, 9.12, -0.1, 1.25])] v: [f64; 4],
+    ) {
+        let mut rng = StdRng::seed_from_u64(42);
+        let v = Cartesian::from(v);
+
+        for _ in 0..10_000 {
+            let dv: DoubleVersor = rng.random();
+            assert_relative_eq!(
+                dv.rotate(&v),
+                RotationMatrix::from(dv).rotate(&v),
+                epsilon = 1e-14,
+            );
+        }
+    }
+
+    #[rstest]
+    fn random_combine_invert_roundtrip(
+        #[values([0.0, 0.0, 0.0, 0.0], [1.0,2.0,3.0,4.0], [-3.0, 9.12, -0.1, 1.25])] v: [f64; 4],
+    ) {
+        let mut rng = StdRng::seed_from_u64(42);
+        let v = Cartesian::from(v);
+
+        for _ in 0..10_000 {
+            let a: DoubleVersor = rng.random();
+            let b: DoubleVersor = rng.random();
+            let c: DoubleVersor = rng.random();
+
+            // Inverting undoes the rotation
+            assert_relative_eq!(a.inverted().rotate(&a.rotate(&v)), v, epsilon = 1e-14);
+
+            // Combine then invert gives identity
+            let combined = a.combine(&b);
+            assert_relative_eq!(
+                combined.inverted().combine(&combined).rotate(&v),
+                v,
+                epsilon = 1e-12,
+            );
+
+            // Combine is associative
+            assert_relative_eq!(
+                a.combine(&b).combine(&c).rotate(&v),
+                a.combine(&b.combine(&c)).rotate(&v),
+                epsilon = 1e-12,
+            );
+        }
+    }
 }

From ab34576cb7dd516fd7a1b002c10896ad1f39cc47 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 21:00:37 -0400
Subject: [PATCH 21/37] Remove redundant test

---
 hoomd-vector/src/double_quaternion.rs | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index c8ef5959d..43c0a70ab 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -268,19 +268,6 @@ mod tests {
 
     use crate::{Cartesian, DoubleVersor, Rotate, Rotation, RotationMatrix, Versor};
 
-    #[test]
-    fn rotation_matrix_matches_direct_rotation() {
-        let l = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into().unwrap(), PI / 3.0);
-        let r = Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), PI / 5.0);
-        let dv = DoubleVersor::from((l, r));
-        let v = Cartesian::from([1.0, 2.0, -3.0, 4.0]);
-
-        let direct = dv.rotate(&v);
-        let via_matrix = RotationMatrix::from(dv).rotate(&v);
-
-        assert_relative_eq!(direct, via_matrix);
-    }
-
     #[rstest]
     fn random_rotations_match_matrix(
         #[values([0.0, 0.0, 0.0, 0.0], [1.0,2.0,3.0,4.0], [-3.0, 9.12, -0.1, 1.25])] v: [f64; 4],

From 4b9daceb8ed7d219d29516202774e38bf18e9ebd Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 21:03:53 -0400
Subject: [PATCH 22/37] Fix combination convention

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 43c0a70ab..59c014796 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -212,7 +212,7 @@ impl Rotation for DoubleVersor {
     fn combine(&self, other: &Self) -> Self {
         Self {
             l: self.l.combine(&other.l),
-            r: self.r.combine(&other.r),
+            r: other.r.combine(&self.r),
         }
     }
 

From 1129da7c7d5520c6efa5e8f125f21239be4b7528 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Thu, 11 Jun 2026 21:04:27 -0400
Subject: [PATCH 23/37] Test combine then rotate and invert().invert()

---
 hoomd-vector/src/double_quaternion.rs | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 59c014796..735417b67 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -264,7 +264,6 @@ mod tests {
     use approxim::assert_relative_eq;
     use rand::{RngExt, SeedableRng, rngs::StdRng};
     use rstest::rstest;
-    use std::f64::consts::PI;
 
     use crate::{Cartesian, DoubleVersor, Rotate, Rotation, RotationMatrix, Versor};
 
@@ -314,6 +313,16 @@ mod tests {
                 a.combine(&b.combine(&c)).rotate(&v),
                 epsilon = 1e-12,
             );
+
+            // Double inverse is the original rotation
+            assert_relative_eq!(a.inverted().inverted(), a);
+
+            // Combine then rotate matches applying sequentially
+            assert_relative_eq!(
+                a.combine(&b).rotate(&v),
+                a.rotate(&b.rotate(&v)),
+                epsilon = 1e-12,
+            );
         }
     }
 }

From d3ed37f98ca50756d1b7ce48ef7e3b79d549c5cf Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 09:56:42 -0400
Subject: [PATCH 24/37] missing import

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 735417b67..3218fa843 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -239,7 +239,7 @@ impl Rotation for DoubleVersor {
     /// # Example
     ///
     /// ```
-    /// use hoomd_vector::{DoubleVersor, Rotation};
+    /// use hoomd_vector::{DoubleVersor, Rotation, Versor};
     ///
     /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
     /// let v = DoubleVersor::from_left_isoclinic(Versor::from_axis_angle(

From c211b09faf0c4ee059ba12315f6842691254fbea Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 10:47:40 -0400
Subject: [PATCH 25/37] Implement Metric for double versor

---
 hoomd-vector/src/double_quaternion.rs | 77 +++++++++++++++++++++++++--
 1 file changed, 74 insertions(+), 3 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 3218fa843..54a7965d3 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -12,7 +12,9 @@ use rand::{Rng, RngExt};
 use rand_distr::{Distribution, StandardUniform};
 use serde::{Deserialize, Serialize};
 
-use crate::{Cartesian, Quaternion, Rotate, Rotation, RotationMatrix, Versor};
+use crate::{
+    Cartesian, InnerProduct, Metric, Quaternion, Rotate, Rotation, RotationMatrix, Versor,
+};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///
@@ -125,6 +127,11 @@ impl From<DoubleVersor> for RotationMatrix<4> {
     }
 }
 
+fn quaternion_as_cartesian4(q: &Quaternion) -> Cartesian<4> {
+    let [x, y, z] = q.vector.coordinates;
+    [q.scalar, x, y, z].into()
+}
+
 impl Rotate<Cartesian<4>> for DoubleVersor {
     type Matrix = RotationMatrix<4>;
 
@@ -161,8 +168,7 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     #[inline]
     fn rotate(&self, vector: &Cartesian<4>) -> Cartesian<4> {
         let q = *self.l.get() * Quaternion::from(vector.coordinates) * *self.r.get();
-        let [x, y, z] = q.vector.coordinates;
-        [q.scalar, x, y, z].into()
+        quaternion_as_cartesian4(&q)
     }
 }
 
@@ -259,6 +265,71 @@ impl Rotation for DoubleVersor {
     }
 }
 
+impl Metric for DoubleVersor {
+    #[inline]
+    fn distance_squared(&self, other: &Self) -> f64 {
+        // We choose a form based on the dot product of pairs of left and right
+        // quaternions -- as expected, this is exactly the same as the standard
+        // Frobenius metric on SO(N), $`||A-B||_F^2`$. The following Mathematica code
+        // symbolically proves the expressions are equivalent
+        // (* See RotationMatrix::from<DoubleVersor> *)
+        // LMat[{a_, b_, c_, d_}] := {{a, -b, -c, -d}, {b, a, -d, c}, {c, d, a, -b}, {d, -c, b, a}}
+        // RMat[{p_, q_, r_, s_}] := {{p, -q, -r, -s}, {q, p, s, -r}, {r, -s, p, q}, {s, r, -q, p}}
+
+        // qL1 = {a1, b1, c1, d1}; qL2 = {a2, b2, c2, d2};
+        // qR1 = {p1, q1, r1, s1}; qR2 = {p2, q2, r2, s2};
+
+        // (* Construct the explicit SO(4) matrices *)
+        // matA = LMat[qL1] . RMat[qR1];
+        // matB = LMat[qL2] . RMat[qR2];
+
+        // (* Frobenius Norm-based metric *)
+        // (* Note we use ComplexExpand, as otherwise the symbolic algebra gets stuck resolving sqrts *)
+        // explicitFrobeniusSquared = ComplexExpand[Norm[matA - matB, "Frobenius"]^2];
+
+        // (* Quaternion form *)
+        // dotL = qL1 . qL2;
+        // dotR = qR1 . qR2;
+        // algebraicForm = 8 - 8 * dotL * dotR;
+
+        // (* Contrain our solutions to the manifold *)
+        // constraints = {
+        //   a1^2 + b1^2 + c1^2 + d1^2 == 1,
+        //   a2^2 + b2^2 + c2^2 + d2^2 == 1,
+        //   p1^2 + q1^2 + r1^2 + s1^2 == 1,
+        //   p2^2 + q2^2 + r2^2 + s2^2 == 1
+        // };
+
+        // (* Validate *)
+        // Print["Norm[A-B, \"Frobenius\"]^2 == 8 - 8*(qL1.qL2)*(qR1.qR2)"];
+        // isEquivalent = FullSimplify[explicitFrobeniusSquared == algebraicForm, constraints];
+        // Print["Result: ", isEquivalent];
+        let left_dot =
+            quaternion_as_cartesian4(self.l.get()).dot(&quaternion_as_cartesian4(other.l.get()));
+        let right_dot =
+            quaternion_as_cartesian4(self.r.get()).dot(&quaternion_as_cartesian4(other.r.get()));
+        8.0 * (1.0 - left_dot * right_dot)
+    }
+
+    #[inline]
+    /// The dimension of the manifold of SO(N) is $`\frac{N(N-1)}{2}`$, or 6 for SO(4).
+    fn n_dimensions(&self) -> usize {
+        6
+    }
+    #[inline]
+    /// The distance between two [`DoubleQuaternion`] points.
+    ///
+    /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$ on the
+    /// manifold SO(4):
+    ///
+    /// ```math
+    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8.0 * (1.0 - (u_l \cdot v_l * u_r \cdot v_r)}
+    /// ```
+    fn distance(&self, other: &Self) -> f64 {
+        (self.distance_squared(other)).max(0.0).sqrt()
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use approxim::assert_relative_eq;

From 254a3273e5f3b05fe01fbb92a7080d0ddd11db3f Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 10:48:53 -0400
Subject: [PATCH 26/37] More docs

---
 hoomd-vector/src/double_quaternion.rs | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 54a7965d3..9d896e099 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -325,6 +325,9 @@ impl Metric for DoubleVersor {
     /// ```math
     /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8.0 * (1.0 - (u_l \cdot v_l * u_r \cdot v_r)}
     /// ```
+    ///
+    /// This is equivalent to the matrix form $`||R_u - R_v||_F`$, but does not require
+    /// forming the matrix representation of each rotation.
     fn distance(&self, other: &Self) -> f64 {
         (self.distance_squared(other)).max(0.0).sqrt()
     }

From 3f721eac327ed20904864702172505170a707685 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 10:49:14 -0400
Subject: [PATCH 27/37] Remove decimals from math

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 9d896e099..fd039d0ca 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -323,7 +323,7 @@ impl Metric for DoubleVersor {
     /// manifold SO(4):
     ///
     /// ```math
-    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8.0 * (1.0 - (u_l \cdot v_l * u_r \cdot v_r)}
+    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8 * (1 - (u_l \cdot v_l * u_r \cdot v_r)}
     /// ```
     ///
     /// This is equivalent to the matrix form $`||R_u - R_v||_F`$, but does not require

From f4e68bd6146313e8ed48012311ae6dbd11b1e8c7 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 10:53:05 -0400
Subject: [PATCH 28/37] Swap internal helper to a method on quaternion

---
 hoomd-vector/src/double_quaternion.rs | 36 ++++++++++++---------------
 hoomd-vector/src/quaternion.rs        |  8 ++++++
 2 files changed, 24 insertions(+), 20 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index fd039d0ca..e5d1a54b7 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -127,11 +127,6 @@ impl From<DoubleVersor> for RotationMatrix<4> {
     }
 }
 
-fn quaternion_as_cartesian4(q: &Quaternion) -> Cartesian<4> {
-    let [x, y, z] = q.vector.coordinates;
-    [q.scalar, x, y, z].into()
-}
-
 impl Rotate<Cartesian<4>> for DoubleVersor {
     type Matrix = RotationMatrix<4>;
 
@@ -168,7 +163,7 @@ impl Rotate<Cartesian<4>> for DoubleVersor {
     #[inline]
     fn rotate(&self, vector: &Cartesian<4>) -> Cartesian<4> {
         let q = *self.l.get() * Quaternion::from(vector.coordinates) * *self.r.get();
-        quaternion_as_cartesian4(&q)
+        q.embed_in_cartesian_4()
     }
 }
 
@@ -280,34 +275,35 @@ impl Metric for DoubleVersor {
         // qR1 = {p1, q1, r1, s1}; qR2 = {p2, q2, r2, s2};
 
         // (* Construct the explicit SO(4) matrices *)
-        // matA = LMat[qL1] . RMat[qR1];
-        // matB = LMat[qL2] . RMat[qR2];
+        // A = LMat[qL1] . RMat[qR1]; B = LMat[qL2] . RMat[qR2];
 
         // (* Frobenius Norm-based metric *)
         // (* Note we use ComplexExpand, as otherwise the symbolic algebra gets stuck resolving sqrts *)
-        // explicitFrobeniusSquared = ComplexExpand[Norm[matA - matB, "Frobenius"]^2];
+        // explicitFrobeniusSquared = ComplexExpand[Norm[A - B, "Frobenius"]^2];
 
         // (* Quaternion form *)
-        // dotL = qL1 . qL2;
-        // dotR = qR1 . qR2;
-        // algebraicForm = 8 - 8 * dotL * dotR;
+        // algebraicForm = 8 - 8 * (qL1 . qL2) * (qR1 . qR2);
 
         // (* Contrain our solutions to the manifold *)
         // constraints = {
-        //   a1^2 + b1^2 + c1^2 + d1^2 == 1,
-        //   a2^2 + b2^2 + c2^2 + d2^2 == 1,
-        //   p1^2 + q1^2 + r1^2 + s1^2 == 1,
-        //   p2^2 + q2^2 + r2^2 + s2^2 == 1
+        //   a1^2 + b1^2 + c1^2 + d1^2 == 1, a2^2 + b2^2 + c2^2 + d2^2 == 1,
+        //   p1^2 + q1^2 + r1^2 + s1^2 == 1, p2^2 + q2^2 + r2^2 + s2^2 == 1
         // };
 
         // (* Validate *)
         // Print["Norm[A-B, \"Frobenius\"]^2 == 8 - 8*(qL1.qL2)*(qR1.qR2)"];
         // isEquivalent = FullSimplify[explicitFrobeniusSquared == algebraicForm, constraints];
         // Print["Result: ", isEquivalent];
-        let left_dot =
-            quaternion_as_cartesian4(self.l.get()).dot(&quaternion_as_cartesian4(other.l.get()));
-        let right_dot =
-            quaternion_as_cartesian4(self.r.get()).dot(&quaternion_as_cartesian4(other.r.get()));
+        let left_dot = self
+            .l
+            .get()
+            .embed_in_cartesian_4()
+            .dot(&other.l.get().embed_in_cartesian_4());
+        let right_dot = self
+            .r
+            .get()
+            .embed_in_cartesian_4()
+            .dot(&other.r.get().embed_in_cartesian_4());
         8.0 * (1.0 - left_dot * right_dot)
     }
 
diff --git a/hoomd-vector/src/quaternion.rs b/hoomd-vector/src/quaternion.rs
index ed0ecfd4e..c19d02e92 100644
--- a/hoomd-vector/src/quaternion.rs
+++ b/hoomd-vector/src/quaternion.rs
@@ -284,6 +284,14 @@ impl Quaternion {
     pub fn to_versor_unchecked(self) -> Versor {
         Versor(self / self.norm())
     }
+
+    /// Embed this quaternion as a point in R^4.
+    #[inline]
+    #[must_use]
+    pub fn embed_in_cartesian_4(&self) -> Cartesian<4> {
+        let [x, y, z] = self.vector.coordinates;
+        [self.scalar, x, y, z].into()
+    }
 }
 
 impl From<[f64; 4]> for Quaternion {

From 3fe2c1700b93bc4e712df75b3da9ab3c9a0ac279 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:19:01 -0400
Subject: [PATCH 29/37] Properly move to chordal distance

---
 hoomd-vector/src/double_quaternion.rs | 104 +++++++++++++++-----------
 1 file changed, 60 insertions(+), 44 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index e5d1a54b7..0132d6e9b 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -86,6 +86,65 @@ impl DoubleVersor {
             r,
         }
     }
+
+    #[inline]
+    #[must_use]
+    /// The distance between two [`DoubleQuaternion`] points.
+    ///
+    /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$
+    /// *along a chord through the manifold*. This is NOT equivalent to the intrinsic
+    /// distance, defined with [`DoubleQuaternion::Metric`].
+    ///
+    /// ```math
+    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8 * (1 - (u_l \cdot v_l * u_r \cdot v_r)}
+    /// ```
+    ///
+    /// This is equivalent to the matrix form $`||R_u - R_v||_F`$, but does not require
+    /// forming the matrix representation of each rotation.
+    pub fn chordal_distanc(&self, other: &Self) -> f64 {
+        // We choose a form based on the dot product of pairs of left and right
+        // quaternions -- as expected, this is exactly the same as the standard
+        // Frobenius metric on SO(N), $`||A-B||_F^2`$. The following Mathematica code
+        // symbolically proves the expressions are equivalent
+        // (* See RotationMatrix::from<DoubleVersor> *)
+        // LMat[{a_, b_, c_, d_}] := {{a, -b, -c, -d}, {b, a, -d, c}, {c, d, a, -b}, {d, -c, b, a}}
+        // RMat[{p_, q_, r_, s_}] := {{p, -q, -r, -s}, {q, p, s, -r}, {r, -s, p, q}, {s, r, -q, p}}
+
+        // qL1 = {a1, b1, c1, d1}; qL2 = {a2, b2, c2, d2};
+        // qR1 = {p1, q1, r1, s1}; qR2 = {p2, q2, r2, s2};
+
+        // (* Construct the explicit SO(4) matrices *)
+        // A = LMat[qL1] . RMat[qR1]; B = LMat[qL2] . RMat[qR2];
+
+        // (* Frobenius Norm-based metric *)
+        // (* Note we use ComplexExpand, as otherwise the symbolic algebra gets stuck resolving sqrts *)
+        // explicitFrobeniusSquared = ComplexExpand[Norm[A - B, "Frobenius"]^2];
+
+        // (* Quaternion form *)
+        // algebraicForm = 8 - 8 * (qL1 . qL2) * (qR1 . qR2);
+
+        // (* Contrain our solutions to the manifold *)
+        // constraints = {
+        //   a1^2 + b1^2 + c1^2 + d1^2 == 1, a2^2 + b2^2 + c2^2 + d2^2 == 1,
+        //   p1^2 + q1^2 + r1^2 + s1^2 == 1, p2^2 + q2^2 + r2^2 + s2^2 == 1
+        // };
+
+        // (* Validate *)
+        // Print["Norm[A-B, \"Frobenius\"]^2 == 8 - 8*(qL1.qL2)*(qR1.qR2)"];
+        // isEquivalent = FullSimplify[explicitFrobeniusSquared == algebraicForm, constraints];
+        // Print["Result: ", isEquivalent];
+        let left_dot = self
+            .l
+            .get()
+            .embed_in_cartesian_4()
+            .dot(&other.l.get().embed_in_cartesian_4());
+        let right_dot = self
+            .r
+            .get()
+            .embed_in_cartesian_4()
+            .dot(&other.r.get().embed_in_cartesian_4());
+        (8.0 * (1.0 - left_dot * right_dot)).max(0.0).sqrt()
+    }
 }
 
 impl From<DoubleVersor> for RotationMatrix<4> {
@@ -262,50 +321,7 @@ impl Rotation for DoubleVersor {
 
 impl Metric for DoubleVersor {
     #[inline]
-    fn distance_squared(&self, other: &Self) -> f64 {
-        // We choose a form based on the dot product of pairs of left and right
-        // quaternions -- as expected, this is exactly the same as the standard
-        // Frobenius metric on SO(N), $`||A-B||_F^2`$. The following Mathematica code
-        // symbolically proves the expressions are equivalent
-        // (* See RotationMatrix::from<DoubleVersor> *)
-        // LMat[{a_, b_, c_, d_}] := {{a, -b, -c, -d}, {b, a, -d, c}, {c, d, a, -b}, {d, -c, b, a}}
-        // RMat[{p_, q_, r_, s_}] := {{p, -q, -r, -s}, {q, p, s, -r}, {r, -s, p, q}, {s, r, -q, p}}
-
-        // qL1 = {a1, b1, c1, d1}; qL2 = {a2, b2, c2, d2};
-        // qR1 = {p1, q1, r1, s1}; qR2 = {p2, q2, r2, s2};
-
-        // (* Construct the explicit SO(4) matrices *)
-        // A = LMat[qL1] . RMat[qR1]; B = LMat[qL2] . RMat[qR2];
-
-        // (* Frobenius Norm-based metric *)
-        // (* Note we use ComplexExpand, as otherwise the symbolic algebra gets stuck resolving sqrts *)
-        // explicitFrobeniusSquared = ComplexExpand[Norm[A - B, "Frobenius"]^2];
-
-        // (* Quaternion form *)
-        // algebraicForm = 8 - 8 * (qL1 . qL2) * (qR1 . qR2);
-
-        // (* Contrain our solutions to the manifold *)
-        // constraints = {
-        //   a1^2 + b1^2 + c1^2 + d1^2 == 1, a2^2 + b2^2 + c2^2 + d2^2 == 1,
-        //   p1^2 + q1^2 + r1^2 + s1^2 == 1, p2^2 + q2^2 + r2^2 + s2^2 == 1
-        // };
-
-        // (* Validate *)
-        // Print["Norm[A-B, \"Frobenius\"]^2 == 8 - 8*(qL1.qL2)*(qR1.qR2)"];
-        // isEquivalent = FullSimplify[explicitFrobeniusSquared == algebraicForm, constraints];
-        // Print["Result: ", isEquivalent];
-        let left_dot = self
-            .l
-            .get()
-            .embed_in_cartesian_4()
-            .dot(&other.l.get().embed_in_cartesian_4());
-        let right_dot = self
-            .r
-            .get()
-            .embed_in_cartesian_4()
-            .dot(&other.r.get().embed_in_cartesian_4());
-        8.0 * (1.0 - left_dot * right_dot)
-    }
+    fn distance_squared(&self, other: &Self) -> f64 {}
 
     #[inline]
     /// The dimension of the manifold of SO(N) is $`\frac{N(N-1)}{2}`$, or 6 for SO(4).

From 00f25cccb6a0ec8c8fbb7b2ae998514d1bc5c1c9 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:32:17 -0400
Subject: [PATCH 30/37] Make dot as cartesian pub crate

---
 hoomd-vector/src/quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/quaternion.rs b/hoomd-vector/src/quaternion.rs
index c19d02e92..861c7ed2d 100644
--- a/hoomd-vector/src/quaternion.rs
+++ b/hoomd-vector/src/quaternion.rs
@@ -518,7 +518,7 @@ pub struct Versor(Quaternion);
 impl Versor {
     /// Take the dot product of the Versor as an element of $`\mathbb{R}^4`$.
     #[inline]
-    fn dot_as_cartesian(&self, other: &Self) -> f64 {
+    pub(crate) fn dot_as_cartesian(&self, other: &Self) -> f64 {
         self.get().scalar * other.get().scalar + self.get().vector.dot(&other.get().vector)
     }
     /// Create a [`Versor`] that rotates by an angle (in radians)

From 9abe980b525ddfc090bb692630faadff83d0d745 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:32:30 -0400
Subject: [PATCH 31/37] Correct implementation of intrinsic distance

---
 hoomd-vector/src/double_quaternion.rs | 28 +++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 0132d6e9b..f09feae92 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -6,6 +6,8 @@
 //! efficient than the equivalent matrix representation, but slower when applying
 //! rotations.
 
+use std::f64::consts::{PI, TAU};
+
 use approxim::RelativeEq;
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
 use rand::{Rng, RngExt};
@@ -133,16 +135,8 @@ impl DoubleVersor {
         // Print["Norm[A-B, \"Frobenius\"]^2 == 8 - 8*(qL1.qL2)*(qR1.qR2)"];
         // isEquivalent = FullSimplify[explicitFrobeniusSquared == algebraicForm, constraints];
         // Print["Result: ", isEquivalent];
-        let left_dot = self
-            .l
-            .get()
-            .embed_in_cartesian_4()
-            .dot(&other.l.get().embed_in_cartesian_4());
-        let right_dot = self
-            .r
-            .get()
-            .embed_in_cartesian_4()
-            .dot(&other.r.get().embed_in_cartesian_4());
+        let left_dot = self.l.dot_as_cartesian(&other.l);
+        let right_dot = self.r.dot_as_cartesian(&other.r);
         (8.0 * (1.0 - left_dot * right_dot)).max(0.0).sqrt()
     }
 }
@@ -321,7 +315,17 @@ impl Rotation for DoubleVersor {
 
 impl Metric for DoubleVersor {
     #[inline]
-    fn distance_squared(&self, other: &Self) -> f64 {}
+    fn distance_squared(&self, other: &Self) -> f64 {
+        let left_angle = self.l.dot_as_cartesian(&other.l).clamp(-1.0, 1.0).acos();
+        let right_angle = self.r.dot_as_cartesian(&other.r).clamp(-1.0, 1.0).acos();
+
+        // Alpha is the shorter arc of {l+r, 2π-(l+r)}
+        let alpha = f64::min(left_angle + right_angle, TAU - (left_angle + right_angle));
+        let beta = f64::abs(left_angle - right_angle);
+
+        // Multiply by 2.0 to match the standard ||log(R_A^T R_B)||_F^2 scalin
+        2.0 * (alpha.powi(2) + beta.powi(2))
+    }
 
     #[inline]
     /// The dimension of the manifold of SO(N) is $`\frac{N(N-1)}{2}`$, or 6 for SO(4).
@@ -341,7 +345,7 @@ impl Metric for DoubleVersor {
     /// This is equivalent to the matrix form $`||R_u - R_v||_F`$, but does not require
     /// forming the matrix representation of each rotation.
     fn distance(&self, other: &Self) -> f64 {
-        (self.distance_squared(other)).max(0.0).sqrt()
+        self.distance_squared(other).sqrt()
     }
 }
 

From a353d7ddf100272cdce8fb53114985a4d1fd1ec3 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:32:52 -0400
Subject: [PATCH 32/37] lints

---
 hoomd-vector/src/double_quaternion.rs | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index f09feae92..94c223085 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -6,7 +6,7 @@
 //! efficient than the equivalent matrix representation, but slower when applying
 //! rotations.
 
-use std::f64::consts::{PI, TAU};
+use std::f64::consts::TAU;
 
 use approxim::RelativeEq;
 use hoomd_linear_algebra::{MatMul, matrix::Matrix44};
@@ -14,9 +14,7 @@ use rand::{Rng, RngExt};
 use rand_distr::{Distribution, StandardUniform};
 use serde::{Deserialize, Serialize};
 
-use crate::{
-    Cartesian, InnerProduct, Metric, Quaternion, Rotate, Rotation, RotationMatrix, Versor,
-};
+use crate::{Cartesian, Metric, Quaternion, Rotate, Rotation, RotationMatrix, Versor};
 
 /// A pair of [`Versor`]s that represent a 4D rotation.
 ///

From f308eb2dbc3cc5b308640488d506d95497f56701 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:38:25 -0400
Subject: [PATCH 33/37] Correct equations

---
 hoomd-vector/src/double_quaternion.rs | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 94c223085..45682f6e7 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -337,11 +337,19 @@ impl Metric for DoubleVersor {
     /// manifold SO(4):
     ///
     /// ```math
-    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{8 * (1 - (u_l \cdot v_l * u_r \cdot v_r)}
+    /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{2\theta_1^2 + 2\theta_2^2}
     /// ```
     ///
-    /// This is equivalent to the matrix form $`||R_u - R_v||_F`$, but does not require
-    /// forming the matrix representation of each rotation.
+    /// Where the principal planar angles $\theta_1$ and $\theta_2$ are derived
+    /// from the left and right quaternion arc lengths ($\alpha$ and $\beta$),
+    /// accounting double cover of the manifold:
+    /// * $`\alpha = \arccos(\vec{u}_L \cdot \vec{v}_L)`$
+    /// * $`\beta = \arccos(\vec{u}_R \cdot \vec{v}_R)`$
+    /// * $`\theta_1 = \min(\alpha + \beta, 2\pi - (\alpha + \beta))`$
+    /// * $`\theta_2 = |\alpha - \beta|`$
+    ///
+    /// This is equivalent to the scaled matrix logarithm form $`||\log(R_u^T R_v)||_F`$,
+    /// measuring the shortest path along the curved manifold.
     fn distance(&self, other: &Self) -> f64 {
         self.distance_squared(other).sqrt()
     }

From cbf616b4399df81b3b1d4f0d8fe442daae572034 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Fri, 12 Jun 2026 11:56:08 -0400
Subject: [PATCH 34/37] Tests for orienation trial moves in SO(4)

---
 hoomd-mc/src/rotate/doubleversor.rs | 96 ++++++++++++++++++++++++++++-
 1 file changed, 95 insertions(+), 1 deletion(-)

diff --git a/hoomd-mc/src/rotate/doubleversor.rs b/hoomd-mc/src/rotate/doubleversor.rs
index 4008c4d20..823aadd3a 100644
--- a/hoomd-mc/src/rotate/doubleversor.rs
+++ b/hoomd-mc/src/rotate/doubleversor.rs
@@ -6,7 +6,7 @@ use std::f64::consts::PI;
 use hoomd_microstate::property::Orientation;
 use hoomd_utility::valid::PositiveReal;
 use hoomd_vector::{DoubleVersor, Rotation};
-use rand::{Rng, RngExt};
+use rand::Rng;
 use rand_distr::Distribution;
 
 use crate::{Adjust, LocalTrial, Rotate, rotate::versor::VersorDisplacement};
@@ -78,3 +78,97 @@ impl Adjust for Rotate<DoubleVersor> {
         }
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use assert2::check;
+    use hoomd_microstate::property::OrientedPoint;
+    use hoomd_vector::{Cartesian, DoubleVersor, Metric};
+    use rand::{SeedableRng, rngs::StdRng};
+    use rstest::*;
+
+    /// Number of trial moves to test.
+    const N: usize = 262_144;
+
+    /// Expected geodesic distance for small a:E[χ₆] = sqrt2*Γ(7/2)/Γ(3) = 15sqrt(2π)/16
+    #[inline]
+    fn expected_chi6_mean() -> f64 {
+        15.0 * (2.0 * PI).sqrt() / 16.0
+    }
+
+    #[rstest]
+    fn rotate(#[values(0.1, 0.5)] a: f64) {
+        let mut rng = StdRng::seed_from_u64(1);
+        let body = OrientedPoint {
+            position: Cartesian::from([0.0, 0.0, 0.0]),
+            orientation: DoubleVersor::default(),
+        };
+        let rotate = Rotate::with_maximum_rotation(
+            a.try_into()
+                .expect("hard-coded constant should be a positive real"),
+        );
+
+        let mut delta_thetas = Vec::with_capacity(N);
+        for _ in 0..N {
+            let trial = rotate.propose(&mut rng, body);
+            delta_thetas.push(trial.orientation.distance(&body.orientation));
+        }
+
+        let mean = delta_thetas.iter().sum::<f64>() / N as f64;
+        let expected = expected_chi6_mean() * a;
+        assert!(
+            (mean - expected).abs() < 0.01 * a, // Should be within a few percent
+            "mean = {mean}, expected = {expected}, diff = {}",
+            (mean - expected).abs(),
+        );
+    }
+
+    #[rstest]
+    fn rotate_mean_distance(#[values(1e-3, 0.1)] a: f64) {
+        // Two independent versor displacements give distance:
+        // sqrt(|sl|² + |sr|²) = a*\chi_6 for small displacements, so
+        // E[distance] = a * sqrt(2) * (15*sqrt(pi)/8)/2 = a * 15sqrt(2π)/16.
+        let expected = expected_chi6_mean() * a;
+
+        let mut rng = StdRng::seed_from_u64(1);
+        let body = OrientedPoint {
+            position: Cartesian::from([0.0, 0.0, 0.0]),
+            orientation: DoubleVersor::default(),
+        };
+        let rotate = Rotate::with_maximum_rotation(
+            a.try_into()
+                .expect("hard-coded constant should be a positive real"),
+        );
+
+        let sum: f64 = (0..N)
+            .map(|_| {
+                let trial = rotate.propose(&mut rng, body);
+                trial.orientation.distance(&body.orientation)
+            })
+            .sum();
+
+        let mean = sum / N as f64;
+        assert!(
+            (mean - expected).abs() < 0.001 * a,
+            "mean = {mean}, expected = {expected}, diff = {}",
+            (mean - expected).abs(),
+        );
+    }
+
+    #[test]
+    fn test_adjust() -> anyhow::Result<()> {
+        let mut rotate = Rotate::<DoubleVersor>::with_maximum_rotation(0.5.try_into()?);
+
+        rotate.adjust(2.0.try_into()?);
+        check!(rotate.maximum_rotation().get() == 1.0);
+
+        rotate.adjust(0.5.try_into()?);
+        check!(rotate.maximum_rotation().get() == 0.5);
+
+        rotate.adjust(10.0.try_into()?);
+        check!(rotate.maximum_rotation().get() == PI / 2.0);
+
+        Ok(())
+    }
+}

From ea4a3968879252e105c2ea581f8a023eb3c99780 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 30 Jun 2026 13:24:41 -0400
Subject: [PATCH 35/37] lint

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 45682f6e7..905689bbd 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -361,7 +361,7 @@ mod tests {
     use rand::{RngExt, SeedableRng, rngs::StdRng};
     use rstest::rstest;
 
-    use crate::{Cartesian, DoubleVersor, Rotate, Rotation, RotationMatrix, Versor};
+    use crate::{Cartesian, DoubleVersor, Rotate, Rotation, RotationMatrix};
 
     #[rstest]
     fn random_rotations_match_matrix(

From e598e9ee83f2a12a4168ef8996ab35a3b1a73928 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 30 Jun 2026 13:26:12 -0400
Subject: [PATCH 36/37] clippy

---
 hoomd-vector/src/double_quaternion.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd-vector/src/double_quaternion.rs b/hoomd-vector/src/double_quaternion.rs
index 905689bbd..0162a8e0e 100644
--- a/hoomd-vector/src/double_quaternion.rs
+++ b/hoomd-vector/src/double_quaternion.rs
@@ -340,7 +340,7 @@ impl Metric for DoubleVersor {
     /// d_{SO(4)}(\vec{u}, \vec{v}) = \sqrt{2\theta_1^2 + 2\theta_2^2}
     /// ```
     ///
-    /// Where the principal planar angles $\theta_1$ and $\theta_2$ are derived
+    /// Where the principal planar angles $`\theta_1`$ and $`\theta_2`$ are derived
     /// from the left and right quaternion arc lengths ($\alpha$ and $\beta$),
     /// accounting double cover of the manifold:
     /// * $`\alpha = \arccos(\vec{u}_L \cdot \vec{v}_L)`$

From 2c03cfbb32a01f867866668c3c8505db457e25c6 Mon Sep 17 00:00:00 2001
From: janbridley <bradleyjenj@gmail.com>
Date: Tue, 30 Jun 2026 13:28:50 -0400
Subject: [PATCH 37/37] lint

---
 hoomd-mc/src/rotate/doubleversor.rs | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/hoomd-mc/src/rotate/doubleversor.rs b/hoomd-mc/src/rotate/doubleversor.rs
index 823aadd3a..a27dac2f6 100644
--- a/hoomd-mc/src/rotate/doubleversor.rs
+++ b/hoomd-mc/src/rotate/doubleversor.rs
@@ -1,6 +1,8 @@
 // Copyright (c) 2024-2026 The Regents of the University of Michigan.
 // Part of hoomd-rs, released under the BSD 3-Clause License.
 
+//! Implement [`DoubleVersor`] trial moves for 4-dimensional orientable bodies.
+
 use std::f64::consts::PI;
 
 use hoomd_microstate::property::Orientation;