40#include "accumulators/AutoUpdateAccumulators.hpp"
45#include "collision_detection/CollisionDetection.hpp"
57#include "system/System.hpp"
63#include <boost/mpi/collectives/all_reduce.hpp>
66#include <caliper/cali.h>
84#ifdef WALBERLA_STATIC_ASSERT
85#error "waLberla headers should not be visible to the ESPResSo core"
90volatile std::sig_atomic_t
ctrl_C = 0;
101 assert(m_protocol !=
nullptr);
108 auto &system = get_system();
109 auto &cell_structure = *system.cell_structure;
110 auto &box_geo = *system.box_geo;
112 m_protocol = std::move(protocol);
113 update_box_params(box_geo, system.get_sim_time());
114 system.propagation->recalc_forces =
true;
119 auto &system = get_system();
120 auto &cell_structure = *system.cell_structure;
121 auto &box_geo = *system.box_geo;
122 m_protocol =
nullptr;
124 system.propagation->recalc_forces =
true;
171#ifdef STOKESIAN_DYNAMICS
177 throw std::runtime_error(
"Unknown value for integ_switch");
183 for (
auto &p : cell_structure->local_particles()) {
184 used_propagations |= p.propagation();
187 used_propagations |= propagation->default_propagation;
189 used_propagations = boost::mpi::all_reduce(
::comm_cart, used_propagations,
191 propagation->used_propagations = used_propagations;
194void System::System::integrator_sanity_checks()
const {
195 auto const thermo_switch = thermostat->thermo_switch;
196 if (time_step <= 0.) {
202 <<
"The steepest descent integrator is incompatible with thermostats";
208 "currently active combination of thermostats";
227#ifdef STOKESIAN_DYNAMICS
233 if (lb.is_solver_set() and (propagation->used_propagations &
236 if (thermostat->lb ==
nullptr) {
240 if (bonded_ias->get_n_thermalized_bonds() >= 1 and
241 (thermostat->thermalized_bond ==
nullptr or
244 <<
"Thermalized bonds require the thermalized_bond thermostat";
248 for (
auto const &p : cell_structure->local_particles()) {
250 if (p.can_rotate() and not p.is_virtual() and
251 (p.propagation() & (SYSTEM_DEFAULT | ROT_EULER | ROT_LANGEVIN |
252 ROT_BROWNIAN | ROT_STOKESIAN)) == 0) {
254 <<
"Rotating particles must have a rotation propagation mode enabled";
264 if (time_step <= 0.) {
268 auto const eps =
static_cast<double>(std::numeric_limits<float>::epsilon());
269 if ((tau - time_step) / (tau + time_step) < -eps)
270 throw std::invalid_argument(method +
" tau (" + std::to_string(tau) +
271 ") must be >= MD time_step (" +
272 std::to_string(time_step) +
")");
273 auto const factor = tau / time_step;
274 if (std::fabs(std::round(factor) - factor) / factor > eps)
275 throw std::invalid_argument(method +
" tau (" + std::to_string(tau) +
276 ") must be an integer multiple of the "
278 std::to_string(time_step) +
"). Factor is " +
279 std::to_string(factor));
289 auto const tol = agrid / 1E6;
290 if ((lattice_left - geo_left).norm2() > tol or
291 (lattice_right - geo_right).norm2() > tol) {
293 <<
"left ESPResSo: [" << geo_left <<
"], "
294 <<
"left waLBerla: [" << lattice_left <<
"]"
296 <<
"right ESPResSo: [" << geo_right <<
"], "
297 <<
"right waLBerla: [" << lattice_right <<
"]"
298 <<
"\nfor method: " << method;
299 throw std::runtime_error(
300 "waLBerla and ESPResSo disagree about domain decomposition.");
308 *system.
box_geo, cell_structure.get_le_pos_offset_at_last_resort());
309 if (cell_structure.check_resort_required(offset)) {
315 auto const &propagation = *this->propagation;
316 if ((not thermostat->langevin) or ((propagation.used_propagations &
321 auto const &langevin = *thermostat->langevin;
322 auto const kT = thermostat->kT;
323 for (
auto &p : cell_structure->local_particles()) {
345 auto const kT = thermostat.kT;
346 for (
auto &p : particles) {
384#ifdef STOKESIAN_DYNAMICS
403 for (
auto &p : particles) {
438 CALI_CXX_MARK_FUNCTION;
440 auto &propagation = *this->propagation;
441#ifdef VIRTUAL_SITES_RELATIVE
442 auto const has_vs_rel = [&propagation]() {
447#ifdef BOND_CONSTRAINT
448 auto const n_rigid_bonds = bonded_ias->get_n_rigid_bonds();
452 propagation.update_default_propagation(thermostat->thermo_switch);
453 update_used_propagations();
454 on_integration_start();
463 propagation.recalc_forces)) {
465 CALI_MARK_BEGIN(
"Initial Force Calculation");
467 thermostat->lb_coupling_deactivate();
469#ifdef VIRTUAL_SITES_RELATIVE
476 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
487 CALI_MARK_END(
"Initial Force Calculation");
491 thermostat->lb_coupling_activate();
497 int n_verlet_updates = 0;
501 auto const singleton_mode =
comm_cart.size() == 1;
502 auto caught_sigint =
false;
503 auto caught_error =
false;
505 auto lb_active =
false;
506 auto ek_active =
false;
508 lb_active = lb.is_solver_set();
509 ek_active = ek.is_ready_for_propagation();
511 auto const calc_md_steps_per_tau = [
this](
double tau) {
512 return static_cast<int>(std::round(tau / time_step));
516 CALLGRIND_START_INSTRUMENTATION;
520 CALI_CXX_MARK_LOOP_BEGIN(integration_loop,
"Integration loop");
522 int integrated_steps = 0;
523 for (
int step = 0; step < n_steps; step++) {
525 CALI_CXX_MARK_LOOP_ITERATION(integration_loop, step);
528 auto particles = cell_structure->local_particles();
530#ifdef BOND_CONSTRAINT
535 lees_edwards->update_box_params(*box_geo, sim_time);
541 sim_time += time_step;
544 for (
auto &p : particles) {
557 thermostat->philox_counter_increment();
559#ifdef BOND_CONSTRAINT
566#ifdef VIRTUAL_SITES_RELATIVE
570 cell_structure->update_ghosts_and_resort_particle(
582 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
584 particles = cell_structure->local_particles();
588#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
589 if (thermostat->lb and
601 for (
auto &p : particles) {
605#ifdef BOND_CONSTRAINT
613 if (lb_active and ek_active) {
615 auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
616 auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
618 if (md_steps_per_lb_step != md_steps_per_ek_step) {
620 <<
"LB and EK are active but with different time steps.";
623 assert(propagation.lb_skipped_md_steps ==
624 propagation.ek_skipped_md_steps);
626 propagation.lb_skipped_md_steps += 1;
627 propagation.ek_skipped_md_steps += 1;
628 if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
629 propagation.lb_skipped_md_steps = 0;
630 propagation.ek_skipped_md_steps = 0;
632 lb.ghost_communication_vel();
635 }
else if (lb_active) {
636 auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
637 propagation.lb_skipped_md_steps += 1;
638 if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
639 propagation.lb_skipped_md_steps = 0;
642 }
else if (ek_active) {
643 auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
644 propagation.ek_skipped_md_steps += 1;
645 if (propagation.ek_skipped_md_steps >= md_steps_per_ek_step) {
646 propagation.ek_skipped_md_steps = 0;
650 if (lb_active and (propagation.used_propagations &
652 thermostat->lb->rng_increment();
655#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
656 if (thermostat->lb and
659 lb.ghost_communication_vel();
665#ifdef COLLISION_DETECTION
666 collision_detection->handle_collisions();
668 bond_breakage->process_queue(*
this);
679 if (singleton_mode and ctrl_C == 1) {
680 caught_sigint =
true;
686 lb.ghost_communication();
688 lees_edwards->update_box_params(*box_geo, sim_time);
690 CALI_CXX_MARK_LOOP_END(integration_loop);
694 CALLGRIND_STOP_INSTRUMENTATION;
697#ifdef VIRTUAL_SITES_RELATIVE
704 cell_structure->update_verlet_stats(n_steps, n_verlet_updates);
718 return integrated_steps;
722 bool update_accumulators) {
723 assert(n_steps >= 0);
729 if (not cell_structure->is_verlet_skin_set()) {
731 cell_structure->set_verlet_skin_heuristic();
740 if (not update_accumulators or n_steps == 0) {
741 return integrate(n_steps, reuse_forces);
744 for (
int i = 0; i < n_steps;) {
748 std::min((n_steps - i), auto_update_accumulators->next_update());
750 auto const local_retval = integrate(steps, reuse_forces);
753 std::remove_const_t<
decltype(local_retval)> global_retval;
754 boost::mpi::all_reduce(
comm_cart, local_retval, global_retval,
756 if (global_retval < 0) {
757 return global_retval;
762 (*auto_update_accumulators)(
comm_cart, steps);
772 propagation->recalc_forces =
true;
773 lees_edwards->update_box_params(*box_geo, sim_time);
@ INTEG_METHOD_STEEPEST_DESCENT
Data structures for bonded interactions.
This file contains everything related to the global cell structure / cell system.
void lees_edwards_update(double pos_offset, double shear_velocity)
Update the Lees-Edwards parameters of the box geometry for the current simulation time.
void update_box_params(BoxGeometry &box_geo, double sim_time)
Update the Lees-Edwards parameters of the box geometry for the current simulation time.
void set_protocol(std::shared_ptr< ActiveProtocol > protocol)
Set a new Lees-Edwards protocol.
void unset_protocol()
Delete the currently active Lees-Edwards protocol.
ParticleRangeFiltered< Predicate > filter(Predicate pred) const
void update_default_propagation(int thermo_switch)
bool should_propagate_with(Particle const &p, int mode) const
RAII guard for signal handling.
void thermostat_force_init()
Calculate initial particle forces from active thermostats.
void update_used_propagations()
Update the global propagation bitmask.
void set_sim_time(double value)
Set sim_time.
int integrate_with_signal_handler(int n_steps, int reuse_forces, bool update_accumulators)
int integrate(int n_steps, int reuse_forces)
Integrate equations of motion.
std::shared_ptr< Thermostat::Thermostat > thermostat
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
std::shared_ptr< IsotropicNptThermostat > npt_iso
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
int check_runtime_errors(boost::mpi::communicator const &comm)
Count runtime errors on all nodes.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
static void resort_particles_if_needed(System::System &system)
static bool integrator_step_1(ParticleRange const &particles, Propagation const &propagation, System::System &system, double time_step)
Calls the hook for propagation kernels before the force calculation.
static void integrator_step_2(ParticleRange const &particles, Propagation const &propagation, Thermostat::Thermostat const &thermostat, double time_step)
void walberla_agrid_sanity_checks(std::string method, Utils::Vector3d const &geo_left, Utils::Vector3d const &geo_right, Utils::Vector3d const &lattice_left, Utils::Vector3d const &lattice_right, double agrid)
void walberla_tau_sanity_checks(std::string method, double tau, double time_step)
Molecular dynamics integrator.
#define INTEG_ERROR_RUNTIME
#define INTEG_ERROR_SIGINT
#define INTEG_REUSE_FORCES_NEVER
recalculate forces unconditionally (mostly used for timing)
#define INTEG_REUSE_FORCES_ALWAYS
do not recalculate forces (mostly when reading checkpoints with forces)
void brownian_dynamics_rotator(BrownianThermostat const &brownian, Particle &p, double time_step, double kT)
void brownian_dynamics_propagator(BrownianThermostat const &brownian, Particle &p, double time_step, double kT)
Utils::Vector3d friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT)
Langevin thermostat for particle translational velocities.
Utils::Vector3d friction_thermo_langevin_rotation(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT)
Langevin thermostat for particle angular velocities.
void lb_tracers_propagate(CellStructure &cell_structure, LB::Solver const &lb, double time_step)
void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure, BoxGeometry const &box_geo, LocalBox const &local_box, LB::Solver &lb)
@ DATA_PART_PROPERTIES
Particle::p.
Utils::Vector3d verlet_list_offset(BoxGeometry const &box, double pos_offset_at_last_resort)
double get_shear_velocity(double time, ActiveProtocol const &protocol)
Calculation of current velocity.
double get_pos_offset(double time, ActiveProtocol const &protocol)
@ TRANS_LB_MOMENTUM_EXCHANGE
volatile std::sig_atomic_t ctrl_C
Various procedures concerning interactions between particles.
void synchronize_npt_state()
Synchronizes NpT state such as instantaneous and average pressure.
Exports for the NpT code.
void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Correction of current velocities using RATTLE algorithm.
void save_old_position(const ParticleRange &particles, const ParticleRange &ghost_particles)
copy current position
void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Propagate velocity and position while using SHAKE algorithm for bond constraint.
void vs_relative_update_particles(CellStructure &cell_structure, BoxGeometry const &box_geo)
void convert_initial_torques(const ParticleRange &particles)
Convert torques to the body-fixed frame before the integration loop.
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
bool steepest_descent_step(const ParticleRange &particles)
Steepest descent integrator.
void stokesian_dynamics_step_1(ParticleRangeStokesian const &particles, StokesianThermostat const &stokesian, double time_step, double kT)
void velocity_verlet_rotator_1(Particle &p, double time_step)
void velocity_verlet_propagator_2(Particle &p, double time_step)
Final integration step of the Velocity Verlet integrator.
void velocity_verlet_propagator_1(Particle &p, double time_step)
Propagate the velocities and positions.
void velocity_verlet_rotator_2(Particle &p, double time_step)
void velocity_verlet_npt_step_2(ParticleRangeNPT const &particles, IsotropicNptThermostat const &npt_iso, double time_step)
Final integration step of the Velocity Verlet+NpT integrator.
void velocity_verlet_npt_step_1(ParticleRangeNPT const &particles, IsotropicNptThermostat const &npt_iso, double time_step, System::System &system)
Special propagator for NpT isotropic.