From eb4e3d7e75a0faadbbb9429ba052b94eaf569e47 Mon Sep 17 00:00:00 2001 From: pmblanco Date: Wed, 23 Nov 2022 13:21:12 +0100 Subject: [PATCH 1/2] bugfix: reaction_methods always check for empty ids before performing reaction steps --- .../reaction_methods/ReactionAlgorithm.cpp | 27 +++++ .../reaction_methods/ReactionAlgorithm.hpp | 1 + src/core/reaction_methods/WidomInsertion.cpp | 3 + testsuite/python/CMakeLists.txt | 1 + testsuite/python/reaction_bookkeeping.py | 103 ++++++++++++++++++ 5 files changed, 135 insertions(+) create mode 100644 testsuite/python/reaction_bookkeeping.py diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index c344fc95bd1..2f627e0358f 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -56,6 +56,8 @@ void ReactionAlgorithm::do_reaction(int reaction_steps) { // assume that the kinetic part drops out in the process of calculating // ensemble averages (kinetic part may be separated and crossed out) auto current_E_pot = mpi_calculate_potential_energy(); + // Setup the list of empty pids for bookeeping + setup_bookkeeping_of_empty_pids(); for (int i = 0; i < reaction_steps; i++) { int reaction_id = i_random(static_cast(reactions.size())); generic_oneway_reaction(*reactions[reaction_id], current_E_pot); @@ -597,4 +599,29 @@ bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, return false; } +/** + * Cleans the list of empty pids and searches for empty pid in the system + */ +void ReactionAlgorithm::setup_bookkeeping_of_empty_pids() { + + // Clean-up the list of empty pids + m_empty_p_ids_smaller_than_max_seen_particle = {}; + + // Loop from 0 to the maximum particle id in the system + + std::vector existing_ids = get_particle_ids(); + int max_id = *max_element(std::begin(existing_ids), std::end(existing_ids)); + + for (int pid = 0; pid < max_id; pid++) { + // Search if the pid is in the system + + if (std::find(existing_ids.begin(), existing_ids.end(), pid) == + existing_ids.end()) { + // if the pid is not in the system, add it to the list of empty pids + + m_empty_p_ids_smaller_than_max_seen_particle.push_back(pid); + } + } +} + } // namespace ReactionMethods diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 434b54dd5ae..20c01cdb67b 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -159,6 +159,7 @@ class ReactionAlgorithm { return {m_slab_start_z, m_slab_end_z}; } + void setup_bookkeeping_of_empty_pids(); void delete_particle(int p_id); void add_reaction(std::shared_ptr const &new_reaction); void delete_reaction(int reaction_id) { diff --git a/src/core/reaction_methods/WidomInsertion.cpp b/src/core/reaction_methods/WidomInsertion.cpp index e6ce8c15d25..4b17f036688 100644 --- a/src/core/reaction_methods/WidomInsertion.cpp +++ b/src/core/reaction_methods/WidomInsertion.cpp @@ -34,6 +34,9 @@ double WidomInsertion::calculate_particle_insertion_potential_energy( auto const E_pot_old = mpi_calculate_potential_energy(); + // Setup the list of empty pids for bookeeping + setup_bookkeeping_of_empty_pids(); + // make reaction attempt and immediately reverse it auto const change_tracker = make_reaction_attempt(current_reaction); auto const E_pot_new = mpi_calculate_potential_energy(); diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 3e8c6516448..64052a2418e 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -175,6 +175,7 @@ python_test(FILE script_interface.py MAX_NUM_PROC 4) python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 1) python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4) python_test(FILE reaction_complex.py MAX_NUM_PROC 1) +python_test(FILE reaction_bookkeeping.py MAX_NUM_PROC 1) python_test(FILE widom_insertion.py MAX_NUM_PROC 1) python_test(FILE constant_pH.py MAX_NUM_PROC 1) python_test(FILE constant_pH_stats.py MAX_NUM_PROC 4 LABELS long) diff --git a/testsuite/python/reaction_bookkeeping.py b/testsuite/python/reaction_bookkeeping.py new file mode 100644 index 00000000000..324f52add80 --- /dev/null +++ b/testsuite/python/reaction_bookkeeping.py @@ -0,0 +1,103 @@ +# +# Copyright (C) 2013-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 unittest_decorators as utx + +import espressomd +import espressomd.reaction_methods +import numpy as np + + +@utx.skipIfMissingFeatures(["ELECTROSTATICS"]) +class ReactionMethodsBookkeepingTest(ut.TestCase): + """ + Test that two different instances of the reaction methods + do not break the particle id bookkeeping. + """ + pH = 10. + pKa = 7. + exclusion_range = 1. + seed = 12345 + kT = 1. + BOX_LENGTH = 100. + N_SALT = 10 + N_acid = 10 + + types = {"H": 0, "Na": 1, "Cl": 2, "HA": 3, "A": 4} + charges = {"H": 1.0, "Na": 1.0, "Cl": -1.0, "HA": 0.0, "A": -1.0} + + system = espressomd.System(box_l=[BOX_LENGTH, ] * 3) + + cph = espressomd.reaction_methods.ConstantpHEnsemble( + constant_pH=pH, + kT=kT, + exclusion_range=exclusion_range, + seed=seed) + widom = espressomd.reaction_methods.WidomInsertion( + kT=kT, + seed=seed) + + @classmethod + def setUpClass(cls): + cls.system.part.add(type=[cls.types["Na"]] * cls.N_SALT, + pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH, + q=[cls.charges["Na"]] * cls.N_SALT) + cls.system.part.add(type=[cls.types["Cl"]] * cls.N_SALT, + pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH, + q=[cls.charges["Cl"]] * cls.N_SALT) + cls.system.part.add(type=[cls.types["HA"]] * cls.N_acid, + pos=np.random.rand(cls.N_acid, 3) * cls.BOX_LENGTH, + q=[cls.charges["HA"]] * cls.N_acid) + + cls.cph.add_reaction( + gamma=10**(-cls.pKa), + reactant_types=[cls.types["HA"]], + product_types=[cls.types["A"], cls.types["H"]], + default_charges={cls.types["HA"]: cls.charges["HA"], + cls.types["A"]: cls.charges["A"], + cls.types["H"]: cls.charges["H"]} + ) + + cls.widom.add_reaction( + reactant_types=[], + reactant_coefficients=[], + product_types=[cls.types["Na"], cls.types["Cl"]], + product_coefficients=[1, 1], + default_charges={cls.types["Na"]: cls.charges["Na"], + cls.types["Cl"]: cls.charges["Cl"]} + ) + cls.system.setup_type_map(type_list=list(cls.types.values())) + + def test_reaction_bookeeping(self): + self.widom.calculate_particle_insertion_potential_energy(reaction_id=0) + self.cph.reaction(reaction_steps=100) + + # Measure the chemical potential + for _ in range(50): + self.widom.calculate_particle_insertion_potential_energy( + reaction_id=0) + charge = 0. + for p in self.system.part: + charge += p.q + self.assertEqual(charge, 0.) + + +if __name__ == "__main__": + ut.main() From 72a029a501b48a61f3af06898004faeb18ae6d98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 12 Dec 2022 15:22:00 +0100 Subject: [PATCH 2/2] core: Use O(N) algorithm instead of O(N^2) --- .../reaction_methods/ReactionAlgorithm.cpp | 21 +++++++------------ testsuite/python/reaction_bookkeeping.py | 3 ++- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 2f627e0358f..be949fd0e26 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -603,24 +603,17 @@ bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, * Cleans the list of empty pids and searches for empty pid in the system */ void ReactionAlgorithm::setup_bookkeeping_of_empty_pids() { - // Clean-up the list of empty pids - m_empty_p_ids_smaller_than_max_seen_particle = {}; - - // Loop from 0 to the maximum particle id in the system - - std::vector existing_ids = get_particle_ids(); - int max_id = *max_element(std::begin(existing_ids), std::end(existing_ids)); - - for (int pid = 0; pid < max_id; pid++) { - // Search if the pid is in the system - - if (std::find(existing_ids.begin(), existing_ids.end(), pid) == - existing_ids.end()) { - // if the pid is not in the system, add it to the list of empty pids + m_empty_p_ids_smaller_than_max_seen_particle.clear(); + auto particle_ids = get_particle_ids(); + std::sort(particle_ids.begin(), particle_ids.end()); + auto pid1 = -1; + for (auto pid2 : particle_ids) { + for (int pid = pid1 + 1; pid < pid2; ++pid) { m_empty_p_ids_smaller_than_max_seen_particle.push_back(pid); } + pid1 = pid2; } } diff --git a/testsuite/python/reaction_bookkeeping.py b/testsuite/python/reaction_bookkeeping.py index 324f52add80..a10862aec00 100644 --- a/testsuite/python/reaction_bookkeeping.py +++ b/testsuite/python/reaction_bookkeeping.py @@ -58,7 +58,8 @@ class ReactionMethodsBookkeepingTest(ut.TestCase): def setUpClass(cls): cls.system.part.add(type=[cls.types["Na"]] * cls.N_SALT, pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH, - q=[cls.charges["Na"]] * cls.N_SALT) + q=[cls.charges["Na"]] * cls.N_SALT, + id=list(range(20, 20 + cls.N_SALT))) cls.system.part.add(type=[cls.types["Cl"]] * cls.N_SALT, pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH, q=[cls.charges["Cl"]] * cls.N_SALT)