Refactor NpT public interface

src/core/communication.cpp

@@ -537,7 +537,7 @@ void mpi_bcast_nptiso_geom() {
 void mpi_bcast_nptiso_geom_slave(int, int) {
   MPI_Bcast(&nptiso.geometry, 1, MPI_INT, 0, comm_cart);
   MPI_Bcast(&nptiso.dimension, 1, MPI_INT, 0, comm_cart);
-  MPI_Bcast(&nptiso.cubic_box, 1, MPI_INT, 0, comm_cart);
+  MPI_Bcast(&nptiso.cubic_box, 1, MPI_LOGICAL, 0, comm_cart);
   MPI_Bcast(&nptiso.non_const_dim, 1, MPI_INT, 0, comm_cart);
 }
 

src/core/electrostatics_magnetostatics/coulomb.cpp

@@ -334,7 +334,7 @@ void calc_long_range_force(const ParticleRange &particles) {
     if (this_node == 0) {
       p3m_gpu_add_farfield_force();
     }
-    /* there is no NPT handling here as long as we cannot compute energies.g
+    /* there is no NPT handling here as long as we cannot compute energies.
        This is checked in integrator_npt_sanity_checks() when integration
        starts. */
     break;

src/core/global.cpp

@@ -116,9 +116,6 @@ const std::unordered_map<int, Datafield> fields{
     {FIELD_NPTISO_PINST,
      {&nptiso.p_inst, Datafield::Type::DOUBLE, 1,
       "npt_p_inst"}}, /* 24 from pressure.cpp */
-    {FIELD_NPTISO_PINSTAV,
-     {&nptiso.p_inst_av, Datafield::Type::DOUBLE, 1,
-      "npt_p_inst_av"}}, /* 25 from pressure.cpp */
     {FIELD_NPTISO_PDIFF,
      {&nptiso.p_diff, Datafield::Type::DOUBLE, 1,
       "npt_p_diff"}}, /* 26 from pressure.cpp */

src/core/global.hpp

@@ -81,8 +81,6 @@ enum Fields {
   FIELD_NPTISO_PEXT,
   /** index of \ref nptiso_struct::p_inst */
   FIELD_NPTISO_PINST,
-  /** index of \ref nptiso_struct::p_inst_av */
-  FIELD_NPTISO_PINSTAV,
   /** index of \ref nptiso_struct::p_diff */
   FIELD_NPTISO_PDIFF,
   /** index of \ref nptiso_struct::piston */

src/core/integrate.cpp

@@ -120,8 +120,8 @@ void integrator_sanity_checks() {
 /************************************************************/
 
 /** @brief Calls the hook for propagation kernels before the force calculation
-    @return whether or not to stop the integration loop early.
-*/
+ *  @return whether or not to stop the integration loop early.
+ */
 bool integrator_step_1(ParticleRange &particles) {
   switch (integ_switch) {
   case INTEG_METHOD_STEEPEST_DESCENT:
@@ -164,7 +164,7 @@ void integrator_step_2(ParticleRange &particles) {
 void integrate_vv(int n_steps, int reuse_forces) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
 
-  /* Prepare the Integrator */
+  /* Prepare the integrator */
   on_integration_start();
 
   /* if any method vetoes (P3M not initialized), immediately bail out */
@@ -174,13 +174,7 @@ void integrate_vv(int n_steps, int reuse_forces) {
   /* Verlet list criterion */
 
   /* Integration Step: Preparation for first integration step:
-     Calculate forces f(t) as function of positions p(t) ( and velocities v(t) )
-   */
-  /* reuse_forces logic:
-     -1: recalculate forces unconditionally, mostly used for timing
-      0: recalculate forces if recalc_forces is set, meaning it is probably
-     necessary
-      1: do not recalculate forces. Mostly when reading checkpoints with forces
+   * Calculate forces F(t) as function of positions x(t) (and velocities v(t))
    */
   if (reuse_forces == -1 || (recalc_forces && reuse_forces != 1)) {
     ESPRESSO_PROFILER_MARK_BEGIN("Initial Force Calculation");
@@ -229,8 +223,8 @@ void integrate_vv(int n_steps, int reuse_forces) {
 #ifdef BOND_CONSTRAINT
     if (n_rigidbonds)
       save_old_pos(particles, ghost_cells.particles());
-
 #endif
+
     bool early_exit = integrator_step_1(particles);
     if (early_exit)
       break;
@@ -246,8 +240,8 @@ void integrate_vv(int n_steps, int reuse_forces) {
     }
 #endif
 
-// VIRTUAL_SITES pos (and vel for DPD) update for security reason !!!
 #ifdef VIRTUAL_SITES
+    // VIRTUAL_SITES pos (and vel for DPD) update for security reason!!!
     virtual_sites()->update(true);
 #endif
 
@@ -262,15 +256,14 @@ void integrate_vv(int n_steps, int reuse_forces) {
     virtual_sites()->after_force_calc();
 #endif
     integrator_step_2(particles);
-// SHAKE velocity updates
 #ifdef BOND_CONSTRAINT
+    // SHAKE velocity updates
     if (n_rigidbonds) {
       correct_vel_shake();
     }
 #endif
 
     // propagate one-step functionalities
-
     if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
       lb_lbfluid_propagate();
       lb_lbcoupling_propagate();
@@ -284,12 +277,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
 #endif
     }
 
-#ifdef NPT
-    if (integ_switch == INTEG_METHOD_NPT_ISO) {
-      npt_update_instantaneous_pressure();
-    }
-#endif
-
     if (check_runtime_errors(comm_cart))
       break;
 
@@ -305,8 +292,8 @@ void integrate_vv(int n_steps, int reuse_forces) {
   CALLGRIND_STOP_INSTRUMENTATION;
 #endif
 
-// VIRTUAL_SITES update vel
 #ifdef VIRTUAL_SITES
+  // VIRTUAL_SITES update vel
   virtual_sites()->update(false); // Recalc positions = false
 #endif
 
@@ -404,10 +391,10 @@ void integrate_set_nvt() {
 }
 
 #ifdef NPT
-/** Parse integrate npt_isotropic command */
-int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
-                                int ydir, int zdir, bool cubic_box) {
-  nptiso.cubic_box = 0;
+int integrate_set_npt_isotropic(double ext_pressure, double piston,
+                                bool xdir_rescale, bool ydir_rescale,
+                                bool zdir_rescale, bool cubic_box) {
+  nptiso.cubic_box = cubic_box;
   nptiso.p_ext = ext_pressure;
   nptiso.piston = piston;
 
@@ -416,50 +403,31 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
                          "this integrator!\n";
     return ES_ERROR;
   }
-  if (xdir || ydir || zdir) {
-    /* set the geometry to include rescaling specified directions only*/
-    nptiso.geometry = 0;
-    nptiso.dimension = 0;
-    nptiso.non_const_dim = -1;
-    if (xdir) {
-      nptiso.geometry = (nptiso.geometry | NPTGEOM_XDIR);
-      nptiso.dimension += 1;
-      nptiso.non_const_dim = 0;
-    }
-    if (ydir) {
-      nptiso.geometry = (nptiso.geometry | NPTGEOM_YDIR);
-      nptiso.dimension += 1;
-      nptiso.non_const_dim = 1;
-    }
-    if (zdir) {
-      nptiso.geometry = (nptiso.geometry | NPTGEOM_ZDIR);
-      nptiso.dimension += 1;
-      nptiso.non_const_dim = 2;
-    }
-  } else {
-    /* set the geometry to include rescaling in all directions; the default*/
-    nptiso.geometry = 0;
-    nptiso.geometry = (nptiso.geometry | NPTGEOM_XDIR);
-    nptiso.geometry = (nptiso.geometry | NPTGEOM_YDIR);
-    nptiso.geometry = (nptiso.geometry | NPTGEOM_ZDIR);
-    nptiso.dimension = 3;
-    nptiso.non_const_dim = 2;
+  /* set the NpT geometry */
+  nptiso.geometry = 0;
+  nptiso.dimension = 0;
+  nptiso.non_const_dim = -1;
+  if (xdir_rescale) {
+    nptiso.geometry |= NPTGEOM_XDIR;
+    nptiso.dimension += 1;
+    nptiso.non_const_dim = 0;
   }
-
-  if (cubic_box) {
-    /* enable if the volume fluctuations should also apply to dimensions which
-   are switched off by the above flags
-   and which do not contribute to the pressure (3D) / tension (2D, 1D) */
-    nptiso.cubic_box = 1;
+  if (ydir_rescale) {
+    nptiso.geometry |= NPTGEOM_YDIR;
+    nptiso.dimension += 1;
+    nptiso.non_const_dim = 1;
+  }
+  if (zdir_rescale) {
+    nptiso.geometry |= NPTGEOM_ZDIR;
+    nptiso.dimension += 1;
+    nptiso.non_const_dim = 2;
   }
 
-/* Sanity Checks */
+  /* Sanity Checks */
 #ifdef ELECTROSTATICS
   if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.prefactor > 0) {
     runtimeErrorMsg() << "WARNING: If electrostatics is being used you must "
                          "use the cubic box npt.";
-    integ_switch = INTEG_METHOD_NVT;
-    mpi_bcast_parameter(FIELD_INTEG_SWITCH);
     return ES_ERROR;
   }
 #endif
@@ -468,20 +436,14 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
   if (nptiso.dimension < 3 && !nptiso.cubic_box && dipole.prefactor > 0) {
     runtimeErrorMsg() << "WARNING: If magnetostatics is being used you must "
                          "use the cubic box npt.";
-    integ_switch = INTEG_METHOD_NVT;
-    mpi_bcast_parameter(FIELD_INTEG_SWITCH);
     return ES_ERROR;
   }
 #endif
 
   if (nptiso.dimension == 0 || nptiso.non_const_dim == -1) {
     runtimeErrorMsg() << "You must enable at least one of the x y z components "
                          "as fluctuating dimension(s) for box length motion!";
-    runtimeErrorMsg() << "Cannot proceed with npt_isotropic, reverting to nvt "
-                         "integration... \n";
-    integ_switch = INTEG_METHOD_NVT;
-    mpi_bcast_parameter(FIELD_INTEG_SWITCH);
-    return (ES_ERROR);
+    return ES_ERROR;
   }
 
   /* set integrator switch */
@@ -490,8 +452,8 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
   mpi_bcast_parameter(FIELD_NPTISO_PISTON);
   mpi_bcast_parameter(FIELD_NPTISO_PEXT);
 
-  /* broadcast npt geometry information to all nodes */
+  /* broadcast NpT geometry information to all nodes */
   mpi_bcast_nptiso_geom();
-  return (ES_OK);
+  return ES_OK;
 }
 #endif

src/core/integrate.hpp

@@ -36,7 +36,7 @@
 /************************************************************/
 /*@{*/
 
-/** Switch determining which Integrator to use. */
+/** Switch determining which integrator to use. */
 extern int integ_switch;
 
 /** incremented if a Verlet update is done, aka particle resorting. */
@@ -61,8 +61,7 @@ extern bool skin_set;
 
 /** If true, the forces will be recalculated before the next integration. */
 extern bool recalc_forces;
-/** Average number of integration steps the Verlet list has been re
-    used. */
+/** Average number of integration steps the Verlet list has been re-using. */
 extern double verlet_reuse;
 
 /** Communicate signal handling to the Python interpreter */
@@ -77,45 +76,75 @@ extern bool set_py_interrupt;
 /** check sanity of integrator params */
 void integrator_sanity_checks();
 
-/** integrate equations of motion
-    \param n_steps number of steps to integrate.
-    \param reuse_forces if nonzero, blindly trust
-    the forces still stored with the particles for the first time step.
-
-    @details This function calls two hooks for propagation kernels such as
-    velocity verlet, velocity verlet + npt box changes, and steepest_descent.
-    One hook is called before and one after the force calculation.
-    It is up to the propagation kernels to increment the simulation time.
-
-    This function propagates the system according to the choice of integrator
-    stored in @ref integ_switch. The general structure is:
-    - if reuse_forces is zero, recalculate the forces based on the current
-      state of the system
-    - Loop over the number of simulation steps:
-      -# initialization (e.g., RATTLE)
-      -# First hook for propagation kernels
-      -# Update dependent particles and properties (RATTLE, virtual sites)
-      -# Calculate forces for the current state of the system. This includes
-   forces added by the Langevin thermostat and the Lattice-Boltzmann-particle
-   coupling
-      -# Second hook for propagation kernels
-      -# Update dependent properties (Virtual sites, RATTLE)
-      -# Run single step algorithms (Lattice-Boltzmann propagation, collision
-   detection, NPT update)
-    - Final update of dependent properties and statistics/counters
-
-    High-level documentation of the integration and thermostatting schemes
-    can be found in doc/sphinx/system_setup.rst and /doc/sphinx/running.rst
-
+/** Integrate equations of motion
+ *  @param n_steps       Number of integration steps, can be zero
+ *  @param reuse_forces  Decide when to re-calculate forces:
+ *                       - -1: recalculate forces unconditionally
+ *                         (mostly used for timing)
+ *                       - 0: recalculate forces if @ref recalc_forces is set,
+ *                         meaning it is probably necessary
+ *                       - 1: do not recalculate forces (mostly when reading
+ *                         checkpoints with forces)
+ *
+ *  @details This function calls two hooks for propagation kernels such as
+ *  velocity verlet, velocity verlet + npt box changes, and steepest_descent.
+ *  One hook is called before and one after the force calculation.
+ *  It is up to the propagation kernels to increment the simulation time.
+ *
+ *  This function propagates the system according to the choice of integrator
+ *  stored in @ref integ_switch. The general structure is:
+ *  - if reuse_forces is zero, recalculate the forces based on the current
+ *    state of the system
+ *  - Loop over the number of simulation steps:
+ *    -# initialization (e.g., RATTLE)
+ *    -# First hook for propagation kernels
+ *    -# Update dependent particles and properties (RATTLE, virtual sites)
+ *    -# Calculate forces for the current state of the system. This includes
+ *       forces added by the Langevin thermostat and the
+ *       Lattice-Boltzmann-particle coupling
+ *    -# Second hook for propagation kernels
+ *    -# Update dependent properties (Virtual sites, RATTLE)
+ *    -# Run single step algorithms (Lattice-Boltzmann propagation, collision
+ *       detection, NpT update)
+ *  - Final update of dependent properties and statistics/counters
+ *
+ *  High-level documentation of the integration and thermostatting schemes
+ *  can be found in doc/sphinx/system_setup.rst and /doc/sphinx/running.rst
+ *
  */
 void integrate_vv(int n_steps, int reuse_forces);
 
 /*@}*/
 
+/** @brief Run the integration loop. Can be interrupted with Ctrl+C.
+ *
+ *  @param n_steps        Number of integration steps, can be zero
+ *  @param recalc_forces  Whether to recalculate forces
+ *  @param reuse_forces   Whether to re-use forces
+ *  @retval ES_OK on success
+ *  @retval ES_ERROR on error
+ */
 int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces);
 
+/** @brief Set the NVT integrator. */
 void integrate_set_nvt();
-int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
-                                int ydir, int zdir, bool cubic_box);
+
+/** @brief Set the NpT isotropic integrator.
+ *
+ *  @param ext_pressure  Reference pressure
+ *  @param piston        Piston mass
+ *  @param xdir_rescale  Enable box rescaling in the *x*-direction
+ *  @param ydir_rescale  Enable box rescaling in the *y*-direction
+ *  @param zdir_rescale  Enable box rescaling in the *z*-direction
+ *  @param cubic_box     Determines if the volume fluctuations should also
+ *                       apply to dimensions which are switched off by the
+ *                       above flags and which do not contribute to the
+ *                       pressure (3D) or tension (2D, 1D)
+ *  @retval ES_OK on success
+ *  @retval ES_ERROR on error
+ */
+int integrate_set_npt_isotropic(double ext_pressure, double piston,
+                                bool xdir_rescale, bool ydir_rescale,
+                                bool zdir_rescale, bool cubic_box);
 
 #endif