Skip to content

Commit

Permalink
core: Give MC functions more descriptive names
Browse files Browse the repository at this point in the history
Co-authored-by: Pablo Miguel Blanco Andrés <[email protected]>
  • Loading branch information
jngrad and pm-blanco committed Feb 13, 2023
1 parent 82cda68 commit 4cfc4ce
Show file tree
Hide file tree
Showing 6 changed files with 39 additions and 30 deletions.
10 changes: 5 additions & 5 deletions src/core/reaction_methods/ReactionAlgorithm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -233,7 +233,7 @@ ReactionAlgorithm::get_particle_numbers(SingleReaction const &reaction) const {
}

std::optional<double>
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;
Expand All @@ -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(
Expand Down
28 changes: 16 additions & 12 deletions src/core/reaction_methods/ReactionAlgorithm.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> generic_oneway_reaction_part_1(int reaction_id);
std::optional<double> 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);
Expand Down
7 changes: 4 additions & 3 deletions src/python/espressomd/reaction_methods.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions src/script_interface/reaction_methods/ReactionAlgorithm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>(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<double>(params, "bf");
auto const E_pot_old = get_value<double>(params, "E_pot_old");
auto const E_pot_new = get_value<double>(params, "E_pot_new");
auto const reaction_id = get_value<int>(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;
}
Expand Down
10 changes: 5 additions & 5 deletions src/script_interface/tests/ReactionEnsemble_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;

Expand All @@ -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);

Expand All @@ -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
Expand All @@ -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
Expand Down
4 changes: 4 additions & 0 deletions testsuite/python/reaction_methods_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.)
Expand Down Expand Up @@ -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):
Expand Down

0 comments on commit 4cfc4ce

Please sign in to comment.