Merge #3164 #3245

3164: Separate quaternion algebra from particle rotation r=fweik a=jngrad Description of changes: - extract quaternion algebra from the `particle_data.hpp` and `rotation.hpp` files - removes code duplication caused by a (now resolved) circular dependency (see #3157) - makes it possible to replace the quaternion code by a dedicated library, e.g. [boost:qvm::quat](https://www.boost.org/doc/libs/1_68_0/libs/qvm/doc/index.html) or [boost::math::quaternion](https://www.boost.org/doc/libs/1_62_0/libs/math/doc/html/quaternions.html) in the core (see #2289), [rowan](https://rowan.readthedocs.io/en/latest/) in the interface (see #2964) - simplify code around call sites using Vector4d arithmetic, `std::tuple`, Particle references - documentation cleanup 3245: Fix thermostat and integrator checkpointing r=KaiSzuttor a=jngrad The checkpointing mechanism silently broke in 4.1 for the SD and NPT integrators and LB and NPT thermostats. This was fixed, and now all integrators and thermostats checkpoints are tested in CI. Co-authored-by: Jean-Noël Grad <[email protected]> Co-authored-by: Kai Szuttor <[email protected]>
espressomd · Oct 15, 2019 · 054183b · 054183b
3 parents 0b7b357 + 4b92283 + fb3c8fe
commit 054183b
Show file tree

Hide file tree

Showing 23 changed files with 479 additions and 322 deletions.
diff --git a/doc/sphinx/particles.rst b/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::
 

diff --git a/src/core/collision.cpp b/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;

diff --git a/src/core/npt.cpp b/src/core/npt.cpp
@@ -16,20 +16,27 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#include "npt.hpp"
+#include <cmath>
+
 #include "communication.hpp"
 #include "config.hpp"
 #include "errorhandling.hpp"
 #include "integrate.hpp"
+#include "npt.hpp"
 
 #ifdef NPT
 void synchronize_npt_state(int n_steps) {
   nptiso.invalidate_p_vel = false;
   MPI_Bcast(&nptiso.p_inst, 1, MPI_DOUBLE, 0, comm_cart);
   MPI_Bcast(&nptiso.p_diff, 1, MPI_DOUBLE, 0, comm_cart);
   MPI_Bcast(&nptiso.volume, 1, MPI_DOUBLE, 0, comm_cart);
-  if (this_node == 0)
-    nptiso.p_inst_av /= 1.0 * n_steps;
+  if (this_node == 0) {
+    if (n_steps == 0) {
+      nptiso.p_inst_av = std::nan("");
+    } else {
+      nptiso.p_inst_av /= 1.0 * n_steps;
+    }
+  }
   MPI_Bcast(&nptiso.p_inst_av, 1, MPI_DOUBLE, 0, comm_cart);
 }
 

diff --git a/src/core/particle_data.cpp b/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

diff --git a/src/core/particle_data.hpp b/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

diff --git a/src/core/rotation.cpp b/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