Expose local_particles only thru function calls

src/core/PartCfg.hpp

@@ -23,45 +23,15 @@
 #include "ParticleCache.hpp"
 #include "cells.hpp"
 #include "grid.hpp"
-#include "particle_data.hpp"
 
 #include "serialization/Particle.hpp"
-#include <utils/SkipIterator.hpp>
-
-#include <boost/iterator/indirect_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
 
 /**
  * @brief Proxy class that gets a particle range from #local_particles.
  */
 class GetLocalParts {
-  class SkipIfNullOrGhost {
-  public:
-    bool operator()(Particle const *p_ptr) const {
-      return (p_ptr == nullptr) or (p_ptr->l.ghost);
-    }
-  };
-
-  using skip_it = Utils::SkipIterator<Particle **, SkipIfNullOrGhost>;
-  using iterator = boost::indirect_iterator<skip_it>;
-  using Range = boost::iterator_range<iterator>;
-
 public:
-  Range operator()() const {
-    if (local_particles.empty()) {
-      auto begin = skip_it(nullptr, nullptr, SkipIfNullOrGhost());
-      return {make_indirect_iterator(begin), make_indirect_iterator(begin)};
-    }
-
-    auto begin = skip_it(local_particles.data(),
-                         local_particles.data() + max_seen_particle + 1,
-                         SkipIfNullOrGhost());
-    auto end = skip_it(local_particles.data() + max_seen_particle + 1,
-                       local_particles.data() + max_seen_particle + 1,
-                       SkipIfNullOrGhost());
-
-    return {make_indirect_iterator(begin), make_indirect_iterator(end)};
-  }
+  auto operator()() const { return cell_structure.local_cells().particles(); }
 };
 
 /** Unfold coordinates to physical position. */

src/core/cells.cpp

@@ -228,8 +228,8 @@ unsigned topology_check_resort(int cs, unsigned local_resort) {
     local_particles. */
 static void invalidate_ghosts() {
   for (auto const &p : cell_structure.ghost_cells().particles()) {
-    if (local_particles[p.identity()] == &p) {
-      local_particles[p.identity()] = {};
+    if (get_local_particle_data(p.identity()) == &p) {
+      set_local_particle_data(p.identity(), nullptr);
     }
   }
 

src/core/collision.cpp

@@ -66,7 +66,14 @@ Collision_parameters collision_params;
 
 namespace {
 Particle &get_part(int id) {
-  return assert(local_particles[id]), *local_particles[id];
+  auto const p = get_local_particle_data(id);
+
+  if (not p) {
+    throw std::runtime_error("Could not handle collision because particle " +
+                             std::to_string(id) + " was not found.");
+  }
+
+  return *p;
 }
 } // namespace
 
@@ -382,17 +389,17 @@ void bind_at_poc_create_bond_between_vs(const int current_vs_pid,
     // Create bond between the virtual particles
     const int bondG[] = {collision_params.bond_vs, current_vs_pid - 2};
     // Only add bond if vs was created on this node
-    if (local_particles[current_vs_pid - 1])
+    if (get_local_particle_data(current_vs_pid - 1))
       local_add_particle_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 (local_particles[current_vs_pid - 1])
+    if (get_local_particle_data(current_vs_pid - 1))
       local_add_particle_bond(get_part(current_vs_pid - 1), bondG);
-    if (local_particles[current_vs_pid - 2])
+    if (get_local_particle_data(current_vs_pid - 2))
       local_add_particle_bond(get_part(current_vs_pid - 2), bondG);
     break;
   }
@@ -471,9 +478,9 @@ void three_particle_binding_domain_decomposition(
   for (auto &c : gathered_queue) {
     // If we have both particles, at least as ghosts, Get the corresponding cell
     // indices
-    if ((local_particles[c.pp1]) && (local_particles[c.pp2])) {
-      Particle &p1 = *local_particles[c.pp1];
-      Particle &p2 = *local_particles[c.pp2];
+    if (get_local_particle_data(c.pp1) && get_local_particle_data(c.pp2)) {
+      Particle &p1 = *get_local_particle_data(c.pp1);
+      Particle &p2 = *get_local_particle_data(c.pp2);
       auto cell1 = find_current_cell(p1);
       auto cell2 = find_current_cell(p2);
 
@@ -504,7 +511,7 @@ void handle_collisions() {
   if (bind_centers()) {
     for (auto &c : local_collision_queue) {
       // put the bond to the non-ghost particle; at least one partner always is
-      if (local_particles[c.pp1]->l.ghost) {
+      if (get_local_particle_data(c.pp1)->l.ghost) {
         std::swap(c.pp1, c.pp2);
       }
       int bondG[2];
@@ -528,17 +535,14 @@ void handle_collisions() {
     MPI_Allreduce(MPI_IN_PLACE, &max_seen_particle, 1, MPI_INT, MPI_MAX,
                   comm_cart);
 
-    // Make sure, the local_particles array is long enough
-    realloc_local_particles(max_seen_particle);
-
     int current_vs_pid = max_seen_particle + 1;
 
     // Iterate over global collision queue
     for (auto &c : gathered_queue) {
 
       // Get particle pointers
-      Particle *p1 = local_particles[c.pp1];
-      Particle *p2 = local_particles[c.pp2];
+      Particle *p1 = get_local_particle_data(c.pp1);
+      Particle *p2 = get_local_particle_data(c.pp2);
 
       // Only nodes take part in particle creation and binding
       // that see both particles
@@ -587,8 +591,8 @@ void handle_collisions() {
                                               p->identity());
               // Particle storage locations may have changed due to
               // added particle
-              p1 = local_particles[c.pp1];
-              p2 = local_particles[c.pp2];
+              p1 = get_local_particle_data(c.pp1);
+              p2 = get_local_particle_data(c.pp2);
             } else {
               added_particle(current_vs_pid);
             }
@@ -651,8 +655,8 @@ void handle_collisions() {
                                             attach_vs_to.identity());
             // Particle storage locations may have changed due to
             // added particle
-            p1 = local_particles[c.pp1];
-            p2 = local_particles[c.pp2];
+            p1 = get_local_particle_data(c.pp1);
+            p2 = get_local_particle_data(c.pp2);
             current_vs_pid++;
           } else { // Just update the books
             added_particle(current_vs_pid);

src/core/communication.cpp

@@ -54,6 +54,7 @@
 #include "nonbonded_interactions/nonbonded_tab.hpp"
 #include "npt.hpp"
 #include "partCfg_global.hpp"
+#include "particle_data.hpp"
 #include "pressure.hpp"
 #include "rotation.hpp"
 #include "statistics.hpp"

src/core/debug.cpp

@@ -61,12 +61,12 @@ void check_particle_consistency() {
           errexit();
         }
       }
-      if (local_particles[p.p.identity] != &p) {
+      if (get_local_particle_data(p.p.identity) != &p) {
         fprintf(stderr,
                 "%d: check_particle_consistency: ERROR: address "
                 "mismatch for part id %d: local: %p cell: %p in cell %d\n",
                 this_node, p.p.identity,
-                static_cast<void *>(local_particles[p.p.identity]),
+                static_cast<void *>(get_local_particle_data(p.p.identity)),
                 static_cast<void const *>(&p), c);
         errexit();
       }
@@ -86,13 +86,13 @@ void check_particle_consistency() {
 
   /* checks: local particle id */
   for (n = 0; n < max_seen_particle + 1; n++) {
-    if (local_particles[n] != nullptr) {
+    if (get_local_particle_data(n) != nullptr) {
       local_part_cnt++;
-      if (local_particles[n]->p.identity != n) {
+      if (get_local_particle_data(n)->p.identity != n) {
         fprintf(stderr,
                 "%d: check_particle_consistency: ERROR: "
                 "local_particles part %d has corrupted id %d\n",
-                this_node, n, local_particles[n]->p.identity);
+                this_node, n, get_local_particle_data(n)->p.identity);
         errexit();
       }
     }
@@ -119,7 +119,7 @@ void check_particle_consistency() {
     }
 
     for (int p = 0; p < n_part; p++)
-      if (local_particles[p])
+      if (get_local_particle_data(p))
         fprintf(stderr, "%d: got particle %d in local_particles\n", this_node,
                 p);
 

src/core/domain_decomposition.cpp

@@ -651,7 +651,7 @@ void move_if_local(ParticleList &src, ParticleList &rest) {
   for (int i = 0; i < src.n; i++) {
     auto &part = src.part[i];
 
-    assert(local_particles[src.part[i].p.identity] == nullptr);
+    assert(get_local_particle_data(src.part[i].p.identity) == nullptr);
 
     auto target_cell = dd_save_position_to_cell(part.r.p);
 
@@ -684,7 +684,7 @@ void move_left_or_right(ParticleList &src, ParticleList &left,
   for (int i = 0; i < src.n; i++) {
     auto &part = src.part[i];
 
-    assert(local_particles[src.part[i].p.identity] == nullptr);
+    assert(get_local_particle_data(src.part[i].p.identity) == nullptr);
 
     if (get_mi_coord(part.r.p[dir], local_geo.my_left()[dir],
                      box_geo.length()[dir], box_geo.periodic(dir)) < 0.0) {

src/core/energy_inline.hpp

@@ -61,6 +61,7 @@
 #include "bonded_interactions/bonded_coulomb_sr.hpp"
 #include "electrostatics_magnetostatics/coulomb_inline.hpp"
 #endif
+#include "particle_data.hpp"
 #include "statistics.hpp"
 
 #include "energy.hpp"
@@ -215,7 +216,7 @@ inline void add_bonded_energy(Particle const *const p1) {
     int n_partners = iaparams.num;
 
     /* fetch particle 2, which is always needed */
-    Particle const *const p2 = local_particles[p1->bl.e[i++]];
+    Particle const *const p2 = get_local_particle_data(p1->bl.e[i++]);
     if (!p2) {
       runtimeErrorMsg() << "bond broken between particles " << p1->p.identity
                         << " and " << p1->bl.e[i - 1]
@@ -225,7 +226,7 @@ inline void add_bonded_energy(Particle const *const p1) {
 
     /* fetch particle 3 eventually */
     if (n_partners >= 2) {
-      p3 = local_particles[p1->bl.e[i++]];
+      p3 = get_local_particle_data(p1->bl.e[i++]);
       if (!p3) {
         runtimeErrorMsg() << "bond broken between particles " << p1->p.identity
                           << ", " << p1->bl.e[i - 2] << " and "
@@ -237,7 +238,7 @@ inline void add_bonded_energy(Particle const *const p1) {
 
     /* fetch particle 4 eventually */
     if (n_partners >= 3) {
-      p4 = local_particles[p1->bl.e[i++]];
+      p4 = get_local_particle_data(p1->bl.e[i++]);
       if (!p4) {
         runtimeErrorMsg() << "bond broken between particles " << p1->p.identity
                           << ", " << p1->bl.e[i - 3] << ", " << p1->bl.e[i - 2]

src/core/event.cpp

@@ -46,6 +46,7 @@
 #include "npt.hpp"
 #include "nsquare.hpp"
 #include "partCfg_global.hpp"
+#include "particle_data.hpp"
 #include "pressure.hpp"
 #include "random.hpp"
 #include "rattle.hpp"

src/core/forces_inline.hpp

@@ -412,15 +412,15 @@ inline void add_bonded_force(Particle *const p1) {
     bool bond_broken = true;
 
     if (n_partners) {
-      p2 = local_particles[p1->bl.e[i++]];
+      p2 = get_local_particle_data(p1->bl.e[i++]);
       if (!p2) {
         runtimeErrorMsg() << "bond broken between particles " << p1->p.identity;
         return;
       }
 
       /* fetch particle 3 eventually */
       if (n_partners >= 2) {
-        p3 = local_particles[p1->bl.e[i++]];
+        p3 = get_local_particle_data(p1->bl.e[i++]);
         if (!p3) {
           runtimeErrorMsg()
               << "bond broken between particles " << p1->p.identity << ", "
@@ -432,7 +432,7 @@ inline void add_bonded_force(Particle *const p1) {
 
       /* fetch particle 4 eventually */
       if (n_partners >= 3) {
-        p4 = local_particles[p1->bl.e[i++]];
+        p4 = get_local_particle_data(p1->bl.e[i++]);
         if (!p4) {
           runtimeErrorMsg()
               << "bond broken between particles " << p1->p.identity << ", "

src/core/ghosts.cpp

@@ -212,8 +212,8 @@ static void put_recv_buffer(CommBuf &recv_buffer,
       for (Particle &part : part_list->particles()) {
         if (data_parts & GHOSTTRANS_PROPRTS) {
           archiver >> part.p;
-          if (local_particles[part.p.identity] == nullptr) {
-            local_particles[part.p.identity] = ∂
+          if (get_local_particle_data(part.p.identity) == nullptr) {
+            set_local_particle_data(part.p.identity, &part);
           }
         }
         if (data_parts & GHOSTTRANS_POSITION) {

src/core/grid_based_algorithms/lbgpu.cpp

@@ -36,6 +36,7 @@
 #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

@@ -168,15 +168,15 @@ void ImmersedBoundaries::calc_volumes() {
           if (type == BONDED_IA_IBM_TRIEL) {
             // Our particle is the leading particle of a triel
             // Get second and third particle of the triangle
-            Particle const *const p2 = local_particles[p1.bl.e[j + 1]];
+            Particle const *const p2 = get_local_particle_data(p1.bl.e[j + 1]);
             if (!p2) {
               runtimeErrorMsg()
                   << "{IBM_calc_volumes: 078 bond broken between particles "
                   << p1.p.identity << " and " << p1.bl.e[j + 1]
                   << " (particles not stored on the same node)} ";
               return;
             }
-            Particle const *const p3 = local_particles[p1.bl.e[j + 2]];
+            Particle const *const p3 = get_local_particle_data(p1.bl.e[j + 2]);
             if (!p3) {
               runtimeErrorMsg()
                   << "{IBM_calc_volumes: 078 bond broken between particles "
@@ -277,8 +277,8 @@ void ImmersedBoundaries::calc_volume_force() {
           if (type == BONDED_IA_IBM_TRIEL) {
             // Our particle is the leading particle of a triel
             // Get second and third particle of the triangle
-            Particle &p2 = *local_particles[p1.bl.e[j + 1]];
-            Particle &p3 = *local_particles[p1.bl.e[j + 2]];
+            Particle &p2 = *get_local_particle_data(p1.bl.e[j + 1]);
+            Particle &p3 = *get_local_particle_data(p1.bl.e[j + 2]);
 
             // Unfold position of first node.
             // This is to get a continuous trajectory with no jumps when box

src/core/io/mpiio/mpiio.cpp

@@ -377,9 +377,6 @@ void mpi_mpiio_common_read(const char *filename, unsigned fields) {
   // Determine nlocalpart (prefix of rank+1 - own prefix) on every node.
   read_prefs(fnam + ".pref", rank, size, nglobalpart, &pref, &nlocalpart);
 
-  // Prepare ESPResSo data structures
-  local_particles.resize(nglobalpart);
-  std::fill(local_particles.begin(), local_particles.end(), nullptr);
   n_part = nglobalpart;
   max_seen_particle = nglobalpart;
 
@@ -410,7 +407,7 @@ void mpi_mpiio_common_read(const char *filename, unsigned fields) {
                           MPI_INT);
 
     for (int i = 0; i < nlocalpart; ++i)
-      local_particles[id[i]]->p.type = type[i];
+      get_local_particle_data(id[i])->p.type = type[i];
   }
 
   if (fields & MPIIO_OUT_VEL) {
@@ -422,7 +419,7 @@ void mpi_mpiio_common_read(const char *filename, unsigned fields) {
 
     for (int i = 0; i < nlocalpart; ++i)
       for (int k = 0; k < 3; ++k)
-        local_particles[id[i]]->m.v[k] = vel[3 * i + k];
+        get_local_particle_data(id[i])->m.v[k] = vel[3 * i + k];
   }
 
   if (fields & MPIIO_OUT_BND) {
@@ -445,7 +442,7 @@ void mpi_mpiio_common_read(const char *filename, unsigned fields) {
 
     for (int i = 0; i < nlocalpart; ++i) {
       int blen = boff[i + 1] - boff[i];
-      auto &bl = local_particles[id[i]]->bl;
+      auto &bl = get_local_particle_data(id[i])->bl;
       bl.resize(blen);
       std::copy_n(bond.begin() + boff[i], blen, bl.begin());
     }

src/core/io/writer/h5md_core.cpp

@@ -24,6 +24,7 @@
 #include "communication.hpp"
 #include "grid.hpp"
 #include "integrate.hpp"
+#include "particle_data.hpp"
 #include "version.hpp"
 
 #include <fstream>

src/core/layered.cpp

@@ -29,6 +29,7 @@
 #include "errorhandling.hpp"
 #include "ghosts.hpp"
 #include "global.hpp"
+#include "particle_data.hpp"
 
 #include <mpi.h>