Add option to skip ghost force reduction

src/core/cell_system/CellStructure.cpp

@@ -264,9 +264,9 @@ void CellStructure::set_atom_decomposition(boost::mpi::communicator const &comm,
 
 void CellStructure::set_regular_decomposition(
     boost::mpi::communicator const &comm, double range, BoxGeometry const &box,
-    LocalBox<double> &local_geo) {
-  set_particle_decomposition(
-      std::make_unique<RegularDecomposition>(comm, range, box, local_geo));
+    LocalBox<double> &local_geo, bool without_ghost_force_reduction) {
+  set_particle_decomposition(std::make_unique<RegularDecomposition>(
+      comm, range, box, local_geo, without_ghost_force_reduction));
   m_type = CellStructureType::CELL_STRUCTURE_REGULAR;
   local_geo.set_cell_structure_type(m_type);
 }

src/core/cell_system/CellStructure.hpp

@@ -527,10 +527,12 @@ struct CellStructure {
    * @param range Interaction range.
    * @param box Box Geometry.
    * @param local_geo Geometry of the local box.
+   * @param without_ghost_force_reduction remove the ghost force reduction.
    */
   void set_regular_decomposition(boost::mpi::communicator const &comm,
                                  double range, BoxGeometry const &box,
-                                 LocalBox<double> &local_geo);
+                                 LocalBox<double> &local_geo,
+                                 bool without_ghost_force_reduction);
 
   /**
    * @brief Set the particle decomposition to @ref HybridDecomposition.

src/core/cell_system/HybridDecomposition.cpp

@@ -40,14 +40,17 @@
 #include <set>
 #include <utility>
 
+// TODO: check if we support without_ghost_force_reduction in
+// HybridDecomposition too...
+
 HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
                                          double cutoff_regular,
                                          BoxGeometry const &box_geo,
                                          LocalBox<double> const &local_box,
                                          std::set<int> n_square_types)
     : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular),
       m_regular_decomposition(RegularDecomposition(
-          m_comm, cutoff_regular + skin, m_box, local_box)),
+          m_comm, cutoff_regular + skin, m_box, local_box, false)),
       m_n_square(AtomDecomposition(m_comm, m_box)),
       m_n_square_types(std::move(n_square_types)) {
 

src/core/cell_system/RegularDecomposition.cpp

@@ -483,6 +483,7 @@ void RegularDecomposition::init_cell_interactions() {
 
           auto cell = &cells.at(
               get_linear_index(local_index(neighbor), ghost_cell_grid));
+
           if (ind2 > ind1) {
             red_neighbors.push_back(cell);
           } else {
@@ -633,8 +634,10 @@ GhostCommunicator RegularDecomposition::prepare_comm() {
 RegularDecomposition::RegularDecomposition(boost::mpi::communicator comm,
                                            double range,
                                            BoxGeometry const &box_geo,
-                                           LocalBox<double> const &local_geo)
-    : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo) {
+                                           LocalBox<double> const &local_geo,
+                                           bool without_ghost_force_reduction)
+    : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo),
+      m_without_ghost_force_reduction(without_ghost_force_reduction) {
   /* set up new regular decomposition cell structure */
   create_cell_grid(range);
 

src/core/cell_system/RegularDecomposition.hpp

@@ -85,11 +85,17 @@ struct RegularDecomposition : public ParticleDecomposition {
   std::vector<Cell *> m_ghost_cells;
   GhostCommunicator m_exchange_ghosts_comm;
   GhostCommunicator m_collect_ghost_force_comm;
+  bool m_without_ghost_force_reduction;
 
 public:
   RegularDecomposition(boost::mpi::communicator comm, double range,
                        BoxGeometry const &box_geo,
-                       LocalBox<double> const &local_geo);
+                       LocalBox<double> const &local_geo,
+                       bool without_ghost_force_reduction);
+
+  bool get_without_ghost_force_reduction() const noexcept {
+    return m_without_ghost_force_reduction;
+  }
 
   GhostCommunicator const &exchange_ghosts_comm() const override {
     return m_exchange_ghosts_comm;

src/core/cells.cpp

@@ -203,12 +203,24 @@ void set_hybrid_decomposition(std::set<int> n_square_types,
   on_cell_structure_change();
 }
 
+void set_regular_decomposition(bool without_ghost_force_reduction) {
+  cell_structure.set_regular_decomposition(comm_cart, interaction_range(),
+                                           box_geo, local_geo,
+                                           without_ghost_force_reduction);
+  on_cell_structure_change();
+}
+
 void cells_re_init(CellStructureType new_cs) {
   switch (new_cs) {
-  case CellStructureType::CELL_STRUCTURE_REGULAR:
-    cell_structure.set_regular_decomposition(comm_cart, interaction_range(),
-                                             box_geo, local_geo);
+  case CellStructureType::CELL_STRUCTURE_REGULAR: {
+    auto &current_regular_decomposition =
+        dynamic_cast<RegularDecomposition const &>(
+            std::as_const(cell_structure).decomposition());
+    cell_structure.set_regular_decomposition(
+        comm_cart, interaction_range(), box_geo, local_geo,
+        current_regular_decomposition.get_without_ghost_force_reduction());
     break;
+  }
   case CellStructureType::CELL_STRUCTURE_NSQUARE:
     cell_structure.set_atom_decomposition(comm_cart, box_geo, local_geo);
     break;

src/core/cells.hpp

@@ -84,6 +84,11 @@ extern CellStructure cell_structure;
 void set_hybrid_decomposition(std::set<int> n_square_types,
                               double cutoff_regular);
 
+/** Initialize cell structure @ref RegularDecomposition
+ *  @param without_ghost_force_reduction   Remove the ghost force reduction
+ */
+void set_regular_decomposition(bool without_ghost_force_reduction);
+
 /** Reinitialize the cell structures.
  *  @param new_cs The new topology to use afterwards.
  */

src/core/event.cpp

@@ -71,7 +71,8 @@ void on_program_start() {
   init_node_grid();
 
   /* initially go for regular decomposition */
-  cells_re_init(CellStructureType::CELL_STRUCTURE_REGULAR);
+  set_regular_decomposition(false);
+  // cells_re_init(CellStructureType::CELL_STRUCTURE_REGULAR);
 
   /* make sure interaction 0<->0 always exists */
   make_particle_type_exist(0);

src/script_interface/cell_system/CellSystem.cpp

@@ -125,6 +125,15 @@ CellSystem::CellSystem() {
       {"max_cut_nonbonded", AutoParameter::read_only, maximal_cutoff_nonbonded},
       {"max_cut_bonded", AutoParameter::read_only, maximal_cutoff_bonded},
       {"interaction_range", AutoParameter::read_only, interaction_range},
+      {"without_ghost_force_reduction", AutoParameter::read_only,
+       []() {
+         if (::cell_structure.decomposition_type() !=
+             CellStructureType::CELL_STRUCTURE_REGULAR) {
+           return Variant{none};
+         }
+         auto const rd = get_regular_decomposition();
+         return Variant{rd.get_without_ghost_force_reduction()};
+       }},
   });
 }
 
@@ -255,6 +264,10 @@ void CellSystem::initialize(CellStructureType const &cs_type,
         get_value_or<std::vector<int>>(params, "n_square_types", {});
     auto n_square_types = std::set<int>{ns_types.begin(), ns_types.end()};
     set_hybrid_decomposition(std::move(n_square_types), cutoff_regular);
+  } else if (cs_type == CellStructureType::CELL_STRUCTURE_REGULAR) {
+    auto const without_ghost_force_reduction =
+        get_value_or<bool>(params, "without_ghost_force_reduction", false);
+    set_regular_decomposition(without_ghost_force_reduction);
   } else {
     cells_re_init(cs_type);
   }