Local particles (#3625)

The goal here is to make it transparent which call sites actually need access to
the cells, and which just need to iterate the particles. There are a quite a lot of
changes her, but most of them are just textual replacements, so should be fast
to check. The final goal here is to be able to hide the cells from client code of
the particle storage. From my side #2899 can be closed with this.

Description of changes:
- Add functions {local_particles, ghost_particles} to CellStructure to directly get the local/ghost particle ranges
- Replaced all occurrences of particles access thru the {local, ghost} cell ranges by
  direct calls to the aforementioned functions
- Replaced CellPList by the now functionally equivalent Utils::Span

src/core/CellStructure.hpp

@@ -53,25 +53,16 @@ enum Resort : unsigned {
 };
 }
 
-/** List of cell pointers. */
-struct CellPList {
-  ParticleRange particles() const {
-    /* Find first non-empty cell */
-    auto first = std::find_if(cell, cell + n,
-                              [](const Cell *c) { return not c->empty(); });
-
-    return {CellParticleIterator(first, cell + n),
-            CellParticleIterator(cell + n)};
-  }
-
-  Cell **begin() { return cell; }
-  Cell **end() { return cell + n; }
-
-  Cell *operator[](int i) { return assert(i < n), cell[i]; }
+namespace Cells {
+inline ParticleRange particles(Utils::Span<Cell *> cells) {
+  /* Find first non-empty cell */
+  auto first_non_empty = std::find_if(
+      cells.begin(), cells.end(), [](const Cell *c) { return not c->empty(); });
 
-  Cell **cell = nullptr;
-  int n = 0;
-};
+  return {CellParticleIterator(first_non_empty, cells.end()),
+          CellParticleIterator(cells.end())};
+}
+} // namespace Cells
 
 /** Describes a cell structure / cell system. Contains information
  *  about the communication of cell contents (particles, ghosts, ...)
@@ -199,12 +190,15 @@ struct CellStructure {
   double min_range;
 
   /** Return the global local_cells */
-  CellPList local_cells() {
-    return {m_local_cells.data(), static_cast<int>(m_local_cells.size())};
+  Utils::Span<Cell *> local_cells() {
+    return {m_local_cells.data(), m_local_cells.size()};
+  }
+
+  ParticleRange local_particles() {
+    return Cells::particles(Utils::make_span(m_local_cells));
   }
-  /** Return the global ghost_cells */
-  CellPList ghost_cells() {
-    return {m_ghost_cells.data(), static_cast<int>(m_ghost_cells.size())};
+  ParticleRange ghost_particles() {
+    return Cells::particles(Utils::make_span(m_ghost_cells));
   }
 
   /** Communicator to exchange ghost particles. */

src/core/EspressoSystemInterface.cpp

@@ -33,7 +33,7 @@ void EspressoSystemInterface::gatherParticles() {
 #ifdef CUDA
   if (m_gpu) {
     if (gpu_get_global_particle_vars_pointer_host()->communication_enabled) {
-      copy_part_data_to_gpu(cell_structure.local_cells().particles());
+      copy_part_data_to_gpu(cell_structure.local_particles());
       reallocDeviceMemory(gpu_get_particle_pointer().size());
       if (m_splitParticleStructGpu && (this_node == 0))
         split_particle_struct();

src/core/cells.cpp

@@ -198,7 +198,7 @@ unsigned topology_check_resort(int cs, unsigned local_resort) {
 /** Go through ghost cells and remove the ghost entries from the
     local particle index. */
 static void invalidate_ghosts() {
-  for (auto const &p : cell_structure.ghost_cells().particles()) {
+  for (auto const &p : cell_structure.ghost_particles()) {
     if (cell_structure.get_local_particle(p.identity()) == &p) {
       cell_structure.update_particle_index(p.identity(), nullptr);
     }
@@ -233,9 +233,7 @@ void cells_re_init(int new_cs, double range) {
 
   clear_particle_node();
 
-  for (auto &p : CellPList{old_local_cells.data(),
-                           static_cast<int>(old_local_cells.size())}
-                     .particles()) {
+  for (auto &p : Cells::particles(Utils::make_span(old_local_cells))) {
     cell_structure.add_particle(std::move(p));
   }
 
@@ -281,7 +279,8 @@ void fold_and_reset(Particle &p) {
  *
  * @returns List of Particles that do not belong on this node.
  */
-ParticleList sort_and_fold_parts(const CellStructure &cs, CellPList cells,
+ParticleList sort_and_fold_parts(const CellStructure &cs,
+                                 Utils::Span<Cell *> cells,
                                  std::vector<Cell *> &modified_cells) {
   ParticleList displaced_parts;
 
@@ -379,7 +378,7 @@ void cells_resort_particles(int global_flag) {
 
   displaced_parts.clear();
 
-  on_resort_particles(cell_structure.local_cells().particles());
+  on_resort_particles(cell_structure.local_particles());
 }
 
 /*************************************************/
@@ -404,14 +403,13 @@ void cells_on_geometry_change(int flags) {
 void check_resort_particles() {
   const double skin2 = Utils::sqr(skin / 2.0);
 
-  auto const level =
-      (std::any_of(cell_structure.local_cells().particles().begin(),
-                   cell_structure.local_cells().particles().end(),
-                   [&skin2](Particle const &p) {
-                     return (p.r.p - p.l.p_old).norm2() > skin2;
-                   }))
-          ? Cells::RESORT_LOCAL
-          : Cells::RESORT_NONE;
+  auto const level = (std::any_of(cell_structure.local_particles().begin(),
+                                  cell_structure.local_particles().end(),
+                                  [&skin2](Particle const &p) {
+                                    return (p.r.p - p.l.p_old).norm2() > skin2;
+                                  }))
+                         ? Cells::RESORT_LOCAL
+                         : Cells::RESORT_NONE;
 
   cell_structure.set_resort_particles(level);
 }
@@ -439,7 +437,7 @@ void cells_update_ghosts(unsigned data_parts) {
 
     /* Add the ghost particles to the index if we don't already
      * have them. */
-    for (auto &part : cell_structure.ghost_cells().particles()) {
+    for (auto &part : cell_structure.ghost_particles()) {
       if (cell_structure.get_local_particle(part.p.identity) == nullptr) {
         cell_structure.update_particle_index(part.identity(), &part);
       }

src/core/communication.cpp

@@ -598,8 +598,7 @@ void mpi_bcast_max_mu() { mpi_call_all(calc_mu_max); }
 
 /***** GALILEI TRANSFORM AND ASSOCIATED FUNCTIONS ****/
 void mpi_kill_particle_motion_slave(int rotation) {
-  local_kill_particle_motion(rotation,
-                             cell_structure.local_cells().particles());
+  local_kill_particle_motion(rotation, cell_structure.local_particles());
   on_particle_change();
 }
 
@@ -610,7 +609,7 @@ void mpi_kill_particle_motion(int rotation) {
 }
 
 void mpi_kill_particle_forces_slave(int torque) {
-  local_kill_particle_forces(torque, cell_structure.local_cells().particles());
+  local_kill_particle_forces(torque, cell_structure.local_particles());
   on_particle_change();
 }
 

src/core/debug.cpp

@@ -29,12 +29,12 @@
 #include <stdexcept>
 
 void check_particle_consistency() {
-  auto local_cells = cell_structure.local_cells();
-  auto const cell_part_cnt = local_cells.particles().size();
+  auto local_particles = cell_structure.local_particles();
+  auto const cell_part_cnt = local_particles.size();
 
   auto const max_id = cell_structure.get_max_local_particle_id();
 
-  for (auto const &p : local_cells.particles()) {
+  for (auto const &p : local_particles) {
     auto const id = p.identity();
 
     if (id < 0 || id > max_id) {

src/core/electrostatics_magnetostatics/coulomb.cpp

@@ -185,8 +185,8 @@ void integrate_sanity_check() {
 }
 
 void update_dependent_particles() {
-  iccp3m_iteration(cell_structure.local_cells().particles(),
-                   cell_structure.ghost_cells().particles());
+  iccp3m_iteration(cell_structure.local_particles(),
+                   cell_structure.ghost_particles());
 }
 
 void on_observable_calc() {

src/core/electrostatics_magnetostatics/mdlc_correction.cpp

@@ -39,7 +39,7 @@ DLC_struct dlc_params = {1e100, 0, 0, 0, 0};
 static double mu_max;
 
 void calc_mu_max() {
-  auto local_particles = cell_structure.local_cells().particles();
+  auto local_particles = cell_structure.local_particles();
   mu_max = std::accumulate(
       local_particles.begin(), local_particles.end(), 0.0,
       [](double mu, Particle const &p) { return std::max(mu, p.p.dipm); });

src/core/electrostatics_magnetostatics/p3m-dipolar.cpp

@@ -1795,7 +1795,7 @@ void dp3m_count_magnetic_particles() {
     tot_sums[i] = 0.0;
   }
 
-  for (auto const &p : cell_structure.local_cells().particles()) {
+  for (auto const &p : cell_structure.local_particles()) {
     if (p.p.dipm != 0.0) {
       node_sums[0] += p.calc_dip().norm2();
       node_sums[1] += 1.0;

src/core/electrostatics_magnetostatics/p3m.cpp

@@ -1697,7 +1697,7 @@ void p3m_count_charged_particles() {
     tot_sums[i] = 0.0;
   }
 
-  for (auto const &p : cell_structure.local_cells().particles()) {
+  for (auto const &p : cell_structure.local_particles()) {
     if (p.p.q != 0.0) {
       node_sums[0] += 1.0;
       node_sums[1] += Utils::sqr(p.p.q);

src/core/electrostatics_magnetostatics/scafacos.cpp

@@ -86,7 +86,7 @@ int ScafacosData::update_particle_data() {
 #endif
   }
 
-  for (auto const &p : cell_structure.local_cells().particles()) {
+  for (auto const &p : cell_structure.local_particles()) {
     auto const pos = folded_position(p.r.p, box_geo);
     positions.push_back(pos[0]);
     positions.push_back(pos[1]);
@@ -113,7 +113,7 @@ void ScafacosData::update_particle_forces() const {
   if (positions.empty())
     return;
 
-  for (auto &p : cell_structure.local_cells().particles()) {
+  for (auto &p : cell_structure.local_particles()) {
     if (!dipolar()) {
       p.f.f += coulomb.prefactor * p.p.q *
                Utils::Vector3d(Utils::Span<const double>(&(fields[it]), 3));

src/core/energy.cpp

@@ -97,9 +97,9 @@ void energy_calc(double *result, const double time) {
         add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2);
       });
 
-  calc_long_range_energies(cell_structure.local_cells().particles());
+  calc_long_range_energies(cell_structure.local_particles());
 
-  auto local_parts = cell_structure.local_cells().particles();
+  auto local_parts = cell_structure.local_particles();
   Constraints::constraints.add_energy(local_parts, time, energy);
 
 #ifdef CUDA

src/core/forces.cpp

@@ -68,7 +68,7 @@ void init_forces(const ParticleRange &particles) {
   /* initialize ghost forces with zero
      set torque to zero for all and rescale quaternions
   */
-  for (auto &p : cell_structure.ghost_cells().particles()) {
+  for (auto &p : cell_structure.ghost_particles()) {
     p.f = init_ghost_force(p);
   }
 }
@@ -88,10 +88,10 @@ void force_calc(CellStructure &cell_structure) {
   prepare_local_collision_queue();
 #endif
 
-  auto particles = cell_structure.local_cells().particles();
-  auto ghost_particles = cell_structure.ghost_cells().particles();
+  auto particles = cell_structure.local_particles();
+  auto ghost_particles = cell_structure.ghost_particles();
 #ifdef ELECTROSTATICS
-  iccp3m_iteration(particles, cell_structure.ghost_cells().particles());
+  iccp3m_iteration(particles, cell_structure.ghost_particles());
 #endif
   init_forces(particles);
 

src/core/galilei.cpp

@@ -56,8 +56,8 @@ void local_kill_particle_forces(int torque, const ParticleRange &particles) {
 /* Calculate the CMS of the system */
 std::pair<Utils::Vector3d, double> local_system_CMS() {
   return boost::accumulate(
-      cell_structure.local_cells().particles(),
-      std::pair<Utils::Vector3d, double>{}, [](auto sum, const Particle &p) {
+      cell_structure.local_particles(), std::pair<Utils::Vector3d, double>{},
+      [](auto sum, const Particle &p) {
         if (not p.p.is_virtual) {
           return std::pair<Utils::Vector3d, double>{
               sum.first +
@@ -70,8 +70,8 @@ std::pair<Utils::Vector3d, double> local_system_CMS() {
 
 std::pair<Utils::Vector3d, double> local_system_CMS_velocity() {
   return boost::accumulate(
-      cell_structure.local_cells().particles(),
-      std::pair<Utils::Vector3d, double>{}, [](auto sum, const Particle &p) {
+      cell_structure.local_particles(), std::pair<Utils::Vector3d, double>{},
+      [](auto sum, const Particle &p) {
         if (not p.p.is_virtual) {
           return std::pair<Utils::Vector3d, double>{
               sum.first + p.p.mass * p.m.v, sum.second + p.p.mass};
@@ -82,7 +82,7 @@ std::pair<Utils::Vector3d, double> local_system_CMS_velocity() {
 
 /* Remove the CMS velocity */
 void local_galilei_transform(const Utils::Vector3d &cmsvel) {
-  for (auto &p : cell_structure.local_cells().particles()) {
+  for (auto &p : cell_structure.local_particles()) {
     p.m.v -= cmsvel;
   }
 }

src/core/immersed_boundary/ImmersedBoundaries.cpp

@@ -124,7 +124,7 @@ void ImmersedBoundaries::calc_volumes() {
   std::vector<double> tempVol(MaxNumIBM);
 
   // Loop over all particles on local node
-  for (auto const &p1 : cell_structure.local_cells().particles()) {
+  for (auto const &p1 : cell_structure.local_particles()) {
     // Check if particle has a BONDED_IA_IBM_TRIEL and a
     // BONDED_IA_IBM_VOLUME_CONSERVATION Basically this loops over all
     // triangles, not all particles First round to check for volume
@@ -226,7 +226,7 @@ void ImmersedBoundaries::calc_volumes() {
 /** Calculate and add the volume force to each node */
 void ImmersedBoundaries::calc_volume_force() {
   // Loop over all particles on local node
-  for (auto &p1 : cell_structure.local_cells().particles()) {
+  for (auto &p1 : cell_structure.local_particles()) {
 
     // Check if particle has a BONDED_IA_IBM_TRIEL and a
     // BONDED_IA_IBM_VOLUME_CONSERVATION. Basically this loops over all

src/core/integrate.cpp

@@ -192,7 +192,7 @@ int integrate(int n_steps, int reuse_forces) {
 
     if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
 #ifdef ROTATION
-      convert_initial_torques(cell_structure.local_cells().particles());
+      convert_initial_torques(cell_structure.local_particles());
 #endif
     }
 
@@ -216,11 +216,11 @@ int integrate(int n_steps, int reuse_forces) {
   for (int step = 0; step < n_steps; step++) {
     ESPRESSO_PROFILER_CXX_MARK_LOOP_ITERATION(integration_loop, step);
 
-    auto particles = cell_structure.local_cells().particles();
+    auto particles = cell_structure.local_particles();
 
 #ifdef BOND_CONSTRAINT
     if (n_rigidbonds)
-      save_old_pos(particles, cell_structure.ghost_cells().particles());
+      save_old_pos(particles, cell_structure.ghost_particles());
 #endif
 
     bool early_exit = integrator_step_1(particles);
@@ -245,7 +245,7 @@ int integrate(int n_steps, int reuse_forces) {
     // Communication step: distribute ghost positions
     cells_update_ghosts(global_ghost_flags());
 
-    particles = cell_structure.local_cells().particles();
+    particles = cell_structure.local_particles();
 
     force_calc(cell_structure);
 

src/core/particle_data.cpp

@@ -404,10 +404,10 @@ void mpi_who_has_slave(int, int) {
 
   sendbuf.resize(n_part);
 
-  auto end = std::transform(cell_structure.local_cells().particles().begin(),
-                            cell_structure.local_cells().particles().end(),
-                            sendbuf.data(),
-                            [](Particle const &p) { return p.p.identity; });
+  auto end =
+      std::transform(cell_structure.local_particles().begin(),
+                     cell_structure.local_particles().end(), sendbuf.data(),
+                     [](Particle const &p) { return p.p.identity; });
 
   auto npart = std::distance(sendbuf.data(), end);
   MPI_Send(sendbuf.data(), npart, MPI_INT, 0, SOME_TAG, comm_cart);
@@ -444,9 +444,7 @@ void mpi_who_has(const ParticleRange &particles) {
 /**
  * @brief Rebuild the particle index.
  */
-void build_particle_node() {
-  mpi_who_has(cell_structure.local_cells().particles());
-}
+void build_particle_node() { mpi_who_has(cell_structure.local_particles()); }
 
 /**
  *  @brief Get the mpi rank which owns the particle with id.
@@ -890,7 +888,7 @@ Particle *local_place_particle(int id, const Utils::Vector3d &pos, int _new) {
 }
 
 void local_rescale_particles(int dir, double scale) {
-  for (auto &p : cell_structure.local_cells().particles()) {
+  for (auto &p : cell_structure.local_particles()) {
     if (dir < 3)
       p.r.p[dir] *= scale;
     else {
@@ -902,7 +900,7 @@ void local_rescale_particles(int dir, double scale) {
 #ifdef EXCLUSIONS
 void local_change_exclusion(int part1, int part2, int _delete) {
   if (part1 == -1 && part2 == -1) {
-    for (auto &p : cell_structure.local_cells().particles()) {
+    for (auto &p : cell_structure.local_particles()) {
       p.el.clear();
     }
 

src/core/pressure.cpp

@@ -125,7 +125,7 @@ void pressure_calc(double *result, double *result_t, double *result_nb,
   /* rescale kinetic energy (=ideal contribution) */
   virials.data.e[0] /= (3.0 * volume * time_step * time_step);
 
-  calc_long_range_virials(cell_structure.local_cells().particles());
+  calc_long_range_virials(cell_structure.local_particles());
 
 #ifdef VIRTUAL_SITES
   {