Cleanup virtual sites

src/core/forces.cpp

@@ -227,7 +227,6 @@ void System::System::calculate_forces(double kT) {
 #endif
 #endif // CUDA
 
-// VIRTUAL_SITES distribute forces
 #ifdef VIRTUAL_SITES_RELATIVE
   vs_relative_back_transfer_forces_and_torques(*cell_structure);
 #endif

src/core/integrate.cpp

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2010-2022 The ESPResSo project
+ * Copyright (C) 2010-2023 The ESPResSo project
  * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
  *   Max-Planck-Institute for Polymer Research, Theory Group
  *
@@ -35,6 +35,7 @@
 
 #include "BoxGeometry.hpp"
 #include "ParticleRange.hpp"
+#include "PropagationMode.hpp"
 #include "accumulators.hpp"
 #include "bond_breakage/bond_breakage.hpp"
 #include "bonded_interactions/rigid_bond.hpp"
@@ -441,6 +442,12 @@ int System::integrate(int n_steps, int reuse_forces) {
 #ifdef CALIPER
   CALI_CXX_MARK_FUNCTION;
 #endif
+#ifdef VIRTUAL_SITES_RELATIVE
+  auto const has_vs_rel = []() {
+    return ::used_propagations & (PropagationMode::ROT_VS_RELATIVE |
+                                  PropagationMode::TRANS_VS_RELATIVE);
+  };
+#endif
 
   // Prepare particle structure and run sanity checks of all active algorithms
   update_default_propagation();
@@ -460,7 +467,9 @@ int System::integrate(int n_steps, int reuse_forces) {
     lb_lbcoupling_deactivate();
 
 #ifdef VIRTUAL_SITES_RELATIVE
-    vs_relative_update_particles(*cell_structure);
+    if (has_vs_rel()) {
+      vs_relative_update_particles(*cell_structure, *box_geo);
+    }
 #endif
 
     // Communication step: distribute ghost positions
@@ -548,15 +557,14 @@ 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)) {
+    if (has_vs_rel()) {
 #ifdef NPT
       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);
+      vs_relative_update_particles(*cell_structure, *box_geo);
     }
 #endif // VIRTUAL_SITES_RELATIVE
 
@@ -574,7 +582,7 @@ int System::integrate(int n_steps, int reuse_forces) {
     if (::used_propagations & PropagationMode::TRANS_LB_TRACER) {
       lb_tracers_add_particle_force_to_fluid(*cell_structure, time_step);
     }
-#endif // VIRTUAL_SITES_INERTIALESS_TRACERS
+#endif
     integrator_step_2(particles, temperature, time_step);
     if (integ_switch == INTEG_METHOD_BD) {
       resort_particles_if_needed(*this);
@@ -627,8 +635,10 @@ int System::integrate(int n_steps, int reuse_forces) {
         }
       }
 
-#ifdef VIRTUAL_SITES
-      lb_tracers_propagate(*cell_structure, time_step);
+#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
+      if (::used_propagations & PropagationMode::TRANS_LB_TRACER) {
+        lb_tracers_propagate(*cell_structure, time_step);
+      }
 #endif
 
 #ifdef COLLISION_DETECTION
@@ -661,7 +671,9 @@ int System::integrate(int n_steps, int reuse_forces) {
 #endif
 
 #ifdef VIRTUAL_SITES_RELATIVE
-  vs_relative_update_particles(*cell_structure);
+  if (has_vs_rel()) {
+    vs_relative_update_particles(*cell_structure, *box_geo);
+  }
 #endif
 
   // Verlet list statistics

src/core/pressure.cpp

@@ -110,8 +110,7 @@ std::shared_ptr<Observable_stat> System::calculate_pressure() {
 
 #ifdef VIRTUAL_SITES_RELATIVE
   if (!obs_pressure.virtual_sites.empty()) {
-    auto const vs_pressure =
-        vs_relative_pressure_tensor(cell_structure->local_particles());
+    auto const vs_pressure = vs_relative_pressure_tensor(*cell_structure);
     boost::copy(Utils::flatten(vs_pressure),
                 obs_pressure.virtual_sites.begin());
   }

src/core/system/System.cpp

@@ -281,7 +281,7 @@ void System::on_particle_charge_change() {
 void System::update_dependent_particles() {
 #ifdef VIRTUAL_SITES
 #ifdef VIRTUAL_SITES_RELATIVE
-  vs_relative_update_particles(*cell_structure);
+  vs_relative_update_particles(*cell_structure, *box_geo);
 #endif
   cell_structure->update_ghosts_and_resort_particle(global_ghost_flags());
 #endif

src/core/virtual_sites/relative.cpp

@@ -28,16 +28,15 @@
 #include "errorhandling.hpp"
 #include "forces.hpp"
 #include "integrate.hpp"
+#include "lees_edwards/lees_edwards.hpp"
 #include "rotation.hpp"
-#include "system/System.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/math/quaternion.hpp>
 #include <utils/math/tensor_product.hpp>
 #include <utils/quaternion.hpp>
 
-#include "integrate.hpp"
-#include "lees_edwards/lees_edwards.hpp"
+#include <cassert>
 
 /**
  * @brief Vector pointing from the real particle to the virtual site.
@@ -79,14 +78,14 @@ static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs) {
  * @param p Virtual site.
  * @return Pointer to real particle.
  */
-static Particle *get_reference_particle(Particle const &p) {
+static Particle *get_reference_particle(CellStructure &cell_structure,
+                                        Particle const &p) {
   auto const &vs_rel = p.vs_relative();
   if (vs_rel.to_particle_id == -1) {
     runtimeErrorMsg() << "Particle with id " << p.id()
                       << " is a dangling virtual site";
     return nullptr;
   }
-  auto &cell_structure = *System::get_system().cell_structure;
   auto p_ref_ptr = cell_structure.get_local_particle(vs_rel.to_particle_id);
   if (!p_ref_ptr) {
     runtimeErrorMsg() << "No real particle with id " << vs_rel.to_particle_id
@@ -109,20 +108,19 @@ static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) {
   return tensor_product(-p_vs.force(), connection_vector(p_ref, p_vs));
 }
 
-static bool is_vs_relative_trans(const Particle &p) {
+static bool is_vs_relative_trans(Particle const &p) {
   return p.propagation() & PropagationMode::TRANS_VS_RELATIVE;
 }
-static bool is_vs_relative_rot(const Particle &p) {
+static bool is_vs_relative_rot(Particle const &p) {
   return p.propagation() & PropagationMode::ROT_VS_RELATIVE;
 }
 
-static bool is_vs_relative(const Particle &p) {
+static bool is_vs_relative(Particle const &p) {
   return (is_vs_relative_trans(p) or is_vs_relative_rot(p));
 }
 
-void vs_relative_update_particles(CellStructure &cell_structure) {
-  auto const &system = System::get_system();
-  auto const &box_geo = *system.box_geo;
+void vs_relative_update_particles(CellStructure &cell_structure,
+                                  BoxGeometry const &box_geo) {
   cell_structure.ghosts_update(Cells::DATA_PART_POSITION |
                                Cells::DATA_PART_MOMENTUM);
 
@@ -132,7 +130,7 @@ void vs_relative_update_particles(CellStructure &cell_structure) {
       continue;
     }
 
-    auto const *p_ref_ptr = get_reference_particle(p);
+    auto const *p_ref_ptr = get_reference_particle(cell_structure, p);
     if (!p_ref_ptr)
       continue;
 
@@ -176,9 +174,8 @@ void vs_relative_back_transfer_forces_and_torques(
     if (!is_vs_relative(p))
       continue;
 
-    auto *p_ref_ptr = get_reference_particle(p);
-    if (!p_ref_ptr)
-      continue;
+    auto *p_ref_ptr = get_reference_particle(cell_structure, p);
+    assert(p_ref_ptr != nullptr);
 
     auto &p_ref = *p_ref_ptr;
     if (is_vs_relative_trans(p)) {
@@ -194,13 +191,13 @@ void vs_relative_back_transfer_forces_and_torques(
 
 // Rigid body contribution to scalar pressure and pressure tensor
 Utils::Matrix<double, 3, 3>
-vs_relative_pressure_tensor(const ParticleRange &particles) {
+vs_relative_pressure_tensor(CellStructure const &cell_structure) {
   Utils::Matrix<double, 3, 3> pressure_tensor = {};
 
-  for (auto const &p : particles) {
+  for (auto const &p : cell_structure.local_particles()) {
     if (is_vs_relative_trans(p)) {
-      auto const *p_ref_ptr = get_reference_particle(p);
-      if (p_ref_ptr) {
+      if (auto const *p_ref_ptr = cell_structure.get_local_particle(
+              p.vs_relative().to_particle_id)) {
         pressure_tensor += constraint_stress(*p_ref_ptr, p);
       }
     }

src/core/virtual_sites/relative.hpp

@@ -24,16 +24,17 @@
 
 #ifdef VIRTUAL_SITES_RELATIVE
 
-#include "ParticleRange.hpp"
+#include "BoxGeometry.hpp"
 #include "cell_system/CellStructure.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/matrix.hpp>
 
-void vs_relative_update_particles(CellStructure &cell_structure);
+void vs_relative_update_particles(CellStructure &cell_structure,
+                                  BoxGeometry const &box_geo);
 void vs_relative_back_transfer_forces_and_torques(
     CellStructure &cell_structure);
 Utils::Matrix<double, 3, 3>
-vs_relative_pressure_tensor(const ParticleRange &particles);
+vs_relative_pressure_tensor(CellStructure const &cell_structure);
 
 #endif // VIRTUAL_SITES_RELATIVE

src/script_interface/system/System.cpp

@@ -111,7 +111,7 @@ System::System() : m_instance{}, m_leaves{std::make_shared<Leaves>()} {
            m_instance->on_boxl_change();
          });
        },
-       []() { return ::System::get_system().box_geo->length(); }},
+       [this]() { return m_instance->box_geo->length(); }},
       {"periodicity",
        [this](Variant const &v) {
          auto const periodicity = get_value<Utils::Vector3b>(v);