tests: Check virtual sites coordinate folding

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/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()