WIP: re-introduction virtual sites at COM

src/config/features.def

@@ -70,6 +70,7 @@ VIRTUAL_SITES
 VIRTUAL_SITES_RELATIVE          implies VIRTUAL_SITES
 VIRTUAL_SITES_RELATIVE          implies ROTATION
 VIRTUAL_SITES_INERTIALESS_TRACERS implies VIRTUAL_SITES
+VIRTUAL_SITES_CENTER_OF_MASS implies VIRTUAL_SITES
 
 /* DPD features */
 DPD

src/core/virtual_sites/CMakeLists.txt

@@ -22,4 +22,5 @@ target_sources(
   PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lb_inertialess_tracers.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/lb_inertialess_tracers_cuda_interface.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesInertialessTracers.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesRelative.cpp)
+          ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesRelative.cpp
+	  ${CMAKE_CURRENT_SOURCE_DIR}/VirtualSitesCenterOfMass.cpp)

src/core/virtual_sites/VirtualSitesCenterOfMass.cpp

@@ -0,0 +1,115 @@
+/*
+ * Copyright (C) 2010-2022 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "VirtualSitesCenterOfMass.hpp"
+
+// #ifdef VIRTUAL_SITES_CENTER_OF_MASS
+
+#include "Particle.hpp"
+#include "cells.hpp"
+#include "errorhandling.hpp"
+#include "forces.hpp"
+#include "grid.hpp"
+#include "integrate.hpp"
+#include "rotation.hpp"
+#include "communication.hpp"
+
+#include <utils/Vector.hpp>
+#include <utils/math/quaternion.hpp>
+#include <utils/math/tensor_product.hpp>
+#include <utils/quaternion.hpp>
+
+#include <boost/mpi/collectives/all_reduce.hpp>
+
+#include <functional>
+#include <unordered_map>
+
+void VirtualSitesCenterOfMass::update() {
+
+  // com_by_mol_id initialization
+  for (const auto &[mol_id, vs_id] : vitual_site_id_for_mol_id) {
+    com_by_mol_id.emplace(std::make_pair(mol_id, std::make_shared<ComInfo>()));
+  }
+
+  auto const particles = cell_structure.local_particles();
+
+  // Update com_by_mol_id
+  for (const auto &p : particles) {
+    if (com_by_mol_id.find(p.mol_id()) != com_by_mol_id.end()) {
+      com_by_mol_id[p.mol_id()]->total_mass += p.mass();
+      com_by_mol_id[p.mol_id()]->weighted_position_sum +=
+          p.mass() * p.pos(); // Are these the unforlded positions?
+    }
+  }
+
+  // Reduction operation
+  for (auto &kv : com_by_mol_id) {
+    auto const tot_mass =
+        boost::mpi::all_reduce(comm_cart, kv.second->total_mass, std::plus());
+    kv.second->total_mass = tot_mass;
+    auto const weighted_position_sum = boost::mpi::all_reduce(
+        comm_cart, kv.second->weighted_position_sum, std::plus());
+    kv.second->weighted_position_sum = weighted_position_sum;
+  }
+
+  Utils::Vector3d com;
+  int vs_id;
+
+  for (auto &[mol_id, com_info] : com_by_mol_id) {
+    com = com_info->weighted_position_sum /
+          com_info->total_mass; // Is '/' definend in Utils::Vector3d?
+    vs_id = vitual_site_id_for_mol_id[mol_id];
+    auto vs_ptr = cell_structure.get_local_particle(
+        vs_id); // When has the VS been defined?
+    if (vs_ptr == nullptr) {
+      continue;
+    } else {
+      vs_ptr->pos() = com;
+      vs_ptr->mass() = com_info->total_mass;
+
+      auto folded_pos = vs_ptr->pos();
+      auto image_box = Utils::Vector3i{};
+      Utils::Vector3i n_shifts{};
+      fold_position(folded_pos, image_box, box_geo);
+
+      vs_ptr->pos() = folded_pos;
+    }
+  }
+}
+
+// Distribute forces that have accumulated on virtual particles to the
+// associated real particles
+
+void VirtualSitesCenterOfMass::back_transfer_forces() {
+
+  // cell_structure.ghosts_reduce_forces();
+  // init_forces_ghosts(cell_structure.ghost_particles());
+
+  //int p_mol_id;
+  //int vs_id;
+  for (auto &p : cell_structure.local_particles()) {
+    //p_mol_id = p.mol_id();
+    //vs_id = vitual_site_id_for_mol_id[p.mol_id()];
+    auto vs_ptr = cell_structure.get_local_particle(vitual_site_id_for_mol_id[p.mol_id()]);
+    p.force() +=
+        (p.mass() / com_by_mol_id[p.mol_id()]->total_mass) * vs_ptr->force();
+  }
+}
+
+// #endif // VIRTUAL_SITES_CENTER_OF_MASS

src/core/virtual_sites/VirtualSitesCenterOfMass.hpp

@@ -0,0 +1,87 @@
+/*
+ * Copyright (C) 2010-2022 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+
+#ifndef VIRTUAL_SITES_VIRTUAL_SITES_CENTER_OF_MASS_HPP
+#define VIRTUAL_SITES_VIRTUAL_SITES_CENTER_OF_MASS_HPP
+
+/** \file
+ *  This file contains routine to handle virtual sites at the center of mass of
+ * a bunch of other particles (say, a molecule). Forces acting on this center of
+ * mass are distributed back onto the constituents. The position/velocity/mass
+ * of the virtual site at center of mass is calculated from the
+ * positions/velocities/masses of many particles.
+ *
+ *  Virtual sites are like particles, but they will not be integrated.
+ *  Step performed for virtual sites:
+ *  - update virtual sites
+ *  - calculate forces
+ *  - distribute forces
+ *  - move non-virtual particles
+ *  - update virtual sites
+ */
+
+#include <iostream>
+#include "config/config.hpp"
+
+//#ifndef VIRTUAL_SITES_CENTER_OF_MASS
+
+#include "VirtualSites.hpp"
+#include <map>
+#include <memory>
+#include <utils/Vector.hpp>
+#include <utils/matrix.hpp>
+
+/** @brief Base class for virtual sites implementations */
+class VirtualSitesCenterOfMass : public VirtualSites {
+public:
+  VirtualSitesCenterOfMass(const std::unordered_map<int, int> & mid_for_vs) : 
+    vitual_site_id_for_mol_id(mid_for_vs) {}
+  /**
+   * @brief Update positions and velocities of virtual sites.
+   */
+  void update();
+
+  /** @brief Back-transfer forces to non-virtual particles. */
+  void back_transfer_forces();
+
+  /**  @brief Store (mol_id, virtual_site_particle_id) pairs */
+  std::unordered_map<int, int> vitual_site_id_for_mol_id = {};
+
+  auto const &get_mid_for_vs() const { return vitual_site_id_for_mol_id; }
+
+  void set_mid_for_vs( std::unordered_map<int, int> const &vitual_site_id_for_mol_id_ ) {
+    vitual_site_id_for_mol_id = vitual_site_id_for_mol_id_;
+  }
+
+private:
+  struct ComInfo {
+    double total_mass = 0.0;
+    Utils::Vector3d weighted_position_sum = {0., 0., 0.};
+  };
+
+  /**  @brief Store (mol_id, com_info) pairs */
+  std::unordered_map<int, std::shared_ptr<ComInfo>> com_by_mol_id;
+};
+
+//#endif  // VIRTUAL_SITES_CENTER_OF_MASS
+
+#endif

src/python/espressomd/virtual_sites.py

@@ -63,3 +63,18 @@ class VirtualSitesRelative(ScriptInterfaceHelper):
 
         """
         _so_name = "VirtualSites::VirtualSitesRelative"
+
+
+if has_features("VIRTUAL_SITES_CENTER_OF_MASS"):
+    @script_interface_register
+    class VirtualSitesCenterOfMass(ScriptInterfaceHelper):
+
+        """Virtual site implementation for applying forces for center of mass of 
+        group of particles.
+
+        """
+        _so_name = "VirtualSites::VirtualSitesCenterOfMass"
+        _so_creation_policy = "GLOBAL"
+        _so_bind_methods = (
+            "back_transfer_forces",
+        )

src/script_interface/virtual_sites/VirtualSitesCenterOfMass.hpp

@@ -0,0 +1,63 @@
+/*
+ * Copyright (C) 2010-2022 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * 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 SCRIPT_INTERFACE_VIRTUAL_SITES_VIRTUAL_SITES_CENTER_OF_MASS_HPP
+#define SCRIPT_INTERFACE_VIRTUAL_SITES_VIRTUAL_SITES_CENTER_OF_MASS_HPP
+
+#include "config/config.hpp"
+
+#ifdef VIRTUAL_SITES_CENTER_OF_MASS
+
+#include "VirtualSites.hpp"
+
+#include "core/virtual_sites/VirtualSitesCenterOfMass.hpp"
+
+#include <memory>
+#include <map>
+
+namespace ScriptInterface {
+namespace VirtualSites {
+
+class VirtualSitesCenterOfMass : public VirtualSites {
+public:
+
+  void do_construct(VariantMap const &args) override {
+    auto mid_for_vs = get_value<std::unordered_map<int, int>>(args, "mid_for_vs");
+    m_virtual_sites = std::make_shared<::VirtualSitesCenterOfMass>(mid_for_vs);
+
+    add_parameters({
+        {"mid_for_vs", AutoParameter::read_only,
+         [this]() { return make_unordered_map_of_variants(m_virtual_sites->get_mid_for_vs()); }}
+    });
+  }
+
+  std::shared_ptr<::VirtualSites> virtual_sites() override {
+    return m_virtual_sites;
+      } 
+
+private:
+  std::shared_ptr<::VirtualSitesCenterOfMass> m_virtual_sites;
+};
+
+} /* namespace VirtualSites */
+} /* namespace ScriptInterface */
+#endif // VIRTUAL_SITES_CENTER_OF_MASS
+#endif

src/script_interface/virtual_sites/initialize.cpp

@@ -24,6 +24,7 @@
 #include "VirtualSitesInertialessTracers.hpp"
 #include "VirtualSitesOff.hpp"
 #include "VirtualSitesRelative.hpp"
+#include "VirtualSitesCenterOfMass.hpp"
 
 namespace ScriptInterface {
 namespace VirtualSites {
@@ -37,6 +38,9 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
 #endif
 #ifdef VIRTUAL_SITES_RELATIVE
   om->register_new<VirtualSitesRelative>("VirtualSites::VirtualSitesRelative");
+#endif
+#ifdef VIRTUAL_SITES_CENTER_OF_MASS
+  om->register_new<VirtualSitesCenterOfMass>("VirtualSites::VirtualSitesCenterOfMass");
 #endif
   om->register_new<ActiveVirtualSitesHandle>(
       "VirtualSites::ActiveVirtualSitesHandle");