diff --git a/src/core/collision.cpp b/src/core/collision.cpp index c67aea601c..51e5cc1e83 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -480,7 +480,7 @@ static void three_particle_binding_domain_decomposition( // Handle the collisions stored in the queue void handle_collisions(CellStructure &cell_structure) { - auto const &system = System::get_system(); + auto &system = System::get_system(); auto const &box_geo = *system.box_geo; // Note that the glue to surface mode adds bonds between the centers // but does so later in the process. This is needed to guarantee that @@ -656,6 +656,7 @@ void handle_collisions(CellStructure &cell_structure) { cell_structure.update_ghosts_and_resort_particle( Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS); } + system.update_used_propagations(); } // are we in one of the vs_based methods #endif // defined VIRTUAL_SITES_RELATIVE diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 4f01e75fb8..6ad8653ff1 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -203,9 +203,9 @@ static void update_default_propagation() { } } -static void update_used_propagations(ParticleRange const &particles) { +void System::System::update_used_propagations() { ::used_propagations = PropagationMode::NONE; - for (auto &p : particles) { + for (auto &p : cell_structure->local_particles()) { ::used_propagations |= p.propagation(); } if (::used_propagations & PropagationMode::SYSTEM_DEFAULT) { @@ -444,7 +444,7 @@ int System::integrate(int n_steps, int reuse_forces) { // Prepare particle structure and run sanity checks of all active algorithms update_default_propagation(); - update_used_propagations(cell_structure->local_particles()); + update_used_propagations(); on_integration_start(); // If any method vetoes (e.g. P3M not initialized), immediately bail out @@ -547,16 +547,18 @@ int System::integrate(int n_steps, int reuse_forces) { } #endif +#ifdef VIRTUAL_SITES_RELATIVE + if (::used_propagations & (PropagationMode::TRANS_VS_RELATIVE | + PropagationMode::ROT_VS_RELATIVE)) { #ifdef NPT - if (integ_switch == INTEG_METHOD_NPT_ISO) { - cell_structure->update_ghosts_and_resort_particle( - Cells::DATA_PART_PROPERTIES); + if (integ_switch == INTEG_METHOD_NPT_ISO) { + cell_structure->update_ghosts_and_resort_particle( + Cells::DATA_PART_PROPERTIES); + } +#endif // NPT + vs_relative_update_particles(*cell_structure); } -#endif - -#ifdef VIRTUAL_SITES_RELATIVE - vs_relative_update_particles(*cell_structure); -#endif +#endif // VIRTUAL_SITES_RELATIVE if (cell_structure->get_resort_particles() >= Cells::RESORT_LOCAL) n_verlet_updates++; @@ -568,9 +570,11 @@ int System::integrate(int n_steps, int reuse_forces) { calculate_forces(::temperature); -#ifdef VIRTUAL_SITES - lb_tracers_add_particle_force_to_fluid(*cell_structure, time_step); -#endif +#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS + if (::used_propagations & PropagationMode::TRANS_LB_TRACER) { + lb_tracers_add_particle_force_to_fluid(*cell_structure, time_step); + } +#endif // VIRTUAL_SITES_INERTIALESS_TRACERS integrator_step_2(particles, temperature, time_step); if (integ_switch == INTEG_METHOD_BD) { resort_particles_if_needed(*this); diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index 3c765e155c..32003dbbd7 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -218,10 +218,6 @@ class System : public std::enable_shared_from_this { void on_particle_change(); /** @brief Called every time a particle charge changes. */ void on_particle_charge_change(); - /** @brief Update particles with properties depending on other particles, - * namely virtual sites and ICC charges. - */ - void update_dependent_particles(); /** called before calculating observables, i.e. energy, pressure or * the integrator (forces). Initialize any methods here which are not * initialized immediately (P3M etc.). @@ -229,6 +225,16 @@ class System : public std::enable_shared_from_this { void on_observable_calc(); /**@}*/ + /** + * @brief Update particles with properties depending on other particles, + * namely virtual sites and ICC charges. + */ + void update_dependent_particles(); + /** + * @brief Update the global propagation bitmask. + */ + void update_used_propagations(); + Coulomb::Solver coulomb; Dipoles::Solver dipoles; LB::Solver lb; diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index 7b7a247ddc..bf136b37ce 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -21,7 +21,7 @@ #include "config/config.hpp" -#ifdef VIRTUAL_SITES +#ifdef VIRTUAL_SITES_RELATIVE #include "virtual_sites.hpp" @@ -40,9 +40,6 @@ #include #include -#ifdef VIRTUAL_SITES_RELATIVE - -/** Calculate the rotation quaternion and distance between two particles */ std::tuple, double> calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, BoxGeometry const &box_geo, double min_global_cut, @@ -117,4 +114,3 @@ calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, } #endif // VIRTUAL_SITES_RELATIVE -#endif // VIRTUAL_SITES diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index f8ef445d51..aeafd663fb 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -21,17 +21,17 @@ #include "config/config.hpp" -#ifdef VIRTUAL_SITES -#include "Particle.hpp" - #ifdef VIRTUAL_SITES_RELATIVE #include "BoxGeometry.hpp" +#include "Particle.hpp" +#include "PropagationMode.hpp" #include #include +/** Calculate the rotation quaternion and distance between two particles */ std::tuple, double> calculate_vs_relate_to_params(Particle const &p_current, Particle const &p_relate_to, @@ -50,11 +50,11 @@ inline void vs_relate_to(Particle &p_vs, Particle const &p_relate_to, // Set the particle id of the particle we want to relate to, the distance // and the relative orientation auto &vs_relative = p_vs.vs_relative(); - p_vs.propagation() = PropagationMode::TRANS_VS_RELATIVE; vs_relative.to_particle_id = p_relate_to.id(); std::tie(vs_relative.rel_orientation, vs_relative.distance) = calculate_vs_relate_to_params(p_vs, p_relate_to, box_geo, min_global_cut); + p_vs.propagation() = + PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE; } #endif // VIRTUAL_SITES_RELATIVE -#endif // VIRTUAL_SITES diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index 2a2eaf3e05..16d126febe 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -217,7 +217,8 @@ static void rotate_system(CellStructure &cell_structure, double phi, } cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); - cell_structure.update_ghosts_and_resort_particle(Cells::DATA_PART_PROPERTIES); + cell_structure.update_ghosts_and_resort_particle(Cells::DATA_PART_POSITION | + Cells::DATA_PART_PROPERTIES); } /** Rescale all particle positions in direction @p dir by a factor @p scale. diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index 15ec84be5c..4bdc6d380b 100644 --- a/testsuite/python/integrator_npt.py +++ b/testsuite/python/integrator_npt.py @@ -24,7 +24,7 @@ import espressomd.propagation -@utx.skipIfMissingFeatures(["NPT", "LENNARD_JONES"]) +@utx.skipIfMissingFeatures(["NPT"]) class IntegratorNPT(ut.TestCase): """This tests the NpT integrator interface.""" @@ -39,16 +39,20 @@ def tearDown(self): self.system.part.clear() if espressomd.has_features(["ELECTROSTATICS"]): self.system.electrostatics.clear() - if espressomd.has_features(["DIPOLES"]): - self.system.magnetostatics.clear() + self.reset_rng_counter() + self.system.thermostat.turn_off() + self.system.integrator.set_vv() + if espressomd.has_features(["LENNARD_JONES"]): + self.system.non_bonded_inter[0, 0].lennard_jones.deactivate() + + def reset_rng_counter(self): # reset RNG counter to make tests independent of execution order + self.system.thermostat.set_npt(kT=0., gamma0=0., gammav=1e-6, seed=42) thmst_list = self.system.thermostat.get_state() for thmst in thmst_list: if thmst["type"] == "NPT_ISO": thmst["counter"] = 0 self.system.thermostat.__setstate__(thmst_list) - self.system.thermostat.turn_off() - self.system.integrator.set_vv() def test_integrator_exceptions(self): # invalid parameters should throw exceptions @@ -119,7 +123,29 @@ def calc_trajectory(p, x0): np.testing.assert_allclose(pos, ref_pos, rtol=1e-7) np.testing.assert_allclose(vel, ref_vel, rtol=1e-7) - def test_01_integrator_recovery(self): + @utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE"]) + def test_07__virtual(self): + system = self.system + system.time_step = 0.01 + + virtual = system.part.add(pos=[0, 0, 0], v=[1, 0, 0]) + physical = system.part.add(pos=[0, 0, 0], v=[2, 0, 0]) + virtual.vs_relative = (physical.id, 0.1, (1., 0., 0., 0.)) + + system.thermostat.set_npt(kT=0., gamma0=2., gammav=1e-6, seed=42) + system.integrator.set_isotropic_npt(ext_pressure=0.01, piston=1e6) + + system.integrator.run(1) + + np.testing.assert_almost_equal(np.copy(physical.f), [0., 0., 0.]) + np.testing.assert_almost_equal(np.copy(physical.v), [1.9602, 0., 0.]) + np.testing.assert_almost_equal(np.copy(physical.pos), [0.0198, 0., 0.]) + np.testing.assert_almost_equal(np.copy(virtual.f), [0., 0., 0.]) + np.testing.assert_almost_equal(np.copy(virtual.v), [1.9602, 0., 0.]) + np.testing.assert_almost_equal(np.copy(virtual.pos), [0.0198, 0., 0.1]) + + @utx.skipIfMissingFeatures(["LENNARD_JONES"]) + def test_09_integrator_recovery(self): # the system is still in a valid state after a failure system = self.system np.random.seed(42) @@ -222,7 +248,7 @@ def run_with_p3m(self, container, p3m, method): with self.assertRaisesRegex(Exception, err_msg): system.integrator.run(0, recalc_forces=True) - @utx.skipIfMissingFeatures(["P3M"]) + @utx.skipIfMissingFeatures(["P3M", "LENNARD_JONES"]) def test_npt_p3m_cpu(self): import espressomd.electrostatics p3m = espressomd.electrostatics.P3M( @@ -231,7 +257,7 @@ def test_npt_p3m_cpu(self): self.run_with_p3m(self.system.electrostatics, p3m, "electrostatics") @utx.skipIfMissingGPU() - @utx.skipIfMissingFeatures(["P3M"]) + @utx.skipIfMissingFeatures(["P3M", "LENNARD_JONES"]) def test_npt_p3m_gpu(self): import espressomd.electrostatics p3m = espressomd.electrostatics.P3MGPU(