Source code for espressomd.reaction_methods

# Copyright (C) 2010-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 <http://www.gnu.org/licenses/>.

import numpy as np
import warnings
import math
import sys
from .script_interface import ScriptInterfaceHelper, script_interface_register
from .code_features import has_features
from . import utils


[docs] @script_interface_register class SingleReaction(ScriptInterfaceHelper): _so_name = "ReactionMethods::SingleReaction" _so_creation_policy = "GLOBAL" _so_bind_methods = ("get_acceptance_rate",) def __init__(self, **kwargs): super().__init__(**kwargs) if not 'sip' in kwargs: utils.check_valid_keys(self.valid_keys(), kwargs.keys()) self.reaction_types = tuple( sorted(set(tuple(self.reactant_types) + tuple(self.product_types))))
[docs] def valid_keys(self): return self.required_keys()
[docs] def required_keys(self): return {"reactant_types", "reactant_coefficients", "gamma", "product_types", "product_coefficients"}
[docs] def make_backward_reaction(self): return SingleReaction( gamma=1. / self.gamma, reactant_types=self.product_types, reactant_coefficients=self.product_coefficients, product_types=self.reactant_types, product_coefficients=self.reactant_coefficients)
[docs] class ReactionAlgorithm(ScriptInterfaceHelper): """ This class provides the base class for Reaction Algorithms like the Reaction Ensemble algorithm and the constant pH method. Initialize the reaction algorithm by setting the standard pressure, temperature, and the exclusion range. Note: When creating particles the velocities of the new particles are set according the Maxwell-Boltzmann distribution. In this step the mass of the new particle is assumed to equal 1. Parameters ---------- kT : :obj:`float` Thermal energy of the system in simulation units exclusion_range : :obj:`float` Minimal distance from any particle, within which new particles will not be inserted. seed : :obj:`int` 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 ------- remove_constraint() Remove any previously defined constraint. Requires setting the volume using :meth:`set_volume`. set_cylindrical_constraint_in_z_direction() Constrain the reaction moves within a cylinder aligned on the z-axis. Requires setting the volume using :meth:`set_volume`. Parameters ---------- center_x : :obj:`float` x coordinate of center of the cylinder. center_y : :obj:`float` y coordinate of center of the cylinder. radius_of_cylinder : :obj:`float` radius of the cylinder set_wall_constraints_in_z_direction() Restrict the sampling area to a slab in z-direction. Requires setting the volume using :meth:`set_volume`. This constraint is necessary when working with :ref:`Electrostatic Layer Correction (ELC)`. Parameters ---------- slab_start_z : :obj:`float` z coordinate of the bottom wall. slab_end_z : :obj:`float` z coordinate of the top wall. Examples -------- >>> import espressomd >>> import espressomd.shapes >>> import espressomd.electrostatics >>> import espressomd.reaction_methods >>> import numpy as np >>> # setup a charged system >>> box_l = 20 >>> elc_gap = 10 >>> system = espressomd.System(box_l=[box_l, box_l, box_l + elc_gap]) >>> system.time_step = 0.001 >>> system.cell_system.skin = 0.4 >>> types = {"HA": 0, "A-": 1, "H+": 2, "wall": 3} >>> charges = {types["HA"]: 0, types["A-"]: -1, types["H+"]: +1} >>> for i in range(10): ... system.part.add(pos=np.random.random(3) * box_l, type=types["A-"], q=charges[types["A-"]]) ... system.part.add(pos=np.random.random(3) * box_l, type=types["H+"], q=charges[types["H+"]]) >>> for particle_type in charges.keys(): ... system.non_bonded_inter[particle_type, types["wall"]].wca.set_params(epsilon=1.0, sigma=1.0) >>> # add ELC actor >>> p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=1e-2) >>> elc = espressomd.electrostatics.ELC(actor=p3m, maxPWerror=1.0, gap_size=elc_gap) >>> system.actors.add(elc) >>> # add constant pH method >>> RE = espressomd.reaction_methods.ConstantpHEnsemble(kT=1, exclusion_range=1, seed=77) >>> RE.constant_pH = 2 >>> RE.add_reaction(gamma=0.0088, reactant_types=[types["HA"]], ... product_types=[types["A-"], types["H+"]], ... default_charges=charges) >>> # add walls for the ELC gap >>> RE.set_wall_constraints_in_z_direction(0, box_l) >>> RE.set_volume(box_l**3) >>> system.constraints.add(shape=espressomd.shapes.Wall(dist=0, normal=[0, 0, 1]), ... particle_type=types["wall"]) >>> system.constraints.add(shape=espressomd.shapes.Wall(dist=-box_l, normal=[0, 0, -1]), ... particle_type=types["wall"]) get_wall_constraints_in_z_direction() Returns the restrictions of the sampling area in z-direction. set_volume() Set the volume to be used in the acceptance probability of the reaction ensemble. This can be useful when using constraints, if the relevant volume is different from the box volume. If not used the default volume which is used, is the box volume. Parameters ---------- volume : :obj:`float` Volume of the system in simulation units get_volume() Get the volume to be used in the acceptance probability of the reaction ensemble. get_acceptance_rate_configurational_moves() Returns the acceptance rate for the configuration moves. get_acceptance_rate_reaction() Returns the acceptance rate for the given reaction. Parameters ---------- reaction_id : :obj:`int` Reaction id set_non_interacting_type() Sets the particle type for non-interacting particles. Default value: 100. This is used to temporarily hide particles during a reaction trial move, if they are to be deleted after the move is accepted. Please change this value if you intend to use the type 100 for some other particle types with interactions, or if you need improved performance, as the default value of 100 causes some overhead. Please also note that particles in the current implementation of the Reaction Ensemble are only hidden with respect to Lennard-Jones and Coulomb interactions. Hiding of other interactions, for example a magnetic, needs to be implemented in the code. Parameters ---------- type : :obj:`int` Particle type for the hidden particles get_non_interacting_type() Returns the type which is used for hiding particle displacement_mc_move_for_particles_of_type() Performs displacement Monte Carlo moves for particles of a given type. New positions of the displaced particles are chosen from the whole box with a uniform probability distribution and new velocities are sampled from the Maxwell-Boltzmann distribution. The sequence of moves is only accepted if each individual move in the sequence was accepted. Particles are sampled without replacement. Therefore, calling this method once for 10 particles is not equivalent to calling this method 10 times for 1 particle. Parameters ---------- type_mc : :obj:`int` Particle type which should be moved particle_number_to_be_changed : :obj:`int` Number of particles to move, defaults to 1. Particles are selected without replacement. Returns ------- :obj:`bool` Whether all moves were accepted. delete_particle() Deletes the particle of the given p_id and makes sure that the particle range has no holes. This function has some restrictions, as e.g. bonds are not deleted. Therefore only apply this function to simple ions. Parameters ---------- p_id : :obj:`int` Id of the particle to be deleted. change_reaction_constant() Changes the reaction constant of a given reaction (for both the forward and backward reactions). The ``reaction_id`` which is assigned to a reaction depends on the order in which :meth:`add_reaction` was called. The 0th reaction has ``reaction_id=0``, the next added reaction needs to be addressed with ``reaction_id=1``, etc. Parameters ---------- reaction_id : :obj:`int` Identifier of the reaction to modify. Will be multiplied by 2 internally! gamma : :obj:`float` New reaction constant for the forward reaction. """ _so_name = "ReactionMethods::ReactionAlgorithm" _so_creation_policy = "GLOBAL" _so_bind_methods = ("remove_constraint", "get_wall_constraints_in_z_direction", "set_wall_constraints_in_z_direction", "set_cylindrical_constraint_in_z_direction", "set_volume", "get_volume", "get_acceptance_rate_reaction", "set_non_interacting_type", "get_non_interacting_type", "displacement_mc_move_for_particles_of_type", "change_reaction_constant", "delete_particle", ) def __init__(self, **kwargs): if self._so_name == ReactionAlgorithm._so_name: raise RuntimeError( f"Base class '{self.__class__.__name__}' cannot be instantiated") if 'exclusion_radius' in kwargs: raise KeyError( 'the keyword `exclusion_radius` is obsolete. Currently, the equivalent keyword is `exclusion_range`') super().__init__(**kwargs) if not 'sip' in kwargs: utils.check_valid_keys(self.valid_keys(), kwargs.keys()) self._rebuild_reaction_cache()
[docs] def valid_keys(self): return {"kT", "exclusion_range", "seed", "exclusion_radius_per_type", "search_algorithm"}
[docs] def required_keys(self): return {"kT", "exclusion_range", "seed"}
[docs] def add_reaction(self, **kwargs): """ Sets up a reaction in the forward and backward direction. Parameters ---------- gamma : :obj:`float` Equilibrium constant :math:`\\Gamma` of the reaction in simulation units (see section :ref:`Reaction Ensemble` for its definition). reactant_types : list of :obj:`int` List of particle types of reactants in the reaction. reactant_coefficients : list of :obj:`int` List of stoichiometric coefficients of the reactants in the same order as the list of their types. product_types : list of :obj:`int` List of particle types of products in the reaction. product_coefficients : list of :obj:`int` List of stoichiometric coefficients of products of the reaction in the same order as the list of their types default_charges : :obj:`dict` A dictionary of default charges for types that occur in the provided reaction. check_for_electroneutrality : :obj:`bool` Check for electroneutrality of the given reaction. Default is ``True``. """ default_charges = kwargs.pop("default_charges") neutrality_check = kwargs.pop("check_for_electroneutrality", True) if not isinstance(default_charges, dict): raise TypeError("Argument 'default_charges' needs to be a dict") forward_reaction = SingleReaction(**kwargs) backward_reaction = forward_reaction.make_backward_reaction() if neutrality_check: self._check_charge_neutrality( type2charge=default_charges, reaction=forward_reaction) old_default_charges = self.default_charges for ptype, charge in default_charges.items(): if ptype in old_default_charges: if abs(old_default_charges[ptype] - charge) > 1e-10: raise ValueError( f"Cannot override charge on particle type {ptype} (from {old_default_charges[ptype]} to {charge})") for ptype, charge in default_charges.items(): if ptype not in old_default_charges: self.call_method( "set_charge_of_type", type=ptype, charge=charge) self.call_method("add_reaction", reaction=forward_reaction) self.call_method("add_reaction", reaction=backward_reaction) self.check_reaction_method() self._rebuild_reaction_cache()
[docs] def delete_reaction(self, **kwargs): """ Delete a reaction from the set of used reactions (the forward and backward reaction). The ``reaction_id`` which is assigned to a reaction depends on the order in which :meth:`add_reaction` was called. The 0th reaction has ``reaction_id=0``, the next added reaction needs to be addressed with ``reaction_id=1``, etc. After the deletion of a reaction subsequent reactions take the ``reaction_id`` of the deleted reaction. Parameters ---------- reaction_id : :obj:`int` Reaction id """ self.call_method("delete_reaction", **kwargs) self._rebuild_reaction_cache()
[docs] def check_reaction_method(self): if len(self.reactions) == 0: raise RuntimeError("Reaction system not initialized") # charges of all reactive types need to be known if has_features("ELECTROSTATICS"): for reaction in self.reactions: for p_type in reaction.reaction_types: if p_type not in self.default_charges: raise RuntimeError( f"Forgot to assign charge to type {p_type}")
def _check_charge_neutrality(self, type2charge, reaction): charges = np.array(list(type2charge.values())) if np.count_nonzero(charges) == 0: # all particles have zero charge # no need to check electroneutrality return # calculate net change of electrical charge for the reaction net_charge_change = 0.0 for coef, ptype in zip( reaction.reactant_coefficients, reaction.reactant_types): net_charge_change -= coef * type2charge[ptype] for coef, ptype in zip( reaction.product_coefficients, reaction.product_types): net_charge_change += coef * type2charge[ptype] min_abs_nonzero_charge = np.min( np.abs(charges[np.nonzero(charges)[0]])) if abs(net_charge_change) / min_abs_nonzero_charge > 1e-10: raise ValueError("Reaction system is not charge neutral") def _rebuild_reaction_cache(self): """ Maintain a local cache of the Python representation of the script interface reactions. This helps reducing the overhead when calculating the acceptance probability. """ self._reactions_cache = [x for x in self.reactions]
[docs] def reaction(self, steps): """ Performs randomly selected reactions. Parameters ---------- steps : :obj:`int`, optional The number of reactions to be performed at once, defaults to 1. """ self.call_method("setup_bookkeeping_of_empty_pids") E_pot = self.call_method("potential_energy") for _ in range(steps): reaction_id = self.call_method("get_random_reaction_index") E_pot = self.generic_oneway_reaction(reaction_id, E_pot)
[docs] def calculate_acceptance_probability(self, reaction_id, E_pot_diff): """ Calculate the acceptance probability of a Monte Carlo move. Parameters ---------- reaction_id : :obj:`int` Identifier of the reaction that was carried out in the move. E_pot_diff : :obj:`float` The potential energy difference for the move. Returns ------- :obj:`float` The acceptance probability. """ factorial_expr = self.call_method("calculate_factorial_expression") reaction = self._reactions_cache[reaction_id] return self.get_volume()**reaction.nu_bar * reaction.gamma * \ factorial_expr * math.exp(-E_pot_diff / self.kT)
[docs] def generic_oneway_reaction(self, reaction_id, E_pot_old): """ Carry out a generic one-way chemical reaction of the type `A+B+...+G +... --> K+...X + Z +...` You need to use `2A --> B` instead of `A+A --> B` since in the latter you assume distinctness of the particles, however both ways to describe the reaction are equivalent in the thermodynamic limit (large particle numbers). Furthermore, the order of the reactant and product types matters since particles will be replaced in that order! If there are less reactants than products, new product particles are created randomly in the box. Reactants get their type and charge changed to the corresponding type and charge of the products. If there are more reactants than products, excess reactant particles are deleted. Parameters ---------- reaction_id : :obj:`int` Identifier of the reaction to attempt. E_pot_old : :obj:`float` The current potential energy. Returns ------- E_pot_new : :obj:`float` The potential energy after the move. """ try: E_pot_new = self.call_method( "create_new_trial_state", reaction_id=reaction_id) if E_pot_new is None: return E_pot_old E_pot_diff = E_pot_new - E_pot_old bf = self.calculate_acceptance_probability(reaction_id, E_pot_diff) return self.call_method("make_reaction_mc_move_attempt", reaction_id=reaction_id, bf=bf, E_pot_new=E_pot_new, E_pot_old=E_pot_old) except BaseException as err: tb = sys.exc_info()[2] raise RuntimeError( "An exception was raised by a chemical reaction; the particle " "state tracking is no longer guaranteed to be correct! -- " f"{err}").with_traceback(tb)
def _check_reaction_index(self, reaction_index): if reaction_index < 0 or reaction_index >= len(self.reactions): raise IndexError(f"No reaction with id {reaction_index}")
[docs] def get_status(self): """ Returns the status of the reaction ensemble in a dictionary containing the used reactions, the used kT and the used exclusion radius. """ self.check_reaction_method() property_keys = {"reactant_coefficients", "reactant_types", "product_coefficients", "product_types", "gamma"} reactions_list = [{key: getattr(reaction, key) for key in property_keys} for reaction in self.reactions] return {"reactions": reactions_list, "kT": self.kT, "exclusion_range": self.exclusion_range, "exclusion_radius_per_type": self.exclusion_radius_per_type}
[docs] @script_interface_register class ReactionEnsemble(ReactionAlgorithm): """ This class implements the Reaction Ensemble. """ _so_name = "ReactionMethods::ReactionEnsemble"
[docs] @script_interface_register class ConstantpHEnsemble(ReactionAlgorithm): """ This class implements the constant pH Ensemble. When adding an acid-base reaction, the acid and base particle types are always assumed to be at index 0 of the lists passed to arguments ``reactant_types`` and ``product_types``. Attributes ---------- constant_pH : :obj:`float` Constant pH value. """ _so_name = "ReactionMethods::ConstantpHEnsemble"
[docs] def valid_keys(self): return {"kT", "exclusion_range", "seed", "constant_pH", "exclusion_radius_per_type", "search_algorithm"}
[docs] def required_keys(self): return {"kT", "exclusion_range", "seed", "constant_pH"}
[docs] def calculate_acceptance_probability(self, reaction_id, E_pot_diff): """ Calculate the acceptance probability of a Monte Carlo move. Parameters ---------- reaction_id : :obj:`int` Identifier of the reaction that was carried out in the move. E_pot_diff : :obj:`float` The potential energy difference for the move. Returns ------- :obj:`float` The acceptance probability. """ factorial_expr = self.call_method("calculate_factorial_expression") reaction = self._reactions_cache[reaction_id] ln_bf = E_pot_diff - reaction.nu_bar * self.kT * math.log(10.) * ( self.constant_pH + reaction.nu_bar * math.log10(reaction.gamma)) return factorial_expr * math.exp(-ln_bf / self.kT)
[docs] def add_reaction(self, *args, **kwargs): warn_msg = ( "arguments 'reactant_coefficients' and 'product_coefficients' " "are deprecated and are no longer necessary for the constant pH " "ensemble. They are kept for backward compatibility but might " "be deleted in future versions.") err_msg = ("All product and reactant coefficients must equal one in " "the constant pH method as implemented in ESPResSo.") warn_user = False if "reactant_coefficients" in kwargs: if kwargs["reactant_coefficients"][0] != 1: raise ValueError(err_msg) else: warn_user = True else: kwargs["reactant_coefficients"] = [1] if "product_coefficients" in kwargs: if kwargs["product_coefficients"][0] != 1 or kwargs["product_coefficients"][1] != 1: raise ValueError(err_msg) else: warn_user = True else: kwargs["product_coefficients"] = [1, 1] if warn_user: warnings.warn(warn_msg, FutureWarning) if (len(kwargs["product_types"]) != 2 or len( kwargs["reactant_types"]) != 1): raise ValueError( "The constant pH method is only implemented for reactions " "with two product types and one adduct type.") super().add_reaction(*args, **kwargs)
[docs] @script_interface_register class WidomInsertion(ReactionAlgorithm): """ This class implements the Widom insertion method in the canonical ensemble for homogeneous systems, where the excess chemical potential is not depending on the location. """ _so_name = "ReactionMethods::WidomInsertion"
[docs] def required_keys(self): return {"kT", "seed"}
[docs] def valid_keys(self): return {"kT", "seed"}
[docs] def add_reaction(self, **kwargs): kwargs['gamma'] = 1. super().add_reaction(**kwargs)
[docs] def calculate_particle_insertion_potential_energy(self, **kwargs): """ Measures the potential energy when particles are inserted in the system following the reaction provided in ``reaction_id``. Please define the insertion moves by calling the method :meth:`~ReactionAlgorithm.add_reaction` (with only product types specified). Note that although this function does not provide directly the chemical potential, it can be used to calculate it. For an example of such an application please check :file:`/samples/widom_insertion.py`. Parameters ---------- reaction_id : :obj:`int` Reaction identifier. Will be multiplied by 2 internally to skip reverse reactions, i.e. deletion reactions! Returns ------- :obj:`float` The particle insertion potential energy. """ return self.call_method( "calculate_particle_insertion_potential_energy", **kwargs)
[docs] def calculate_excess_chemical_potential(self, **kwargs): """ Given a set of samples of the particle insertion potential energy, calculates the excess chemical potential and its statistical error. Parameters ---------- particle_insertion_potential_energy_samples : array_like of :obj:`float` Samples of the particle insertion potential energy. N_blocks : :obj:`int`, optional Number of bins for binning analysis. Returns ------- mean : :obj:`float` Mean excess chemical potential. error : :obj:`float` Standard error of the mean. """ def do_block_analysis(samples, N_blocks): """ Performs a binning analysis of samples. Divides the samples in ``N_blocks`` equispaced blocks and returns the mean and its uncertainty """ size_block = int(len(samples) / N_blocks) block_list = [] for block in range(N_blocks): block_list.append( np.mean(samples[block * size_block:(block + 1) * size_block])) sample_mean = np.mean(block_list) sample_std = np.std(block_list, ddof=1) sample_uncertainty = sample_std / np.sqrt(N_blocks) return sample_mean, sample_uncertainty kT = self.kT gamma_samples = np.exp(-1.0 * np.array( kwargs["particle_insertion_potential_energy_samples"]) / kT) gamma_mean, gamma_std = do_block_analysis( samples=gamma_samples, N_blocks=kwargs.get("N_blocks", 16)) mu_ex_mean = -kT * np.log(gamma_mean) # full propagation of error mu_ex_Delta = 0.5 * kT * abs(-np.log(gamma_mean + gamma_std) - (-np.log(gamma_mean - gamma_std))) return mu_ex_mean, mu_ex_Delta