From ccd5d0b266e6200c3c35641b56527f30f9143051 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 19 Jul 2024 19:28:27 +0200 Subject: [PATCH] checkpointing --- src/core/forces.cpp | 14 +--- src/core/immersed_boundary/ibm_tribend.cpp | 4 + src/core/immersed_boundary/ibm_tribend.hpp | 4 +- src/core/immersed_boundary/ibm_triel.cpp | 4 + src/core/immersed_boundary/ibm_triel.hpp | 3 +- .../object-in-fluid/oif_global_forces.cpp | 84 +++++++++++-------- .../object-in-fluid/oif_global_forces.hpp | 39 ++++----- .../oif_global_forces_params.hpp | 16 +--- src/core/object-in-fluid/oif_local_forces.hpp | 21 +---- src/core/system/System.cpp | 3 + src/core/system/System.hpp | 2 + src/core/system/System.impl.hpp | 1 + src/python/espressomd/interactions.py | 37 ++------ src/python/espressomd/system.py | 2 +- src/script_interface/ObjectMap.hpp | 6 +- .../bond_breakage/BreakageSpecs.hpp | 3 +- .../interactions/BondedInteraction.hpp | 62 +++++++++----- .../interactions/BondedInteractions.hpp | 2 - .../interactions/NonBondedInteractions.hpp | 9 +- .../particle_data/ParticleList.cpp | 19 +---- src/script_interface/system/System.cpp | 34 +++++--- src/script_interface/tests/ObjectMap_test.cpp | 1 - testsuite/python/CMakeLists.txt | 1 + testsuite/python/ibm.py | 75 ++++++++++++++++- .../interactions_non-bonded_interface.py | 4 + testsuite/python/save_checkpoint.py | 1 + testsuite/python/system.py | 39 +++++++++ testsuite/python/test_checkpoint.py | 8 +- 28 files changed, 287 insertions(+), 211 deletions(-) create mode 100644 testsuite/python/system.py diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 60f0c769e43..163a94fde8a 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -211,19 +211,7 @@ void System::System::calculate_forces() { collision_detection_cutoff}); constraints->add_forces(particles, get_sim_time()); - - for (int i = 0; i < max_oif_objects; i++) { - // There are two global quantities that need to be evaluated: - // object's surface and object's volume. - auto const area_volume = boost::mpi::all_reduce( - comm_cart, calc_oif_global(i, *box_geo, *cell_structure), std::plus()); - auto const oif_part_area = std::abs(area_volume[0]); - auto const oif_part_vol = std::abs(area_volume[1]); - if (oif_part_area < 1e-100 and oif_part_vol < 1e-100) { - break; - } - add_oif_global_forces(area_volume, i, *box_geo, *cell_structure); - } + oif_global->calculate_forces(); // Must be done here. Forces need to be ghost-communicated immersed_boundaries->volume_conservation(*cell_structure); diff --git a/src/core/immersed_boundary/ibm_tribend.cpp b/src/core/immersed_boundary/ibm_tribend.cpp index 0d11326b78b..9815d2df963 100644 --- a/src/core/immersed_boundary/ibm_tribend.cpp +++ b/src/core/immersed_boundary/ibm_tribend.cpp @@ -95,6 +95,9 @@ IBMTribend::calc_forces(BoxGeometry const &box_geo, Particle const &p1, void IBMTribend::initialize(BoxGeometry const &box_geo, CellStructure const &cell_structure) { + if (is_initialized) { + return; + } // Compute theta0 if (flat) { theta0 = 0.; @@ -128,4 +131,5 @@ void IBMTribend::initialize(BoxGeometry const &box_geo, if (desc < 0.) theta0 = 2. * std::numbers::pi - theta0; } + is_initialized = true; } diff --git a/src/core/immersed_boundary/ibm_tribend.hpp b/src/core/immersed_boundary/ibm_tribend.hpp index 076bdc0addc..7a77b72712c 100644 --- a/src/core/immersed_boundary/ibm_tribend.hpp +++ b/src/core/immersed_boundary/ibm_tribend.hpp @@ -48,6 +48,7 @@ struct IBMTribend { std::tuple p_ids; bool flat; + bool is_initialized; double cutoff() const { return 0.; } @@ -61,7 +62,8 @@ struct IBMTribend { CellStructure const &cell_structure); IBMTribend(int ind1, int ind2, int ind3, int ind4, double kb, bool flat) - : kb{kb}, theta0{0.}, p_ids{ind1, ind2, ind3, ind4}, flat{flat} {} + : kb{kb}, theta0{0.}, p_ids{ind1, ind2, ind3, ind4}, flat{flat}, + is_initialized{false} {} /** Calculate the forces * The equations can be found in Appendix C of @cite kruger12a. diff --git a/src/core/immersed_boundary/ibm_triel.cpp b/src/core/immersed_boundary/ibm_triel.cpp index 9a357190cc4..df19cc1a7a4 100644 --- a/src/core/immersed_boundary/ibm_triel.cpp +++ b/src/core/immersed_boundary/ibm_triel.cpp @@ -193,6 +193,9 @@ IBMTriel::calc_forces(Utils::Vector3d const &vec1, void IBMTriel::initialize(BoxGeometry const &box_geo, CellStructure const &cell_structure) { + if (is_initialized) { + return; + } // collect particles from nodes auto const [ind1, ind2, ind3] = p_ids; auto const pos1 = get_ibm_particle_position(cell_structure, ind1); @@ -221,4 +224,5 @@ void IBMTriel::initialize(BoxGeometry const &box_geo, b1 = (l0 * cosPhi0 - lp0) / area2; b2 = -(l0 * cosPhi0) / area2; area0 = 0.5 * area2; + is_initialized = true; } diff --git a/src/core/immersed_boundary/ibm_triel.hpp b/src/core/immersed_boundary/ibm_triel.hpp index 9b0f8a3396d..98450544f29 100644 --- a/src/core/immersed_boundary/ibm_triel.hpp +++ b/src/core/immersed_boundary/ibm_triel.hpp @@ -56,6 +56,7 @@ struct IBMTriel { /** Particle ids */ std::tuple p_ids; + bool is_initialized; double cutoff() const { return maxDist; } @@ -71,7 +72,7 @@ struct IBMTriel { double k1, double k2) : l0{0.}, lp0{0.}, sinPhi0{0.}, cosPhi0{0.}, area0{0.}, a1{0.}, a2{0.}, b1{0.}, b2{0.}, maxDist{maxDist}, elasticLaw{elasticLaw}, k1{k1}, - k2{k2}, p_ids{ind1, ind2, ind3} {} + k2{k2}, p_ids{ind1, ind2, ind3}, is_initialized{false} {} /** Calculate the forces. * The equations can be found in Appendix C of @cite kruger12a. diff --git a/src/core/object-in-fluid/oif_global_forces.cpp b/src/core/object-in-fluid/oif_global_forces.cpp index 13a8613268c..6ac1465b622 100644 --- a/src/core/object-in-fluid/oif_global_forces.cpp +++ b/src/core/object-in-fluid/oif_global_forces.cpp @@ -21,68 +21,68 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" +#include "bonded_interactions/bonded_interaction_data.hpp" #include "cell_system/CellStructure.hpp" +#include "communication.hpp" +#include "oif_global_forces_params.hpp" #include "system/System.hpp" -#include "bonded_interactions/bonded_interaction_data.hpp" - #include #include +#include +#include + +#include +#include +#include #include -int max_oif_objects = 0; +/** Calculate the mesh volume and area. */ +static auto calc_oif_mesh(int molType, BoxGeometry const &box_geo, + CellStructure &cs, + BondedInteractionsMap const &bonded_ias) { -Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo, - CellStructure &cs) { - auto const &bonded_ias = *::System::get_system().bonded_ias; - // first-fold-then-the-same approach - double partArea = 0.0; - // z volume - double VOL_partVol = 0.; + double area = 0.; + double volume = 0.; - cs.bond_loop([&partArea, &VOL_partVol, &box_geo, &bonded_ias, molType]( + cs.bond_loop([&area, &volume, &box_geo, &bonded_ias, molType]( Particle &p1, int bond_id, std::span partners) { if (p1.mol_id() != molType) return false; - if (boost::get(bonded_ias.at(bond_id).get()) != - nullptr) { - // remaining neighbors fetched + if (boost::get(bonded_ias.at(bond_id).get())) { auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box()); auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11); auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11); - // unfolded positions correct auto const VOL_A = Utils::area_triangle(p11, p22, p33); - partArea += VOL_A; + area += VOL_A; auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33); auto const VOL_dn = VOL_norm.norm(); - auto const VOL_hz = 1.0 / 3.0 * (p11[2] + p22[2] + p33[2]); - VOL_partVol -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz; + auto const VOL_hz = (1.0 / 3.0) * (p11[2] + p22[2] + p33[2]); + volume -= VOL_A * VOL_norm[2] / VOL_dn * VOL_hz; } return false; }); - return {{partArea, VOL_partVol}}; + return Utils::Vector2d{{area, volume}}; } -void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, - BoxGeometry const &box_geo, CellStructure &cs) { - auto const &bonded_ias = *::System::get_system().bonded_ias; - // first-fold-then-the-same approach - double area = area_volume[0]; - double VOL_volume = area_volume[1]; +/** Distribute the OIF global forces to all particles in the mesh. */ +static void add_oif_global_forces(double area, double volume, int molType, + BoxGeometry const &box_geo, CellStructure &cs, + BondedInteractionsMap const &bonded_ias) { - cs.bond_loop([&box_geo, &bonded_ias, area, VOL_volume, molType]( + cs.bond_loop([&box_geo, &bonded_ias, area, volume, molType]( Particle &p1, int bond_id, std::span partners) { if (p1.mol_id() != molType) return false; - if (auto const *iaparams = - boost::get(bonded_ias.at(bond_id).get())) { + auto const *bond_ptr = bonded_ias.at(bond_id).get(); + if (auto const *bond = boost::get(bond_ptr)) { auto const p11 = box_geo.unfolded_position(p1.pos(), p1.image_box()); auto const p22 = p11 + box_geo.get_mi_vector(partners[0]->pos(), p11); auto const p33 = p11 + box_geo.get_mi_vector(partners[1]->pos(), p11); @@ -91,14 +91,13 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, // starting code from volume force auto const VOL_norm = Utils::get_n_triangle(p11, p22, p33).normalize(); auto const VOL_A = Utils::area_triangle(p11, p22, p33); - auto const VOL_vv = (VOL_volume - iaparams->V0) / iaparams->V0; + auto const VOL_vv = (volume - bond->V0) / bond->V0; - auto const VOL_force = - (1.0 / 3.0) * iaparams->kv * VOL_vv * VOL_A * VOL_norm; + auto const VOL_force = (1. / 3.) * bond->kv * VOL_vv * VOL_A * VOL_norm; auto const h = (1. / 3.) * (p11 + p22 + p33); - auto const deltaA = (area - iaparams->A0_g) / iaparams->A0_g; + auto const deltaA = (area - bond->A0_g) / bond->A0_g; auto const m1 = h - p11; auto const m2 = h - p22; @@ -108,7 +107,7 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, auto const m2_length = m2.norm(); auto const m3_length = m3.norm(); - auto const fac = iaparams->ka_g * VOL_A * deltaA / + auto const fac = bond->ka_g * VOL_A * deltaA / (m1_length * m1_length + m2_length * m2_length + m3_length * m3_length); @@ -120,3 +119,22 @@ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, return false; }); } + +void OifGlobal::run_force_kernel() { + auto &system = get_system(); + auto &box_geo = *system.box_geo; + auto &bonded_ias = *system.bonded_ias; + auto &cell_structure = *system.cell_structure; + for (int i = 0; i < max_oif_objects; ++i) { + // There are two global quantities that need to be evaluated: + // object's surface and object's volume. + auto const local = calc_oif_mesh(i, box_geo, cell_structure, bonded_ias); + auto const global = boost::mpi::all_reduce(comm_cart, local, std::plus()); + auto const area = std::abs(global[0]); + auto const volume = std::abs(global[1]); + if (area < 1e-100 and volume < 1e-100) { + break; + } + add_oif_global_forces(area, volume, i, box_geo, cell_structure, bonded_ias); + } +} diff --git a/src/core/object-in-fluid/oif_global_forces.hpp b/src/core/object-in-fluid/oif_global_forces.hpp index 908a349e136..e74d045198f 100644 --- a/src/core/object-in-fluid/oif_global_forces.hpp +++ b/src/core/object-in-fluid/oif_global_forces.hpp @@ -16,33 +16,26 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef CORE_OBJECT_IN_FLUID_OIF_GLOBAL_FORCES_HPP -#define CORE_OBJECT_IN_FLUID_OIF_GLOBAL_FORCES_HPP + +#pragma once + /** \file - * Routines to calculate the OIF global forces energy or/and and force + * Routines to calculate the OIF global forces * for a particle triple (triangle from mesh). See @cite dupin07a. */ -#include "BoxGeometry.hpp" -#include "cell_system/CellStructure.hpp" -#include "oif_global_forces_params.hpp" - -#include - -/** Calculate the OIF global force. - * Called in force_calc() from within forces.cpp - * - calculates the global area and global volume for a cell before the forces - * are handled - * - MPI synchronization with all reduce - * - !!! loop over particles from regular_decomposition !!! - */ -Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo, - CellStructure &cs); +#include "system/Leaf.hpp" -/** Distribute the OIF global forces to all particles in the mesh. */ -void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, - BoxGeometry const &box_geo, CellStructure &cs); +class OifGlobal : public System::Leaf { +public: + int max_oif_objects = 0; -extern int max_oif_objects; + void calculate_forces() { + if (max_oif_objects > 0) { + run_force_kernel(); + } + } -#endif +private: + void run_force_kernel(); +}; diff --git a/src/core/object-in-fluid/oif_global_forces_params.hpp b/src/core/object-in-fluid/oif_global_forces_params.hpp index 8906b759a9b..3f2e1d578d8 100644 --- a/src/core/object-in-fluid/oif_global_forces_params.hpp +++ b/src/core/object-in-fluid/oif_global_forces_params.hpp @@ -16,10 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef OIF_GLOBAL_FORCES_PARAMS_HPP -#define OIF_GLOBAL_FORCES_PARAMS_HPP -#include +#pragma once /** Parameters for OIF global forces * @@ -46,16 +44,4 @@ struct OifGlobalForcesBond { this->V0 = V0; this->kv = kv; } - -private: - friend boost::serialization::access; - template - void serialize(Archive &ar, long int /* version */) { - ar & A0_g; - ar & ka_g; - ar & V0; - ar & kv; - } }; - -#endif diff --git a/src/core/object-in-fluid/oif_local_forces.hpp b/src/core/object-in-fluid/oif_local_forces.hpp index 3f741d1d75e..763a00f44bd 100644 --- a/src/core/object-in-fluid/oif_local_forces.hpp +++ b/src/core/object-in-fluid/oif_local_forces.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP -#define CORE_OBJECT_IN_FLUID_OIF_LOCAL_FORCES_HPP + +#pragma once /** \file * Routines to calculate the OIF local forces for a particle quadruple @@ -80,21 +80,6 @@ struct OifLocalForcesBond { std::tuple calc_forces(BoxGeometry const &box_geo, Particle const &p2, Particle const &p1, Particle const &p3, Particle const &p4) const; - -private: - friend boost::serialization::access; - template - void serialize(Archive &ar, long int /* version */) { - ar & r0; - ar & ks; - ar & kslin; - ar & phi0; - ar & kb; - ar & A01; - ar & A02; - ar & kal; - ar & kvisc; - } }; /** Compute the OIF local forces. @@ -252,5 +237,3 @@ OifLocalForcesBond::calc_forces(BoxGeometry const &box_geo, Particle const &p2, } return std::make_tuple(force2, force1, force3, force4); } - -#endif diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index 3e0d31fa2c8..e9083ff749b 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -72,6 +72,7 @@ System::System(Private) { nonbonded_ias = std::make_shared(); comfixed = std::make_shared(); galilei = std::make_shared(); + oif_global = std::make_shared(); immersed_boundaries = std::make_shared(); #ifdef COLLISION_DETECTION collision_detection = std::make_shared(); @@ -96,6 +97,8 @@ void System::initialize() { bonded_ias->bind_system(handle); thermostat->bind_system(handle); nonbonded_ias->bind_system(handle); + oif_global->bind_system(handle); + immersed_boundaries->bind_system(handle); #ifdef COLLISION_DETECTION collision_detection->bind_system(handle); #endif diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index d30892551ac..bc576045c55 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -48,6 +48,7 @@ class Thermostat; class ComFixed; class Galilei; class Observable_stat; +class OifGlobal; class ImmersedBoundaries; class CollisionDetection; namespace BondBreakage { @@ -280,6 +281,7 @@ class System : public std::enable_shared_from_this { std::shared_ptr thermostat; std::shared_ptr comfixed; std::shared_ptr galilei; + std::shared_ptr oif_global; std::shared_ptr immersed_boundaries; #ifdef COLLISION_DETECTION std::shared_ptr collision_detection; diff --git a/src/core/system/System.impl.hpp b/src/core/system/System.impl.hpp index f6fca811a99..24a3f8c9c0e 100644 --- a/src/core/system/System.impl.hpp +++ b/src/core/system/System.impl.hpp @@ -37,6 +37,7 @@ #include "integrators/Propagation.hpp" #include "lees_edwards/lees_edwards.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "object-in-fluid/oif_global_forces.hpp" #include "thermostat.hpp" #include "BoxGeometry.hpp" diff --git a/src/python/espressomd/interactions.py b/src/python/espressomd/interactions.py index 44d6ce9cdfc..9b95b78fefa 100644 --- a/src/python/espressomd/interactions.py +++ b/src/python/espressomd/interactions.py @@ -739,39 +739,16 @@ def __init__(self, **kwargs): feature = self.__class__.__dict__.get("_so_feature") if feature is not None: code_features.assert_features(feature) - - if "sip" not in kwargs: - assert "bond_id" not in kwargs - # create a new script interface object from bond parameters - # (normal bond creation and checkpointing constructor #2) + if "sip" in kwargs: + super().__init__(**kwargs) + self._ctor_params = self.get_params() + self._bond_id = -1 + else: params = self.get_default_params() params.update(kwargs) super().__init__(**params) self._ctor_params = params self._bond_id = -1 - else: - # create a new bond based on a bond in the script interface - # (checkpointing constructor #3) - super().__init__(**kwargs) - self._bond_id = -1 - self._ctor_params = self.get_params() - - def __reduce__(self): - if self._bond_id != -1: - # checkpointing constructor #1 - return (BondedInteraction._restore_object, - (self.__class__, {"bond_id": self._bond_id})) - else: - # checkpointing constructor #2 - return (BondedInteraction._restore_object, - (self.__class__, self._serialize())) - - def _serialize(self): - return self._ctor_params.copy() - - @classmethod - def _restore_object(cls, derived_class, kwargs): - return derived_class(**kwargs) def __setattr__(self, attr, value): super().__setattr__(attr, value) @@ -1224,7 +1201,7 @@ def get_default_params(self): """Gets default values of optional parameters. """ - return {"k2": 0} + return {"k2": 0., "is_initialized": False, "_cache": None} @script_interface_register @@ -1254,7 +1231,7 @@ def get_default_params(self): """Gets default values of optional parameters. """ - return {"refShape": "Flat"} + return {"refShape": "Flat", "theta0": 0., "is_initialized": False} @script_interface_register diff --git a/src/python/espressomd/system.py b/src/python/espressomd/system.py index 8bcb76eb4ca..ed21a3b760f 100644 --- a/src/python/espressomd/system.py +++ b/src/python/espressomd/system.py @@ -136,7 +136,7 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self._setup_atexit() return - super().__init__(_normal_construction=True, **kwargs) + super().__init__(_regular_constructor=True, **kwargs) if has_features("CUDA"): self.cuda_init_handle = cuda_init.CudaInitHandle() if has_features("WALBERLA"): diff --git a/src/script_interface/ObjectMap.hpp b/src/script_interface/ObjectMap.hpp index 1dfa3aade8b..0b4adea93a6 100644 --- a/src/script_interface/ObjectMap.hpp +++ b/src/script_interface/ObjectMap.hpp @@ -53,7 +53,6 @@ class ObjectMap : public ObjectContainer { virtual void insert_in_core(KeyType const &key, std::shared_ptr const &obj_ptr) = 0; virtual void erase_in_core(KeyType const &key) = 0; - virtual void before_do_construct() = 0; public: ObjectMap() { @@ -63,10 +62,7 @@ class ObjectMap : public ObjectContainer { }); } - void do_construct(VariantMap const ¶ms) override { - before_do_construct(); - restore_from_checkpoint(params); - } + void do_construct(VariantMap const ¶ms) override = 0; /** * @brief Add an element to the map. diff --git a/src/script_interface/bond_breakage/BreakageSpecs.hpp b/src/script_interface/bond_breakage/BreakageSpecs.hpp index f05bcf61fe4..35574840737 100644 --- a/src/script_interface/bond_breakage/BreakageSpecs.hpp +++ b/src/script_interface/bond_breakage/BreakageSpecs.hpp @@ -47,8 +47,9 @@ class BreakageSpecs private: std::shared_ptr<::BondBreakage::BondBreakage> m_bond_breakage; - void before_do_construct() override { + void do_construct(VariantMap const ¶ms) override { m_bond_breakage = std::make_shared<::BondBreakage::BondBreakage>(); + restore_from_checkpoint(params); } void on_bind_system(::System::System &system) override { diff --git a/src/script_interface/interactions/BondedInteraction.hpp b/src/script_interface/interactions/BondedInteraction.hpp index 50d6dd8cead..6350152cf9e 100644 --- a/src/script_interface/interactions/BondedInteraction.hpp +++ b/src/script_interface/interactions/BondedInteraction.hpp @@ -47,6 +47,7 @@ #include #include #include +#include #include namespace ScriptInterface { @@ -79,12 +80,12 @@ class BondedInteraction : public AutoParameters { void check_valid_parameters(VariantMap const ¶ms) const { auto const valid_keys = get_valid_parameters(); for (auto const &key : valid_keys) { - if (params.count(std::string(key)) == 0) { + if (not params.contains(std::string(key))) { throw std::runtime_error("Parameter '" + key + "' is missing"); } } for (auto const &kv : params) { - if (valid_keys.count(kv.first) == 0) { + if (not valid_keys.contains(kv.first)) { throw std::runtime_error("Parameter '" + kv.first + "' is not recognized"); } @@ -450,6 +451,14 @@ class IBMTriel : public BondedInteractionImpl<::IBMTriel> { } return std::string("Skalak"); }}, + {"is_initialized", AutoParameter::read_only, + [this]() { return get_struct().is_initialized; }}, + {"_cache", AutoParameter::read_only, + [this]() { + auto &s = get_struct(); + return std::vector{{s.l0, s.lp0, s.sinPhi0, s.cosPhi0, + s.area0, s.a1, s.a2, s.b1, s.b2}}; + }}, }); } @@ -465,12 +474,26 @@ class IBMTriel : public BondedInteractionImpl<::IBMTriel> { throw std::invalid_argument( "Invalid value for parameter 'elasticLaw': '" + law_name + "'"); } - m_bonded_ia = - std::make_shared<::Bonded_IA_Parameters>(CoreBondedInteraction( - get_value(params, "ind1"), get_value(params, "ind2"), - get_value(params, "ind3"), - get_value(params, "maxDist"), elastic_law, - get_value(params, "k1"), get_value(params, "k2"))); + auto bond = CoreBondedInteraction( + get_value(params, "ind1"), get_value(params, "ind2"), + get_value(params, "ind3"), get_value(params, "maxDist"), + elastic_law, get_value(params, "k1"), + get_value(params, "k2")); + if (get_value_or(params, "is_initialized", false)) { + auto const cache = get_value>(params, "_cache"); + assert(cache.size() == 9ul); + bond.l0 = cache[0]; + bond.lp0 = cache[1]; + bond.sinPhi0 = cache[2]; + bond.cosPhi0 = cache[3]; + bond.area0 = cache[4]; + bond.a1 = cache[5]; + bond.a2 = cache[6]; + bond.b1 = cache[7]; + bond.b2 = cache[8]; + bond.is_initialized = true; + } + m_bonded_ia = std::make_shared<::Bonded_IA_Parameters>(std::move(bond)); } }; @@ -520,6 +543,8 @@ class IBMTribend : public BondedInteractionImpl<::IBMTribend> { }}, {"theta0", AutoParameter::read_only, [this]() { return get_struct().theta0; }}, + {"is_initialized", AutoParameter::read_only, + [this]() { return get_struct().is_initialized; }}, }); } @@ -535,18 +560,15 @@ class IBMTribend : public BondedInteractionImpl<::IBMTribend> { throw std::invalid_argument("Invalid value for parameter 'refShape': '" + shape_name + "'"); } - m_bonded_ia = - std::make_shared<::Bonded_IA_Parameters>(CoreBondedInteraction( - get_value(params, "ind1"), get_value(params, "ind2"), - get_value(params, "ind3"), get_value(params, "ind4"), - get_value(params, "kb"), flat)); - } - - std::set get_valid_parameters() const override { - auto names = - BondedInteractionImpl::get_valid_parameters(); - names.erase("theta0"); - return names; + auto bond = CoreBondedInteraction( + get_value(params, "ind1"), get_value(params, "ind2"), + get_value(params, "ind3"), get_value(params, "ind4"), + get_value(params, "kb"), flat); + if (get_value_or(params, "is_initialized", false)) { + bond.theta0 = get_value(params, "theta0"); + bond.is_initialized = true; + } + m_bonded_ia = std::make_shared<::Bonded_IA_Parameters>(std::move(bond)); } }; diff --git a/src/script_interface/interactions/BondedInteractions.hpp b/src/script_interface/interactions/BondedInteractions.hpp index 319af2701d3..a24f11ba702 100644 --- a/src/script_interface/interactions/BondedInteractions.hpp +++ b/src/script_interface/interactions/BondedInteractions.hpp @@ -57,8 +57,6 @@ class BondedInteractions : public BondedInteractionsBase_t { std::shared_ptr<::BondedInteractionsMap> m_handle; std::unique_ptr m_params; - void before_do_construct() override {} - void do_construct(VariantMap const ¶ms) override { m_handle = std::make_shared<::BondedInteractionsMap>(); m_params = std::make_unique(params); diff --git a/src/script_interface/interactions/NonBondedInteractions.hpp b/src/script_interface/interactions/NonBondedInteractions.hpp index 18da8f818b8..0b0027c6ce9 100644 --- a/src/script_interface/interactions/NonBondedInteractions.hpp +++ b/src/script_interface/interactions/NonBondedInteractions.hpp @@ -67,16 +67,12 @@ class NonBondedInteractions : public System::Leaf { m_handle->bind_system(m_system.lock()); m_handle->on_non_bonded_ia_change(); *m_notify_cutoff_change = [this]() { - if (m_handle) { + if (m_handle and not m_system.expired()) { m_handle->on_non_bonded_ia_change(); } }; } - void on_detach_system(::System::System &) override { - *m_notify_cutoff_change = []() {}; - } - std::pair get_key(Variant const &key) const { try { auto const types = get_value>(key); @@ -96,9 +92,6 @@ class NonBondedInteractions : public System::Leaf { protected: Variant do_call_method(std::string const &method, VariantMap const ¶ms) override { - if (method == "get_n_types") { - return Variant{m_handle->get_max_seen_particle_type() + 1}; - } if (method == "reset") { if (not context()->is_head_node()) { return {}; diff --git a/src/script_interface/particle_data/ParticleList.cpp b/src/script_interface/particle_data/ParticleList.cpp index be47d333a43..9bb11034303 100644 --- a/src/script_interface/particle_data/ParticleList.cpp +++ b/src/script_interface/particle_data/ParticleList.cpp @@ -49,23 +49,6 @@ namespace ScriptInterface { namespace Particles { -#ifdef EXCLUSIONS -static void set_exclusions(ParticleHandle &p, Variant const &exclusions) { - p.call_method("set_exclusions", {{"p_ids", exclusions}}); -} -#endif // EXCLUSIONS - -static void set_bonds(ParticleHandle &p, Variant const &bonds) { - auto const bond_list_flat = get_value>>(bonds); - for (auto const &bond_flat : bond_list_flat) { - auto const bond_id = bond_flat[0]; - auto const part_id = - std::vector{bond_flat.begin() + 1, bond_flat.end()}; - p.call_method("add_bond", - {{"bond_id", bond_id}, {"part_id", std::move(part_id)}}); - } -} - #ifdef EXCLUSIONS /** * @brief Use the bond topology to automatically add exclusions between @@ -208,7 +191,7 @@ Variant ParticleList::do_call_method(std::string const &name, context()->make_shared("Particles::ParticleHandle", local_params)); #ifdef EXCLUSIONS if (params.count("exclusions")) { - set_exclusions(*so, params.at("exclusions")); + so->call_method("set_exclusions", {{"p_ids", params.at("exclusions")}}); } #endif // EXCLUSIONS return so->get_parameter("id"); diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index 39ef51241e3..214d91ee642 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -181,7 +181,11 @@ System::System() : m_instance{}, m_leaves{std::make_shared()} { }); }, [this]() { return m_instance->get_min_global_cut(); }}, - {"max_oif_objects", ::max_oif_objects}, + {"max_oif_objects", + [this](Variant const &v) { + m_instance->oif_global->max_oif_objects = get_value(v); + }, + [this]() { return m_instance->oif_global->max_oif_objects; }}, }); // note: the order of leaves matters! e.g. bonds depend on thermostats, // and thus a thermostat object must be instantiated before the bonds @@ -226,13 +230,15 @@ void System::do_construct(VariantMap const ¶ms) { * runtime errors about the local geometry being smaller * than the interaction range would be raised. */ - if (not params.contains("box_l")) { - throw std::domain_error("Required argument 'box_l' not provided."); - } - if (params.contains("_normal_construction") and system_created) { - throw std::runtime_error( - "You can only have one instance of the system class at a time"); - } + context()->parallel_try_catch([&]() { + if (not params.contains("box_l")) { + throw std::domain_error("Required argument 'box_l' not provided."); + } + if (params.contains("_regular_constructor") and system_created) { + throw std::runtime_error( + "You can only have one instance of the system class at a time"); + } + }); m_instance = ::System::System::create(); ::System::set_system(m_instance); @@ -244,17 +250,19 @@ void System::do_construct(VariantMap const ¶ms) { m_instance->lb.bind_system(m_instance); m_instance->ek.bind_system(m_instance); - if (params.contains("_normal_construction")) { + if (params.contains("_regular_constructor")) { std::set const setable_properties = { "box_l", "min_global_cut", "periodicity", "time", "time_step", "force_cap", - "max_oif_objects", "_normal_construction"}; + "max_oif_objects", "_regular_constructor"}; for (auto const &kv : params) { if (not setable_properties.contains(kv.first)) { - throw std::domain_error( - "Property '" + kv.first + - "' can not be set via argument to System class"); + context()->parallel_try_catch([&kv]() { + throw std::domain_error( + "Property '" + kv.first + + "' cannot be set via argument to System class"); + }); } } for (std::string attr : diff --git a/src/script_interface/tests/ObjectMap_test.cpp b/src/script_interface/tests/ObjectMap_test.cpp index 8c54ec4f47a..e735dd68c29 100644 --- a/src/script_interface/tests/ObjectMap_test.cpp +++ b/src/script_interface/tests/ObjectMap_test.cpp @@ -43,7 +43,6 @@ struct ObjectMapImpl : ObjectMap { std::unordered_map mock_core; private: - void before_do_construct() override {} void insert_in_core(KeyType const &key, ObjectRef const &obj_ptr) override { next_key = std::max(next_key, key + 1); mock_core[key] = obj_ptr; diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 73618834b1c..e09ee314478 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -375,6 +375,7 @@ python_test(FILE elc.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE elc_vs_analytic.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE rotation.py MAX_NUM_PROC 1) python_test(FILE shapes.py MAX_NUM_PROC 1) +python_test(FILE system.py MAX_NUM_PROC 2) python_test(FILE h5md.py MAX_NUM_PROC 2) python_test(FILE h5md.py MAX_NUM_PROC 1 SUFFIX 1_core) python_test(FILE p3m_fft.py MAX_NUM_PROC 6) diff --git a/testsuite/python/ibm.py b/testsuite/python/ibm.py index 36affc2791f..e73696fac29 100644 --- a/testsuite/python/ibm.py +++ b/testsuite/python/ibm.py @@ -18,6 +18,7 @@ # import numpy as np import itertools +import pickle import unittest as ut import espressomd @@ -46,6 +47,10 @@ def compute_dihedral_angle(self, pos0, pos1, pos2, pos3): cos_alpha = min(1, np.dot(n1, n2)) alpha = np.arccos(cos_alpha) + desc = np.dot((pos1 - pos0), np.cross(n1, n2)) + if desc > 0.: + alpha = 2. * np.pi - alpha + return alpha def test_tribend(self): @@ -55,7 +60,7 @@ def test_tribend(self): system = self.system system.thermostat.set_langevin(kT=0, gamma=10, seed=1) - # Add four particles + # Add four particles arranged in a U-shape, dihedral angle is zero p0 = system.part.add(pos=[4, 4, 4]) p1 = system.part.add(pos=[4, 4, 5]) p2 = system.part.add(pos=[4, 5, 5]) @@ -84,8 +89,9 @@ def test_tribend(self): # Perform integration system.integrator.run(200) - angle = self.compute_dihedral_angle(p0.pos, p1.pos, p2.pos, p3.pos) - self.assertLess(angle, 2E-2) + theta = self.compute_dihedral_angle(p0.pos, p1.pos, p2.pos, p3.pos) + theta = np.fmod(theta + 0.1, 2. * np.pi) - 0.1 + self.assertLess(np.abs(theta), 2E-2) # IBM doesn't implement energy and pressure kernels. energy = self.system.analysis.energy() @@ -230,6 +236,69 @@ def test_volcons(self): mesh_center = np.mean(partcls.pos, axis=0) np.testing.assert_allclose(mesh_center, mesh_center_ref, rtol=1e-3) + def test_tribend_checkpointing(self): + system = self.system + # Add four particles arranged in a U-shape, dihedral angle is non-zero + p1 = system.part.add(pos=[3.5, 3.5, 4.0]) + p2 = system.part.add(pos=[4.5, 3.5, 4.0]) + p3 = system.part.add(pos=[4.5, 4.5, 4.0]) + p4 = system.part.add(pos=[3.5, 4.5, 4.2]) + theta0 = self.compute_dihedral_angle(p1.pos, p2.pos, p3.pos, p4.pos) + assert abs(theta0 - 6.) < 0.1 + tribend_original = espressomd.interactions.IBM_Tribend( + ind1=p1.id, ind2=p2.id, ind3=p3.id, ind4=p4.id, kb=1., refShape="Initial") + tribend_unpickled = pickle.loads(pickle.dumps(tribend_original)) + for tribend in [tribend_original, tribend_unpickled]: + self.assertFalse(tribend.is_initialized) + self.assertEqual(tribend.theta0, 0.) + self.assertEqual(tribend.ind1, p1.id) + self.assertEqual(tribend.ind2, p2.id) + self.assertEqual(tribend.ind3, p3.id) + self.assertEqual(tribend.ind4, p4.id) + system.bonded_inter.add(tribend_original) + system.bonded_inter.add(tribend_unpickled) + for tribend in [tribend_original, tribend_unpickled]: + self.assertTrue(tribend.is_initialized) + self.assertAlmostEqual(tribend.theta0, theta0, delta=1e-6) + self.assertEqual(tribend.ind1, p1.id) + self.assertEqual(tribend.ind2, p2.id) + self.assertEqual(tribend.ind3, p3.id) + self.assertEqual(tribend.ind4, p4.id) + + def test_triel_checkpointing(self): + system = self.system + p1 = system.part.add(pos=[1, 4, 4]) + p2 = system.part.add(pos=[1, 4, 5]) + p3 = system.part.add(pos=[1, 5, 5]) + sqrt2 = np.sqrt(2.) + cache = [sqrt2, 1., sqrt2 / 2., sqrt2 / 2., 0.5, -1., 1., 0., -1.] + triel_original = espressomd.interactions.IBM_Triel( + ind1=p1.id, ind2=p2.id, ind3=p3.id, elasticLaw="Skalak", k1=15., + k2=0.5, maxDist=2.4) + triel_unpickled = pickle.loads(pickle.dumps(triel_original)) + for triel in [triel_original, triel_unpickled]: + self.assertFalse(triel.is_initialized) + np.testing.assert_allclose(np.copy(triel._cache), 0.) + self.assertEqual(triel.ind1, p1.id) + self.assertEqual(triel.ind2, p2.id) + self.assertEqual(triel.ind3, p3.id) + self.assertEqual(triel.elasticLaw, "Skalak") + self.assertEqual(triel.k1, 15.) + self.assertEqual(triel.k2, 0.5) + self.assertEqual(triel.maxDist, 2.4) + system.bonded_inter.add(triel_original) + system.bonded_inter.add(triel_unpickled) + for triel in [triel_original, triel_unpickled]: + self.assertTrue(triel.is_initialized) + np.testing.assert_allclose(np.copy(triel._cache), cache) + self.assertEqual(triel.ind1, p1.id) + self.assertEqual(triel.ind2, p2.id) + self.assertEqual(triel.ind3, p3.id) + self.assertEqual(triel.elasticLaw, "Skalak") + self.assertEqual(triel.k1, 15.) + self.assertEqual(triel.k2, 0.5) + self.assertEqual(triel.maxDist, 2.4) + if __name__ == "__main__": ut.main() diff --git a/testsuite/python/interactions_non-bonded_interface.py b/testsuite/python/interactions_non-bonded_interface.py index 755df31ec1d..f6a119b1d8f 100644 --- a/testsuite/python/interactions_non-bonded_interface.py +++ b/testsuite/python/interactions_non-bonded_interface.py @@ -188,6 +188,10 @@ def test_exceptions(self): with self.assertRaisesRegex(ValueError, "LJ parameter 'shift' has to be 'auto' or a float"): self.system.non_bonded_inter[0, 0].lennard_jones.set_params( epsilon=1., sigma=2., cutoff=3., shift="automatic") + with self.assertRaisesRegex(ValueError, r"NonBondedInteractions\[\] expects two particle types as indices"): + self.system.non_bonded_inter[0, 0, 1].lennard_jones + with self.assertRaisesRegex(ValueError, r"NonBondedInteractions\[\] expects two particle types as indices"): + self.system.non_bonded_inter[-1, 0].lennard_jones skin = self.system.cell_system.skin box_l = self.system.box_l diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 417055baad3..c2c72434005 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -307,6 +307,7 @@ ibm_triel_bond = espressomd.interactions.IBM_Triel( ind1=p1.id, ind2=p2.id, ind3=p3.id, k1=1.1, k2=1.2, maxDist=1.6, elasticLaw="NeoHookean") +system.bonded_inter.add(ibm_tribend_bond) break_spec = espressomd.bond_breakage.BreakageSpec( breakage_length=5., action_type="delete_bond") system.bond_breakage[strong_harmonic_bond._bond_id] = break_spec diff --git a/testsuite/python/system.py b/testsuite/python/system.py new file mode 100644 index 00000000000..8cd0ac86255 --- /dev/null +++ b/testsuite/python/system.py @@ -0,0 +1,39 @@ +# +# Copyright (C) 2024 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import unittest as ut +import espressomd + + +class Test(ut.TestCase): + + def test_system(self): + with self.assertRaisesRegex(ValueError, "Required argument 'box_l' not provided"): + espressomd.System() + with self.assertRaisesRegex(ValueError, "Property 'unknown' cannot be set via argument to System class"): + espressomd.System(box_l=[1., 1., 1.], unknown=1) + system = espressomd.System(box_l=[1., 1., 1.], min_global_cut=0.01) + self.assertEqual(system.min_global_cut, 0.01) + self.assertEqual(system.time_step, -1.) + with self.assertRaisesRegex(RuntimeError, "You can only have one instance of the system class at a time"): + espressomd.System(box_l=[1., 1., 1.]) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index cc4e28dc84a..7b49c2cb90a 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -696,12 +696,12 @@ def test_bonded_inter(self): ibm_volcons_bond.params, {'softID': 15, 'kappaV': 0.01}) self.assertEqual( {**ibm_tribend_bond.params, **{'theta0': 0.}}, - {'kb': 2., 'theta0': 0., 'refShape': 'Initial', - 'ind1': 0, 'ind2': 1, 'ind3': 3, 'ind4': 4}) + {'kb': 2., 'refShape': 'Initial', 'is_initialized': True, + 'theta0': 0., 'ind1': 0, 'ind2': 1, 'ind3': 3, 'ind4': 4}) self.assertEqual( - ibm_triel_bond.params, + {k: v for k, v in ibm_triel_bond.params.items() if k[0] != "_"}, {'k1': 1.1, 'k2': 1.2, 'maxDist': 1.6, 'elasticLaw': 'NeoHookean', - 'ind1': 0, 'ind2': 1, 'ind3': 3}) + 'ind1': 0, 'ind2': 1, 'ind3': 3, 'is_initialized': False}) # check new bonds can be added if not has_lb_mode: new_bond = espressomd.interactions.HarmonicBond(r_0=0.2, k=1.)