Simplify VS_Relative, fix VS support for Lees Edwards (espressomd#4564)

* Remove the less-than-transparent shifting logic which prevents posiiton folding of vs relative. this is no longer needed, since we now use minimum_image_distance everywhere.
* Apply Lees Edwards position shift to VS on update
* Always fold shear plane normal coordinates in Lees Edwards minimum image distance

src/core/BoxGeometry.hpp

@@ -187,7 +187,14 @@ 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) {
-      return lees_edwards_bc().distance(a - b, length(), length_half(),
+      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(
+          a_tmp[shear_plane_normal], length()[shear_plane_normal]);
+      b_tmp[shear_plane_normal] = Algorithm::periodic_fold(
+          b_tmp[shear_plane_normal], length()[shear_plane_normal]);
+      return lees_edwards_bc().distance(a_tmp - b_tmp, length(), length_half(),
                                         length_inv(), m_periodic);
     }
     assert(type() == BoxType::CUBOID);

src/core/integrate.cpp

@@ -135,7 +135,7 @@ void unset_protocol() {
 
 template <class Kernel> void run_kernel() {
   if (box_geo.type() == BoxType::LEES_EDWARDS) {
-    auto const kernel = Kernel{box_geo, time_step};
+    auto const kernel = Kernel{box_geo};
     auto const particles = cell_structure.local_particles();
     std::for_each(particles.begin(), particles.end(),
                   [&kernel](auto &p) { kernel(p); });

src/core/lees_edwards/LeesEdwardsBC.hpp

@@ -39,12 +39,16 @@ struct LeesEdwardsBC {
 
     double n_le_crossings =
         std::round(res[shear_plane_normal] * l_inv[shear_plane_normal]);
-    res[shear_direction] += n_le_crossings * pos_offset;
+    if (n_le_crossings >= 1)
+      res[shear_direction] += pos_offset;
+    if (n_le_crossings <= -1)
+      res[shear_direction] -= pos_offset;
 
     for (int i : {0, 1, 2}) {
-      if (periodic[i])
+      if (periodic[i]) {
         n_jumps[i] = std::round(res[i] * l_inv[i]);
-      res[i] -= n_jumps[i] * l[i];
+        res[i] -= n_jumps[i] * l[i];
+      }
     }
 
     return res;

src/core/lees_edwards/lees_edwards.hpp

@@ -30,33 +30,32 @@ namespace LeesEdwards {
 class UpdateOffset {
 protected:
   LeesEdwardsBC const &m_le;
-  double const m_half_time_step;
 
 public:
-  UpdateOffset(BoxGeometry const &box, double time_step)
-      : m_le{box.lees_edwards_bc()}, m_half_time_step{0.5 * time_step} {}
+  UpdateOffset(BoxGeometry const &box) : m_le{box.lees_edwards_bc()} {}
 
-  void operator()(Particle &p) const {
-    p.lees_edwards_offset() -= m_half_time_step *
-                               p.image_box()[m_le.shear_plane_normal] *
-                               m_le.shear_velocity;
+  void operator()(Particle &p, double pos_prefactor = 1.0) const {
+    // Disabled as long as we do not use a two step LE update
+    //    p.lees_edwards_offset() -= pos_prefactor *
+    //                               static_cast<double>(p.lees_edwards_flag())
+    //                               * m_le.pos_offset / 2;
   }
 };
 
 class Push : public UpdateOffset {
-  Utils::Vector3d const &m_box_l;
+  BoxGeometry const &m_box;
 
 public:
-  Push(BoxGeometry const &box, double time_step)
-      : UpdateOffset{box, time_step}, m_box_l{box.length()} {}
+  Push(BoxGeometry const &box) : UpdateOffset{box}, m_box{box} {}
 
-  void operator()(Particle &p) const {
+  void operator()(Particle &p, double pos_prefactor = 1.0) const {
     // Lees-Edwards acts at the zero step for the velocity half update and at
     // the midstep for the position.
     //
     // The update of the velocity at the end of the time step is triggered by
     // the flag and occurs in @ref propagate_vel_finalize_p_inst
-    if (p.pos()[m_le.shear_plane_normal] >= m_box_l[m_le.shear_plane_normal]) {
+    if (p.pos()[m_le.shear_plane_normal] >=
+        m_box.length()[m_le.shear_plane_normal]) {
       p.lees_edwards_flag() = -1; // perform a negative half velocity shift
     } else if (p.pos()[m_le.shear_plane_normal] < 0) {
       p.lees_edwards_flag() = 1; // perform a positive half velocity shift
@@ -66,9 +65,9 @@ class Push : public UpdateOffset {
 
     auto const dir = static_cast<double>(p.lees_edwards_flag());
     p.v()[m_le.shear_direction] += dir * m_le.shear_velocity;
-    p.pos()[m_le.shear_direction] += dir * m_le.pos_offset;
-    p.lees_edwards_offset() -= dir * m_le.pos_offset;
-    UpdateOffset::operator()(p);
+    p.pos()[m_le.shear_direction] += pos_prefactor * dir * m_le.pos_offset;
+    p.lees_edwards_offset() -= pos_prefactor * dir * m_le.pos_offset;
+    fold_position(p.pos(), p.image_box(), m_box);
   }
 };
 

src/core/unit_tests/grid_test.cpp

@@ -155,7 +155,7 @@ BOOST_AUTO_TEST_CASE(lees_edwards_mi_vector) {
 
     for (int i = 0; i < 3; i++) {
       auto const expected = get_mi_coord(a[i], b[i], box_l[i], box.periodic(i));
-      BOOST_CHECK_SMALL(std::abs(expected - result[i]), epsilon<double>);
+      BOOST_CHECK_SMALL(std::abs(expected - result[i]), 5. * epsilon<double>);
     }
   }
 
@@ -170,14 +170,17 @@ BOOST_AUTO_TEST_CASE(lees_edwards_mi_vector) {
     auto const result = box.get_mi_vector(a, b);
 
     for (int i = 0; i < 3; i++) {
-      auto const expected = get_mi_coord(a[i], b[i], box_l[i], box.periodic(i));
-      BOOST_CHECK_SMALL(std::abs(expected - result[i]), epsilon<double>);
+      auto expected = get_mi_coord(a[i], b[i], box_l[i], box.periodic(i));
+      if (i == le.shear_direction) {
+        expected -= le.pos_offset = 1.;
+      }
+      BOOST_CHECK_SMALL(std::abs(expected - result[i]), 5. * epsilon<double>);
     }
   }
   // LE pos offset and distance > box in shear plane normal direction
   {
-    auto const a = Vector3d{1.1, 12.2, -13};
-    auto const b = Vector3d{-11, 8.8, -13};
+    auto const a = Vector3d{1.1, 12.2, -13.};
+    auto const b = Vector3d{-11., 8.8, -13.};
     auto const le_jumps = std::round((a[0] - b[0]) / box.length()[0]);
 
     auto const result = box.get_mi_vector(a, b);
@@ -188,7 +191,7 @@ BOOST_AUTO_TEST_CASE(lees_edwards_mi_vector) {
             box.length()[2] // Manually apply minimum image convention
     }};
     for (int i : {0, 1, 2}) {
-      BOOST_CHECK_CLOSE(expected[i], result[i], epsilon<double>);
+      BOOST_CHECK_CLOSE(expected[i], result[i], 100. * epsilon<double>);
     }
   }
 
@@ -197,7 +200,7 @@ BOOST_AUTO_TEST_CASE(lees_edwards_mi_vector) {
   box.set_periodic(0, true);
   box.set_periodic(1, true);
   box.set_periodic(2, true);
-  box.set_length(Vector3d{5, 5, 5});
+  box.set_length(Vector3d{5., 5., 5.});
   box.set_lees_edwards_bc(LeesEdwardsBC{2.98, 0., 0, 1});
   auto const result =
       box.get_mi_vector(Vector3d{2.5, 1., 2.5}, Vector3d{4.52, 4., 2.5});

src/core/unit_tests/lees_edwards_test.cpp

@@ -48,24 +48,27 @@ BOOST_AUTO_TEST_CASE(test_shear_direction) {
 }
 
 BOOST_AUTO_TEST_CASE(test_update_offset) {
+  auto const prefactor = 2.5;
+  auto const old_offset = 1.5;
   Particle p;
   p.image_box() = {2, 4, -1};
-  p.lees_edwards_offset() = 1.5;
+  p.lees_edwards_offset() = old_offset;
+  p.lees_edwards_flag() = 1;
 
   LeesEdwardsBC le{0., 3.5, 1, 2};
   BoxGeometry box;
   box.set_lees_edwards_bc(le);
-  UpdateOffset(box, 3.5)(p);
-  auto const expected_offset = 1.5 - le.shear_velocity * 0.5 * 3.5 *
-                                         (p.image_box()[le.shear_plane_normal]);
+  UpdateOffset{box}(p, prefactor);
+  // Disabled as long as we do not use a two step LE update
+  auto const expected_offset = old_offset; // - prefactor * le.pos_offset / 2.
   BOOST_CHECK_CLOSE(p.lees_edwards_offset(), expected_offset, tol);
 }
 
 BOOST_AUTO_TEST_CASE(test_push) {
   auto const old_offset = -1.2;
   auto const shear_l = 6.;
   auto const shear_normal_l = 4.5;
-  auto const dt = 0.8;
+  auto const prefactor = 1.5;
   auto const old_pos = Utils::Vector3d{{3., shear_normal_l * 1.1, 10.}};
   auto const old_vel = Utils::Vector3d{{-1.2, 2., 4.1}};
 
@@ -85,23 +88,23 @@ BOOST_AUTO_TEST_CASE(test_push) {
   box.set_lees_edwards_bc(le);
 
   // Test transition in one direction
-  Push(box, dt)(p);
-  auto expected_pos = old_pos - shear_direction(box) * le.pos_offset;
+  Push{box}(p, prefactor);
+  auto expected_pos =
+      old_pos - prefactor * shear_direction(box) * le.pos_offset;
   auto expected_vel = old_vel - shear_direction(box) * le.shear_velocity;
-  auto expected_offset =
-      old_offset + le.pos_offset -
-      le.shear_velocity * 0.5 * dt * p.image_box()[le.shear_plane_normal];
+  auto expected_offset = old_offset + prefactor * le.pos_offset;
+  fold_position(expected_pos, p.image_box(), box);
   BOOST_CHECK_SMALL((p.pos() - expected_pos).norm(), eps);
   BOOST_CHECK_SMALL((p.v() - expected_vel).norm(), eps);
   BOOST_CHECK_CLOSE(p.lees_edwards_offset(), expected_offset, tol);
 
   // Test transition in the other direction
   p.pos()[le.shear_plane_normal] = -1;
-  Push(box, dt)(p);
+  Push{box}(p, prefactor);
   expected_pos = {old_pos[0], -1., old_pos[2]};
+  fold_position(expected_pos, p.image_box(), box);
   expected_vel = old_vel;
-  expected_offset = old_offset - le.shear_velocity * dt *
-                                     p.image_box()[le.shear_plane_normal];
+  expected_offset = old_offset;
   BOOST_CHECK_SMALL((p.pos() - expected_pos).norm(), eps);
   BOOST_CHECK_SMALL((p.v() - expected_vel).norm(), eps);
   BOOST_CHECK_CLOSE(p.lees_edwards_offset(), expected_offset, tol);

src/core/virtual_sites/VirtualSitesRelative.cpp

@@ -34,6 +34,9 @@
 #include <utils/math/tensor_product.hpp>
 #include <utils/quaternion.hpp>
 
+#include "integrate.hpp"
+#include "lees_edwards/lees_edwards.hpp"
+
 /**
  * @brief Vector pointing from the real particle to the virtual site.
  *
@@ -112,26 +115,12 @@ void VirtualSitesRelative::update() const {
       continue;
 
     auto const &p_ref = *p_ref_ptr;
-    auto new_pos = p_ref.pos() + connection_vector(p_ref, p);
-    /* The shift has to respect periodic boundaries: if the reference
-     * particles is not in the same image box, we potentially avoid shifting
-     * to the other side of the box. */
-    auto shift = box_geo.get_mi_vector(new_pos, p.pos());
-    p.pos() += shift;
-    Utils::Vector3i image_shift{};
-    fold_position(shift, image_shift, box_geo);
-    p.image_box() = p_ref.image_box() - image_shift;
-
+    p.pos() = p_ref.pos() + connection_vector(p_ref, p);
     p.v() = velocity(p_ref, p);
 
     if (box_geo.type() == BoxType::LEES_EDWARDS) {
-      auto const &lebc = box_geo.lees_edwards_bc();
-      auto const shear_dir = lebc.shear_direction;
-      auto const shear_normal = lebc.shear_plane_normal;
-      auto const le_vel = lebc.shear_velocity;
-      Utils::Vector3i n_shifts{};
-      fold_position(new_pos, n_shifts, box_geo);
-      p.v()[shear_dir] -= n_shifts[shear_normal] * le_vel;
+      auto push = LeesEdwards::Push(box_geo);
+      push(p, -1);
     }
 
     if (have_quaternions())

testsuite/python/lees_edwards.py

@@ -34,6 +34,13 @@
 params_osc = {'initial_pos_offset': 0.1, 'time_0': -2.1, 'amplitude': 2.3,
               'omega': 2.51}
 lin_protocol = espressomd.lees_edwards.LinearShear(**params_lin)
+
+
+def get_lin_pos_offset(time, initial_pos_offset=None,
+                       time_0=None, shear_velocity=None):
+    return initial_pos_offset + (time - time_0) * shear_velocity
+
+
 osc_protocol = espressomd.lees_edwards.OscillatoryShear(**params_osc)
 off_protocol = espressomd.lees_edwards.Off()
 
@@ -284,29 +291,28 @@ def test_boundary_crossing_off(self):
 
     def test_trajectory_reconstruction(self):
         system = self.system
+        system.time = 3.4
 
-        protocol = espressomd.lees_edwards.LinearShear(
-            shear_velocity=1., initial_pos_offset=0.0, time_0=0.0)
         system.lees_edwards.set_boundary_conditions(
-            shear_direction="x", shear_plane_normal="y", protocol=protocol)
+            shear_direction="x", shear_plane_normal="y", protocol=lin_protocol)
 
         pos = system.box_l - 0.01
         vel = np.array([0, 1, 0])
         p = system.part.add(pos=pos, v=vel)
+        old_x = p.pos_folded[0]
 
         system.integrator.run(1)
+        new_x = p.pos[0]
 
         np.testing.assert_almost_equal(
-            p.lees_edwards_flag * 1.0 * system.time_step * 0.5,
-            p.lees_edwards_offset)
+            p.lees_edwards_offset, 
+            -(new_x - old_x))
         np.testing.assert_almost_equal(p.lees_edwards_flag, -1)
 
-        offset1 = p.lees_edwards_flag * 1.0 * system.time_step * 0.5
-
-        system.integrator.run(1)
-
+        system.integrator.run(1)  # no boundary crossing
         np.testing.assert_almost_equal(
-            offset1 - 1.0 * 0.5, p.lees_edwards_offset)
+            p.lees_edwards_offset, 
+            -(new_x - old_x))  # unchanged
         np.testing.assert_almost_equal(p.lees_edwards_flag, 0)
 
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
@@ -421,20 +427,33 @@ def test_virt_sites(self):
 
         # Construct pair of VS across normal boundary
         system.lees_edwards.protocol = None
-        p1 = system.part.add(pos=(2.5, 0.0, 2.5), rotation=[False] * 3, id=0)
+        p1 = system.part.add(pos=(2.5, 0.0, 2.5), rotation=[
+                             False] * 3, id=0, v=np.array((-1, 2, 3)))
         p2 = system.part.add(pos=(2.5, 1.0, 2.5))
         p2.vs_auto_relate_to(p1)
         p3 = system.part.add(pos=(2.5, 4.0, 2.5))
         p3.vs_auto_relate_to(p1)
-        system.integrator.run(1)
 
         system.lees_edwards.set_boundary_conditions(
             shear_direction="x", shear_plane_normal="y", protocol=lin_protocol)
-        system.integrator.run(1)
+        # Test position and velocity of VS with Le shift
+        old_p3_pos = p3.pos
+        expected_p3_pos = old_p3_pos - \
+            np.array((get_lin_pos_offset(system.time, **params_lin), 0, 0))
+        system.integrator.run(0, recalc_forces=True)
+        np.testing.assert_allclose(np.copy(p3.pos_folded), expected_p3_pos)
+        np.testing.assert_allclose(
+            p3.v, p1.v + np.array((params_lin["shear_velocity"], 0, 0)))
+
+        # Check distances
         np.testing.assert_allclose(
             np.copy(system.distance_vec(p3, p2)), [0, 2, 0], atol=tol)
+        np.testing.assert_allclose(
+            np.copy(system.distance_vec(p2, p3)), [0, -2, 0], atol=tol)
         np.testing.assert_allclose(
             np.copy(system.velocity_difference(p3, p2)), [0, 0, 0], atol=tol)
+        np.testing.assert_allclose(
+            np.copy(system.velocity_difference(p2, p3)), [0, 0, 0], atol=tol)
         system.integrator.run(0)
         np.testing.assert_allclose(
             np.copy(system.distance_vec(p3, p2)), [0, 2, 0], atol=tol)
@@ -526,7 +545,7 @@ def test_virt_sites_interaction(self):
         ["EXTERNAL_FORCES", "VIRTUAL_SITES_RELATIVE", "COLLISION_DETECTION"])
     def test_le_colldet(self):
         system = self.system
-        system.min_global_cut = 1.0
+        system.min_global_cut = 1.2
         system.time = 0
         protocol = espressomd.lees_edwards.LinearShear(
             shear_velocity=-1.0, initial_pos_offset=0.0)
@@ -615,7 +634,7 @@ def test_le_colldet(self):
     @utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE"])
     def test_le_breaking_bonds(self):
         system = self.system
-        system.min_global_cut = 1.0
+        system.min_global_cut = 1.2
         system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
         protocol = espressomd.lees_edwards.LinearShear(
             shear_velocity=-1.0, initial_pos_offset=0.0)