Skip to content

Commit

Permalink
core: Extract exclusion range neighbor search
Browse files Browse the repository at this point in the history
  • Loading branch information
jngrad committed Jul 28, 2023
1 parent a770f49 commit a0359cb
Show file tree
Hide file tree
Showing 18 changed files with 474 additions and 195 deletions.
38 changes: 27 additions & 11 deletions doc/sphinx/reaction_methods.rst
Original file line number Diff line number Diff line change
Expand Up @@ -353,14 +353,30 @@ The Monte Carlo (MC) sampling of the reaction can be coupled with a configurati
For non-interacting systems this coupling is not an issue, but for interacting systems the insertion of new particles
can lead to instabilities in the MD integration ultimately leading to a crash of the simulation.

This integration instabilities can be avoided by defining a distance around the particles which already exist in the system
where new particles will not be inserted, which is defined by the required keyword ``exclusion_range``.
This prevents big overlaps with the newly inserted particles, avoiding too big forces between particles, which prevents the MD integration from crashing.
The value of the exclusion range does not affect the limiting result and it only affects the convergence and the stability of the integration. For interacting systems,
it is usually a good practice to choose the exclusion range such that it is comparable to the diameter of the particles.

If particles with significantly different sizes are present, it is desired to define a different exclusion range for each pair of particle types. This can be done by
defining an exclusion radius per particle type by using the optional argument ``exclusion_radius_per_type``. Then, their exclusion range is calculated using
the Lorentz-Berthelot combination rule, *i.e.* ``exclusion_range = exclusion_radius_per_type[particle_type_1] + exclusion_radius_per_type[particle_type_2]``.
If the exclusion radius of one particle type is not defined, the value of the parameter provided in ``exclusion_range`` is used by default.
If the value in ``exclusion_radius_per_type`` is equal to 0, then the exclusion range of that particle type with any other particle is 0.
This integration instabilities can be avoided by defining a distance around
the particles which already exist in the system where new particles will not
be inserted, which is defined by the required keyword ``exclusion_range``.
This prevents big overlaps with the newly inserted particles, which would
otherwise cause large forces between particles and crash the MD integrator.
The value of the exclusion range does not affect the limiting result
and it only affects the convergence and the stability of the integration.
For interacting systems, it is usually a good practice to choose the exclusion
range such that it is comparable to the diameter of the particles.

If particles with significantly different sizes are present, it is desirable
to define a different exclusion range for each pair of particle types.
This can be done by defining an exclusion radius per particle type via
the optional argument ``exclusion_radius_per_type``.
Then, their exclusion range is calculated using the Lorentz-Berthelot
combination rule, *i.e.*

.. code-block:: python
exclusion_range = exclusion_radius_per_type[particle_type_1] + \
exclusion_radius_per_type[particle_type_2]
If the exclusion radius of one particle type is not defined, the value
of the parameter provided in ``exclusion_range`` is used by default.
If the value in ``exclusion_radius_per_type`` is equal to 0,
then the exclusion range of that particle type with any other particle is 0.
For more detail, see :class:`espressomd.cell_system.ExclusionRadius`.
6 changes: 4 additions & 2 deletions src/core/reaction_methods/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,10 @@
#

target_sources(
espresso_core PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp
${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp)
espresso_core
PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp
${CMAKE_CURRENT_SOURCE_DIR}/ExclusionRadius.cpp
${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp)

if(ESPRESSO_BUILD_TESTS)
add_subdirectory(tests)
Expand Down
12 changes: 6 additions & 6 deletions src/core/reaction_methods/ConstantpHEnsemble.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,11 @@
#ifndef REACTION_METHODS_CONSTANT_PH_ENSEMBLE_HPP
#define REACTION_METHODS_CONSTANT_PH_ENSEMBLE_HPP

#include "reaction_methods/ExclusionRadius.hpp"
#include "reaction_methods/ReactionAlgorithm.hpp"

#include <map>
#include <memory>

namespace ReactionMethods {

Expand All @@ -40,12 +42,10 @@ namespace ReactionMethods {
*/
class ConstantpHEnsemble : public ReactionAlgorithm {
public:
ConstantpHEnsemble(
boost::mpi::communicator const &comm, int seed, double kT,
double exclusion_range, double constant_pH,
const std::unordered_map<int, double> &exclusion_radius_per_type)
: ReactionAlgorithm(comm, seed, kT, exclusion_range,
exclusion_radius_per_type),
ConstantpHEnsemble(boost::mpi::communicator const &comm, int seed, double kT,
std::shared_ptr<ExclusionRadius> exclusion,
double constant_pH)
: ReactionAlgorithm(comm, seed, kT, exclusion),
m_constant_pH(constant_pH) {}
double m_constant_pH;
};
Expand Down
154 changes: 154 additions & 0 deletions src/core/reaction_methods/ExclusionRadius.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
/*
* Copyright (C) 2022-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 <http://www.gnu.org/licenses/>.
*/

#include "reaction_methods/ExclusionRadius.hpp"

#include "cells.hpp"
#include "event.hpp"
#include "grid.hpp"
#include "particle_node.hpp"

#include <boost/mpi/collectives/all_reduce.hpp>
#include <boost/mpi/collectives/broadcast.hpp>
#include <boost/mpi/communicator.hpp>

#include <algorithm>
#include <cassert>
#include <functional>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>

static auto get_real_particle(boost::mpi::communicator const &comm, int p_id) {
assert(p_id >= 0);
auto ptr = ::cell_structure.get_local_particle(p_id);
if (ptr != nullptr and ptr->is_ghost()) {
ptr = nullptr;
}
assert(boost::mpi::all_reduce(comm, static_cast<int>(ptr != nullptr),
std::plus<>()) == 1);
return ptr;
}

void ExclusionRadius::set_exclusion_range(double range) {
if (range < 0.) {
throw std::domain_error("Invalid value for exclusion range");
}
exclusion_range = range;
recalc_derived_parameters();
}

void ExclusionRadius::set_exclusion_radius_per_type(map_type const &map) {
for (auto const &[type, exclusion_radius] : map) {
if (exclusion_radius < 0.) {
throw std::domain_error("Invalid exclusion radius for type " +
std::to_string(type) + ": radius " +
std::to_string(exclusion_radius));
}
}
exclusion_radius_per_type = map;
recalc_derived_parameters();
}

void ExclusionRadius::recalc_derived_parameters() {
m_max_exclusion_range = exclusion_range;
for (auto const &item : exclusion_radius_per_type) {
auto const radius = item.second;
m_max_exclusion_range = std::max(m_max_exclusion_range, 2. * radius);
}
}

/**
* Check if an inserted particle is too close to neighboring particles.
*/
bool ExclusionRadius::check_exclusion_range(int p_id, int p_type) {

/* Check the exclusion radius of the inserted particle */
if (exclusion_radius_per_type.count(p_type) != 0) {
if (exclusion_radius_per_type[p_type] == 0.) {
return false;
}
}

auto p1_ptr = get_real_particle(m_comm, p_id);

std::vector<int> particle_ids;
if (neighbor_search_order_n) {
auto all_ids = get_particle_ids_parallel();
/* remove the inserted particle id */
all_ids.erase(std::remove(all_ids.begin(), all_ids.end(), p_id),
all_ids.end());
particle_ids = all_ids;
} else {
on_observable_calc();
auto const local_ids =
get_short_range_neighbors(p_id, m_max_exclusion_range);
assert(p1_ptr == nullptr or !!local_ids);
if (local_ids) {
particle_ids = std::move(*local_ids);
}
}

auto within_exclusion_range = false;
if (p1_ptr != nullptr) {
auto &p1 = *p1_ptr;

/* Check if the inserted particle within any exclusion radius */
for (auto const p2_id : particle_ids) {
if (auto const p2_ptr = ::cell_structure.get_local_particle(p2_id)) {
auto const &p2 = *p2_ptr;
double excluded_distance;
if (exclusion_radius_per_type.count(p_type) == 0 or
exclusion_radius_per_type.count(p2.type()) == 0) {
excluded_distance = exclusion_range;
} else if (exclusion_radius_per_type[p2.type()] == 0.) {
continue;
} else {
excluded_distance = exclusion_radius_per_type[p_type] +
exclusion_radius_per_type[p2.type()];
}

auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm();

if (d_min < excluded_distance) {
within_exclusion_range = true;
break;
}
}
}
if (m_comm.rank() != 0) {
m_comm.send(0, 1, within_exclusion_range);
}
} else if (m_comm.rank() == 0) {
m_comm.recv(boost::mpi::any_source, 1, within_exclusion_range);
}
boost::mpi::broadcast(m_comm, within_exclusion_range, 0);
return within_exclusion_range;
}

bool ExclusionRadius::check_exclusion_range(int pid) {
int type_local = 0;
if (auto p = get_real_particle(m_comm, pid)) {
type_local = p->type();
}
auto const type =
boost::mpi::all_reduce(m_comm, type_local, std::plus<int>());
return check_exclusion_range(pid, type);
}
43 changes: 43 additions & 0 deletions src/core/reaction_methods/ExclusionRadius.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
/*
* Copyright (C) 2022-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 <http://www.gnu.org/licenses/>.
*/

#pragma once

#include <boost/mpi/communicator.hpp>

#include <unordered_map>

class ExclusionRadius {
boost::mpi::communicator const &m_comm;
double m_max_exclusion_range = 0.;
void recalc_derived_parameters();

public:
using map_type = std::unordered_map<int, double>;
explicit ExclusionRadius(boost::mpi::communicator const &comm)
: m_comm{comm} {}
void set_exclusion_range(double range);
void set_exclusion_radius_per_type(map_type const &map);
bool check_exclusion_range(int p_id, int p_type);
bool check_exclusion_range(int pid);

bool neighbor_search_order_n = true;
double exclusion_range = 0.;
map_type exclusion_radius_per_type{};
};
64 changes: 2 additions & 62 deletions src/core/reaction_methods/ReactionAlgorithm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -305,68 +305,8 @@ void ReactionAlgorithm::hide_particle(int p_id, int p_type) const {
* Check if the inserted particle is too close to neighboring particles.
*/
void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) {

/* Check the exclusion radius of the inserted particle */
if (exclusion_radius_per_type.count(p_type) != 0) {
if (exclusion_radius_per_type[p_type] == 0.) {
return;
}
}

auto p1_ptr = get_real_particle(p_id);

std::vector<int> particle_ids;
if (neighbor_search_order_n) {
auto all_ids = get_particle_ids_parallel();
/* remove the inserted particle id */
all_ids.erase(std::remove(all_ids.begin(), all_ids.end(), p_id),
all_ids.end());
particle_ids = all_ids;
} else {
on_observable_calc();
auto const local_ids =
get_short_range_neighbors(p_id, m_max_exclusion_range);
assert(p1_ptr == nullptr or !!local_ids);
if (local_ids) {
particle_ids = std::move(*local_ids);
}
}

if (p1_ptr != nullptr) {
auto &p1 = *p1_ptr;

/* Check if the inserted particle within the exclusion radius of any other
* particle */
for (auto const p2_id : particle_ids) {
if (auto const p2_ptr = ::cell_structure.get_local_particle(p2_id)) {
auto const &p2 = *p2_ptr;
double excluded_distance;
if (exclusion_radius_per_type.count(p_type) == 0 ||
exclusion_radius_per_type.count(p2.type()) == 0) {
excluded_distance = exclusion_range;
} else if (exclusion_radius_per_type[p2.type()] == 0.) {
continue;
} else {
excluded_distance = exclusion_radius_per_type[p_type] +
exclusion_radius_per_type[p2.type()];
}

auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm();

if (d_min < excluded_distance) {
particle_inside_exclusion_range_touched = true;
break;
}
}
}
if (m_comm.rank() != 0) {
m_comm.send(0, 1, particle_inside_exclusion_range_touched);
}
} else if (m_comm.rank() == 0) {
m_comm.recv(boost::mpi::any_source, 1,
particle_inside_exclusion_range_touched);
}
boost::mpi::broadcast(m_comm, particle_inside_exclusion_range_touched, 0);
particle_inside_exclusion_range_touched |=
m_exclusion->check_exclusion_range(p_id, p_type);
}

/**
Expand Down
Loading

0 comments on commit a0359cb

Please sign in to comment.