From 4cfc4cec6f21bbbe8e9c36a334221f11e971c983 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 13 Feb 2023 18:15:31 +0100 Subject: [PATCH] core: Give MC functions more descriptive names MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Pablo Miguel Blanco Andrés --- .../reaction_methods/ReactionAlgorithm.cpp | 10 +++---- .../reaction_methods/ReactionAlgorithm.hpp | 28 +++++++++++-------- src/python/espressomd/reaction_methods.py | 7 +++-- .../reaction_methods/ReactionAlgorithm.cpp | 10 +++---- .../tests/ReactionEnsemble_test.cpp | 10 +++---- .../python/reaction_methods_interface.py | 4 +++ 6 files changed, 39 insertions(+), 30 deletions(-) diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 78cedbab7ba..f723af7a90d 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -233,7 +233,7 @@ ReactionAlgorithm::get_particle_numbers(SingleReaction const &reaction) const { } std::optional -ReactionAlgorithm::generic_oneway_reaction_part_1(int reaction_id) { +ReactionAlgorithm::create_new_trial_state(int reaction_id) { auto &reaction = *reactions[reaction_id]; reaction.tried_moves++; particle_inside_exclusion_range_touched = false; @@ -253,10 +253,10 @@ ReactionAlgorithm::generic_oneway_reaction_part_1(int reaction_id) { return {E_pot_new}; } -double ReactionAlgorithm::generic_oneway_reaction_part_2(int reaction_id, - double bf, - double E_pot_old, - double E_pot_new) { +double ReactionAlgorithm::make_reaction_mc_move_attempt(int reaction_id, + double bf, + double E_pot_old, + double E_pot_new) { auto const exponential = std::exp(-(E_pot_new - E_pot_old) / kT); auto &reaction = *reactions[reaction_id]; reaction.accumulator_potential_energy_difference_exponential( diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 9d61219a3aa..e617dbdebca 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -177,39 +177,43 @@ class ReactionAlgorithm { public: /** - * Attempt a reaction move and calculate the new potential energy. + * Carry out a reaction MC move and calculate the new potential energy. * Particles are selected without replacement. + * The previous state of the system is cached. * @returns Potential energy of the system after the move. */ - std::optional generic_oneway_reaction_part_1(int reaction_id); + std::optional create_new_trial_state(int reaction_id); /** - * Accept or reject moves made by @ref generic_oneway_reaction_part_1 based - * on a probability acceptance @c bf. + * Accept or reject a reaction MC move made by @ref create_new_trial_state + * based on a probability acceptance @c bf. + * The previous state of the system is either restored from the cache if + * the move is rejected, or cleared from the cache if the move is accepted. * @returns Potential energy of the system after the move was accepted or * rejected. */ - double generic_oneway_reaction_part_2(int reaction_id, double bf, - double E_pot_old, double E_pot_new); + double make_reaction_mc_move_attempt(int reaction_id, double bf, + double E_pot_old, double E_pot_new); /** - * Attempt a global MC move for particles of a given type. + * Attempt displacement MC moves for particles of a given type. * Particles are selected without replacement. * @param type Type of particles to move. * @param n_particles Number of particles of that type to move. * @returns true if all moves were accepted. */ bool make_displacement_mc_move_attempt(int type, int n_particles); + + /** @brief Compute the system potential energy. */ + double calculate_potential_energy() const; + +protected: /** - * Perform a global MC move for particles of a given type. + * Carry out displacement MC moves for particles of a given type. * Particles are selected without replacement. * @param type Type of particles to move. * @param n_particles Number of particles of that type to move. */ void displacement_mc_move(int type, int n_particles); - /** @brief Compute the system potential energy. */ - double calculate_potential_energy() const; - -protected: /** @brief Carry out a chemical reaction and save the old system state. */ void make_reaction_attempt(::ReactionMethods::SingleReaction const &reaction, ParticleChanges &bookkeeping); diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index c02445736e7..7d6a45fe510 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -476,12 +476,12 @@ def generic_oneway_reaction(self, reaction_id, E_pot_old): """ try: E_pot_new = self.call_method( - "generic_oneway_reaction_part_1", reaction_id=reaction_id) + "create_new_trial_state", reaction_id=reaction_id) if E_pot_new is None: return E_pot_old E_pot_diff = E_pot_new - E_pot_old bf = self.calculate_acceptance_probability(reaction_id, E_pot_diff) - return self.call_method("generic_oneway_reaction_part_2", + return self.call_method("make_reaction_mc_move_attempt", reaction_id=reaction_id, bf=bf, E_pot_new=E_pot_new, E_pot_old=E_pot_old) except BaseException as err: @@ -515,7 +515,8 @@ def get_status(self): reactions_list.append(reaction) return {"reactions": reactions_list, "kT": self.kT, - "exclusion_range": self.exclusion_range} + "exclusion_range": self.exclusion_range, + "exclusion_radius_per_type": self.exclusion_radius_per_type} @script_interface_register diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp index fbe51a7a707..d11c875a4bf 100644 --- a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp @@ -130,26 +130,26 @@ Variant ReactionAlgorithm::do_call_method(std::string const &name, if (name == "potential_energy") { return RE()->calculate_potential_energy(); } - if (name == "generic_oneway_reaction_part_1") { + if (name == "create_new_trial_state") { auto const reaction_id = get_value(params, "reaction_id"); Variant result{}; context()->parallel_try_catch([&]() { - auto const optional = RE()->generic_oneway_reaction_part_1(reaction_id); + auto const optional = RE()->create_new_trial_state(reaction_id); if (optional) { result = *optional; } }); return result; } - if (name == "generic_oneway_reaction_part_2") { + if (name == "make_reaction_mc_move_attempt") { auto const bf = get_value(params, "bf"); auto const E_pot_old = get_value(params, "E_pot_old"); auto const E_pot_new = get_value(params, "E_pot_new"); auto const reaction_id = get_value(params, "reaction_id"); Variant result; context()->parallel_try_catch([&]() { - result = RE()->generic_oneway_reaction_part_2(reaction_id, bf, E_pot_old, - E_pot_new); + result = RE()->make_reaction_mc_move_attempt(reaction_id, bf, E_pot_old, + E_pot_new); }); return result; } diff --git a/src/script_interface/tests/ReactionEnsemble_test.cpp b/src/script_interface/tests/ReactionEnsemble_test.cpp index 664d819fcbf..dc29e554ef5 100644 --- a/src/script_interface/tests/ReactionEnsemble_test.cpp +++ b/src/script_interface/tests/ReactionEnsemble_test.cpp @@ -182,7 +182,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { espresso::system->set_box_l(box_l); // without reactants, no reaction will take place - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(not result.has_value()); // the reaction was updated @@ -200,7 +200,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { // system auto const energy_ref = 0.; - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(result.has_value()); auto const energy_move = *result; @@ -218,7 +218,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { auto const bf = r_algo_si->calculate_acceptance_probability( reaction, energy_move, {{type_D, 1}, {type_E, 0}}); - auto const energy_end = r_algo.generic_oneway_reaction_part_2( + auto const energy_end = r_algo.make_reaction_mc_move_attempt( reaction_id, bf, 0., energy_move); BOOST_CHECK_CLOSE(energy_end, energy_ref, tol); @@ -231,7 +231,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { } { // attempt a second reaction - auto const result = r_algo.generic_oneway_reaction_part_1(reaction_id); + auto const result = r_algo.create_new_trial_state(reaction_id); BOOST_REQUIRE(result.has_value()); // verify bookkeeping @@ -247,7 +247,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) { // force move to be rejected auto const energy_reject = - r_algo.generic_oneway_reaction_part_2(reaction_id, 0., 0.2, 0.1); + r_algo.make_reaction_mc_move_attempt(reaction_id, 0., 0.2, 0.1); BOOST_CHECK_CLOSE(energy_reject, 0.2, tol); // the reaction was updated diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py index 2f677d08f2e..8f2bfcf4276 100644 --- a/testsuite/python/reaction_methods_interface.py +++ b/testsuite/python/reaction_methods_interface.py @@ -93,6 +93,7 @@ def count_by_type(types): list(method.exclusion_radius_per_type.keys()), [2]) self.assertAlmostEqual( method.exclusion_radius_per_type[2], 0.2, delta=1e-10) + exclusion_radius_per_type = {2: 0.2} self.assertAlmostEqual( method.get_volume(), self.system.volume(), delta=1e-10) method.set_volume(volume=1.) @@ -125,6 +126,9 @@ def count_by_type(types): status = method.get_status() self.assertEqual(status['kT'], kT) self.assertEqual(status['exclusion_range'], exclusion_range) + self.assertEqual( + status['exclusion_radius_per_type'], + exclusion_radius_per_type) self.assertEqual(len(status['reactions']), 2) for reaction_flat, params in zip( status['reactions'], reaction_parameters):