Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactored MC displacement move of particles #4599

Open
wants to merge 1 commit into
base: python
Choose a base branch
from
Open

Refactored MC displacement move of particles #4599

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 43 additions & 86 deletions src/core/reaction_methods/ReactionAlgorithm.cpp
View file
Open in desktop
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 = ids_to_move[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
Open in desktop
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
154 changes: 47 additions & 107 deletions src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp
View file
Open in desktop
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,92 +145,59 @@ 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
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);
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;
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);
}

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

auto const box_l = 1.;
std::vector<Utils::Vector3d> ref_positions{{0.1, 0.2, 0.3},
{0.4, 0.5, 0.6}};
place_particle(0, ref_positions[0]);
place_particle(1, ref_positions[1]);
set_particle_type(0, type_A);
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,
&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));
BOOST_CHECK_THROW(displacement_move(type_A, -2), std::domain_error);
BOOST_CHECK_THROW(displacement_move(-2, 1), std::domain_error);

// check that the particle moves if there is no restriction
r_algo.do_particle_displacement_MC_move(1, {});
BOOST_CHECK_GE((get_particle_data(0).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_B});

BOOST_CHECK_EQUAL((get_particle_data(0).pos() - ref_positions[0]).norm(),
0);
BOOST_CHECK_GE((get_particle_data(1).pos() - ref_positions[1]).norm(), tol);

// check that no particle moves if they do not match types in
// particle_types_to_move

place_particle(0, ref_positions[0]);
place_particle(1, ref_positions[1]);
r_algo.do_particle_displacement_MC_move(10, {3});
BOOST_CHECK_EQUAL((get_particle_data(0).pos() - ref_positions[0]).norm(),
0);
BOOST_CHECK_EQUAL((get_particle_data(1).pos() - 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;
r_algo.particle_inside_exclusion_range_touched = false;
BOOST_REQUIRE(!r_algo.displacement_move_for_particles_of_type(type_A, 2));
// check none of the particles moved
for (auto const pid : {0, 1}) {
auto const ref_old_pos = ref_positions[pid];
auto const &p = get_particle_data(pid);
auto const &new_pos = p.pos();
BOOST_CHECK_LE((new_pos - ref_old_pos).norm(), tol);
}
// force a MC move to be accepted by using a constant Hamiltonian
r_algo.exclusion_range = 0.;
r_algo.particle_inside_exclusion_range_touched = false;
BOOST_REQUIRE(r_algo.displacement_move_for_particles_of_type(type_A, 1));
std::vector<double> distances(2);
// check that only one particle moved
for (auto const pid : {0, 1}) {
auto const &p = get_particle_data(pid);
distances[pid] = (ref_positions[pid] - p.pos()).norm();
}
BOOST_CHECK_LE(std::min(distances[0], distances[1]), tol);
BOOST_CHECK_GE(std::max(distances[0], distances[1]), 0.1);
// cleanup
remove_particle(0);
remove_particle(1);

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

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

BOOST_CHECK_EQUAL((get_particle_data(0).pos() - ref_positions[0]).norm(),
0);
BOOST_CHECK_EQUAL((get_particle_data(1).pos() - ref_positions[1]).norm(),
0);
}

// check random positions generator
Expand Down Expand Up @@ -312,33 +279,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