From e185d900f9a8405f0f9c3c7bab6e0cf86a06ec93 Mon Sep 17 00:00:00 2001 From: Julian Hossbach Date: Thu, 12 May 2022 11:55:33 +0200 Subject: [PATCH 1/2] core: Efficient search algorithm for pair interactions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implement parallel kernels to find particle pairs within a cutoff distance of a specific particle and to calculate the short-range non-bonded energy of a specific particle. Co-authored-by: Jean-Noël Grad Co-authored-by: Rudolf Weeber --- src/core/cell_system/CellStructure.hpp | 52 ++++++ src/core/cells.cpp | 44 +++++ src/core/cells.hpp | 11 ++ src/core/energy.cpp | 33 ++++ src/core/energy.hpp | 11 ++ src/python/espressomd/analyze.pxd | 5 + src/python/espressomd/analyze.pyx | 19 ++ src/python/espressomd/cell_system.py | 33 ++++ .../cell_system/CellSystem.hpp | 21 +++ testsuite/python/CMakeLists.txt | 2 + testsuite/python/analyze_energy.py | 12 ++ testsuite/python/get_neighbors.py | 164 ++++++++++++++++++ 12 files changed, 407 insertions(+) create mode 100644 testsuite/python/get_neighbors.py 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/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/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/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/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() From a257acf072330519ed79eaa3340f5588d990c9c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 18 May 2022 16:34:11 +0200 Subject: [PATCH 2/2] core: Parallel search for reaction methods --- .../reaction_methods/ReactionAlgorithm.cpp | 50 ++++++++++--------- .../reaction_methods/ReactionAlgorithm.hpp | 14 ++++-- src/python/espressomd/reaction_methods.py | 12 ++++- .../reaction_methods/ConstantpHEnsemble.hpp | 3 ++ .../reaction_methods/ReactionAlgorithm.hpp | 18 +++++++ .../reaction_methods/ReactionEnsemble.hpp | 3 ++ .../reaction_methods/WidomInsertion.hpp | 9 ++++ testsuite/python/constant_pH_stats.py | 29 +++++------ testsuite/python/reaction_ensemble.py | 6 +-- .../python/reaction_methods_interface.py | 20 +++++--- 10 files changed, 110 insertions(+), 54 deletions(-) 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/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/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/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/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"):