Remove virtual flag

samples/immersed_boundary/addSoft.py

@@ -31,7 +31,7 @@ def AddSoft(system, comX, comY, comZ, k1, k2):
             X = float(line[0]) + comX
             Y = float(line[1]) + comY
             Z = float(line[2]) + comZ
-            system.part.add(id=i, pos=[X, Y, Z], virtual=True)
+            system.part.add(id=i, pos=[X, Y, Z])
 
     # triangles
     import espressomd.interactions

samples/immersed_boundary/sampleImmersedBoundary.py

@@ -77,7 +77,7 @@
     agrid=1, density=1, kinematic_viscosity=1, tau=system.time_step,
     ext_force_density=[force, 0, 0])
 system.lb = lbf
-system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0, act_on_virtual=False)
+system.thermostat.set_lb(LB_fluid=lbf, gamma=1.0)
 
 # Setup boundaries
 wall_shapes = [None] * 2

src/core/Particle.hpp

@@ -16,8 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef ESPRESSO_CORE_PARTICLE_HPP
-#define ESPRESSO_CORE_PARTICLE_HPP
+
+#pragma once
 
 #include "config/config.hpp"
 
@@ -79,13 +79,6 @@ struct ParticleProperties {
   /** determines which propagation schemes should be applied to the particle **/
   int propagation = PropagationMode::TRANS_SYSTEM_DEFAULT;
 
-#ifdef VIRTUAL_SITES
-  /** is particle virtual */
-  bool is_virtual = false;
-#else  // VIRTUAL_SITES
-  static constexpr bool is_virtual = false;
-#endif // VIRTUAL_SITES
-
 #ifdef ROTATION
   /** Bitfield for the particle axes of rotation.
    *  Values:
@@ -234,12 +227,9 @@ struct ParticleProperties {
 #ifdef DIPOLE_FIELD_TRACKING
     ar &dip_fld;
 #endif
-#ifdef VIRTUAL_SITES
-    ar &is_virtual;
 #ifdef VIRTUAL_SITES_RELATIVE
     ar &vs_relative;
 #endif
-#endif // VIRTUAL_SITES
 
 #ifdef THERMOSTAT_PER_PARTICLE
     ar &gamma;
@@ -524,18 +514,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 &virtual_flag() { return p.is_virtual; }
-  auto const &virtual_flag() const { return p.is_virtual; }
-  auto is_virtual() const { return p.is_virtual; }
-  void set_virtual(bool const virt_flag) { p.is_virtual = virt_flag; }
+  auto is_virtual() const {
+    return (p.propagation & (PropagationMode::TRANS_VS_RELATIVE |
+                             PropagationMode::ROT_VS_RELATIVE |
+                             PropagationMode::TRANS_LB_TRACER)) != 0;
+  }
 #ifdef VIRTUAL_SITES_RELATIVE
   auto const &vs_relative() const { return p.vs_relative; }
   auto &vs_relative() { return p.vs_relative; }
 #endif // VIRTUAL_SITES_RELATIVE
-#else
-  constexpr auto is_virtual() const { return p.is_virtual; }
-#endif
 #ifdef THERMOSTAT_PER_PARTICLE
   auto const &gamma() const { return p.gamma; }
   auto &gamma() { return p.gamma; }
@@ -633,5 +620,3 @@ BOOST_IS_BITWISE_SERIALIZABLE(ParticleRattle)
 #ifdef VIRTUAL_SITES_RELATIVE
 BOOST_IS_BITWISE_SERIALIZABLE(decltype(ParticleProperties::vs_relative))
 #endif
-
-#endif

src/core/ParticleRange.hpp

@@ -16,8 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef ESPRESSO_SRC_CORE_PARTICLE_RANGE_HPP
-#define ESPRESSO_SRC_CORE_PARTICLE_RANGE_HPP
+
+#pragma once
 
 #include "CellParticleIterator.hpp"
 #include "Particle.hpp"
@@ -51,14 +51,11 @@ class ParticleRange : public boost::iterator_range<CellParticleIterator> {
 
   template <typename Predicate>
   ParticleRangeFiltered<Predicate> filter(Predicate pred) const {
-    return {boost::make_filter_iterator(PropagationPredicate<Predicate>(pred),
-                                        begin(), end()),
-            boost::make_filter_iterator(PropagationPredicate<Predicate>(pred),
-                                        end(), end())};
+    auto const predicate = PropagationPredicate<Predicate>(pred);
+    return {boost::make_filter_iterator(predicate, begin(), end()),
+            boost::make_filter_iterator(predicate, end(), end())};
   }
 
 private:
   base_type::difference_type mutable m_size = -1;
 };
-
-#endif

src/core/PropagationModes.hpp

@@ -19,6 +19,9 @@
 
 #pragma once
 
+/**
+ * @brief Bitfield for propagation modes.
+ */
 enum PropagationMode {
   NONE = 0,
   TRANS_SYSTEM_DEFAULT = 1 << 0,
@@ -34,6 +37,9 @@ enum PropagationMode {
   ROT_BROWNIAN = 1 << 10
 };
 
+/**
+ * @brief Check allowlist of valid propagation modes combinations.
+ */
 inline bool is_valid_propagation_combination(int propagation) {
   if (propagation == NONE)
     return true;

src/core/collision.cpp

@@ -356,7 +356,6 @@ static void place_vs_and_relate_to_particle(CellStructure &cell_structure,
   auto p_vs = cell_structure.add_particle(std::move(new_part));
   vs_relate_to(*p_vs, get_part(cell_structure, relate_to), box_geo,
                min_global_cut);
-  p_vs->set_virtual(true);
   p_vs->type() = collision_params.vs_particle_type;
 }
 

src/core/forces.cpp

@@ -96,14 +96,22 @@ inline ParticleForce thermostat_force(Particle const &p, double time_step,
     return {};
   }
 
+  auto const prop = p.propagation();
+  auto const translation = prop & (PropagationMode::TRANS_LANGEVIN |
+                                   PropagationMode::TRANS_SYSTEM_DEFAULT);
 #ifdef ROTATION
-  return {friction_thermo_langevin(langevin, p, time_step, kT),
-          p.can_rotate() ? convert_vector_body_to_space(
-                               p, friction_thermo_langevin_rotation(
-                                      langevin, p, time_step, kT))
-                         : Utils::Vector3d{}};
+  auto const rotation = prop & (PropagationMode::ROT_LANGEVIN |
+                                PropagationMode::TRANS_SYSTEM_DEFAULT) and
+                        p.can_rotate();
+  return {translation ? friction_thermo_langevin(langevin, p, time_step, kT)
+                      : Utils::Vector3d{},
+          rotation ? convert_vector_body_to_space(
+                         p, friction_thermo_langevin_rotation(langevin, p,
+                                                              time_step, kT))
+                   : Utils::Vector3d{}};
 #else
-  return friction_thermo_langevin(langevin, p, time_step, kT);
+  return translation ? friction_thermo_langevin(langevin, p, time_step, kT)
+                     : Utils::Vector3d{};
 #endif
 }
 
@@ -254,7 +262,7 @@ void System::System::calculate_forces(double kT) {
   immersed_boundaries.volume_conservation(*cell_structure);
 
   if (lb.is_solver_set()) {
-    LB::couple_particles(thermo_virtual, particles, ghost_particles, time_step);
+    LB::couple_particles(particles, ghost_particles, time_step);
   }
 
 #ifdef CUDA

src/core/ghosts.cpp

@@ -126,11 +126,7 @@ serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts,
                      BoxGeometry const &box_geo,
                      Utils::Vector3d const *ghost_shift) {
   if (data_parts & GHOSTTRANS_PROPRTS) {
-    ar &p.id() & p.mol_id() & p.type();
-#ifdef VIRTUAL_SITES
-    ar &p.virtual_flag();
-#endif
-    ar &p.propagation();
+    ar &p.id() & p.mol_id() & p.type() & p.propagation();
 #ifdef ROTATION
     ar &p.rotation();
 #ifdef ROTATIONAL_INERTIA

src/core/integrate.cpp

@@ -332,7 +332,7 @@ static bool integrator_step_1(ParticleRange const &particles, double kT,
   }
 #endif // STOKESIAN_DYNAMICS
 
-  increment_sim_time(time_step);
+  sim_time += time_step;
   return false;
 }
 
@@ -671,8 +671,6 @@ double get_time_step() { return System::get_system().get_time_step(); }
 
 double get_sim_time() { return sim_time; }
 
-void increment_sim_time(double amount) { sim_time += amount; }
-
 void set_time(double value) {
   ::sim_time = value;
   ::recalc_forces = true;

src/core/integrate.hpp

@@ -85,9 +85,6 @@ double get_time_step();
 /** Get simulation time */
 double get_sim_time();
 
-/** Increase simulation time (only on head node) */
-void increment_sim_time(double amount);
-
 /** @brief Set the simulation time. */
 void set_time(double value);
 

src/core/integrators/brownian_inline.hpp

@@ -19,8 +19,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef INTEGRATORS_BROWNIAN_INLINE_HPP
-#define INTEGRATORS_BROWNIAN_INLINE_HPP
+#pragma once
 
 #include "config/config.hpp"
 
@@ -33,13 +32,10 @@
 inline void brownian_dynamics_propagator(BrownianThermostat const &brownian,
                                          Particle &p, double time_step,
                                          double kT) {
-  // Don't propagate translational degrees of freedom of vs
-  if (!p.is_virtual() or thermo_virtual) {
-    p.pos() += bd_drag(brownian.gamma, p, time_step);
-    p.v() = bd_drag_vel(brownian.gamma, p);
-    p.pos() += bd_random_walk(brownian, p, time_step, kT);
-    p.v() += bd_random_walk_vel(brownian, p);
-  }
+  p.pos() += bd_drag(brownian.gamma, p, time_step);
+  p.v() = bd_drag_vel(brownian.gamma, p);
+  p.pos() += bd_random_walk(brownian, p, time_step, kT);
+  p.v() += bd_random_walk_vel(brownian, p);
 }
 
 #ifdef ROTATION
@@ -55,4 +51,3 @@ inline void brownian_dynamics_rotator(BrownianThermostat const &brownian,
   p.omega() += bd_random_walk_vel_rot(brownian, p);
 }
 #endif // ROTATION
-#endif // INTEGRATORS_BROWNIAN_INLINE_HPP

src/core/integrators/steepest_descent.cpp

@@ -55,19 +55,16 @@ bool steepest_descent_step(const ParticleRange &particles) {
     for (unsigned int j = 0; j < 3; j++) {
       // Skip, if coordinate is fixed
       if (!p.is_fixed_along(j)) {
-        // Skip positional increments of virtual particles
-        if (!p.is_virtual()) {
-          // Square of force on particle
-          f += Utils::sqr(p.force()[j]);
-
-          // Positional increment, crop to maximum allowed by user
-          auto const dp =
-              std::clamp(params.gamma * p.force()[j], -params.max_displacement,
-                         params.max_displacement);
-
-          // Move particle
-          p.pos()[j] += dp;
-        }
+        // Square of force on particle
+        f += Utils::sqr(p.force()[j]);
+
+        // Positional increment, crop to maximum allowed by user
+        auto const dp =
+            std::clamp(params.gamma * p.force()[j], -params.max_displacement,
+                       params.max_displacement);
+
+        // Move particle
+        p.pos()[j] += dp;
       }
     }
 #ifdef ROTATION

src/core/integrators/stokesian_dynamics_inline.hpp

@@ -16,12 +16,13 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef STOKESIAN_DYNAMICS_INLINE_HPP
-#define STOKESIAN_DYNAMICS_INLINE_HPP
+
+#pragma once
 
 #include "config/config.hpp"
 
 #ifdef STOKESIAN_DYNAMICS
+
 #include "ParticleRange.hpp"
 #include "communication.hpp"
 #include "integrate.hpp"
@@ -62,4 +63,3 @@ inline void stokesian_dynamics_step_1(const ParticleRange &particles,
 }
 
 #endif // STOKESIAN_DYNAMICS
-#endif

src/core/integrators/velocity_verlet_inline.hpp

@@ -16,8 +16,8 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef INTEGRATORS_VELOCITY_VERLET_HPP
-#define INTEGRATORS_VELOCITY_VERLET_HPP
+
+#pragma once
 
 #include "config/config.hpp"
 
@@ -34,9 +34,6 @@
  */
 inline void velocity_verlet_propagate_vel_pos_par(Particle &p,
                                                   double time_step) {
-  // Don't propagate translational degrees of freedom of vs
-  if (p.is_virtual())
-    return;
   for (unsigned int j = 0; j < 3; j++) {
     if (!p.is_fixed_along(j)) {
       /* Propagate velocities: v(t+0.5*dt) = v(t) + 0.5 * dt * a(t) */
@@ -54,16 +51,10 @@ inline void velocity_verlet_propagate_vel_pos_par(Particle &p,
  */
 inline void velocity_verlet_propagate_vel_final_par(Particle &p,
                                                     double time_step) {
-  // Virtual sites are not propagated during integration
-  if (p.is_virtual())
-    return;
-
   for (unsigned int j = 0; j < 3; j++) {
     if (!p.is_fixed_along(j)) {
       /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */
       p.v()[j] += 0.5 * time_step * p.force()[j] / p.mass();
     }
   }
 }
-
-#endif

src/core/integrators/velocity_verlet_npt.cpp

@@ -51,9 +51,6 @@ velocity_verlet_npt_propagate_vel_final(ParticleRangeNPT const &particles,
   nptiso.p_vel = {};
 
   for (auto &p : particles) {
-    // Virtual sites are not propagated during integration
-    if (p.is_virtual())
-      continue;
     auto const noise = friction_therm0_nptiso<2>(npt_iso, p.v(), p.id());
     for (unsigned int j = 0; j < 3; j++) {
       if (!p.is_fixed_along(j)) {
@@ -127,8 +124,6 @@ static void velocity_verlet_npt_propagate_pos(ParticleRangeNPT const &particles,
 
   /* propagate positions while rescaling positions and velocities */
   for (auto &p : particles) {
-    if (p.is_virtual())
-      continue;
     for (unsigned int j = 0; j < 3; j++) {
       if (!p.is_fixed_along(j)) {
         if (nptiso.geometry & ::nptgeom_dir[j]) {
@@ -170,10 +165,6 @@ static void velocity_verlet_npt_propagate_vel(ParticleRangeNPT const &particles,
   nptiso.p_vel = {};
 
   for (auto &p : particles) {
-
-    // Don't propagate translational degrees of freedom of vs
-    if (p.is_virtual())
-      continue;
     for (unsigned int j = 0; j < 3; j++) {
       if (!p.is_fixed_along(j)) {
         auto const noise = friction_therm0_nptiso<1>(npt_iso, p.v(), p.id());