Rewrite coordinates folding of virtual sites relative

src/core/BoxGeometry.hpp

@@ -187,7 +187,7 @@ class BoxGeometry {
   Utils::Vector<T, 3> get_mi_vector(const Utils::Vector<T, 3> &a,
                                     const Utils::Vector<T, 3> &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(

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<double> rel_orientation =
@@ -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 <class Archive> 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 <class Archive> 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; }

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 <utils/Span.hpp>
 #include <utils/serialization/memcpy_archive.hpp>
@@ -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();

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
@@ -115,12 +120,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())

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

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)]
 

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(

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')

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 <http://www.gnu.org/licenses/>.
+#
+
+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()