Monte Carlo re-implementation in Python (#4765)

Description of changes:
- rewrite Monte Carlo methods in Python
- allow user-defined MC acceptance criteria
- API change: reaction methods now take a `system` argument

Co-authored-by: Pablo Miguel Blanco Andrés <blancoapa@natur.cuni.cz>

Author: jngrad

CMakeLists.txt

@@ -35,7 +35,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 # project info
 project(
   ESPResSo
-  VERSION "5.0.1"
+  VERSION "5.1.0"
   LANGUAGES C CXX
   HOMEPAGE_URL "https://espressomd.org"
   DESCRIPTION

doc/sphinx/particles.rst

@@ -734,23 +734,13 @@ Particle number counting feature
 
 
 Knowing the number of particles of a certain type in simulations where
-particle numbers can fluctuate is of interest.
-Particle ids can be stored in a map for each individual type::
+particle numbers can fluctuate is of interest::
 
     import espressomd
     system = espressomd.System(box_l=[1, 1, 1])
-    system.setup_type_map([_type])
     system.number_of_particles(_type)
 
-If you want to keep track of particle ids of a certain type you have to
-initialize the method by calling  ::
-
-    system.setup_type_map([_type])
-
-After that the system will keep track of particle ids of that type. Keeping
-track of particles of a given type is not enabled by default since it requires
-memory. The keyword ``number_of_particles`` as argument will return the number
-of particles which have the given type. For counting the number of particles
+For counting the number of particles
 of a given type you could also use
 :meth:`ParticleList.select() <espressomd.particle_data.ParticleList.select>` ::
 
@@ -763,9 +753,8 @@ of a given type you could also use
     print(len(system.part.select(type=2)))
 
 However calling ``select(type=type)`` results in looping over all particles,
-which is slow. In contrast, the system
-:meth:`~espressomd.system.System.number_of_particles` method can return the
-number of particles with that type.
+and creating temporary particle objects, which is slow. In contrast, the system
+:meth:`~espressomd.system.System.number_of_particles` method has low overhead.
 
 .. _Self-propelled swimmers:
 

doc/sphinx/reaction_methods.rst

@@ -35,12 +35,6 @@ Please keep in mind the following remarks:
   :meth:`~espressomd.reaction_methods.ReactionAlgorithm.set_non_interacting_type`
   in all reaction method classes.
 
-* Some of the functionality requires particle book-keeping. If your simulation
-  script raises runtime errors about "provided particle type X is currently not
-  tracked by the system", use :meth:`system.setup_type_map(type_list=[X])
-  <espressomd.system.System.setup_type_map>` where ``X`` is the particle
-  type to track.
-
 Thermodynamic ensembles
 -----------------------
 
@@ -196,7 +190,8 @@ As before in the reaction ensemble, one can define multiple reactions (e.g. for
 
 .. code-block:: python
 
-    cpH = reaction_methods.ConstantpHEnsemble(kT=1, exclusion_range=1., seed=77)
+    cpH = reaction_methods.ConstantpHEnsemble(
+        kT=1., seed=77, exclusion_range=1., system=system)
     cpH.add_reaction(gamma=K_diss,
                      reactant_types=[0],
                      reactant_coefficients=[1],
@@ -252,7 +247,7 @@ For this one has to provide the following reaction to the Widom method:
 .. code-block:: python
 
     type_B = 1
-    widom = reaction_methods.WidomInsertion(kT=1, seed=77)
+    widom = reaction_methods.WidomInsertion(kT=1., seed=77, system=system)
     widom.add_reaction(reactant_types=[],
                        reactant_coefficients=[],
                        product_types=[type_B],
@@ -353,14 +348,44 @@ 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.
+These integration instabilities can be avoided by defining a distance
+``exclusion_range`` around particles in the system, within which new particles
+will not be inserted (and old particles will not be removed, to maintain detailed balance).
+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 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.reaction_methods.ExclusionRadius`.
+
+.. _Writing new Monte Carlo methods:
+
+Writing new Monte Carlo methods
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Most of the logic for reaction methods is implemented at the Python level.
+The C++ core is only used for performance-relevant operations on particles.
+Hence, one can prototype new reaction methods with relative ease.
+For example, the acceptance probability for a Monte Carlo trial move
+is exposed in :meth:`ReactionAlgorithm.calculate_acceptance_probability()
+<espressomd.reaction_methods.ReactionAlgorithm.calculate_acceptance_probability>`.
+Reaction method classes override this function with their custom expression
+for the acceptance probability.

doc/tutorials/boltzmann_inversion/boltzmann_inversion.ipynb

@@ -230,7 +230,6 @@
     "system = espressomd.System(box_l=[BOX_LENGTH]*3)\n",
     "system.time_step = DT\n",
     "system.cell_system.skin = SKIN\n",
-    "system.setup_type_map(type_list=list(types.values()))\n",
     "np.random.seed(42)\n",
     "\n",
     "# Add the colloidal particles\n",

doc/tutorials/constant_pH/constant_pH.ipynb

@@ -623,6 +623,7 @@
     "# SOLUTION CELL\n",
     "exclusion_range = PARTICLE_SIZE_REDUCED if USE_WCA else 0.0\n",
     "RE = espressomd.reaction_methods.ConstantpHEnsemble(\n",
+    "    system=system,\n",
     "    kT=KT_REDUCED,\n",
     "    exclusion_range=exclusion_range,\n",
     "    seed=77,\n",

doc/tutorials/grand_canonical_monte_carlo/grand_canonical_monte_carlo.ipynb

@@ -597,7 +597,8 @@
     "        c_salt_res_sim.magnitude))\n",
     "RE = espressomd.reaction_methods.ReactionEnsemble(kT=KT,\n",
     "                                                  exclusion_range=1.0,\n",
-    "                                                  seed=reaction_seed)\n",
+    "                                                  seed=reaction_seed,\n",
+    "                                                  system=system)\n",
     "RE.add_reaction(gamma=K_XX,\n",
     "                reactant_types=[],\n",
     "                reactant_coefficients=[],\n",

doc/tutorials/widom_insertion/widom_insertion.ipynb

@@ -338,7 +338,7 @@
    "source": [
     "# SOLUTION CELL\n",
     "# add reaction for Widom particle insertion\n",
-    "widom = espressomd.reaction_methods.WidomInsertion(kT=KT, seed=42)\n",
+    "widom = espressomd.reaction_methods.WidomInsertion(kT=KT, seed=42, system=system)\n",
     "widom.add_reaction(reactant_types=[], reactant_coefficients=[],\n",
     "                   product_types=[types[\"Xplus\"], types[\"Xminus\"]], product_coefficients=[1, 1],\n",
     "                   default_charges={types[\"Xplus\"]: charges[\"Xplus\"], types[\"Xminus\"]: charges[\"Xminus\"]})\n",

maintainer/benchmarks/mc_acid_base_reservoir.py

@@ -226,6 +226,7 @@ def calc_donnan_coefficient(c_acid, I_res, charge=-1):
 # ### Set up the constant pH ensemble using the reaction ensemble module
 exclusion_range = PARTICLE_SIZE_REDUCED
 RE = espressomd.reaction_methods.ReactionEnsemble(
+    system=system,
     kT=KT_REDUCED,
     exclusion_range=exclusion_range,
     seed=77

samples/grand_canonical.py

@@ -89,13 +89,12 @@
             epsilon=wca_eps, sigma=wca_sig)
 
 RE = espressomd.reaction_methods.ReactionEnsemble(
-    kT=temperature, exclusion_range=wca_sig, seed=3)
+    system=system, kT=temperature, exclusion_range=wca_sig, seed=3)
 RE.add_reaction(
     gamma=cs_bulk**2 * np.exp(excess_chemical_potential_pair / temperature),
     reactant_types=[], reactant_coefficients=[], product_types=[1, 2],
     product_coefficients=[1, 1], default_charges={1: -1, 2: +1})
 print(RE.get_status())
-system.setup_type_map(type_list=[0, 1, 2])
 
 # Set the hidden particle type to the lowest possible number to speed
 # up the simulation

samples/reaction_ensemble_complex_reaction.py

@@ -82,7 +82,7 @@
 
 # use an exclusion radius of 0 to simulate an ideal gas
 RE = espressomd.reaction_methods.ReactionEnsemble(
-    kT=1, exclusion_range=0, seed=4)
+    kT=1., exclusion_range=0., seed=4, system=system)
 
 
 RE.add_reaction(