diff --git a/doc/tutorials/constant_pH/constant_pH.ipynb b/doc/tutorials/constant_pH/constant_pH.ipynb index e8bba2a0c8c..16541c6faa0 100644 --- a/doc/tutorials/constant_pH/constant_pH.ipynb +++ b/doc/tutorials/constant_pH/constant_pH.ipynb @@ -488,7 +488,7 @@ "```python\n", "exclusion_radius = 1.0 if USE_WCA else 0.0\n", "RE = espressomd.reaction_ensemble.ConstantpHEnsemble(\n", - " temperature=KT.to(\"sim_energy\").magnitude,\n", + " kT=KT.to(\"sim_energy\").magnitude,\n", " exclusion_radius=exclusion_radius,\n", " seed=77\n", ")\n", diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 28f515d5c20..35ca8290f32 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -89,7 +89,7 @@ epsilon=wca_eps, sigma=wca_sig) RE = espressomd.reaction_ensemble.ReactionEnsemble( - temperature=temperature, exclusion_radius=wca_sig, seed=3) + kT=temperature, exclusion_radius=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], diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index bdab9dbe625..cd0002c1933 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -72,19 +72,23 @@ RE = None if args.mode == "reaction_ensemble": RE = espressomd.reaction_ensemble.ReactionEnsemble( - temperature=1, + kT=1, exclusion_radius=1, seed=77) + RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], + product_types=[1, 2], product_coefficients=[1, 1], + default_charges={0: 0, 1: -1, 2: +1}) elif args.mode == "constant_pH_ensemble": RE = espressomd.reaction_ensemble.ConstantpHEnsemble( - temperature=1, exclusion_radius=1, seed=77) + kT=1, exclusion_radius=1, seed=77) RE.constant_pH = 2 + RE.add_reaction(gamma=K_diss, reactant_types=[0], + product_types=[1, 2], + default_charges={0: 0, 1: -1, 2: +1}) else: raise RuntimeError( "Please provide either --reaction_ensemble or --constant_pH_ensemble as argument ") -RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], - product_types=[1, 2], product_coefficients=[1, 1], - default_charges={0: 0, 1: -1, 2: +1}) + print(RE.get_status()) system.setup_type_map([0, 1, 2]) diff --git a/samples/reaction_ensemble_complex_reaction.py b/samples/reaction_ensemble_complex_reaction.py index 6c2bc95edb1..f8f78c9f70f 100644 --- a/samples/reaction_ensemble_complex_reaction.py +++ b/samples/reaction_ensemble_complex_reaction.py @@ -82,7 +82,7 @@ # use an exclusion radius of 0 to simulate an ideal gas RE = espressomd.reaction_ensemble.ReactionEnsemble( - temperature=1, exclusion_radius=0, seed=4) + kT=1, exclusion_radius=0, seed=4) RE.add_reaction( diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 463e30eac1c..86b42e272ea 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -103,7 +103,7 @@ system.thermostat.set_langevin(kT=temperature, gamma=1.0, seed=42) widom = espressomd.reaction_ensemble.WidomInsertion( - temperature=temperature, seed=77) + kT=temperature, seed=77) # add insertion reaction insertion_reaction_id = 0 diff --git a/src/core/reaction_methods/ConstantpHEnsemble.cpp b/src/core/reaction_methods/ConstantpHEnsemble.cpp index d071ae7a1bb..3741a6c2e7d 100644 --- a/src/core/reaction_methods/ConstantpHEnsemble.cpp +++ b/src/core/reaction_methods/ConstantpHEnsemble.cpp @@ -36,7 +36,7 @@ namespace ReactionMethods { double ConstantpHEnsemble::calculate_acceptance_probability( SingleReaction const ¤t_reaction, double E_pot_old, double E_pot_new, std::map const &old_particle_numbers) const { - auto const beta = 1.0 / temperature; + auto const beta = 1.0 / kT; auto const pKa = -current_reaction.nu_bar * log10(current_reaction.gamma); auto const ln_bf = (E_pot_new - E_pot_old) - current_reaction.nu_bar / beta * log(10) * diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index af78b2b7c2d..a395ffc8ea5 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -85,13 +85,13 @@ void ReactionAlgorithm::check_reaction_method() const { throw std::runtime_error("Reaction system not initialized"); } - if (temperature < 0) { - throw std::runtime_error("Temperatures cannot be negative. Please provide " - "a temperature (in k_B T) to the simulation. " - "Normally it should be 1.0. This will be used " - "directly to calculate beta:=1/(k_B T) which " + if (kT < 0) { + throw std::runtime_error("kT cannot be negative," + "normally it should be 1.0. This will be used" + "directly to calculate beta:=1/(k_B T) which" "occurs in the exp(-beta*E)\n"); } + #ifdef ELECTROSTATICS // check for the existence of default charges for all types that take part in // the reactions @@ -336,9 +336,9 @@ void ReactionAlgorithm::generic_oneway_reaction( auto const bf = calculate_acceptance_probability( current_reaction, E_pot_old, E_pot_new, old_particle_numbers); - std::vector exponential = { - exp(-1.0 / temperature * (E_pot_new - E_pot_old))}; - current_reaction.accumulator_exponentials(exponential); + std::vector exponential = {exp(-1.0 / kT * (E_pot_new - E_pot_old))}; + current_reaction.accumulator_potential_energy_difference_exponential( + exponential); if (m_uniform_real_distribution(m_generator) < bf) { // accept @@ -500,7 +500,7 @@ int ReactionAlgorithm::create_particle(int desired_type) { } // we use mass=1 for all particles, think about adapting this - move_particle(p_id, get_random_position_in_box(), std::sqrt(temperature)); + move_particle(p_id, get_random_position_in_box(), std::sqrt(kT)); set_particle_type(p_id, desired_type); #ifdef ELECTROSTATICS set_particle_q(p_id, charges_of_types[desired_type]); @@ -540,7 +540,7 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) { auto const &p = get_particle_data(p_id); old_positions.emplace_back(std::pair{p_id, p.r.p}); // write new position - auto const prefactor = std::sqrt(temperature / p.p.mass); + auto const prefactor = std::sqrt(kT / p.p.mass); auto const new_pos = get_random_position_in_box(); move_particle(p_id, new_pos, prefactor); check_exclusion_radius(p_id); @@ -573,7 +573,7 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type( ? std::numeric_limits::max() : calculate_current_potential_energy_of_system(); - double beta = 1.0 / temperature; + double beta = 1.0 / kT; // Metropolis algorithm since proposal density is symmetric auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old))); diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp index 6d04654d83e..cc84c25381e 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.hpp +++ b/src/core/reaction_methods/ReactionAlgorithm.hpp @@ -53,7 +53,7 @@ class ReactionAlgorithm { std::vector reactions; std::map charges_of_types; - double temperature = -10.0; + double kT = -10.0; /** * Hard sphere radius. If particles are closer than this value, * it is assumed that their interaction energy gets approximately diff --git a/src/core/reaction_methods/ReactionEnsemble.cpp b/src/core/reaction_methods/ReactionEnsemble.cpp index 09400bb1f84..0f66934b711 100644 --- a/src/core/reaction_methods/ReactionEnsemble.cpp +++ b/src/core/reaction_methods/ReactionEnsemble.cpp @@ -36,7 +36,7 @@ double ReactionEnsemble::calculate_acceptance_probability( const double factorial_expr = calculate_factorial_expression(current_reaction, old_particle_numbers); - const double beta = 1.0 / temperature; + const double beta = 1.0 / kT; // calculate Boltzmann factor return std::pow(volume, current_reaction.nu_bar) * current_reaction.gamma * factorial_expr * exp(-beta * (E_pot_new - E_pot_old)); diff --git a/src/core/reaction_methods/SingleReaction.hpp b/src/core/reaction_methods/SingleReaction.hpp index f2ad8c118f7..ffc89835db2 100644 --- a/src/core/reaction_methods/SingleReaction.hpp +++ b/src/core/reaction_methods/SingleReaction.hpp @@ -51,7 +51,8 @@ struct SingleReaction { double gamma = {}; // calculated values that are stored for performance reasons int nu_bar = {}; ///< change in particle numbers for the reaction - Utils::Accumulator accumulator_exponentials = Utils::Accumulator(1); + Utils::Accumulator accumulator_potential_energy_difference_exponential = + Utils::Accumulator(1); int tried_moves = 0; int accepted_moves = 0; double get_acceptance_rate() const { diff --git a/src/core/reaction_methods/WidomInsertion.cpp b/src/core/reaction_methods/WidomInsertion.cpp index ff10cdf50ca..113c5e79c65 100644 --- a/src/core/reaction_methods/WidomInsertion.cpp +++ b/src/core/reaction_methods/WidomInsertion.cpp @@ -62,17 +62,15 @@ std::pair WidomInsertion::measure_excess_chemical_potential( restore_properties(hidden_particles_properties, number_of_saved_properties); // 3) restore previously changed reactant particles restore_properties(changed_particles_properties, number_of_saved_properties); - std::vector exponential = { - exp(-1.0 / temperature * (E_pot_new - E_pot_old))}; - current_reaction.accumulator_exponentials(exponential); + std::vector exponential = {exp(-1.0 / kT * (E_pot_new - E_pot_old))}; + current_reaction.accumulator_potential_energy_difference_exponential( + exponential); // calculate mean excess chemical potential and standard error of the mean - std::pair result = std::make_pair( - -temperature * log(current_reaction.accumulator_exponentials.mean()[0]), - std::abs(-temperature / - current_reaction.accumulator_exponentials.mean()[0] * - current_reaction.accumulator_exponentials.std_error()[0])); - return result; + auto const &accumulator = + current_reaction.accumulator_potential_energy_difference_exponential; + return {-kT * log(accumulator.mean()[0]), + std::abs(-kT / accumulator.mean()[0] * accumulator.std_error()[0])}; } } // namespace ReactionMethods diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp index 09813744ce7..d110205fcaa 100644 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp @@ -51,7 +51,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { constexpr double tol = 100 * std::numeric_limits::epsilon(); ConstantpHEnsembleTest r_algo(42); - r_algo.temperature = 20.; + r_algo.kT = 20.; r_algo.m_constant_pH = 1.; // exception if no reaction was added @@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) { // acceptance = f_expr * exp(- E / T + nu_bar * log(10) * (pH - nu_bar * // pKa)) auto const acceptance_ref = - f_expr * std::exp(energy / r_algo.temperature + + f_expr * std::exp(energy / r_algo.kT + std::log(10.) * (r_algo.m_constant_pH + std::log10(reaction.gamma))); auto const acceptance = r_algo.calculate_acceptance_probability( diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index 537c2eec4c8..640f72dd6df 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -107,11 +107,11 @@ BOOST_AUTO_TEST_CASE(ReactionAlgorithm_test) { BOOST_CHECK_EQUAL(probability, -10.); } - // exception if temperature is negative + // exception if kT is negative BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); - // set temperature - r_algo.temperature = 1.; + // set kT + r_algo.kT = 1.; #ifdef ELECTROSTATICS // exception if reactant types have no charge information diff --git a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp index 6b78d0e023c..2d761f8642a 100644 --- a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp @@ -52,7 +52,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { ReactionEnsembleTest r_algo(42); r_algo.volume = 10.; - r_algo.temperature = 20.; + r_algo.kT = 20.; // exception if no reaction was added BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error); @@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(ReactionEnsemble_test) { // acceptance = V^{nu_bar} * gamma * f_expr * exp(- E / T) auto const acceptance_ref = std::pow(r_algo.volume, reaction.nu_bar) * reaction.gamma * f_expr * - std::exp(energy / r_algo.temperature); + std::exp(energy / r_algo.kT); auto const acceptance = r_algo.calculate_acceptance_probability( reaction, energy, 0., p_numbers); BOOST_CHECK_CLOSE(acceptance, acceptance_ref, 5 * tol); diff --git a/src/python/espressomd/reaction_ensemble.pxd b/src/python/espressomd/reaction_ensemble.pxd index 0afb10082e2..48c701c2e16 100644 --- a/src/python/espressomd/reaction_ensemble.pxd +++ b/src/python/espressomd/reaction_ensemble.pxd @@ -48,7 +48,7 @@ cdef extern from "reaction_methods/ReactionAlgorithm.hpp" namespace "ReactionMet vector[SingleReaction] reactions map[int, double] charges_of_types - double temperature + double kT double exclusion_radius double volume bool box_is_cylindric_around_z_axis diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx index 30f545da207..045f4684c9a 100644 --- a/src/python/espressomd/reaction_ensemble.pyx +++ b/src/python/espressomd/reaction_ensemble.pyx @@ -19,6 +19,7 @@ from libcpp.vector cimport vector from libcpp.memory cimport unique_ptr from cython.operator cimport dereference as deref import numpy as np +import warnings cdef class ReactionAlgorithm: @@ -36,8 +37,8 @@ cdef class ReactionAlgorithm: Parameters ---------- - temperature : :obj:`float` - The temperature at which the reaction is performed. + kT : :obj:`float` + Thermal energy of the system in simulation units exclusion_radius : :obj:`float` Minimal distance from any particle, within which new particle will not be inserted. This is useful to avoid integrator failures if particles @@ -53,13 +54,13 @@ cdef class ReactionAlgorithm: cdef CReactionAlgorithm * RE def _valid_keys(self): - return "temperature", "exclusion_radius", "seed" + return "kT", "exclusion_radius", "seed" def _required_keys(self): - return "temperature", "exclusion_radius", "seed" + return "kT", "exclusion_radius", "seed" def _set_params_in_es_core(self): - deref(self.RE).temperature = self._params["temperature"] + deref(self.RE).kT = self._params["kT"] # setting a volume is a side effect, sets the default volume of the # reaction ensemble as the volume of the cuboid simulation box. this # volume can be altered by the command "reaction ensemble volume @@ -191,8 +192,9 @@ cdef class ReactionAlgorithm: self._params["check_for_electroneutrality"] = True for k in self._required_keys_add(): if k not in kwargs: - raise ValueError("At least the following keys have to be given as keyword arguments: " + - self._required_keys_add().__str__() + " got " + kwargs.__str__()) + raise ValueError( + f"At least the following keys have to be given as keyword " + f"arguments: {self._required_keys()}, got {kwargs}") self._params[k] = kwargs[k] for k in self._valid_keys_add(): @@ -203,7 +205,9 @@ cdef class ReactionAlgorithm: self._set_params_in_es_core_add() def _valid_keys_add(self): - return "gamma", "reactant_types", "reactant_coefficients", "product_types", "product_coefficients", "default_charges", "check_for_electroneutrality" + return ("gamma", "reactant_types", "reactant_coefficients", + "product_types", "product_coefficients", "default_charges", + "check_for_electroneutrality") def _required_keys_add(self): return ["gamma", "reactant_types", "reactant_coefficients", @@ -301,7 +305,7 @@ cdef class ReactionAlgorithm: def get_status(self): """ Returns the status of the reaction ensemble in a dictionary containing - the used reactions, the used temperature and the used exclusion radius. + the used reactions, the used kT and the used exclusion radius. """ deref(self.RE).check_reaction_method() @@ -336,8 +340,8 @@ cdef class ReactionAlgorithm: "gamma": deref(self.RE).reactions[single_reaction_i].gamma} reactions.append(reaction) - return {"reactions": reactions, "temperature": deref( - self.RE).temperature, "exclusion_radius": deref(self.RE).exclusion_radius} + return {"reactions": reactions, "kT": deref( + self.RE).kT, "exclusion_radius": deref(self.RE).exclusion_radius} def delete_particle(self, p_id): """ @@ -406,12 +410,13 @@ cdef class ReactionEnsemble(ReactionAlgorithm): cdef unique_ptr[CReactionEnsemble] REptr def __init__(self, *args, **kwargs): - self._params = {"temperature": 1, + self._params = {"kT": 1, "exclusion_radius": 0} for k in self._required_keys(): if k not in kwargs: raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys().__str__() + " got " + kwargs.__str__()) + f"At least the following keys have to be given as keyword " + f"arguments: {self._required_keys()}, got {kwargs}") self._params[k] = kwargs[k] self.REptr.reset(new CReactionEnsemble(int(self._params["seed"]))) @@ -437,12 +442,12 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm): cdef unique_ptr[CConstantpHEnsemble] constpHptr def __init__(self, *args, **kwargs): - self._params = {"temperature": 1, - "exclusion_radius": 0} + self._params = {"kT": 1, "exclusion_radius": 0} for k in self._required_keys(): if k not in kwargs: raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys().__str__() + " got " + kwargs.__str__()) + f"At least the following keys have to be given as keyword " + f"arguments: {self._required_keys()}, got {kwargs}") self._params[k] = kwargs[k] self.constpHptr.reset(new CConstantpHEnsemble(int(self._params["seed"]))) @@ -457,12 +462,39 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm): self._set_params_in_es_core() 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.") - if(kwargs["reactant_coefficients"][0] != 1 or kwargs["product_coefficients"][0] != 1 or kwargs["product_coefficients"][1] != 1): - raise ValueError( - "All product and reactant coefficients must equal one in the constant pH method as implemented in ESPResSo.") + "The constant pH method is only implemented for reactions " + "with two product types and one adduct type.") + super().add_reaction(*args, **kwargs) property constant_pH: @@ -490,29 +522,30 @@ cdef class WidomInsertion(ReactionAlgorithm): cdef unique_ptr[CWidomInsertion] WidomInsertionPtr def _required_keys(self): - return "temperature", "seed" + return ("kT", "seed") def _valid_keys(self): - return "temperature", "seed" + return ("kT", "seed") def _valid_keys_add(self): - return "reactant_types", "reactant_coefficients", "product_types", "product_coefficients", "default_charges", "check_for_electroneutrality" + return ("reactant_types", "reactant_coefficients", "product_types", + "product_coefficients", "default_charges", + "check_for_electroneutrality") def _required_keys_add(self): - return ["reactant_types", "reactant_coefficients", - "product_types", "product_coefficients", "default_charges"] + return ("reactant_types", "reactant_coefficients", + "product_types", "product_coefficients", "default_charges") def __init__(self, *args, **kwargs): - self._params = {"temperature": 1} + self._params = {"kT": 1} for k in self._required_keys(): if k not in kwargs: raise ValueError( - "At least the following keys have to be given as keyword arguments: " + self._required_keys().__str__() + " got " + kwargs.__str__()) + f"At least the following keys have to be given as keyword " + f"arguments: {self._required_keys()}, got {kwargs}") self._params[k] = kwargs[k] - self._params[ - "exclusion_radius"] = 0.0 # this is not used by the widom insertion method - self._params[ - "gamma"] = 1.0 # this is not used by the widom insertion method + self._params["exclusion_radius"] = 0.0 # not used by this method + self._params["gamma"] = 1.0 # not used by this method self.WidomInsertionPtr.reset(new CWidomInsertion(int(self._params["seed"]))) self.RE = self.WidomInsertionPtr.get() diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index 1f7a67c02ab..accca5e22f7 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -45,15 +45,13 @@ def test_ideal_alpha(self): type=N0 * [type_A, type_H]) RE = espressomd.reaction_ensemble.ConstantpHEnsemble( - temperature=1.0, + kT=1.0, exclusion_radius=1, seed=44) RE.add_reaction( gamma=10**(-pKa), reactant_types=[type_HA], - reactant_coefficients=[1], product_types=[type_A, type_H], - product_coefficients=[1, 1], default_charges={type_HA: 0, type_A: -1, type_H: +1}) RE.constant_pH = pH diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index dcf8355102d..34e24bd6a79 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -45,7 +45,7 @@ class ReactionEnsembleTest(ut.TestCase): system.cell_system.skin = 0.4 system.time_step = 0.01 RE = espressomd.reaction_ensemble.ConstantpHEnsemble( - temperature=1.0, exclusion_radius=1, seed=44) + kT=1.0, exclusion_radius=1, seed=44) @classmethod def setUpClass(cls): @@ -56,9 +56,7 @@ def setUpClass(cls): cls.RE.add_reaction( gamma=cls.Ka, reactant_types=[cls.type_HA], - reactant_coefficients=[1], product_types=[cls.type_A, cls.type_H], - product_coefficients=[1, 1], default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}) cls.RE.constant_pH = cls.pH diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index a1fd42f0fd3..00471545038 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -61,7 +61,7 @@ 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_ensemble.ReactionEnsemble( - temperature=temperature, + kT=temperature, exclusion_radius=exclusion_radius, seed=12) @classmethod @@ -151,9 +151,9 @@ def test_reaction_system(self): self.assertAlmostEqual( ReactionEnsembleTest.temperature, - RE_status["temperature"], + RE_status["kT"], places=9, - msg="reaction ensemble temperature not set correctly.") + msg="reaction ensemble kT not set correctly.") self.assertAlmostEqual( ReactionEnsembleTest.exclusion_radius, RE_status["exclusion_radius"], diff --git a/testsuite/python/widom_insertion.py b/testsuite/python/widom_insertion.py index cc9b69e53a3..4ee45129d58 100644 --- a/testsuite/python/widom_insertion.py +++ b/testsuite/python/widom_insertion.py @@ -73,7 +73,7 @@ class WidomInsertionTest(ut.TestCase): volume = system.volume() Widom = espressomd.reaction_ensemble.WidomInsertion( - temperature=TEMPERATURE, seed=1) + kT=TEMPERATURE, seed=1) def setUp(self): self.system.part.add(pos=0.5 * self.system.box_l, type=self.TYPE_HA)