Fix virtual sites relative regressions

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
 

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);

src/core/system/System.hpp

@@ -218,17 +218,23 @@ class System : public std::enable_shared_from_this<System> {
   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.).
    */
   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;

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 <cstdio>
 #include <tuple>
 
-#ifdef VIRTUAL_SITES_RELATIVE
-
-/** Calculate the rotation quaternion and distance between two particles */
 std::tuple<Utils::Quaternion<double>, 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

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 <utils/quaternion.hpp>
 
 #include <tuple>
 
+/** Calculate the rotation quaternion and distance between two particles */
 std::tuple<Utils::Quaternion<double>, 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

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.

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(