Separate quaternion algebra from particle rotation

doc/sphinx/particles.rst

@@ -369,8 +369,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
@@ -382,21 +384,23 @@ 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.
 
-#. To update the positions of all virtual sites, call
+#. To update the positions of all virtual sites, call::
 
-   system.integrator.run(0,recalc_forces=True)
+       system.integrator.run(0, recalc_forces=True)
 
 Please note:
 
 -  The relative position of the virtual site is defined by its distance
    from the non-virtual particle, the id of the non-virtual particle and
    a quaternion which defines the vector from non-virtual particle to
    virtual site in the non-virtual particles body-fixed frame. This
-   information is saved in the virtual site's `espressomd.particle_data.ParticleHandle.vs_relative` attribute.
+   information is saved in the virtual site's :attr:`espressomd.particle_data.ParticleHandle.vs_relative` attribute.
    Take care, not to overwrite it after using ``vs_auto_relate``.
 
 -  Virtual sites can not be placed relative to other virtual sites, as
@@ -405,13 +409,14 @@ Please note:
    placed in the center of mass of the rigid arrangement of particles.
 
 -  In case you know the correct quaternions, you can also setup a
-   virtual site using its :attr:`espressomd.particle_data.ParticleHandle.vs_relative` and :attr:`espressomd.particle_data.ParticleHandle.virtual` attributes.
+   virtual site using its :attr:`espressomd.particle_data.ParticleHandle.vs_relative`
+   and :attr:`espressomd.particle_data.ParticleHandle.virtual` attributes.
 
 -  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
@@ -428,9 +433,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

@@ -156,7 +156,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
@@ -924,24 +924,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

@@ -40,6 +40,8 @@
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
 
+#include "quaternion.hpp"
+
 #include <memory>
 
 /************************************************
@@ -96,8 +98,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};
@@ -120,8 +122,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. */
@@ -137,26 +138,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) {
@@ -166,7 +167,6 @@ struct ParticleProperties {
       ar &quat;
     }
   } vs_relative;
-
 #endif
 #else  /* VIRTUAL_SITES */
   static constexpr const bool is_virtual = false;
@@ -179,7 +179,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.;
@@ -210,31 +210,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 */
   inline const 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 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;
@@ -730,9 +729,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/quaternion.hpp

@@ -0,0 +1,94 @@
+/*
+  Copyright (C) 2010-2018 The ESPResSo project
+  Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+    Max-Planck-Institute for Polymer Research, Theory Group
+
+  This file is part of ESPResSo.
+
+  ESPResSo is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  ESPResSo is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef QUATERNION_H
+#define QUATERNION_H
+/** \file
+ *  Quaternion algebra.
+ */
+
+#ifdef ROTATION
+
+#include <utils/Vector.hpp>
+#include <utils/constants.hpp>
+
+/** Multiply two quaternions */
+template <class T>
+Utils::Vector<T, 4> multiply_quaternions(Utils::Vector<T, 4> const &a,
+                                         Utils::Vector<T, 4> const &b) {
+  // Formula from http://www.j3d.org/matrix_faq/matrfaq_latest.html
+  return {a[0] * b[0] - a[1] * b[1] - a[2] * b[2] - a[3] * b[3],
+          a[0] * b[1] + a[1] * b[0] + a[2] * b[3] - a[3] * b[2],
+          a[0] * b[2] + a[2] * b[0] + a[3] * b[1] - a[1] * b[3],
+          a[0] * b[3] + a[3] * b[0] + a[1] * b[2] - a[2] * b[1]};
+}
+
+/** Convert quaternion to director */
+template <class T>
+Utils::Vector<T, 3> quaternion_to_director(Utils::Vector<T, 4> const &quat) {
+  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]};
+}
+
+/** Convert director to quaternion
+ *  @param d  Director
+ *  @return A quaternion from the director, or {1, 0, 0, 0}
+ *  if the director is the null vector.
+ */
+template <class T>
+Utils::Vector<T, 4>
+convert_director_to_quaternion(Utils::Vector<T, 3> const &d) {
+
+  auto const dm = d.norm();
+
+  // null vectors cannot be converted to quaternions
+  if (dm < ROUND_ERROR_PREC) {
+    return {1, 0, 0, 0};
+  }
+
+  // Calculate angles
+  auto const d_xy = sqrt(d[0] * d[0] + d[1] * d[1]);
+  double theta2, phi2;
+  if (d_xy == 0) {
+    // Here the director is co-linear with the z-azis
+    // We need to distinguish between (0, 0, +d_z) and (0, 0, -d_z)
+    theta2 = (d[2] > 0) ? 0 : Utils::pi() / 2;
+    phi2 = 0;
+  } else {
+    // Here we take care of all other directions
+    // We suppose that theta2 = theta/2 and phi2 = (phi - pi/2)/2,
+    // where angles theta and phi are in spherical coordinates
+    theta2 = acos(d[2] / dm) / 2;
+    phi2 = ((d[1] > 0) ? 1 : -1) * acos(d[0] / d_xy) / 2 - Utils::pi() / 4;
+  }
+
+  // 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);
+  return {cos_theta2 * cos_phi2, -sin_theta2 * cos_phi2, -sin_theta2 * sin_phi2,
+          cos_theta2 * sin_phi2};
+}
+
+#endif // ROTATION
+#endif