Merge branch 'python' into propagators

doc/sphinx/installation.rst

@@ -442,15 +442,15 @@ General features
 -  ``THERMOSTAT_PER_PARTICLE`` Allows setting a per-particle friction
    coefficient for the Langevin and Brownian thermostats.
 
--  ``ROTATIONAL_INERTIA`` Allows to set the eigenvalues for rotational inertia matrix of particles
+-  ``ROTATIONAL_INERTIA`` Allows particles to have individual rotational inertia matrix eigenvalues.
+   When not built in, all eigenvalues are unity in simulation units.
 
 -  ``EXTERNAL_FORCES`` Allows to define an arbitrary constant force for each particle
    individually. Also allows to fix individual coordinates of particles,
    keep them at a fixed position or within a plane.
 
--  ``MASS`` Allows particles to have individual masses. Note that some analysis
-   procedures have not yet been adapted to take the masses into account
-   correctly.
+-  ``MASS`` Allows particles to have individual masses.
+   When not built in, all masses are unity in simulation units.
 
    .. seealso:: :attr:`espressomd.particle_data.ParticleHandle.mass`
 

doc/sphinx/integration.rst

@@ -199,6 +199,7 @@ Notes:
 * The particle forces :math:`F` include interactions as well as a friction (:math:`\gamma^0`) and noise term (:math:`\sqrt{k_B T \gamma^0 dt} \overline{\eta}`) analogous to the terms in the :ref:`Langevin thermostat`.
 * The particle forces are only calculated in step 5 and then reused in step 1 of the next iteration. See :ref:`Velocity Verlet Algorithm` for the implications of that.
 * The NpT algorithm doesn't support :ref:`Lees-Edwards boundary conditions`.
+* The NpT algorithm doesn't support propagation of angular velocities.
 
 .. _Steepest descent:
 

src/core/electrostatics/p3m.cpp

@@ -492,6 +492,8 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
     /* sqrt(-1)*k differentiation */
     int j[3];
     int ind = 0;
+    /* compute electric field */
+    // Eq. (3.49) @cite deserno00b
     for (j[0] = 0; j[0] < p3m.fft.plan[3].new_mesh[0]; j[0]++) {
       for (j[1] = 0; j[1] < p3m.fft.plan[3].new_mesh[1]; j[1]++) {
         for (j[2] = 0; j[2] < p3m.fft.plan[3].new_mesh[2]; j[2]++) {
@@ -535,6 +537,7 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
         boost::combine(p_q_range, p_force_range));
 
     // add dipole forces
+    // Eq. (3.19) @cite deserno00b
     if (box_dipole) {
       auto const dm = prefactor * pref * box_dipole.value();
       for (auto zipped : boost::combine(p_q_range, p_force_range)) {
@@ -550,6 +553,7 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
     auto node_energy = 0.;
     for (int i = 0; i < p3m.fft.plan[3].new_size; i++) {
       // Use the energy optimized influence function for energy!
+      // Eq. (3.40) @cite deserno00b
       node_energy += p3m.g_energy[i] * (Utils::sqr(p3m.rs_mesh[2 * i]) +
                                         Utils::sqr(p3m.rs_mesh[2 * i + 1]));
     }
@@ -559,11 +563,14 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
     boost::mpi::reduce(comm_cart, node_energy, energy, std::plus<>(), 0);
     if (this_node == 0) {
       /* self energy correction */
+      // Eq. (3.8) @cite deserno00b
       energy -= p3m.sum_q2 * p3m.params.alpha * Utils::sqrt_pi_i();
       /* net charge correction */
+      // Eq. (3.11) @cite deserno00b
       energy -= p3m.square_sum_q * Utils::pi() /
                 (2. * volume * Utils::sqr(p3m.params.alpha));
       /* dipole correction */
+      // Eq. (3.9) @cite deserno00b
       if (box_dipole) {
         energy += pref * box_dipole.value().norm2();
       }

src/core/electrostatics/p3m.hpp

@@ -202,6 +202,7 @@ struct CoulombP3M : public Coulomb::Actor<CoulombP3M> {
   }
 
   /** Calculate real-space contribution of Coulomb pair energy. */
+  // Eq. (3.6) @cite deserno00b
   double pair_energy(double q1q2, double dist) const {
     if ((q1q2 == 0.) || dist >= p3m.params.r_cut || dist <= 0.) {
       return {};

src/core/energy_inline.hpp

@@ -305,7 +305,7 @@ inline double translational_kinetic_energy(Particle const &p) {
  */
 inline double rotational_kinetic_energy(Particle const &p) {
 #ifdef ROTATION
-  return (not p.is_virtual() and p.can_rotate())
+  return (p.can_rotate() and not p.is_virtual())
              ? 0.5 * (hadamard_product(p.omega(), p.omega()) * p.rinertia())
              : 0.0;
 #else

src/core/integrators/velocity_verlet_npt.cpp

@@ -30,7 +30,6 @@
 #include "errorhandling.hpp"
 #include "integrate.hpp"
 #include "npt.hpp"
-#include "rotation.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
 #include "thermostats/npt_inline.hpp"

src/python/espressomd/particle_data.py

@@ -179,7 +179,7 @@ class ParticleHandle(ScriptInterfaceHelper):
     rinertia: (3,) array_like of :obj:`float`
         The particle rotational inertia.
 
-        Sets the diagonal elements of this particles rotational inertia
+        Sets the diagonal elements of this particle's rotational inertia
         tensor. These correspond with the inertial moments along the
         coordinate axes in the particle's co-rotating coordinate system.
         When the particle's quaternions are set to ``[1, 0, 0, 0,]``, the

testsuite/python/elc.py

@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2010-2022 The ESPResSo project
+# Copyright (C) 2010-2023 The ESPResSo project
 #
 # This file is part of ESPResSo.
 #
@@ -22,15 +22,13 @@
 import espressomd.electrostatics
 
 import numpy as np
+import itertools
 
-GAP = np.array([0., 0., 3.])
-BOX_L = np.array(3 * [10.]) + GAP
 TIME_STEP = 1e-100
-POTENTIAL_DIFFERENCE = -3.
 
 
 class ElcTest:
-    system = espressomd.System(box_l=BOX_L, time_step=TIME_STEP)
+    system = espressomd.System(box_l=[1.] * 3, time_step=TIME_STEP)
     system.cell_system.skin = 0.0
 
     def tearDown(self):
@@ -40,6 +38,12 @@ def tearDown(self):
     def test_finite_potential_drop(self):
         system = self.system
 
+        GAP = np.array([0., 0., 3.])
+        BOX_L = np.array(3 * [10.]) + GAP
+        POTENTIAL_DIFFERENCE = -3.
+
+        system.box_l = BOX_L
+
         p1 = system.part.add(pos=[0, 0, 1], q=+1)
         p2 = system.part.add(pos=[0, 0, 9], q=-1)
 
@@ -90,6 +94,48 @@ def test_finite_potential_drop(self):
         with self.assertRaisesRegex(Exception, 'entered ELC gap region'):
             self.system.integrator.run(2)
 
+    def test_elc_p3m_madelung(self):
+        system = self.system
+
+        n_pairs = 6
+
+        BOX_L = 2 * n_pairs
+        ELC_GAP = 0.4 * BOX_L
+        system.box_l = [BOX_L, BOX_L, BOX_L + ELC_GAP]
+
+        for j, k, l in itertools.product(range(2 * n_pairs), repeat=3):
+            system.part.add(pos=[j + 0.5, k + 0.5, l + 0.5],
+                            q=(-1)**(j + k + l))
+
+        p3m = self.p3m_class(
+            prefactor=1,
+            mesh=[60, 60, 84],
+            cao=7,
+            r_cut=3.314,
+            alpha=1.18,
+            accuracy=3e-7,
+            tune=False,
+        )
+        elc = espressomd.electrostatics.ELC(
+            actor=p3m,
+            gap_size=ELC_GAP,
+            maxPWerror=3e-7,
+            delta_mid_top=-1,
+            delta_mid_bot=-1,
+            const_pot=True,
+            pot_diff=0,
+        )
+
+        system.electrostatics.solver = elc
+
+        MADELUNG = -1.74756459463318219
+        U_expected = MADELUNG
+
+        U_elc = 2. * \
+            system.analysis.energy()['coulomb'] / len(system.part.all())
+
+        np.testing.assert_allclose(U_elc, U_expected, atol=0., rtol=1e-6)
+
 
 @utx.skipIfMissingFeatures(["P3M"])
 class ElcTestCPU(ElcTest, ut.TestCase):

testsuite/python/integrator_npt.py

@@ -207,8 +207,8 @@ def run_with_p3m(self, container, p3m, method):
         system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
         system.integrator.run(10)
         # check runtime warnings
-        self.system.thermostat.turn_off()
-        self.system.integrator.set_vv()
+        system.thermostat.turn_off()
+        system.integrator.set_vv()
         err_msg = f"If {method} is being used you must use the cubic box NpT"
         with self.assertRaisesRegex(RuntimeError, err_msg):
             system.integrator.set_isotropic_npt(**npt_kwargs_rectangular)