Introduce boost qvm and a Quaternion type (#4031)

Fixes #2289

src/core/Particle.hpp

@@ -25,6 +25,7 @@
 
 #include <utils/Vector.hpp>
 #include <utils/math/quaternion.hpp>
+#include <utils/quaternion.hpp>
 
 #include <boost/serialization/vector.hpp>
 
@@ -128,9 +129,10 @@ struct ParticleProperties {
     int to_particle_id = 0;
     double distance = 0;
     /** Relative position of the virtual site. */
-    Utils::Vector4d rel_orientation = {0., 0., 0., 0.};
+    Utils::Quaternion<double> rel_orientation =
+        Utils::Quaternion<double>::zero();
     /** Orientation of the virtual particle in the body fixed frame. */
-    Utils::Vector4d quat = {0., 0., 0., 0.};
+    Utils::Quaternion<double> quat = Utils::Quaternion<double>::zero();
 
     template <class Archive> void serialize(Archive &ar, long int) {
       ar &to_particle_id;
@@ -248,7 +250,7 @@ struct ParticlePosition {
 
 #ifdef ROTATION
   /** quaternion to define particle orientation */
-  Utils::Vector4d quat = {1., 0., 0., 0.};
+  Utils::Quaternion<double> quat = Utils::Quaternion<double>::identity();
   /** unit director calculated from the quaternion */
   Utils::Vector3d calc_director() const {
     return Utils::convert_quaternion_to_director(quat);

src/core/observables/CylindricalLBProfileObservable.hpp

@@ -55,7 +55,8 @@ class CylindricalLBProfileObservable : public CylindricalProfileObservable {
       // z-axis symmetry.
       std::tie(theta, rotation_axis) =
           Utils::rotation_params(Utils::Vector3d{{0.0, 0.0, 1.0}}, axis);
-      p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart);
+      if (theta > std::numeric_limits<double>::epsilon())
+        p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart);
       p = p_cart + center;
     }
   }

src/core/particle_data.cpp

@@ -147,7 +147,7 @@ using UpdatePropertyMessage = boost::variant
 using UpdatePositionMessage = boost::variant
         < UpdatePosition<Utils::Vector3d, &ParticlePosition::p>
 #ifdef ROTATION
-        , UpdatePosition<Utils::Vector4d, &ParticlePosition::quat>
+        , UpdatePosition<Utils::Quaternion<double>, &ParticlePosition::quat>
 #endif
         >;
 
@@ -750,7 +750,7 @@ void set_particle_dipm(int part, double dipm) {
 }
 
 void set_particle_dip(int part, double const *const dip) {
-  Utils::Vector4d quat;
+  Utils::Quaternion<double> quat;
   double dipm;
   std::tie(quat, dipm) =
       convert_dip_to_quat(Utils::Vector3d({dip[0], dip[1], dip[2]}));
@@ -769,7 +769,8 @@ void set_particle_virtual(int part, bool is_virtual) {
 #endif
 
 #ifdef VIRTUAL_SITES_RELATIVE
-void set_particle_vs_quat(int part, Utils::Vector4d const &vs_relative_quat) {
+void set_particle_vs_quat(int part,
+                          Utils::Quaternion<double> const &vs_relative_quat) {
   auto vs_relative = get_particle_data(part).p.vs_relative;
   vs_relative.quat = vs_relative_quat;
 
@@ -779,7 +780,7 @@ void set_particle_vs_quat(int part, Utils::Vector4d const &vs_relative_quat) {
 }
 
 void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance,
-                              Utils::Vector4d const &rel_ori) {
+                              Utils::Quaternion<double> const &rel_ori) {
   ParticleProperties::VirtualSitesRelativeParameters vs_relative;
   vs_relative.distance = vs_distance;
   vs_relative.to_particle_id = vs_relative_to;
@@ -837,9 +838,9 @@ void set_particle_mol_id(int part, int mid) {
 
 #ifdef ROTATION
 void set_particle_quat(int part, double *quat) {
-  mpi_update_particle<ParticlePosition, &Particle::r, Utils::Vector4d,
-                      &ParticlePosition::quat>(part,
-                                               Utils::Vector4d(quat, quat + 4));
+  mpi_update_particle<ParticlePosition, &Particle::r, Utils::Quaternion<double>,
+                      &ParticlePosition::quat>(
+      part, Utils::Quaternion<double>{quat[0], quat[1], quat[2], quat[3]});
 }
 
 void set_particle_omega_lab(int part, const Utils::Vector3d &omega_lab) {

src/core/particle_data.hpp

@@ -41,6 +41,7 @@
 
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
+#include <utils/quaternion.hpp>
 
 #include <cstddef>
 #include <memory>
@@ -235,9 +236,10 @@ void set_particle_dipm(int part, double dipm);
 void set_particle_virtual(int part, bool is_virtual);
 #endif
 #ifdef VIRTUAL_SITES_RELATIVE
-void set_particle_vs_quat(int part, Utils::Vector4d const &vs_relative_quat);
+void set_particle_vs_quat(int part,
+                          Utils::Quaternion<double> const &vs_relative_quat);
 void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance,
-                              Utils::Vector4d const &rel_ori);
+                              Utils::Quaternion<double> const &rel_ori);
 #endif
 
 #if defined(LANGEVIN_PER_PARTICLE) || defined(BROWNIAN_PER_PARTICLE)

src/core/rotation.cpp

@@ -54,8 +54,8 @@
  *                   Lagrange parameter lambda
  *  @param[out] Wd   Angular acceleration of the particle
  */
-static void define_Qdd(Particle const &p, Utils::Vector4d &Qd,
-                       Utils::Vector4d &Qdd, Utils::Vector3d &S,
+static void define_Qdd(Particle const &p, Utils::Quaternion<double> &Qd,
+                       Utils::Quaternion<double> &Qdd, Utils::Vector3d &S,
                        Utils::Vector3d &Wd) {
   /* calculate the first derivative of the quaternion */
   /* Eq. (4) @cite sonnenschein85a */
@@ -107,7 +107,7 @@ static void define_Qdd(Particle const &p, Utils::Vector4d &Qd,
       p.r.quat[3] * S1;
 
   S[0] = S1;
-  S[1] = Qd * Qdd;
+  S[1] = Utils::dot(Qd, Qdd);
   S[2] = Qdd.norm2();
 }
 
@@ -128,7 +128,7 @@ void propagate_omega_quat_particle(Particle &p) {
   if (p.p.rotation == ROTATION_FIXED)
     return;
 
-  Utils::Vector4d Qd{}, Qdd{};
+  Utils::Quaternion<double> Qd{}, Qdd{};
   Utils::Vector3d S{}, Wd{};
 
   // Clear rotational velocity for blocked rotation axes.

src/core/rotation.hpp

@@ -35,6 +35,7 @@
 #include <utils/mask.hpp>
 #include <utils/math/quaternion.hpp>
 #include <utils/math/rotation_matrix.hpp>
+#include <utils/quaternion.hpp>
 
 #include <cmath>
 #include <utility>
@@ -55,7 +56,7 @@ void convert_initial_torques(const ParticleRange &particles);
 // Frame conversion routines
 inline Utils::Vector3d
 convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec) {
-  return rotation_matrix(p.r.quat) * vec;
+  return p.r.quat * vec;
 }
 
 inline Utils::Vector3d convert_vector_space_to_body(const Particle &p,
@@ -89,7 +90,7 @@ auto convert_body_to_space(const Particle &p, const Utils::Matrix<T, 3, 3> &A) {
 #ifdef DIPOLES
 
 /** convert a dipole moment to quaternions and dipolar strength  */
-inline std::pair<Utils::Vector4d, double>
+inline std::pair<Utils::Quaternion<double>, double>
 convert_dip_to_quat(const Utils::Vector3d &dip) {
   auto quat = Utils::convert_director_to_quaternion(dip);
   return {quat, dip.norm()};
@@ -100,7 +101,7 @@ convert_dip_to_quat(const Utils::Vector3d &dip) {
 /** Rotate the particle p around the body-frame defined NORMALIZED axis
  *  @p aBodyFrame by amount @p phi.
  */
-inline Utils::Vector4d
+inline Utils::Quaternion<double>
 local_rotate_particle_body(Particle const &p,
                            const Utils::Vector3d &axis_body_frame,
                            const double phi) {
@@ -112,13 +113,8 @@ local_rotate_particle_body(Particle const &p,
 
   // Convert rotation axis to body-fixed frame
   axis = mask(p.p.rotation, axis).normalize();
-
-  auto const s = std::sin(phi / 2);
-  auto const q =
-      Utils::Vector4d{cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]}
-          .normalize();
-
-  return Utils::multiply_quaternions(p.r.quat, q);
+  Utils::Quaternion<double> q = boost::qvm::rot_quat(axis, phi);
+  return p.r.quat * q;
 }
 
 /** Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi

src/core/thermostats/brownian_inline.hpp

@@ -316,7 +316,7 @@ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian,
  *  @param[in]     p              %Particle
  *  @param[in]     dt             Time interval
  */
-inline Utils::Vector4d
+inline Utils::Quaternion<double>
 bd_drag_rot(Thermostat::GammaType const &brownian_gamma_rotation, Particle &p,
             double dt) {
   Thermostat::GammaType gamma;
@@ -395,8 +395,9 @@ bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation,
  *  @param[in]     p              %Particle
  *  @param[in]     dt             Time interval
  */
-inline Utils::Vector4d bd_random_walk_rot(BrownianThermostat const &brownian,
-                                          Particle const &p, double dt) {
+inline Utils::Quaternion<double>
+bd_random_walk_rot(BrownianThermostat const &brownian, Particle const &p,
+                   double dt) {
   // first, set defaults
   Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation;
 

src/core/unit_tests/random_test.hpp

@@ -23,6 +23,7 @@
 /* Helper functions to compute random numbers covariance in a single pass */
 
 #include <utils/Vector.hpp>
+#include <utils/quaternion.hpp>
 
 #include <boost/accumulators/accumulators.hpp>
 #include <boost/accumulators/statistics.hpp>
@@ -40,7 +41,8 @@
 #include <vector>
 
 namespace Utils {
-using VariantVectorXd = boost::variant<double, Vector2d, Vector3d, Vector4d>;
+using VariantVectorXd =
+    boost::variant<double, Vector2d, Vector3d, Vector4d, Quaternion<double>>;
 } // namespace Utils
 
 using Utils::VariantVectorXd;
@@ -54,6 +56,7 @@ class visitor_size : public boost::static_visitor<size_t> {
   template <size_t N> size_t operator()(Vector<double, N> const &v) const {
     return v.size();
   }
+  size_t operator()(Utils::Quaternion<double> const &q) const { return 4; }
   size_t operator()(double v) const { return 1; }
 };
 
@@ -63,6 +66,9 @@ class visitor_get : public boost::static_visitor<double> {
   double operator()(Vector<double, N> const &v, size_t i) const {
     return v[i];
   }
+  double operator()(Utils::Quaternion<double> const &q, size_t i) const {
+    return q[i];
+  }
   double operator()(double v, size_t i) const {
     assert(i == 0);
     return v;

src/core/virtual_sites.cpp

@@ -54,7 +54,7 @@ void set_virtual_sites(std::shared_ptr<VirtualSites> const &v) {
 #ifdef VIRTUAL_SITES_RELATIVE
 
 /** Calculate the rotation quaternion and distance between two particles */
-inline std::tuple<Utils::Vector4d, double>
+inline std::tuple<Utils::Quaternion<double>, double>
 calculate_vs_relate_to_params(Particle const &p_current,
                               Particle const &p_relate_to) {
   // get the distance between the particles
@@ -85,38 +85,38 @@ calculate_vs_relate_to_params(Particle const &p_current,
   //   = quat_(obtained from desired director)
   // Resolving this for the quat_(virtual particle)
 
-  Utils::Vector4d quat;
+  Utils::Quaternion<double> quat;
   // If the distance between real & virtual particle is 0 we just set the
   // relative orientation to {1 0 0 0}, as it is irrelevant but needs to be
   // a valid quaternion
   if (dist == 0) {
-    quat = {1, 0, 0, 0};
+    quat = Utils::Quaternion<double>::identity();
   } else {
     d.normalize();
 
     // Obtain quaternion from desired director
-    Utils::Vector4d quat_director = Utils::convert_director_to_quaternion(d);
+    Utils::Quaternion<double> quat_director =
+        Utils::convert_director_to_quaternion(d);
 
     // Define quaternion as described above
-    quat = {p_relate_to.r.quat * quat_director,
-            -quat_director[0] * p_relate_to.r.quat[1] +
-                quat_director[1] * p_relate_to.r.quat[0] +
-                quat_director[2] * p_relate_to.r.quat[3] -
-                quat_director[3] * p_relate_to.r.quat[2],
-            p_relate_to.r.quat[1] * quat_director[3] +
-                p_relate_to.r.quat[0] * quat_director[2] -
-                p_relate_to.r.quat[3] * quat_director[1] -
-                p_relate_to.r.quat[2] * quat_director[0],
-            quat_director[3] * p_relate_to.r.quat[0] -
-                p_relate_to.r.quat[3] * quat_director[0] +
-                p_relate_to.r.quat[2] * quat_director[1] -
-                p_relate_to.r.quat[1] * quat_director[2]};
-    auto const norm = p_relate_to.r.quat * p_relate_to.r.quat;
-    quat /= norm;
+    quat =
+        Utils::Quaternion<double>{Utils::dot(p_relate_to.r.quat, quat_director),
+                                  -quat_director[0] * p_relate_to.r.quat[1] +
+                                      quat_director[1] * p_relate_to.r.quat[0] +
+                                      quat_director[2] * p_relate_to.r.quat[3] -
+                                      quat_director[3] * p_relate_to.r.quat[2],
+                                  p_relate_to.r.quat[1] * quat_director[3] +
+                                      p_relate_to.r.quat[0] * quat_director[2] -
+                                      p_relate_to.r.quat[3] * quat_director[1] -
+                                      p_relate_to.r.quat[2] * quat_director[0],
+                                  quat_director[3] * p_relate_to.r.quat[0] -
+                                      p_relate_to.r.quat[3] * quat_director[0] +
+                                      p_relate_to.r.quat[2] * quat_director[1] -
+                                      p_relate_to.r.quat[1] * quat_director[2]};
+    quat /= p_relate_to.r.quat.norm2();
 
     // Verify result
-    Utils::Vector4d qtemp =
-        Utils::multiply_quaternions(p_relate_to.r.quat, quat);
+    Utils::Quaternion<double> qtemp = p_relate_to.r.quat * quat;
     for (int i = 0; i < 4; i++)
       if (fabs(qtemp[i] - quat_director[i]) > 1E-9)
         fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n", i,
@@ -140,7 +140,7 @@ void vs_relate_to(int part_num, int relate_to) {
   auto const &p_current = get_particle_data(part_num);
   auto const &p_relate_to = get_particle_data(relate_to);
 
-  Utils::Vector4d quat;
+  Utils::Quaternion<double> quat;
   double dist;
   std::tie(quat, dist) = calculate_vs_relate_to_params(p_current, p_relate_to);
 

src/core/virtual_sites/VirtualSitesRelative.cpp

@@ -32,6 +32,7 @@
 #include <utils/math/quaternion.hpp>
 #include <utils/math/sqr.hpp>
 #include <utils/math/tensor_product.hpp>
+#include <utils/quaternion.hpp>
 
 #include <stdexcept>
 
@@ -42,10 +43,10 @@ namespace {
  * @param vs_rel Parameters for the virtual site.
  * @return Orientation quaternion of the virtual site.
  */
-Utils::Vector4d
+Utils::Quaternion<double>
 orientation(Particle const *p_ref,
             const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) {
-  return multiply_quaternions(p_ref->r.quat, vs_rel.quat);
+  return p_ref->r.quat * vs_rel.quat;
 }
 
 /**
@@ -62,10 +63,9 @@ Utils::Vector3d connection_vector(
   // This is obtained, by multiplying the quaternion representing the director
   // of the real particle with the quaternion of the virtual particle, which
   // specifies the relative orientation.
-  auto const director =
-      Utils::convert_quaternion_to_director(
-          Utils::multiply_quaternions(p_ref->r.quat, vs_rel.rel_orientation))
-          .normalize();
+  auto const director = Utils::convert_quaternion_to_director(
+                            p_ref->r.quat * vs_rel.rel_orientation)
+                            .normalize();
 
   return vs_rel.distance * director;
 }

src/python/espressomd/particle_data.pxd

@@ -17,7 +17,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Here we create something to handle particles
-from .utils cimport Vector4d, Vector3d, Vector3i, Span
+from .utils cimport Vector4d, Vector3d, Vector3i, Span, Quaternion
 from libcpp cimport bool
 from libcpp.vector cimport vector  # import std::vector as vector
 from libc cimport stdint
@@ -151,8 +151,8 @@ cdef extern from "particle_data.hpp":
     IF VIRTUAL_SITES_RELATIVE:
         void pointer_to_vs_relative(const particle * P, const int * & res1, const double * & res2, const double * & res3)
         void pointer_to_vs_quat(const particle * P, const double * & res)
-        void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance, Vector4d rel_ori)
-        void set_particle_vs_quat(int part, Vector4d vs_quat)
+        void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance, Quaternion[double] rel_ori)
+        void set_particle_vs_quat(int part, Quaternion[double] vs_quat)
 
     void pointer_to_q(const particle * P, const double * & res)
 

src/python/espressomd/particle_data.pyx

@@ -634,7 +634,7 @@ cdef class ParticleHandle:
             def __set__(self, q):
                 check_type_or_throw_except(
                     q, 4, float, "vs_quat has to be an array-like of length 4")
-                cdef Vector4d _q
+                cdef Quaternion[double] _q
                 for i in range(4):
                     _q[i] = q[i]
                 set_particle_vs_quat(self._id, _q)
@@ -668,7 +668,7 @@ cdef class ParticleHandle:
                 _relto = x[0]
                 _dist = x[1]
                 q = x[2]
-                cdef Vector4d _q
+                cdef Quaternion[double] _q
                 for i in range(4):
                     _q[i] = q[i]