From d61db70c88a5622533e49196b120af9df7b19c5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 7 Dec 2021 19:54:22 +0100 Subject: [PATCH 1/3] tests: Check Verlet list update --- src/core/EspressoSystemStandAlone.cpp | 7 + src/core/EspressoSystemStandAlone.hpp | 1 + src/core/unit_tests/CMakeLists.txt | 2 + src/core/unit_tests/Verlet_list_test.cpp | 246 +++++++++++++++++++++++ 4 files changed, 256 insertions(+) create mode 100644 src/core/unit_tests/Verlet_list_test.cpp diff --git a/src/core/EspressoSystemStandAlone.cpp b/src/core/EspressoSystemStandAlone.cpp index 8ff2e8d4a2e..d57fd206a25 100644 --- a/src/core/EspressoSystemStandAlone.cpp +++ b/src/core/EspressoSystemStandAlone.cpp @@ -56,6 +56,13 @@ void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const { mpi_set_box_length(box_l); } +void EspressoSystemStandAlone::set_node_grid( + Utils::Vector3i const &node_grid) const { + if (!head_node) + return; + mpi_set_node_grid(node_grid); +} + void EspressoSystemStandAlone::set_time_step(double time_step) const { if (!head_node) return; diff --git a/src/core/EspressoSystemStandAlone.hpp b/src/core/EspressoSystemStandAlone.hpp index 049861aea93..82cb94fdadf 100644 --- a/src/core/EspressoSystemStandAlone.hpp +++ b/src/core/EspressoSystemStandAlone.hpp @@ -28,6 +28,7 @@ class EspressoSystemStandAlone { public: EspressoSystemStandAlone(int argc, char **argv); void set_box_l(Utils::Vector3d const &box_l) const; + void set_node_grid(Utils::Vector3i const &node_grid) const; void set_time_step(double time_step) const; void set_skin(double new_skin) const; diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 6f9c701753c..263f303305f 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -47,6 +47,8 @@ unit_test(NAME BoxGeometry_test SRC BoxGeometry_test.cpp DEPENDS EspressoCore) unit_test(NAME LocalBox_test SRC LocalBox_test.cpp DEPENDS EspressoCore) unit_test(NAME Lattice_test SRC Lattice_test.cpp DEPENDS EspressoCore) unit_test(NAME lb_exceptions SRC lb_exceptions.cpp DEPENDS EspressoCore) +unit_test(NAME Verlet_list_test SRC Verlet_list_test.cpp DEPENDS EspressoCore + NUM_PROC 4) unit_test(NAME thermostats_test SRC thermostats_test.cpp DEPENDS EspressoCore) unit_test(NAME random_test SRC random_test.cpp DEPENDS EspressoUtils Random123) unit_test(NAME BondList_test SRC BondList_test.cpp DEPENDS EspressoCore) diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp new file mode 100644 index 00000000000..c5b9e202a93 --- /dev/null +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -0,0 +1,246 @@ +/* + * Copyright (C) 2021 The ESPResSo project + * + * 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 . + */ + +#define BOOST_TEST_MODULE Verlet list update test + +#include "config.hpp" + +#ifdef LENNARD_JONES + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_NO_MAIN +#define BOOST_TEST_ALTERNATIVE_INIT_API +#include +#include +#include +namespace utf = boost::unit_test; +namespace bdata = boost::unit_test::data; + +#include "EspressoSystemStandAlone.hpp" +#include "Particle.hpp" +#include "ParticleFactory.hpp" +#include "integrate.hpp" +#include "nonbonded_interactions/lj.hpp" +#include "particle_data.hpp" +#include "thermostat.hpp" + +#include + +#include + +#include +#include +#include +#include +#include + +namespace espresso { +// ESPResSo system instance +std::unique_ptr system; +} // namespace espresso + +/** Decorator to run a unit test only on the head node. */ +struct if_head_node { + boost::test_tools::assertion_result operator()(utf::test_unit_id) { + return world.rank() == 0; + } + +private: + boost::mpi::communicator world; +}; + +namespace Testing { +/** + * Helper class to setup an integrator and particle properties such that the + * particle is on a collision course with another particle on the x-axis. + */ +struct IntegratorHelper { + IntegratorHelper() = default; + virtual ~IntegratorHelper() = default; + /** Set integrator parameters. */ + virtual void set_integrator() const = 0; + /** Set particle to move along the x-axis. */ + virtual void set_particle_properties(int) const = 0; + virtual char const *name() const = 0; + friend auto operator<<(std::ostream &os, IntegratorHelper const &obj) + -> std::ostream & { + return os << obj.name(); + } +}; + +struct : public IntegratorHelper { + void set_integrator() const override { + mpi_set_thermo_switch(THERMO_OFF); + integrate_set_steepest_descent(0., 0.01, 100.); + } + void set_particle_properties(int pid) const override { + set_particle_ext_force(pid, {20., 0., 0.}); + } + char const *name() const override { return "SteepestDescent"; } +} steepest_descent; + +struct : public IntegratorHelper { + void set_integrator() const override { + mpi_set_thermo_switch(THERMO_OFF); + integrate_set_nvt(); + } + void set_particle_properties(int pid) const override { + set_particle_v(pid, {20., 0., 0.}); + } + char const *name() const override { return "VelocityVerlet"; } +} velocity_verlet; + +#ifdef NPT +struct : public IntegratorHelper { + void set_integrator() const override { + integrate_set_npt_isotropic(1., 1e9, true, true, true, true); + mpi_set_temperature(1.); + npt_iso_set_rng_seed(0); + mpi_set_thermo_switch(thermo_switch | THERMO_NPT_ISO); + mpi_set_nptiso_gammas(0., 0.); // disable box volume change + } + void set_particle_properties(int pid) const override { + set_particle_v(pid, {20., 0., 0.}); + } + char const *name() const override { return "VelocityVerletNpT"; } +} velocity_verlet_npt; +#endif // NPT +} // namespace Testing + +inline double get_dist_from_last_verlet_update(Particle const &p) { + return (p.r.p - p.l.p_old).norm(); +} + +inline double get_dist_from_pair(Particle const &p1, Particle const &p2) { + return (p1.r.p - p2.r.p).norm(); +} + +auto const node_grids = std::vector{{4, 1, 1}, {2, 2, 1}}; +auto const propagators = + std::vector>{ + Testing::velocity_verlet, +#ifdef NPT + Testing::velocity_verlet_npt, +#endif + Testing::steepest_descent}; + +BOOST_TEST_DECORATOR(*utf::precondition(if_head_node())) +BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update, + bdata::make(node_grids) * bdata::make(propagators), + node_grid, integration_helper) { + constexpr auto tol = 100. * std::numeric_limits::epsilon(); + boost::mpi::communicator world; + + auto const box_l = 8.; + auto const box_center = box_l / 2.; + espresso::system->set_box_l(Utils::Vector3d::broadcast(box_l)); + espresso::system->set_node_grid(node_grid); + + // particle properties + auto const pid1 = 9; + auto const pid2 = 2; + auto const pid3 = 5; + + // distance between particles + auto const dist = 0.2; + // set up LJ potential + auto const eps = 1.0; + auto const sig = 0.05; + auto const shift = 0.0; + auto const offset = 0.1; + auto const min = 0.0; + auto const r_off = dist - offset; + auto const cut = r_off + 1e-3; + lennard_jones_set_params(0, 1, eps, sig, cut, shift, offset, min); + + // set up velocity-Verlet integrator + auto const time_step = 0.01; + auto const skin = 0.1; + espresso::system->set_time_step(time_step); + espresso::system->set_skin(skin); + integration_helper.get().set_integrator(); + + // If the Verlet list is not updated, two particles initially placed in + // different cells will never see each other, even when closer than the + // interaction cutoff. To reproduce that bug, we check the case where + // we have 4 cells in one direction and particles in cells #0 and #2 + // (i.e. the non-neighboring cell case) and the case where we have + // 2 cells in one direction (i.e. the neighboring cell case). + create_particle({2. - 0.10, 1., 1.}, pid1, 0); + create_particle({4. + 0.15, 1., 1.}, pid2, 1); + create_particle({2. - 0.10, 5., 1.}, pid3, 0); + BOOST_REQUIRE_EQUAL(get_particle_node(pid1), 0); + BOOST_REQUIRE_EQUAL(get_particle_node(pid2), 2); + BOOST_REQUIRE_EQUAL(get_particle_node(pid3), (node_grid[0] == 4) ? 0 : 1); + + // check that particles in different cells will eventually interact during + // normal integration (the integration helper sets a collision course) + { + integration_helper.get().set_particle_properties(pid1); + + // integrate until both particles are closer than cutoff + { + mpi_integrate(11, 0); + auto const &p1 = get_particle_data(pid1); + auto const &p2 = get_particle_data(pid2); + BOOST_REQUIRE_LT(get_dist_from_pair(p1, p2), cut); + } + + // check forces and Verlet update + { + mpi_integrate(1, 0); + auto const &p1 = get_particle_data(pid1); + auto const &p2 = get_particle_data(pid2); + BOOST_CHECK_CLOSE(p1.f.f[0] - p1.p.ext_force[0], 480., 1e-9); + BOOST_CHECK_CLOSE(p1.f.f[1], 0., tol); + BOOST_CHECK_CLOSE(p1.f.f[2], 0., tol); + BOOST_TEST(p1.f.f - p1.p.ext_force == -p2.f.f, + boost::test_tools::per_element()); + BOOST_CHECK_LT(get_dist_from_last_verlet_update(p1), skin / 2.); + } + } + + // check that particles in different cells will interact when manually + // placed next to each other (@c place_particle() resorts particles) + { + place_particle(pid3, {4. + 0.10, 1., 1.0}); + { + auto const &p2 = get_particle_data(pid2); + auto const &p3 = get_particle_data(pid3); + BOOST_REQUIRE_LT(get_dist_from_pair(p2, p3), cut); + BOOST_CHECK_GT(get_dist_from_last_verlet_update(p3), skin / 2.); + } + { + mpi_integrate(0, 0); + auto const &p3 = get_particle_data(pid3); + BOOST_CHECK_LT(get_dist_from_last_verlet_update(p3), skin / 2.); + } + } +} + +int main(int argc, char **argv) { + espresso::system = std::make_unique(argc, argv); + // the test case only works for 4 MPI ranks + boost::mpi::communicator world; + if (world.size() == 4) + return boost::unit_test::unit_test_main(init_unit_test, argc, argv); +} +#else // ifdef LENNARD_JONES +int main(int argc, char **argv) {} +#endif From 2a23bd13738d1f0e797a75dc6daf132a52c9e96d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 7 Dec 2021 20:42:09 +0100 Subject: [PATCH 2/3] core: Centralize particle resort logic Also take virtual particles into account. Co-authored-by: Rudolf Weeber --- src/core/CellStructure.hpp | 16 ++++++++++++++++ src/core/cells.cpp | 9 ++------- src/core/integrate.cpp | 15 ++++++++++++--- src/core/integrators/brownian_inline.hpp | 4 ---- .../integrators/stokesian_dynamics_inline.hpp | 14 ++------------ src/core/integrators/velocity_verlet_inline.hpp | 8 -------- src/core/virtual_sites/VirtualSitesRelative.cpp | 10 ++++++---- .../virtual_sites/lb_inertialess_tracers.cpp | 14 +++++--------- 8 files changed, 43 insertions(+), 47 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index 3d33dc356dd..f5dfc34c33d 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -34,6 +34,8 @@ #include "bond_error.hpp" #include "ghosts.hpp" +#include + #include #include #include @@ -361,6 +363,20 @@ struct CellStructure { */ void clear_resort_particles() { m_resort_particles = Cells::RESORT_NONE; } + /** + * @brief Check whether a particle has moved further than the skin, + * thus requiring a resort. + * @param particles Particles to check + * @param skin Skin + * @return Whether a resort is needed. + */ + bool check_resort_required(ParticleRange const &particles, double skin) { + auto const lim = Utils::sqr(skin / 2.); + return std::any_of( + particles.begin(), particles.end(), + [lim](const auto &p) { return ((p.r.p - p.l.p_old).norm2() > lim); }); + } + /** * @brief Synchronize number of ghosts. */ diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 63b961d8075..ffc74396290 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -216,13 +216,8 @@ REGISTER_CALLBACK(cells_re_init) void mpi_bcast_cell_structure(int cs) { mpi_call_all(cells_re_init, cs); } void check_resort_particles() { - const double skin2 = Utils::sqr(skin / 2.0); - - auto const level = (std::any_of(cell_structure.local_particles().begin(), - cell_structure.local_particles().end(), - [&skin2](Particle const &p) { - return (p.r.p - p.l.p_old).norm2() > skin2; - })) + auto const level = (cell_structure.check_resort_required( + cell_structure.local_particles(), skin)) ? Cells::RESORT_LOCAL : Cells::RESORT_NONE; diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 011ae63ec9b..7cacbcbe261 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -133,17 +133,24 @@ void integrator_sanity_checks() { } } +static void resort_particles_if_needed(ParticleRange &particles) { + if (cell_structure.check_resort_required(particles, skin)) { + cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + } +} + /** @brief Calls the hook for propagation kernels before the force calculation * @return whether or not to stop the integration loop early. */ bool integrator_step_1(ParticleRange &particles) { + bool early_exit = false; switch (integ_switch) { case INTEG_METHOD_STEEPEST_DESCENT: - if (steepest_descent_step(particles)) - return true; // early exit + early_exit = steepest_descent_step(particles); break; case INTEG_METHOD_NVT: velocity_verlet_step_1(particles, time_step); + resort_particles_if_needed(particles); break; #ifdef NPT case INTEG_METHOD_NPT_ISO: @@ -157,12 +164,13 @@ bool integrator_step_1(ParticleRange &particles) { #ifdef STOKESIAN_DYNAMICS case INTEG_METHOD_SD: stokesian_dynamics_step_1(particles, time_step); + resort_particles_if_needed(particles); break; #endif // STOKESIAN_DYNAMICS default: throw std::runtime_error("Unknown value for integ_switch"); } - return false; + return early_exit; } /** Calls the hook of the propagation kernels after force calculation */ @@ -182,6 +190,7 @@ void integrator_step_2(ParticleRange &particles, double kT) { case INTEG_METHOD_BD: // the Ermak-McCammon's Brownian Dynamics requires a single step brownian_dynamics_propagator(brownian, particles, time_step, kT); + resort_particles_if_needed(particles); break; #ifdef STOKESIAN_DYNAMICS case INTEG_METHOD_SD: diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp index 04060514872..ae3fd6055e5 100644 --- a/src/core/integrators/brownian_inline.hpp +++ b/src/core/integrators/brownian_inline.hpp @@ -25,7 +25,6 @@ #include "config.hpp" #include "ParticleRange.hpp" -#include "cells.hpp" #include "integrate.hpp" #include "rotation.hpp" #include "thermostat.hpp" @@ -43,9 +42,6 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, p.m.v = bd_drag_vel(brownian.gamma, p); p.r.p += bd_random_walk(brownian, p, time_step, kT); p.m.v += bd_random_walk_vel(brownian, p); - /* Verlet criterion check */ - if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } #ifdef ROTATION if (!p.p.rotation) diff --git a/src/core/integrators/stokesian_dynamics_inline.hpp b/src/core/integrators/stokesian_dynamics_inline.hpp index 1d3ced0a6e2..238da65e55b 100644 --- a/src/core/integrators/stokesian_dynamics_inline.hpp +++ b/src/core/integrators/stokesian_dynamics_inline.hpp @@ -22,20 +22,14 @@ #include "config.hpp" #ifdef STOKESIAN_DYNAMICS -#include "CellStructure.hpp" #include "ParticleRange.hpp" -#include "cells.hpp" #include "communication.hpp" #include "integrate.hpp" #include "rotation.hpp" #include "stokesian_dynamics/sd_interface.hpp" -#include -#include - inline void stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles, double time_step) { - auto const skin2 = Utils::sqr(0.5 * skin); // Compute new (translational and rotational) velocities propagate_vel_pos_sd(particles, comm_cart, time_step); @@ -46,15 +40,11 @@ inline void stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles, p.r.p += p.m.v * time_step; // Perform rotation - double norm = p.m.omega.norm(); + auto const norm = p.m.omega.norm(); if (norm != 0) { - Utils::Vector3d omega_unit = (1 / norm) * p.m.omega; + auto const omega_unit = (1. / norm) * p.m.omega; local_rotate_particle(p, omega_unit, norm * time_step); } - - // Verlet criterion check - if ((p.r.p - p.l.p_old).norm2() > skin2) - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } } diff --git a/src/core/integrators/velocity_verlet_inline.hpp b/src/core/integrators/velocity_verlet_inline.hpp index ee68daf1916..2ea31122ad4 100644 --- a/src/core/integrators/velocity_verlet_inline.hpp +++ b/src/core/integrators/velocity_verlet_inline.hpp @@ -24,12 +24,9 @@ #include "CellStructure.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" -#include "cells.hpp" #include "integrate.hpp" #include "rotation.hpp" -#include - /** Propagate the velocities and positions. Integration steps before force * calculation of the Velocity Verlet integrator:
\f[ v(t+0.5 \Delta t) = * v(t) + 0.5 \Delta t f(t)/m \f]
\f[ p(t+\Delta t) = p(t) + \Delta t @@ -38,7 +35,6 @@ inline void velocity_verlet_propagate_vel_pos(const ParticleRange &particles, double time_step) { - auto const skin2 = Utils::sqr(0.5 * skin); for (auto &p : particles) { #ifdef ROTATION propagate_omega_quat_particle(p, time_step); @@ -57,10 +53,6 @@ inline void velocity_verlet_propagate_vel_pos(const ParticleRange &particles, p.r.p[j] += time_step * p.m.v[j]; } } - - /* Verlet criterion check*/ - if ((p.r.p - p.l.p_old).norm2() > skin2) - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } } diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index c3e038d2955..85f0bc1ca55 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -153,7 +153,8 @@ void VirtualSitesRelative::update() const { cell_structure.ghosts_update(Cells::DATA_PART_POSITION | Cells::DATA_PART_MOMENTUM); - for (auto &p : cell_structure.local_particles()) { + auto const particles = cell_structure.local_particles(); + for (auto &p : particles) { if (!p.p.is_virtual) continue; @@ -169,10 +170,11 @@ void VirtualSitesRelative::update() const { if (get_have_quaternion()) p.r.quat = orientation(p_ref, p.p.vs_relative); + } - if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); - } // namespace + if (cell_structure.check_resort_required(particles, skin)) { + cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + } } // Distribute forces that have accumulated on virtual particles to the diff --git a/src/core/virtual_sites/lb_inertialess_tracers.cpp b/src/core/virtual_sites/lb_inertialess_tracers.cpp index bce83e17163..ff36a6c6473 100644 --- a/src/core/virtual_sites/lb_inertialess_tracers.cpp +++ b/src/core/virtual_sites/lb_inertialess_tracers.cpp @@ -100,9 +100,7 @@ void IBM_UpdateParticlePositions(ParticleRange const &particles, ParticleVelocitiesFromLB_GPU(particles); #endif - // Do update: Euler - const double skin2 = Utils::sqr(0.5 * skin); - // Loop over particles in local cells + // Euler integrator for (auto &p : particles) { if (p.p.is_virtual) { #ifdef EXTERNAL_FORCES @@ -117,14 +115,12 @@ void IBM_UpdateParticlePositions(ParticleRange const &particles, if (!(p.p.ext_flag & 8)) #endif p.r.p[2] += p.m.v[2] * time_step; - - // Check if the particle might have crossed a box border - const double dist2 = (p.r.p - p.l.p_old).norm2(); - if (dist2 > skin2) { - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); - } } } + + if (cell_structure.check_resort_required(particles, skin)) { + cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + } } /** Put the momentum of a given particle into the LB fluid. */ From 1576e469598d0b724ba550f6bf0f9c43f2688281 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 7 Dec 2021 20:53:51 +0100 Subject: [PATCH 3/3] core: Update Verlet list documentation --- src/core/CellStructure.hpp | 4 ++-- src/core/Particle.hpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index f5dfc34c33d..d792ae99a55 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -364,8 +364,8 @@ struct CellStructure { void clear_resort_particles() { m_resort_particles = Cells::RESORT_NONE; } /** - * @brief Check whether a particle has moved further than the skin, - * thus requiring a resort. + * @brief Check whether a particle has moved further than half the skin + * since the last Verlet list update, thus requiring a resort. * @param particles Particles to check * @param skin Skin * @return Whether a resort is needed. diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 4bdbe713138..004715e240a 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -355,7 +355,7 @@ struct ParticleMomentum { struct ParticleLocal { /** is particle a ghost particle. */ bool ghost = false; - /** position in the last time step before last Verlet list update. */ + /** position from the last Verlet list update. */ Utils::Vector3d p_old = {0, 0, 0}; /** index of the simulation box image where the particle really sits. */ Utils::Vector3i i = {0, 0, 0};