...
 
Commits (2)
......@@ -21,9 +21,10 @@ set(EspressoCore_SRC
immersed_boundaries.cpp
event.cpp
integrate.cpp
integrators/velocity_verlet_npt.cpp
layered.cpp
metadynamics.cpp
minimize_energy.cpp
integrators/steepest_descent.cpp
npt.cpp
nsquare.cpp
partCfg_global.cpp
......
......@@ -48,8 +48,8 @@
#include "grid_based_algorithms/lb_interpolation.hpp"
#include "grid_based_algorithms/lb_particle_coupling.hpp"
#include "integrate.hpp"
#include "integrators/steepest_descent.hpp"
#include "io/mpiio/mpiio.hpp"
#include "minimize_energy.hpp"
#include "nonbonded_interactions/nonbonded_tab.hpp"
#include "npt.hpp"
#include "partCfg_global.hpp"
......
......@@ -151,8 +151,9 @@ void on_integration_start() {
recalc_forces = true;
}
/* Ensemble preparation: NVT or NPT */
integrate_ensemble_init();
#ifdef NPT
npt_ensemble_init(box_geo);
#endif
/* Update particle and observable information for routines in statistics.cpp
*/
......@@ -465,9 +466,6 @@ void on_ghost_flags_change() {
void update_dependent_particles() {
#ifdef VIRTUAL_SITES
if (virtual_sites()->is_relative()) {
ghost_communicator(&cell_structure.update_ghost_pos_comm);
}
virtual_sites()->update();
#endif
cells_update_ghosts();
......
......@@ -54,9 +54,7 @@ void init_forces(const ParticleRange &particles) {
thermodynamic ensemble */
#ifdef NPT
/* reset virial part of instantaneous pressure */
if (integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[0] = nptiso.p_vir[1] = nptiso.p_vir[2] = 0.0;
npt_reset_instantaneous_virials();
#endif
/* initialize forces with Langevin thermostat forces
......@@ -165,10 +163,6 @@ void force_calc(CellStructure &cell_structure) {
// VIRTUAL_SITES distribute forces
#ifdef VIRTUAL_SITES
if (virtual_sites()->is_relative()) {
ghost_communicator(&cell_structure.collect_ghost_force_comm);
init_forces_ghosts(ghost_particles);
}
virtual_sites()->back_transfer_forces_and_torques();
#endif
......
......@@ -293,11 +293,7 @@ inline void add_non_bonded_pair_force(Particle &p1, Particle &p2,
/* but nothing afterwards */
/*********************************************************************/
#ifdef NPT
if (integ_switch == INTEG_METHOD_NPT_ISO) {
for (int j = 0; j < 3; j++) {
nptiso.p_vir[j] += force[j] * d[j];
}
}
npt_add_virial_contribution(force, d);
#endif
/***********************************************/
......@@ -470,10 +466,7 @@ inline void add_bonded_force(Particle *const p1) {
}
#ifdef NPT
if (integ_switch == INTEG_METHOD_NPT_ISO) {
for (int j = 0; j < 3; j++)
nptiso.p_vir[j] += force1[j] * dx[j];
}
npt_add_virial_contribution(force1, dx);
#endif
switch (type) {
......
This diff is collapsed.
......@@ -77,16 +77,36 @@ extern bool set_py_interrupt;
/** check sanity of integrator params */
void integrator_sanity_checks();
/** check sanity of npt params */
void integrator_npt_sanity_checks();
/** Initialize the used thermodynamic Ensemble (NVT or NPT) */
void integrate_ensemble_init();
/** integrate with velocity Verlet integrator.
/** 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
*/
void integrate_vv(int n_steps, int reuse_forces);
......
......@@ -19,7 +19,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "minimize_energy.hpp"
#include "integrators/steepest_descent.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include "event.hpp"
......
#ifndef INTEGRATORS_VELOCITY_VERLET_HPP
#define INTEGRATORS_VELOCITY_VERLET_HPP
#include "ParticleRange.hpp"
#include "config.hpp"
#include "particle_data.hpp"
#include "rotation.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
* v(t+0.5 \Delta t) \f]
*/
inline void velocity_verlet_propagate_vel_pos(const ParticleRange &particles) {
auto const skin2 = Utils::sqr(0.5 * skin);
for (auto &p : particles) {
#ifdef ROTATION
propagate_omega_quat_particle(p);
#endif
// Don't propagate translational degrees of freedom of vs
#ifdef VIRTUAL_SITES
if (p.p.is_virtual)
continue;
#endif
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p.p.ext_flag & COORD_FIXED(j)))
#endif
{
/* 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;
/* Propagate positions (only NVT): p(t + dt) = p(t) + dt *
* v(t+0.5*dt) */
p.r.p[j] += time_step * p.m.v[j];
}
}
/* Verlet criterion check*/
if (Utils::sqr(p.r.p[0] - p.l.p_old[0]) +
Utils::sqr(p.r.p[1] - p.l.p_old[1]) +
Utils::sqr(p.r.p[2] - p.l.p_old[2]) >
skin2)
set_resort_particles(Cells::RESORT_LOCAL);
}
}
/** Final integration step of the Velocity Verlet integrator
* \f[ v(t+\Delta t) = v(t+0.5 \Delta t) + 0.5 \Delta t f(t+\Delta t)/m \f]
*/
inline void
velocity_verlet_propagate_vel_final(const ParticleRange &particles) {
for (auto &p : particles) {
#ifdef VIRTUAL_SITES
// Virtual sites are not propagated during integration
if (p.p.is_virtual)
continue;
#endif
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p.p.ext_flag & COORD_FIXED(j))) {
#endif
/* 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
}
#endif
}
}
}
inline void velocity_verlet_step_1(const ParticleRange &particles) {
velocity_verlet_propagate_vel_pos(particles);
sim_time += time_step;
}
inline void velocity_verlet_step_2(const ParticleRange &particles) {
velocity_verlet_propagate_vel_final(particles);
#ifdef ROTATION
convert_torques_propagate_omega(particles);
#endif
}
#endif
#include "config.hpp"
#ifdef NPT
#include "ParticleRange.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include "grid.hpp"
#include "integrate.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "npt.hpp"
#include "particle_data.hpp"
#include "thermostat.hpp"
#include "utils/math/sqr.hpp"
void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) {
nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0;
for (auto &p : particles) {
#ifdef VIRTUAL_SITES
// Virtual sites are not propagated during integration
if (p.p.is_virtual)
continue;
#endif
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p.p.ext_flag & COORD_FIXED(j))) {
#endif
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;
} else
/* 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
}
#endif
}
}
}
/** Scale and communicate instantaneous NPT pressure */
void velocity_verlet_npt_finalize_p_inst() {
double p_tmp = 0.0;
int i;
/* finalize derivation of p_inst */
nptiso.p_inst = 0.0;
for (i = 0; i < 3; i++) {
if (nptiso.geometry & nptiso.nptgeom_dir[i]) {
nptiso.p_vel[i] /= Utils::sqr(time_step);
nptiso.p_inst += nptiso.p_vir[i] + nptiso.p_vel[i];
}
}
MPI_Reduce(&nptiso.p_inst, &p_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, comm_cart);
if (this_node == 0) {
nptiso.p_inst = p_tmp / (nptiso.dimension * nptiso.volume);
nptiso.p_diff = nptiso.p_diff +
(nptiso.p_inst - nptiso.p_ext) * 0.5 * time_step +
friction_thermV_nptiso(nptiso.p_diff);
}
}
void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) {
double scal[3] = {0., 0., 0.}, L_new = 0.0;
/* finalize derivation of p_inst */
velocity_verlet_npt_finalize_p_inst();
/* adjust \ref nptiso_struct::nptiso.volume; prepare pos- and
* vel-rescaling
*/
if (this_node == 0) {
nptiso.volume += nptiso.inv_piston * nptiso.p_diff * 0.5 * time_step;
scal[2] = Utils::sqr(box_geo.length()[nptiso.non_const_dim]) /
pow(nptiso.volume, 2.0 / nptiso.dimension);
nptiso.volume += nptiso.inv_piston * nptiso.p_diff * 0.5 * time_step;
if (nptiso.volume < 0.0) {
runtimeErrorMsg()
<< "your choice of piston= " << nptiso.piston << ", dt= " << time_step
<< ", p_diff= " << nptiso.p_diff
<< " just caused the volume to become negative, decrease dt";
nptiso.volume =
box_geo.length()[0] * box_geo.length()[1] * box_geo.length()[2];
scal[2] = 1;
}
L_new = pow(nptiso.volume, 1.0 / nptiso.dimension);
scal[1] = L_new / box_geo.length()[nptiso.non_const_dim];
scal[0] = 1 / scal[1];
}
MPI_Bcast(scal, 3, MPI_DOUBLE, 0, comm_cart);
/* propagate positions while rescaling positions and velocities */
for (auto &p : particles) {
#ifdef VIRTUAL_SITES
if (p.p.is_virtual)
continue;
#endif
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p.p.ext_flag & COORD_FIXED(j))) {
#endif
if (nptiso.geometry & nptiso.nptgeom_dir[j]) {
{
p.r.p[j] = scal[1] * (p.r.p[j] + scal[2] * p.m.v[j] * time_step);
p.l.p_old[j] *= scal[1];
p.m.v[j] *= scal[0];
}
} else {
p.r.p[j] += p.m.v[j] * time_step;
}
#ifdef EXTERNAL_FORCES
}
#endif
}
}
set_resort_particles(Cells::RESORT_LOCAL);
/* Apply new volume to the box-length, communicate it, and account for
* necessary adjustments to the cell geometry */
if (this_node == 0) {
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) {
new_box[i] = L_new;
}
}
box_geo.set_length(new_box);
}
MPI_Bcast(box_geo.m_length.data(), 3, MPI_DOUBLE, 0, comm_cart);
/* fast box length update */
grid_changed_box_l(box_geo);
recalc_maximal_cutoff();
cells_on_geometry_change(CELL_FLAG_FAST);
}
void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) {
#ifdef NPT
nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0;
#endif
for (auto &p : particles) {
#ifdef ROTATION
propagate_omega_quat_particle(p);
#endif
// Don't propagate translational degrees of freedom of vs
#ifdef VIRTUAL_SITES
if (p.p.is_virtual)
continue;
#endif
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p.p.ext_flag & COORD_FIXED(j)))
#endif
{
#ifdef NPT
if (integ_switch == INTEG_METHOD_NPT_ISO &&
(nptiso.geometry & nptiso.nptgeom_dir[j])) {
p.m.v[j] += p.f.f[j] * 0.5 * time_step / p.p.mass +
friction_therm0_nptiso(p.m.v[j]) / p.p.mass;
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) */
p.m.v[j] += 0.5 * time_step * p.f.f[j] / p.p.mass;
}
}
}
}
void velocity_verlet_npt_step_1(const ParticleRange &particles) {
velocity_verlet_npt_propagate_vel(particles);
velocity_verlet_npt_propagate_pos(particles);
sim_time += time_step;
}
void velocity_verlet_npt_step_2(const ParticleRange &particles) {
velocity_verlet_npt_propagate_vel_final(particles);
#ifdef ROTATION
convert_torques_propagate_omega(particles);
#endif
velocity_verlet_npt_finalize_p_inst();
}
#endif
#ifndef INTEGRATORS_VELOCITY_VERLET_NPT_HPP
#define INTEGRATORS_VELOCITY_VERLET_NPT_HPP
#include "config.hpp"
#ifdef NPT
#include "ParticleRange.hpp"
#include "particle_data.hpp"
/** Special propagator for NPT ISOTROPIC
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
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:<br>
\f[ v(t+\Delta t) = v(t+0.5 \Delta t) + 0.5 \Delta t f(t+\Delta t)/m \f] */
void velocity_verlet_npt_step_2(const ParticleRange &particles);
#endif
#endif
......@@ -17,4 +17,68 @@ You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "npt.hpp"
#include "communication.hpp"
#include "config.hpp"
#include "errorhandling.hpp"
#include "integrate.hpp"
#ifdef NPT
void synchronize_npt_state(int n_steps) {
nptiso.invalidate_p_vel = false;
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)
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 "
"dimension not yet set. this should not happen. ");
}
nptiso.volume = pow(box.length()[nptiso.non_const_dim], nptiso.dimension);
if (recalc_forces) {
nptiso.p_inst = 0.0;
nptiso.p_vir[0] = nptiso.p_vir[1] = nptiso.p_vir[2] = 0.0;
nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0;
}
}
}
void integrator_npt_sanity_checks() {
if (integ_switch == INTEG_METHOD_NPT_ISO) {
if (nptiso.piston <= 0.0) {
runtimeErrorMsg() << "npt on, but piston mass not set";
}
}
}
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)
nptiso.p_vir[0] = nptiso.p_vir[1] = nptiso.p_vir[2] = 0.0;
}
void npt_add_virial_contribution(const Utils::Vector3d &force,
const Utils::Vector3d &d) {
if (integ_switch == INTEG_METHOD_NPT_ISO) {
for (int j = 0; j < 3; j++) {
nptiso.p_vir[j] += force[j] * d[j];
}
}
}
#endif // NPT
......@@ -25,6 +25,7 @@
#ifndef NPT_H
#define NPT_H
#include "BoxGeometry.hpp"
/************************************************
* data types
************************************************/
......@@ -88,4 +89,14 @@ extern nptiso_struct nptiso;
#define NPTGEOM_YDIR 2
#define NPTGEOM_ZDIR 4
/** @brief Synchronizes npt state such as instantaneous and average pressure
* @param n_steps Number of interation steps since the last sync
*/
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);
#endif
......@@ -160,7 +160,8 @@ void app_pos_correction(const ParticleRange &particles) {
} // for i loop
}
void correct_pos_shake(const ParticleRange &particles) {
void correct_pos_shake(ParticleRange const &particles) {
cells_update_ghosts();
int repeat_, cnt = 0;
int repeat = 1;
......@@ -169,7 +170,7 @@ void correct_pos_shake(const ParticleRange &particles) {
repeat_ = 0;
compute_pos_corr_vec(&repeat_, cell_structure.local_cells().particles());
ghost_communicator(&cell_structure.collect_ghost_force_comm);
app_pos_correction(particles);
app_pos_correction(cell_structure.local_cells().particles());
/**Ghost Positions Update*/
ghost_communicator(&cell_structure.update_ghost_pos_comm);
if (this_node == 0)
......@@ -273,7 +274,9 @@ void revert_force(const ParticleRange &particles,
revert(p);
}
void correct_vel_shake(CellStructure &cell_structure) {
void correct_vel_shake() {
ghost_communicator(&cell_structure.update_ghost_pos_comm);
int repeat_, repeat = 1, cnt = 0;
/**transfer the current forces to r.p_old of the particle structure so that
velocity corrections can be stored temporarily at the f.f[3] of the particle
......
......@@ -47,7 +47,7 @@ void save_old_pos(const ParticleRange &particles,
void correct_pos_shake(const ParticleRange &particles);
/** Correction of current velocities using RATTLE algorithm*/
void correct_vel_shake(CellStructure &cell_structure);
void correct_vel_shake();
/** set the parameter for a rigid, aka RATTLE bond */
int rigid_bond_set_params(int bond_type, double d, double p_tol, double v_tol);
......
......@@ -63,13 +63,6 @@ public:
m_have_quaternion = have_quaternion;
};
bool get_have_quaternion() const { return m_have_quaternion; };
/** Is a ghost comm needed before a velocity update */
virtual bool need_ghost_comm_before_vel_update() const = 0;
/**
* @brief true if the implementation depends on the position of
* real particles.
*/
virtual bool is_relative() const = 0;
virtual ~VirtualSites() = default;
......
......@@ -32,8 +32,6 @@ class VirtualSitesInertialessTracers : public VirtualSites {
void back_transfer_forces_and_torques() const override{};
void after_force_calc() override;
void after_lb_propagation() override;
bool need_ghost_comm_before_vel_update() const override { return false; };
bool is_relative() const override { return false; }
};
#endif
......
......@@ -27,8 +27,6 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
class VirtualSitesOff : public VirtualSites {
void update(bool) const override{};
void back_transfer_forces_and_torques() const override{};
bool need_ghost_comm_before_vel_update() const override { return false; };
bool is_relative() const override { return false; }
};
#endif
......
......@@ -18,6 +18,7 @@
*/
#include "VirtualSitesRelative.hpp"
#include "forces_inline.hpp"
#include <utils/math/sqr.hpp>
#ifdef VIRTUAL_SITES_RELATIVE
......@@ -30,7 +31,12 @@
#include "rotation.hpp"
void VirtualSitesRelative::update(bool recalc_positions) const {
// Ghost update logic
if (n_nodes > 0) {
if (recalc_positions or get_have_velocity()) {
ghost_communicator(&cell_structure.update_ghost_pos_comm);
}
}
for (auto &p : local_cells.particles()) {
if (!p.p.is_virtual)
continue;
......@@ -128,6 +134,9 @@ void VirtualSitesRelative::update_vel(Particle &p) const {
// Distribute forces that have accumulated on virtual particles to the
// associated real particles
void VirtualSitesRelative::back_transfer_forces_and_torques() const {
ghost_communicator(&cell_structure.collect_ghost_force_comm);
init_forces_ghosts(cell_structure.ghost_cells().particles());
// Iterate over all the particles in the local cells
for (auto &p : local_cells.particles()) {
// We only care about virtual particles
......
......@@ -35,17 +35,12 @@ public:
void update(bool recalc_positions) const override;
/** @copydoc VirtualSites::back_transfer_forces_and_torques */
void back_transfer_forces_and_torques() const override;
/** @copydoc VirtualSites::need_ghost_comm_before_vel_update */
bool need_ghost_comm_before_vel_update() const override {
return (n_nodes > 1) && get_have_velocity();
};
/** @copydoc VirtualSites::n_pressure_contribs */
int n_pressure_contribs() const override { return 1; };
/** @copydoc VirtualSites::pressure_and_stress_tensor_contribution */
void
pressure_and_stress_tensor_contribution(double *pressure,
double *stress_tensor) const override;
bool is_relative() const override { return true; }
private:
void update_pos(Particle &p) const;
......
......@@ -20,14 +20,20 @@ from libcpp.string cimport string # import std::string as string
from libcpp.vector cimport vector # import std::vector as vector
from libcpp cimport bool as cbool
include "myconfig.pxi"
cdef extern from "config.hpp":
pass
cdef extern from "integrate.hpp" nogil:
cdef int python_integrate(int n_steps, cbool recalc_forces, int reuse_forces)
cdef void integrate_set_nvt()
cdef int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir, int ydir, int zdir, int cubic_box)
cdef extern cbool skin_set
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 inline int _integrate(int nSteps, cbool recalc_forces, int reuse_forces):
with nogil:
return python_integrate(nSteps, recalc_forces, reuse_forces)
......@@ -48,7 +54,7 @@ cdef extern from "RuntimeError.hpp" namespace "ErrorHandling":
cdef extern from "errorhandling.hpp" namespace "ErrorHandling":
cdef vector[RuntimeError]mpi_gather_runtime_errors()
cdef extern from "minimize_energy.hpp":
cdef extern from "integrators/steepest_descent.hpp":
void minimize_energy_init(const double f_max, const double gamma, const int max_steps, const double max_displacement)
cdef extern from "communication.hpp":
int mpi_minimize_energy()
......
......@@ -145,18 +145,20 @@ cdef class Integrator:
If this optional parameter is true, a cubic box is assumed.
"""
self._method = "NPT"
self._isotropic_npt_params['ext_pressure'] = ext_pressure
self._isotropic_npt_params['piston'] = piston
self._isotropic_npt_params['direction'] = direction
self._isotropic_npt_params['cubic_box'] = cubic_box
if "NPT" not in espressomd.code_info.features():
raise Exception("NPT is not compiled in")
check_type_or_throw_except(
ext_pressure, 1, float, "NPT parameter ext_pressure must be a float")
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")
if (integrate_set_npt_isotropic(ext_pressure, piston, direction[0], direction[1], direction[2], cubic_box)):
handle_errors("Encountered errors setting up the NPT integrator")
IF NPT:
self._method = "NPT"
self._isotropic_npt_params['ext_pressure'] = ext_pressure
self._isotropic_npt_params['piston'] = piston
self._isotropic_npt_params['direction'] = direction
self._isotropic_npt_params['cubic_box'] = cubic_box
check_type_or_throw_except(
ext_pressure, 1, float, "NPT parameter ext_pressure must be a float")
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")
if (integrate_set_npt_isotropic(ext_pressure, piston, direction[0], direction[1], direction[2], cubic_box)):
handle_errors(
"Encountered errors setting up the NPT integrator")
ELSE:
raise Exception("NPT not compiled in.")
......@@ -22,7 +22,7 @@ from espressomd.system cimport *
cimport numpy as np
from espressomd.utils cimport *
cdef extern from "minimize_energy.hpp":
cdef extern from "integrators/steepest_descent.hpp":
cdef void minimize_energy_init(const double f_max, const double gamma, const int max_steps, const double max_displacement)
cdef extern from "communication.hpp":
......