Implement Lees-Edwards boundary conditions (#4457)

Fixes #2051
Fixes #2083

doc/sphinx/advanced_methods.rst

@@ -119,14 +119,6 @@ The following limitations currently apply for the collision detection:
   between virtual sites
 
 
-.. _Lees-Edwards boundary conditions:
-
-Lees-Edwards boundary conditions
---------------------------------
-
-Lees-Edwards boundary conditions are not available in the current version of |es|.
-
-
 .. _Immersed Boundary Method for soft elastic objects:
 
 Immersed Boundary Method for soft elastic objects

doc/sphinx/running.rst

@@ -198,6 +198,7 @@ Notes:
 * For the instantaneous pressure, the same limitations of applicability hold as described in :ref:`Pressure`.
 * The particle forces :math:`F` include interactions as well as a friction (:math:`\gamma^0`) and noise term (:math:`\sqrt{k_B T \gamma^0 dt} \overline{\eta}`) analogous to the terms in the :ref:`Langevin thermostat`.
 * The particle forces are only calculated in step 5 and then reused in step 1 of the next iteration. See :ref:`Velocity Verlet Algorithm` for the implications of that.
+* The NpT algorithm doesn't support :ref:`Lees-Edwards boundary conditions`.
 
 .. _Steepest descent:
 

doc/sphinx/system_setup.rst

@@ -40,6 +40,7 @@ for further calculations, you should explicitly make a copy e.g. via
   dimension non-periodic does not hinder particles from leaving the box in
   this direction. In this case for keeping particles in the simulation box
   a constraint has to be set.
+  For more details, see :ref:`Boundary conditions`.
 
 * :py:attr:`~espressomd.system.System.time_step`
 
@@ -81,6 +82,131 @@ or by calling the corresponding ``get_state()`` methods like::
     gamma = system.thermostat.get_state()[0]['gamma']
     gamma_rot = system.thermostat.get_state()[0]['gamma_rotation']
 
+.. _Boundary conditions:
+
+Boundary conditions
+-------------------
+
+.. _Periodic boundaries:
+
+Periodic boundaries
+~~~~~~~~~~~~~~~~~~~
+
+With periodic boundary conditions, particles interact with periodic
+images of all particles in the system. This is the default behavior.
+When particles cross a box boundary, their position are folded and
+their image box counter are incremented.
+
+From the Python interface, the folded position is accessed with
+:attr:`~espressomd.particle_data.ParticleHandle.pos_folded` and the image
+box counter with :attr:`~espressomd.particle_data.ParticleHandle.image_box`.
+Note that :attr:`~espressomd.particle_data.ParticleHandle.pos` gives the
+unfolded particle position.
+
+Example::
+
+    import espressomd
+    system = espressomd.System(box_l=[5.0, 5.0, 5.0], periodicity=[True, True, True])
+    system.time_step = 0.1
+    system.cell_system.skin = 0.0
+    p = system.part.add(pos=[4.9, 0.0, 0.0], v=[0.1, 0.0, 0.0])
+    system.integrator.run(20)
+    print(f"pos        = {p.pos}")
+    print(f"pos_folded = {p.pos_folded}")
+    print(f"image_box  = {p.image_box}")
+
+Output:
+
+.. code-block:: none
+
+    pos        = [5.1 0.  0. ]
+    pos_folded = [0.1 0.  0. ]
+    image_box  = [1 0 0]
+
+.. _Open boundaries:
+
+Open boundaries
+~~~~~~~~~~~~~~~
+
+With open boundaries, particles can leave the simulation box.
+What happens in this case depends on which algorithm is used.
+Some algorithms may require open boundaries,
+such as :ref:`Stokesian Dynamics`.
+
+Example::
+
+    import espressomd
+    system = espressomd.System(box_l=[5.0, 5.0, 5.0], periodicity=[False, False, False])
+    system.time_step = 0.1
+    system.cell_system.skin = 0.0
+    p = system.part.add(pos=[4.9, 0.0, 0.0], v=[0.1, 0.0, 0.0])
+    system.integrator.run(20)
+    print(f"pos        = {p.pos}")
+    print(f"pos_folded = {p.pos_folded}")
+    print(f"image_box  = {p.image_box}")
+
+Output:
+
+.. code-block:: none
+
+    pos        = [5.1 0.  0. ]
+    pos_folded = [5.1 0.  0. ]
+    image_box  = [0 0 0]
+
+.. _Lees-Edwards boundary conditions:
+
+Lees--Edwards boundary conditions
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Lees--Edwards boundary conditions (LEbc) are special periodic boundary
+conditions to simulation systems under shear stress :cite:`lees72a`.
+Periodic images of particles across the shear boundary appear with a
+time-dependent position offset. When a particle crosses the shear boundary,
+it appears to the opposite side of the simulation box with a position offset
+and a shear velocity :cite:`bindgen21a`.
+
+LEbc require a fully periodic system and are configured with
+:class:`~espressomd.lees_edwards.LinearShear` and
+:class:`~espressomd.lees_edwards.OscillatoryShear`.
+To temporarily disable LEbc, use :class:`~espressomd.lees_edwards.Off`.
+To completely disable LEbc and reinitialize the box geometry, do
+``system.lees_edwards.protocol = None``.
+
+Example::
+
+    import espressomd
+    import espressomd.lees_edwards
+    system = espressomd.System(box_l=[5.0, 5.0, 5.0])
+    system.time_step = 0.1
+    system.cell_system.skin = 0.0
+    system.cell_system.set_n_square(use_verlet_lists=True)
+    le_protocol = espressomd.lees_edwards.LinearShear(
+        shear_velocity=-0.1, initial_pos_offset=0.0, time_0=-0.1)
+    system.lees_edwards.protocol = le_protocol
+    system.lees_edwards.shear_direction = 1 # shear along y-axis
+    system.lees_edwards.shear_plane_normal = 0 # shear when crossing the x-boundary
+    p = system.part.add(pos=[4.9, 0.0, 0.0], v=[0.1, 0.0, 0.0])
+    system.integrator.run(20)
+    print(f"pos        = {p.pos}")
+    print(f"pos_folded = {p.pos_folded}")
+    print(f"image_box  = {p.image_box}")
+    print(f"velocity   = {p.v}")
+
+Output:
+
+.. code-block:: none
+
+    pos        = [5.1 0.2 0. ]
+    pos_folded = [0.1 0.2 0. ]
+    image_box  = [1 0 0]
+    velocity   = [0.1 0.1 0. ]
+
+Particles inserted outside the box boundaries will be wrapped around
+using the normal periodic boundary rules, i.e. they will not be sheared,
+even though their :attr:`~espressomd.particle_data.ParticleHandle.image_box`
+is *not* zero.
+
+
 .. _Cellsystems:
 
 Cellsystems

doc/sphinx/zrefs.bib

@@ -67,6 +67,17 @@ @ARTICLE{ballenegger09a
   publisher = {AIP},
 }
 
+@Article{bindgen21a,
+  author = {Bindgen, Sebastian and Weik, Florian and Weeber, Rudolf and Koos, Erin and de Buyl, Pierre},
+  title = {Lees-Edwards boundary conditions for translation invariant shear flow: implementation and transport properties},
+  journal = {Physics of Fluids},
+  year = {2021},
+  volume = {33},
+  number = {8},
+  doi = {10.1063/5.0055396},
+  pages = {083615},
+}
+
 @ARTICLE{brodka04a,
   author = {Br\'{o}dka, A.},
   title = {{E}wald summation method with electrostatic layer correction for interactions
@@ -296,6 +307,16 @@ @ARTICLE{humphrey96a
   doi = {10.1016/0263-7855(96)00018-5},
 }
 
+@Article{lees72a,
+  author    = {Lees, A. W. and Edwards, S. F.},
+  title     = {The computer study of transport processes under extreme conditions},
+  journal   = {Journal of Physics C: Solid State Physics},
+  year      = {1972},
+  volume    = {5},
+  number    = {15},
+  doi       = {10.1088/0022-3719/5/15/006},
+}
+
 @ARTICLE{tyagi08a,
   author = {Sandeep Tyagi and Axel Arnold and Christian Holm},
   title = {Electrostatic layer correction with image charges: A linear scaling

src/core/AtomDecomposition.hpp

@@ -101,6 +101,8 @@ class AtomDecomposition : public ParticleDecomposition {
     return m_box;
   }
 
+  BoxGeometry box() const override { return m_box; };
+
 private:
   /**
    * @brief Find cell for id.

src/core/BoxGeometry.hpp

@@ -19,9 +19,11 @@
 #ifndef CORE_BOX_GEOMETRY_HPP
 #define CORE_BOX_GEOMETRY_HPP
 
+#include "LeesEdwardsBC.hpp"
 #include "algorithm/periodic_fold.hpp"
 
 #include <utils/Vector.hpp>
+#include <utils/math/sgn.hpp>
 
 #include <bitset>
 #include <cassert>
@@ -68,8 +70,28 @@ template <typename T> T get_mi_coord(T a, T b, T box_length, bool periodic) {
 }
 } // namespace detail
 
+enum class BoxType { CUBOID = 0, LEES_EDWARDS = 1 };
+
 class BoxGeometry {
+public:
+  BoxGeometry() {
+    set_length(Utils::Vector3d{1., 1., 1.});
+    set_periodic(0, true);
+    set_periodic(1, true);
+    set_periodic(2, true);
+    set_type(BoxType::CUBOID);
+  }
+  BoxGeometry(BoxGeometry const &rhs) {
+    m_type = rhs.type();
+    set_length(rhs.length());
+    set_periodic(0, rhs.periodic(0));
+    set_periodic(1, rhs.periodic(1));
+    set_periodic(2, rhs.periodic(2));
+    m_lees_edwards_bc = rhs.m_lees_edwards_bc;
+  }
+
 private:
+  BoxType m_type = BoxType::CUBOID;
   /** Flags for all three dimensions whether pbc are applied (default). */
   std::bitset<3> m_periodic = 0b111;
   /** Side lengths of the box */
@@ -79,6 +101,9 @@ class BoxGeometry {
   /** Half side lengths of the box */
   Utils::Vector3d m_length_half = {0.5, 0.5, 0.5};
 
+  /** Lees-Edwards boundary conditions */
+  LeesEdwardsBC m_lees_edwards_bc;
+
 public:
   /**
    * @brief Set periodicity for direction
@@ -161,9 +186,37 @@ class BoxGeometry {
   template <typename T>
   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 clees_edwards_bc().distance(a - b, length(), length_half(),
+                                         length_inv(), m_periodic);
+    }
+    assert(type() == BoxType::CUBOID);
     return {get_mi_coord(a[0], b[0], 0), get_mi_coord(a[1], b[1], 1),
             get_mi_coord(a[2], b[2], 2)};
   }
+
+  BoxType type() const { return m_type; }
+  void set_type(BoxType type) { m_type = type; }
+
+  LeesEdwardsBC &lees_edwards_bc() { return m_lees_edwards_bc; }
+  LeesEdwardsBC const &clees_edwards_bc() const { return m_lees_edwards_bc; }
+  void set_lees_edwards_bc(LeesEdwardsBC bc) { m_lees_edwards_bc = bc; }
+
+  /** Calculate the velocity difference including the Lees-Edwards velocity */
+  Utils::Vector3d velocity_difference(Utils::Vector3d const &x,
+                                      Utils::Vector3d const &y,
+                                      Utils::Vector3d const &u,
+                                      Utils::Vector3d const &v) const {
+    auto ret = u - v;
+    if (type() == BoxType::LEES_EDWARDS) {
+      auto const &le = m_lees_edwards_bc;
+      auto const dy = x[le.shear_plane_normal] - y[le.shear_plane_normal];
+      if (fabs(dy) > 0.5 * length_half()[le.shear_plane_normal]) {
+        ret[le.shear_direction] -= Utils::sgn(dy) * le.shear_velocity;
+      }
+    }
+    return ret;
+  }
 };
 
 /** @brief Fold a coordinate to primary simulation box.

src/core/CellStructure.cpp

@@ -24,6 +24,8 @@
 #include "AtomDecomposition.hpp"
 #include "CellStructureType.hpp"
 #include "RegularDecomposition.hpp"
+#include "grid.hpp"
+#include "lees_edwards.hpp"
 
 #include <utils/contains.hpp>
 
@@ -221,7 +223,7 @@ struct UpdateParticleIndexVisitor {
 };
 } // namespace
 
-void CellStructure::resort_particles(int global_flag) {
+void CellStructure::resort_particles(int global_flag, BoxGeometry const &box) {
   invalidate_ghosts();
 
   static std::vector<ParticleChange> diff;
@@ -234,6 +236,7 @@ void CellStructure::resort_particles(int global_flag) {
   }
 
   m_rebuild_verlet_list = true;
+  m_le_pos_offset_at_last_resort = box.clees_edwards_bc().pos_offset;
 
 #ifdef ADDITIONAL_CHECKS
   check_particle_index();