From 8a7fdcf64c9152fc91e1d847353cc20bed3c3390 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 5 Sep 2022 16:34:51 +0200 Subject: [PATCH] Simplify VS_Relative, fix VS support for Lees Edwards --- src/core/BoxGeometry.hpp | 15 +++++- src/core/integrate.cpp | 2 +- src/core/lees_edwards/LeesEdwardsBC.hpp | 10 ++-- src/core/lees_edwards/lees_edwards.hpp | 33 ++++++------ .../virtual_sites/VirtualSitesRelative.cpp | 34 ++++++------- testsuite/python/lees_edwards.py | 51 +++++++++++++------ 6 files changed, 90 insertions(+), 55 deletions(-) diff --git a/src/core/BoxGeometry.hpp b/src/core/BoxGeometry.hpp index 7fdc9c13a64..9970ae27207 100644 --- a/src/core/BoxGeometry.hpp +++ b/src/core/BoxGeometry.hpp @@ -187,7 +187,20 @@ class BoxGeometry { Utils::Vector get_mi_vector(const Utils::Vector &a, const Utils::Vector &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); diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index a424f7c24d4..a0a149549d0 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -136,7 +136,7 @@ void unset_protocol() { template 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); }); diff --git a/src/core/lees_edwards/LeesEdwardsBC.hpp b/src/core/lees_edwards/LeesEdwardsBC.hpp index 29beaceb16f..e78d2059857 100644 --- a/src/core/lees_edwards/LeesEdwardsBC.hpp +++ b/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; diff --git a/src/core/lees_edwards/lees_edwards.hpp b/src/core/lees_edwards/lees_edwards.hpp index 90e926a5f64..abadb0966d5 100644 --- a/src/core/lees_edwards/lees_edwards.hpp +++ b/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(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(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; + fold_position(p.pos(), p.image_box(), m_box); + // UpdateOffset::operator()(p,pos_prefactor); } }; diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index cdeb15a4f34..6f04dc00870 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -34,6 +34,9 @@ #include #include +#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()) diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index b635724284b..0faf7a0f0f4 100644 --- a/testsuite/python/lees_edwards.py +++ b/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() @@ -284,29 +292,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 +428,34 @@ 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( + 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 +547,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 +636,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)