From 2f3a935c0567036ee6bf319d85499c8121d36a57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 15 Oct 2019 20:10:37 +0200 Subject: [PATCH 1/7] Refactor NPT interface Allow bool values in the script interface. The default value for box dimension fluctuations is now [True, True, True] instead of the counter-intuitive [0, 0, 0] that was implicitly converted to [1, 1, 1] in the core. --- src/core/integrate.cpp | 49 ++++++++++++----------------- src/core/integrate.hpp | 5 +-- src/python/espressomd/integrate.pxd | 4 ++- src/python/espressomd/integrate.pyx | 17 +++++----- 4 files changed, 36 insertions(+), 39 deletions(-) diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 95afbbe04e1..3b5b3b2d15f 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -405,8 +405,9 @@ 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) { +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 = 0; nptiso.p_ext = ext_pressure; nptiso.piston = piston; @@ -416,33 +417,23 @@ 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; + /* set the geometry to include rescaling specified directions only */ + 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 (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; } diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index d832ef8e7ff..581f0eabf78 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -115,7 +115,8 @@ void integrate_vv(int n_steps, int reuse_forces); int python_integrate(int n_steps, bool recalc_forces, bool reuse_forces); void integrate_set_nvt(); -int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir, - int ydir, int zdir, bool cubic_box); +int integrate_set_npt_isotropic(double ext_pressure, double piston, + bool xdir_rescale, bool ydir_rescale, + bool zdir_rescale, bool cubic_box); #endif diff --git a/src/python/espressomd/integrate.pxd b/src/python/espressomd/integrate.pxd index 1eff4454a4e..51903947fb5 100644 --- a/src/python/espressomd/integrate.pxd +++ b/src/python/espressomd/integrate.pxd @@ -32,7 +32,9 @@ cdef extern from "integrate.hpp" nogil: IF NPT: cdef extern from "integrate.hpp" nogil: - cdef int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir, int ydir, int zdir, int cubic_box) + cdef int integrate_set_npt_isotropic(double ext_pressure, double piston, + cbool xdir_rescale, cbool ydir_rescale, + cbool zdir_rescale, cbool cubic_box) cdef inline int _integrate(int nSteps, cbool recalc_forces, int reuse_forces): with nogil: diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx index 9d04138b2ad..b2e50285cc9 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.pyx @@ -132,8 +132,8 @@ cdef class Integrator: self._method = "NVT" integrate_set_nvt() - def set_isotropic_npt(self, ext_pressure, piston, direction=[0, 0, 0], - cubic_box=False): + def set_isotropic_npt(self, ext_pressure, piston, + direction=(True, True, True), cubic_box=False): """ Set the integration method to NPT. @@ -143,16 +143,17 @@ cdef class Integrator: The external pressure. piston : :obj:`float` The mass of the applied piston. - direction : (3,) array_like of :obj:`int`, optional + direction : (3,) array_like of :obj:`bool`, optional Select which dimensions are allowed to fluctuate by assigning - them to ``1``. In the special case where all dimensions are set - to ``0`` (default), they are all set to ``1`` in the core. + them to ``True``. cubic_box : :obj:`bool`, optional - If this optional parameter is true, a cubic box is assumed. + If this optional parameter is ``True``, a cubic box is assumed. """ IF NPT: self._method = "NPT" + if isinstance(direction, np.ndarray): + direction = list(map(int, direction)) self._isotropic_npt_params['ext_pressure'] = ext_pressure self._isotropic_npt_params['piston'] = piston self._isotropic_npt_params['direction'] = direction @@ -162,7 +163,9 @@ cdef class Integrator: check_type_or_throw_except( piston, 1, float, "NPT parameter piston must be a float") check_type_or_throw_except( - direction, 3, int, "NPT parameter direction must be an array-like of three ints") + direction, 3, int, "NPT parameter direction must be an array-like of three bools") + check_type_or_throw_except( + cubic_box, 1, int, "NPT parameter cubic_box must be a bool") if (integrate_set_npt_isotropic(ext_pressure, piston, direction[0], direction[1], direction[2], cubic_box)): handle_errors( From 277b8d1cb61eba0c6377243726dfd51ddf902aba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 15 Oct 2019 21:25:29 +0200 Subject: [PATCH 2/7] Update integrator documentation --- src/core/integrate.cpp | 24 ++--- src/core/integrate.hpp | 94 +++++++++++++------- src/core/integrators/velocity_verlet_npt.cpp | 8 +- src/core/integrators/velocity_verlet_npt.hpp | 25 +++--- src/core/npt.hpp | 44 +++++---- src/python/espressomd/integrate.pyx | 3 +- 6 files changed, 109 insertions(+), 89 deletions(-) diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 3b5b3b2d15f..d89006923a2 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -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,8 +256,8 @@ 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(); } @@ -305,8 +299,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,7 +398,6 @@ void integrate_set_nvt() { } #ifdef NPT -/** Parse integrate npt_isotropic command */ int integrate_set_npt_isotropic(double ext_pressure, double piston, bool xdir_rescale, bool ydir_rescale, bool zdir_rescale, bool cubic_box) { @@ -438,13 +431,10 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, } 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; } -/* 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 " diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index 581f0eabf78..1483443130d 100644 --- a/src/core/integrate.hpp +++ b/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,44 +76,73 @@ 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(); + +/** @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); diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 3e7a0615c3b..278e93d1401 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -49,7 +49,7 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { p.m.v[j] += 0.5 * time_step / p.p.mass * p.f.f[j] + friction_therm0_nptiso(p.m.v[j]) / p.p.mass; } else - /* Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) */ + // Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) p.m.v[j] += 0.5 * time_step * p.f.f[j] / p.p.mass; #ifdef EXTERNAL_FORCES } @@ -58,7 +58,7 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { } } -/** Scale and communicate instantaneous NPT pressure */ +/** Scale and communicate instantaneous NpT pressure */ void velocity_verlet_npt_finalize_p_inst() { double p_tmp = 0.0; int i; @@ -173,8 +173,8 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { propagate_omega_quat_particle(p); #endif -// Don't propagate translational degrees of freedom of vs #ifdef VIRTUAL_SITES + // Don't propagate translational degrees of freedom of vs if (p.p.is_virtual) continue; #endif @@ -191,7 +191,7 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { nptiso.p_vel[j] += Utils::sqr(p.m.v[j] * time_step) * p.p.mass; } else #endif - /* Propagate velocities: v(t+0.5*dt) = v(t) + 0.5*dt * a(t) */ + // Propagate velocities: v(t+0.5*dt) = v(t) + 0.5*dt * a(t) p.m.v[j] += 0.5 * time_step * p.f.f[j] / p.p.mass; } } diff --git a/src/core/integrators/velocity_verlet_npt.hpp b/src/core/integrators/velocity_verlet_npt.hpp index 1b4769590f1..b1a9a8b3b18 100644 --- a/src/core/integrators/velocity_verlet_npt.hpp +++ b/src/core/integrators/velocity_verlet_npt.hpp @@ -25,19 +25,22 @@ #include "Particle.hpp" #include "ParticleRange.hpp" -/** Special propagator for NPT ISOTROPIC - Propagate the velocities and positions. Integration steps before force - calculation of the Velocity Verlet integrator:
\f[ v(t+0.5 \Delta t) = - v(t) + 0.5 \Delta t f(t)/m \f]
\f[ p(t+\Delta t) = p(t) + \Delta t - v(t+0.5 \Delta t) \f] - - Propagate pressure, box_length (2 times) and positions, rescale - positions and velocities and check Verlet list criterion (only NPT) */ +/** Special propagator for NpT isotropic. + * Propagate the velocities and positions. Integration steps before force + * calculation of the Velocity Verlet integrator: + * \f[ v(t+0.5 \Delta t) = v(t) + 0.5 \Delta t \cdot F(t)/m \f] + * \f[ x(t+\Delta t) = x(t) + \Delta t \cdot v(t+0.5 \Delta t) \f] + * + * Propagate pressure, box_length (2 times) and positions, rescale + * positions and velocities and check Verlet list criterion (only NpT). + */ void velocity_verlet_npt_step_1(const ParticleRange &particles); -/** Final integration step of the Velocity Verlet+NPT integrator and finalize - instantaneous pressure calculation:
- \f[ v(t+\Delta t) = v(t+0.5 \Delta t) + 0.5 \Delta t f(t+\Delta t)/m \f] */ +/** Final integration step of the Velocity Verlet+NpT integrator. + * Finalize instantaneous pressure calculation: + * \f[ v(t+\Delta t) = v(t+0.5 \Delta t) + * + 0.5 \Delta t \cdot F(t+\Delta t)/m \f] + */ void velocity_verlet_npt_step_2(const ParticleRange &particles); #endif diff --git a/src/core/npt.hpp b/src/core/npt.hpp index 421096b1176..84001cba4e1 100644 --- a/src/core/npt.hpp +++ b/src/core/npt.hpp @@ -19,30 +19,24 @@ * along with this program. If not, see . */ /** \file - * Exports for the NPT code, which otherwise is really spread all over... + * Exports for the NpT code. */ #ifndef NPT_H #define NPT_H #include "BoxGeometry.hpp" -/************************************************ - * data types - ************************************************/ -/** Structure to hold all variables related to the isotropic NpT-integration - * scheme. - */ +/** Parameters of the isotropic NpT-integration scheme. */ typedef struct { /** mass of a virtual piston representing the shaken box */ double piston; - /** inverse of piston */ + /** inverse of \ref piston */ double inv_piston; - /** isotropic volume. Note that we use the term volume throughout + /** isotropic volume. Note that we use the term volume throughout, * although for a 2d or 1d system we mean Area and Length respectively */ double volume; - /** desired pressure to which the algorithm strives to */ double p_ext; /** instantaneous pressure the system currently has */ @@ -59,38 +53,42 @@ typedef struct { * in offline pressure calculations such as 'analyze p_inst' */ bool invalidate_p_vel; - /** geometry information for the npt integrator. Holds the vector - * where a positive value for dir indicates that + /** geometry information for the NpT integrator. Holds the vector + * \< dir, dir, dir \> where a positive value for dir indicates that * box movement is allowed in that direction. To check whether a - * given direction is turned on use bitwise comparison with \ref + * given direction is turned on, use bitwise comparison with \ref * nptgeom_dir */ int geometry; - /** bitwise comparison values corresponding to different directions*/ + /** bitwise comparison values corresponding to different directions */ int nptgeom_dir[3]; - /** The number of dimensions in which npt boxlength motion is coupled to - * particles */ + /** The number of dimensions in which NpT boxlength motion is coupled to + * particles */ int dimension; /** Set this flag if you want all box dimensions to be identical. Needed for - * electrostatics and magnetostatics. If the value of dimension is less than - * 3 then box length motion in one or more directions will be decoupled from - * the particle motion + * electrostatics and magnetostatics. If the value of \ref dimension is + * less than 3, then box length motion in one or more directions will + * be decoupled from the particle motion */ int cubic_box; - /** An index to one of the non_constant dimensions. handy if you just want + /** An index to one of the non-constant dimensions. Handy if you just want * the variable box_l */ int non_const_dim; } nptiso_struct; extern nptiso_struct nptiso; -/** Allowable values for nptiso.geometry */ +/** @name NpT geometry bitmasks. + * Allowed values for @ref nptiso_struct::geometry. + */ +/*@{*/ #define NPTGEOM_XDIR 1 #define NPTGEOM_YDIR 2 #define NPTGEOM_ZDIR 4 +/*@}*/ -/** @brief Synchronizes npt state such as instantaneous and average pressure - * @param n_steps Number of integration steps since the last sync +/** @brief Synchronizes NpT state such as instantaneous and average pressure + * @param n_steps Number of integration steps since the last sync */ void synchronize_npt_state(int n_steps); void npt_ensemble_init(const BoxGeometry &box); diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx index b2e50285cc9..7d825a5fd37 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.pyx @@ -147,7 +147,8 @@ cdef class Integrator: Select which dimensions are allowed to fluctuate by assigning them to ``True``. cubic_box : :obj:`bool`, optional - If this optional parameter is ``True``, a cubic box is assumed. + If ``True``, a cubic box is assumed and the value of ``direction`` + will be ignored when rescaling the box. """ IF NPT: From ea3e4e54d2035f243206bcdcbd69cdbd5eb9ebb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 15 Oct 2019 21:34:57 +0200 Subject: [PATCH 3/7] Make NpT cubic_box a boolean flag --- src/core/communication.cpp | 2 +- src/core/integrate.cpp | 6 +----- src/core/integrators/velocity_verlet_npt.cpp | 6 ++---- src/core/npt.hpp | 2 +- 4 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/core/communication.cpp b/src/core/communication.cpp index 366271bba52..6eadbf54039 100644 --- a/src/core/communication.cpp +++ b/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); } diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index d89006923a2..936ea10744e 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -401,7 +401,7 @@ void integrate_set_nvt() { 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 = 0; + nptiso.cubic_box = cubic_box; nptiso.p_ext = ext_pressure; nptiso.piston = piston; @@ -430,10 +430,6 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, nptiso.non_const_dim = 2; } - if (cubic_box) { - nptiso.cubic_box = 1; - } - /* Sanity Checks */ #ifdef ELECTROSTATICS if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.prefactor > 0) { diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 278e93d1401..ad8c09cbbec 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -44,7 +44,7 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { #ifdef EXTERNAL_FORCES if (!(p.p.ext_flag & COORD_FIXED(j))) { #endif - if ((nptiso.geometry & nptiso.nptgeom_dir[j])) { + if (nptiso.geometry & nptiso.nptgeom_dir[j]) { nptiso.p_vel[j] += Utils::sqr(p.m.v[j] * time_step) * p.p.mass; p.m.v[j] += 0.5 * time_step / p.p.mass * p.f.f[j] + friction_therm0_nptiso(p.m.v[j]) / p.p.mass; @@ -145,9 +145,7 @@ void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) { Utils::Vector3d new_box = box_geo.length(); for (int i = 0; i < 3; i++) { - if (nptiso.geometry & nptiso.nptgeom_dir[i]) { - new_box[i] = L_new; - } else if (nptiso.cubic_box) { + if (nptiso.geometry & nptiso.nptgeom_dir[i] || nptiso.cubic_box) { new_box[i] = L_new; } } diff --git a/src/core/npt.hpp b/src/core/npt.hpp index 84001cba4e1..496007f8bb3 100644 --- a/src/core/npt.hpp +++ b/src/core/npt.hpp @@ -70,7 +70,7 @@ typedef struct { * less than 3, then box length motion in one or more directions will * be decoupled from the particle motion */ - int cubic_box; + bool cubic_box; /** An index to one of the non-constant dimensions. Handy if you just want * the variable box_l */ From 09348abf8f2955f4edccce9522ccc0dc02dcfd2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 15 Oct 2019 21:53:54 +0200 Subject: [PATCH 4/7] Remove NpT parameter p_inst_av Average instantaneous pressure: instantaenous pressure divided by the number of steps in the integration (NaN if nsteps=0). Unused in the core and observables, inaccessible from the Python interface. --- src/core/global.cpp | 3 --- src/core/global.hpp | 2 -- src/core/integrate.cpp | 6 ------ src/core/npt.cpp | 16 ---------------- src/core/npt.hpp | 3 --- src/core/pressure.cpp | 1 - src/python/espressomd/globals.pxd | 1 - 7 files changed, 32 deletions(-) diff --git a/src/core/global.cpp b/src/core/global.cpp index e749e70a344..96718473543 100644 --- a/src/core/global.cpp +++ b/src/core/global.cpp @@ -116,9 +116,6 @@ const std::unordered_map 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 */ diff --git a/src/core/global.hpp b/src/core/global.hpp index 6a288c39295..c08e025e233 100644 --- a/src/core/global.hpp +++ b/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 */ diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 936ea10744e..2abd278a20b 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -278,12 +278,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; diff --git a/src/core/npt.cpp b/src/core/npt.cpp index ac9b2d42316..0798543c164 100644 --- a/src/core/npt.cpp +++ b/src/core/npt.cpp @@ -16,8 +16,6 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#include - #include "communication.hpp" #include "config.hpp" #include "errorhandling.hpp" @@ -30,21 +28,12 @@ void synchronize_npt_state(int n_steps) { MPI_Bcast(&nptiso.p_inst, 1, MPI_DOUBLE, 0, comm_cart); MPI_Bcast(&nptiso.p_diff, 1, MPI_DOUBLE, 0, comm_cart); MPI_Bcast(&nptiso.volume, 1, MPI_DOUBLE, 0, comm_cart); - if (this_node == 0) { - if (n_steps == 0) { - nptiso.p_inst_av = std::nan(""); - } else { - nptiso.p_inst_av /= 1.0 * n_steps; - } - } - MPI_Bcast(&nptiso.p_inst_av, 1, MPI_DOUBLE, 0, comm_cart); } void npt_ensemble_init(const BoxGeometry &box) { if (integ_switch == INTEG_METHOD_NPT_ISO) { /* prepare NpT-integration */ nptiso.inv_piston = 1 / (1.0 * nptiso.piston); - nptiso.p_inst_av = 0.0; if (nptiso.dimension == 0) { throw std::runtime_error( "%d: INTERNAL ERROR: npt integrator was called but " @@ -69,11 +58,6 @@ void integrator_npt_sanity_checks() { } } -void npt_update_instantaneous_pressure() { - if ((this_node == 0) && (integ_switch == INTEG_METHOD_NPT_ISO)) - nptiso.p_inst_av += nptiso.p_inst; -} - /** reset virial part of instantaneous pressure */ void npt_reset_instantaneous_virials() { if (integ_switch == INTEG_METHOD_NPT_ISO) diff --git a/src/core/npt.hpp b/src/core/npt.hpp index 496007f8bb3..745e494bfab 100644 --- a/src/core/npt.hpp +++ b/src/core/npt.hpp @@ -41,8 +41,6 @@ typedef struct { double p_ext; /** instantaneous pressure the system currently has */ double p_inst; - /** instantaneous pressure averaged over current integration cycle */ - double p_inst_av; /** difference between \ref p_ext and \ref p_inst */ double p_diff; /** virial (short-range) components of \ref p_inst */ @@ -93,7 +91,6 @@ extern nptiso_struct nptiso; void synchronize_npt_state(int n_steps); void npt_ensemble_init(const BoxGeometry &box); void integrator_npt_sanity_checks(); -void npt_update_instantaneous_pressure(); void npt_reset_instantaneous_virials(); void npt_add_virial_contribution(const Utils::Vector3d &force, const Utils::Vector3d &d); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index f5fb7e7cad3..c26ac6690bf 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -53,7 +53,6 @@ nptiso_struct nptiso = {0.0, 0.0, 0.0, 0.0, - 0.0, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, true, diff --git a/src/python/espressomd/globals.pxd b/src/python/espressomd/globals.pxd index 19a520758b1..d4f0b8b07dd 100644 --- a/src/python/espressomd/globals.pxd +++ b/src/python/espressomd/globals.pxd @@ -120,7 +120,6 @@ cdef extern from "npt.hpp": ctypedef struct nptiso_struct: double p_ext double p_inst - double p_inst_av double p_diff double piston extern nptiso_struct nptiso From 3bf96e489899470e9af6088421fb1d9e7c7166e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 09:48:09 +0200 Subject: [PATCH 5/7] Don't decay NpT integrator to NVT integrator When the NpT integrator is initialized with incorrect values, halt script execution with an exception instead of initializing NVT. --- src/core/integrate.cpp | 12 ++---------- src/python/espressomd/integrate.pyx | 3 ++- 2 files changed, 4 insertions(+), 11 deletions(-) diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 2abd278a20b..f3cf41f6f6e 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -429,8 +429,6 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, 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 @@ -439,8 +437,6 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, 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 @@ -448,11 +444,7 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, 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 */ @@ -463,6 +455,6 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, /* broadcast npt geometry information to all nodes */ mpi_bcast_nptiso_geom(); - return (ES_OK); + return ES_OK; } #endif diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx index 7d825a5fd37..562f5d80aee 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.pyx @@ -148,7 +148,8 @@ cdef class Integrator: them to ``True``. cubic_box : :obj:`bool`, optional If ``True``, a cubic box is assumed and the value of ``direction`` - will be ignored when rescaling the box. + will be ignored when rescaling the box. This is required e.g. for + electrostatics and magnetostatics. """ IF NPT: From 857d3eeab101593b6004da70d85796e1618f4896 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 10:01:57 +0200 Subject: [PATCH 6/7] Cleanup comments --- src/core/electrostatics_magnetostatics/coulomb.cpp | 2 +- src/core/integrate.cpp | 11 +++++------ src/core/npt.cpp | 2 +- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 8d8b9588cf4..79ea280267b 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/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; diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index f3cf41f6f6e..517bb524644 100644 --- a/src/core/integrate.cpp +++ b/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: @@ -241,7 +241,7 @@ void integrate_vv(int n_steps, int reuse_forces) { #endif #ifdef VIRTUAL_SITES - // VIRTUAL_SITES pos (and vel for DPD) update for security reason !!! + // VIRTUAL_SITES pos (and vel for DPD) update for security reason!!! virtual_sites()->update(true); #endif @@ -264,7 +264,6 @@ void integrate_vv(int n_steps, int reuse_forces) { #endif // propagate one-step functionalities - if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { lb_lbfluid_propagate(); lb_lbcoupling_propagate(); @@ -404,7 +403,7 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, "this integrator!\n"; return ES_ERROR; } - /* set the geometry to include rescaling specified directions only */ + /* set the NpT geometry */ nptiso.geometry = 0; nptiso.dimension = 0; nptiso.non_const_dim = -1; @@ -453,7 +452,7 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, 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; } diff --git a/src/core/npt.cpp b/src/core/npt.cpp index 0798543c164..fd1078e2246 100644 --- a/src/core/npt.cpp +++ b/src/core/npt.cpp @@ -16,11 +16,11 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ +#include "npt.hpp" #include "communication.hpp" #include "config.hpp" #include "errorhandling.hpp" #include "integrate.hpp" -#include "npt.hpp" #ifdef NPT void synchronize_npt_state(int n_steps) { From 98e274936106fd0f035d6eadc701a34a760c9563 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 16 Oct 2019 11:39:31 +0200 Subject: [PATCH 7/7] Fix regression --- src/core/pressure.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index c26ac6690bf..b41d03cc567 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -59,7 +59,7 @@ nptiso_struct nptiso = {0.0, 0, {NPTGEOM_XDIR, NPTGEOM_YDIR, NPTGEOM_ZDIR}, 0, - 0, + false, 0}; /************************************************************/