From bf4890824cc7f3fb1a63506b1eb012cbac531aa4 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 8 Dec 2023 15:51:45 +0100 Subject: [PATCH 1/4] testsuite: ELC-IC: added Madelung-energy testcase Co-authored-by: Alexander Schlaich --- testsuite/python/elc.py | 56 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 51 insertions(+), 5 deletions(-) diff --git a/testsuite/python/elc.py b/testsuite/python/elc.py index ec35086b96..f13d940e2f 100644 --- a/testsuite/python/elc.py +++ b/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): From c759ee8a66f3ba3bf0529f607d572704eb5f8f10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 11 Dec 2023 15:05:40 +0100 Subject: [PATCH 2/4] Remove kinetic energy contribution from virtual sites Virtual sites no longer contribute a kinetic energy. The regression was introduced in 4.2.0 by a4d5b7680cf06ee90cbd29341f7516c180134e08. --- doc/sphinx/installation.rst | 8 ++++---- src/core/energy_inline.hpp | 2 +- src/core/unit_tests/energy_test.cpp | 19 ++++++++++++++++++- src/python/espressomd/particle_data.py | 2 +- 4 files changed, 24 insertions(+), 7 deletions(-) diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 0a5edce42f..46ba6fbe7e 100644 --- a/doc/sphinx/installation.rst +++ b/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` diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 657a1ce33c..99d41334cf 100644 --- a/src/core/energy_inline.hpp +++ b/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 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 diff --git a/src/core/unit_tests/energy_test.cpp b/src/core/unit_tests/energy_test.cpp index f04aa0ace3..cd9e9dbf45 100644 --- a/src/core/unit_tests/energy_test.cpp +++ b/src/core/unit_tests/energy_test.cpp @@ -17,7 +17,7 @@ * along with this program. If not, see . */ -#define BOOST_TEST_MODULE tests +#define BOOST_TEST_MODULE energy calculation #define BOOST_TEST_DYN_LINK #include @@ -69,6 +69,23 @@ BOOST_AUTO_TEST_CASE(rotational_kinetic_energy_) { 0.5 * (hadamard_product(p.omega(), p.omega()) * p.rinertia()); BOOST_CHECK_EQUAL(rotational_kinetic_energy(p), expected); } + + // virtual particle + { +#ifdef VIRTUAL_SITES + + Particle p; +#ifdef ROTATIONAL_INERTIA + p.rinertia() = {1., 2., 3.}; +#endif + p.set_virtual(true); + p.omega() = {3., 4., 5.}; + p.set_can_rotate_all_axes(); + + auto const expected = 0.; + BOOST_CHECK_EQUAL(rotational_kinetic_energy(p), expected); +#endif + } #endif } diff --git a/src/python/espressomd/particle_data.py b/src/python/espressomd/particle_data.py index 547fdae8cd..d9b72500f1 100644 --- a/src/python/espressomd/particle_data.py +++ b/src/python/espressomd/particle_data.py @@ -178,7 +178,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 From dd34132da2253da1692e1aa6c875e7d1760c4e97 Mon Sep 17 00:00:00 2001 From: RiccardoFrenner Date: Mon, 4 Dec 2023 21:44:32 +0100 Subject: [PATCH 3/4] Add equation references to P3M code --- src/core/electrostatics/p3m.cpp | 7 +++++++ src/core/electrostatics/p3m.hpp | 1 + 2 files changed, 8 insertions(+) diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index 03758b444a..e62d5ff772 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/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(); } diff --git a/src/core/electrostatics/p3m.hpp b/src/core/electrostatics/p3m.hpp index cdf4dedd26..82506a9359 100644 --- a/src/core/electrostatics/p3m.hpp +++ b/src/core/electrostatics/p3m.hpp @@ -202,6 +202,7 @@ struct CoulombP3M : public Coulomb::Actor { } /** 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 {}; From b61516a874f3c067c0a2fce81d76bd50d6bdd72d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 18 Dec 2023 11:51:30 +0100 Subject: [PATCH 4/4] Disable NpT propagation of angular velocities --- doc/sphinx/integration.rst | 1 + src/core/integrators/velocity_verlet_npt.cpp | 9 ++++---- testsuite/python/integrator_npt.py | 23 +++++++++----------- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/doc/sphinx/integration.rst b/doc/sphinx/integration.rst index cf521a64d6..9e2a4e3967 100644 --- a/doc/sphinx/integration.rst +++ b/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: diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 21a3c9fab0..f73592cb93 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/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" @@ -170,7 +169,10 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles, for (auto &p : particles) { #ifdef ROTATION - propagate_omega_quat_particle(p, time_step); + if (p.can_rotate()) { + runtimeErrorMsg() << "The isotropic NpT integrator doesn't propagate " + "angular velocities"; + } #endif // Don't propagate translational degrees of freedom of vs @@ -201,9 +203,6 @@ void velocity_verlet_npt_step_1(const ParticleRange &particles, void velocity_verlet_npt_step_2(const ParticleRange &particles, double time_step) { velocity_verlet_npt_propagate_vel_final(particles, time_step); -#ifdef ROTATION - convert_torques_propagate_omega(particles, time_step); -#endif velocity_verlet_npt_finalize_p_inst(time_step); } #endif // NPT diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index 5c86768694..919c926d36 100644 --- a/testsuite/python/integrator_npt.py +++ b/testsuite/python/integrator_npt.py @@ -18,7 +18,6 @@ # import unittest as ut import unittest_decorators as utx -import tests_common import numpy as np import espressomd import espressomd.integrate @@ -123,6 +122,14 @@ def test_00_integrator_recovery(self): np.copy(system.part.all().pos), positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]])) + if espressomd.has_features("ROTATION"): + # propagator doesn't handle angular velocities + system.part.all().rotation = [3 * [False], [True, False, False]] + system.thermostat.set_npt(kT=1., gamma0=2., gammav=0.04, seed=42) + system.integrator.set_isotropic_npt(**npt_params) + with self.assertRaisesRegex(Exception, "The isotropic NpT integrator doesn't propagate angular velocities"): + system.integrator.run(1) + def run_with_p3m(self, container, p3m, method): system = self.system npt_kwargs = {"ext_pressure": 0.001, "piston": 0.001} @@ -135,8 +142,6 @@ def run_with_p3m(self, container, p3m, method): pos=np.random.uniform(0, system.box_l[0], (11, 3))) if espressomd.has_features("P3M"): partcl.q = np.sign(np.arange(-5, 6)) - if espressomd.has_features("DP3M"): - partcl.dip = tests_common.random_dipoles(11) system.non_bonded_inter[0, 0].lennard_jones.set_params( epsilon=1, sigma=1, cutoff=2**(1 / 6), shift=0.25) system.integrator.set_steepest_descent( @@ -150,8 +155,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) @@ -164,14 +169,6 @@ def run_with_p3m(self, container, p3m, method): with self.assertRaisesRegex(Exception, err_msg): system.integrator.run(0, recalc_forces=True) - @utx.skipIfMissingFeatures(["DP3M"]) - def test_npt_dp3m_cpu(self): - import espressomd.magnetostatics - dp3m = espressomd.magnetostatics.DipolarP3M( - prefactor=1.0, accuracy=1e-2, mesh=3 * [36], cao=7, r_cut=1.0, - alpha=2.995, tune=False) - self.run_with_p3m(self.system.magnetostatics, dp3m, "magnetostatics") - @utx.skipIfMissingFeatures(["P3M"]) def test_npt_p3m_cpu(self): import espressomd.electrostatics