diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp
index ca0bb4f7a23..7d03b8d865c 100644
--- a/src/core/cell_system/CellStructure.hpp
+++ b/src/core/cell_system/CellStructure.hpp
@@ -692,6 +692,58 @@ struct CellStructure {
return particle_to_cell(p);
}
+
+ /**
+ * @brief Run kernel on all particles inside local cell and its neighbors.
+ *
+ * @param p Particle to find cell for
+ * @param kernel Function with signature double(Particle const&,
+ * Particle const&, Utils::Vector3d const&)
+ * @return false if cell is not found, otherwise true
+ */
+ template
+ bool run_on_particle_short_range_neighbors(Particle const &p,
+ Kernel &kernel) {
+ auto const cell = find_current_cell(p);
+
+ if (cell == nullptr) {
+ return false;
+ }
+
+ auto const maybe_box = decomposition().minimum_image_distance();
+
+ if (maybe_box) {
+ auto const distance_function = detail::MinimalImageDistance{box_geo};
+ short_range_neighbor_loop(p, cell, kernel, distance_function);
+ } else {
+ auto const distance_function = detail::EuclidianDistance{};
+ short_range_neighbor_loop(p, cell, kernel, distance_function);
+ }
+ return true;
+ }
+
+private:
+ template
+ void short_range_neighbor_loop(Particle const &p1, Cell *const cell,
+ Kernel &kernel, DistanceFunc const &df) {
+ /* Iterate over particles inside cell */
+ for (auto const &p2 : cell->particles()) {
+ if (p1.id() != p2.id()) {
+ auto const vec = df(p1, p2).vec21;
+ kernel(p1, p2, vec);
+ }
+ }
+ /* Iterate over all neighbors */
+ for (auto const neighbor : cell->neighbors().all()) {
+ /* Iterate over particles in neighbors */
+ if (neighbor != cell) {
+ for (auto const &p2 : neighbor->particles()) {
+ auto const vec = df(p1, p2).vec21;
+ kernel(p1, p2, vec);
+ }
+ }
+ }
+ }
};
#endif
diff --git a/src/core/cells.cpp b/src/core/cells.cpp
index 0e98eb31d3c..dd0216ce5c5 100644
--- a/src/core/cells.cpp
+++ b/src/core/cells.cpp
@@ -33,6 +33,7 @@
#include "Particle.hpp"
#include "communication.hpp"
+#include "errorhandling.hpp"
#include "event.hpp"
#include "grid.hpp"
#include "integrate.hpp"
@@ -113,8 +114,51 @@ static void search_distance_sanity_check(double const distance) {
std::to_string(range));
}
}
+static void search_neighbors_sanity_check(double const distance) {
+ search_distance_sanity_check(distance);
+ if (cell_structure.decomposition_type() ==
+ CellStructureType::CELL_STRUCTURE_HYBRID) {
+ throw std::runtime_error("Cannot search for neighbors in the hybrid "
+ "decomposition cell system");
+ }
+}
} // namespace detail
+boost::optional>
+mpi_get_short_range_neighbors_local(int const pid, double const distance,
+ bool run_sanity_checks) {
+
+ if (run_sanity_checks) {
+ detail::search_neighbors_sanity_check(distance);
+ }
+ on_observable_calc();
+
+ auto const p = cell_structure.get_local_particle(pid);
+ if (not p or p->is_ghost()) {
+ return {};
+ }
+
+ std::vector ret;
+ auto const cutoff2 = distance * distance;
+ auto kernel = [&ret, cutoff2](Particle const &p1, Particle const &p2,
+ Utils::Vector3d const &vec) {
+ if (vec.norm2() < cutoff2) {
+ ret.emplace_back(p2.id());
+ }
+ };
+ cell_structure.run_on_particle_short_range_neighbors(*p, kernel);
+ return ret;
+}
+
+REGISTER_CALLBACK_ONE_RANK(mpi_get_short_range_neighbors_local)
+
+std::vector mpi_get_short_range_neighbors(int const pid,
+ double const distance) {
+ detail::search_neighbors_sanity_check(distance);
+ return mpi_call(::Communication::Result::one_rank,
+ mpi_get_short_range_neighbors_local, pid, distance, false);
+}
+
std::vector> get_pairs(double const distance) {
detail::search_distance_sanity_check(distance);
auto pairs =
diff --git a/src/core/cells.hpp b/src/core/cells.hpp
index 99b1f6452b3..372ce2d466e 100644
--- a/src/core/cells.hpp
+++ b/src/core/cells.hpp
@@ -56,6 +56,8 @@
#include "Particle.hpp"
+#include
+
#include
#include
@@ -116,6 +118,15 @@ get_pairs_of_types(double distance, std::vector const &types);
/** Check if a particle resorting is required. */
void check_resort_particles();
+/**
+ * @brief Get ids of particles that are within a certain distance
+ * of another particle.
+ */
+std::vector mpi_get_short_range_neighbors(int pid, double distance);
+boost::optional>
+mpi_get_short_range_neighbors_local(int pid, double distance,
+ bool run_sanity_checks);
+
/**
* @brief Find the cell in which a particle is stored.
*
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index 73eb3fdcc36..e5710b90491 100644
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -128,3 +128,36 @@ double observable_compute_energy() {
auto const obs_energy = calculate_energy();
return obs_energy->accumulate(0);
}
+
+double particle_short_range_energy_contribution_local(int pid) {
+ double ret = 0.0;
+
+ if (cell_structure.get_resort_particles()) {
+ cells_update_ghosts(global_ghost_flags());
+ }
+
+ auto const p = cell_structure.get_local_particle(pid);
+
+ if (p) {
+ auto kernel = [&ret](Particle const &p, Particle const &p1,
+ Utils::Vector3d const &vec) {
+#ifdef EXCLUSIONS
+ if (not do_nonbonded(p, p1))
+ return;
+#endif
+ auto const &ia_params = *get_ia_param(p.type(), p1.type());
+ // Add energy for current particle pair to result
+ ret += calc_non_bonded_pair_energy(p, p1, ia_params, vec, vec.norm());
+ };
+ cell_structure.run_on_particle_short_range_neighbors(*p, kernel);
+ }
+ return ret;
+}
+
+REGISTER_CALLBACK_REDUCTION(particle_short_range_energy_contribution_local,
+ std::plus())
+
+double particle_short_range_energy_contribution(int pid) {
+ return mpi_call(Communication::Result::reduction, std::plus(),
+ particle_short_range_energy_contribution_local, pid);
+}
diff --git a/src/core/energy.hpp b/src/core/energy.hpp
index d6e40ca947e..2a4b5fd1eb3 100644
--- a/src/core/energy.hpp
+++ b/src/core/energy.hpp
@@ -42,4 +42,15 @@ double calculate_current_potential_energy_of_system();
/** Helper function for @ref Observables::Energy. */
double observable_compute_energy();
+/**
+ * @brief Compute short-range energy of a particle.
+ *
+ * Iterates through particles inside cell and neighboring cells and computes
+ * energy contribution for a specificparticle.
+ *
+ * @param pid Particle id
+ * @return Non-bonded energy of the particle.
+ */
+double particle_short_range_energy_contribution(int pid);
+
#endif
diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp
index c7839ce1f5e..d8c33fb1507 100644
--- a/src/core/reaction_methods/ReactionAlgorithm.cpp
+++ b/src/core/reaction_methods/ReactionAlgorithm.cpp
@@ -21,6 +21,7 @@
#include "reaction_methods/ReactionAlgorithm.hpp"
+#include "cells.hpp"
#include "energy.hpp"
#include "grid.hpp"
#include "partCfg_global.hpp"
@@ -381,45 +382,46 @@ void ReactionAlgorithm::check_exclusion_range(int inserted_particle_id) {
auto const &inserted_particle = get_particle_data(inserted_particle_id);
- /* Check the excluded radius of the inserted particle */
-
+ /* Check the exclusion radius of the inserted particle */
if (exclusion_radius_per_type.count(inserted_particle.type()) != 0) {
- if (exclusion_radius_per_type[inserted_particle.type()] == 0) {
+ if (exclusion_radius_per_type[inserted_particle.type()] == 0.) {
return;
}
}
- auto particle_ids = get_particle_ids();
- /* remove the inserted particle id*/
- particle_ids.erase(std::remove(particle_ids.begin(), particle_ids.end(),
- inserted_particle_id),
- particle_ids.end());
-
- /* Check if the inserted particle within the excluded_range of any other
- * particle*/
- double excluded_distance;
- for (const auto &particle_id : particle_ids) {
- auto const &already_present_particle = get_particle_data(particle_id);
+ std::vector particle_ids;
+ if (neighbor_search_order_n) {
+ particle_ids = get_particle_ids();
+ /* remove the inserted particle id */
+ particle_ids.erase(std::remove(particle_ids.begin(), particle_ids.end(),
+ inserted_particle_id),
+ particle_ids.end());
+ } else {
+ particle_ids = mpi_get_short_range_neighbors(inserted_particle.identity(),
+ m_max_exclusion_range);
+ }
+
+ /* Check if the inserted particle within the exclusion radius of any other
+ * particle */
+ for (auto const &particle_id : particle_ids) {
+ auto const &p = get_particle_data(particle_id);
+ double excluded_distance;
if (exclusion_radius_per_type.count(inserted_particle.type()) == 0 ||
- exclusion_radius_per_type.count(inserted_particle.type()) == 0) {
+ exclusion_radius_per_type.count(p.type()) == 0) {
excluded_distance = exclusion_range;
- } else if (exclusion_radius_per_type[already_present_particle.type()] ==
- 0.) {
+ } else if (exclusion_radius_per_type[p.type()] == 0.) {
continue;
} else {
- excluded_distance =
- exclusion_radius_per_type[inserted_particle.type()] +
- exclusion_radius_per_type[already_present_particle.type()];
+ excluded_distance = exclusion_radius_per_type[inserted_particle.type()] +
+ exclusion_radius_per_type[p.type()];
}
auto const d_min =
- box_geo
- .get_mi_vector(already_present_particle.r.p, inserted_particle.r.p)
- .norm();
+ box_geo.get_mi_vector(p.pos(), inserted_particle.pos()).norm();
if (d_min < excluded_distance) {
particle_inside_exclusion_range_touched = true;
- return;
+ break;
}
}
}
diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp
index ffc6248a13f..69e12a6b033 100644
--- a/src/core/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/core/reaction_methods/ReactionAlgorithm.hpp
@@ -50,8 +50,9 @@ class ReactionAlgorithm {
public:
ReactionAlgorithm(
int seed, double kT, double exclusion_range,
- const std::unordered_map &exclusion_radius_per_type)
- : m_generator(Random::mt19937(std::seed_seq({seed, seed, seed}))),
+ std::unordered_map const &exclusion_radius_per_type)
+ : kT{kT}, exclusion_range{exclusion_range},
+ m_generator(Random::mt19937(std::seed_seq({seed, seed, seed}))),
m_normal_distribution(0.0, 1.0), m_uniform_real_distribution(0.0, 1.0) {
if (kT < 0.) {
throw std::domain_error("Invalid value for 'kT'");
@@ -59,8 +60,6 @@ class ReactionAlgorithm {
if (exclusion_range < 0.) {
throw std::domain_error("Invalid value for 'exclusion_range'");
}
- this->kT = kT;
- this->exclusion_range = exclusion_range;
set_exclusion_radius_per_type(exclusion_radius_per_type);
update_volume();
}
@@ -100,6 +99,7 @@ class ReactionAlgorithm {
void update_volume();
void
set_exclusion_radius_per_type(std::unordered_map const &map) {
+ auto max_exclusion_range = exclusion_range;
for (auto const &item : map) {
auto const type = item.first;
auto const exclusion_radius = item.second;
@@ -108,8 +108,11 @@ class ReactionAlgorithm {
std::to_string(type) + ": radius " +
std::to_string(exclusion_radius));
}
+ max_exclusion_range =
+ std::max(max_exclusion_range, 2. * exclusion_radius);
}
exclusion_radius_per_type = map;
+ m_max_exclusion_range = max_exclusion_range;
}
virtual int do_reaction(int reaction_steps);
@@ -134,6 +137,7 @@ class ReactionAlgorithm {
int particle_number_of_type);
bool particle_inside_exclusion_range_touched = false;
+ bool neighbor_search_order_n = true;
protected:
std::vector m_empty_p_ids_smaller_than_max_seen_particle;
@@ -178,7 +182,6 @@ class ReactionAlgorithm {
bool
all_reactant_particles_exist(SingleReaction const ¤t_reaction) const;
-protected:
virtual double
calculate_acceptance_probability(SingleReaction const &, double, double,
std::map const &) const {
@@ -210,6 +213,7 @@ class ReactionAlgorithm {
double m_cyl_y = -10.0;
double m_slab_start_z = -10.0;
double m_slab_end_z = -10.0;
+ double m_max_exclusion_range = 0.;
protected:
Utils::Vector3d get_random_position_in_box();
diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd
index f93dd857166..af70f233f38 100644
--- a/src/python/espressomd/analyze.pxd
+++ b/src/python/espressomd/analyze.pxd
@@ -40,6 +40,10 @@ cdef extern from "" namespace "std" nogil:
array2() except+
double & operator[](size_t)
+cdef extern from "particle_data.hpp":
+ ctypedef struct particle "Particle"
+ const particle & get_particle_data(int pid) except +
+
cdef extern from "PartCfg.hpp":
cppclass PartCfg:
pass
@@ -92,6 +96,7 @@ cdef extern from "pressure.hpp":
cdef extern from "energy.hpp":
cdef shared_ptr[Observable_stat] calculate_energy()
double calculate_current_potential_energy_of_system()
+ double particle_short_range_energy_contribution(int pid)
cdef extern from "dpd.hpp":
Vector9d dpd_stress()
diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx
index 0c64825d432..d7523cff457 100644
--- a/src/python/espressomd/analyze.pyx
+++ b/src/python/espressomd/analyze.pyx
@@ -29,6 +29,7 @@ from .utils cimport Vector3i, Vector3d, Vector9d
from .utils cimport make_Vector3d
from .utils cimport make_array_locked
from .utils cimport create_nparray_from_double_array
+from .particle_data cimport particle, ParticleHandle
def autocorrelation(time_series):
@@ -344,6 +345,24 @@ class Analysis:
utils.handle_errors("calculate_energy() failed")
return obs
+ def particle_energy(self, particle):
+ """
+ Calculate the non-bonded energy of a single given particle.
+
+ Parameters
+ ----------
+ particle : :class:`~espressomd.particle_data.ParticleHandle`
+
+ Returns
+ -------
+ :obj: `float`
+ non-bonded energy of that particle
+
+ """
+ energy_contribution = particle_short_range_energy_contribution(
+ particle.id)
+ return energy_contribution
+
def calc_re(self, chain_start=None, number_of_chains=None,
chain_length=None):
"""
diff --git a/src/python/espressomd/cell_system.py b/src/python/espressomd/cell_system.py
index 786f3d7ec04..f72b911f8d8 100644
--- a/src/python/espressomd/cell_system.py
+++ b/src/python/espressomd/cell_system.py
@@ -17,6 +17,7 @@
# along with this program. If not, see .
#
from . import utils
+from . import particle_data
from .script_interface import ScriptInterfaceHelper, script_interface_register
@@ -177,6 +178,38 @@ def get_pairs(self, distance, types='all'):
pairs = self.call_method("get_pairs", distance=distance, types=types)
return [tuple(pair) for pair in pairs]
+ def get_neighbors(self, particle, distance):
+ """
+ Get neighbors of a given particle up to a certain distance.
+
+ The choice of :ref:`cell systems | ` has an impact on
+ how far the algorithm can detect particle pairs:
+
+ * N-square: no restriction on the search distance, no double counting
+ if search distance is larger than the box size
+ * regular decomposition: the search distance is bounded by half
+ the local cell geometry
+ * hybrid decomposition: not supported
+
+ Parameters
+ ----------
+ particle : :class:`~espressomd.particle_data.ParticleHandle`
+ distance : :obj:`float`
+ Pairs of particles closer than ``distance`` are found.
+
+ Returns
+ -------
+ (N,) array_like of :obj:`int`
+ The list of neighbor particles surrounding the particle
+
+ """
+ utils.check_type_or_throw_except(
+ particle, 1, particle_data.ParticleHandle, "Parameter 'particle' must be a particle")
+ utils.check_type_or_throw_except(
+ distance, 1, float, "Parameter 'distance' must be a float")
+ return self.call_method(
+ "get_neighbors", distance=distance, pid=particle.id)
+
def non_bonded_loop_trace(self):
pairs = self.call_method("non_bonded_loop_trace")
res = []
diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py
index 1e237a22640..66a678cc418 100644
--- a/src/python/espressomd/reaction_methods.py
+++ b/src/python/espressomd/reaction_methods.py
@@ -72,6 +72,13 @@ class ReactionAlgorithm(ScriptInterfaceHelper):
Initial counter value (or seed) of the Mersenne Twister RNG.
exclusion_radius_per_type : :obj:`dict`, optional
Mapping of particle types to exclusion radii.
+ search_algorithm : :obj:`str`
+ Pair search algorithm. Default is ``"order_n"``, which evaluates the
+ distance between the inserted particle and all other particles in the
+ system, which scales with O(N). For MPI-parallel simulations, the
+ ``"parallel"`` method is faster. The ``"parallel"`` method is not
+ recommended for simulations on 1 MPI rank, since it comes with the
+ overhead of a ghost particle update.
Methods
-------
@@ -286,7 +293,8 @@ def __init__(self, **kwargs):
utils.check_valid_keys(self.valid_keys(), kwargs.keys())
def valid_keys(self):
- return {"kT", "exclusion_range", "seed", "exclusion_radius_per_type"}
+ return {"kT", "exclusion_range", "seed",
+ "exclusion_radius_per_type", "search_algorithm"}
def required_keys(self):
return {"kT", "exclusion_range", "seed"}
@@ -409,7 +417,7 @@ class ConstantpHEnsemble(ReactionAlgorithm):
def valid_keys(self):
return {"kT", "exclusion_range", "seed",
- "constant_pH", "exclusion_radius_per_type"}
+ "constant_pH", "exclusion_radius_per_type", "search_algorithm"}
def required_keys(self):
return {"kT", "exclusion_range", "seed", "constant_pH"}
diff --git a/src/script_interface/cell_system/CellSystem.hpp b/src/script_interface/cell_system/CellSystem.hpp
index c2f6d49df4a..8346bf44d14 100644
--- a/src/script_interface/cell_system/CellSystem.hpp
+++ b/src/script_interface/cell_system/CellSystem.hpp
@@ -217,6 +217,27 @@ class CellSystem : public AutoParameters {
});
return out;
}
+ if (name == "get_neighbors") {
+ std::vector> neighbors_global;
+ context()->parallel_try_catch([&neighbors_global, ¶ms]() {
+ auto const dist = get_value(params, "distance");
+ auto const pid = get_value(params, "pid");
+ auto const ret = mpi_get_short_range_neighbors_local(pid, dist, true);
+ std::vector neighbors_local;
+ if (ret) {
+ neighbors_local = *ret;
+ }
+ boost::mpi::gather(comm_cart, neighbors_local, neighbors_global, 0);
+ });
+ std::vector neighbors;
+ for (auto const &neighbors_local : neighbors_global) {
+ if (not neighbors_local.empty()) {
+ neighbors = neighbors_local;
+ break;
+ }
+ }
+ return neighbors;
+ }
if (name == "non_bonded_loop_trace") {
std::vector out;
auto const pair_list = non_bonded_loop_trace();
diff --git a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
index 47c0a672e3f..93bde19c9c4 100644
--- a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
@@ -56,6 +56,9 @@ class ConstantpHEnsemble : public ReactionAlgorithm {
get_value(params, "constant_pH"),
get_value_or>(
params, "exclusion_radius_per_type", {}));
+ do_set_parameter("search_algorithm",
+ Variant{get_value_or(
+ params, "search_algorithm", "order_n")});
}
private:
diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
index a21d125ea6f..eb35edf59f3 100644
--- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
@@ -67,6 +67,24 @@ class ReactionAlgorithm : public AutoParameters {
return out;
}},
{"kT", AutoParameter::read_only, [this]() { return RE()->get_kT(); }},
+ {"search_algorithm",
+ [this](Variant const &v) {
+ auto const key = get_value(v);
+ if (key == "order_n") {
+ RE()->neighbor_search_order_n = true;
+ } else if (key == "parallel") {
+ RE()->neighbor_search_order_n = false;
+ } else {
+ throw std::invalid_argument("Unknown search algorithm '" + key +
+ "'");
+ }
+ },
+ [this]() {
+ if (RE()->neighbor_search_order_n) {
+ return std::string("order_n");
+ }
+ return std::string("parallel");
+ }},
{"exclusion_range", AutoParameter::read_only,
[this]() { return RE()->get_exclusion_range(); }},
{"exclusion_radius_per_type",
diff --git a/src/script_interface/reaction_methods/ReactionEnsemble.hpp b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
index 20478e1e8b6..57bc0be625b 100644
--- a/src/script_interface/reaction_methods/ReactionEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
@@ -45,6 +45,9 @@ class ReactionEnsemble : public ReactionAlgorithm {
get_value(params, "exclusion_range"),
get_value_or>(
params, "exclusion_radius_per_type", {}));
+ do_set_parameter("search_algorithm",
+ Variant{get_value_or(
+ params, "search_algorithm", "order_n")});
}
private:
diff --git a/src/script_interface/reaction_methods/WidomInsertion.hpp b/src/script_interface/reaction_methods/WidomInsertion.hpp
index 343679a0184..4dab2dcf31d 100644
--- a/src/script_interface/reaction_methods/WidomInsertion.hpp
+++ b/src/script_interface/reaction_methods/WidomInsertion.hpp
@@ -40,6 +40,15 @@ class WidomInsertion : public ReactionAlgorithm {
return m_re;
}
+ WidomInsertion() {
+ add_parameters({{"search_algorithm",
+ [](Variant const &) {
+ throw std::runtime_error(
+ "No search algorithm for WidomInsertion");
+ },
+ []() { return none; }}});
+ }
+
void do_construct(VariantMap const ¶ms) override {
m_re = std::make_shared<::ReactionMethods::WidomInsertion>(
get_value(params, "seed"), get_value(params, "kT"), 0.,
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index 3c020d018e9..0d34b2bd196 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -107,6 +107,8 @@ checkpoint_test(MODES therm_sdm__int_sdm)
python_test(FILE bond_breakage.py MAX_NUM_PROC 4)
python_test(FILE cell_system.py MAX_NUM_PROC 4)
+python_test(FILE get_neighbors.py MAX_NUM_PROC 4)
+python_test(FILE get_neighbors.py MAX_NUM_PROC 3 SUFFIX 3_cores)
python_test(FILE tune_skin.py MAX_NUM_PROC 1)
python_test(FILE constraint_homogeneous_magnetic_field.py MAX_NUM_PROC 4)
python_test(FILE cutoffs.py MAX_NUM_PROC 4)
diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py
index 6e1e50cc589..4fb1aa9e474 100644
--- a/testsuite/python/analyze_energy.py
+++ b/testsuite/python/analyze_energy.py
@@ -81,6 +81,9 @@ def test_non_bonded(self):
self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
+ # Test the single particle energy function
+ self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
+ [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
# add another pair of particles
self.system.part.add(pos=[3, 2, 2], type=1)
self.system.part.add(pos=[4, 2, 2], type=1)
@@ -93,6 +96,9 @@ def test_non_bonded(self):
self.assertAlmostEqual(energy["non_bonded", 0, 0]
+ energy["non_bonded", 0, 1]
+ energy["non_bonded", 1, 1], energy["total"], delta=1e-7)
+ # Test the single particle energy function
+ self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
+ [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
def test_bonded(self):
p0, p1 = self.system.part.all()
@@ -134,6 +140,8 @@ def test_all(self):
self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
self.assertAlmostEqual(energy["bonded"], 3. / 2., delta=1e-7)
self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
+ self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
+ [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
# two bonds
p1.add_bond((self.harmonic, p0))
energy = self.system.analysis.energy()
@@ -141,6 +149,8 @@ def test_all(self):
self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
self.assertAlmostEqual(energy["bonded"], 3., delta=1e-7)
self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
+ self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
+ [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
# add another pair of particles
self.system.part.add(pos=[1, 5, 5], type=1)
self.system.part.add(pos=[2, 5, 5], type=1)
@@ -150,6 +160,8 @@ def test_all(self):
self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
self.assertAlmostEqual(energy["bonded"], 3., delta=1e-7)
self.assertAlmostEqual(energy["non_bonded"], 1. + 1., delta=1e-7)
+ self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
+ [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
@utx.skipIfMissingFeatures(["ELECTROSTATICS", "P3M"])
def test_electrostatics(self):
diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py
index 07ed7fdd418..a3a3fc56884 100644
--- a/testsuite/python/constant_pH_stats.py
+++ b/testsuite/python/constant_pH_stats.py
@@ -23,7 +23,7 @@
import espressomd.reaction_methods
-class ReactionEnsembleTest(ut.TestCase):
+class Test(ut.TestCase):
"""Test the core implementation of the constant pH reaction ensemble."""
@@ -39,20 +39,21 @@ class ReactionEnsembleTest(ut.TestCase):
types["A-"]: -1,
types["H+"]: +1,
}
- temperature = 1.0
+ temperature = 1.
# choose target alpha not too far from 0.5 to get good statistics in a
# small number of steps
pKa_minus_pH = -0.2
- pH = 2
+ pH = 2.
pKa = pKa_minus_pH + pH
- Ka = 10**(-pKa)
- box_l = (N0 / c0)**(1.0 / 3.0)
+ Ka = 10.**(-pKa)
+ box_l = np.cbrt(N0 / c0)
system = espressomd.System(box_l=[box_l, box_l, box_l])
np.random.seed(69) # make reaction code fully deterministic
system.cell_system.skin = 0.4
system.time_step = 0.01
RE = espressomd.reaction_methods.ConstantpHEnsemble(
- kT=1.0, exclusion_range=1, seed=44, constant_pH=pH)
+ kT=1., exclusion_range=1., seed=44, constant_pH=pH,
+ search_algorithm="parallel")
@classmethod
def setUpClass(cls):
@@ -68,13 +69,13 @@ def setUpClass(cls):
@classmethod
def ideal_alpha(cls, pH):
- return 1.0 / (1 + 10**(cls.pKa - pH))
+ return 1. / (1. + 10.**(cls.pKa - pH))
def test_ideal_titration_curve(self):
- N0 = ReactionEnsembleTest.N0
- types = ReactionEnsembleTest.types
- system = ReactionEnsembleTest.system
- RE = ReactionEnsembleTest.RE
+ RE = self.RE
+ N0 = self.N0
+ types = self.types
+ system = self.system
# Set the hidden particle type to the lowest possible number to speed
# up the simulation
@@ -100,9 +101,9 @@ def test_ideal_titration_curve(self):
# note you cannot calculate the pH via -log10(/volume) in the
# constant pH ensemble, since the volume is totally arbitrary and does
# not influence the average number of protons
- pH = ReactionEnsembleTest.pH
- pKa = ReactionEnsembleTest.pKa
- target_alpha = ReactionEnsembleTest.ideal_alpha(pH)
+ pH = self.pH
+ pKa = self.pKa
+ target_alpha = Test.ideal_alpha(pH)
rel_error_alpha = abs(average_alpha - target_alpha) / target_alpha
# relative error
self.assertLess(
diff --git a/testsuite/python/get_neighbors.py b/testsuite/python/get_neighbors.py
new file mode 100644
index 00000000000..67b904d4b9b
--- /dev/null
+++ b/testsuite/python/get_neighbors.py
@@ -0,0 +1,164 @@
+# Copyright (C) 2022 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 .
+import unittest as ut
+import numpy as np
+import espressomd
+import espressomd.lees_edwards
+
+
+class Test(ut.TestCase):
+ system = espressomd.System(box_l=[1., 1., 1.])
+ system.cell_system.skin = 0.01
+ system.time_step = 0.01
+ n_nodes = system.cell_system.get_state()["n_nodes"]
+ node_grid = system.cell_system.node_grid
+
+ def setUp(self):
+ self.system.cell_system.set_regular_decomposition()
+ self.system.cell_system.node_grid = self.node_grid
+ self.system.periodicity = [True, True, True]
+
+ def tearDown(self):
+ self.system.part.clear()
+ self.system.lees_edwards.protocol = None
+
+ def get_neighbors_list(self, p1, dist):
+ result = []
+ for p2 in self.system.part.all():
+ if p1.id != p2.id and self.system.distance(p1, p2) <= dist:
+ result.append(p2.id)
+ return result
+
+ def check_neighbors(self, dist):
+ for p in self.system.part.all():
+ list_ref = np.sort(self.get_neighbors_list(p, dist))
+ list_core = np.sort(self.system.cell_system.get_neighbors(p, dist))
+ np.testing.assert_array_equal(list_core, list_ref)
+
+ @ut.skipIf(n_nodes == 1, "only runs for 2 or more MPI ranks")
+ def test_get_neighbors_different_domains(self):
+ system = self.system
+ get_neighbors = system.cell_system.get_neighbors
+ pos_red = system.box_l - 0.01
+ pos_black = np.array(3 * [0.01])
+ with self.subTest(msg='no ghost update required'):
+ # check particles on neighboring nodes
+ p1, p2 = system.part.add(id=[1, 2], pos=[pos_black, pos_red])
+ assert p1.node != p2.node
+ np.testing.assert_array_equal(get_neighbors(p1, 0.1), [2])
+ np.testing.assert_array_equal(get_neighbors(p2, 0.1), [1])
+ system.part.clear()
+ with self.subTest(msg='ghost update required'):
+ # check particles on the same node
+ p1, p2 = system.part.add(id=[1, 2], pos=[pos_black, pos_black])
+ assert p1.node == p2.node
+ np.testing.assert_array_equal(get_neighbors(p1, 0.1), [2])
+ np.testing.assert_array_equal(get_neighbors(p2, 0.1), [1])
+ # check after move to the same node
+ p2.pos = pos_black + [0.12, 0., 0.]
+ assert p1.node == p2.node
+ np.testing.assert_array_equal(get_neighbors(p1, 0.1), [])
+ np.testing.assert_array_equal(get_neighbors(p2, 0.1), [])
+ # check after move to a neighboring node (ghost update required)
+ p2.pos = pos_red
+ assert p1.node == p2.node # ghosts not updated yet
+ np.testing.assert_array_equal(get_neighbors(p1, 0.1), [2])
+ np.testing.assert_array_equal(get_neighbors(p2, 0.1), [1])
+ # check after move to a neighboring node + forced ghost update
+ system.integrator.run(0, recalc_forces=False)
+ assert p1.node != p2.node
+ np.testing.assert_array_equal(get_neighbors(p1, 0.1), [2])
+ np.testing.assert_array_equal(get_neighbors(p2, 0.1), [1])
+
+ @ut.skipIf(n_nodes == 3, "Only runs if number of MPI cores is not 3")
+ def test_get_neighbors_random_positions(self):
+ system = self.system
+ system.part.add(pos=np.random.random((20, 3)) * system.box_l)
+ for _ in range(10):
+ system.part.all().pos = np.random.random((20, 3)) * system.box_l
+ self.check_neighbors(0.2)
+ system.part.clear()
+
+ def test_n_square(self):
+ """
+ Check the N-square system is supported. The search distance is not
+ bounded by the box size. No double counting occurs.
+ """
+ system = self.system
+ system.cell_system.set_n_square()
+ p1 = system.part.add(pos=[0., 0., 0.])
+ p2 = system.part.add(pos=[0., 0., 0.49])
+ self.assertEqual(system.cell_system.get_neighbors(p1, 0.5), [p2.id])
+ self.assertEqual(system.cell_system.get_neighbors(p1, 50.), [p2.id])
+
+ @ut.skipIf(n_nodes == 1, "Only runs if number of MPI cores is 2 or more")
+ def test_domain_decomposition(self):
+ """
+ Check the neighbor search distance is bounded by the regular
+ decomposition maximal range and that aperiodic systems are supported.
+ """
+ system = self.system
+ get_neighbors = system.cell_system.get_neighbors
+ local_box_l = np.copy(system.box_l / system.cell_system.node_grid)
+ max_range = np.min(local_box_l) / 2.
+ pos = np.array([0.01, 0.01, 0.01])
+ p1, p2 = system.part.add(pos=[pos, -pos])
+ self.check_neighbors(0.9999 * max_range)
+ with np.testing.assert_raises_regex(ValueError, "pair search distance .* bigger than the decomposition range"):
+ get_neighbors(p1, 1.0001 * max_range)
+
+ # change periodicity
+ pair_dist = np.linalg.norm(system.distance_vec(p1, p2))
+ self.assertEqual(get_neighbors(p1, 1.1 * pair_dist), [p2.id])
+ system.periodicity = [False, True, True]
+ self.assertEqual(get_neighbors(p1, 1.1 * pair_dist), [])
+
+ def test_hybrid_decomposition(self):
+ """
+ Set up colloids in a N-square cell system coupled to ions on a
+ regular decomposition cell system.
+ """
+ system = self.system
+ system.cell_system.set_hybrid_decomposition(
+ n_square_types={1}, cutoff_regular=0.1)
+ p_ion = system.part.add(pos=[0., 0., 0.], type=0)
+ p_colloid = system.part.add(pos=[0., 0., 0.], type=1)
+
+ msg = "Cannot search for neighbors in the hybrid decomposition cell system"
+ with np.testing.assert_raises_regex(RuntimeError, msg):
+ system.cell_system.get_neighbors(p_ion, 0.05)
+ with np.testing.assert_raises_regex(RuntimeError, msg):
+ system.cell_system.get_neighbors(p_colloid, 0.05)
+
+ def test_lees_edwards(self):
+ """
+ Check the Lees-Edwards position offset is taken into account
+ in the distance calculation.
+ """
+ system = self.system
+ system.lees_edwards.protocol = espressomd.lees_edwards.LinearShear(
+ initial_pos_offset=0.1, time_0=0., shear_velocity=1.0)
+ system.lees_edwards.shear_direction = 0
+ system.lees_edwards.shear_plane_normal = 2
+ p1 = system.part.add(pos=[0., 0., 0.])
+ p2 = system.part.add(pos=[0., 0., 0.99])
+ self.assertEqual(system.cell_system.get_neighbors(p1, 0.10), [])
+ self.assertEqual(system.cell_system.get_neighbors(p1, 0.15), [p2.id])
+
+
+if __name__ == "__main__":
+ ut.main()
diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py
index b205e98bbd8..3daa92f4625 100644
--- a/testsuite/python/reaction_ensemble.py
+++ b/testsuite/python/reaction_ensemble.py
@@ -58,7 +58,7 @@ class ReactionEnsembleTest(ut.TestCase):
product_types = [types["A-"], types["H+"]]
product_coefficients = [1, 1]
nubar = 1
- system = espressomd.System(box_l=np.ones(3) * (N0 / c0)**(1.0 / 3.0))
+ system = espressomd.System(box_l=np.ones(3) * np.cbrt(N0 / c0))
np.random.seed(69) # make reaction code fully deterministic
system.cell_system.skin = 0.4
volume = system.volume()
@@ -68,8 +68,8 @@ class ReactionEnsembleTest(ut.TestCase):
# degree of dissociation alpha = N_A / N_HA = N_H / N_0
gamma = target_alpha**2 / (1. - target_alpha) * N0 / (volume**nubar)
RE = espressomd.reaction_methods.ReactionEnsemble(
- kT=temperature,
- exclusion_range=exclusion_range, seed=12)
+ seed=12, kT=temperature,
+ exclusion_range=exclusion_range, search_algorithm="parallel")
@classmethod
def setUpClass(cls):
diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py
index dc05cb85b68..59fa1e1335f 100644
--- a/testsuite/python/reaction_methods_interface.py
+++ b/testsuite/python/reaction_methods_interface.py
@@ -34,7 +34,7 @@ def tearDown(self):
self.system.part.clear()
def check_interface(self, method, kT, exclusion_range,
- gamma, exclusion_radius_per_type):
+ gamma, exclusion_radius_per_type, search_algorithm):
def check_reaction_parameters(reactions, parameters):
for reaction, params in zip(reactions, parameters):
for key in reaction.required_keys():
@@ -72,6 +72,7 @@ def check_reaction_parameters(reactions, parameters):
method.exclusion_range,
exclusion_range,
delta=1e-10)
+ self.assertEqual(method.search_algorithm, search_algorithm)
if not isinstance(method, espressomd.reaction_methods.WidomInsertion):
self.assertEqual(
list(method.exclusion_radius_per_type.keys()), [1])
@@ -162,18 +163,20 @@ def test_interface(self):
# reaction ensemble
method = espressomd.reaction_methods.ReactionEnsemble(
- kT=1.4, seed=12, **params)
- self.check_interface(method, kT=1.4, gamma=1.2, **params)
+ kT=1.4, seed=12, search_algorithm="order_n", **params)
+ self.check_interface(method, kT=1.4, gamma=1.2,
+ search_algorithm="order_n", **params)
# constant pH ensemble
method = espressomd.reaction_methods.ConstantpHEnsemble(
- kT=1.5, seed=14, constant_pH=10, **params)
- self.check_interface(method, kT=1.5, gamma=1.2, **params)
+ kT=1.5, seed=14, search_algorithm="parallel", constant_pH=10., **params)
+ self.check_interface(method, kT=1.5, gamma=1.2,
+ search_algorithm="parallel", **params)
# Widom insertion
method = espressomd.reaction_methods.WidomInsertion(kT=1.6, seed=16)
self.check_interface(method, kT=1.6, gamma=1., exclusion_range=0.,
- exclusion_radius_per_type={})
+ exclusion_radius_per_type={}, search_algorithm=None)
def test_exceptions(self):
single_reaction_params = {
@@ -257,6 +260,8 @@ def test_exceptions(self):
method.delete_particle(p_id=0)
with self.assertRaisesRegex(RuntimeError, "Trying to remove some non-existing particles from the system via the inverse Widom scheme"):
widom.calculate_particle_insertion_potential_energy(reaction_id=0)
+ with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"):
+ widom.search_algorithm = "order_n"
# check other exceptions
with self.assertRaisesRegex(ValueError, "Invalid value for 'volume'"):
@@ -287,6 +292,9 @@ def test_exceptions(self):
with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 1: radius -0.10"):
espressomd.reaction_methods.ReactionEnsemble(
kT=1., seed=12, exclusion_range=1., exclusion_radius_per_type={1: -0.1})
+ with self.assertRaisesRegex(ValueError, "Unknown search algorithm 'unknown'"):
+ espressomd.reaction_methods.ReactionEnsemble(
+ kT=1., seed=12, exclusion_range=1., search_algorithm="unknown")
method = espressomd.reaction_methods.ReactionEnsemble(
kT=1., exclusion_range=1., seed=12, exclusion_radius_per_type={1: 0.1})
with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 2: radius -0.10"):
|