Virtual Sites Cleanup (espressomd#3564)

This is some intermediate cleanup of virtual sites. Biggest change
is that the velocities are now always updated. This has actually very
little performance impact, and is less error prone. 

Description of changes:
 - Simplified interface of `VirtualSites`
 - Some cleanup of the pressure calc
 - Removed `VirtualSites::have_velocity()` mechanism

doc/sphinx/particles.rst

@@ -309,12 +309,11 @@ To switch the active scheme, the attribute :attr:`espressomd.system.System.virtu
     from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative
 
     system = espressomd.System()
-    system.virtual_sites = VirtualSitesRelative(have_velocity=True, have_quaternion=False)
+    system.virtual_sites = VirtualSitesRelative(have_quaternion=False)
     # or
     system.virtual_sites = VirtualSitesOff()
 
 By default, :class:`espressomd.virtual_sites.VirtualSitesOff` is selected. This means that virtual particles are not touched during integration.
-The ``have_velocity`` parameter determines whether or not the velocity of virtual sites is calculated, which carries a performance cost.
 The ``have_quaternion`` parameter determines whether the quaternion of the virtual particle is updated (useful in combination with the
 :attr:`espressomd.particle_data.ParticleHandle.vs_quat` property of the virtual particle which defines the orientation of the virtual particle
 in the body fixed frame of the related real particle.

doc/tutorials/05-raspberry_electrophoresis/05-raspberry_electrophoresis.ipynb

@@ -294,7 +294,7 @@
    "outputs": [],
    "source": [
     "# Select the desired implementation for virtual sites\n",
-    "system.virtual_sites = VirtualSitesRelative(have_velocity=True)\n",
+    "system.virtual_sites = VirtualSitesRelative()\n",
     "# Setting min_global_cut is necessary when there is no interaction defined with a range larger than\n",
     "# the colloid such that the virtual particles are able to communicate their forces to the real particle\n",
     "# at the center of the colloid\n",

samples/drude_bmimpf6.py

@@ -75,7 +75,7 @@
 print("\n-->Ion pairs:", n_ionpairs, "Box size:", box_l)
 
 system = espressomd.System(box_l=[box_l, box_l, box_l])
-system.virtual_sites = VirtualSitesRelative(have_velocity=True)
+system.virtual_sites = VirtualSitesRelative()
 
 if args.visu:
     d_scale = 0.988 * 0.5

samples/rigid_body.py

@@ -32,8 +32,7 @@
 
 
 system = espressomd.System(box_l=[10.0] * 3)
-system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(
-    have_velocity=True)
+system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
 system.time_step = 0.01
 system.thermostat.set_langevin(kT=1.0, gamma=20.0, seed=42)
 

src/core/event.cpp

@@ -423,7 +423,7 @@ unsigned global_ghost_flags() {
 
 void update_dependent_particles() {
 #ifdef VIRTUAL_SITES
-  virtual_sites()->update(true);
+  virtual_sites()->update();
   cells_update_ghosts(global_ghost_flags());
 #endif
 

src/core/ghosts.cpp

@@ -320,6 +320,9 @@ static bool is_poststorable(GhostCommunication const &ghost_comm) {
 }
 
 void ghost_communicator(GhostCommunicator *gcr, unsigned int data_parts) {
+  if (GHOSTTRANS_NONE == data_parts)
+    return;
+
   static CommBuf send_buffer, recv_buffer;
 
   for (auto it = gcr->comm.begin(); it != gcr->comm.end(); ++it) {

src/core/integrate.cpp

@@ -182,7 +182,7 @@ int integrate(int n_steps, int reuse_forces) {
     lb_lbcoupling_deactivate();
 
 #ifdef VIRTUAL_SITES
-    virtual_sites()->update(true);
+    virtual_sites()->update();
 #endif
 
     // Communication step: distribute ghost positions
@@ -239,7 +239,7 @@ int integrate(int n_steps, int reuse_forces) {
 #endif
 
 #ifdef VIRTUAL_SITES
-    virtual_sites()->update(true);
+    virtual_sites()->update();
 #endif
 
     // Communication step: distribute ghost positions
@@ -292,8 +292,7 @@ int integrate(int n_steps, int reuse_forces) {
 #endif
 
 #ifdef VIRTUAL_SITES
-  // VIRTUAL_SITES update vel
-  virtual_sites()->update(false); // Recalc positions = false
+  virtual_sites()->update();
 #endif
 
   /* verlet list statistics */

src/core/integrators/brownian_inline.hpp

@@ -507,11 +507,8 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian,
                                          const ParticleRange &particles) {
   for (auto &p : particles) {
     // Don't propagate translational degrees of freedom of vs
-#ifdef VIRTUAL_SITES
     extern bool thermo_virtual;
-    if (!(p.p.is_virtual) or thermo_virtual)
-#endif
-    {
+    if (!(p.p.is_virtual) or thermo_virtual) {
       p.r.p += bd_drag(brownian.gamma, p, time_step);
       p.m.v = bd_drag_vel(brownian.gamma, p);
       p.r.p += bd_random_walk(brownian, p, time_step);

src/core/integrators/velocity_verlet_inline.hpp

@@ -67,11 +67,10 @@ inline void
 velocity_verlet_propagate_vel_final(const ParticleRange &particles) {
 
   for (auto &p : particles) {
-#ifdef VIRTUAL_SITES
     // Virtual sites are not propagated during integration
     if (p.p.is_virtual)
       continue;
-#endif
+
     for (int j = 0; j < 3; j++) {
       if (!(p.p.ext_flag & COORD_FIXED(j))) {
         /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */

src/core/pressure.cpp

@@ -29,6 +29,7 @@
 #include "npt.hpp"
 #include "pressure_inline.hpp"
 #include "virtual_sites.hpp"
+#include <boost/range/algorithm/copy.hpp>
 
 #include "short_range_loop.hpp"
 
@@ -127,8 +128,12 @@ void pressure_calc(double *result, double *result_t, double *result_nb,
   calc_long_range_virials(cell_structure.local_cells().particles());
 
 #ifdef VIRTUAL_SITES
-  virtual_sites()->pressure_and_stress_tensor_contribution(
-      virials.virtual_sites, p_tensor.virtual_sites);
+  {
+    auto const vs_stress = virtual_sites()->stress_tensor();
+
+    *virials.virtual_sites += trace(vs_stress);
+    boost::copy(flatten(vs_stress), p_tensor.virtual_sites);
+  }
 #endif
 
   for (int n = 1; n < virials.data.n; n++)
@@ -190,7 +195,7 @@ void init_virials(Observable_stat *stat) {
   Dipole::pressure_n();
 #endif
 #ifdef VIRTUAL_SITES
-  n_vs = virtual_sites()->n_pressure_contribs();
+  n_vs = 1;
 #endif
 
   // Allocate memory for the data
@@ -228,7 +233,7 @@ void init_p_tensor(Observable_stat *stat) {
   auto const n_dipolar = 0;
 #endif
 #ifdef VIRTUAL_SITES
-  n_vs = virtual_sites()->n_pressure_contribs();
+  n_vs = 1;
 #endif
 
   obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded,

src/core/pressure_inline.hpp

@@ -73,17 +73,17 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2,
   }
 
 #ifdef ELECTROSTATICS
-  /* real space Coulomb */
-  auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist);
+  {
+    /* real space Coulomb */
+    auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist);
 
-  for (int i = 0; i < 3; i++) {
-    for (int j = 0; j < 3; j++) {
-      p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j];
+    for (int i = 0; i < 3; i++) {
+      for (int j = 0; j < 3; j++) {
+        p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j];
+      }
     }
-  }
 
-  for (int i = 0; i < 3; i++) {
-    virials.coulomb[0] += p_coulomb[i][i];
+    virials.coulomb[0] += trace(p_coulomb);
   }
 #endif /*ifdef ELECTROSTATICS */
 

src/core/virtual_sites.hpp

@@ -26,6 +26,8 @@
 
 #include "virtual_sites/VirtualSites.hpp"
 
+#include <memory>
+
 /** @brief get active virtual sites implementation */
 const std::shared_ptr<VirtualSites> &virtual_sites();
 

src/core/virtual_sites/VirtualSites.hpp

@@ -33,42 +33,33 @@
  */
 
 #ifdef VIRTUAL_SITES
-#include <memory>
+#include <utils/Vector.hpp>
 
 /** @brief Base class for virtual sites implementations */
 class VirtualSites {
 public:
-  VirtualSites() : m_have_velocity(true), m_have_quaternion(false){};
-  /** @brief Update positions and/or velocities of virtual sites.
-   *  Velocities are only updated if get_have_velocity() returns true.
-   *  @param recalc_positions  Skip the recalculation of positions if false.
+  VirtualSites() = default;
+  virtual ~VirtualSites() = default;
+
+  /**
+   * @brief Update positions and velocities of virtual sites.
    */
-  virtual void update(bool recalc_positions) const = 0;
+  virtual void update() const {}
   /** Back-transfer forces (and torques) to non-virtual particles. */
-  virtual void back_transfer_forces_and_torques() const = 0;
+  virtual void back_transfer_forces_and_torques() const {}
   /** @brief Called after force calculation (and before rattle/shake) */
   virtual void after_force_calc(){};
   virtual void after_lb_propagation(){};
-  /** @brief Number of pressure contributions */
-  virtual int n_pressure_contribs() const { return 0; };
   /** @brief Pressure contribution. */
-  virtual void
-  pressure_and_stress_tensor_contribution(double *pressure,
-                                          double *stress_tensor) const {};
-  /** @brief Enable/disable velocity calculations for vs. */
-  void set_have_velocity(const bool &v) { m_have_velocity = v; };
-  const bool &get_have_velocity() const { return m_have_velocity; };
+  virtual Utils::Matrix<double, 3, 3> stress_tensor() const { return {}; };
   /** @brief Enable/disable quaternion calculations for vs.*/
   void set_have_quaternion(const bool &have_quaternion) {
     m_have_quaternion = have_quaternion;
   };
   bool get_have_quaternion() const { return m_have_quaternion; };
 
-  virtual ~VirtualSites() = default;
-
 private:
-  bool m_have_velocity;
-  bool m_have_quaternion;
+  bool m_have_quaternion = false;
 };
 
 #endif

src/core/virtual_sites/VirtualSitesInertialessTracers.hpp

@@ -28,8 +28,6 @@
  * instantaneously transferred to the fluid
  */
 class VirtualSitesInertialessTracers : public VirtualSites {
-  void update(bool recalc_positions) const override{};
-  void back_transfer_forces_and_torques() const override{};
   void after_force_calc() override;
   void after_lb_propagation() override;
 };

src/core/virtual_sites/VirtualSitesOff.hpp

@@ -24,10 +24,7 @@
 #include "VirtualSites.hpp"
 
 /** @brief Do-nothing virtual-sites implementation */
-class VirtualSitesOff : public VirtualSites {
-  void update(bool) const override{};
-  void back_transfer_forces_and_torques() const override{};
-};
+class VirtualSitesOff : public VirtualSites {};
 
 #endif
 #endif