From a0359cb9dcac937bb342804d22fb85f6b38eaf03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 28 Jul 2023 23:46:08 +0200 Subject: [PATCH] core: Extract exclusion range neighbor search --- doc/sphinx/reaction_methods.rst | 38 +++-- src/core/reaction_methods/CMakeLists.txt | 6 +- .../reaction_methods/ConstantpHEnsemble.hpp | 12 +- src/core/reaction_methods/ExclusionRadius.cpp | 154 ++++++++++++++++++ src/core/reaction_methods/ExclusionRadius.hpp | 43 +++++ .../reaction_methods/ReactionAlgorithm.cpp | 64 +------- .../reaction_methods/ReactionAlgorithm.hpp | 36 +--- .../reaction_methods/ReactionEnsemble.hpp | 9 +- src/core/reaction_methods/WidomInsertion.hpp | 12 +- src/python/espressomd/reaction_methods.py | 60 ++++++- .../reaction_methods/ConstantpHEnsemble.hpp | 11 +- .../reaction_methods/ExclusionRadius.hpp | 109 +++++++++++++ .../reaction_methods/ReactionAlgorithm.cpp | 45 +---- .../reaction_methods/ReactionAlgorithm.hpp | 10 ++ .../reaction_methods/ReactionEnsemble.hpp | 9 +- .../reaction_methods/WidomInsertion.hpp | 34 ++-- .../reaction_methods/initialize.cpp | 2 + .../python/reaction_methods_interface.py | 15 +- 18 files changed, 474 insertions(+), 195 deletions(-) create mode 100644 src/core/reaction_methods/ExclusionRadius.cpp create mode 100644 src/core/reaction_methods/ExclusionRadius.hpp create mode 100644 src/script_interface/reaction_methods/ExclusionRadius.hpp diff --git a/doc/sphinx/reaction_methods.rst b/doc/sphinx/reaction_methods.rst index 5b2a1b012f8..9a0d7140eb7 100644 --- a/doc/sphinx/reaction_methods.rst +++ b/doc/sphinx/reaction_methods.rst @@ -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`. diff --git a/src/core/reaction_methods/CMakeLists.txt b/src/core/reaction_methods/CMakeLists.txt index 88a046e9b79..fe0c8971060 100644 --- a/src/core/reaction_methods/CMakeLists.txt +++ b/src/core/reaction_methods/CMakeLists.txt @@ -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) diff --git a/src/core/reaction_methods/ConstantpHEnsemble.hpp b/src/core/reaction_methods/ConstantpHEnsemble.hpp index ea526aa0613..66012572fa6 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.hpp @@ -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 +#include namespace ReactionMethods { @@ -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 &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 exclusion, + double constant_pH) + : ReactionAlgorithm(comm, seed, kT, exclusion), m_constant_pH(constant_pH) {} double m_constant_pH; }; diff --git a/src/core/reaction_methods/ExclusionRadius.cpp b/src/core/reaction_methods/ExclusionRadius.cpp new file mode 100644 index 00000000000..edd3d594ecf --- /dev/null +++ b/src/core/reaction_methods/ExclusionRadius.cpp @@ -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 . + */ + +#include "reaction_methods/ExclusionRadius.hpp" + +#include "cells.hpp" +#include "event.hpp" +#include "grid.hpp" +#include "particle_node.hpp" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +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(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 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()); + return check_exclusion_range(pid, type); +} diff --git a/src/core/reaction_methods/ExclusionRadius.hpp b/src/core/reaction_methods/ExclusionRadius.hpp new file mode 100644 index 00000000000..70013120ef9 --- /dev/null +++ b/src/core/reaction_methods/ExclusionRadius.hpp @@ -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 . + */ + +#pragma once + +#include + +#include + +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; + 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{}; +}; diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index f19f91b6d40..9e364577736 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -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 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); } /** diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index e617dbdebca..c62770cafaf 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -21,6 +21,7 @@ #include "config/config.hpp" +#include "ExclusionRadius.hpp" #include "SingleReaction.hpp" #include "Particle.hpp" @@ -50,20 +51,14 @@ class ReactionAlgorithm { boost::mpi::communicator const &m_comm; public: - ReactionAlgorithm( - boost::mpi::communicator const &comm, int seed, double kT, - double exclusion_range, - std::unordered_map const &exclusion_radius_per_type) - : m_comm{comm}, kT{kT}, exclusion_range{exclusion_range}, + ReactionAlgorithm(boost::mpi::communicator const &comm, int seed, double kT, + std::shared_ptr exclusion) + : m_comm{comm}, kT{kT}, m_exclusion{exclusion}, 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'"); } - if (exclusion_range < 0.) { - throw std::domain_error("Invalid value for 'exclusion_range'"); - } - set_exclusion_radius_per_type(exclusion_radius_per_type); update_volume(); } @@ -78,8 +73,7 @@ class ReactionAlgorithm { * infinite, therefore these configurations do not contribute * to the partition function and ensemble averages. */ - double exclusion_range; - std::unordered_map exclusion_radius_per_type; + std::shared_ptr m_exclusion; double volume; int non_interacting_type = 100; @@ -91,7 +85,6 @@ class ReactionAlgorithm { } auto get_kT() const { return kT; } - auto get_exclusion_range() const { return exclusion_range; } auto get_volume() const { return volume; } void set_volume(double new_volume) { if (new_volume <= 0.) { @@ -100,23 +93,6 @@ class ReactionAlgorithm { volume = new_volume; } 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; - if (exclusion_radius < 0.) { - throw std::domain_error("Invalid excluded_radius value for type " + - 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; - } void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; } void set_cyl_constraint(double center_x, double center_y, double radius); @@ -136,7 +112,6 @@ class ReactionAlgorithm { } 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; @@ -261,7 +236,6 @@ 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.; Particle *get_real_particle(int p_id) const; Particle *get_local_particle(int p_id) const; diff --git a/src/core/reaction_methods/ReactionEnsemble.hpp b/src/core/reaction_methods/ReactionEnsemble.hpp index f6ef81dd7c1..50d430fa919 100644 --- a/src/core/reaction_methods/ReactionEnsemble.hpp +++ b/src/core/reaction_methods/ReactionEnsemble.hpp @@ -21,8 +21,6 @@ #include "reaction_methods/ReactionAlgorithm.hpp" -#include - namespace ReactionMethods { /** Reaction ensemble method. @@ -36,12 +34,7 @@ namespace ReactionMethods { */ class ReactionEnsemble : public ReactionAlgorithm { public: - ReactionEnsemble( - boost::mpi::communicator const &comm, int seed, double kT, - double exclusion_radius, - const std::unordered_map &exclusion_radius_per_type) - : ReactionAlgorithm(comm, seed, kT, exclusion_radius, - exclusion_radius_per_type) {} + using ReactionAlgorithm::ReactionAlgorithm; }; } // namespace ReactionMethods diff --git a/src/core/reaction_methods/WidomInsertion.hpp b/src/core/reaction_methods/WidomInsertion.hpp index d176b84a12b..bf9ae4dda1f 100644 --- a/src/core/reaction_methods/WidomInsertion.hpp +++ b/src/core/reaction_methods/WidomInsertion.hpp @@ -19,22 +19,16 @@ #ifndef REACTION_METHODS_WIDOM_INSERTION_HPP #define REACTION_METHODS_WIDOM_INSERTION_HPP -#include "ReactionAlgorithm.hpp" +#include "reaction_methods/ReactionAlgorithm.hpp" -#include -#include +#include namespace ReactionMethods { /** Widom insertion method */ class WidomInsertion : public ReactionAlgorithm { public: - WidomInsertion( - boost::mpi::communicator const &comm, int seed, double kT, - double exclusion_range, - const std::unordered_map &exclusion_radius_per_type) - : ReactionAlgorithm(comm, seed, kT, exclusion_range, - exclusion_radius_per_type) {} + using ReactionAlgorithm::ReactionAlgorithm; double calculate_particle_insertion_potential_energy(int reaction_id) { auto &reaction = *reactions[reaction_id]; diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index 7cf5eb9ca17..768b1a2072e 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -52,6 +52,64 @@ def make_backward_reaction(self): product_coefficients=self.reactant_coefficients) +@script_interface_register +class ExclusionRadius(ScriptInterfaceHelper): + """ + Neighbor search algorithm that detects when a particle enters the exclusion + zone of another particle. The exclusion radii are particle type-dependent. + + During the neighbor search, the following cases can arise: + + * the central particle per-type exclusion radius is zero: return ``False`` + * the neighbor particle per-type exclusion radius is zero: return ``False`` + * the central and neighbor particles per-type exclusion radii are non-zero: + return ``True`` if the inter-particle distance is smaller than the sum of + their respective exclusion radii, ``False`` otherwise + * either the central particle type or the neighbor particle type is not in + ``exclusion_radius_per_type``: return ``True`` if the inter-particle + distance is smaller than ``exclusion_range``, ``False`` otherwise + + Attributes + ---------- + exclusion_radius_per_type : :obj:`dict`, optional + Mapping of particle types to exclusion radii. + exclusion_range : :obj:`float` + Minimal distance from any particle whose type + is not in ``exclusion_radius_per_type``. + search_algorithm : :obj:`str` + Pair search algorithm. Default is ``"order_n"``, which evaluates the + distance between the queried particle and all other particles in the + system, and 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 + ------- + check_exclusion_range() + Check the neighborhood of a central particle and detect if any neighbor + is too close. + + Parameters + ----------- + pid : :obj:`int` + Particle id. + ptype : :obj:`int`, optional + Particle type. If not provided, it will be read from the particle + and communicated to all MPI ranks. + + Returns + ------- + :obj:`bool` : + Whether the particle is within the exclusion radius + of another particle. + + """ + _so_name = "ReactionMethods::ExclusionRadius" + _so_creation_policy = "GLOBAL" + _so_bind_methods = ("check_exclusion_range",) + + class ReactionAlgorithm(ScriptInterfaceHelper): """ @@ -280,7 +338,7 @@ def __init__(self, **kwargs): if 'exclusion_radius' in kwargs: raise KeyError( 'the keyword `exclusion_radius` is obsolete. Currently, the equivalent keyword is `exclusion_range`') - super().__init__(**kwargs) + super().__init__(exclusion=ExclusionRadius(**kwargs), **kwargs) if not 'sip' in kwargs: utils.check_valid_keys(self.valid_keys(), kwargs.keys()) self._rebuild_reaction_cache() diff --git a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp index 4156a98ef89..77e91f9c9fa 100644 --- a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp +++ b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp @@ -55,18 +55,13 @@ class ConstantpHEnsemble : public ReactionAlgorithm { } void do_construct(VariantMap const ¶ms) override { + setup_neighbor_search(params); context()->parallel_try_catch([&]() { m_re = std::make_shared<::ReactionMethods::ConstantpHEnsemble>( context()->get_comm(), get_value(params, "seed"), - get_value(params, "kT"), - get_value(params, "exclusion_range"), - get_value(params, "constant_pH"), - get_value_or>( - params, "exclusion_radius_per_type", {})); + get_value(params, "kT"), m_exclusion->get_instance(), + get_value(params, "constant_pH")); }); - do_set_parameter("search_algorithm", - Variant{get_value_or( - params, "search_algorithm", "order_n")}); } protected: diff --git a/src/script_interface/reaction_methods/ExclusionRadius.hpp b/src/script_interface/reaction_methods/ExclusionRadius.hpp new file mode 100644 index 00000000000..616f17c50e3 --- /dev/null +++ b/src/script_interface/reaction_methods/ExclusionRadius.hpp @@ -0,0 +1,109 @@ +/* + * Copyright (C) 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 . + */ + +#pragma once + +#include "script_interface/ScriptInterface.hpp" +#include "script_interface/auto_parameters/AutoParameters.hpp" + +#include "core/reaction_methods/ExclusionRadius.hpp" + +#include +#include +#include + +namespace ScriptInterface { +namespace ReactionMethods { + +class ExclusionRadius : public AutoParameters { + std::shared_ptr<::ExclusionRadius> m_obj; + +public: + ExclusionRadius() { + add_parameters({{"search_algorithm", + [this](Variant const &v) { + context()->parallel_try_catch([&]() { + auto const key = get_value(v); + if (key == "order_n") { + m_obj->neighbor_search_order_n = true; + } else if (key == "parallel") { + m_obj->neighbor_search_order_n = false; + } else { + throw std::invalid_argument( + "Unknown search algorithm '" + key + "'"); + } + }); + }, + [this]() { + if (m_obj->neighbor_search_order_n) { + return std::string("order_n"); + } + return std::string("parallel"); + }}, + {"exclusion_range", + [this](Variant const &v) { + context()->parallel_try_catch([&]() { + m_obj->set_exclusion_range(get_value(v)); + }); + }, + [this]() { return m_obj->exclusion_range; }}, + {"exclusion_radius_per_type", + [this](Variant const &v) { + context()->parallel_try_catch([&]() { + m_obj->set_exclusion_radius_per_type( + get_value<::ExclusionRadius::map_type>(v)); + }); + }, + [this]() { + return make_unordered_map_of_variants( + m_obj->exclusion_radius_per_type); + }}}); + } + + void do_construct(VariantMap const ¶ms) override { + context()->parallel_try_catch([&]() { + m_obj = std::make_shared<::ExclusionRadius>(context()->get_comm()); + }); + if (params.count("exclusion_range")) { + do_set_parameter("exclusion_range", params.at("exclusion_range")); + } + if (params.count("exclusion_radius_per_type")) { + do_set_parameter("exclusion_radius_per_type", + params.at("exclusion_radius_per_type")); + } + } + + Variant do_call_method(std::string const &name, + VariantMap const ¶ms) override { + if (name == "check_exclusion_range") { + auto const pid = get_value(params, "pid"); + if (params.count("ptype")) { + auto const ptype = get_value(params, "ptype"); + return m_obj->check_exclusion_range(pid, ptype); + } + return m_obj->check_exclusion_range(pid); + } + return {}; + } + + auto get_instance() const { return m_obj; } +}; + +} // namespace ReactionMethods +} // namespace ScriptInterface diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp index d11c875a4bf..7723624e27b 100644 --- a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp @@ -60,24 +60,9 @@ ReactionAlgorithm::ReactionAlgorithm() { {"kT", AutoParameter::read_only, [this]() { return RE()->get_kT(); }}, {"search_algorithm", [this](Variant const &v) { - context()->parallel_try_catch([&]() { - 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 + - "'"); - } - }); + m_exclusion->do_set_parameter("search_algorithm", v); }, - [this]() { - if (RE()->neighbor_search_order_n) { - return std::string("order_n"); - } - return std::string("parallel"); - }}, + [this]() { return m_exclusion->get_parameter("search_algorithm"); }}, {"particle_inside_exclusion_range_touched", [this](Variant const &v) { RE()->particle_inside_exclusion_range_touched = get_value(v); @@ -87,32 +72,20 @@ ReactionAlgorithm::ReactionAlgorithm() { [this]() { return make_unordered_map_of_variants(RE()->charges_of_types); }}, - {"exclusion_range", AutoParameter::read_only, - [this]() { return RE()->get_exclusion_range(); }}, + {"exclusion_range", + [this](Variant const &v) { + m_exclusion->do_set_parameter("exclusion_range", v); + }, + [this]() { return m_exclusion->get_parameter("exclusion_range"); }}, {"exclusion_radius_per_type", [this](Variant const &v) { - context()->parallel_try_catch([&]() { - RE()->set_exclusion_radius_per_type( - get_value>(v)); - }); + m_exclusion->do_set_parameter("exclusion_radius_per_type", v); }, [this]() { - return make_unordered_map_of_variants( - RE()->exclusion_radius_per_type); + return m_exclusion->get_parameter("exclusion_radius_per_type"); }}}); } -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(ptr != nullptr), - std::plus<>()) == 1); - return ptr; -} - Variant ReactionAlgorithm::do_call_method(std::string const &name, VariantMap const ¶ms) { if (name == "calculate_factorial_expression") { diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp index ae3fa822fa3..f21e3a2bde2 100644 --- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp +++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp @@ -20,6 +20,7 @@ #ifndef SCRIPT_INTERFACE_REACTION_METHODS_REACTION_ALGORITHM_HPP #define SCRIPT_INTERFACE_REACTION_METHODS_REACTION_ALGORITHM_HPP +#include "ExclusionRadius.hpp" #include "SingleReaction.hpp" #include "script_interface/ScriptInterface.hpp" @@ -41,6 +42,7 @@ class ReactionAlgorithm : public AutoParameters { protected: /** Keep track of the script interface pointer of each reaction. */ std::vector> m_reactions; + std::shared_ptr m_exclusion; /** * Check reaction id is within the reaction container bounds. * Since each reaction has a corresponding backward reaction, @@ -58,6 +60,14 @@ class ReactionAlgorithm : public AutoParameters { return index; } + void setup_neighbor_search(VariantMap const ¶ms) { + auto so_ptr = get_value(params, "exclusion"); + m_exclusion = std::dynamic_pointer_cast(so_ptr); + m_exclusion->do_set_parameter("search_algorithm", + Variant{get_value_or( + params, "search_algorithm", "order_n")}); + } + public: virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() = 0; virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const diff --git a/src/script_interface/reaction_methods/ReactionEnsemble.hpp b/src/script_interface/reaction_methods/ReactionEnsemble.hpp index 6d2ea320dd2..82d4b4524ec 100644 --- a/src/script_interface/reaction_methods/ReactionEnsemble.hpp +++ b/src/script_interface/reaction_methods/ReactionEnsemble.hpp @@ -44,17 +44,12 @@ class ReactionEnsemble : public ReactionAlgorithm { } void do_construct(VariantMap const ¶ms) override { + setup_neighbor_search(params); context()->parallel_try_catch([&]() { m_re = std::make_shared<::ReactionMethods::ReactionEnsemble>( context()->get_comm(), get_value(params, "seed"), - get_value(params, "kT"), - get_value(params, "exclusion_range"), - get_value_or>( - params, "exclusion_radius_per_type", {})); + get_value(params, "kT"), m_exclusion->get_instance()); }); - 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 bd660f64187..d03d11e6b07 100644 --- a/src/script_interface/reaction_methods/WidomInsertion.hpp +++ b/src/script_interface/reaction_methods/WidomInsertion.hpp @@ -45,22 +45,31 @@ class WidomInsertion : public ReactionAlgorithm { } WidomInsertion() { - add_parameters({{"search_algorithm", - [this](Variant const &) { - if (context()->is_head_node()) { - throw std::runtime_error( - "No search algorithm for WidomInsertion"); - } - }, - []() { return none; }}}); + add_parameters( + {{"search_algorithm", + [this](Variant const &) { throw_on_exclusion_change(); }, + []() { return none; }}, + {"exclusion_range", + [this](Variant const &) { throw_on_exclusion_change(); }, + [this]() { return m_exclusion->get_parameter("exclusion_range"); }}, + {"exclusion_radius_per_type", + [this](Variant const &) { throw_on_exclusion_change(); }, + [this]() { + return m_exclusion->get_parameter("exclusion_radius_per_type"); + }}}); } void do_construct(VariantMap const ¶ms) override { + setup_neighbor_search(params); context()->parallel_try_catch([&]() { + auto exclusion = m_exclusion->get_instance(); + if (exclusion->exclusion_range != 0. or + not exclusion->exclusion_radius_per_type.empty()) { + throw std::runtime_error("No search algorithm for WidomInsertion"); + } m_re = std::make_shared<::ReactionMethods::WidomInsertion>( context()->get_comm(), get_value(params, "seed"), - get_value(params, "kT"), 0., - std::unordered_map{}); + get_value(params, "kT"), m_exclusion->get_instance()); }); } @@ -80,6 +89,11 @@ class WidomInsertion : public ReactionAlgorithm { private: std::shared_ptr<::ReactionMethods::WidomInsertion> m_re; + void throw_on_exclusion_change() const { + if (context()->is_head_node()) { + throw std::runtime_error("No search algorithm for WidomInsertion"); + } + } }; } /* namespace ReactionMethods */ diff --git a/src/script_interface/reaction_methods/initialize.cpp b/src/script_interface/reaction_methods/initialize.cpp index dcf692d228b..ad289de5ad2 100644 --- a/src/script_interface/reaction_methods/initialize.cpp +++ b/src/script_interface/reaction_methods/initialize.cpp @@ -19,6 +19,7 @@ #include "initialize.hpp" +#include "ExclusionRadius.hpp" #include "SingleReaction.hpp" #include "ConstantpHEnsemble.hpp" @@ -34,6 +35,7 @@ void initialize(Utils::Factory *om) { om->register_new("ReactionMethods::WidomInsertion"); om->register_new("ReactionMethods::ReactionEnsemble"); om->register_new("ReactionMethods::ConstantpHEnsemble"); + om->register_new("ReactionMethods::ExclusionRadius"); } } // namespace ReactionMethods } // namespace ScriptInterface diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py index 8f2bfcf4276..29e5409b2f1 100644 --- a/testsuite/python/reaction_methods_interface.py +++ b/testsuite/python/reaction_methods_interface.py @@ -94,6 +94,9 @@ def count_by_type(types): self.assertAlmostEqual( method.exclusion_radius_per_type[2], 0.2, delta=1e-10) exclusion_radius_per_type = {2: 0.2} + method.exclusion_range = 0.8 + self.assertAlmostEqual(method.exclusion_range, 0.8, delta=1e-10) + method.exclusion_range = 0.8 self.assertAlmostEqual( method.get_volume(), self.system.volume(), delta=1e-10) method.set_volume(volume=1.) @@ -229,7 +232,7 @@ def test_exceptions(self): } widom = espressomd.reaction_methods.WidomInsertion(kT=1., seed=12) method = espressomd.reaction_methods.ReactionEnsemble( - kT=1.5, exclusion_range=0.8, seed=12, exclusion_radius_per_type={1: 0.1}) + kT=1.5, exclusion_range=0.1, seed=12, exclusion_radius_per_type={1: 0.1}) method.add_reaction(**reaction_params) widom.add_reaction(**reaction_params) @@ -310,6 +313,10 @@ def test_exceptions(self): widom.calculate_particle_insertion_potential_energy(reaction_id=0) with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"): widom.search_algorithm = "order_n" + with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"): + widom.exclusion_range = 1. + with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"): + widom.exclusion_radius_per_type = {1: 2.} # check other exceptions with self.assertRaisesRegex(ValueError, "Invalid value for 'volume'"): @@ -346,10 +353,10 @@ def test_exceptions(self): method.set_non_interacting_type(type=-1) # check invalid exclusion ranges and radii - with self.assertRaisesRegex(ValueError, "Invalid value for 'exclusion_range'"): + with self.assertRaisesRegex(ValueError, "Invalid value for exclusion range"): espressomd.reaction_methods.ReactionEnsemble( kT=1., seed=12, exclusion_range=-1.) - with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 1: radius -0.10"): + with self.assertRaisesRegex(ValueError, "Invalid exclusion radius 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'"): @@ -357,7 +364,7 @@ def test_exceptions(self): 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"): + with self.assertRaisesRegex(ValueError, "Invalid exclusion radius for type 2: radius -0.10"): method.exclusion_radius_per_type = {2: -0.1} self.assertEqual(list(method.exclusion_radius_per_type.keys()), [1])