From c7ea55777de1d90eefa6ba7a14112cb2e3e84561 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 13 Apr 2023 17:21:58 +0200 Subject: [PATCH 1/4] core: Fix particle folding bug for vs_relative Co-authored-by: Rudolf Weeber --- src/core/BoxGeometry.hpp | 2 +- src/core/virtual_sites/VirtualSitesRelative.cpp | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/core/BoxGeometry.hpp b/src/core/BoxGeometry.hpp index a8af7fc8fef..81503f9976d 100644 --- a/src/core/BoxGeometry.hpp +++ b/src/core/BoxGeometry.hpp @@ -187,7 +187,7 @@ class BoxGeometry { Utils::Vector get_mi_vector(const Utils::Vector &a, const Utils::Vector &b) const { if (type() == BoxType::LEES_EDWARDS) { - auto const &shear_plane_normal = lees_edwards_bc().shear_plane_normal; + auto const shear_plane_normal = lees_edwards_bc().shear_plane_normal; auto a_tmp = a; auto b_tmp = b; a_tmp[shear_plane_normal] = Algorithm::periodic_fold( diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 04f64ededb7..7a7b154cb89 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -115,12 +115,15 @@ void VirtualSitesRelative::update() const { continue; auto const &p_ref = *p_ref_ptr; + p.image_box() = p_ref.image_box(); p.pos() = p_ref.pos() + connection_vector(p_ref, p); p.v() = velocity(p_ref, p); if (box_geo.type() == BoxType::LEES_EDWARDS) { auto push = LeesEdwards::Push(box_geo); - push(p, -1); + push(p, -1); // includes a position fold + } else { + fold_position(p.pos(), p.image_box(), box_geo); } if (have_quaternions()) From 144a6cf3922ced8e13f40a0699413dc4adb63bc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 13 Apr 2023 17:26:37 +0200 Subject: [PATCH 2/4] core: Fix ghost communication Particles tracked by relative virtual sites need to communicate their image_box to ensure proper coordinates folding when the virtual site gets the real particle information from a ghost layer. --- src/core/Particle.hpp | 10 +++++----- src/core/ghosts.cpp | 6 ++++++ 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index b57c88552ce..53aca298245 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -265,6 +265,8 @@ struct ParticleProperties { struct ParticlePosition { /** periodically folded position. */ Utils::Vector3d p = {0., 0., 0.}; + /** index of the simulation box image where the particle really sits. */ + Utils::Vector3i i = {0, 0, 0}; #ifdef ROTATION /** quaternion to define particle orientation */ @@ -282,6 +284,7 @@ struct ParticlePosition { template void serialize(Archive &ar, long int /* version */) { ar &p; + ar &i; #ifdef ROTATION ar &quat; #endif @@ -363,8 +366,6 @@ struct ParticleLocal { /** is particle a ghost particle. */ bool ghost = false; short int lees_edwards_flag = 0; - /** index of the simulation box image where the particle really sits. */ - Utils::Vector3i i = {0, 0, 0}; /** position from the last Verlet list update. */ Utils::Vector3d p_old = {0., 0., 0.}; /** Accumulated applied Lees-Edwards offset. */ @@ -373,7 +374,6 @@ struct ParticleLocal { template void serialize(Archive &ar, long int /* version */) { ar &ghost; ar &lees_edwards_flag; - ar &i; ar &p_old; ar &lees_edwards_offset; } @@ -446,8 +446,8 @@ struct Particle { // NOLINT(bugprone-exception-escape) void set_ghost(bool const ghost_flag) { l.ghost = ghost_flag; } auto &pos_at_last_verlet_update() { return l.p_old; } auto const &pos_at_last_verlet_update() const { return l.p_old; } - auto const &image_box() const { return l.i; } - auto &image_box() { return l.i; } + auto const &image_box() const { return r.i; } + auto &image_box() { return r.i; } auto const &lees_edwards_offset() const { return l.lees_edwards_offset; } auto &lees_edwards_offset() { return l.lees_edwards_offset; } auto const &lees_edwards_flag() const { return l.lees_edwards_flag; } diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 55ccc6f3e08..900f59b7759 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -29,7 +29,9 @@ * - a "GhostCommunication" is always named "ghost_comm". */ #include "ghosts.hpp" + #include "Particle.hpp" +#include "grid.hpp" #include #include @@ -168,9 +170,13 @@ serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts, if (direction == SerializationDirection::SAVE and ghost_shift != nullptr) { /* ok, this is not nice, but perhaps fast */ auto pos = p.pos() + *ghost_shift; + auto img = p.image_box(); + fold_position(pos, img, ::box_geo); ar &pos; + ar &img; } else { ar &p.pos(); + ar &p.image_box(); } #ifdef ROTATION ar &p.quat(); From 7dd8db258aa43e3fb7a24da27aafe558d967fd23 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 13 Apr 2023 17:42:25 +0200 Subject: [PATCH 3/4] core: Check for dangling virtual sites --- src/core/Particle.hpp | 2 +- src/core/virtual_sites/VirtualSitesRelative.cpp | 5 +++++ testsuite/python/particle_slice.py | 6 +++--- testsuite/python/test_checkpoint.py | 11 ++++++++++- testsuite/python/virtual_sites_relative.py | 9 ++++++++- 5 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 53aca298245..7b0520a65b1 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -161,7 +161,7 @@ struct ParticleProperties { * quaternion attribute. */ struct VirtualSitesRelativeParameters { - int to_particle_id = 0; + int to_particle_id = -1; double distance = 0.; /** Relative position of the virtual site. */ Utils::Quaternion rel_orientation = diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 7a7b154cb89..7da96d1221b 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -82,6 +82,11 @@ static Particle *get_reference_particle(Particle const &p) { return nullptr; } auto const &vs_rel = p.vs_relative(); + if (vs_rel.to_particle_id == -1) { + runtimeErrorMsg() << "Particle with id " << p.id() + << " is a dangling virtual site"; + return nullptr; + } auto p_ref_ptr = cell_structure.get_local_particle(vs_rel.to_particle_id); if (!p_ref_ptr) { runtimeErrorMsg() << "No real particle with id " << vs_rel.to_particle_id diff --git a/testsuite/python/particle_slice.py b/testsuite/python/particle_slice.py index 4cd2aaf0a8a..4a9332047e0 100644 --- a/testsuite/python/particle_slice.py +++ b/testsuite/python/particle_slice.py @@ -336,9 +336,9 @@ def test_vs_relative(self): self.assertEqual(repr(all_partcls.vs_relative), repr([[1, 1.0, np.array([1., 1., 1., 1.])], - [0, 0.0, np.array([1., 0., 0., 0.])], - [0, 0.0, np.array([1., 0., 0., 0.])], - [0, 0.0, np.array([1., 0., 0., 0.])]])) + [-1, 0.0, np.array([1., 0., 0., 0.])], + [-1, 0.0, np.array([1., 0., 0., 0.])], + [-1, 0.0, np.array([1., 0., 0., 0.])]])) all_partcls.vs_relative = [1, 1.0, (1.0, 1.0, 1.0, 1.0)] diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 532cdc32cc0..2174d3b5a53 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -496,12 +496,21 @@ def test_drude_helpers(self): @utx.skipIfMissingFeatures(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']) def test_virtual_sites(self): - self.assertTrue(system.part.by_id(1).virtual) self.assertIsInstance( system.virtual_sites, espressomd.virtual_sites.VirtualSitesRelative) self.assertTrue(system.virtual_sites.have_quaternion) self.assertTrue(system.virtual_sites.override_cutoff_check) + p_real = system.part.by_id(0) + p_virt = system.part.by_id(1) + self.assertTrue(p_virt.virtual) + self.assertFalse(p_real.virtual) + self.assertEqual(p_real.vs_relative[0], -1) + self.assertEqual(p_virt.vs_relative[0], p_real.id) + self.assertEqual(p_real.vs_relative[1], 0.) + self.assertEqual(p_virt.vs_relative[1], np.sqrt(2.)) + np.testing.assert_allclose( + p_real.vs_relative[2], [1., 0., 0., 0.], atol=1e-10) def test_mean_variance_calculator(self): np.testing.assert_array_equal( diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index eda97ebafe0..6a575542d83 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -57,7 +57,7 @@ def director_from_quaternion(self, quat): - quat[2] * quat[2] + quat[3] * quat[3]))) def verify_vs(self, vs, verify_velocity=True): - """Verify vs position and (if compiled in) velocity.""" + """Verify virtual site position and velocity.""" self.assertTrue(vs.virtual) vs_r = vs.vs_relative @@ -156,6 +156,13 @@ def test_vs_exceptions(self): p1 = system.part.add(pos=[0.0, 0.0, 0.0], rotation=3 * [True], id=1) p2 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=3 * [True], id=2) p3 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=3 * [True], id=3) + p4 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=3 * [True], id=4) + # dangling virtual sites are not allowed + with self.assertRaisesRegex(Exception, "Particle with id 4 is a dangling virtual site"): + p4.virtual = True + self.assertEqual(p4.vs_relative[0], -1) + system.integrator.run(0, recalc_forces=True) + p4.remove() # relating to anything else other than a particle or id is not allowed with self.assertRaisesRegex(ValueError, "Argument of 'vs_auto_relate_to' has to be of type ParticleHandle or int"): p2.vs_auto_relate_to('0') From 86bf9c7c2dab6b9eb1ab5e7baa9a04c03a4a7f8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 14 Apr 2023 00:07:00 +0200 Subject: [PATCH 4/4] tests: Check virtual sites coordinate folding --- testsuite/python/CMakeLists.txt | 1 + .../python/virtual_sites_relative_pbc.py | 198 ++++++++++++++++++ 2 files changed, 199 insertions(+) create mode 100644 testsuite/python/virtual_sites_relative_pbc.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index ca149768c89..1072822819e 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -267,6 +267,7 @@ python_test(FILE brownian_dynamics_stats.py MAX_NUM_PROC 1 LABELS long) python_test(FILE lees_edwards.py MAX_NUM_PROC 4) python_test(FILE nsquare.py MAX_NUM_PROC 4) python_test(FILE virtual_sites_relative.py MAX_NUM_PROC 2) +python_test(FILE virtual_sites_relative_pbc.py MAX_NUM_PROC 2) python_test(FILE virtual_sites_tracers.py MAX_NUM_PROC 2 DEPENDENCIES virtual_sites_tracers_common.py) python_test(FILE virtual_sites_tracers_gpu.py MAX_NUM_PROC 2 GPU_SLOTS 1 diff --git a/testsuite/python/virtual_sites_relative_pbc.py b/testsuite/python/virtual_sites_relative_pbc.py new file mode 100644 index 00000000000..745dc6c7bd6 --- /dev/null +++ b/testsuite/python/virtual_sites_relative_pbc.py @@ -0,0 +1,198 @@ +# +# Copyright (C) 2023 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 . +# + +import unittest as ut +import unittest_decorators as utx +import espressomd +import espressomd.virtual_sites +import espressomd.lees_edwards +import numpy as np + + +@utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE", "ROTATION"]) +class Test(ut.TestCase): + """ + This test case checks the behavior of a virtual and real particle pair + as they move towards, and eventually cross, a periodic boundary. + Folded and unfolded coordinates are checked, as well as rotations, + with and without Lees-Edwards boundary conditions. + + If this test fails, there must be a bug in the virtual site update method, + the Lees-Edwards update method, or the particle position ghost communicator. + """ + vs_dist = 1. + system = espressomd.System(box_l=[20., 20., 20.]) + system.time_step = 0.01 + system.min_global_cut = 1.5 * vs_dist + system.cell_system.skin = 0.1 + + def setUp(self): + self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() + + def tearDown(self): + self.system.part.clear() + self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesOff() + self.system.lees_edwards.protocol = None + + def check_pbc(self): + system = self.system + vs_dist = self.vs_dist + box_l = system.box_l[0] + director = np.array([1, 0, 0]) + unit_cell = np.array([0, 0, 0]) + eps = 1e-4 + v = (vs_dist / 2. + 1.5 * eps) / system.time_step + le_offset = 0. + if system.lees_edwards.protocol is not None: + le_offset = system.lees_edwards.protocol.initial_pos_offset + + # make sure the periodic boundary (x-axis) uses the ghost communicator + if np.prod(system.cell_system.node_grid) != 1: + assert system.cell_system.node_grid[0] != 1 + # make sure the periodic boundary (x-axis) is the shear boundary + if system.lees_edwards.protocol is not None: + assert system.lees_edwards.shear_plane_normal == "x" + assert system.lees_edwards.shear_direction == "y" + + def check_dist(real_part, virt_part, check_unfolded_dist=True): + dist_folded = system.distance(real_part, virt_part) + self.assertAlmostEqual(dist_folded, vs_dist) + if check_unfolded_dist: + dist_unfolded = np.linalg.norm(real_part.pos - virt_part.pos) + self.assertAlmostEqual(dist_unfolded, vs_dist) + + def check_pos(p, ref_image_box, ref_pos): + np.testing.assert_allclose(np.copy(p.pos), ref_pos, atol=1e-10) + np.testing.assert_array_equal(np.copy(p.image_box), ref_image_box) + + for le_dir in (-1., +1.): + if le_dir == -1: + # set one particle near to the left periodic boundary and the + # other one further away by one unit of vs distance + start = 0. + elif le_dir == +1: + # set one particle near to the right periodic boundary and the + # other one further away by one unit of vs distance + start = box_l + le_vec = le_dir * np.array([0., le_offset, 0.]) + image_box = le_dir * director + # trajectory of the particle nearest to the boundary at each step + pos_near = np.zeros((3, 3)) + pos_near[0] = [start - le_dir * + (eps * 1.0 + vs_dist * 0.0), 1., 1.] + pos_near[1] = [start + le_dir * + (eps * 0.5 + vs_dist * 0.5), 1., 1.] + pos_near[2] = [start + le_dir * + (eps * 2.0 + vs_dist * 1.0), 1. - le_dir * le_offset, 1.] + pos_away = pos_near - [le_dir * vs_dist, 0., 0.] + real_kwargs = {"v": le_dir * v * director, "director": director} + virt_kwargs = {"virtual": True} + # case 1: virtual site goes through boundary before real particle + # case 2: virtual site goes through boundary after real particle + # In the second time step, the virtual site and real particle are + # in the same image box. + real_part1 = system.part.add(pos=pos_away[0], **real_kwargs) + virt_part1 = system.part.add(pos=pos_near[0], **virt_kwargs) + real_part2 = system.part.add(pos=pos_near[0], **real_kwargs) + virt_part2 = system.part.add(pos=pos_away[0], **virt_kwargs) + virt_part1.vs_auto_relate_to(real_part1) + virt_part2.vs_auto_relate_to(real_part2) + system.integrator.run(0) + check_dist(real_part1, virt_part1) + check_dist(real_part2, virt_part2) + check_pos(real_part1, unit_cell, pos_away[0]) + check_pos(virt_part1, unit_cell, pos_near[0]) + check_pos(real_part2, unit_cell, pos_near[0]) + check_pos(virt_part2, unit_cell, pos_away[0]) + system.integrator.run(1) + check_dist(real_part1, virt_part1, check_unfolded_dist=False) + check_dist(real_part2, virt_part2, check_unfolded_dist=False) + check_pos(real_part1, unit_cell, pos_away[1] + 0 * le_vec) + check_pos(virt_part1, image_box, pos_near[1] + 1 * le_vec) + check_pos(real_part2, image_box, pos_near[1] - 1 * le_vec) + check_pos(virt_part2, unit_cell, pos_away[1] - 2 * le_vec) + system.integrator.run(1) + check_dist(real_part1, virt_part1) + check_dist(real_part2, virt_part2) + check_pos(real_part1, image_box, pos_away[2]) + check_pos(virt_part1, image_box, pos_near[2]) + check_pos(real_part2, image_box, pos_near[2]) + check_pos(virt_part2, image_box, pos_away[2]) + + system.part.clear() + + # case 1: virtual site goes through boundary before real particle + # case 2: virtual site goes through boundary after real particle + # Afterwards, the real particle director changes direction to bring + # the virtual site in the same image box as the real particle. + vs_offset = le_dir * director * vs_dist + real_part1 = system.part.add(pos=pos_away[0], **real_kwargs) + virt_part1 = system.part.add(pos=pos_near[0], **virt_kwargs) + real_part2 = system.part.add(pos=pos_near[0], **real_kwargs) + virt_part2 = system.part.add(pos=pos_away[0], **virt_kwargs) + virt_part1.vs_auto_relate_to(real_part1) + virt_part2.vs_auto_relate_to(real_part2) + system.integrator.run(1) + # flip director + real_part1.director = -real_part1.director + real_part2.director = -real_part2.director + # all virtual sites rotate by pi around the real particle + system.integrator.run(0) + check_dist(real_part1, virt_part1) + check_dist(real_part2, virt_part2) + check_pos(real_part1, unit_cell, pos_away[1]) + check_pos(virt_part1, unit_cell, pos_away[1] - vs_offset) + check_pos(real_part2, image_box, pos_near[1] - le_vec) + check_pos(virt_part2, image_box, pos_near[1] - le_vec + vs_offset) + + system.part.clear() + + def test_pbc(self): + system = self.system + + with self.subTest(msg='without Lees-Edwards boundary conditions'): + self.check_pbc() + + system.part.clear() + + # add Lees-Edwards boundary conditions + protocol = espressomd.lees_edwards.LinearShear( + initial_pos_offset=0.1, time_0=0., shear_velocity=0.) + system.lees_edwards.set_boundary_conditions( + shear_direction="y", shear_plane_normal="x", protocol=protocol) + + with self.subTest(msg='with Lees-Edwards boundary conditions'): + self.check_pbc() + + def test_image_box_update(self): + system = self.system + # virtual site image box is overriden by the real particle image box + real_part = system.part.add(pos=[0., 0., +self.vs_dist / 2.]) + virt_part = system.part.add( + pos=[0., 0., -self.vs_dist / 2. + 4. * system.box_l[2]], virtual=True) + virt_part.vs_auto_relate_to(real_part) + np.testing.assert_array_equal(np.copy(real_part.image_box), [0, 0, 0]) + np.testing.assert_array_equal(np.copy(virt_part.image_box), [0, 0, 3]) + system.integrator.run(0) + np.testing.assert_array_equal(np.copy(real_part.image_box), [0, 0, 0]) + np.testing.assert_array_equal(np.copy(virt_part.image_box), [0, 0, -1]) + + +if __name__ == "__main__": + ut.main()