Merge branch 'python' into rework_particle_access

src/core/CellStructure.hpp

@@ -34,6 +34,8 @@
 #include "bond_error.hpp"
 #include "ghosts.hpp"
 
+#include <utils/math/sqr.hpp>
+
 #include <boost/algorithm/cxx11/any_of.hpp>
 #include <boost/container/static_vector.hpp>
 #include <boost/iterator/indirect_iterator.hpp>
@@ -361,6 +363,20 @@ struct CellStructure {
    */
   void clear_resort_particles() { m_resort_particles = Cells::RESORT_NONE; }
 
+  /**
+   * @brief Check whether a particle has moved further than half the skin
+   * since the last Verlet list update, thus requiring a resort.
+   * @param particles Particles to check
+   * @param skin Skin
+   * @return Whether a resort is needed.
+   */
+  bool check_resort_required(ParticleRange const &particles, double skin) {
+    auto const lim = Utils::sqr(skin / 2.);
+    return std::any_of(
+        particles.begin(), particles.end(),
+        [lim](const auto &p) { return ((p.r.p - p.l.p_old).norm2() > lim); });
+  }
+
   /**
    * @brief Synchronize number of ghosts.
    */

src/core/EspressoSystemStandAlone.cpp

@@ -56,6 +56,13 @@ void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const {
   mpi_set_box_length(box_l);
 }
 
+void EspressoSystemStandAlone::set_node_grid(
+    Utils::Vector3i const &node_grid) const {
+  if (!head_node)
+    return;
+  mpi_set_node_grid(node_grid);
+}
+
 void EspressoSystemStandAlone::set_time_step(double time_step) const {
   if (!head_node)
     return;

src/core/EspressoSystemStandAlone.hpp

@@ -28,6 +28,7 @@ class EspressoSystemStandAlone {
 public:
   EspressoSystemStandAlone(int argc, char **argv);
   void set_box_l(Utils::Vector3d const &box_l) const;
+  void set_node_grid(Utils::Vector3i const &node_grid) const;
   void set_time_step(double time_step) const;
   void set_skin(double new_skin) const;
 

src/core/Particle.hpp

@@ -355,7 +355,7 @@ struct ParticleMomentum {
 struct ParticleLocal {
   /** is particle a ghost particle. */
   bool ghost = false;
-  /** position in the last time step before last Verlet list update. */
+  /** position from the last Verlet list update. */
   Utils::Vector3d p_old = {0, 0, 0};
   /** index of the simulation box image where the particle really sits. */
   Utils::Vector3i i = {0, 0, 0};

src/core/cells.cpp

@@ -216,13 +216,8 @@ REGISTER_CALLBACK(cells_re_init)
 void mpi_bcast_cell_structure(int cs) { mpi_call_all(cells_re_init, cs); }
 
 void check_resort_particles() {
-  const double skin2 = Utils::sqr(skin / 2.0);
-
-  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;
-                                  }))
+  auto const level = (cell_structure.check_resort_required(
+                         cell_structure.local_particles(), skin))
                          ? Cells::RESORT_LOCAL
                          : Cells::RESORT_NONE;
 

src/core/integrate.cpp

@@ -133,17 +133,24 @@ void integrator_sanity_checks() {
   }
 }
 
+static void resort_particles_if_needed(ParticleRange &particles) {
+  if (cell_structure.check_resort_required(particles, skin)) {
+    cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
+  }
+}
+
 /** @brief Calls the hook for propagation kernels before the force calculation
  *  @return whether or not to stop the integration loop early.
  */
 bool integrator_step_1(ParticleRange &particles) {
+  bool early_exit = false;
   switch (integ_switch) {
   case INTEG_METHOD_STEEPEST_DESCENT:
-    if (steepest_descent_step(particles))
-      return true; // early exit
+    early_exit = steepest_descent_step(particles);
     break;
   case INTEG_METHOD_NVT:
     velocity_verlet_step_1(particles, time_step);
+    resort_particles_if_needed(particles);
     break;
 #ifdef NPT
   case INTEG_METHOD_NPT_ISO:
@@ -157,12 +164,13 @@ bool integrator_step_1(ParticleRange &particles) {
 #ifdef STOKESIAN_DYNAMICS
   case INTEG_METHOD_SD:
     stokesian_dynamics_step_1(particles, time_step);
+    resort_particles_if_needed(particles);
     break;
 #endif // STOKESIAN_DYNAMICS
   default:
     throw std::runtime_error("Unknown value for integ_switch");
   }
-  return false;
+  return early_exit;
 }
 
 /** Calls the hook of the propagation kernels after force calculation */
@@ -182,6 +190,7 @@ void integrator_step_2(ParticleRange &particles, double kT) {
   case INTEG_METHOD_BD:
     // the Ermak-McCammon's Brownian Dynamics requires a single step
     brownian_dynamics_propagator(brownian, particles, time_step, kT);
+    resort_particles_if_needed(particles);
     break;
 #ifdef STOKESIAN_DYNAMICS
   case INTEG_METHOD_SD:

src/core/integrators/brownian_inline.hpp

@@ -25,7 +25,6 @@
 #include "config.hpp"
 
 #include "ParticleRange.hpp"
-#include "cells.hpp"
 #include "integrate.hpp"
 #include "rotation.hpp"
 #include "thermostat.hpp"
@@ -43,9 +42,6 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian,
       p.m.v = bd_drag_vel(brownian.gamma, p);
       p.r.p += bd_random_walk(brownian, p, time_step, kT);
       p.m.v += bd_random_walk_vel(brownian, p);
-      /* Verlet criterion check */
-      if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin))
-        cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
     }
 #ifdef ROTATION
     if (!p.p.rotation)

src/core/integrators/stokesian_dynamics_inline.hpp

@@ -22,20 +22,14 @@
 #include "config.hpp"
 
 #ifdef STOKESIAN_DYNAMICS
-#include "CellStructure.hpp"
 #include "ParticleRange.hpp"
-#include "cells.hpp"
 #include "communication.hpp"
 #include "integrate.hpp"
 #include "rotation.hpp"
 #include "stokesian_dynamics/sd_interface.hpp"
 
-#include <utils/Vector.hpp>
-#include <utils/math/sqr.hpp>
-
 inline void stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles,
                                                  double time_step) {
-  auto const skin2 = Utils::sqr(0.5 * skin);
 
   // Compute new (translational and rotational) velocities
   propagate_vel_pos_sd(particles, comm_cart, time_step);
@@ -46,15 +40,11 @@ inline void stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles,
     p.r.p += p.m.v * time_step;
 
     // Perform rotation
-    double norm = p.m.omega.norm();
+    auto const norm = p.m.omega.norm();
     if (norm != 0) {
-      Utils::Vector3d omega_unit = (1 / norm) * p.m.omega;
+      auto const omega_unit = (1. / norm) * p.m.omega;
       local_rotate_particle(p, omega_unit, norm * time_step);
     }
-
-    // Verlet criterion check
-    if ((p.r.p - p.l.p_old).norm2() > skin2)
-      cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
   }
 }
 

src/core/integrators/velocity_verlet_inline.hpp

@@ -24,12 +24,9 @@
 #include "CellStructure.hpp"
 #include "Particle.hpp"
 #include "ParticleRange.hpp"
-#include "cells.hpp"
 #include "integrate.hpp"
 #include "rotation.hpp"
 
-#include <utils/math/sqr.hpp>
-
 /** Propagate the velocities and positions. Integration steps before force
  *  calculation of the Velocity Verlet integrator: <br> \f[ v(t+0.5 \Delta t) =
  *  v(t) + 0.5 \Delta t f(t)/m \f] <br> \f[ p(t+\Delta t) = p(t) + \Delta t
@@ -38,7 +35,6 @@
 inline void velocity_verlet_propagate_vel_pos(const ParticleRange &particles,
                                               double time_step) {
 
-  auto const skin2 = Utils::sqr(0.5 * skin);
   for (auto &p : particles) {
 #ifdef ROTATION
     propagate_omega_quat_particle(p, time_step);
@@ -57,10 +53,6 @@ inline void velocity_verlet_propagate_vel_pos(const ParticleRange &particles,
         p.r.p[j] += time_step * p.m.v[j];
       }
     }
-
-    /* Verlet criterion check*/
-    if ((p.r.p - p.l.p_old).norm2() > skin2)
-      cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
   }
 }
 

src/core/unit_tests/CMakeLists.txt

@@ -47,6 +47,8 @@ unit_test(NAME BoxGeometry_test SRC BoxGeometry_test.cpp DEPENDS EspressoCore)
 unit_test(NAME LocalBox_test SRC LocalBox_test.cpp DEPENDS EspressoCore)
 unit_test(NAME Lattice_test SRC Lattice_test.cpp DEPENDS EspressoCore)
 unit_test(NAME lb_exceptions SRC lb_exceptions.cpp DEPENDS EspressoCore)
+unit_test(NAME Verlet_list_test SRC Verlet_list_test.cpp DEPENDS EspressoCore
+          NUM_PROC 4)
 unit_test(NAME thermostats_test SRC thermostats_test.cpp DEPENDS EspressoCore)
 unit_test(NAME random_test SRC random_test.cpp DEPENDS EspressoUtils Random123)
 unit_test(NAME BondList_test SRC BondList_test.cpp DEPENDS EspressoCore)