WIP: Simplify VS_Relative, fix VS support for Lees Edwards

src/core/BoxGeometry.hpp

@@ -187,7 +187,20 @@ 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;
+      //      if (a[shear_plane_normal] <0 or a[shear_plane_normal]
+      //      >=length()[shear_plane_normal]
+      //          or b[shear_plane_normal] <0 or b[shear_plane_normal]
+      //          >=length()[shear_plane_normal]) throw
+      //          std::runtime_error("Unfolded position in shear plane normal
+      //          direction in LE get_mi_vector will lead to wrong result");
+      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

@@ -136,7 +136,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

@@ -26,37 +26,36 @@
 #include <cmath>
 #include <memory>
 
+#include <iostream>
 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 {
+    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,12 +65,10 @@ 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;
-    // TODO: clarify whether the fold is needed
-    //  p.pos()[m_le.shear_direction] = Algorithm::periodic_fold(
-    //      p.pos()[m_le.shear_direction], m_box_l()[m_le.shear_direction]);
-    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 / 2;
+    fold_position(p.pos(), p.image_box(), m_box);
+    //    UpdateOffset::operator()(p,pos_prefactor);
   }
 };
 

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,23 @@ 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);
+    //    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.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,14 @@
 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)
+
+
+# pass in **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()
 
@@ -294,11 +302,13 @@ def test_trajectory_reconstruction(self):
         vel = np.array([0, 1, 0])
         p = system.part.add(pos=pos, v=vel)
 
+        print(system.time, system.lees_edwards.pos_offset)      
         system.integrator.run(1)
+        print(system.time, system.lees_edwards.pos_offset)      
 
         np.testing.assert_almost_equal(
-            p.lees_edwards_flag * 1.0 * system.time_step * 0.5,
-            p.lees_edwards_offset)
+            p.lees_edwards_offset, 
+            get_lin_pos_offset(system.time - system.time_step / 2, **params_lin))
         np.testing.assert_almost_equal(p.lees_edwards_flag, -1)
 
         offset1 = p.lees_edwards_flag * 1.0 * system.time_step * 0.5
@@ -421,16 +431,26 @@ 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.integrator.run(0, recalc_forces=True)
 
         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(