ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
core/system/System.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "config/config.hpp"
21
22#include "System.hpp"
23#include "System.impl.hpp"
24
25#include "BoxGeometry.hpp"
26#include "LocalBox.hpp"
27#include "PropagationMode.hpp"
28#include "accumulators/AutoUpdateAccumulators.hpp"
34#include "collision_detection/CollisionDetection.hpp"
35#include "communication.hpp"
37#include "errorhandling.hpp"
38#include "npt.hpp"
39#include "particle_node.hpp"
40#include "thermostat.hpp"
42
43#include <utils/Vector.hpp>
45
46#include <boost/mpi/collectives/all_reduce.hpp>
47
48#include <algorithm>
49#include <cstddef>
50#include <functional>
51#include <memory>
52#include <stdexcept>
53#include <utility>
54
55namespace System {
56
57static std::shared_ptr<System> instance = System::create();
58
59std::shared_ptr<System> System::create() {
60 auto handle = std::make_shared<System>(Private());
61 handle->initialize();
62 return handle;
63}
64
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
78 collision_detection =
79 std::make_shared<CollisionDetection::CollisionDetection>();
80#endif
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>();
86 reinit_thermo = true;
87 time_step = -1.;
88 sim_time = 0.;
89 force_cap = 0.;
90 min_global_cut = INACTIVE_CUTOFF;
91}
92
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);
105#endif
106 auto_update_accumulators->bind_system(handle);
107 constraints->bind_system(handle);
108#ifdef CUDA
109 gpu.bind_system(handle);
110 gpu.initialize();
111#endif
112 lb.bind_system(handle);
113 ek.bind_system(handle);
114}
115
116void reset_system() { instance.reset(); }
117
118void set_system(std::shared_ptr<System> new_instance) {
119 instance = new_instance;
120}
121
123
124void System::set_time_step(double value) {
125 if (value <= 0.)
126 throw std::domain_error("time_step must be > 0.");
127 if (lb.is_solver_set()) {
128 lb.veto_time_step(value);
129 }
130 if (ek.is_solver_set()) {
131 ek.veto_time_step(value);
132 }
133 time_step = value;
134 on_timestep_change();
135}
136
137void System::check_kT(double value) const {
138 if (lb.is_solver_set()) {
139 lb.veto_kT(value);
140 }
141 if (ek.is_solver_set()) {
142 ek.veto_kT(value);
143 }
144}
145
146void System::set_force_cap(double value) {
147 force_cap = value;
148 propagation->recalc_forces = true;
149}
150
151void System::set_min_global_cut(double value) {
152 min_global_cut = value;
153 on_verlet_skin_change();
154}
155
157 if (topology == CellStructureType::REGULAR) {
158 if (cell_structure->decomposition_type() == CellStructureType::REGULAR) {
159 // get fully connected info from exising regular decomposition
160 auto &old_regular_decomposition =
161 dynamic_cast<RegularDecomposition const &>(
162 std::as_const(*cell_structure).decomposition());
163 cell_structure->set_regular_decomposition(
164 get_interaction_range(),
165 old_regular_decomposition.fully_connected_boundary());
166 } else { // prev. decomposition is not a regular decomposition
167 cell_structure->set_regular_decomposition(get_interaction_range(), {});
168 }
169 } else if (topology == CellStructureType::NSQUARE) {
170 cell_structure->set_atom_decomposition();
171 } else {
172 assert(topology == CellStructureType::HYBRID);
173 /* Get current HybridDecomposition to extract n_square_types */
174 auto &old_hybrid_decomposition = dynamic_cast<HybridDecomposition const &>(
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());
179 }
180}
181
183 set_cell_structure_topology(cell_structure->decomposition_type());
184}
185
186void System::on_boxl_change(bool skip_method_adaption) {
187 update_local_geo();
188 rebuild_cell_structure();
189
190 /* Now give methods a chance to react to the change in box length */
191 if (not skip_method_adaption) {
192 lb.on_boxl_change();
193 ek.on_boxl_change();
194#ifdef ELECTROSTATICS
195 coulomb.on_boxl_change();
196#endif
197#ifdef DIPOLES
198 dipoles.on_boxl_change();
199#endif
200 }
201 constraints->on_boxl_change();
202}
203
204void System::veto_boxl_change(bool skip_particle_checks) const {
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<>());
208 if (n_part > 0ul) {
209 throw std::runtime_error(
210 "Cannot reset the box length when particles are present");
211 }
212 }
213 constraints->veto_boxl_change();
214 lb.veto_boxl_change();
215 ek.veto_boxl_change();
216}
217
219 update_local_geo();
220 lb.on_node_grid_change();
221 ek.on_node_grid_change();
222#ifdef ELECTROSTATICS
223 coulomb.on_node_grid_change();
224#endif
225#ifdef DIPOLES
226 dipoles.on_node_grid_change();
227#endif
228 rebuild_cell_structure();
229}
230
232#ifdef ELECTROSTATICS
233 coulomb.on_periodicity_change();
234#endif
235
236#ifdef DIPOLES
237 dipoles.on_periodicity_change();
238#endif
239
240#ifdef STOKESIAN_DYNAMICS
241 if (propagation->integ_switch == INTEG_METHOD_SD) {
242 if (box_geo->periodic(0u) or box_geo->periodic(1u) or box_geo->periodic(2u))
243 runtimeErrorMsg() << "Stokesian Dynamics requires periodicity "
244 << "(False, False, False)\n";
245 }
246#endif
247 on_verlet_skin_change();
248}
249
252 lb.on_cell_structure_change();
253 ek.on_cell_structure_change();
254#ifdef ELECTROSTATICS
255 coulomb.on_cell_structure_change();
256#endif
257#ifdef DIPOLES
258 dipoles.on_cell_structure_change();
259#endif
260}
261
262void System::on_thermostat_param_change() { reinit_thermo = true; }
263
265 rebuild_cell_structure();
266#ifdef ELECTROSTATICS
267 coulomb.on_coulomb_change();
268#endif
269#ifdef DIPOLES
270 dipoles.on_dipoles_change();
271#endif
272 on_short_range_ia_change();
273}
274
276 lb.on_temperature_change();
277 ek.on_temperature_change();
278}
279
281 lb.on_timestep_change();
282 ek.on_timestep_change();
283 on_thermostat_param_change();
284}
285
287 rebuild_cell_structure();
288 propagation->recalc_forces = true;
289}
290
292 nonbonded_ias->recalc_maximal_cutoffs();
293 rebuild_cell_structure();
294 propagation->recalc_forces = true;
295}
296
298#ifdef ELECTROSTATICS
299 coulomb.on_coulomb_change();
300#endif
301 on_short_range_ia_change();
302}
303
305#ifdef DIPOLES
306 dipoles.on_dipoles_change();
307#endif
308 on_short_range_ia_change();
309}
310
311void System::on_constraint_change() { propagation->recalc_forces = true; }
312
314 propagation->recalc_forces = true;
315}
316
318 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
319 propagation->recalc_forces = true;
320}
321
323 if (cell_structure->decomposition_type() == CellStructureType::HYBRID) {
324 cell_structure->set_resort_particles(Cells::RESORT_GLOBAL);
325 } else {
326 cell_structure->set_resort_particles(Cells::RESORT_LOCAL);
327 }
328#ifdef ELECTROSTATICS
329 coulomb.on_particle_change();
330#endif
331#ifdef DIPOLES
332 dipoles.on_particle_change();
333#endif
334 propagation->recalc_forces = true;
335
336 /* the particle information is no longer valid */
338}
339
341#ifdef ELECTROSTATICS
342 coulomb.on_particle_change();
343#endif
344}
345
347#ifdef VIRTUAL_SITES
348#ifdef VIRTUAL_SITES_RELATIVE
349 vs_relative_update_particles(*cell_structure, *box_geo);
350#endif
351 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
352#endif
353
354#ifdef ELECTROSTATICS
355 update_icc_particles();
356#endif
357
358 // Here we initialize volume conservation
359 // This function checks if the reference volumes have been set and if
360 // necessary calculates them
361 immersed_boundaries->init_volume_conservation(*cell_structure);
362}
363
365 /* Prepare particle structure: Communication step: number of ghosts and ghost
366 * information */
367 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
368 update_dependent_particles();
369
370#ifdef ELECTROSTATICS
371 coulomb.on_observable_calc();
372#endif
373
374#ifdef DIPOLES
375 dipoles.on_observable_calc();
376#endif
377
379}
380
381void System::on_lees_edwards_change() { lb.on_lees_edwards_change(); }
382
388
390 auto max_cut = INACTIVE_CUTOFF;
391 max_cut = std::max(max_cut, get_min_global_cut());
392#ifdef ELECTROSTATICS
393 max_cut = std::max(max_cut, coulomb.cutoff());
394#endif
395#ifdef DIPOLES
396 max_cut = std::max(max_cut, dipoles.cutoff());
397#endif
398 if (::communicator.size > 1) {
399 // If there is just one node, the bonded cutoff can be omitted
400 // because bond partners are always on the local node.
401 max_cut = std::max(max_cut, bonded_ias->maximal_cutoff());
402 }
403 max_cut = std::max(max_cut, nonbonded_ias->maximal_cutoff());
404
405#ifdef COLLISION_DETECTION
406 max_cut = std::max(max_cut, collision_detection->cutoff());
407#endif
408 return max_cut;
409}
410
412 try {
413#ifdef ELECTROSTATICS
414 coulomb.sanity_checks();
415#endif
416#ifdef DIPOLES
417 dipoles.sanity_checks();
418#endif
419 } catch (std::runtime_error const &err) {
420 runtimeErrorMsg() << err.what();
421 return true;
422 }
423 return false;
424}
425
427 auto const max_cut = maximal_cutoff();
428 auto const verlet_skin = cell_structure->get_verlet_skin();
429 /* Consider skin only if there are actually interactions */
430 return (max_cut > 0.) ? max_cut + verlet_skin : INACTIVE_CUTOFF;
431}
432
434 box_geo->set_length(box_l);
435 on_boxl_change();
436}
437
439 // sanity checks
440 integrator_sanity_checks();
441#ifdef NPT
443#endif
444 long_range_interactions_sanity_checks();
445 lb.sanity_checks();
446 ek.sanity_checks();
447
448 /* Prepare the thermostat */
449 if (reinit_thermo) {
450 thermostat->recalc_prefactors(time_step);
451 reinit_thermo = false;
452 propagation->recalc_forces = true;
453 }
454
455#ifdef NPT
456 if (propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
457 npt_ensemble_init(box_geo->length(), propagation->recalc_forces);
458 }
459#endif
460
462
463#ifdef ADDITIONAL_CHECKS
464 if (!Utils::Mpi::all_compare(::comm_cart, cell_structure->use_verlet_list)) {
465 runtimeErrorMsg() << "Nodes disagree about use of verlet lists.";
466 }
467#ifdef ELECTROSTATICS
468 {
469 auto const &actor = coulomb.impl->solver;
470 if (not Utils::Mpi::all_compare(::comm_cart, static_cast<bool>(actor)) or
471 (actor and not Utils::Mpi::all_compare(::comm_cart, (*actor).index())))
472 runtimeErrorMsg() << "Nodes disagree about Coulomb long-range method";
473 }
474#endif
475#ifdef DIPOLES
476 {
477 auto const &actor = dipoles.impl->solver;
478 if (not Utils::Mpi::all_compare(::comm_cart, static_cast<bool>(actor)) or
479 (actor and not Utils::Mpi::all_compare(::comm_cart, (*actor).index())))
480 runtimeErrorMsg() << "Nodes disagree about dipolar long-range method";
481 }
482#endif
483#endif /* ADDITIONAL_CHECKS */
484
485 on_observable_calc();
486}
487
488/**
489 * @brief Returns the ghost flags required for running pair
490 * kernels for the global state, e.g. the force calculation.
491 * @return Required data parts;
492 */
494 /* Position and Properties are always requested. */
496
497 if (lb.is_solver_set())
498 data_parts |= Cells::DATA_PART_MOMENTUM;
499
500 if (thermostat->thermo_switch & THERMO_DPD)
501 data_parts |= Cells::DATA_PART_MOMENTUM;
502
503 if (thermostat->thermo_switch & THERMO_BOND) {
504 data_parts |= Cells::DATA_PART_MOMENTUM;
505 data_parts |= Cells::DATA_PART_BONDS;
506 }
507
508#ifdef COLLISION_DETECTION
509 if (not collision_detection->is_off()) {
510 data_parts |= Cells::DATA_PART_BONDS;
511 }
512#endif
513
514 return data_parts;
515}
516
517} // namespace System
CellStructureType
Cell structure topology.
@ NSQUARE
Atom decomposition (N-square).
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
@ INTEG_METHOD_SD
@ INTEG_METHOD_NPT_ISO
@ THERMO_BOND
@ THERMO_DPD
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)
Definition LocalBox.hpp:69
Main system class.
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.
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 set_force_cap(double value)
Set force_cap.
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_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 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 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.
System & get_system()
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)
Definition npt.cpp:112
void integrator_npt_sanity_checks()
Definition npt.cpp:122
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)
Definition relative.cpp:122
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.