ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
script_interface/system/System.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2013-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 "System.hpp"
21
22#include "config/config.hpp"
23
24#include "core/BoxGeometry.hpp"
25#include "core/Particle.hpp"
29#include "core/cells.hpp"
31#include "core/exclusions.hpp"
35#include "core/propagation.hpp"
36#include "core/rotation.hpp"
39
58
59#include <utils/Vector.hpp>
60#include <utils/demangle.hpp>
63
64#include <boost/mpi/collectives.hpp>
65
66#include <algorithm>
67#include <array>
68#include <cassert>
69#include <cmath>
70#include <functional>
71#include <initializer_list>
72#include <memory>
73#include <stdexcept>
74#include <string>
75#include <type_traits>
76#include <vector>
77
78namespace ScriptInterface {
79namespace System {
80
81static bool system_created = false;
82
83#ifdef EXCLUSIONS
85 Variant const &exclusions) {
86 p.call_method("set_exclusions", {{"p_ids", exclusions}});
87}
88#endif // EXCLUSIONS
89
90static void set_bonds(Particles::ParticleHandle &p, Variant const &bonds) {
92 for (auto const &bond_flat : bond_list_flat) {
93 auto const bond_id = bond_flat[0];
94 auto const part_id =
95 std::vector<int>{bond_flat.begin() + 1, bond_flat.end()};
96 p.call_method("add_bond",
97 {{"bond_id", bond_id}, {"part_id", std::move(part_id)}});
98 }
99}
100
101/** @brief Container for leaves of the system class. */
103 Leaves() = default;
104 std::shared_ptr<CellSystem::CellSystem> cell_system;
105 std::shared_ptr<Integrators::IntegratorHandle> integrator;
106 std::shared_ptr<Interactions::BondedInteractions> bonded_interactions;
107#ifdef COLLISION_DETECTION
108 std::shared_ptr<CollisionDetection::CollisionDetection> collision_detection;
109#endif
110 std::shared_ptr<Thermostat::Thermostat> thermostat;
111 std::shared_ptr<Analysis::Analysis> analysis;
112 std::shared_ptr<Galilei::ComFixed> comfixed;
113 std::shared_ptr<Galilei::Galilei> galilei;
114 std::shared_ptr<BondBreakage::BreakageSpecs> bond_breakage;
115 std::shared_ptr<LeesEdwards::LeesEdwards> lees_edwards;
116 std::shared_ptr<Accumulators::AutoUpdateAccumulators>
118 std::shared_ptr<Constraints::Constraints> constraints;
119 std::shared_ptr<Interactions::NonBondedInteractions> non_bonded_inter;
120#ifdef ELECTROSTATICS
121 std::shared_ptr<Coulomb::Container> electrostatics;
122#endif
123#ifdef DIPOLES
124 std::shared_ptr<Dipoles::Container> magnetostatics;
125#endif
126 std::shared_ptr<Particles::ParticleList> part;
127};
128
129System::System() : m_instance{}, m_leaves{std::make_shared<Leaves>()} {
130 auto const add_parameter =
131 [this, ptr = m_leaves.get()](std::string key, auto(Leaves::*member)) {
133 key.c_str(),
134 [this, ptr, member, key](Variant const &val) {
135 auto &dst = ptr->*member;
136 if (dst != nullptr) {
137 throw WriteError(key);
138 }
139 dst = get_value<std::remove_reference_t<decltype(dst)>>(val);
140 dst->bind_system(m_instance);
141 },
142 [ptr, member]() { return ptr->*member; })});
143 };
144
145 add_parameters({
146 {"box_l",
147 [this](Variant const &v) {
148 context()->parallel_try_catch([&]() {
151 throw std::domain_error("Attribute 'box_l' must be > 0");
152 }
153 m_instance->veto_boxl_change();
154 m_instance->box_geo->set_length(new_value);
155 m_instance->on_boxl_change();
156 });
157 },
158 [this]() { return m_instance->box_geo->length(); }},
159 {"periodicity",
160 [this](Variant const &v) {
162 for (unsigned int i = 0u; i < 3u; ++i) {
163 m_instance->box_geo->set_periodic(i, periodicity[i]);
164 }
165 context()->parallel_try_catch(
166 [&]() { m_instance->on_periodicity_change(); });
167 },
168 [this]() {
169 return Utils::Vector3b{m_instance->box_geo->periodic(0),
170 m_instance->box_geo->periodic(1),
171 m_instance->box_geo->periodic(2)};
172 }},
173 {"min_global_cut",
174 [this](Variant const &v) {
175 context()->parallel_try_catch([&]() {
176 auto const new_value = get_value<double>(v);
178 throw std::domain_error("Attribute 'min_global_cut' must be >= 0");
179 }
180 m_instance->set_min_global_cut(new_value);
181 });
182 },
183 [this]() { return m_instance->get_min_global_cut(); }},
184 {"max_oif_objects",
185 [this](Variant const &v) {
186 m_instance->oif_global->max_oif_objects = get_value<int>(v);
187 },
188 [this]() { return m_instance->oif_global->max_oif_objects; }},
189 });
190 // note: the order of leaves matters! e.g. bonds depend on thermostats,
191 // and thus a thermostat object must be instantiated before the bonds
192 add_parameter("cell_system", &Leaves::cell_system);
193 add_parameter("integrator", &Leaves::integrator);
194 add_parameter("thermostat", &Leaves::thermostat);
195 add_parameter("analysis", &Leaves::analysis);
196 add_parameter("comfixed", &Leaves::comfixed);
197 add_parameter("galilei", &Leaves::galilei);
198 add_parameter("bonded_inter", &Leaves::bonded_interactions);
199#ifdef COLLISION_DETECTION
200 add_parameter("collision_detection", &Leaves::collision_detection);
201#endif
202 add_parameter("bond_breakage", &Leaves::bond_breakage);
203 add_parameter("lees_edwards", &Leaves::lees_edwards);
204 add_parameter("auto_update_accumulators", &Leaves::auto_update_accumulators);
205 add_parameter("constraints", &Leaves::constraints);
206 add_parameter("non_bonded_inter", &Leaves::non_bonded_inter);
207#ifdef ELECTROSTATICS
208 add_parameter("electrostatics", &Leaves::electrostatics);
209#endif
210#ifdef DIPOLES
211 add_parameter("magnetostatics", &Leaves::magnetostatics);
212#endif
213 add_parameter("part", &Leaves::part);
214}
215
216template <typename LeafType>
217void System::do_set_default_parameter(std::string const &name) {
218 assert(context()->is_head_node());
219 auto const so_name = Utils::demangle<LeafType>().substr(17);
220 set_parameter(name, Variant{context()->make_shared(so_name, {})});
221}
222
223void System::do_construct(VariantMap const &params) {
224 /* When reloading the system state from a checkpoint file,
225 * the order of global variables instantiation matters.
226 * The @c box_l must be set before any other global variable.
227 * All these globals re-initialize the cell system, and we
228 * cannot use the default-constructed @c box_geo when e.g.
229 * long-range interactions exist in the system, otherwise
230 * runtime errors about the local geometry being smaller
231 * than the interaction range would be raised.
232 */
233 context()->parallel_try_catch([&]() {
234 if (not params.contains("box_l")) {
235 throw std::domain_error("Required argument 'box_l' not provided.");
236 }
237 if (params.contains("_regular_constructor") and system_created) {
238 throw std::runtime_error(
239 "You can only have one instance of the system class at a time");
240 }
241 });
242 m_instance = ::System::System::create();
243 ::System::set_system(m_instance);
244
245 // domain decomposition can only be set after box_l is set
246 m_instance->set_cell_structure_topology(CellStructureType::NSQUARE);
247 do_set_parameter("box_l", params.at("box_l"));
248 m_instance->set_cell_structure_topology(CellStructureType::REGULAR);
249
250 m_instance->lb.bind_system(m_instance);
251 m_instance->ek.bind_system(m_instance);
252
253 if (params.contains("_regular_constructor")) {
254 std::set<std::string> const setable_properties = {
255 "box_l", "min_global_cut",
256 "periodicity", "time",
257 "time_step", "force_cap",
258 "max_oif_objects", "_regular_constructor"};
259 for (auto const &kv : params) {
260 if (not setable_properties.contains(kv.first)) {
261 context()->parallel_try_catch([&kv]() {
262 throw std::domain_error(
263 "Property '" + kv.first +
264 "' cannot be set via argument to System class");
265 });
266 }
267 }
268 for (std::string attr :
269 {"min_global_cut", "periodicity", "max_oif_objects"}) {
270 if (params.contains(attr)) {
271 do_set_parameter(attr, params.at(attr));
272 }
273 }
274 if (not context()->is_head_node()) {
275 return;
276 }
277 auto integrator = std::dynamic_pointer_cast<Integrators::IntegratorHandle>(
278 context()->make_shared("Integrators::IntegratorHandle", {}));
279 set_parameter("integrator", integrator);
280 for (std::string attr : {"time", "time_step", "force_cap"}) {
281 if (params.contains(attr)) {
282 integrator->set_parameter(attr, params.at(attr));
283 }
284 }
285 // note: the order of leaves matters! e.g. bonds depend on thermostats,
286 // and thus a thermostat object must be instantiated before the bonds
290#ifdef COLLISION_DETECTION
292 "collision_detection");
293#endif
300 "auto_update_accumulators");
303 "non_bonded_inter");
304#ifdef ELECTROSTATICS
306#endif
307#ifdef DIPOLES
309#endif
311 } else {
312 for (auto const &key : get_parameter_insertion_order()) {
313 if (key != "box_l" and params.contains(key)) {
314 do_set_parameter(key, params.at(key));
315 }
316 }
317 }
318 if (not context()->is_head_node()) {
319 return;
320 }
321 call_method("internal_attach_leaves", {});
322}
323
324static void rotate_system(CellStructure &cell_structure, double phi,
325 double theta, double alpha) {
326 auto const particles = cell_structure.local_particles();
327
328 // Calculate center of mass
330 double local_mass = 0.0;
331
332 for (auto const &p : particles) {
333 if (not p.is_virtual()) {
334 local_com += p.mass() * p.pos();
335 local_mass += p.mass();
336 }
337 }
338
339 auto const total_mass =
340 boost::mpi::all_reduce(comm_cart, local_mass, std::plus<>());
341 auto const com =
342 boost::mpi::all_reduce(comm_cart, local_com, std::plus<>()) / total_mass;
343
344 // Rotation axis in Cartesian coordinates
345 Utils::Vector3d axis;
346 axis[0] = std::sin(theta) * std::cos(phi);
347 axis[1] = std::sin(theta) * std::sin(phi);
348 axis[2] = std::cos(theta);
349
350 // Rotate particle coordinates
351 for (auto &p : particles) {
352 // Move the center of mass of the system to the origin
353 p.pos() = com + Utils::vec_rotate(axis, alpha, p.pos() - com);
354#ifdef ROTATION
355 local_rotate_particle(p, axis, alpha);
356#endif
357 }
358
362}
363
364/** Rescale all particle positions in direction @p dir by a factor @p scale.
365 * @param cell_structure cell structure
366 * @param dir direction to scale (0/1/2 = x/y/z, 3 = x+y+z isotropically)
367 * @param scale factor by which to rescale (>1: stretch, <1: contract)
368 */
369static void rescale_particles(CellStructure &cell_structure, int dir,
370 double scale) {
371 for (auto &p : cell_structure.local_particles()) {
372 if (dir < 3)
373 p.pos()[dir] *= scale;
374 else {
375 p.pos() *= scale;
376 }
377 }
378}
379
380Variant System::do_call_method(std::string const &name,
381 VariantMap const &parameters) {
382 if (name == "lock_system_creation") {
383 system_created = true;
384 return {};
385 }
386 if (name == "rescale_boxl") {
387 auto &box_geo = *m_instance->box_geo;
388 auto const coord = get_value<int>(parameters, "coord");
389 auto const length = get_value<double>(parameters, "length");
390 assert(coord >= 0);
391 assert(coord != 3 or ((box_geo.length()[0] == box_geo.length()[1]) and
392 (box_geo.length()[1] == box_geo.length()[2])));
393 auto const scale = (coord == 3) ? length * box_geo.length_inv()[0]
394 : length * box_geo.length_inv()[coord];
395 context()->parallel_try_catch([&]() {
396 if (length <= 0.) {
397 throw std::domain_error("Parameter 'd_new' must be > 0");
398 }
399 m_instance->veto_boxl_change(true);
400 });
401 auto new_value = Utils::Vector3d{};
402 if (coord == 3) {
404 } else {
405 new_value = box_geo.length();
406 new_value[static_cast<unsigned>(coord)] = length;
407 }
408 // when shrinking, rescale the particles first
409 if (scale <= 1.) {
410 rescale_particles(*m_instance->cell_structure, coord, scale);
411 m_instance->on_particle_change();
412 }
413 m_instance->box_geo->set_length(new_value);
414 m_instance->on_boxl_change();
415 if (scale > 1.) {
416 rescale_particles(*m_instance->cell_structure, coord, scale);
417 m_instance->on_particle_change();
418 }
419 return {};
420 }
421 if (name == "setup_type_map") {
422 auto const types = get_value<std::vector<int>>(parameters, "type_list");
423 for (auto const type : types) {
424 ::init_type_map(type);
425 }
426 return {};
427 }
428 if (name == "number_of_particles") {
429 auto const type = get_value<int>(parameters, "type");
430 return ::number_of_particles_with_type(type);
431 }
432 if (name == "velocity_difference") {
433 auto const pos1 = get_value<Utils::Vector3d>(parameters, "pos1");
434 auto const pos2 = get_value<Utils::Vector3d>(parameters, "pos2");
435 auto const v1 = get_value<Utils::Vector3d>(parameters, "v1");
436 auto const v2 = get_value<Utils::Vector3d>(parameters, "v2");
437 return m_instance->box_geo->velocity_difference(pos2, pos1, v2, v1);
438 }
439 if (name == "distance_vec") {
440 auto const pos1 = get_value<Utils::Vector3d>(parameters, "pos1");
441 auto const pos2 = get_value<Utils::Vector3d>(parameters, "pos2");
442 return m_instance->box_geo->get_mi_vector(pos2, pos1);
443 }
444 if (name == "rotate_system") {
445 rotate_system(*m_instance->cell_structure,
448 get_value<double>(parameters, "alpha"));
449 m_instance->on_particle_change();
450 m_instance->update_dependent_particles();
451 return {};
452 }
453 if (name == "get_propagation_modes_enum") {
455 }
456 if (name == "session_shutdown") {
457 if (m_instance) {
458 if (&::System::get_system() == m_instance.get()) {
460 }
461 assert(m_instance.use_count() == 1l);
462 m_instance.reset();
463 }
464 return {};
465 }
466 if (name == "internal_attach_leaves") {
467 m_leaves->part->attach(m_leaves->cell_system,
468 m_leaves->bonded_interactions);
469#ifdef COLLISION_DETECTION
470 m_leaves->collision_detection->attach(m_leaves->bonded_interactions);
471#endif
472 return {};
473 }
474 return {};
475}
476
477/**
478 * @brief Serialize particles.
479 * Particles need to be serialized here to reduce overhead,
480 * and also to guarantee particles get instantiated after the cell structure
481 * was instantiated (since they store a weak pointer to it).
482 */
483std::string System::get_internal_state() const {
484 auto const p_ids = get_particle_ids();
485 std::vector<std::string> object_states(p_ids.size());
486
487 std::ranges::transform(p_ids, object_states.begin(), [this](auto const p_id) {
488 auto p_obj = context()->make_shared(
489 "Particles::ParticleHandle",
490 {{"id", p_id}, {"__cell_structure", m_leaves->cell_system}});
491 auto &p_handle = dynamic_cast<Particles::ParticleHandle &>(*p_obj);
492 auto const packed_state = p_handle.serialize();
493 // custom particle serialization
494 auto state = Utils::unpack<ObjectState>(packed_state);
495 state.name = "Particles::ParticleHandle";
496 auto const bonds_view = p_handle.call_method("get_bonds_view", {});
497 state.params.emplace_back(std::string{"bonds"}, pack(bonds_view));
498#ifdef EXCLUSIONS
499 auto const exclusions = p_handle.call_method("get_exclusions", {});
500 state.params.emplace_back(std::string{"exclusions"}, pack(exclusions));
501#endif // EXCLUSIONS
502 state.params.emplace_back(std::string{"__cpt_sentinel"}, pack(None{}));
503 return Utils::pack(state);
504 });
505
507}
508
509void System::set_internal_state(std::string const &state) {
510 auto const object_states = Utils::unpack<std::vector<std::string>>(state);
511#ifdef EXCLUSIONS
512 std::unordered_map<int, Variant> exclusions = {};
513#endif // EXCLUSIONS
514 std::unordered_map<int, Variant> bonds = {};
515
516 for (auto const &packed_object : object_states) {
517 auto state = Utils::unpack<ObjectState>(packed_object);
518 VariantMap params = {};
519 for (auto const &kv : state.params) {
520 params[kv.first] = unpack(kv.second, {});
521 }
522 auto const p_id = get_value<int>(params.at("id"));
523 bonds[p_id] = params.extract("bonds").mapped();
524#ifdef EXCLUSIONS
525 exclusions[p_id] = params.extract("exclusions").mapped();
526#endif // EXCLUSIONS
527 params["__cell_structure"] = get_parameter("cell_system");
528 context()->make_shared("Particles::ParticleHandle", params);
529 }
530
531 for (auto const p_id : get_particle_ids()) {
532 auto p_obj = context()->make_shared(
533 "Particles::ParticleHandle",
534 {{"id", p_id}, {"__cell_structure", m_leaves->cell_system}});
535 auto &p_handle = dynamic_cast<Particles::ParticleHandle &>(*p_obj);
536 set_bonds(p_handle, bonds[p_id]);
537#ifdef EXCLUSIONS
538 set_exclusions(p_handle, exclusions[p_id]);
539#endif // EXCLUSIONS
540 }
541}
542
543} // namespace System
544} // namespace ScriptInterface
@ NSQUARE
Atom decomposition (N-square).
@ REGULAR
Regular decomposition.
static int coord(std::string const &s)
Vector implementation and trait types for boost qvm interoperability.
This file contains everything related to the global cell structure / cell system.
void add_parameters(std::vector< AutoParameter > &&params)
Type to indicate no value in Variant.
std::string serialize() const
Variant call_method(const std::string &name, const VariantMap &params)
Call a method on the object.
static std::shared_ptr< System > create()
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Definition Vector.hpp:110
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
This file contains the asynchronous MPI communication.
@ DATA_PART_PROPERTIES
Particle::p.
@ DATA_PART_POSITION
Particle::r.
static void set_bonds(Particles::ParticleHandle &p, Variant const &bonds)
static void set_exclusions(Particles::ParticleHandle &p, Variant const &exclusions)
static void rotate_system(CellStructure &cell_structure, double phi, double theta, double alpha)
static void rescale_particles(CellStructure &cell_structure, int dir, double scale)
Rescale all particle positions in direction dir by a factor scale.
PackedVariant pack(const Variant &v)
Transform a Variant to a PackedVariant.
T get_value(Variant const &v)
Extract value of specific type T from a Variant.
std::unordered_map< std::string, Variant > VariantMap
Definition Variant.hpp:69
auto make_unordered_map_of_variants(std::unordered_map< K, V > const &v)
Definition Variant.hpp:80
boost::make_recursive_variant< None, bool, int, std::size_t, double, std::string, ObjectRef, Utils::Vector3b, Utils::Vector3i, Utils::Vector2d, Utils::Vector3d, Utils::Vector4d, std::vector< int >, std::vector< double >, std::vector< boost::recursive_variant_ >, std::unordered_map< int, boost::recursive_variant_ >, std::unordered_map< std::string, boost::recursive_variant_ > >::type Variant
Possible types for parameters.
Definition Variant.hpp:67
Variant unpack(const PackedVariant &v, std::unordered_map< ObjectId, ObjectRef > const &objects)
Unpack a PackedVariant.
System & get_system()
void set_system(std::shared_ptr< System > new_instance)
Vector3d vec_rotate(const Vector3d &axis, double angle, const Vector3d &vector)
Rotate a vector around an axis.
std::string pack(T const &v)
Pack a serialize type into a string.
Definition pack.hpp:38
Various procedures concerning interactions between particles.
constexpr double INACTIVE_CUTOFF
Cutoff for deactivated interactions.
Routines to calculate the OIF global forces for a particle triple (triangle from mesh).
void init_type_map(int type)
std::vector< int > get_particle_ids()
Get all particle ids.
Particles creation and deletion.
std::unordered_map< std::string, int > propagation_flags_map()
Convert PropagationMode::PropagationMode to name/value pairs.
This file contains all subroutines required to process rotational motion.
void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame, const double phi)
Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi.
Definition rotation.hpp:133
static SteepestDescentParameters params
Currently active steepest descent instance.
Describes a cell structure / cell system.
void update_ghosts_and_resort_particle(unsigned data_parts)
Update ghost particles, with particle resort if needed.
void set_resort_particles(Cells::Resort level)
Increase the local resort level at least to level.
ParticleRange local_particles() const
Description and getter/setter for a parameter.
Container for leaves of the system class.
std::shared_ptr< Thermostat::Thermostat > thermostat
std::shared_ptr< Integrators::IntegratorHandle > integrator
std::shared_ptr< Dipoles::Container > magnetostatics
std::shared_ptr< CellSystem::CellSystem > cell_system
std::shared_ptr< Particles::ParticleList > part
std::shared_ptr< BondBreakage::BreakageSpecs > bond_breakage
std::shared_ptr< Constraints::Constraints > constraints
std::shared_ptr< Interactions::BondedInteractions > bonded_interactions
std::shared_ptr< CollisionDetection::CollisionDetection > collision_detection
std::shared_ptr< Interactions::NonBondedInteractions > non_bonded_inter
std::shared_ptr< Coulomb::Container > electrostatics
std::shared_ptr< Accumulators::AutoUpdateAccumulators > auto_update_accumulators
std::shared_ptr< LeesEdwards::LeesEdwards > lees_edwards