Factor out code duplication

src/core/Particle.hpp

@@ -513,11 +513,15 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   auto const &mu_E() const { return p.mu_E; }
   auto &mu_E() { return p.mu_E; }
 #endif
+#ifdef VIRTUAL_SITES
   auto is_virtual() const {
     return (p.propagation & (PropagationMode::TRANS_VS_RELATIVE |
                              PropagationMode::ROT_VS_RELATIVE |
                              PropagationMode::TRANS_LB_TRACER)) != 0;
   }
+#else
+  constexpr auto is_virtual() const { return false; }
+#endif // VIRTUAL_SITES
 #ifdef VIRTUAL_SITES_RELATIVE
   auto const &vs_relative() const { return p.vs_relative; }
   auto &vs_relative() { return p.vs_relative; }

src/core/PropagationMode.hpp

@@ -23,7 +23,7 @@ namespace PropagationMode {
 /**
  * @brief Flags to create bitmasks for propagation modes.
  */
-enum PropagationMode {
+enum PropagationMode : int {
   NONE = 0,
   SYSTEM_DEFAULT = 1 << 0,
   TRANS_NEWTON = 1 << 1,

src/core/integrate.cpp

@@ -326,22 +326,14 @@ static bool should_propagate_with(Particle const &p, int mode) {
 }
 
 void System::System::thermostats_force_init(double kT) {
-  int bitmask = PropagationMode::TRANS_LANGEVIN | PropagationMode::ROT_LANGEVIN;
-  if (::used_propagations & bitmask) {
-    if (default_propagation & bitmask) {
-      bitmask |= PropagationMode::SYSTEM_DEFAULT;
-    }
-    for (auto &p : cell_structure->local_particles()) {
-      if (p.propagation() & bitmask) {
-        if (should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
-          p.force() += friction_thermo_langevin(langevin, p, time_step, kT);
+  for (auto &p : cell_structure->local_particles()) {
+    if (should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
+      p.force() += friction_thermo_langevin(langevin, p, time_step, kT);
 #ifdef ROTATION
-        if (should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
-          p.torque() += convert_vector_body_to_space(
-              p, friction_thermo_langevin_rotation(langevin, p, time_step, kT));
+    if (should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
+      p.torque() += convert_vector_body_to_space(
+          p, friction_thermo_langevin_rotation(langevin, p, time_step, kT));
 #endif
-      }
-    }
   }
 }
 
@@ -356,8 +348,11 @@ static bool integrator_step_1(ParticleRange const &particles, double kT,
     return steepest_descent_step(particles);
 
   for (auto &p : particles) {
-    if (should_propagate_with(p, TRANS_VS_RELATIVE))
+#ifdef VIRTUAL_SITES
+    // virtual sites are updated later in the integration loop
+    if (p.is_virtual())
       continue;
+#endif
     if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE))
       velocity_verlet_propagator_1(p, time_step);
     if (should_propagate_with(p, TRANS_NEWTON))
@@ -409,8 +404,11 @@ static void integrator_step_2(ParticleRange const &particles, double kT,
   if (integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
     return;
   for (auto &p : particles) {
-    if (should_propagate_with(p, TRANS_VS_RELATIVE))
+#ifdef VIRTUAL_SITES
+    // virtual sites are updated later in the integration loop
+    if (p.is_virtual())
       continue;
+#endif
     if (should_propagate_with(p, TRANS_LB_MOMENTUM_EXCHANGE))
       velocity_verlet_propagator_2(p, time_step);
     if (should_propagate_with(p, TRANS_NEWTON))
@@ -551,8 +549,8 @@ int System::integrate(int n_steps, int reuse_forces) {
 
 #ifdef NPT
     if (integ_switch == INTEG_METHOD_NPT_ISO) {
-      // TODO: workaround to force a ghost update
-      cell_structure->update_ghosts_and_resort_particle(0xffff);
+      cell_structure->update_ghosts_and_resort_particle(
+          Cells::DATA_PART_PROPERTIES);
     }
 #endif
 

src/core/unit_tests/Particle_test.cpp

@@ -24,6 +24,7 @@
 #include <boost/test/unit_test.hpp>
 
 #include "Particle.hpp"
+#include "PropagationMode.hpp"
 #include "config/config.hpp"
 
 #include <utils/Span.hpp>
@@ -36,6 +37,7 @@
 #include <algorithm>
 #include <array>
 #include <sstream>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -290,4 +292,8 @@ BOOST_AUTO_TEST_CASE(particle_bitfields) {
   p.set_cannot_rotate_all_axes();
   BOOST_CHECK(not p.can_rotate());
 #endif
+
+  static_assert(
+      std::is_same_v<std::underlying_type_t<PropagationMode::PropagationMode>,
+                     decltype(ParticleProperties::propagation)>);
 }

src/python/espressomd/particle_data.py

@@ -372,7 +372,6 @@ class ParticleHandle(ScriptInterfaceHelper):
     is_virtual()
         Whether the particle is a virtual site.
 
-
     """
 
     _so_name = "Particles::ParticleHandle"

src/script_interface/particle_data/ParticleHandle.cpp

@@ -593,14 +593,10 @@ Variant ParticleHandle::do_call_method(std::string const &name,
       remove_particle(m_pid);
     });
   } else if (name == "is_virtual") {
-#ifdef VIRTUAL_SITES
     if (not context()->is_head_node()) {
       return {};
     }
     return get_particle_data(m_pid).is_virtual();
-#else
-    return false;
-#endif // VIRTUAL_SITES
 #ifdef VIRTUAL_SITES_RELATIVE
   } else if (name == "vs_relate_to") {
     if (not context()->is_head_node()) {

src/script_interface/system/System.cpp

@@ -59,7 +59,6 @@
 #include <stdexcept>
 #include <string>
 #include <type_traits>
-#include <unordered_map>
 #include <vector>
 
 namespace ScriptInterface {
@@ -218,8 +217,7 @@ static void rotate_system(CellStructure &cell_structure, double phi,
   }
 
   cell_structure.set_resort_particles(Cells::RESORT_GLOBAL);
-  // TODO: workaround to force a ghost update
-  cell_structure.update_ghosts_and_resort_particle(0xffff);
+  cell_structure.update_ghosts_and_resort_particle(Cells::DATA_PART_PROPERTIES);
 }
 
 /** Rescale all particle positions in direction @p dir by a factor @p scale.

testsuite/python/integrator_brownian.py

@@ -181,23 +181,27 @@ def test_propagation(self):
         for various combinations of propagation modes.
         """
         Propagation = espressomd.propagation.Propagation
-        gamma_trans = 1.2
-        gamma_rot = 1.5
-        x0 = [0., 1., 2.]
+        aniso = espressomd.has_features("PARTICLE_ANISOTROPY")
+        gamma_trans = np.array([1.2, 1.2, 1.2]) if aniso else 1.2
+        gamma_rot = np.array([1.5, 1.5, 1.5]) if aniso else 1.5
+        x0 = np.array([0., 1., 2.])
         v0 = np.array([1., 2., 3.])
         o0 = np.array([2., 3., 4.])
         ext_force = np.array([-1., +2., -4.])
         ext_torque = np.array([+1., -3., +5.])
 
+        def has_gamma(gamma):
+            return gamma[0] > 0. if np.shape(gamma) else gamma > 0.
+
         def calc_trajectory(p):
             t = self.system.time
             gamma_t = gamma_trans
             gamma_r = gamma_rot
             if espressomd.has_features("THERMOSTAT_PER_PARTICLE") and \
                     p.propagation & Propagation.SYSTEM_DEFAULT:
-                if p.gamma[0] > 0.:
+                if has_gamma(p.gamma):
                     gamma_t = p.gamma
-                if p.gamma_rot[0] > 0.:
+                if has_gamma(p.gamma_rot):
                     gamma_r = p.gamma_rot
             if (p.propagation & (Propagation.SYSTEM_DEFAULT |
                                  Propagation.TRANS_BROWNIAN)):

testsuite/python/integrator_langevin.py

@@ -174,22 +174,18 @@ def test_03__friction_rot(self):
         """Test the rotational friction-only part of the thermostat."""
 
         system = self.system
+        aniso = espressomd.has_features("PARTICLE_ANISOTROPY")
         o0 = np.array([5., 5., 5.])
-        if espressomd.has_features("PARTICLE_ANISOTROPY"):
-            gamma_t = np.array([0.5, 2., 1.5])
-            gamma_r = np.array([1.5, 0.7, 1.2])
-        else:
-            gamma_t = 2.
-            gamma_r = 3.
+        gamma_t = np.array([0.5, 2., 1.5]) if aniso else 2.
+        gamma_r = np.array([1.5, 0.7, 1.2]) if aniso else 3.
+        rinertia = np.array([1., 1., 1.])
 
         system.time = 0.
         system.time_step = 0.0001
         p = system.part.add(pos=(0, 0, 0), omega_body=o0, rotation=3 * [True])
         if espressomd.has_features("ROTATIONAL_INERTIA"):
-            rinertia = np.array((2, 2, 2))
+            rinertia = np.array([2., 2., 2.])
             p.rinertia = rinertia
-        else:
-            rinertia = np.array((1, 1, 1))
         system.thermostat.set_langevin(
             kT=0, gamma=gamma_t, gamma_rotation=gamma_r, seed=41)
 
@@ -317,23 +313,27 @@ def test_propagation(self):
         for various combinations of propagation modes.
         """
         Propagation = espressomd.propagation.Propagation
-        gamma_trans = 1.2
-        gamma_rot = 0.6
-        x0 = [0., 1., 2.]
+        aniso = espressomd.has_features("PARTICLE_ANISOTROPY")
+        gamma_trans = np.array([1.2, 1.2, 1.2]) if aniso else 1.2
+        gamma_rot = np.array([0.6, 0.6, 0.6]) if aniso else 0.6
+        x0 = np.array([0., 1., 2.])
         v0 = np.array([1., 2., 3.])
         o0 = np.array([1., 2., 1.])
         ext_force = np.array([-2., +2., -3.])
         ext_torque = np.array([+0.25, -0.5, +1.])
 
+        def has_gamma(gamma):
+            return gamma[0] > 0. if np.shape(gamma) else gamma > 0.
+
         def calc_trajectory(p):
             t = self.system.time
             gamma_t = gamma_trans
             gamma_r = gamma_rot
             if espressomd.has_features("THERMOSTAT_PER_PARTICLE") and \
                     p.propagation & Propagation.SYSTEM_DEFAULT:
-                if p.gamma[0] > 0.:
+                if has_gamma(p.gamma):
                     gamma_t = p.gamma
-                if p.gamma_rot[0] > 0.:
+                if has_gamma(p.gamma_rot):
                     gamma_r = p.gamma_rot
             if (p.propagation & (Propagation.SYSTEM_DEFAULT |
                                  Propagation.TRANS_LANGEVIN)):

testsuite/python/integrator_langevin_stats.py

@@ -117,6 +117,10 @@ def verify_diffusion(self, p, corr, kT, gamma):
         acf = c.result()
         tau = c.lag_times()
 
+        gamma = np.copy(gamma)
+        if not gamma.shape:
+            gamma = gamma * np.ones(3)
+
         # Integrate with trapezoidal rule
         for i in range(3):
             I = np.trapz(acf[:, p.id, i], tau)
@@ -130,21 +134,20 @@ def test_06__diffusion(self):
         kT = 1.37
         dt = 0.1
         system.time_step = dt
+        aniso = espressomd.has_features("PARTICLE_ANISOTROPY")
 
         # Translational gamma. We cannot test per-component, if rotation is on,
         # because body and space frames become different.
         gamma = 3.1
 
         # Rotational gamma
-        gamma_rot_i = 4.7
-        gamma_rot_a = [4.2, 1, 1.2]
+        gamma_rot = [4.2, 1, 1.2] if aniso else 4.7
 
         # If we have langevin per particle:
         # Translation
         per_part_gamma = 1.63
         # Rotational
-        per_part_gamma_rot_i = 2.6
-        per_part_gamma_rot_a = [2.4, 3.8, 1.1]
+        per_part_gamma_rot = [2.4, 3.8, 1.1] if aniso else 2.6
 
         # Particle with global thermostat params
         p_global = system.part.add(pos=(0, 0, 0))
@@ -155,25 +158,14 @@ def test_06__diffusion(self):
         if espressomd.has_features("THERMOSTAT_PER_PARTICLE"):
             p_gamma = system.part.add(pos=(0, 0, 0))
             self.setup_diff_mass_rinertia(p_gamma)
-            if espressomd.has_features("PARTICLE_ANISOTROPY"):
-                p_gamma.gamma = 3 * [per_part_gamma]
-                if espressomd.has_features("ROTATION"):
-                    p_gamma.gamma_rot = per_part_gamma_rot_a
-            else:
-                p_gamma.gamma = per_part_gamma
-                if espressomd.has_features("ROTATION"):
-                    p_gamma.gamma_rot = per_part_gamma_rot_i
+            p_gamma.gamma = 3 * [per_part_gamma] if aniso else per_part_gamma
+            if espressomd.has_features("ROTATION"):
+                p_gamma.gamma_rot = per_part_gamma_rot
 
         # Thermostat setup
         if espressomd.has_features("ROTATION"):
-            if espressomd.has_features("PARTICLE_ANISOTROPY"):
-                # particle anisotropy and rotation
-                system.thermostat.set_langevin(
-                    kT=kT, gamma=gamma, gamma_rotation=gamma_rot_a, seed=41)
-            else:
-                # Rotation without particle anisotropy
-                system.thermostat.set_langevin(
-                    kT=kT, gamma=gamma, gamma_rotation=gamma_rot_i, seed=41)
+            system.thermostat.set_langevin(
+                kT=kT, gamma=gamma, gamma_rotation=gamma_rot, seed=41)
         else:
             # No rotation
             system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=41)
@@ -215,31 +207,15 @@ def test_06__diffusion(self):
             corr_omega.finalize()
 
         # Verify diffusion
-        # Translation
-        # Cast gammas to vector, to make checks independent of
-        # PARTICLE_ANISOTROPY
-        gamma = np.ones(3) * gamma
-        per_part_gamma = np.ones(3) * per_part_gamma
         self.verify_diffusion(p_global, corr_vel, kT, gamma)
         if espressomd.has_features("THERMOSTAT_PER_PARTICLE"):
             self.verify_diffusion(p_gamma, corr_vel, kT, per_part_gamma)
 
-        # Rotation
         if espressomd.has_features("ROTATION"):
-            # Decide on effective gamma rotation, since for rotation it is
-            # direction dependent
-            eff_gamma_rot = None
-            if espressomd.has_features("PARTICLE_ANISOTROPY"):
-                eff_gamma_rot = gamma_rot_a
-                eff_per_part_gamma_rot = per_part_gamma_rot_a
-            else:
-                eff_gamma_rot = gamma_rot_i * np.ones(3)
-                eff_per_part_gamma_rot = per_part_gamma_rot_i * np.ones(3)
-
-            self.verify_diffusion(p_global, corr_omega, kT, eff_gamma_rot)
+            self.verify_diffusion(p_global, corr_omega, kT, gamma_rot)
             if espressomd.has_features("THERMOSTAT_PER_PARTICLE"):
                 self.verify_diffusion(
-                    p_gamma, corr_omega, kT, eff_per_part_gamma_rot)
+                    p_gamma, corr_omega, kT, per_part_gamma_rot)
 
     def test_08__noise_correlation(self):
         """Checks that the Langevin noise is uncorrelated"""

testsuite/python/integrator_npt.py

@@ -101,6 +101,7 @@ def calc_trajectory(p, x0):
         modes_trans = [
             Propagation.NONE,
             Propagation.SYSTEM_DEFAULT,
+            Propagation.TRANS_NEWTON,
             Propagation.TRANS_LANGEVIN_NPT,
         ]
         for mode_trans in modes_trans: