Split local particle index from rest of particle data (#3583)

First baby steps to split the front end code from the parallel code.

Changes:
- Move things related to the local particle index into separate files
- Cleaned up/move some of the bond code

src/core/CMakeLists.txt

@@ -24,6 +24,7 @@ set(EspressoCore_SRC
     nsquare.cpp
     partCfg_global.cpp
     particle_data.cpp
+    particle_index.cpp
     polymer.cpp
     polynom.cpp
     pressure.cpp
@@ -41,7 +42,7 @@ set(EspressoCore_SRC
     thermostat.cpp
     tuning.cpp
     virtual_sites.cpp
-)
+    exclusions.cpp)
 
 if(CUDA)
   set(EspressoCuda_SRC

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -19,7 +19,9 @@
 #include "bonded_interaction_data.hpp"
 #include "communication.hpp"
 
+#include <boost/range/algorithm/copy.hpp>
 #include <boost/range/numeric.hpp>
+
 #include <utils/constants.hpp>
 
 std::vector<Bonded_ia_parameters> bonded_ia_params;
@@ -129,3 +131,65 @@ int virtual_set_params(int bond_type) {
 
   return ES_OK;
 }
+
+void remove_all_bonds_to(Particle &p, int id) {
+  IntList *bl = &p.bl;
+  int i, j, partners;
+
+  for (i = 0; i < bl->n;) {
+    partners = bonded_ia_params[bl->e[i]].num;
+    for (j = 1; j <= partners; j++)
+      if (bl->e[i + j] == id)
+        break;
+    if (j <= partners) {
+      bl->erase(bl->begin() + i, bl->begin() + i + 1 + partners);
+    } else
+      i += 1 + partners;
+  }
+  assert(i == bl->n);
+}
+
+void add_bond(Particle &p, Utils::Span<const int> bond) {
+  boost::copy(bond, std::back_inserter(p.bl));
+}
+
+int delete_bond(Particle *part, const int *bond) {
+  IntList *bl = &part->bl;
+  int i, j, type, partners;
+
+  // Empty bond means: delete all bonds
+  if (!bond) {
+    bl->clear();
+
+    return ES_OK;
+  }
+
+  // Go over the bond list to find the bond to delete
+  for (i = 0; i < bl->n;) {
+    type = bl->e[i];
+    partners = bonded_ia_params[type].num;
+
+    // If the bond type does not match the one, we want to delete, skip
+    if (type != bond[0])
+      i += 1 + partners;
+    else {
+      // Go over the bond partners
+      for (j = 1; j <= partners; j++) {
+        // Leave the loop early, if the bond to delete and the bond with in the
+        // particle don't match
+        if (bond[j] != bl->e[i + j])
+          break;
+      }
+      // If we did not exit from the loop early, all parameters matched
+      // and we go on with deleting
+      if (j > partners) {
+        // New length of bond list
+        bl->erase(bl->begin() + i, bl->begin() + i + 1 + partners);
+
+        return ES_OK;
+      }
+      i += 1 + partners;
+    }
+  }
+  return ES_ERROR;
+}

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -19,11 +19,12 @@
 #ifndef _BONDED_INTERACTION_DATA_HPP
 #define _BONDED_INTERACTION_DATA_HPP
 
-#include <boost/optional.hpp>
-
 #include "Particle.hpp"
 #include "TabulatedPotential.hpp"
-#include <utils/Counter.hpp>
+
+#include <utils/Span.hpp>
+
+#include <boost/optional.hpp>
 
 /** @file
  *  Data structures for bonded interactions.
@@ -476,6 +477,24 @@ inline bool pair_bond_enum_exists_between(Particle const &p1,
          (p2.bl.n > 0 && pair_bond_enum_exists_on(p2, p1, bond));
 }
 
+/** @brief Add bond to local particle.
+ *  @param p     identity of principal atom of the bond.
+ *  @param bond  field containing the bond type number and the identity
+ *               of all bond partners (secondary atoms of the bond).
+ */
+void add_bond(Particle &p, Utils::Span<const int> bond);
+
+/** Remove bond from particle. */
+int delete_bond(Particle *part, const int *bond);
+
+/**
+ * @brief Remove all bonds on particle involving another particle.
+ *
+ * @param p Particle whose bond list is modified.
+ * @param id Bonds involving this id are removed.
+ */
+void remove_all_bonds_to(Particle &p, int id);
+
 /** Calculate the maximal cutoff of bonded interactions, required to
  *  determine the cell size for communication.
  *

src/core/cells.cpp

@@ -38,15 +38,14 @@
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "nsquare.hpp"
 #include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <utils/NoOp.hpp>
 #include <utils/mpi/gather_buffer.hpp>
 
 #include <boost/iterator/indirect_iterator.hpp>
 
 #include <cstdio>
-#include <cstdlib>
-#include <cstring>
 
 /** list of all cells. */
 std::vector<Cell> cells;

src/core/collision.cpp

@@ -27,8 +27,7 @@
 #include "event.hpp"
 #include "grid.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "particle_data.hpp"
-#include "rotation.hpp"
+#include "particle_index.hpp"
 #include "virtual_sites/VirtualSitesRelative.hpp"
 
 #include <utils/mpi/all_compare.hpp>
@@ -362,7 +361,7 @@ void coldet_do_three_particle_bond(Particle &p, Particle &p1, Particle &p2) {
   // First, fill bond data structure
   const Utils::Vector3i bondT = {bond_id, p1.p.identity, p2.p.identity};
 
-  local_add_particle_bond(p, bondT);
+  add_bond(p, bondT);
 }
 
 #ifdef VIRTUAL_SITES_RELATIVE
@@ -389,17 +388,17 @@ void bind_at_poc_create_bond_between_vs(const int current_vs_pid,
     const int bondG[] = {collision_params.bond_vs, current_vs_pid - 2};
     // Only add bond if vs was created on this node
     if (get_local_particle_data(current_vs_pid - 1))
-      local_add_particle_bond(get_part(current_vs_pid - 1), bondG);
+      add_bond(get_part(current_vs_pid - 1), bondG);
     break;
   }
   case 2: {
     // Create 1st bond between the virtual particles
     const int bondG[] = {collision_params.bond_vs, c.pp1, c.pp2};
     // Only add bond if vs was created on this node
     if (get_local_particle_data(current_vs_pid - 1))
-      local_add_particle_bond(get_part(current_vs_pid - 1), bondG);
+      add_bond(get_part(current_vs_pid - 1), bondG);
     if (get_local_particle_data(current_vs_pid - 2))
-      local_add_particle_bond(get_part(current_vs_pid - 2), bondG);
+      add_bond(get_part(current_vs_pid - 2), bondG);
     break;
   }
   }
@@ -413,9 +412,9 @@ void glue_to_surface_bind_part_to_vs(const Particle *const p1,
   const int bondG[] = {collision_params.bond_vs, vs_pid_plus_one - 1};
 
   if (p1->p.type == collision_params.part_type_after_glueing) {
-    local_add_particle_bond(get_part(p1->p.identity), bondG);
+    add_bond(get_part(p1->p.identity), bondG);
   } else {
-    local_add_particle_bond(get_part(p2->p.identity), bondG);
+    add_bond(get_part(p2->p.identity), bondG);
   }
 }
 
@@ -514,7 +513,7 @@ void handle_collisions() {
       int bondG[2];
       bondG[0] = collision_params.bond_centers;
       bondG[1] = c.pp2;
-      local_add_particle_bond(get_part(c.pp1), bondG);
+      add_bond(get_part(c.pp1), bondG);
     }
   }
 
@@ -629,7 +628,7 @@ void handle_collisions() {
             int bondG[2];
             bondG[0] = collision_params.bond_centers;
             bondG[1] = c.pp2;
-            local_add_particle_bond(get_part(c.pp1), bondG);
+            add_bond(get_part(c.pp1), bondG);
           }
 
           // Change type of particle being attached, to make it inert

src/core/communication.cpp

@@ -276,13 +276,10 @@ void mpi_remove_particle(int pnode, int part) {
 
 void mpi_remove_particle_slave(int pnode, int part) {
   if (part != -1) {
-    if (pnode == this_node) {
-      local_remove_particle(part);
-    } else {
-      remove_all_bonds_to(part);
-    }
-  } else
+    local_remove_particle(part);
+  } else {
     local_remove_all_particles();
+  }
 
   on_particle_change();
 }

src/core/debug.cpp

@@ -25,7 +25,7 @@
 #include "debug.hpp"
 
 #include "cells.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <stdexcept>
 

src/core/domain_decomposition.cpp

@@ -28,7 +28,7 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 #include "grid.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include "serialization/ParticleList.hpp"
 #include <utils/index.hpp>

src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp

@@ -28,8 +28,6 @@
 #include "electrostatics_magnetostatics/dipole.hpp"
 #include "errorhandling.hpp"
 #include "grid.hpp"
-#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "particle_data.hpp"
 
 #include <utils/constants.hpp>
 

src/core/energy_inline.hpp

@@ -61,7 +61,8 @@
 #include "bonded_interactions/bonded_coulomb_sr.hpp"
 #include "electrostatics_magnetostatics/coulomb_inline.hpp"
 #endif
-#include "particle_data.hpp"
+#include "exclusions.hpp"
+#include "particle_index.hpp"
 #include "statistics.hpp"
 
 #include "energy.hpp"

src/core/exclusions.cpp

@@ -0,0 +1,19 @@
+#include "exclusions.hpp"
+
+#ifdef EXCLUSIONS
+void add_exclusion(Particle *part, int part2) {
+  for (int i = 0; i < part->el.n; i++)
+    if (part->el.e[i] == part2)
+      return;
+
+  part->el.push_back(part2);
+}
+
+void delete_exclusion(Particle *part, int part2) {
+  IntList &el = part->el;
+
+  if (!el.empty()) {
+    el.erase(std::remove(el.begin(), el.end(), part2), el.end());
+  };
+}
+#endif

src/core/exclusions.hpp

@@ -0,0 +1,23 @@
+#ifndef ESPRESSO_EXCLUSIONS_HPP
+#define ESPRESSO_EXCLUSIONS_HPP
+
+#include "Particle.hpp"
+
+#ifdef EXCLUSIONS
+/** Determine if the non-bonded interactions between @p p1 and @p p2 should be
+ *  calculated.
+ */
+inline bool do_nonbonded(Particle const &p1, Particle const &p2) {
+  /* check for particle 2 in particle 1's exclusion list. The exclusion list is
+   * symmetric, so this is sufficient. */
+  return std::none_of(p1.el.begin(), p1.el.end(),
+                      [&p2](int id) { return p2.p.identity == id; });
+}
+
+/** Remove exclusion from particle if possible */
+void delete_exclusion(Particle *part, int part2);
+
+/** Insert an exclusion if not already set */
+void add_exclusion(Particle *part, int part2);
+#endif
+#endif // ESPRESSO_EXCLUSIONS_HPP

src/core/forces_inline.hpp

@@ -36,6 +36,7 @@
 #include "bonded_interactions/thermalized_bond.hpp"
 #include "bonded_interactions/umbrella.hpp"
 #include "errorhandling.hpp"
+#include "exclusions.hpp"
 #include "forces.hpp"
 #include "immersed_boundary/ibm_tribend.hpp"
 #include "immersed_boundary/ibm_triel.hpp"
@@ -59,7 +60,7 @@
 #include "npt.hpp"
 #include "object-in-fluid/oif_global_forces.hpp"
 #include "object-in-fluid/oif_local_forces.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 #include "rotation.hpp"
 #include "thermostat.hpp"
 

src/core/ghosts.cpp

@@ -31,7 +31,7 @@
 #include "ghosts.hpp"
 #include "Particle.hpp"
 #include "communication.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <mpi.h>
 #include <utils/Span.hpp>

src/core/global.cpp

@@ -33,7 +33,6 @@
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "npt.hpp"
 #include "object-in-fluid/oif_global_forces.hpp"
-#include "particle_data.hpp"
 #include "rattle.hpp"
 #include "thermostat.hpp"
 #include "tuning.hpp"

src/core/grid_based_algorithms/lbgpu.cpp

@@ -36,7 +36,6 @@
 #include "grid_based_algorithms/lbgpu.hpp"
 #include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "particle_data.hpp"
 #include "statistics.hpp"
 
 #include <utils/constants.hpp>

src/core/immersed_boundary/ImmersedBoundaries.cpp

@@ -25,7 +25,7 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 #include "grid.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <utils/constants.hpp>
 

src/core/io/mpiio/mpiio.cpp

@@ -55,17 +55,16 @@
 #include "cells.hpp"
 #include "errorhandling.hpp"
 #include "event.hpp"
-#include "integrate.hpp"
 #include "mpiio.hpp"
 #include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <mpi.h>
 
 #include <cerrno>
 #include <cstring>
 #include <string>
 #include <sys/stat.h>
-#include <unistd.h>
 #include <vector>
 
 namespace Mpiio {

src/core/nsquare.cpp

@@ -27,7 +27,7 @@
 #include "communication.hpp"
 #include "constraints.hpp"
 #include "ghosts.hpp"
-#include "particle_data.hpp"
+#include "particle_index.hpp"
 
 #include <mpi.h>