From a1b5610963068620009c7a6aa1fe14961f7d84b8 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 12:20:00 +0100 Subject: [PATCH 01/16] core: virtual_sites: Default all the things --- src/core/virtual_sites/VirtualSites.hpp | 4 ++-- src/core/virtual_sites/VirtualSitesInertialessTracers.hpp | 2 -- src/core/virtual_sites/VirtualSitesOff.hpp | 5 +---- 3 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index 454b416af09..d36e6deb818 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -43,9 +43,9 @@ class VirtualSites { * Velocities are only updated if get_have_velocity() returns true. * @param recalc_positions Skip the recalculation of positions if false. */ - virtual void update(bool recalc_positions) const = 0; + virtual void update(bool recalc_positions) const {} /** Back-transfer forces (and torques) to non-virtual particles. */ - virtual void back_transfer_forces_and_torques() const = 0; + virtual void back_transfer_forces_and_torques() const {} /** @brief Called after force calculation (and before rattle/shake) */ virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp index 3eb8973a5a4..c22a66687f1 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.hpp @@ -28,8 +28,6 @@ * instantaneously transferred to the fluid */ class VirtualSitesInertialessTracers : public VirtualSites { - void update(bool recalc_positions) const override{}; - void back_transfer_forces_and_torques() const override{}; void after_force_calc() override; void after_lb_propagation() override; }; diff --git a/src/core/virtual_sites/VirtualSitesOff.hpp b/src/core/virtual_sites/VirtualSitesOff.hpp index ba4ddc8959b..75918ba616f 100644 --- a/src/core/virtual_sites/VirtualSitesOff.hpp +++ b/src/core/virtual_sites/VirtualSitesOff.hpp @@ -24,10 +24,7 @@ #include "VirtualSites.hpp" /** @brief Do-nothing virtual-sites implementation */ -class VirtualSitesOff : public VirtualSites { - void update(bool) const override{}; - void back_transfer_forces_and_torques() const override{}; -}; +class VirtualSitesOff : public VirtualSites {}; #endif #endif From 53914907078ffe02e5816b35a10d398af110f5d5 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 13:51:13 +0100 Subject: [PATCH 02/16] core: virtual_sites: Allways have pressure contribution --- src/core/pressure.cpp | 4 ++-- src/core/virtual_sites/VirtualSites.hpp | 2 -- src/core/virtual_sites/VirtualSitesRelative.cpp | 1 + src/core/virtual_sites/VirtualSitesRelative.hpp | 2 -- testsuite/python/virtual_sites_relative.py | 8 +++++--- 5 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 70b4c40e3a1..192ceca2739 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -190,7 +190,7 @@ void init_virials(Observable_stat *stat) { Dipole::pressure_n(); #endif #ifdef VIRTUAL_SITES - n_vs = virtual_sites()->n_pressure_contribs(); + n_vs = 1; #endif // Allocate memory for the data @@ -228,7 +228,7 @@ void init_p_tensor(Observable_stat *stat) { auto const n_dipolar = 0; #endif #ifdef VIRTUAL_SITES - n_vs = virtual_sites()->n_pressure_contribs(); + n_vs = 1; #endif obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded, diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index d36e6deb818..4164ede18d0 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -49,8 +49,6 @@ class VirtualSites { /** @brief Called after force calculation (and before rattle/shake) */ virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; - /** @brief Number of pressure contributions */ - virtual int n_pressure_contribs() const { return 0; }; /** @brief Pressure contribution. */ virtual void pressure_and_stress_tensor_contribution(double *pressure, diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 5dc59824aa8..a77931c567e 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -67,6 +67,7 @@ void VirtualSitesRelative::update_virtual_particle_quaternion( "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " "particle associated with virtual site.\n"); } + p.r.quat = multiply_quaternions(p_real->r.quat, p.p.vs_relative.quat); } diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 975cecee5de..7b6c1f018b8 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -35,8 +35,6 @@ class VirtualSitesRelative : public VirtualSites { void update(bool recalc_positions) const override; /** @copydoc VirtualSites::back_transfer_forces_and_torques */ void back_transfer_forces_and_torques() const override; - /** @copydoc VirtualSites::n_pressure_contribs */ - int n_pressure_contribs() const override { return 1; }; /** @copydoc VirtualSites::pressure_and_stress_tensor_contribution */ void pressure_and_stress_tensor_contribution(double *pressure, diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index b7a125c7ed7..68b4461c6e0 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -320,10 +320,12 @@ def test_zz_stress_tensor(self): system.time_step = 0.01 system.cell_system.skin = 0.1 system.min_global_cut = 0.2 - # Should not have one if vs are turned off + # Should not have a pressure system.virtual_sites = VirtualSitesOff() - self.assertNotIn("virtual_sites", system.analysis.pressure()) - self.assertNotIn("virtual_sites", system.analysis.stress_tensor()) + stress_vs = system.analysis.stress_tensor()["virtual_sites", 0] + p_vs = system.analysis.pressure()["virtual_sites", 0] + np.testing.assert_allclose(stress_vs, 0., atol=1e-10) + np.testing.assert_allclose(p_vs, 0., atol=1e-10) # vs relative contrib system.virtual_sites = VirtualSitesRelative() From c445ab93b3659cac91f4924aa328fa7f109d5843 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 14:18:29 +0100 Subject: [PATCH 03/16] utils: trace(Matrix) + test --- src/utils/include/utils/Vector.hpp | 19 +++++++++++++++++++ src/utils/tests/Vector_test.cpp | 8 ++++++++ 2 files changed, 27 insertions(+) diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index 8e0467b09a9..5b819e9c8ea 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -155,6 +155,25 @@ using Vector3i = VectorXi<3>; template using Matrix = Vector, N>; +/** + * @brief Trace of a matrix. + * + * Returns the sum of the diagonal elements + * of a square matrix. + * + * @tparam T Arithmetic type + * @tparam N Matrix dimension + * @param m Input matrix + * @return Trace of matrix. + */ +template T trace(Matrix const &m) { + T tr = T{}; + for (size_t i = 0; i < N; i++) + tr += m[i][i]; + + return tr; +} + namespace detail { template auto binary_op(Vector const &a, Vector const &b, Op op) { diff --git a/src/utils/tests/Vector_test.cpp b/src/utils/tests/Vector_test.cpp index 421bad9f825..635c26344d2 100644 --- a/src/utils/tests/Vector_test.cpp +++ b/src/utils/tests/Vector_test.cpp @@ -381,3 +381,11 @@ BOOST_AUTO_TEST_CASE(diag_matrix) { for (int j = 0; j < 3; j++) BOOST_CHECK_EQUAL(result[i][j], (i == j) ? v[i] : 0); } + +BOOST_AUTO_TEST_CASE(trace_) { + auto const A = Utils::Matrix{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}; + auto const result = trace(A); + auto const expected = A[0][0] + A[1][1] + A[2][2]; + + BOOST_CHECK_EQUAL(expected, result); +} \ No newline at end of file From 4ef7de30eb12a48d406f7401d1f7512d119afeb2 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 15:19:58 +0100 Subject: [PATCH 04/16] utils: flatten + test --- src/utils/include/utils/Vector.hpp | 17 +++++++++++++++++ src/utils/tests/Vector_test.cpp | 10 +++++++++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index 5b819e9c8ea..91a46baff29 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -174,6 +174,23 @@ template T trace(Matrix const &m) { return tr; } +/** + * @brief Flatten a matrix to a linear vector. + * + * @param m Input Matrix + * @return Flat vector with elements of @ref m. + */ +template +Vector flatten(Matrix const &m) { + Vector ret; + + for (size_t i = 0; i < N; i++) + for (size_t j = 0; j < M; j++) + ret[i * M + j] = m[j][i]; + + return ret; +} + namespace detail { template auto binary_op(Vector const &a, Vector const &b, Op op) { diff --git a/src/utils/tests/Vector_test.cpp b/src/utils/tests/Vector_test.cpp index 635c26344d2..689faaff662 100644 --- a/src/utils/tests/Vector_test.cpp +++ b/src/utils/tests/Vector_test.cpp @@ -388,4 +388,12 @@ BOOST_AUTO_TEST_CASE(trace_) { auto const expected = A[0][0] + A[1][1] + A[2][2]; BOOST_CHECK_EQUAL(expected, result); -} \ No newline at end of file +} + +BOOST_AUTO_TEST_CASE(flatten_) { + auto const A = Utils::Matrix{{1, 2}, {3, 4}}; + auto const result = flatten(A); + auto const expected = Utils::Vector{1, 3, 2, 4}; + + BOOST_CHECK(result == expected); +} From 7c53ddfc0bd946888e5efc65f3d44018a47f13a4 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 15:20:48 +0100 Subject: [PATCH 05/16] core: virtual_sites: Return stress tensor from VS --- src/core/pressure.cpp | 10 ++++++++-- src/core/virtual_sites/VirtualSites.hpp | 9 ++++++--- src/core/virtual_sites/VirtualSitesRelative.cpp | 17 +++++++---------- src/core/virtual_sites/VirtualSitesRelative.hpp | 5 ++--- 4 files changed, 23 insertions(+), 18 deletions(-) diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 192ceca2739..c2ee42ca16e 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -29,6 +29,7 @@ #include "npt.hpp" #include "pressure_inline.hpp" #include "virtual_sites.hpp" +#include #include "short_range_loop.hpp" @@ -127,8 +128,13 @@ void pressure_calc(double *result, double *result_t, double *result_nb, calc_long_range_virials(cell_structure.local_cells().particles()); #ifdef VIRTUAL_SITES - virtual_sites()->pressure_and_stress_tensor_contribution( - virials.virtual_sites, p_tensor.virtual_sites); + { + auto const vs_stress = + virtual_sites()->pressure_and_stress_tensor_contribution(); + + *virials.virtual_sites += trace(vs_stress); + boost::copy(flatten(vs_stress), p_tensor.virtual_sites); + } #endif for (int n = 1; n < virials.data.n; n++) diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index 4164ede18d0..4797db842a9 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -33,6 +33,8 @@ */ #ifdef VIRTUAL_SITES +#include + #include /** @brief Base class for virtual sites implementations */ @@ -50,9 +52,10 @@ class VirtualSites { virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; /** @brief Pressure contribution. */ - virtual void - pressure_and_stress_tensor_contribution(double *pressure, - double *stress_tensor) const {}; + virtual Utils::Matrix + pressure_and_stress_tensor_contribution() const { + return {}; + }; /** @brief Enable/disable velocity calculations for vs. */ void set_have_velocity(const bool &v) { m_have_velocity = v; }; const bool &get_have_velocity() const { return m_have_velocity; }; diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index a77931c567e..3c4b1303bb1 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -158,11 +158,13 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { } // Rigid body contribution to scalar pressure and stress tensor -void VirtualSitesRelative::pressure_and_stress_tensor_contribution( - double *pressure, double *stress_tensor) const { +Utils::Matrix +VirtualSitesRelative::pressure_and_stress_tensor_contribution() const { // Division by 3 volume is somewhere else. (pressure.cpp after all pressure // calculations) Iterate over all the particles in the local cells + Utils::Matrix stress_tensor = {}; + for (auto &p : cell_structure.local_cells().particles()) { if (!p.p.is_virtual) continue; @@ -178,15 +180,10 @@ void VirtualSitesRelative::pressure_and_stress_tensor_contribution( // conditions auto const d = get_mi_vector(p_real->r.p, p.r.p, box_geo); - // Stress tensor contribution - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - stress_tensor[k * 3 + l] += p.f.f[k] * d[l]; - - // Pressure = 1/3 trace of stress tensor - // but the 1/3 is applied somewhere else. - *pressure += p.f.f * d; + stress_tensor += tensor_product(d, p.f.f); } + + return stress_tensor; } #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 7b6c1f018b8..263fb53bda3 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -36,9 +36,8 @@ class VirtualSitesRelative : public VirtualSites { /** @copydoc VirtualSites::back_transfer_forces_and_torques */ void back_transfer_forces_and_torques() const override; /** @copydoc VirtualSites::pressure_and_stress_tensor_contribution */ - void - pressure_and_stress_tensor_contribution(double *pressure, - double *stress_tensor) const override; + Utils::Matrix + pressure_and_stress_tensor_contribution() const override; private: /** Update the position of the given virtual particle as defined by the real From 8a61ce59e5ce6bebf0faa1a16f1d85f508d9dcd0 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 15:58:44 +0100 Subject: [PATCH 06/16] core: virtual_sites: Cleanup --- src/core/pressure_inline.hpp | 16 ++++++++-------- src/core/virtual_sites/VirtualSitesRelative.cpp | 6 +----- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 4d6a3d9c1ba..fb4a8ee0a58 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -73,17 +73,17 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, } #ifdef ELECTROSTATICS - /* real space Coulomb */ - auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist); + { + /* real space Coulomb */ + auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist); - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; + } } - } - for (int i = 0; i < 3; i++) { - virials.coulomb[0] += p_coulomb[i][i]; + virials.coulomb[0] += trace(p_coulomb); } #endif /*ifdef ELECTROSTATICS */ diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 3c4b1303bb1..06f7963175e 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -19,6 +19,7 @@ #include "VirtualSitesRelative.hpp" #include "forces_inline.hpp" +#include #include #include @@ -160,17 +161,12 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // Rigid body contribution to scalar pressure and stress tensor Utils::Matrix VirtualSitesRelative::pressure_and_stress_tensor_contribution() const { - // Division by 3 volume is somewhere else. (pressure.cpp after all pressure - // calculations) Iterate over all the particles in the local cells - Utils::Matrix stress_tensor = {}; for (auto &p : cell_structure.local_cells().particles()) { if (!p.p.is_virtual) continue; - update_pos(p); - // First obtain the real particle responsible for this virtual particle: const Particle *p_real = get_local_particle_data(p.p.vs_relative.to_particle_id); From a5ea4b456a3d0d34d64d015a10fbea4d55046bf7 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 16:24:27 +0100 Subject: [PATCH 07/16] core: ghosts: Skip ghost comm if there is no data to transmit --- src/core/ghosts.cpp | 3 +++ src/core/virtual_sites/VirtualSitesRelative.cpp | 4 +--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index e01fe2f0563..74f3fe72a4c 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -320,6 +320,9 @@ static bool is_poststorable(GhostCommunication const &ghost_comm) { } void ghost_communicator(GhostCommunicator *gcr, unsigned int data_parts) { + if (GHOSTTRANS_NONE == data_parts) + return; + static CommBuf send_buffer, recv_buffer; for (auto it = gcr->comm.begin(); it != gcr->comm.end(); ++it) { diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 06f7963175e..d0368f3b1a9 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -40,9 +40,7 @@ void VirtualSitesRelative::update(bool recalc_positions) const { (get_have_velocity() ? (GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM) : GHOSTTRANS_NONE); - if (recalc_positions or get_have_velocity()) { - ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); - } + ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); } for (auto &p : cell_structure.local_cells().particles()) { if (!p.p.is_virtual) From cc745223cb751c6f3184a2eaea08e93bfb0326ed Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 17:23:43 +0100 Subject: [PATCH 08/16] core: header cleanup --- src/core/virtual_sites.hpp | 2 ++ src/core/virtual_sites/VirtualSites.hpp | 2 -- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index 8aab54a63ad..3bb3dd719a7 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -26,6 +26,8 @@ #include "virtual_sites/VirtualSites.hpp" +#include + /** @brief get active virtual sites implementation */ const std::shared_ptr &virtual_sites(); diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index 4797db842a9..f855afc9467 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -35,8 +35,6 @@ #ifdef VIRTUAL_SITES #include -#include - /** @brief Base class for virtual sites implementations */ class VirtualSites { public: From 50d79059a46e2bc9f9a26f4294443d68e4438af7 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 17:30:08 +0100 Subject: [PATCH 09/16] core: virtual_sites: Better name for stress calculation --- src/core/pressure.cpp | 3 +-- src/core/virtual_sites/VirtualSites.hpp | 3 +-- src/core/virtual_sites/VirtualSitesRelative.cpp | 2 +- src/core/virtual_sites/VirtualSitesRelative.hpp | 5 ++--- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index c2ee42ca16e..1b63fc88c9b 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -129,8 +129,7 @@ void pressure_calc(double *result, double *result_t, double *result_nb, #ifdef VIRTUAL_SITES { - auto const vs_stress = - virtual_sites()->pressure_and_stress_tensor_contribution(); + auto const vs_stress = virtual_sites()->stress_tensor(); *virials.virtual_sites += trace(vs_stress); boost::copy(flatten(vs_stress), p_tensor.virtual_sites); diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index f855afc9467..bbe33d98db6 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -50,8 +50,7 @@ class VirtualSites { virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; /** @brief Pressure contribution. */ - virtual Utils::Matrix - pressure_and_stress_tensor_contribution() const { + virtual Utils::Matrix stress_tensor() const { return {}; }; /** @brief Enable/disable velocity calculations for vs. */ diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index d0368f3b1a9..945cc930d44 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -158,7 +158,7 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // Rigid body contribution to scalar pressure and stress tensor Utils::Matrix -VirtualSitesRelative::pressure_and_stress_tensor_contribution() const { +VirtualSitesRelative::stress_tensor() const { Utils::Matrix stress_tensor = {}; for (auto &p : cell_structure.local_cells().particles()) { diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 263fb53bda3..4b15043cb57 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -35,9 +35,8 @@ class VirtualSitesRelative : public VirtualSites { void update(bool recalc_positions) const override; /** @copydoc VirtualSites::back_transfer_forces_and_torques */ void back_transfer_forces_and_torques() const override; - /** @copydoc VirtualSites::pressure_and_stress_tensor_contribution */ - Utils::Matrix - pressure_and_stress_tensor_contribution() const override; + /** @copydoc VirtualSites::stress_tensor */ + Utils::Matrix stress_tensor() const override; private: /** Update the position of the given virtual particle as defined by the real From 46df169aa43315940e18778b7d371f7a842f88b2 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 18:55:52 +0100 Subject: [PATCH 10/16] core: Removed redundant if guards --- src/core/integrators/brownian_inline.hpp | 5 +---- src/core/integrators/velocity_verlet_inline.hpp | 3 +-- src/core/virtual_sites/VirtualSites.hpp | 4 +--- src/core/virtual_sites/VirtualSitesRelative.cpp | 3 +-- 4 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp index 902fe3c2505..dcf246398b4 100644 --- a/src/core/integrators/brownian_inline.hpp +++ b/src/core/integrators/brownian_inline.hpp @@ -507,11 +507,8 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, const ParticleRange &particles) { for (auto &p : particles) { // Don't propagate translational degrees of freedom of vs -#ifdef VIRTUAL_SITES extern bool thermo_virtual; - if (!(p.p.is_virtual) or thermo_virtual) -#endif - { + if (!(p.p.is_virtual) or thermo_virtual) { p.r.p += bd_drag(brownian.gamma, p, time_step); p.m.v = bd_drag_vel(brownian.gamma, p); p.r.p += bd_random_walk(brownian, p, time_step); diff --git a/src/core/integrators/velocity_verlet_inline.hpp b/src/core/integrators/velocity_verlet_inline.hpp index 7d51915f621..8dafc9a1e50 100644 --- a/src/core/integrators/velocity_verlet_inline.hpp +++ b/src/core/integrators/velocity_verlet_inline.hpp @@ -67,11 +67,10 @@ inline void velocity_verlet_propagate_vel_final(const ParticleRange &particles) { for (auto &p : particles) { -#ifdef VIRTUAL_SITES // Virtual sites are not propagated during integration if (p.p.is_virtual) continue; -#endif + for (int j = 0; j < 3; j++) { if (!(p.p.ext_flag & COORD_FIXED(j))) { /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */ diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index bbe33d98db6..96eef3fd50b 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -50,9 +50,7 @@ class VirtualSites { virtual void after_force_calc(){}; virtual void after_lb_propagation(){}; /** @brief Pressure contribution. */ - virtual Utils::Matrix stress_tensor() const { - return {}; - }; + virtual Utils::Matrix stress_tensor() const { return {}; }; /** @brief Enable/disable velocity calculations for vs. */ void set_have_velocity(const bool &v) { m_have_velocity = v; }; const bool &get_have_velocity() const { return m_have_velocity; }; diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 945cc930d44..d523bbbd926 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -157,8 +157,7 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { } // Rigid body contribution to scalar pressure and stress tensor -Utils::Matrix -VirtualSitesRelative::stress_tensor() const { +Utils::Matrix VirtualSitesRelative::stress_tensor() const { Utils::Matrix stress_tensor = {}; for (auto &p : cell_structure.local_cells().particles()) { From 9815c8ea64b236c0d4f1fe094934158e0c4a433f Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 19:39:06 +0100 Subject: [PATCH 11/16] core: virtual_sites: Refactor --- .../virtual_sites/VirtualSitesRelative.cpp | 148 ++++++++---------- .../virtual_sites/VirtualSitesRelative.hpp | 14 -- 2 files changed, 68 insertions(+), 94 deletions(-) diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index d523bbbd926..31f2e05e8a3 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -19,7 +19,6 @@ #include "VirtualSitesRelative.hpp" #include "forces_inline.hpp" -#include #include #include @@ -32,6 +31,58 @@ #include "integrate.hpp" #include "rotation.hpp" +namespace { +Utils::Vector4d updated_quaternion( + Particle const *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return multiply_quaternions(p_real->r.quat, vs_rel.quat); +} + +Utils::Vector3d +updated_pos(Particle const *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + // Calculate the quaternion defining the orientation of the vector connecting + // the virtual site and the real particle + // This is obtained, by multiplying the quaternion representing the director + // of the real particle with the quaternion of the virtual particle, which + // specifies the relative orientation. + auto const director = + Utils::convert_quaternion_to_director( + Utils::multiply_quaternions(p_real->r.quat, vs_rel.rel_orientation)) + .normalize(); + + return p_real->r.p + director * vs_rel.distance; +} + +Utils::Vector3d +updated_vel(const Particle *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + auto const d = + Utils::convert_quaternion_to_director( + Utils::multiply_quaternions(p_real->r.quat, vs_rel.rel_orientation)) + .normalize() * + vs_rel.distance; + + // Get omega of real particle in space-fixed frame + auto const omega_space_frame = + convert_vector_body_to_space(*p_real, p_real->m.omega); + // Obtain velocity from v=v_real particle + omega_real_particle \times + // director + return vector_product(omega_space_frame, d) + p_real->m.v; +} + +Particle *get_real_particle( + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + auto p_real = get_local_particle_data(vs_rel.to_particle_id); + if (!p_real) { + throw std::runtime_error("No real particle associated with virtual site."); + } + + return p_real; +} + +} // namespace + void VirtualSitesRelative::update(bool recalc_positions) const { // Ghost update logic if (n_nodes > 1) { @@ -46,85 +97,25 @@ void VirtualSitesRelative::update(bool recalc_positions) const { if (!p.p.is_virtual) continue; - if (recalc_positions) - update_pos(p); + const Particle *p_real = get_real_particle(p.p.vs_relative); - if (get_have_velocity()) - update_vel(p); + if (recalc_positions) { + auto const new_pos = updated_pos(p_real, p.p.vs_relative); + /* The shift has to respect periodic boundaries: if the reference + * particles is not in the same image box, we potentially avoid to shift + * to the other side of the box. */ + p.r.p += get_mi_vector(new_pos, p.r.p, box_geo); - if (get_have_quaternion()) - update_virtual_particle_quaternion(p); - } -} - -void VirtualSitesRelative::update_virtual_particle_quaternion( - Particle &p) const { - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - if (!p_real) { - throw std::runtime_error( - "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"); - } - - p.r.quat = multiply_quaternions(p_real->r.quat, p.p.vs_relative.quat); -} - -void VirtualSitesRelative::update_pos(Particle &p) const { - // First obtain the real particle responsible for this virtual particle: - // Find the 1st real particle in the topology for the virtual particle's - // mol_id - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); - // Check, if a real particle was found - if (!p_real) { - runtimeErrorMsg() - << "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"; - return; - } - - // Calculate the quaternion defining the orientation of the vector connecting - // the virtual site and the real particle - // This is obtained, by multiplying the quaternion representing the director - // of the real particle with the quaternion of the virtual particle, which - // specifies the relative orientation. - auto const director = - Utils::convert_quaternion_to_director( - Utils::multiply_quaternions(p_real->r.quat, - p.p.vs_relative.rel_orientation)) - .normalize(); - - auto const new_pos = p_real->r.p + director * p.p.vs_relative.distance; - /* The shift has to respect periodic boundaries: if the reference particles - * is not in the same image box, we potentially avoid to shift to the other - * side of the box. */ - auto const shift = get_mi_vector(new_pos, p.r.p, box_geo); - p.r.p += shift; + if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) + set_resort_particles(Cells::RESORT_LOCAL); + } - if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) - set_resort_particles(Cells::RESORT_LOCAL); -} + if (get_have_velocity()) + p.m.v = updated_vel(p_real, p.p.vs_relative); -void VirtualSitesRelative::update_vel(Particle &p) const { - // First obtain the real particle responsible for this virtual particle: - Particle *p_real = get_local_particle_data(p.p.vs_relative.to_particle_id); - // Check, if a real particle was found - if (!p_real) { - runtimeErrorMsg() - << "virtual_sites_relative.cpp - update_mol_pos_particle(): No real " - "particle associated with virtual site.\n"; - return; + if (get_have_quaternion()) + p.r.quat = updated_quaternion(p_real, p.p.vs_relative); } - - auto const d = get_mi_vector(p.r.p, p_real->r.p, box_geo); - - // Get omega of real particle in space-fixed frame - Utils::Vector3d omega_space_frame = - convert_vector_body_to_space(*p_real, p_real->m.omega); - // Obtain velocity from v=v_real particle + omega_real_particle \times - // director - p.m.v = vector_product(omega_space_frame, d) + p_real->m.v; } // Distribute forces that have accumulated on virtual particles to the @@ -139,8 +130,7 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // We only care about virtual particles if (p.p.is_virtual) { // First obtain the real particle responsible for this virtual particle: - Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); + Particle *p_real = get_real_particle(p.p.vs_relative); // The rules for transferring forces are: // F_realParticle +=F_virtualParticle @@ -165,8 +155,7 @@ Utils::Matrix VirtualSitesRelative::stress_tensor() const { continue; // First obtain the real particle responsible for this virtual particle: - const Particle *p_real = - get_local_particle_data(p.p.vs_relative.to_particle_id); + const Particle *p_real = get_real_particle(p.p.vs_relative); // Get distance vector pointing from real to virtual particle, respecting // periodic boundary i @@ -178,5 +167,4 @@ Utils::Matrix VirtualSitesRelative::stress_tensor() const { return stress_tensor; } - #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 4b15043cb57..1f2f5b02cec 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -37,20 +37,6 @@ class VirtualSitesRelative : public VirtualSites { void back_transfer_forces_and_torques() const override; /** @copydoc VirtualSites::stress_tensor */ Utils::Matrix stress_tensor() const override; - -private: - /** Update the position of the given virtual particle as defined by the real - * particles in the same molecule - */ - void update_pos(Particle &p) const; - /** Update the velocity of the given virtual particle as defined by the real - * particles in the same molecule - */ - void update_vel(Particle &p) const; - /** @brief Update the orientation of the virtual particles with respect - * to the real particle. - */ - void update_virtual_particle_quaternion(Particle &p) const; }; #endif From 1f6f2de776dba86966792980fc39b45597a15e4e Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Tue, 11 Feb 2020 20:09:32 +0100 Subject: [PATCH 12/16] core: virtual_sites: Better names and docstrings --- .../virtual_sites/VirtualSitesRelative.cpp | 87 +++++++++++++------ 1 file changed, 62 insertions(+), 25 deletions(-) diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 31f2e05e8a3..0b425298db4 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -32,15 +32,27 @@ #include "rotation.hpp" namespace { -Utils::Vector4d updated_quaternion( - Particle const *p_real, - const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { +/** + * @brief Orientation of the virtual site. + * @param p_real Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Orientation quaternion of the virtual site. + */ +Utils::Vector4d +orientation(Particle const *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { return multiply_quaternions(p_real->r.quat, vs_rel.quat); } -Utils::Vector3d -updated_pos(Particle const *p_real, - const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { +/** + * @brief Vector pointing from the real particle + * to the virtual site. + * + * @return Relative distance. + */ +Utils::Vector3d connection_vector( + Particle const *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { // Calculate the quaternion defining the orientation of the vector connecting // the virtual site and the real particle // This is obtained, by multiplying the quaternion representing the director @@ -51,17 +63,31 @@ updated_pos(Particle const *p_real, Utils::multiply_quaternions(p_real->r.quat, vs_rel.rel_orientation)) .normalize(); - return p_real->r.p + director * vs_rel.distance; + return vs_rel.distance * director; } +/** + * @brief Position of the virtual site + * @param p_real Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Position of the virtual site. + */ Utils::Vector3d -updated_vel(const Particle *p_real, - const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { - auto const d = - Utils::convert_quaternion_to_director( - Utils::multiply_quaternions(p_real->r.quat, vs_rel.rel_orientation)) - .normalize() * - vs_rel.distance; +position(Particle const *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return p_real->r.p + connection_vector(p_real, vs_rel); +} + +/** + * @brief Velocity of the virtual site + * @param p_real Reference particle for the virtual site. + * @param vs_rel Parameters for the virtual site. + * @return Velocity of the virtual site. + */ +Utils::Vector3d +velocity(const Particle *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + auto const d = connection_vector(p_real, vs_rel); // Get omega of real particle in space-fixed frame auto const omega_space_frame = @@ -81,6 +107,22 @@ Particle *get_real_particle( return p_real; } +/** + * @brief Constraint force to hold the particles at its prescribed position. + * + * @param f Force on the virtual site. + * @param p_real Reference particle. + * @param vs_real Parameters. + * @return Constraint force. + */ +auto constraint_stress( + const Utils::Vector3d &f, const Particle *p_real, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + /* The constraint force is minus the force on the particle, make it force + * free. The counter force is translated by the connection vector to the + * real particle, hence the virial stress is */ + return tensor_product(connection_vector(p_real, vs_rel), -f); +} } // namespace void VirtualSitesRelative::update(bool recalc_positions) const { @@ -100,7 +142,7 @@ void VirtualSitesRelative::update(bool recalc_positions) const { const Particle *p_real = get_real_particle(p.p.vs_relative); if (recalc_positions) { - auto const new_pos = updated_pos(p_real, p.p.vs_relative); + auto const new_pos = position(p_real, p.p.vs_relative); /* The shift has to respect periodic boundaries: if the reference * particles is not in the same image box, we potentially avoid to shift * to the other side of the box. */ @@ -111,11 +153,11 @@ void VirtualSitesRelative::update(bool recalc_positions) const { } if (get_have_velocity()) - p.m.v = updated_vel(p_real, p.p.vs_relative); + p.m.v = velocity(p_real, p.p.vs_relative); if (get_have_quaternion()) - p.r.quat = updated_quaternion(p_real, p.p.vs_relative); - } + p.r.quat = orientation(p_real, p.p.vs_relative); + } // namespace } // Distribute forces that have accumulated on virtual particles to the @@ -139,7 +181,7 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // Add forces and torques p_real->f.torque += - vector_product(get_mi_vector(p.r.p, p_real->r.p, box_geo), p.f.f) + + vector_product(connection_vector(p_real, p.p.vs_relative), p.f.f) + p.f.torque; p_real->f.f += p.f.f; } @@ -157,12 +199,7 @@ Utils::Matrix VirtualSitesRelative::stress_tensor() const { // First obtain the real particle responsible for this virtual particle: const Particle *p_real = get_real_particle(p.p.vs_relative); - // Get distance vector pointing from real to virtual particle, respecting - // periodic boundary i - // conditions - auto const d = get_mi_vector(p_real->r.p, p.r.p, box_geo); - - stress_tensor += tensor_product(d, p.f.f); + stress_tensor += constraint_stress(p.f.f, p_real, p.p.vs_relative); } return stress_tensor; From da62e07c54f320593b3bb622a5b2c1e45216de25 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 4 Mar 2020 16:07:48 +0100 Subject: [PATCH 13/16] core: virtual_sites: Better names --- .../virtual_sites/VirtualSitesRelative.cpp | 106 ++++++++++-------- .../virtual_sites/VirtualSitesRelative.hpp | 6 +- 2 files changed, 61 insertions(+), 51 deletions(-) diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 0b425298db4..6d666730441 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -18,30 +18,30 @@ */ #include "VirtualSitesRelative.hpp" -#include "forces_inline.hpp" -#include -#include #ifdef VIRTUAL_SITES_RELATIVE #include "cells.hpp" -#include "config.hpp" -#include "errorhandling.hpp" +#include "forces_inline.hpp" #include "grid.hpp" #include "integrate.hpp" +#include "particle_data.hpp" #include "rotation.hpp" +#include +#include + namespace { /** * @brief Orientation of the virtual site. - * @param p_real Reference particle for the virtual site. + * @param p_ref Reference particle for the virtual site. * @param vs_rel Parameters for the virtual site. * @return Orientation quaternion of the virtual site. */ Utils::Vector4d -orientation(Particle const *p_real, +orientation(Particle const *p_ref, const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { - return multiply_quaternions(p_real->r.quat, vs_rel.quat); + return multiply_quaternions(p_ref->r.quat, vs_rel.quat); } /** @@ -51,7 +51,7 @@ orientation(Particle const *p_real, * @return Relative distance. */ Utils::Vector3d connection_vector( - Particle const *p_real, + Particle const *p_ref, const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { // Calculate the quaternion defining the orientation of the vector connecting // the virtual site and the real particle @@ -60,7 +60,7 @@ Utils::Vector3d connection_vector( // specifies the relative orientation. auto const director = Utils::convert_quaternion_to_director( - Utils::multiply_quaternions(p_real->r.quat, vs_rel.rel_orientation)) + Utils::multiply_quaternions(p_ref->r.quat, vs_rel.rel_orientation)) .normalize(); return vs_rel.distance * director; @@ -68,81 +68,99 @@ Utils::Vector3d connection_vector( /** * @brief Position of the virtual site - * @param p_real Reference particle for the virtual site. + * @param p_ref Reference particle for the virtual site. * @param vs_rel Parameters for the virtual site. * @return Position of the virtual site. */ Utils::Vector3d -position(Particle const *p_real, +position(Particle const *p_ref, const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { - return p_real->r.p + connection_vector(p_real, vs_rel); + return p_ref->r.p + connection_vector(p_ref, vs_rel); } /** * @brief Velocity of the virtual site - * @param p_real Reference particle for the virtual site. + * @param p_ref Reference particle for the virtual site. * @param vs_rel Parameters for the virtual site. * @return Velocity of the virtual site. */ Utils::Vector3d -velocity(const Particle *p_real, +velocity(const Particle *p_ref, const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { - auto const d = connection_vector(p_real, vs_rel); + auto const d = connection_vector(p_ref, vs_rel); // Get omega of real particle in space-fixed frame auto const omega_space_frame = - convert_vector_body_to_space(*p_real, p_real->m.omega); + convert_vector_body_to_space(*p_ref, p_ref->m.omega); // Obtain velocity from v=v_real particle + omega_real_particle \times // director - return vector_product(omega_space_frame, d) + p_real->m.v; + return vector_product(omega_space_frame, d) + p_ref->m.v; } -Particle *get_real_particle( +/** + * @brief Get reference particle. + * + * @param vs_rel Parameters to get the reference particle for. + * @return Pointer to reference particle. + */ +Particle *get_reference_particle( const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { - auto p_real = get_local_particle_data(vs_rel.to_particle_id); - if (!p_real) { + auto p_ref = get_local_particle_data(vs_rel.to_particle_id); + if (!p_ref) { throw std::runtime_error("No real particle associated with virtual site."); } - return p_real; + return p_ref; +} + +/** + * @brief Contraint force on the real particle. + * + * Calculates the force exerted by the constraint on the + * reference particle. + */ +ParticleForce constraint_force( + const ParticleForce &f, const Particle *p_ref, + const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { + return {f.f, + vector_product(connection_vector(p_ref, vs_rel), f.f) + f.torque}; } /** * @brief Constraint force to hold the particles at its prescribed position. * * @param f Force on the virtual site. - * @param p_real Reference particle. + * @param p_ref Reference particle. * @param vs_real Parameters. * @return Constraint force. */ auto constraint_stress( - const Utils::Vector3d &f, const Particle *p_real, + const Utils::Vector3d &f, const Particle *p_ref, const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) { /* The constraint force is minus the force on the particle, make it force * free. The counter force is translated by the connection vector to the * real particle, hence the virial stress is */ - return tensor_product(connection_vector(p_real, vs_rel), -f); + return tensor_product(connection_vector(p_ref, vs_rel), -f); } } // namespace void VirtualSitesRelative::update(bool recalc_positions) const { // Ghost update logic - if (n_nodes > 1) { - auto const data_parts = - (recalc_positions ? GHOSTTRANS_POSITION : GHOSTTRANS_NONE) | - (get_have_velocity() ? (GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM) - : GHOSTTRANS_NONE); + auto const data_parts = + (recalc_positions ? GHOSTTRANS_POSITION : GHOSTTRANS_NONE) | + (get_have_velocity() ? (GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM) + : GHOSTTRANS_NONE); + + ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); - ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); - } for (auto &p : cell_structure.local_cells().particles()) { if (!p.p.is_virtual) continue; - const Particle *p_real = get_real_particle(p.p.vs_relative); + const Particle *p_ref = get_reference_particle(p.p.vs_relative); if (recalc_positions) { - auto const new_pos = position(p_real, p.p.vs_relative); + auto const new_pos = position(p_ref, p.p.vs_relative); /* The shift has to respect periodic boundaries: if the reference * particles is not in the same image box, we potentially avoid to shift * to the other side of the box. */ @@ -153,10 +171,10 @@ void VirtualSitesRelative::update(bool recalc_positions) const { } if (get_have_velocity()) - p.m.v = velocity(p_real, p.p.vs_relative); + p.m.v = velocity(p_ref, p.p.vs_relative); if (get_have_quaternion()) - p.r.quat = orientation(p_real, p.p.vs_relative); + p.r.quat = orientation(p_ref, p.p.vs_relative); } // namespace } @@ -172,18 +190,10 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const { // We only care about virtual particles if (p.p.is_virtual) { // First obtain the real particle responsible for this virtual particle: - Particle *p_real = get_real_particle(p.p.vs_relative); - - // The rules for transferring forces are: - // F_realParticle +=F_virtualParticle - // T_realParticle +=f_realParticle \times - // (r_virtualParticle-r_realParticle) + Particle *p_ref = get_reference_particle(p.p.vs_relative); // Add forces and torques - p_real->f.torque += - vector_product(connection_vector(p_real, p.p.vs_relative), p.f.f) + - p.f.torque; - p_real->f.f += p.f.f; + p_ref->f += constraint_force(p.f, p_ref, p.p.vs_relative); } } } @@ -197,9 +207,9 @@ Utils::Matrix VirtualSitesRelative::stress_tensor() const { continue; // First obtain the real particle responsible for this virtual particle: - const Particle *p_real = get_real_particle(p.p.vs_relative); + const Particle *p_ref = get_reference_particle(p.p.vs_relative); - stress_tensor += constraint_stress(p.f.f, p_real, p.p.vs_relative); + stress_tensor += constraint_stress(p.f.f, p_ref, p.p.vs_relative); } return stress_tensor; diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 1f2f5b02cec..9c1a138a2a3 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -23,9 +23,9 @@ #include "config.hpp" #ifdef VIRTUAL_SITES_RELATIVE -#include "Particle.hpp" -#include "communication.hpp" -#include "virtual_sites.hpp" +#include "VirtualSites.hpp" + +#include /** @brief Virtual sites implementation for rigid bodies */ class VirtualSitesRelative : public VirtualSites { From b36714016e9c45c7b1ae38d30ad3bbff64a63c8f Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 4 Mar 2020 16:45:49 +0100 Subject: [PATCH 14/16] core: virtual_sites: Always update velocities --- doc/sphinx/particles.rst | 3 +- .../05-raspberry_electrophoresis.ipynb | 2 +- samples/drude_bmimpf6.py | 2 +- samples/rigid_body.py | 3 +- src/core/event.cpp | 2 +- src/core/integrate.cpp | 7 ++--- src/core/virtual_sites/VirtualSites.hpp | 19 +++++-------- .../virtual_sites/VirtualSitesRelative.cpp | 28 ++++++++----------- .../virtual_sites/VirtualSitesRelative.hpp | 2 +- src/python/espressomd/virtual_sites.py | 9 ------ .../virtual_sites/VirtualSites.hpp | 7 +---- testsuite/python/save_checkpoint.py | 3 +- testsuite/python/virtual_sites_relative.py | 16 ++--------- .../python/virtual_sites_tracers_common.py | 1 - 14 files changed, 32 insertions(+), 72 deletions(-) diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index 7f542dee264..7dbc6cfd7b5 100644 --- a/doc/sphinx/particles.rst +++ b/doc/sphinx/particles.rst @@ -309,12 +309,11 @@ To switch the active scheme, the attribute :attr:`espressomd.system.System.virtu from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative system = espressomd.System() - system.virtual_sites = VirtualSitesRelative(have_velocity=True, have_quaternion=False) + system.virtual_sites = VirtualSitesRelative(have_quaternion=False) # or system.virtual_sites = VirtualSitesOff() By default, :class:`espressomd.virtual_sites.VirtualSitesOff` is selected. This means that virtual particles are not touched during integration. -The ``have_velocity`` parameter determines whether or not the velocity of virtual sites is calculated, which carries a performance cost. The ``have_quaternion`` parameter determines whether the quaternion of the virtual particle is updated (useful in combination with the :attr:`espressomd.particle_data.ParticleHandle.vs_quat` property of the virtual particle which defines the orientation of the virtual particle in the body fixed frame of the related real particle. diff --git a/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb b/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb index d10ff376b34..a4035b00950 100644 --- a/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb +++ b/doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb @@ -294,7 +294,7 @@ "outputs": [], "source": [ "# Select the desired implementation for virtual sites\n", - "system.virtual_sites = VirtualSitesRelative(have_velocity=True)\n", + "system.virtual_sites = VirtualSitesRelative()\n", "# Setting min_global_cut is necessary when there is no interaction defined with a range larger than\n", "# the colloid such that the virtual particles are able to communicate their forces to the real particle\n", "# at the center of the colloid\n", diff --git a/samples/drude_bmimpf6.py b/samples/drude_bmimpf6.py index b0714d7c752..86cc66ddc9a 100644 --- a/samples/drude_bmimpf6.py +++ b/samples/drude_bmimpf6.py @@ -75,7 +75,7 @@ print("\n-->Ion pairs:", n_ionpairs, "Box size:", box_l) system = espressomd.System(box_l=[box_l, box_l, box_l]) -system.virtual_sites = VirtualSitesRelative(have_velocity=True) +system.virtual_sites = VirtualSitesRelative() if args.visu: d_scale = 0.988 * 0.5 diff --git a/samples/rigid_body.py b/samples/rigid_body.py index 841e0f771db..b3b70cd9ac2 100644 --- a/samples/rigid_body.py +++ b/samples/rigid_body.py @@ -32,8 +32,7 @@ system = espressomd.System(box_l=[10.0] * 3) -system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_velocity=True) +system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative() system.time_step = 0.01 system.thermostat.set_langevin(kT=1.0, gamma=20.0, seed=42) diff --git a/src/core/event.cpp b/src/core/event.cpp index 1c202e3374a..5df85f482e2 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -423,7 +423,7 @@ unsigned global_ghost_flags() { void update_dependent_particles() { #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); cells_update_ghosts(global_ghost_flags()); #endif diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 9b39347faed..7989c8312c2 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -182,7 +182,7 @@ int integrate(int n_steps, int reuse_forces) { lb_lbcoupling_deactivate(); #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); #endif // Communication step: distribute ghost positions @@ -239,7 +239,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - virtual_sites()->update(true); + virtual_sites()->update(); #endif // Communication step: distribute ghost positions @@ -292,8 +292,7 @@ int integrate(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - // VIRTUAL_SITES update vel - virtual_sites()->update(false); // Recalc positions = false + virtual_sites()->update(); #endif /* verlet list statistics */ diff --git a/src/core/virtual_sites/VirtualSites.hpp b/src/core/virtual_sites/VirtualSites.hpp index 96eef3fd50b..9e63d826668 100644 --- a/src/core/virtual_sites/VirtualSites.hpp +++ b/src/core/virtual_sites/VirtualSites.hpp @@ -38,12 +38,13 @@ /** @brief Base class for virtual sites implementations */ class VirtualSites { public: - VirtualSites() : m_have_velocity(true), m_have_quaternion(false){}; - /** @brief Update positions and/or velocities of virtual sites. - * Velocities are only updated if get_have_velocity() returns true. - * @param recalc_positions Skip the recalculation of positions if false. + VirtualSites() = default; + virtual ~VirtualSites() = default; + + /** + * @brief Update positions and velocities of virtual sites. */ - virtual void update(bool recalc_positions) const {} + virtual void update() const {} /** Back-transfer forces (and torques) to non-virtual particles. */ virtual void back_transfer_forces_and_torques() const {} /** @brief Called after force calculation (and before rattle/shake) */ @@ -51,20 +52,14 @@ class VirtualSites { virtual void after_lb_propagation(){}; /** @brief Pressure contribution. */ virtual Utils::Matrix stress_tensor() const { return {}; }; - /** @brief Enable/disable velocity calculations for vs. */ - void set_have_velocity(const bool &v) { m_have_velocity = v; }; - const bool &get_have_velocity() const { return m_have_velocity; }; /** @brief Enable/disable quaternion calculations for vs.*/ void set_have_quaternion(const bool &have_quaternion) { m_have_quaternion = have_quaternion; }; bool get_have_quaternion() const { return m_have_quaternion; }; - virtual ~VirtualSites() = default; - private: - bool m_have_velocity; - bool m_have_quaternion; + bool m_have_quaternion = false; }; #endif diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 6d666730441..07396b8361e 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -144,12 +144,9 @@ auto constraint_stress( } } // namespace -void VirtualSitesRelative::update(bool recalc_positions) const { +void VirtualSitesRelative::update() const { // Ghost update logic - auto const data_parts = - (recalc_positions ? GHOSTTRANS_POSITION : GHOSTTRANS_NONE) | - (get_have_velocity() ? (GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM) - : GHOSTTRANS_NONE); + auto const data_parts = GHOSTTRANS_POSITION | GHOSTTRANS_MOMENTUM; ghost_communicator(&cell_structure.exchange_ghosts_comm, data_parts); @@ -159,22 +156,19 @@ void VirtualSitesRelative::update(bool recalc_positions) const { const Particle *p_ref = get_reference_particle(p.p.vs_relative); - if (recalc_positions) { - auto const new_pos = position(p_ref, p.p.vs_relative); - /* The shift has to respect periodic boundaries: if the reference - * particles is not in the same image box, we potentially avoid to shift - * to the other side of the box. */ - p.r.p += get_mi_vector(new_pos, p.r.p, box_geo); + auto const new_pos = position(p_ref, p.p.vs_relative); + /* The shift has to respect periodic boundaries: if the reference + * particles is not in the same image box, we potentially avoid to shift + * to the other side of the box. */ + p.r.p += get_mi_vector(new_pos, p.r.p, box_geo); - if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) - set_resort_particles(Cells::RESORT_LOCAL); - } - - if (get_have_velocity()) - p.m.v = velocity(p_ref, p.p.vs_relative); + p.m.v = velocity(p_ref, p.p.vs_relative); if (get_have_quaternion()) p.r.quat = orientation(p_ref, p.p.vs_relative); + + if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) + set_resort_particles(Cells::RESORT_LOCAL); } // namespace } diff --git a/src/core/virtual_sites/VirtualSitesRelative.hpp b/src/core/virtual_sites/VirtualSitesRelative.hpp index 9c1a138a2a3..cfe4e512928 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.hpp +++ b/src/core/virtual_sites/VirtualSitesRelative.hpp @@ -32,7 +32,7 @@ class VirtualSitesRelative : public VirtualSites { public: VirtualSitesRelative() = default; /** @copydoc VirtualSites::update */ - void update(bool recalc_positions) const override; + void update() const override; /** @copydoc VirtualSites::back_transfer_forces_and_torques */ void back_transfer_forces_and_torques() const override; /** @copydoc VirtualSites::stress_tensor */ diff --git a/src/python/espressomd/virtual_sites.py b/src/python/espressomd/virtual_sites.py index ac2e721fe8f..2987a054d8a 100644 --- a/src/python/espressomd/virtual_sites.py +++ b/src/python/espressomd/virtual_sites.py @@ -58,14 +58,5 @@ class VirtualSitesRelative(ScriptInterfaceHelper): """Virtual sites implementation placing virtual sites relative to other particles. See :ref:`Rigid arrangements of particles` for details. - Attributes can be set on the instance or passed to the constructor as - keyword arguments. - - Attributes - ---------- - have_velocity : :obj:`bool` - Determines whether the velocity of the virtual sites is calculated. - This carries a performance cost. - """ _so_name = "VirtualSites::VirtualSitesRelative" diff --git a/src/script_interface/virtual_sites/VirtualSites.hpp b/src/script_interface/virtual_sites/VirtualSites.hpp index 5bccbc3bd80..228b4ffff82 100644 --- a/src/script_interface/virtual_sites/VirtualSites.hpp +++ b/src/script_interface/virtual_sites/VirtualSites.hpp @@ -34,12 +34,7 @@ class VirtualSites : public AutoParameters { public: VirtualSites() { add_parameters( - {{"have_velocity", - [this](const Variant &v) { - virtual_sites()->set_have_velocity(get_value(v)); - }, - [this]() { return virtual_sites()->get_have_velocity(); }}, - {"have_quaternion", + {{"have_quaternion", [this](const Variant &v) { virtual_sites()->set_have_quaternion(get_value(v)); }, diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index c2c9b7253ab..84df33a61df 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -160,8 +160,7 @@ max_displacement=0.01) if espressomd.has_features(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']): - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( - have_velocity=True, have_quaternion=True) + system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(have_quaternion=True) system.part[1].vs_auto_relate_to(0) if espressomd.has_features(['LENNARD_JONES']) and 'LJ' in modes: diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py index 68b4461c6e0..35c72e6f038 100644 --- a/testsuite/python/virtual_sites_relative.py +++ b/testsuite/python/virtual_sites_relative.py @@ -77,9 +77,8 @@ def test_aa_method_switching(self): self.assertIsInstance(self.system.virtual_sites, VirtualSitesOff) # Switch implementation - self.system.virtual_sites = VirtualSitesRelative(have_velocity=False) + self.system.virtual_sites = VirtualSitesRelative() self.assertIsInstance(self.system.virtual_sites, VirtualSitesRelative) - self.assertFalse(self.system.virtual_sites.have_velocity) def test_vs_quat(self): self.system.part.clear() @@ -128,7 +127,7 @@ def test_vs_quat(self): def test_pos_vel_forces(self): system = self.system system.cell_system.skin = 0.3 - system.virtual_sites = VirtualSitesRelative(have_velocity=True) + system.virtual_sites = VirtualSitesRelative() system.box_l = [10, 10, 10] system.part.clear() system.time_step = 0.004 @@ -203,21 +202,12 @@ def test_pos_vel_forces(self): # Check self.assertLessEqual(np.linalg.norm(t_exp - t), 1E-6) - # Check virtual sites without velocity - system.virtual_sites.have_velocity = False - - # Velocity should not change - v2 = np.copy(system.part[2].v) - system.part[1].v = [17, -13.5, 2] - system.integrator.run(0, recalc_forces=True) - np.testing.assert_array_equal(v2, np.copy(system.part[2].v)) - def run_test_lj(self): """This fills the system with vs-based dumbells, adds a lj potential, integrates and verifies forces. This is to make sure that no pairs get lost or are outdated in the short range loop""" system = self.system - system.virtual_sites = VirtualSitesRelative(have_velocity=True) + system.virtual_sites = VirtualSitesRelative() # Parameters n = 90 phi = 0.6 diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index 4ebab12362a..9328e7cf179 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -64,7 +64,6 @@ def test_aa_method_switching(self): self.system.virtual_sites = VirtualSitesInertialessTracers() self.assertIsInstance( self.system.virtual_sites, VirtualSitesInertialessTracers) - self.assertTrue(self.system.virtual_sites.have_velocity) def test_advection(self): self.reset_lb(ext_force_density=[0.1, 0, 0]) From c75698dd389f6ea31ea13a17efd1920886ff1ca3 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 4 Mar 2020 19:58:03 +0100 Subject: [PATCH 15/16] Formatting --- src/core/virtual_sites/VirtualSitesRelative.cpp | 2 +- src/utils/include/utils/Vector.hpp | 2 +- testsuite/python/save_checkpoint.py | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index 07396b8361e..660e404e718 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -131,7 +131,7 @@ ParticleForce constraint_force( * * @param f Force on the virtual site. * @param p_ref Reference particle. - * @param vs_real Parameters. + * @param vs_rel Parameters. * @return Constraint force. */ auto constraint_stress( diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index 91a46baff29..8b4a23aff3c 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -178,7 +178,7 @@ template T trace(Matrix const &m) { * @brief Flatten a matrix to a linear vector. * * @param m Input Matrix - * @return Flat vector with elements of @ref m. + * @return Flat vector with elements of the matrix. */ template Vector flatten(Matrix const &m) { diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 84df33a61df..7d5c524af99 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -160,7 +160,8 @@ max_displacement=0.01) if espressomd.has_features(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']): - system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(have_quaternion=True) + system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative( + have_quaternion=True) system.part[1].vs_auto_relate_to(0) if espressomd.has_features(['LENNARD_JONES']) and 'LJ' in modes: From 3b81f054f9745b1683d7c977a7c8b5f27df795b4 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 4 Mar 2020 20:43:47 +0100 Subject: [PATCH 16/16] utils: Fixed compiler warning --- src/utils/include/utils/Vector.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index 8b4a23aff3c..bdebb05c061 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -167,7 +167,7 @@ template using Matrix = Vector, N>; * @return Trace of matrix. */ template T trace(Matrix const &m) { - T tr = T{}; + auto tr = T{}; for (size_t i = 0; i < N; i++) tr += m[i][i];