From fc1c08fe2728809df88b5f8182e30ffab3f71a95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 24 Aug 2023 11:56:05 +0200 Subject: [PATCH 1/6] core: Throw on box resize when constraints are present Fixes a regression introduced in e2294223d589dba3e1b1cef21b01f7b9d0f2f1c5 --- src/core/constraints/Constraints.hpp | 2 +- src/core/event.cpp | 3 +++ testsuite/python/constraint_shape_based.py | 14 ++++++++++++++ testsuite/python/integrator_steepest_descent.py | 1 + 4 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/core/constraints/Constraints.hpp b/src/core/constraints/Constraints.hpp index 754301b86a9..ead490ea029 100644 --- a/src/core/constraints/Constraints.hpp +++ b/src/core/constraints/Constraints.hpp @@ -102,7 +102,7 @@ template class Constraints { } void on_boxl_change() const { - if (not this->empty()) { + if (not m_constraints.empty()) { throw std::runtime_error("The box size can not be changed because there " "are active constraints."); } diff --git a/src/core/event.cpp b/src/core/event.cpp index 14b6957fda4..7cc4f59d809 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -31,6 +31,7 @@ #include "collision.hpp" #include "communication.hpp" #include "config/config.hpp" +#include "constraints.hpp" #include "cuda/init.hpp" #include "cuda/utils.hpp" #include "electrostatics/coulomb.hpp" @@ -241,6 +242,8 @@ void on_boxl_change(bool skip_method_adaption) { system.dipoles.on_boxl_change(); #endif } + + Constraints::constraints.on_boxl_change(); } void on_cell_structure_change() { diff --git a/testsuite/python/constraint_shape_based.py b/testsuite/python/constraint_shape_based.py index 07baeb4c5a2..738ef27606c 100644 --- a/testsuite/python/constraint_shape_based.py +++ b/testsuite/python/constraint_shape_based.py @@ -36,6 +36,9 @@ class ShapeBasedConstraintTest(ut.TestCase): box_l = 30. system = espressomd.System(box_l=3 * [box_l]) + def setUp(self): + self.system.box_l = 3 * [self.box_l] + def tearDown(self): self.system.part.clear() self.system.constraints.clear() @@ -1067,6 +1070,17 @@ def test_torus(self): system.non_bonded_inter[0, 1].lennard_jones.set_params( epsilon=0.0, sigma=0.0, cutoff=0.0, shift=0) + def test_exceptions(self): + system = self.system + wall = espressomd.shapes.Wall(normal=[0., 1., 0.], dist=0.) + constraint = espressomd.constraints.ShapeBasedConstraint( + shape=wall, particle_type=1) + system.constraints.add(constraint) + with self.assertRaisesRegex(RuntimeError, "there are active constraints"): + system.box_l = 0.5 * system.box_l + system.constraints.remove(constraint) + system.box_l = 0.75 * system.box_l + if __name__ == "__main__": ut.main() diff --git a/testsuite/python/integrator_steepest_descent.py b/testsuite/python/integrator_steepest_descent.py index dfb5cf91965..486fbf40899 100644 --- a/testsuite/python/integrator_steepest_descent.py +++ b/testsuite/python/integrator_steepest_descent.py @@ -55,6 +55,7 @@ def setUp(self): def tearDown(self): self.system.part.clear() + self.system.constraints.clear() self.system.integrator.set_vv() def test_relaxation_integrator(self): From 60da45645a0e07057b5ec93e107ea889100501a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 25 Aug 2023 15:13:51 +0200 Subject: [PATCH 2/6] core: Prevent duplicate objects in ObjectList --- src/core/accumulators.cpp | 35 ++++++++++++------- src/core/accumulators.hpp | 1 + src/core/constraints/Constraints.hpp | 14 ++++---- src/script_interface/ObjectList.hpp | 14 ++++++++ .../accumulators/AutoUpdateAccumulators.hpp | 5 +++ .../constraints/Constraints.hpp | 3 ++ src/script_interface/shapes/Union.hpp | 3 ++ .../tests/ObjectList_test.cpp | 8 +++-- src/script_interface/walberla/EKContainer.hpp | 3 ++ src/script_interface/walberla/EKReactions.hpp | 3 ++ src/shapes/include/shapes/Union.hpp | 4 +++ src/shapes/unit_tests/Union_test.cpp | 4 +++ testsuite/python/script_interface.py | 16 +++++++++ 13 files changed, 90 insertions(+), 23 deletions(-) diff --git a/src/core/accumulators.cpp b/src/core/accumulators.cpp index 6f707afd0a3..bc2b8306033 100644 --- a/src/core/accumulators.cpp +++ b/src/core/accumulators.cpp @@ -60,24 +60,33 @@ int auto_update_next_update() { }); } +namespace detail { +struct MatchPredicate { + AccumulatorBase const *m_acc; + template bool operator()(T const &a) const { + return a.acc == m_acc; + } +}; +} // namespace detail + void auto_update_add(AccumulatorBase *acc) { - assert(acc); - assert(std::find_if(auto_update_accumulators.begin(), - auto_update_accumulators.end(), [acc](auto const &a) { - return a.acc == acc; - }) == auto_update_accumulators.end()); + assert(not auto_update_contains(acc)); auto_update_accumulators.emplace_back(acc); } + void auto_update_remove(AccumulatorBase *acc) { - assert(std::find_if(auto_update_accumulators.begin(), - auto_update_accumulators.end(), [acc](auto const &a) { - return a.acc == acc; - }) != auto_update_accumulators.end()); + assert(auto_update_contains(acc)); + auto const beg = auto_update_accumulators.begin(); + auto const end = auto_update_accumulators.end(); auto_update_accumulators.erase( - boost::remove_if( - auto_update_accumulators, - [acc](AutoUpdateAccumulator const &au) { return au.acc == acc; }), - auto_update_accumulators.end()); + std::remove_if(beg, end, detail::MatchPredicate{acc}), end); +} + +bool auto_update_contains(AccumulatorBase const *acc) noexcept { + assert(acc); + auto const beg = auto_update_accumulators.begin(); + auto const end = auto_update_accumulators.end(); + return std::find_if(beg, end, detail::MatchPredicate{acc}) != end; } } // namespace Accumulators diff --git a/src/core/accumulators.hpp b/src/core/accumulators.hpp index 68a9da09a05..bc8199b05d6 100644 --- a/src/core/accumulators.hpp +++ b/src/core/accumulators.hpp @@ -31,6 +31,7 @@ namespace Accumulators { */ void auto_update(boost::mpi::communicator const &comm, int steps); int auto_update_next_update(); +bool auto_update_contains(AccumulatorBase const *) noexcept; void auto_update_add(AccumulatorBase *); void auto_update_remove(AccumulatorBase *); diff --git a/src/core/constraints/Constraints.hpp b/src/core/constraints/Constraints.hpp index ead490ea029..8c19dbf6678 100644 --- a/src/core/constraints/Constraints.hpp +++ b/src/core/constraints/Constraints.hpp @@ -49,22 +49,20 @@ template class Constraints { container_type m_constraints; public: + bool contains(std::shared_ptr const &constraint) const noexcept { + return std::find(begin(), end(), constraint) != end(); + } void add(std::shared_ptr const &constraint) { if (not constraint->fits_in_box(box_geo.length())) { throw std::runtime_error("Constraint not compatible with box size."); } - assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) == - m_constraints.end()); - + assert(not contains(constraint)); m_constraints.emplace_back(constraint); on_constraint_change(); } void remove(std::shared_ptr const &constraint) { - assert(std::find(m_constraints.begin(), m_constraints.end(), constraint) != - m_constraints.end()); - m_constraints.erase( - std::remove(m_constraints.begin(), m_constraints.end(), constraint), - m_constraints.end()); + assert(contains(constraint)); + m_constraints.erase(std::remove(begin(), end(), constraint), end()); on_constraint_change(); } diff --git a/src/script_interface/ObjectList.hpp b/src/script_interface/ObjectList.hpp index 72b7bb8f9b7..f9ad6f61093 100644 --- a/src/script_interface/ObjectList.hpp +++ b/src/script_interface/ObjectList.hpp @@ -52,6 +52,8 @@ class ObjectList : public ObjectContainer { virtual void add_in_core(const std::shared_ptr &obj_ptr) = 0; virtual void remove_in_core(const std::shared_ptr &obj_ptr) = 0; + virtual bool + has_in_core(const std::shared_ptr &obj_ptr) const = 0; public: ObjectList() { @@ -74,6 +76,12 @@ class ObjectList : public ObjectContainer { * @param element The element to add. */ void add(std::shared_ptr const &element) { + if (has_in_core(element)) { + if (Base::context()->is_head_node()) { + throw std::runtime_error("This object is already present in the list"); + } + throw Exception(""); + } add_in_core(element); m_elements.push_back(element); } @@ -84,6 +92,12 @@ class ObjectList : public ObjectContainer { * @param element The element to remove. */ void remove(std::shared_ptr const &element) { + if (not has_in_core(element)) { + if (Base::context()->is_head_node()) { + throw std::runtime_error("This object is absent from the list"); + } + throw Exception(""); + } remove_in_core(element); m_elements.erase(std::remove(m_elements.begin(), m_elements.end(), element), m_elements.end()); diff --git a/src/script_interface/accumulators/AutoUpdateAccumulators.hpp b/src/script_interface/accumulators/AutoUpdateAccumulators.hpp index 842f2e899bb..a5ada2cbd98 100644 --- a/src/script_interface/accumulators/AutoUpdateAccumulators.hpp +++ b/src/script_interface/accumulators/AutoUpdateAccumulators.hpp @@ -29,6 +29,11 @@ namespace ScriptInterface { namespace Accumulators { class AutoUpdateAccumulators : public ObjectList { + bool + has_in_core(std::shared_ptr const &obj_ptr) const override { + return ::Accumulators::auto_update_contains(obj_ptr->accumulator().get()); + } + void add_in_core(std::shared_ptr const &obj_ptr) override { ::Accumulators::auto_update_add(obj_ptr->accumulator().get()); } diff --git a/src/script_interface/constraints/Constraints.hpp b/src/script_interface/constraints/Constraints.hpp index f2b9d537c37..a8b304c872a 100644 --- a/src/script_interface/constraints/Constraints.hpp +++ b/src/script_interface/constraints/Constraints.hpp @@ -32,6 +32,9 @@ namespace ScriptInterface { namespace Constraints { class Constraints : public ObjectList { + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return ::Constraints::constraints.contains(obj_ptr->constraint()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { ::Constraints::constraints.add(obj_ptr->constraint()); } diff --git a/src/script_interface/shapes/Union.hpp b/src/script_interface/shapes/Union.hpp index 1227ec0b504..475e899ef24 100644 --- a/src/script_interface/shapes/Union.hpp +++ b/src/script_interface/shapes/Union.hpp @@ -39,6 +39,9 @@ class Union : public ObjectList { Union() : m_core_shape(std::make_shared<::Shapes::Union>()) {} private: + bool has_in_core(const std::shared_ptr &obj_ptr) const override { + return m_core_shape->contains(obj_ptr->shape()); + } void add_in_core(const std::shared_ptr &obj_ptr) override { m_core_shape->add(obj_ptr->shape()); } diff --git a/src/script_interface/tests/ObjectList_test.cpp b/src/script_interface/tests/ObjectList_test.cpp index 2645ffdd87e..4298d37fa1b 100644 --- a/src/script_interface/tests/ObjectList_test.cpp +++ b/src/script_interface/tests/ObjectList_test.cpp @@ -43,6 +43,10 @@ struct ObjectListImpl : ObjectList { std::vector mock_core; private: + bool has_in_core(const ObjectRef &obj_ptr) const override { + return std::find(mock_core.begin(), mock_core.end(), obj_ptr) != + mock_core.end(); + } void add_in_core(const ObjectRef &obj_ptr) override { mock_core.push_back(obj_ptr); } @@ -82,8 +86,8 @@ BOOST_AUTO_TEST_CASE(removing_elements) { BOOST_AUTO_TEST_CASE(clearing_elements) { // A cleared list is empty. ObjectListImpl list; - list.add(ObjectRef{}); - list.add(ObjectRef{}); + list.add(std::make_shared()); + list.add(std::make_shared()); list.clear(); BOOST_CHECK(list.elements().empty()); BOOST_CHECK(list.mock_core.empty()); diff --git a/src/script_interface/walberla/EKContainer.hpp b/src/script_interface/walberla/EKContainer.hpp index 5f33f50c28e..26a6c617ae4 100644 --- a/src/script_interface/walberla/EKContainer.hpp +++ b/src/script_interface/walberla/EKContainer.hpp @@ -61,6 +61,9 @@ class EKContainer : public ObjectList { std::shared_ptr<::EK::EKWalberla::ek_container_type> m_ek_container; bool m_is_active; + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return m_ek_container->contains(obj_ptr->get_ekinstance()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { context()->parallel_try_catch( [this, &obj_ptr]() { m_ek_container->add(obj_ptr->get_ekinstance()); }); diff --git a/src/script_interface/walberla/EKReactions.hpp b/src/script_interface/walberla/EKReactions.hpp index 84104049454..12f0fe81695 100644 --- a/src/script_interface/walberla/EKReactions.hpp +++ b/src/script_interface/walberla/EKReactions.hpp @@ -40,6 +40,9 @@ namespace ScriptInterface::walberla { class EKReactions : public ObjectList { std::shared_ptr<::EK::EKWalberla::ek_reactions_type> m_ek_reactions; + bool has_in_core(std::shared_ptr const &obj_ptr) const override { + return m_ek_reactions->contains(obj_ptr->get_instance()); + } void add_in_core(std::shared_ptr const &obj_ptr) override { m_ek_reactions->add(obj_ptr->get_instance()); } diff --git a/src/shapes/include/shapes/Union.hpp b/src/shapes/include/shapes/Union.hpp index 492a7c42cc4..c719ac37826 100644 --- a/src/shapes/include/shapes/Union.hpp +++ b/src/shapes/include/shapes/Union.hpp @@ -33,6 +33,10 @@ namespace Shapes { class Union : public Shape { public: + bool contains(std::shared_ptr const &shape) const noexcept { + return std::find(m_shapes.begin(), m_shapes.end(), shape) != m_shapes.end(); + } + void add(std::shared_ptr const &shape) { m_shapes.emplace_back(shape); } diff --git a/src/shapes/unit_tests/Union_test.cpp b/src/shapes/unit_tests/Union_test.cpp index eb2357f5386..dfacb5ecfa8 100644 --- a/src/shapes/unit_tests/Union_test.cpp +++ b/src/shapes/unit_tests/Union_test.cpp @@ -42,8 +42,12 @@ BOOST_AUTO_TEST_CASE(dist_function) { wall2->d() = -10.0; Shapes::Union uni; + BOOST_CHECK(not uni.contains(wall1)); + BOOST_CHECK(not uni.contains(wall2)); uni.add(wall1); uni.add(wall2); + BOOST_CHECK(uni.contains(wall1)); + BOOST_CHECK(uni.contains(wall2)); auto check_union = [&wall1, &wall2, &uni](Utils::Vector3d const &pos) { double wall1_dist; diff --git a/testsuite/python/script_interface.py b/testsuite/python/script_interface.py index 7665a682514..8848e45e18c 100644 --- a/testsuite/python/script_interface.py +++ b/testsuite/python/script_interface.py @@ -102,6 +102,22 @@ def test_autoparameter_exceptions(self): with self.assertRaisesRegex(AttributeError, "Object 'HarmonicBond' has no attribute 'unknown'"): bond.unknown + def test_objectlist_exceptions(self): + """Check ObjectList framework""" + wall = espressomd.shapes.Wall(normal=[-1, 0, 0]) + constraint = espressomd.constraints.ShapeBasedConstraint(shape=wall) + constraints = espressomd.constraints.Constraints() + self.assertEqual(len(constraints), 0) + constraints.add(constraint) + self.assertEqual(len(constraints), 1) + with self.assertRaisesRegex(RuntimeError, "This object is already present in the list"): + constraints.add(constraint) + self.assertEqual(len(constraints), 1) + constraints.remove(constraint) + self.assertEqual(len(constraints), 0) + with self.assertRaisesRegex(RuntimeError, "This object is absent from the list"): + constraints.remove(constraint) + def test_feature_exceptions(self): """Check feature verification""" all_features = set(espressomd.code_info.all_features()) From a9e0e9b154e286180ccfd9791a815b98e00e325b Mon Sep 17 00:00:00 2001 From: David Beyer Date: Wed, 9 Aug 2023 10:35:32 +0200 Subject: [PATCH 3/6] Removed overhead from make_reaction_attempt --- src/core/event.cpp | 6 ++++++ src/core/event.hpp | 3 +++ .../reaction_methods/ReactionAlgorithm.cpp | 19 ++++++++++++++----- 3 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/core/event.cpp b/src/core/event.cpp index 7cc4f59d809..b78142796ac 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -160,6 +160,12 @@ void on_particle_charge_change() { partCfg().invalidate(); } +void on_particle_local_change() { + cells_update_ghosts(global_ghost_flags()); + + recalc_forces = true; +} + void on_particle_change() { if (cell_structure.decomposition_type() == CellStructureType::CELL_STRUCTURE_HYBRID) { diff --git a/src/core/event.hpp b/src/core/event.hpp index 3fda162f9ad..bec60711a08 100644 --- a/src/core/event.hpp +++ b/src/core/event.hpp @@ -66,6 +66,9 @@ void on_particle_change(); /** called every time the charge of a particle has changed. */ void on_particle_charge_change(); +/** called every time that local properties of a particle have changed. */ +void on_particle_local_change(); + /** called every time the Coulomb parameters are changed. all Coulomb methods have a short range part, aka near field diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index f19f91b6d40..d2e10732aae 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -149,6 +149,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, auto const random_index = i_random(number_of_particles_with_type(type)); return get_random_p_id(type, random_index); }; + auto only_local_changes = true; for (int i = 0; i < std::min(n_product_types, n_reactant_types); i++) { auto const n_product_coef = reaction.product_coefficients[i]; auto const n_reactant_coef = reaction.reactant_coefficients[i]; @@ -156,6 +157,11 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, // particles of reactant_types(i) to product_types(i) auto const old_type = reaction.reactant_types[i]; auto const new_type = reaction.product_types[i]; +#ifdef ELECTROSTATICS + if (charges_of_types.at(new_type) != charges_of_types.at(old_type)) { + only_local_changes = false; + } +#endif for (int j = 0; j < std::min(n_product_coef, n_reactant_coef); j++) { auto const p_id = get_random_p_id_of_type(old_type); on_particle_type_change(p_id, old_type, new_type); @@ -167,7 +173,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, } bookkeeping.changed.emplace_back(p_id, old_type); } - on_particle_change(); // create product_coefficients(i)-reactant_coefficients(i) many product // particles iff product_coefficients(i)-reactant_coefficients(i)>0, // iff product_coefficients(i)-reactant_coefficients(i)<0, hide this number @@ -180,7 +185,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); bookkeeping.created.emplace_back(p_id); } - on_particle_change(); + only_local_changes = false; } else if (delta_n < 0) { auto const type = reaction.reactant_types[i]; for (int j = 0; j < -delta_n; j++) { @@ -189,7 +194,7 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); hide_particle(p_id, type); } - on_particle_change(); + only_local_changes = false; } } // create or hide particles of types with noncorresponding replacement types @@ -204,7 +209,6 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); hide_particle(p_id, type); } - on_particle_change(); } else { // create additional product_types particles auto const type = reaction.product_types[i]; @@ -213,9 +217,14 @@ void ReactionAlgorithm::make_reaction_attempt(SingleReaction const &reaction, check_exclusion_range(p_id, type); bookkeeping.created.emplace_back(p_id); } - on_particle_change(); } } + // determine which fine-grained event to trigger + if (n_product_types == n_reactant_types and only_local_changes) { + on_particle_local_change(); + } else { + on_particle_change(); + } } std::unordered_map From fcb3947ec6bd4188613d819207636d972320e556 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 31 Aug 2023 19:58:51 +0200 Subject: [PATCH 4/6] Add test case --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/reaction_trivial.py | 78 ++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+) create mode 100644 testsuite/python/reaction_trivial.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index ce792497bc6..a3b20673a03 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -252,6 +252,7 @@ python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1) python_test(FILE script_interface.py MAX_NUM_PROC 4) python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 2) python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4) +python_test(FILE reaction_trivial.py MAX_NUM_PROC 2) python_test(FILE reaction_complex.py MAX_NUM_PROC 1) python_test(FILE reaction_bookkeeping.py MAX_NUM_PROC 1) python_test(FILE widom_insertion.py MAX_NUM_PROC 1) diff --git a/testsuite/python/reaction_trivial.py b/testsuite/python/reaction_trivial.py new file mode 100644 index 00000000000..7229928b11f --- /dev/null +++ b/testsuite/python/reaction_trivial.py @@ -0,0 +1,78 @@ +# +# Copyright (C) 2023 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 espressomd +import espressomd.reaction_methods +import unittest as ut +import numpy as np + + +class Test(ut.TestCase): + + """Test the trivial reaction 1A <-> 1B""" + + system = espressomd.System(box_l=[1., 1., 1.]) + system.time_step = 0.1 + system.cell_system.skin = 0.1 + np.random.seed(42) + + def test_equilibration(self): + system = self.system + N0 = 40 + c0 = 0.001 + types = {"A": 0, "B": 1} + system.box_l = np.ones(3) * np.cbrt(N0 / c0) + RE = espressomd.reaction_methods.ReactionEnsemble( + seed=42, kT=1., exclusion_range=1., search_algorithm="parallel") + RE.set_non_interacting_type(type=max(types.values()) + 1) + system.part.add( + pos=np.random.random((N0, 3)) * system.box_l, + type=N0 * [types["A"]]) + RE.add_reaction( + gamma=1., + reactant_types=[types["A"]], + reactant_coefficients=[1], + product_types=[types["B"]], + product_coefficients=[1], + default_charges={types["A"]: 0, types["B"]: 0}) + + # warmup + system.thermostat.set_langevin(kT=1., gamma=3., seed=42) + system.integrator.run(200) + RE.reaction(steps=5 * N0) + + # sampling + average_NA = 0.0 + num_samples = 100 + for _ in range(num_samples): + RE.reaction(steps=10) + system.integrator.run(20) + average_NA += system.number_of_particles(type=types["A"]) + average_NA /= num_samples + + alpha = average_NA / float(N0) + rate0 = RE.get_acceptance_rate_reaction(reaction_id=0) + rate1 = RE.get_acceptance_rate_reaction(reaction_id=1) + self.assertAlmostEqual(alpha, 0.50, delta=0.01) + self.assertAlmostEqual(rate0, 0.85, delta=0.05) + self.assertAlmostEqual(rate1, 0.85, delta=0.05) + + +if __name__ == "__main__": + ut.main() From da3e66f54eb8f0be519bccd5c84223d6a462f620 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 14 Sep 2023 21:12:48 +0200 Subject: [PATCH 5/6] core: Fix energy/pressure observable bug --- src/core/Observable_stat.cpp | 8 +++---- src/core/Observable_stat.hpp | 20 +++++++++-------- src/core/constraints/ShapeBasedConstraint.cpp | 6 +++-- src/core/energy_inline.hpp | 2 +- src/core/pressure_inline.hpp | 6 ++--- .../EspressoSystemStandAlone_test.cpp | 12 ++++++---- testsuite/python/analyze_energy.py | 22 +++++++++++++++---- 7 files changed, 48 insertions(+), 28 deletions(-) diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 252296760c4..f7fd6f5abbf 100644 --- a/src/core/Observable_stat.cpp +++ b/src/core/Observable_stat.cpp @@ -58,8 +58,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size) constexpr std::size_t n_ext_fields = 1; // reduction over all fields constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions - auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb + - n_dipolar + n_vs + n_ext_fields; + auto const n_elements = n_kinetic + n_bonded + 2ul * n_non_bonded + + n_coulomb + n_dipolar + n_vs + n_ext_fields; m_data = std::vector(m_chunk_size * n_elements); // spans for the different contributions @@ -78,8 +78,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size) } Utils::Span -Observable_stat::non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const { +Observable_stat::get_non_bonded_contribution(Utils::Span base_pointer, + int type1, int type2) const { auto const offset = static_cast(Utils::upper_triangular( std::min(type1, type2), std::max(type1, type2), max_seen_particle_type)); return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size}; diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index b34157004bc..b84bcfdda00 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -38,8 +38,9 @@ class Observable_stat { std::size_t m_chunk_size; /** Get contribution from a non-bonded interaction */ - Utils::Span non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const; + Utils::Span + get_non_bonded_contribution(Utils::Span base_pointer, int type1, + int type2) const; public: explicit Observable_stat(std::size_t chunk_size); @@ -91,29 +92,30 @@ class Observable_stat { return {bonded.data() + offset, m_chunk_size}; } - void add_non_bonded_contribution(int type1, int type2, + void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, Utils::Span data) { assert(data.size() == m_chunk_size); - auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter; - auto const dest = non_bonded_contribution(source, type1, type2); + auto const span = (molid1 == molid2) ? non_bonded_intra : non_bonded_inter; + auto const dest = get_non_bonded_contribution(span, type1, type2); boost::transform(dest, data, dest.begin(), std::plus<>{}); } - void add_non_bonded_contribution(int type1, int type2, double data) { - add_non_bonded_contribution(type1, type2, {&data, 1}); + void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2, + double data) { + add_non_bonded_contribution(type1, type2, molid1, molid2, {&data, 1}); } /** Get contribution from a non-bonded intramolecular interaction */ Utils::Span non_bonded_intra_contribution(int type1, int type2) const { - return non_bonded_contribution(non_bonded_intra, type1, type2); + return get_non_bonded_contribution(non_bonded_intra, type1, type2); } /** Get contribution from a non-bonded intermolecular interaction */ Utils::Span non_bonded_inter_contribution(int type1, int type2) const { - return non_bonded_contribution(non_bonded_inter, type1, type2); + return get_non_bonded_contribution(non_bonded_inter, type1, type2); } /** MPI reduction. */ diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 8a3a1b85542..9b4e2fe2107 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -164,7 +164,9 @@ void ShapeBasedConstraint::add_energy(const Particle &p, } } // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) - if (part_rep.type() >= 0) - obs_energy.add_non_bonded_contribution(p.type(), part_rep.type(), energy); + if (part_rep.type() >= 0) { + obs_energy.add_non_bonded_contribution( + p.type(), part_rep.type(), p.mol_id(), part_rep.mol_id(), energy); + } } } // namespace Constraints diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index eaa24054fd3..f4c09b93a7e 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -183,7 +183,7 @@ inline void add_non_bonded_pair_energy( if (do_nonbonded(p1, p2)) #endif obs_energy.add_non_bonded_contribution( - p1.type(), p2.type(), + p1.type(), p2.type(), p1.mol_id(), p2.mol_id(), calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist, coulomb_kernel)); diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index a1f6354eb89..7a6ef6589aa 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -70,10 +70,8 @@ inline void add_non_bonded_pair_virials( .f + calc_non_central_force(p1, p2, ia_params, d, dist).f; auto const stress = Utils::tensor_product(d, force); - - auto const type1 = p1.mol_id(); - auto const type2 = p2.mol_id(); - obs_pressure.add_non_bonded_contribution(type1, type2, flatten(stress)); + obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(), + p2.mol_id(), flatten(stress)); } #ifdef ELECTROSTATICS diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index abd19fb80f8..d5d23a38574 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -114,6 +114,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { BOOST_REQUIRE_GE(get_particle_node_parallel(pid2), rank ? -1 : 1); BOOST_REQUIRE_GE(get_particle_node_parallel(pid3), rank ? -1 : 1); } + set_particle_property(pid1, &Particle::mol_id, type_a); + set_particle_property(pid2, &Particle::mol_id, type_b); + set_particle_property(pid3, &Particle::mol_id, type_b); auto const reset_particle_positions = [&start_positions]() { for (auto const &kv : start_positions) { @@ -225,10 +228,10 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { on_non_bonded_ia_change(); // matrix indices and reference energy value - auto const max_type = type_b + 1; - auto const n_pairs = Utils::upper_triangular(type_b, type_b, max_type) + 1; - auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, max_type); - auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, max_type); + auto const size = std::max(type_a, type_b) + 1; + auto const n_pairs = Utils::upper_triangular(type_b, type_b, size) + 1; + auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, size); + auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, size); auto const frac6 = Utils::int_pow<6>(sig / r_off); auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift); @@ -236,6 +239,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { auto const obs_energy = calculate_energy(); if (rank == 0) { for (int i = 0; i < n_pairs; ++i) { + // particles were set up with type == mol_id auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10); diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py index cbe2e829e16..9b1eeff8063 100644 --- a/testsuite/python/analyze_energy.py +++ b/testsuite/python/analyze_energy.py @@ -49,8 +49,8 @@ def setUpClass(cls): cls.system.bonded_inter[5] = cls.harmonic def setUp(self): - self.system.part.add(pos=[1, 2, 2], type=0) - self.system.part.add(pos=[5, 2, 2], type=0) + self.system.part.add(pos=[1, 2, 2], type=0, mol_id=6) + self.system.part.add(pos=[5, 2, 2], type=0, mol_id=6) def tearDown(self): self.system.part.clear() @@ -86,12 +86,14 @@ def test_non_bonded(self): self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7) self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7) self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7) + self.assertAlmostEqual(energy["non_bonded_intra"], 1., delta=1e-7) + self.assertAlmostEqual(energy["non_bonded_inter"], 0., delta=1e-7) # Test the single particle energy function self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum( [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7) # add another pair of particles - self.system.part.add(pos=[3, 2, 2], type=1) - self.system.part.add(pos=[4, 2, 2], type=1) + self.system.part.add(pos=[3, 2, 2], type=1, mol_id=7) + self.system.part.add(pos=[4, 2, 2], type=1, mol_id=7) energy = self.system.analysis.energy() self.assertAlmostEqual(energy["total"], 3., delta=1e-7) self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7) @@ -101,6 +103,18 @@ def test_non_bonded(self): self.assertAlmostEqual(energy["non_bonded", 0, 0] + energy["non_bonded", 0, 1] + energy["non_bonded", 1, 1], energy["total"], delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 0, 0], 1., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 1, 1], 1., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_intra", 0, 1], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 0, 0], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 1, 1], 0., delta=1e-7) + self.assertAlmostEqual( + energy["non_bonded_inter", 0, 1], 1., delta=1e-7) # Test the single particle energy function self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum( [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7) From ffac3db63b4e97f69ef057be10fcc4510baf160f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 19 Sep 2023 17:12:11 +0200 Subject: [PATCH 6/6] doc: Update installation instructions --- doc/sphinx/installation.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 57644a6d958..83b894d1b70 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -8,10 +8,10 @@ This chapter will describe how to get, compile and run the software. |es| releases are available as source code packages from the homepage [1]_. This is where new users should get the code. The code within release packages is tested and known to run on a number of platforms. -Alternatively, people that want to use the newest features of |es| or that -want to start contributing to the software can instead obtain the +Alternatively, people who want to use the newest features of |es| or +start contributing to the software can instead obtain the current development code via the version control system software [2]_ -from |es|'s project page at Github [3]_. This code might be not as well +from |es|'s project page at GitHub [3]_. This code might be not as well tested and documented as the release code; it is recommended to use this code only if you have already gained some experience in using |es|. @@ -93,7 +93,7 @@ To compile |es| on Ubuntu 22.04 LTS, install the following dependencies: .. code-block:: bash sudo apt install build-essential cmake cython3 python3-pip python3-numpy \ - libboost-all-dev openmpi-common fftw3-dev libhdf5-dev libhdf5-openmpi-dev \ + libboost-all-dev openmpi-common fftw3-dev libfftw3-mpi-dev libhdf5-dev libhdf5-openmpi-dev \ python3-scipy python3-opengl libgsl-dev freeglut3 Optionally the ccmake utility can be installed for easier configuration: @@ -221,7 +221,7 @@ To use Jupyter Notebook, install the following packages: .. code-block:: bash - pip3 install --user nbformat notebook 'jupyter_contrib_nbextensions==0.5.1' + pip3 install --user 'nbformat==5.1.3' 'nbconvert==6.4.5' 'notebook==6.4.8' 'jupyter_contrib_nbextensions==0.5.1' jupyter contrib nbextension install --user jupyter nbextension enable rubberband/main jupyter nbextension enable exercise2/main