28#include "accumulators/AutoUpdateAccumulators.hpp"
34#include "collision_detection/CollisionDetection.hpp"
35#include "communication.hpp"
40#include "thermostat.hpp"
46#include <boost/mpi/collectives/all_reduce.hpp>
60 auto handle = std::make_shared<System>(Private());
66 box_geo = std::make_shared<BoxGeometry>();
67 local_geo = std::make_shared<LocalBox>();
68 cell_structure = std::make_shared<CellStructure>(*box_geo);
69 propagation = std::make_shared<Propagation>();
70 bonded_ias = std::make_shared<BondedInteractionsMap>();
71 thermostat = std::make_shared<Thermostat::Thermostat>();
72 nonbonded_ias = std::make_shared<InteractionsNonBonded>();
73 comfixed = std::make_shared<ComFixed>();
74 galilei = std::make_shared<Galilei>();
75 oif_global = std::make_shared<OifGlobal>();
76 immersed_boundaries = std::make_shared<ImmersedBoundaries>();
77#ifdef COLLISION_DETECTION
79 std::make_shared<CollisionDetection::CollisionDetection>();
81 bond_breakage = std::make_shared<BondBreakage::BondBreakage>();
82 lees_edwards = std::make_shared<LeesEdwards::LeesEdwards>();
83 auto_update_accumulators =
84 std::make_shared<Accumulators::AutoUpdateAccumulators>();
85 constraints = std::make_shared<Constraints::Constraints>();
93void System::initialize() {
94 auto handle = shared_from_this();
95 cell_structure->bind_system(handle);
96 lees_edwards->bind_system(handle);
97 immersed_boundaries->bind_system(handle);
98 bonded_ias->bind_system(handle);
99 thermostat->bind_system(handle);
100 nonbonded_ias->bind_system(handle);
101 oif_global->bind_system(handle);
102 immersed_boundaries->bind_system(handle);
103#ifdef COLLISION_DETECTION
104 collision_detection->bind_system(handle);
106 auto_update_accumulators->bind_system(handle);
107 constraints->bind_system(handle);
109 gpu.bind_system(handle);
112 lb.bind_system(handle);
113 ek.bind_system(handle);
126 throw std::domain_error(
"time_step must be > 0.");
127 if (lb.is_solver_set()) {
128 lb.veto_time_step(value);
130 if (ek.is_solver_set()) {
131 ek.veto_time_step(value);
134 on_timestep_change();
138 if (lb.is_solver_set()) {
141 if (ek.is_solver_set()) {
148 propagation->recalc_forces =
true;
152 min_global_cut = value;
153 on_verlet_skin_change();
160 auto &old_regular_decomposition =
162 std::as_const(*cell_structure).decomposition());
163 cell_structure->set_regular_decomposition(
164 get_interaction_range(),
165 old_regular_decomposition.fully_connected_boundary());
167 cell_structure->set_regular_decomposition(get_interaction_range(), {});
170 cell_structure->set_atom_decomposition();
175 std::as_const(*cell_structure).decomposition());
176 cell_structure->set_hybrid_decomposition(
177 old_hybrid_decomposition.get_cutoff_regular(),
178 old_hybrid_decomposition.get_n_square_types());
183 set_cell_structure_topology(cell_structure->decomposition_type());
188 rebuild_cell_structure();
191 if (not skip_method_adaption) {
195 coulomb.on_boxl_change();
198 dipoles.on_boxl_change();
201 constraints->on_boxl_change();
205 if (not skip_particle_checks) {
206 auto const n_part = boost::mpi::all_reduce(
207 ::comm_cart, cell_structure->local_particles().size(), std::plus<>());
209 throw std::runtime_error(
210 "Cannot reset the box length when particles are present");
213 constraints->veto_boxl_change();
214 lb.veto_boxl_change();
215 ek.veto_boxl_change();
220 lb.on_node_grid_change();
221 ek.on_node_grid_change();
223 coulomb.on_node_grid_change();
226 dipoles.on_node_grid_change();
228 rebuild_cell_structure();
233 coulomb.on_periodicity_change();
237 dipoles.on_periodicity_change();
240#ifdef STOKESIAN_DYNAMICS
242 if (box_geo->periodic(0u) or box_geo->periodic(1u) or box_geo->periodic(2u))
244 <<
"(False, False, False)\n";
247 on_verlet_skin_change();
252 lb.on_cell_structure_change();
253 ek.on_cell_structure_change();
255 coulomb.on_cell_structure_change();
258 dipoles.on_cell_structure_change();
265 rebuild_cell_structure();
267 coulomb.on_coulomb_change();
270 dipoles.on_dipoles_change();
272 on_short_range_ia_change();
276 lb.on_temperature_change();
277 ek.on_temperature_change();
281 lb.on_timestep_change();
282 ek.on_timestep_change();
283 on_thermostat_param_change();
287 rebuild_cell_structure();
288 propagation->recalc_forces =
true;
292 nonbonded_ias->recalc_maximal_cutoffs();
293 rebuild_cell_structure();
294 propagation->recalc_forces =
true;
299 coulomb.on_coulomb_change();
301 on_short_range_ia_change();
306 dipoles.on_dipoles_change();
308 on_short_range_ia_change();
314 propagation->recalc_forces =
true;
318 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
319 propagation->recalc_forces =
true;
329 coulomb.on_particle_change();
332 dipoles.on_particle_change();
334 propagation->recalc_forces =
true;
342 coulomb.on_particle_change();
348#ifdef VIRTUAL_SITES_RELATIVE
351 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
355 update_icc_particles();
361 immersed_boundaries->init_volume_conservation(*cell_structure);
367 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
368 update_dependent_particles();
371 coulomb.on_observable_calc();
375 dipoles.on_observable_calc();
391 max_cut = std::max(max_cut, get_min_global_cut());
393 max_cut = std::max(max_cut, coulomb.cutoff());
396 max_cut = std::max(max_cut, dipoles.cutoff());
401 max_cut = std::max(max_cut, bonded_ias->maximal_cutoff());
403 max_cut = std::max(max_cut, nonbonded_ias->maximal_cutoff());
405#ifdef COLLISION_DETECTION
406 max_cut = std::max(max_cut, collision_detection->cutoff());
414 coulomb.sanity_checks();
417 dipoles.sanity_checks();
419 }
catch (std::runtime_error
const &err) {
427 auto const max_cut = maximal_cutoff();
428 auto const verlet_skin = cell_structure->get_verlet_skin();
434 box_geo->set_length(box_l);
440 integrator_sanity_checks();
444 long_range_interactions_sanity_checks();
450 thermostat->recalc_prefactors(time_step);
451 reinit_thermo =
false;
452 propagation->recalc_forces =
true;
463#ifdef ADDITIONAL_CHECKS
469 auto const &actor = coulomb.impl->solver;
477 auto const &actor = dipoles.impl->solver;
485 on_observable_calc();
497 if (lb.is_solver_set())
508#ifdef COLLISION_DETECTION
509 if (not collision_detection->is_off()) {
CellStructureType
Cell structure topology.
@ NSQUARE
Atom decomposition (N-square).
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
Vector implementation and trait types for boost qvm interoperability.
Data structures for bonded interactions.
Hybrid decomposition cell system.
static LocalBox make_regular_decomposition(Utils::Vector3d const &box_l, Utils::Vector3i const &node_index, Utils::Vector3i const &node_grid)
void on_periodicity_change()
double get_interaction_range() const
Get the interaction range.
double maximal_cutoff() const
Calculate the maximal cutoff of all interactions.
void on_constraint_change()
Called every time a constraint is changed.
void on_cell_structure_change()
unsigned get_global_ghost_flags() const
Returns the ghost flags required for running pair kernels for the global state, e....
void on_boxl_change(bool skip_method_adaption=false)
Called when the box length has changed.
void set_box_l(Utils::Vector3d const &box_l)
Change the box dimensions.
void on_temperature_change()
void on_thermostat_param_change()
void on_lees_edwards_change()
void on_integration_start()
void set_force_cap(double value)
Set force_cap.
void on_non_bonded_ia_change()
void on_verlet_skin_change()
void update_dependent_particles()
Update particles with properties depending on other particles, namely virtual sites and ICC charges.
void on_lb_boundary_conditions_change()
Called when the LB boundary conditions change (geometry, slip velocity, or both).
void veto_boxl_change(bool skip_particle_checks=false) const
void set_time_step(double value)
Set time_step.
void on_particle_change()
Called every time a particle property changes.
bool long_range_interactions_sanity_checks() const
Check electrostatic and magnetostatic methods are properly initialized.
void on_node_grid_change()
void on_particle_charge_change()
Called every time a particle charge changes.
void on_particle_local_change()
Called every time a particle local property changes.
void on_observable_calc()
called before calculating observables, i.e.
void on_timestep_change()
void check_kT(double value) const
Veto temperature change.
static std::shared_ptr< System > create()
void set_min_global_cut(double value)
Set min_global_cut.
void set_cell_structure_topology(CellStructureType topology)
Change cell structure topology.
void on_short_range_ia_change()
void rebuild_cell_structure()
Rebuild cell lists.
Communicator communicator
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
ICC is a method that allows to take into account the influence of arbitrarily shaped dielectric inter...
@ DATA_PART_MOMENTUM
Particle::m.
@ DATA_PART_PROPERTIES
Particle::p.
@ DATA_PART_BONDS
Particle::bonds.
@ DATA_PART_POSITION
Particle::r.
static std::shared_ptr< System > instance
void set_system(std::shared_ptr< System > new_instance)
bool all_compare(boost::mpi::communicator const &comm, T const &value)
Compare values on all nodes.
constexpr double INACTIVE_CUTOFF
Cutoff for deactivated interactions.
void npt_ensemble_init(Utils::Vector3d const &box_l, bool recalc_forces)
void integrator_npt_sanity_checks()
Exports for the NpT code.
void invalidate_fetch_cache()
Invalidate the fetch cache for get_particle_data.
void clear_particle_node()
Invalidate particle_node.
Particles creation and deletion.
void vs_relative_update_particles(CellStructure &cell_structure, BoxGeometry const &box_geo)
Utils::Vector3i calc_node_index() const
Calculate the node index in the Cartesian topology.
Utils::Vector3i node_grid
Regular decomposition cell system.
Routines to thermalize the center of mass and distance of a particle pair.