Skip to content

Commit

Permalink
refactored MC displacement move of particles
Browse files Browse the repository at this point in the history
  • Loading branch information
@pm-blanco
pm-blanco committed Oct 26, 2022
1 parent 3a8d4e5 commit 7f53fc7
Show file tree
Hide file tree
Showing 7 changed files with 124 additions and 227 deletions.
129 changes: 43 additions & 86 deletions src/core/reaction_methods/ReactionAlgorithm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -547,106 +547,63 @@ void ReactionAlgorithm::move_particle(int p_id, Utils::Vector3d const &new_pos,
set_particle_v(p_id, vel);
}

std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>>
ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {

std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>> old_state;
old_state.reserve(n_particles);

// draw particle ids at random without replacement
int p_id = -1;
std::vector<int> drawn_pids{p_id};
for (int i = 0; i < n_particles; i++) {
// draw a new particle id
while (Utils::contains(drawn_pids, p_id)) {
auto const random_index = i_random(number_of_particles_with_type(type));
p_id = get_random_p_id(type, random_index);
}
drawn_pids.emplace_back(p_id);
// store original state
auto const &p = get_particle_data(p_id);
old_state.emplace_back(std::tuple<int, Utils::Vector3d, Utils::Vector3d>{
p_id, p.pos(), p.v()});
// write new position and new velocity
auto const prefactor = std::sqrt(kT / p.mass());
auto const new_pos = get_random_position_in_box();
move_particle(p_id, new_pos, prefactor);
check_exclusion_range(p_id);
if (particle_inside_exclusion_range_touched) {
break;
}
}

return old_state;
}

/**
* Performs a global MC move for particles of a given type.
* Particles are selected without replacement.
* @param type Type of particles to move.
* @param n_part Number of particles of that type to move.
* @returns true if all moves were accepted.
* Performs particle displacement MC moves in the canonical ensemble
* @param mc_steps Number of trial MC steps
* @param particle_types_to_move List of particle types from which particles
* are selected. If empty, any particle can be chosen by default.
*
*/
bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type,
int n_part) {

if (type < 0) {
throw std::domain_error("Parameter 'type_mc' must be >= 0");
}
if (n_part < 0) {
throw std::domain_error(
"Parameter 'particle_number_to_be_changed' must be >= 0");
}
void ReactionAlgorithm::do_particle_displacement_MC_move(
int mc_steps, std::vector<int> particle_types_to_move) {
std::vector<signed int> ids_to_move;

if (n_part == 0) {
// reject
return false;
if (particle_types_to_move.empty()) {
ids_to_move = get_particle_ids();
} else {
for (signed int p_id : get_particle_ids()) {
if (std::find(
particle_types_to_move.begin(), particle_types_to_move.end(),
get_particle_data(p_id).type()) != particle_types_to_move.end()) {
ids_to_move.push_back(p_id);
}
}
}

m_tried_configurational_MC_moves += 1;
particle_inside_exclusion_range_touched = false;

auto const n_particles_of_type = number_of_particles_with_type(type);
if (n_part > n_particles_of_type) {
// reject
return false;
// If there are no particles available to move, return
if (ids_to_move.empty()) {
return;
}

auto const E_pot_old = mpi_calculate_potential_energy();
for (int step = 0; step < mc_steps; step++) {

auto const original_state = generate_new_particle_positions(type, n_part);

auto const E_pot_new = (particle_inside_exclusion_range_touched)
? std::numeric_limits<double>::max()
: mpi_calculate_potential_energy();
N_trial_particle_displacement_MC_moves += 1;
particle_inside_exclusion_range_touched = false;
int p_id = i_random(static_cast<int>(ids_to_move.size()));
auto const particle = get_particle_data(p_id);
auto const old_position = particle.pos();
auto const E_pot_old = mpi_calculate_potential_energy();
auto const new_position = get_random_position_in_box();

auto const beta = 1.0 / kT;
place_particle(p_id, new_position);
check_exclusion_range(p_id);

// Metropolis algorithm since proposal density is symmetric
auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old)));
auto const E_pot_new = (particle_inside_exclusion_range_touched)
? std::numeric_limits<double>::max()
: mpi_calculate_potential_energy();

// // correct for enhanced proposal of small radii by using the
// // Metropolis-Hastings algorithm for asymmetric proposal densities
// double old_radius =
// std::sqrt(std::pow(particle_positions[0][0]-cyl_x,2) +
// std::pow(particle_positions[0][1]-cyl_y,2));
// double new_radius =
// std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2));
// auto const bf = std::min(1.0,
// exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
// Metropolis algorithm
auto const bf = std::min(1.0, exp(-1.0 * (E_pot_new - E_pot_old) / kT));

// Metropolis-Hastings algorithm for asymmetric proposal density
if (m_uniform_real_distribution(m_generator) < bf) {
// accept
m_accepted_configurational_MC_moves += 1;
return true;
}
// reject: restore original particle properties
for (auto const &item : original_state) {
set_particle_v(std::get<0>(item), std::get<2>(item));
place_particle(std::get<0>(item), std::get<1>(item));
if (m_uniform_real_distribution(m_generator) < bf) {
// accept
N_accepted_particle_displacement_MC_moves += 1;
} else {
// reject: restore original particle position
place_particle(p_id, old_position);
}
}
return false;
}

} // namespace ReactionMethods
16 changes: 8 additions & 8 deletions src/core/reaction_methods/ReactionAlgorithm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ class ReactionAlgorithm {
double volume;
int non_interacting_type = 100;

int m_accepted_configurational_MC_moves = 0;
int m_tried_configurational_MC_moves = 0;
double get_acceptance_rate_configurational_moves() const {
return static_cast<double>(m_accepted_configurational_MC_moves) /
static_cast<double>(m_tried_configurational_MC_moves);
int N_trial_particle_displacement_MC_moves = 0;
int N_accepted_particle_displacement_MC_moves = 0;
double get_acceptance_rate_particle_displacement_MC_moves() const {
return static_cast<double>(N_accepted_particle_displacement_MC_moves) /
static_cast<double>(N_trial_particle_displacement_MC_moves);
}

auto get_kT() const { return kT; }
Expand Down Expand Up @@ -133,7 +133,9 @@ class ReactionAlgorithm {
reactions.erase(reactions.begin() + reaction_id);
}

bool displacement_move_for_particles_of_type(int type, int n_part);
void
do_particle_displacement_MC_move(int mc_steps,
std::vector<int> particle_types_to_move);

bool particle_inside_exclusion_range_touched = false;
bool neighbor_search_order_n = true;
Expand Down Expand Up @@ -164,8 +166,6 @@ class ReactionAlgorithm {
std::tuple<std::vector<StoredParticleProperty>, std::vector<int>,
std::vector<StoredParticleProperty>>
make_reaction_attempt(SingleReaction const &current_reaction);
std::vector<std::tuple<int, Utils::Vector3d, Utils::Vector3d>>
generate_new_particle_positions(int type, int n_particles);
void
restore_properties(std::vector<StoredParticleProperty> const &property_list);
/**
Expand Down
121 changes: 51 additions & 70 deletions src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,6 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
class ReactionAlgorithmTest : public ReactionAlgorithm {
public:
using ReactionAlgorithm::calculate_acceptance_probability;
using ReactionAlgorithm::generate_new_particle_positions;
using ReactionAlgorithm::get_random_position_in_box;
using ReactionAlgorithm::ReactionAlgorithm;
};
Expand All @@ -65,12 +64,13 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
ReactionAlgorithmTest r_algo(42, 1., 0., {});
for (int tried_moves = 1; tried_moves < 5; ++tried_moves) {
for (int accepted_moves = 0; accepted_moves < 5; ++accepted_moves) {
r_algo.m_tried_configurational_MC_moves = tried_moves;
r_algo.m_accepted_configurational_MC_moves = accepted_moves;
r_algo.N_trial_particle_displacement_MC_moves = tried_moves;
r_algo.N_accepted_particle_displacement_MC_moves = accepted_moves;
auto const ref_rate = static_cast<double>(accepted_moves) /
static_cast<double>(tried_moves);
BOOST_CHECK_CLOSE(r_algo.get_acceptance_rate_configurational_moves(),
ref_rate, tol);
BOOST_CHECK_CLOSE(
r_algo.get_acceptance_rate_particle_displacement_MC_moves(), ref_rate,
tol);
}
}

Expand Down Expand Up @@ -145,48 +145,54 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
// exception if deleting a non-existent particle
BOOST_CHECK_THROW(r_algo.delete_particle(5), std::runtime_error);

// check particle moves
// check particle displacement Monte Carlo moves
{
// set up particles

// set up one particle

auto const box_l = 1.;
std::vector<std::pair<Utils::Vector3d, Utils::Vector3d>> ref_positions{
{{0.1, 0.2, 0.3}, {+10., +20., +30.}},
{{0.4, 0.5, 0.6}, {-10., -20., -30.}}};
place_particle(0, ref_positions[0].first);
place_particle(1, ref_positions[1].first);
set_particle_v(0, ref_positions[0].second);
set_particle_v(1, ref_positions[1].second);
std::vector<Utils::Vector3d> ref_positions{{0.1, 0.2, 0.3},
{0.4, 0.5, 0.6}};
place_particle(0, ref_positions[0]);
set_particle_type(0, type_A);
set_particle_type(1, type_A);
// update particle positions and velocities
BOOST_CHECK(!r_algo.particle_inside_exclusion_range_touched);
r_algo.particle_inside_exclusion_range_touched = false;

// check that the particle moves if there is no restriction

r_algo.do_particle_displacement_MC_move(1, {type_A});
auto const &p = get_particle_data(0);
auto const &new_pos = p.pos();

BOOST_CHECK_GE((new_pos - ref_positions[0]).norm(), tol);

// check that only particles with types in particle_types_to_move are
// allowed to move

place_particle(0, ref_positions[0]);
place_particle(1, ref_positions[1]);

set_particle_type(1, type_B);
r_algo.do_particle_displacement_MC_move(10, {type_A});

auto const &p_1 = get_particle_data(0);
auto const &new_pos_p1 = p_1.pos();
BOOST_CHECK_GE((new_pos_p1 - ref_positions[0]).norm(), tol);

auto const &p_2 = get_particle_data(1);
auto const &new_pos_p2 = p_2.pos();
BOOST_CHECK_EQUAL((new_pos_p2 - ref_positions[1]).norm(), 0);

// force all MC moves to be rejected by picking particles inside
// their exclusion radius

r_algo.exclusion_range = box_l;
auto const bookkeeping = r_algo.generate_new_particle_positions(0, 2);
BOOST_CHECK(r_algo.particle_inside_exclusion_range_touched);
// check moves and bookkeeping
for (auto const &item : bookkeeping) {
auto const pid = std::get<0>(item);
BOOST_REQUIRE(pid == 0 or pid == 1);
auto const ref_old_pos = ref_positions[pid].first;
auto const ref_old_vel = ref_positions[pid].second;
auto const &p = get_particle_data(pid);
auto const &new_pos = p.pos();
auto const &new_vel = p.v();
BOOST_CHECK_EQUAL(std::get<1>(item), ref_old_pos);
BOOST_CHECK_EQUAL(std::get<2>(item), ref_old_vel);
BOOST_CHECK_GE(new_pos, Utils::Vector3d::broadcast(0.));
BOOST_CHECK_LE(new_pos, Utils::Vector3d::broadcast(box_l));
BOOST_CHECK_GE((new_pos - ref_old_pos).norm(), 0.1);
BOOST_CHECK_GE((new_vel - ref_old_vel).norm(), 10.);
}
// cleanup
remove_particle(0);
remove_particle(1);
}
r_algo.particle_inside_exclusion_range_touched = false;

// check Monte Carlo moves
{
// displacement_move(1)
// auto const &p = get_particle_data(0);
// auto const &new_pos = p.pos();
// BOOST_CHECK_LE((new_pos - ref_positions[0]).norm(), tol);

/*
// set up particles
auto const box_l = 1.;
std::vector<Utils::Vector3d> ref_positions{{0.1, 0.2, 0.3},
Expand All @@ -197,8 +203,9 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
set_particle_type(1, type_A);
// check runtime errors when a MC move cannot be physically performed
auto displacement_move =
std::bind(&ReactionAlgorithm::displacement_move_for_particles_of_type,
std::bind(&ReactionAlgorithm::do_particle_displacement_MC_move,
&r_algo, std::placeholders::_1, std::placeholders::_2);
BOOST_REQUIRE(!displacement_move(type_C, 1));
BOOST_REQUIRE(!displacement_move(type_B, 2));
BOOST_REQUIRE(!displacement_move(type_A, 0));
Expand Down Expand Up @@ -231,6 +238,7 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
// cleanup
remove_particle(0);
remove_particle(1);
*/
}

// check random positions generator
Expand Down Expand Up @@ -312,33 +320,6 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) {
exclusion_radius_per_type[type_A] = 0.1;
exclusion_radius_per_type[type_B] = 1;
ReactionAlgorithmTest r_algo(40, 1., 0, exclusion_radius_per_type);

// the new position will always be in the excluded range since the sum of
// the radii of both particle types is larger than box length. The exclusion
// range value should be ignored

r_algo.generate_new_particle_positions(type_B, 1);

BOOST_REQUIRE(r_algo.particle_inside_exclusion_range_touched);

// the new position will never be in the excluded range because the
// exclusion_radius of the particle is 0

r_algo.exclusion_radius_per_type[type_B] = 0;
r_algo.particle_inside_exclusion_range_touched = false;
r_algo.generate_new_particle_positions(type_B, 1);

BOOST_REQUIRE(!r_algo.particle_inside_exclusion_range_touched);
// the new position will never accepted since the value in exclusion_range
// will be used if the particle does not have a defined excluded radius

r_algo.exclusion_range = 1;
r_algo.exclusion_radius_per_type = {{type_A, 0}};
r_algo.generate_new_particle_positions(type_B, 1);

BOOST_REQUIRE(r_algo.particle_inside_exclusion_range_touched);

//
}
}

Expand Down
Loading

0 comments on commit 7f53fc7

Please sign in to comment.