diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib
index cf12c2459e4..967580da4d2 100644
--- a/doc/doxygen/bibliography.bib
+++ b/doc/doxygen/bibliography.bib
@@ -443,16 +443,37 @@ @Article{sonnenschein85a
}
@book{allen2017,
- title={Computer simulation of liquids},
+ title={{Computer simulation of liquids}},
author={Allen, Michael P and Tildesley, Dominic J},
year={2017},
publisher={Oxford University Press},
- url={https://global.oup.com/academic/product/computer-simulation-of-liquids-9780198803201},
doi={10.1093/oso/9780198803195.001.0001},
isbn={9780198803195},
edition={2nd},
}
+@article{martys99a,
+ title = {{Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions}},
+ author = {Martys, Nicos S. and Mountain, Raymond D.},
+ journal = {Physical Review E},
+ volume = {59},
+ number = {3},
+ pages = {3733--3736},
+ year = {1999},
+ doi = {10.1103/PhysRevE.59.3733},
+}
+
+@article{omelyan98a,
+author = {Omelyan, Igor P.},
+title = {{On the numerical integration of motion for rigid polyatomics: The modified quaternion approach}},
+journal = {Computers in Physics},
+volume = {12},
+number = {1},
+pages = {97-103},
+year = {1998},
+doi = {10.1063/1.168642},
+}
+
@Article{swope92a,
author = {Swope, William C. and Ferguson, David M.},
title = {{Alternative expressions for energies and forces due to angle bending and torsional energy}},
diff --git a/doc/sphinx/system_manipulation.rst b/doc/sphinx/system_manipulation.rst
index d7eba327dd5..4188de7e1f1 100644
--- a/doc/sphinx/system_manipulation.rst
+++ b/doc/sphinx/system_manipulation.rst
@@ -80,8 +80,13 @@ Force capping can be activated via :class:`espressomd.system.System`::
This command will limit the magnitude of the force to :math:`r F_\mathrm{max}`.
Energies are not affected by the capping, so the energy can be used to
-identify the remaining overlap. The force capping is switched off by setting
-:math:`F_\mathrm{max}=0`.
+identify the remaining overlap. Torques are also not affected by the
+capping. Force capping is switched off by setting :math:`F_\mathrm{max}=0`.
+
+For simple systems, it is often more convenient to use the
+:ref:`Steepest descent` algorithm instead of writing a tailored warmup
+loop in Python. The steepest descent algorithm will integrate the system
+while capping both the maximum displacement and maximum rotation.
.. _Galilei Transform and Particle Velocity Manipulation:
@@ -96,14 +101,14 @@ in affecting the velocity of the system. ::
* Particle motion and rotation
-::
+ ::
gt.kill_particle_motion()
-This command halts all particles in the current simulation, setting
-their velocities to zero, as well as their angular momentum if the
-option ``rotation`` is specified and the feature ``ROTATION`` has been
-compiled in.
+ This command halts all particles in the current simulation, setting
+ their velocities to zero, as well as their angular momentum if the
+ option ``rotation`` is specified and the feature ``ROTATION`` has been
+ compiled in.
* Forces and torques acting on the particles
diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp
index 2b5ec9b045f..098fbe80cf0 100644
--- a/src/core/rotation.cpp
+++ b/src/core/rotation.cpp
@@ -21,16 +21,14 @@
/** \file
* Molecular dynamics integrator for rotational motion.
*
- * A velocity Verlet algorithm
- * using quaternions is implemented to tackle rotational motion.
- * See @cite allen2017 for the quaternion components indexing used here.
+ * A velocity Verlet algorithm using quaternions is implemented to tackle
+ * rotational motion. See @cite martys99a for the method and
+ * @cite allen2017 for the quaternion components indexing used here.
* A random torque and a friction
* term are added to provide the constant NVT conditions. Due to this feature
- * all particles are
+ * all particles are
* treated as 3D objects with 3 translational and 3 rotational degrees of
- * freedom if ROTATION
- * flag is set in \ref config.hpp "config.hpp".
+ * freedom if ROTATION is compiled in.
*/
#include "rotation.hpp"
@@ -44,9 +42,11 @@
#include
-/** Calculate the derivatives of the quaternion and angular acceleration
- * for a given particle.
- * See @cite sonnenschein85a
+/** @brief Calculate the derivatives of the quaternion and angular
+ * acceleration for a given particle.
+ * See @cite sonnenschein85a. Please note that ESPResSo uses scalar-first
+ * notation for quaternions, while @cite sonnenschein85a use scalar-last
+ * notation.
* @param[in] p %Particle
* @param[out] Qd First derivative of the particle quaternion
* @param[out] Qdd Second derivative of the particle quaternion
@@ -54,8 +54,9 @@
* Lagrange parameter lambda
* @param[out] Wd Angular acceleration of the particle
*/
-static void define_Qdd(Particle const &p, double Qd[4], double Qdd[4],
- double S[3], double Wd[3]) {
+static void define_Qdd(Particle const &p, Utils::Vector4d &Qd,
+ Utils::Vector4d &Qdd, Utils::Vector3d &S,
+ Utils::Vector3d &Wd) {
/* calculate the first derivative of the quaternion */
/* Eq. (4) @cite sonnenschein85a */
Qd[0] = 0.5 * (-p.r.quat[1] * p.m.omega[0] - p.r.quat[2] * p.m.omega[1] -
@@ -76,22 +77,16 @@ static void define_Qdd(Particle const &p, double Qd[4], double Qdd[4],
Wd[0] = (p.f.torque[0] + p.m.omega[1] * p.m.omega[2] *
(p.p.rinertia[1] - p.p.rinertia[2])) /
p.p.rinertia[0];
- else
- Wd[0] = 0.0;
if (p.p.rotation & ROTATION_Y)
Wd[1] = (p.f.torque[1] + p.m.omega[2] * p.m.omega[0] *
(p.p.rinertia[2] - p.p.rinertia[0])) /
p.p.rinertia[1];
- else
- Wd[1] = 0.0;
if (p.p.rotation & ROTATION_Z)
Wd[2] = (p.f.torque[2] + p.m.omega[0] * p.m.omega[1] *
(p.p.rinertia[0] - p.p.rinertia[1])) /
p.p.rinertia[2];
- else
- Wd[2] = 0.0;
- auto const S1 = Qd[0] * Qd[0] + Qd[1] * Qd[1] + Qd[2] * Qd[2] + Qd[3] * Qd[3];
+ auto const S1 = Qd.norm2();
/* Calculate the second derivative of the quaternion. */
/* Eq. (8) @cite sonnenschein85a */
@@ -112,12 +107,20 @@ static void define_Qdd(Particle const &p, double Qd[4], double Qdd[4],
p.r.quat[3] * S1;
S[0] = S1;
- S[1] = Qd[0] * Qdd[0] + Qd[1] * Qdd[1] + Qd[2] * Qdd[2] + Qd[3] * Qdd[3];
- S[2] = Qdd[0] * Qdd[0] + Qdd[1] * Qdd[1] + Qdd[2] * Qdd[2] + Qdd[3] * Qdd[3];
+ S[1] = Qd * Qdd;
+ S[2] = Qdd.norm2();
}
-/** propagate angular velocities and quaternions
- * \todo implement for fixed_coord_flag
+/**
+ * See @cite omelyan98a. Please note that ESPResSo uses scalar-first
+ * notation for quaternions, while @cite omelyan98a use scalar-last
+ * notation.
+ *
+ * For very high angular velocities (e.g. if the product of @ref time_step
+ * with the largest component of @ref ParticleMomentum::omega "p.m.omega"
+ * is superior to ~2.0), the calculation might fail.
+ *
+ * \todo implement for fixed_coord_flag
*/
void propagate_omega_quat_particle(Particle &p) {
@@ -131,14 +134,15 @@ void propagate_omega_quat_particle(Particle &p) {
// Clear rotational velocity for blocked rotation axes.
p.m.omega = Utils::mask(p.p.rotation, p.m.omega);
- define_Qdd(p, Qd.data(), Qdd.data(), S.data(), Wd.data());
+ define_Qdd(p, Qd, Qdd, S, Wd);
- /* Eq. (12) @cite sonnenschein85a. */
- auto const lambda =
- 1 - S[0] * time_step_squared_half -
- sqrt(1 - time_step_squared *
- (S[0] + time_step * (S[1] + time_step_half / 2. *
- (S[2] - S[0] * S[0]))));
+ /* Eq. (12) @cite omelyan98a. */
+ auto const square =
+ 1 - time_step_squared *
+ (S[0] +
+ time_step * (S[1] + time_step_half / 2. * (S[2] - S[0] * S[0])));
+ assert(square >= 0.);
+ auto const lambda = 1 - S[0] * time_step_squared_half - sqrt(square);
p.m.omega += time_step_half * Wd;
p.r.quat += time_step * (Qd + time_step_half * Qdd) - lambda * p.r.quat;
@@ -152,8 +156,6 @@ void propagate_omega_quat_particle(Particle &p) {
}
}
-/** convert the torques to the body-fixed frames and propagate angular
- * velocities */
void convert_torques_propagate_omega(const ParticleRange &particles) {
for (auto &p : particles) {
// Skip particle if rotation is turned off entirely for it.
@@ -171,7 +173,7 @@ void convert_torques_propagate_omega(const ParticleRange &particles) {
/* if the tensor of inertia is isotropic, the following refinement is not
needed.
Otherwise repeat this loop 2-3 times depending on the required accuracy
- */
+ */
const double rinertia_diff_01 = p.p.rinertia[0] - p.p.rinertia[1];
const double rinertia_diff_12 = p.p.rinertia[1] - p.p.rinertia[2];
@@ -188,7 +190,6 @@ void convert_torques_propagate_omega(const ParticleRange &particles) {
}
}
-/** convert the torques to the body-fixed frames before the integration loop */
void convert_initial_torques(const ParticleRange &particles) {
for (auto &p : particles) {
if (!p.p.rotation)
diff --git a/src/core/rotation.hpp b/src/core/rotation.hpp
index 191269560ae..2472ac776b4 100644
--- a/src/core/rotation.hpp
+++ b/src/core/rotation.hpp
@@ -36,15 +36,17 @@
#include
#include
-/** Propagate angular velocities and update quaternions on a particle */
+/** @brief Propagate angular velocities and update quaternions on a
+ * particle.
+ */
void propagate_omega_quat_particle(Particle &p);
-/** Convert torques to the body-fixed frame and propagate
- angular velocities */
+/** @brief Convert torques to the body-fixed frame and propagate
+ * angular velocities.
+ */
void convert_torques_propagate_omega(const ParticleRange &particles);
-/** Convert torques to the body-fixed frame to start
- the integration loop */
+/** Convert torques to the body-fixed frame before the integration loop. */
void convert_initial_torques(const ParticleRange &particles);
// Frame conversion routines