Virtual sites refactor

based on review and offline discusison
* pressure calculation moved into the VirtualSites implementation classes.
* virtual sites can decide how many pressure/stress contribs they want
* renamed stress/pressure contribution to "virtual_sites"  in core and py interface
* moved VirtualSites and derived classes into src/core/virtual_sites
* functions are in src/core/virtual_stes.?pp
* unneeded lb checking removed
* unneeded type checking in script interface removed

src/core/CMakeLists.txt

@@ -1,6 +1,6 @@
 file(GLOB EspressoCore_SRC
-          "*.cpp"
-          )
+          "*.cpp" "virtual_sites/VirtualSitesRelative.cpp")
+
 
 if( WITH_COVERAGE )
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -Og --coverage -fprofile-arcs -ftest-coverage")

src/core/collision.cpp

@@ -26,7 +26,7 @@
 #include "grid.hpp"
 #include "domain_decomposition.hpp"
 #include "rotation.hpp"
-#include "virtual_sites_relative.hpp"
+#include "virtual_sites/VirtualSitesRelative.hpp"
 
 using namespace std;
 

src/core/collision.hpp

@@ -21,7 +21,6 @@
 
 #include "particle_data.hpp"
 #include "interaction_data.hpp"
-#include "virtual_sites_relative.hpp"
 #include "virtual_sites.hpp"
 #include "integrate.hpp"
 

src/core/forces.cpp

@@ -103,17 +103,7 @@ void force_calc() {
   }
 #endif
 
-#if defined(VIRTUAL_SITES_RELATIVE) && defined(LB)
-  // This is on a workaround stage:
-  // When using virtual sites relative and LB at the same time, it is necessary
-  // to reassemble the cell lists after all position updates, also of virtual
-  // particles.
-  if ((lattice_switch & LATTICE_LB) &&
-      cell_structure.type == CELL_STRUCTURE_DOMDEC && (!cell_structure.use_verlet_list) &&
-      (std::dynamic_pointer_cast<VirtualSitesRelative>(virtual_sites()) != nullptr))
-        cells_update_ghosts();
-#endif
-
+
   espressoSystemInterface.update();
 
 #ifdef COLLISION_DETECTION

src/core/initialize.cpp

@@ -741,6 +741,4 @@ void on_ghost_flags_change() {
   };
 #endif
 
-//  if (old_have_v != ghosts_have_v)
-//    cells_re_init(CELL_STRUCTURE_CURRENT);
 }

src/core/pressure.cpp

@@ -26,7 +26,7 @@
 #include "cells.hpp"
 #include "integrate.hpp"
 #include "initialize.hpp"
-#include "virtual_sites_relative.hpp" 
+#include "virtual_sites.hpp" 
 #include "npt.hpp"
 #include "p3m.hpp"
 #include "p3m-dipolar.hpp"
@@ -122,11 +122,8 @@ void pressure_calc(double *result, double *result_t, double *result_nb, double *
 
   calc_long_range_virials();
 
-#ifdef VIRTUAL_SITES_RELATIVE  
-  std::shared_ptr<VirtualSitesRelative> vsr =std::dynamic_pointer_cast<VirtualSitesRelative>(virtual_sites());
-  if (vsr) {
-    vsr->pressure_and_stress_tensor_contribution(virials.vs_relative,p_tensor.vs_relative);
-  }
+#ifdef VIRTUAL_SITES
+  virtual_sites()->pressure_and_stress_tensor_contribution(virials.virtual_sites,p_tensor.virtual_sites);
 #endif
 
 
@@ -235,14 +232,14 @@ void init_virials(Observable_stat *stat)
 {
     // Determine number of contribution for different interaction types
     // bonded, nonbonded, coulomb, dipolar, rigid bodies
-    int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vsr;
+    int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vs;
 
   n_pre        = 1;
   n_non_bonded = (n_particle_types*(n_particle_types+1))/2;
 
   n_coulomb    = 0;
   n_dipolar    = 0;
-  n_vsr=0;
+  n_vs=0;
 
 #ifdef ELECTROSTATICS
   switch (coulomb.method) {
@@ -263,14 +260,13 @@ void init_virials(Observable_stat *stat)
       break;
   }
 #endif
-#ifdef VIRTUAL_SITES_RELATIVE
-  // rigid bodies 
-  n_vsr=1;
+#ifdef VIRTUAL_SITES
+  n_vs=virtual_sites()->n_pressure_contribs();
 #endif
 
 
   // Allocate memory for the data
-  obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vsr, 1);
+  obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vs, 1);
   stat->init_status = 0;
 }
 
@@ -291,15 +287,15 @@ void init_p_tensor(Observable_stat *stat)
 {
     // Determine number of contribution for different interaction types
     // bonded, nonbonded, coulomb, dipolar, rigid bodies
-    int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vsr;
+    int n_pre, n_non_bonded, n_coulomb, n_dipolar,n_vs;
 
 
   n_pre        = 1;
   n_non_bonded = (n_particle_types*(n_particle_types+1))/2;
 
   n_coulomb = 0;
   n_dipolar = 0;
-  n_vsr=0;
+  n_vs=0;
 
 #ifdef ELECTROSTATICS
   switch (coulomb.method) {
@@ -319,12 +315,11 @@ void init_p_tensor(Observable_stat *stat)
   default: n_dipolar = 0;
   }
 #endif
-#ifdef VIRTUAL_SITES_RELATIVE
-  // rigid bodies 
-  n_vsr=1;
+#ifdef VIRTUAL_SITES
+  n_vs=virtual_sites()->n_pressure_contribs();
 #endif
 
-  obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vsr, 9);
+  obsstat_realloc_and_clear(stat, n_pre, n_bonded_ia, n_non_bonded, n_coulomb, n_dipolar, n_vs, 9);
   stat->init_status = 0;
 }
 

src/core/statistics.cpp

@@ -1184,12 +1184,12 @@ void analyze_activate(PartCfg &partCfg, int ind) {
 
 void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded,
                                int n_non_bonded, int n_coulomb, int n_dipolar,
-                               int n_vsr, int c_size) {
+                               int n_vs, int c_size) {
 
   int i;
   // Number of doubles to store pressure in
   int total = c_size * (n_pre + n_bonded_ia + n_non_bonded + n_coulomb +
-                        n_dipolar + n_vsr);
+                        n_dipolar + n_vs);
 
   // Allocate mem for the double list
   stat->data.resize(total);
@@ -1201,13 +1201,13 @@ void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded,
   stat->n_coulomb = n_coulomb;
   stat->n_dipolar = n_dipolar;
   stat->n_non_bonded = n_non_bonded;
-  stat->n_vs_relative = n_vsr; // virtual sites relative (rigid bodies)
+  stat->n_virtual_sites = n_vs;
   // Pointers to the start of different contributions
   stat->bonded = stat->data.e + c_size * n_pre;
   stat->non_bonded = stat->bonded + c_size * n_bonded_ia;
   stat->coulomb = stat->non_bonded + c_size * n_non_bonded;
   stat->dipolar = stat->coulomb + c_size * n_coulomb;
-  stat->vs_relative = stat->dipolar + c_size * n_dipolar;
+  stat->virtual_sites = stat->dipolar + c_size * n_dipolar;
 
   // Set all obseravables to zero
   for (i = 0; i < total; i++)

src/core/statistics.hpp

@@ -57,7 +57,7 @@ typedef struct {
   /** number of non bonded interactions */
   int n_non_bonded;
   /** Number of virtual sites relative (rigid body) conributions */
-  int n_vs_relative;
+  int n_virtual_sites;
 
   /** start of bonded interactions. Right after the special ones */
   double *bonded;
@@ -68,7 +68,7 @@ typedef struct {
   /** start of observables for coulomb interaction. */
   double *dipolar;
   /** Start of observables for virtual sites relative (rigid bodies) */
-  double *vs_relative;
+  double *virtual_sites;
 
   /** number of doubles per data item */
   int chunk_size;
@@ -437,7 +437,7 @@ void invalidate_obs();
 
 void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded,
                                int n_non_bonded, int n_coulomb, int n_dipolar,
-                               int n_vsr, int chunk_size);
+                               int n_vs, int chunk_size);
 
 void obsstat_realloc_and_clear_non_bonded(Observable_stat_non_bonded *stat_nb,
                                           int n_nonbonded, int chunk_size_nb);

src/core/virtual_sites.cpp

@@ -24,6 +24,7 @@
 #include "initialize.hpp" 
 #include "statistics.hpp" 
 #include "integrate.hpp" 
+#include "rotation.hpp" 
 #ifdef VIRTUAL_SITES
 
 namespace {
@@ -41,4 +42,111 @@ void set_virtual_sites(std::shared_ptr<VirtualSites> v) {
  on_ghost_flags_change();
 }
 
+#ifdef VIRTUAL_SITES_RELATIVE
+int vs_relate_to(int part_num, int relate_to)
+{
+    // Get the data for the particle we act on and the one we wnat to relate
+    // it to.
+    auto p_current = get_particle_data(part_num);
+    auto p_relate_to = get_particle_data(relate_to);
+    if (!p_current || !p_relate_to) {
+        runtimeErrorMsg()  <<"Could not retrieve particle data for the given id";
+      return ES_ERROR;
+    }
+
+    // get the distance between the particles
+    double d[3];
+    get_mi_vector(d, p_current->r.p,p_relate_to->r.p);
+
+
+
+    // Check, if the distance between virtual and non-virtual particles is larger htan minimum global cutoff
+    // If so, warn user
+    double l=sqrt(sqrlen(d));
+    if (l>min_global_cut) {
+        runtimeErrorMsg() << "Warning: The distance between virtual and non-virtual particle (" << l << ") is\nlarger than the minimum global cutoff (" << min_global_cut << "). This may lead to incorrect simulations\nunder certain conditions. Set the \"System()\" class property \"min_global_cut\" to\nincrease the minimum cutoff.\n";
+      return ES_ERROR;
+    }
+
+    // Now, calculate the quaternions which specify the angle between 
+    // the director of the particel we relate to and the vector
+    // (paritlce_we_relate_to - this_particle)
+    double quat[4];
+    // The vs_relative implemnation later obtains the direcotr by multiplying
+    // the quaternions representing the orientation of the real particle
+    // with those in the virtual particle. The re quulting quaternion is then
+    // converted to a director.
+    // Whe have quat_(real particle) *quat_(virtual particle) 
+    // = quat_(obtained from desired director)
+    // Resolving this for the quat_(virtaul particle)
+
+    //Normalize desired director
+    int i;
+
+    // If the distance between real & virtual particle is 0
+    // we just set the relative orientation to 1 0 0 0, as it is irrelevant but
+    // needs to be a valid quaternion
+    if (l!=0)
+    {
+      for (i=0;i<3;i++)
+        d[i]/=l;
+
+      // Obtain quaternions from desired director
+      double quat_director[4];
+      convert_quatu_to_quat(d, quat_director);
+
+      // Define quat as described above:
+      double x=0;
+      for (i=0;i<4;i++)
+       x+=p_relate_to->r.quat[i]*p_relate_to->r.quat[i];
+
+      quat[0]=0;
+      for (i=0;i<4;i++)
+       quat[0] +=p_relate_to->r.quat[i]*quat_director[i];
+
+      quat[1] =-quat_director[0] *p_relate_to->r.quat[1] 
+         +quat_director[1] *p_relate_to->r.quat[0]
+         +quat_director[2] *p_relate_to->r.quat[3]
+         -quat_director[3] *p_relate_to->r.quat[2];
+      quat[2] =p_relate_to->r.quat[1] *quat_director[3] 
+        + p_relate_to->r.quat[0] *quat_director[2] 
+        - p_relate_to->r.quat[3] *quat_director[1] 
+        - p_relate_to->r.quat[2] * quat_director[0];
+      quat[3] =quat_director[3] *p_relate_to->r.quat[0]
+        - p_relate_to->r.quat[3] *quat_director[0] 
+        + p_relate_to->r.quat[2] * quat_director[1] 
+        - p_relate_to->r.quat[1] *quat_director[2];
+      for (i=0;i<4;i++)
+       quat[i]/=x;
+
+
+     // Verify result
+     double qtemp[4];
+     multiply_quaternions(p_relate_to->r.quat,quat,qtemp);
+     for (i=0;i<4;i++)
+       if (fabs(qtemp[i]-quat_director[i])>1E-9)
+         fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n",
+  	       i, qtemp[i], quat_director[i]);
+    }
+    else
+    { 
+     quat[0]=1;
+     quat[1]=quat[2]=quat[3]=0;
+    }
+
+    // Set the particle id of the particle we want to relate to, the distnace
+    // and the relative orientation
+    if (set_particle_vs_relative(part_num, relate_to, l, quat) == ES_ERROR) {
+      runtimeErrorMsg() << "setting the vs_relative attributes failed";
+      return ES_ERROR;
+    }
+    set_particle_virtual(part_num,1);
+
+   return ES_OK;
+}
+
+
+
+
+#endif
 #endif

src/core/virtual_sites.hpp → src/core/virtual_sites/VirtualSites.hpp

@@ -18,8 +18,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 _VIRTUAL_SITES_HPP
-#define _VIRTUAL_SITES_HPP
+#ifndef VIRTUAL_SITES_VIRTUAL_SITES_HPP
+#define VIRTUAL_SITES_VIRTUAL_SITES_HPP
 
 
 /** \file virtual_sites.hpp
@@ -48,6 +48,10 @@ class VirtualSites {
     virtual void update(bool recalc_positions=true) const =0;
     /** Back-transfer forces (and torques) to non-virtual particles */
     virtual void back_transfer_forces_and_torques() const =0;
+    /** @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& have_velocity() const { return m_have_velocity; };
@@ -62,26 +66,5 @@ class VirtualSites {
    };
 
 
-   /** @brief Do-nothing virtual-sites implementation */
-   class VirtualSitesOff : public VirtualSites {
-    /** @brief Update positions and/or velocities of virtual sites 
-
-    * Velocities are only updated update_velocities() return true 
-    * @param recalc_positions can be used to skip the reculation of positions 
-    */
-    void update(bool recalc_positions=true) const override {};
-    /** Back-transfer forces (and torques) to non-virtual particles */
-    void back_transfer_forces_and_torques() const override {};
-    /** @brief Is a ghost communication needed after position updates */
-    bool need_ghost_comm_after_pos_update() const override { return false;} 
-    /** Is a ghost comm needed before a velocity update */
-    bool need_ghost_comm_before_vel_update() const override {return false;};
-    /** Is a ghost comm needed before back_transfer */
-    bool need_ghost_comm_before_back_transfer() const override {return false;};
-   };
-
-const std::shared_ptr<VirtualSites>& virtual_sites();
-
-void set_virtual_sites(std::shared_ptr<VirtualSites> v); 
 #endif
 #endif