Separate quaternion algebra from particle rotation

doc/sphinx/particles.rst

@@ -365,8 +365,10 @@ site is placed at a fixed distance from the non-virtual particle. When
 the non-virtual particle rotates, the virtual sites rotates on an orbit
 around the non-virtual particles center.
 
-To use this implementation of virtual sites, activate the feature ``VIRTUAL_SITES_RELATIVE``. Furthermore, an instance of :class:`espressomd.virtual_sites.VirtualSitesRelative` has to be set as the active virtual sites scheme (see above).
-To set up a virtual site,
+To use this implementation of virtual sites, activate the feature
+``VIRTUAL_SITES_RELATIVE``. Furthermore, an instance of
+:class:`espressomd.virtual_sites.VirtualSitesRelative` has to be set as the
+active virtual sites scheme (see above). To set up a virtual site,
 
 #. Place the particle to which the virtual site should be related. It
    needs to be in the center of mass of the rigid arrangement of
@@ -378,7 +380,9 @@ To set up a virtual site,
        p = system.part.add(pos=(1, 2, 3))
        p.vs_auto_relate_to(<ID>)
 
-   where <ID> is the id of the central particle. This will also set the :attr:`espressomd.particle_data.ParticleHandle.virtual` attribute on the particle to 1.
+   where <ID> is the id of the central particle. This will also set the
+   :attr:`espressomd.particle_data.ParticleHandle.virtual` attribute on
+   the particle to ``True``.
 
 #. Repeat the previous step with more virtual sites, if desired.
 
@@ -407,8 +411,8 @@ Please note:
 -  In a simulation on more than one CPU, the effective cell size needs
    to be larger than the largest distance between a non-virtual particle
    and its associated virtual sites. To this aim, when running on more than one core,
-   you need to set the
-   system's :attr:`espressomd.system.System.min_global_cut` attribute to this largest distance.
+   you need to set the system's :attr:`espressomd.system.System.min_global_cut`
+   attribute to this largest distance.
    An error is generated when this requirement is not met.
 
 -  If the virtual sites represent actual particles carrying a mass, the
@@ -425,9 +429,13 @@ Inertialess lattice Boltzmann tracers
 
 :class:`espressomd.virtual_sites.VirtualSitesInertialessTracers`
 
-When this implementation is selected, the virtual sites follow the motion of a lattice Boltzmann fluid (both, CPU and GPU). This is achieved by integrating their position using the fluid velocity at the virtual sites' position.
-Forces acting on the virtual sites are directly transferred as force density onto the lattice Boltzmann fluid, making the coupling free of inertia.
-The feature stems from the implementation of the :ref:`Immersed Boundary Method for soft elastic objects`, but can be used independently.
+When this implementation is selected, the virtual sites follow the motion of a
+lattice Boltzmann fluid (both, CPU and GPU). This is achieved by integrating
+their position using the fluid velocity at the virtual sites' position.
+Forces acting on the virtual sites are directly transferred as force density
+onto the lattice Boltzmann fluid, making the coupling free of inertia.
+The feature stems from the implementation of the
+:ref:`Immersed Boundary Method for soft elastic objects`, but can be used independently.
 
 For correct results, the LB thermostat has to be deactivated for virtual sites::
 

src/core/collision.cpp

@@ -368,7 +368,7 @@ void place_vs_and_relate_to_particle(const int current_vs_pid,
   new_part.r.p = pos;
   auto p_vs = append_indexed_particle(local_cells.cell[0], std::move(new_part));
 
-  local_vs_relate_to(p_vs, &get_part(relate_to));
+  local_vs_relate_to(*p_vs, get_part(relate_to));
 
   p_vs->p.is_virtual = true;
   p_vs->p.type = collision_params.vs_particle_type;

src/core/particle_data.cpp

@@ -153,7 +153,7 @@ using UpdatePropertyMessage = boost::variant
 #ifdef VIRTUAL_SITES
         , UpdateProperty<bool, &Prop::is_virtual>
 #ifdef VIRTUAL_SITES_RELATIVE
-        , UpdateProperty<ParticleProperties::VirtualSitesRelativeParameteres,
+        , UpdateProperty<ParticleProperties::VirtualSitesRelativeParameters,
                          &Prop::vs_relative>
 #endif
 #endif
@@ -908,24 +908,24 @@ void set_particle_virtual(int part, bool is_virtual) {
 #endif
 
 #ifdef VIRTUAL_SITES_RELATIVE
-void set_particle_vs_quat(int part, double *vs_relative_quat) {
+void set_particle_vs_quat(int part, Utils::Vector4d const &vs_relative_quat) {
   auto vs_relative = get_particle_data(part).p.vs_relative;
-  vs_relative.quat = Utils::Vector4d(vs_relative_quat, vs_relative_quat + 4);
+  vs_relative.quat = vs_relative_quat;
 
   mpi_update_particle_property<
-      ParticleProperties::VirtualSitesRelativeParameteres,
+      ParticleProperties::VirtualSitesRelativeParameters,
       &ParticleProperties::vs_relative>(part, vs_relative);
 }
 
 void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance,
-                              double *rel_ori) {
-  ParticleProperties::VirtualSitesRelativeParameteres vs_relative;
+                              Utils::Vector4d const &rel_ori) {
+  ParticleProperties::VirtualSitesRelativeParameters vs_relative;
   vs_relative.distance = vs_distance;
   vs_relative.to_particle_id = vs_relative_to;
-  vs_relative.rel_orientation = {rel_ori, rel_ori + 4};
+  vs_relative.rel_orientation = rel_ori;
 
   mpi_update_particle_property<
-      ParticleProperties::VirtualSitesRelativeParameteres,
+      ParticleProperties::VirtualSitesRelativeParameters,
       &ParticleProperties::vs_relative>(part, vs_relative);
 }
 #endif

src/core/particle_data.hpp

@@ -39,6 +39,7 @@
 #include <utils/List.hpp>
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
+#include <utils/math/quaternion.hpp>
 
 #include <memory>
 
@@ -93,8 +94,8 @@ struct ParticleProperties {
   /** particle type, used for non bonded interactions. */
   int type = 0;
 
-#ifdef MASS
   /** particle mass */
+#ifdef MASS
   double mass = 1.0;
 #else
   constexpr static double mass{1.0};
@@ -112,8 +113,7 @@ struct ParticleProperties {
   Utils::Vector3d out_direction = {0., 0., 0.};
 #endif
 
-  // Determines, whether a particle's rotational degrees of freedom are
-  // integrated
+  /** bitfield for the particle axes of rotation */
   int rotation = 0;
 
   /** charge. */
@@ -129,26 +129,26 @@ struct ParticleProperties {
 #endif
 
 #ifdef DIPOLES
-  /** dipole moment (absolute value)*/
+  /** dipole moment (absolute value) */
   double dipm = 0.;
 #endif
 
 #ifdef VIRTUAL_SITES
   /** is particle virtual */
   bool is_virtual = false;
 #ifdef VIRTUAL_SITES_RELATIVE
-  /** In case, the "relative" implementation of virtual sites is enabled, the
-  following properties define, with respect to which real particle a virtual
-  site is placed and in what distance. The relative orientation of the vector
-  pointing from real particle to virtual site with respect to the orientation
-  of the real particle is stored in the virtual site's quaternion attribute.
-  */
-  struct VirtualSitesRelativeParameteres {
+  /** The following properties define, with respect to which real particle a
+   *  virtual site is placed and at what distance. The relative orientation of
+   *  the vector pointing from real particle to virtual site with respect to the
+   *  orientation of the real particle is stored in the virtual site's
+   *  quaternion attribute.
+   */
+  struct VirtualSitesRelativeParameters {
     int to_particle_id = 0;
     double distance = 0;
-    // Store relative position of the virtual site.
+    /** Relative position of the virtual site. */
     Utils::Vector4d rel_orientation = {0., 0., 0., 0.};
-    // Store the orientation of the virtual particle in the body fixed frame.
+    /** Orientation of the virtual particle in the body fixed frame. */
     Utils::Vector4d quat = {0., 0., 0., 0.};
 
     template <class Archive> void serialize(Archive &ar, long int) {
@@ -158,7 +158,6 @@ struct ParticleProperties {
       ar &quat;
     }
   } vs_relative;
-
 #endif
 #else  /* VIRTUAL_SITES */
   static constexpr const bool is_virtual = false;
@@ -171,7 +170,7 @@ struct ParticleProperties {
 #else
   Utils::Vector3d gamma = {-1., -1., -1.};
 #endif // PARTICLE_ANISOTROPY
-/* Friction coefficient gamma for rotation */
+/** Friction coefficient gamma for rotation */
 #ifdef ROTATION
 #ifndef PARTICLE_ANISOTROPY
   double gamma_rot = -1.;
@@ -202,31 +201,30 @@ struct ParticleProperties {
 };
 
 /** Positional information on a particle. Information that is
-    communicated to calculate interactions with ghost particles. */
+ *  communicated to calculate interactions with ghost particles.
+ */
 struct ParticlePosition {
   /** periodically folded position. */
   Utils::Vector3d p = {0, 0, 0};
 
 #ifdef ROTATION
-  /** quaternions to define particle orientation */
+  /** quaternion to define particle orientation */
   Utils::Vector4d quat = {1., 0., 0., 0.};
-  /** unit director calculated from the quaternions */
+  /** unit director calculated from the quaternion */
   Utils::Vector3d calc_director() const {
-    return {2 * (quat[1] * quat[3] + quat[0] * quat[2]),
-            2 * (quat[2] * quat[3] - quat[0] * quat[1]),
-            quat[0] * quat[0] - quat[1] * quat[1] - quat[2] * quat[2] +
-                quat[3] * quat[3]};
+    return Utils::convert_quaternion_to_director(quat);
   };
 #endif
 
 #ifdef BOND_CONSTRAINT
-  /**stores the particle position at the previous time step*/
+  /** particle position at the previous time step */
   Utils::Vector3d p_old = {0., 0., 0.};
 #endif
 };
 
 /** Force information on a particle. Forces of ghost particles are
-    collected and added up to the force of the original particle. */
+ *  collected and added up to the force of the original particle.
+ */
 struct ParticleForce {
   ParticleForce() = default;
   ParticleForce(ParticleForce const &) = default;
@@ -676,9 +674,9 @@ 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, double *vs_relative_quat);
+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,
-                              double *rel_ori);
+                              Utils::Vector4d const &rel_ori);
 #endif
 
 #ifdef LANGEVIN_PER_PARTICLE

src/core/rotation.cpp

@@ -46,6 +46,7 @@
 #include "thermostat.hpp"
 
 #include <utils/constants.hpp>
+#include <utils/math/quaternion.hpp>
 #include <utils/math/rotation_matrix.hpp>
 
 #include <cmath>
@@ -62,51 +63,6 @@
 static void define_Qdd(Particle const &p, double Qd[4], double Qdd[4],
                        double S[3], double Wd[3]);
 
-/** Convert director to quaternions */
-int convert_director_to_quat(const Utils::Vector3d &d, Utils::Vector4d &quat) {
-  double theta2, phi2;
-
-  // Calculate magnitude of the given vector
-  auto const dm = d.norm();
-
-  // The vector needs to be != 0 to be converted into a quaternion
-  if (dm < ROUND_ERROR_PREC) {
-    return 1;
-  }
-  // Calculate angles
-  auto const d_xy = sqrt(d[0] * d[0] + d[1] * d[1]);
-  // If dipole points along z axis:
-  if (d_xy == 0) {
-    // We need to distinguish between (0,0,d_z) and (0,0,d_z)
-    if (d[2] > 0)
-      theta2 = 0;
-    else
-      theta2 = Utils::pi() / 2.;
-    phi2 = 0;
-  } else {
-    // Here, we take care of all other directions
-    // Here we suppose that theta2 = 0.5*theta and phi2 = 0.5*(phi -
-    // Utils::pi()/2), where theta and phi - angles are in spherical coordinates
-    theta2 = 0.5 * acos(d[2] / dm);
-    if (d[1] < 0)
-      phi2 = -0.5 * acos(d[0] / d_xy) - Utils::pi() * 0.25;
-    else
-      phi2 = 0.5 * acos(d[0] / d_xy) - Utils::pi() * 0.25;
-  }
-
-  // Calculate the quaternion from the angles
-  auto const cos_theta2 = cos(theta2);
-  auto const sin_theta2 = sin(theta2);
-  auto const cos_phi2 = cos(phi2);
-  auto const sin_phi2 = sin(phi2);
-  quat[0] = cos_theta2 * cos_phi2;
-  quat[1] = -sin_theta2 * cos_phi2;
-  quat[2] = -sin_theta2 * sin_phi2;
-  quat[3] = cos_theta2 * sin_phi2;
-
-  return 0;
-}
-
 void define_Qdd(Particle const &p, double Qd[4], double Qdd[4], double S[3],
                 double Wd[3]) {
   /* calculate the first derivative of the quaternion */
@@ -202,10 +158,7 @@ void propagate_omega_quat_particle(Particle &p) {
                    (S[0] + time_step * (S[1] + time_step_half / 2. *
                                                    (S[2] - S[0] * S[0]))));
 
-  for (int j = 0; j < 3; j++) {
-    p.m.omega[j] += time_step_half * Wd[j];
-  }
-
+  p.m.omega += time_step_half * Wd;
   p.r.quat += time_step * (Qd + time_step_half * Qdd) - lambda * p.r.quat;
 
   /* and rescale quaternion, so it is exactly of unit length */
@@ -363,21 +316,12 @@ void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame,
 
   axis /= l;
 
-  double q[4];
-  q[0] = cos(phi / 2);
-  double tmp = sin(phi / 2);
-  q[1] = tmp * axis[0];
-  q[2] = tmp * axis[1];
-  q[3] = tmp * axis[2];
-
-  // Normalize
-  normalize_quaternion(q);
+  double s = sin(phi / 2);
+  Utils::Vector4d q = {cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]};
+  q.normalize();
 
   // Rotate the particle
-  double qn[4]; // Resulting quaternion
-  multiply_quaternions(p.r.quat, q, qn);
-  for (int k = 0; k < 4; k++)
-    p.r.quat[k] = qn[k];
+  p.r.quat = Utils::multiply_quaternions(p.r.quat, q);
 }
 
-#endif
+#endif // ROTATION