From 0a380c1de9dfeba1da8ce2c1f80eceb550ac3826 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 28 Jul 2023 23:59:18 +0200 Subject: [PATCH] Rapid prototyping of Monte Carlo methods MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Provide an alternative implementation of SingleReaction, ReactionAlgorithm and ConstantpHEnsemble in pure Python for rapid prototyping of new Monte Carlo methods. Co-authored-by: Pablo Miguel Blanco Andrés --- doc/sphinx/reaction_methods.rst | 39 +- samples/monte_carlo.py | 519 ++++++++++++++++++ src/python/espressomd/analyze.py | 10 + src/python/espressomd/reaction_methods.py | 2 + src/script_interface/analysis/Analysis.cpp | 15 + testsuite/scripts/samples/CMakeLists.txt | 4 + testsuite/scripts/samples/test_monte_carlo.py | 47 ++ 7 files changed, 635 insertions(+), 1 deletion(-) create mode 100644 samples/monte_carlo.py create mode 100644 testsuite/scripts/samples/test_monte_carlo.py diff --git a/doc/sphinx/reaction_methods.rst b/doc/sphinx/reaction_methods.rst index 9a0d7140eb..91a93f4dd9 100644 --- a/doc/sphinx/reaction_methods.rst +++ b/doc/sphinx/reaction_methods.rst @@ -379,4 +379,41 @@ 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`. +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() +`. +Reaction method classes override this function with their custom expression +for the acceptance probability. + +Alternatively, the sample script :file:`samples/monte_carlo.py` provides +a re-implementation of the core functionality of reaction methods in Python, +with a focus on the :ref:`constant pH ` and +:ref:`reaction ensemble ` methods. +More specifically, the :class:`~espressomd.reaction_methods.SingleReaction`, +:class:`~espressomd.reaction_methods.ReactionAlgorithm`, +:class:`~espressomd.reaction_methods.ReactionEnsemble`, and +:class:`~espressomd.reaction_methods.ConstantpHEnsemble` +classes are rewritten in Python. + +The goal of this sample is to assist in the rapid prototyping of new Monte Carlo +methods. In particular, the sampling and move generation schemes are expressed +in Python, and can be leveraged by users without C++ programming experience. +The sample is designed to run with the :ref:`kernprof` profiler attached: + +.. code-block:: bash + + pypresso --kernprof monte_carlo.py --mode=core + pypresso --kernprof monte_carlo.py --mode=python + +These Python implementations are roughly four times slower +than their corresponding core implementations. diff --git a/samples/monte_carlo.py b/samples/monte_carlo.py new file mode 100644 index 0000000000..c4c2f1ff48 --- /dev/null +++ b/samples/monte_carlo.py @@ -0,0 +1,519 @@ +# +# 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 . +# + +""" +Rapid prototyping of new Monte Carlo methods in Python. + +This sample provides a re-implementation of the core functionality +of :ref:`reaction methods ` in Python, +with a focus on the :ref:`constant pH ` and +:ref:`reaction ensemble ` methods. +See :ref:`Writing new Monte Carlo methods` for more details. +The sample simulates the acid-based titration of polyelectrolyte chains. + +The sample is designed to run with the :ref:`kernprof` profiler attached: + +.. code-block:: bash + + pypresso --kernprof monte_carlo.py --mode=core + pypresso --kernprof monte_carlo.py --mode=python + +""" + +import numpy as np +import itertools +import argparse +import math +import time +import tqdm +import os + +import espressomd +import espressomd.polymer +import espressomd.electrostatics +import espressomd.reaction_methods + +required_features = ["P3M", "WCA"] +espressomd.assert_features(required_features) + +parser = argparse.ArgumentParser( + prog=f"pypresso --kernprof {os.path.basename(__file__)}", + epilog=__doc__.lstrip().split("\n", 1)[0]) +parser.add_argument("--mode", choices=["core", "python"], default="python", + help="use C++ (core) or Python (python) implementation") +parser.add_argument("--method", choices=["cph", "re"], default="cph", + help="use constant pH (cph) or reaction ensemble (re)") +args = parser.parse_args() + + +if "line_profiler" not in dir(): + def profile(func): + def wrapper(*args, **kwargs): + return func(*args, **kwargs) + return wrapper + + +def factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0): + value = 1. + if nu_i > 0: + value /= math.factorial(N_i0 + nu_i) // math.factorial(N_i0) + elif nu_i < 0: + value *= math.factorial(N_i0) // math.factorial(N_i0 + nu_i) + return value + + +class SingleReaction: + + def __init__(self, **kwargs): + self.reactant_types = kwargs["reactant_types"] + self.reactant_coefficients = kwargs["reactant_coefficients"] + self.product_types = kwargs["product_types"] + self.product_coefficients = kwargs["product_coefficients"] + self.gamma = kwargs["gamma"] + self.accepted_moves = 0 + self.trial_moves = 0 + self.accumulator_potential_energy_difference_exponential = [] + self.nu_bar = sum(self.product_coefficients) - \ + sum(self.reactant_coefficients) + + def get_acceptance_rate(self): + return self.accepted_moves / self.trial_moves + + 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) + + +class ReactionAlgorithm: + + def __init__(self, **kwargs): + self.system = kwargs["system"] + self.kT = kwargs["kT"] + self.non_interacting_type = 100 + self.reactions = [] + self.particle_inside_exclusion_range_touched = False + self.default_charges = {} + self.m_empty_p_ids_smaller_than_max_seen_particle = [] + self.rng = np.random.default_rng(seed=kwargs["seed"]) + self.exclusion = espressomd.reaction_methods.ExclusionRadius(**kwargs) + + def set_non_interacting_type(self, type): + self.non_interacting_type = type + + def add_reaction(self, **kwargs): + """ + Set up a reaction in the forward and backward directions. + """ + self.default_charges.update(kwargs["default_charges"]) + forward_reaction = SingleReaction( + reactant_types=kwargs["reactant_types"], + reactant_coefficients=kwargs["reactant_coefficients"], + product_types=kwargs["product_types"], + product_coefficients=kwargs["product_coefficients"], + gamma=kwargs["gamma"]) + backward_reaction = forward_reaction.make_backward_reaction() + self.reactions.append(forward_reaction) + self.reactions.append(backward_reaction) + + @profile + def reaction(self, steps): + """ + Perform reaction steps. Chemical reactions are selected at random. + """ + E_pot = self.system.analysis.potential_energy() + random = self.rng.choice(len(self.reactions), size=steps, replace=True) + for i in random: + E_pot = self.generic_oneway_reaction(self.reactions[i], E_pot) + + @profile + def generic_oneway_reaction(self, reaction, E_pot_old): + """ + Carry out a generic one-way chemical reaction of the type + `A + B + ... --> C + D + ...` and return the new potential + energy if the trial move is accepted. + """ + reaction.trial_moves += 1 + self.particle_inside_exclusion_range_touched = False + if not self.all_reactant_particles_exist(reaction): + return E_pot_old + + old_particle_numbers = self.save_old_particle_numbers(reaction) + p_properties_old_state = self.make_reaction_attempt(reaction) + + if self.particle_inside_exclusion_range_touched: # reject trial move + self.restore_system(p_properties_old_state) + self.particle_inside_exclusion_range_touched = False + return E_pot_old + + E_pot_new = self.system.analysis.potential_energy() + E_pot_diff = E_pot_new - E_pot_old + bf = self.calculate_acceptance_probability( + reaction, E_pot_diff, old_particle_numbers) + reaction.accumulator_potential_energy_difference_exponential.append( + math.exp(-E_pot_diff / self.kT)) + if self.rng.uniform() < bf: # accept trial move + self.delete_hidden_particles(p_properties_old_state) + reaction.accepted_moves += 1 + return E_pot_new + else: # reject trial move + self.restore_system(p_properties_old_state) + return E_pot_old + + @profile + def make_reaction_attempt(self, reaction): + """ + Carry out a chemical reaction and save the old system state. + """ + minimum_number_of_types = min(len(reaction.reactant_types), + len(reaction.product_types)) + maximum_number_of_types = max(len(reaction.reactant_types), + len(reaction.product_types)) + p_properties_old_state = {"changed_particles": [], + "created_particles": [], + "hidden_particles": []} + + for index in range(minimum_number_of_types): + r_type = reaction.reactant_types[index] + p_type = reaction.product_types[index] + r_charge = self.default_charges[r_type] + p_charge = self.default_charges[p_type] + + # change reactant particles to product particles + size = min(reaction.reactant_coefficients[index], + reaction.product_coefficients[index]) + for random_pid in self.get_random_pids(r_type, size): + p = self.system.part.by_id(random_pid) + p.type = p_type + p.q = p_charge + p_properties_old_state["changed_particles"].append( + {"pid": random_pid, "type": r_type, "charge": r_charge}) + + # measure stoichiometric excess + delta_n = reaction.product_coefficients[index] - \ + reaction.reactant_coefficients[index] + + if delta_n > 0: + # create product particles + for _ in range(delta_n): + pid = self.create_particle(p_type) + self.check_exclusion_range(pid) + p_properties_old_state["created_particles"].append( + {"pid": pid, "type": p_type, "charge": p_charge}) + elif delta_n < 0: + # hide reactant particles + for random_pid in self.get_random_pids(r_type, -delta_n): + p_properties_old_state["hidden_particles"].append( + {"pid": random_pid, "type": r_type, "charge": r_charge}) + self.check_exclusion_range(random_pid) + self.hide_particle(random_pid) + + # create/hide particles with non-corresponding replacement types + for index in range(minimum_number_of_types, maximum_number_of_types): + if len(reaction.product_types) < len(reaction.reactant_types): + r_type = reaction.reactant_types[index] + r_charge = self.default_charges[r_type] + size = reaction.reactant_coefficients[index] + # hide superfluous reactant particles + for random_pid in self.get_random_pids(r_type, size): + p_properties_old_state["hidden_particles"].append( + {"pid": random_pid, "type": r_type, "charge": r_charge}) + self.check_exclusion_range(random_pid) + self.hide_particle(random_pid) + else: + p_type = reaction.product_types[index] + p_charge = self.default_charges[p_type] + # create additional product particles + for _ in range(reaction.product_coefficients[index]): + pid = self.create_particle(p_type) + self.check_exclusion_range(pid) + p_properties_old_state["created_particles"].append( + {"pid": pid, "type": p_type, "charge": p_charge}) + + return p_properties_old_state + + def calculate_acceptance_probability( + self, reaction, E_pot_diff, old_particle_numbers): + raise NotImplementedError("Derived classes must implement this method") + + def check_exclusion_range(self, pid): + self.particle_inside_exclusion_range_touched |= self.exclusion.check_exclusion_range( + pid=pid) + + def all_reactant_particles_exist(self, reaction): + for r_type in reaction.reactant_types: + r_index = reaction.reactant_types.index(r_type) + r_coef = reaction.reactant_coefficients[r_index] + if self.system.number_of_particles(type=r_type) < r_coef: + return False + return True + + def get_random_pids(self, ptype, size): + pids = self.system.analysis.call_method( + "get_pids_of_type", ptype=ptype) + indices = self.rng.choice(len(pids), size=size, replace=False) + return [pids[i] for i in indices] + + def save_old_particle_numbers(self, reaction): + old_particle_numbers = {} + for r_type in reaction.reactant_types + reaction.product_types: + old_particle_numbers[r_type] = self.system.number_of_particles( + type=r_type) + return old_particle_numbers + + def delete_created_particles(self, p_properties_old_state): + for particle_info in p_properties_old_state["created_particles"]: + self.system.part.by_id(particle_info["pid"]).remove() + + def delete_hidden_particles(self, p_properties_old_state): + for particle_info in p_properties_old_state["hidden_particles"]: + self.system.part.by_id(particle_info["pid"]).remove() + + def restore_system(self, p_properties_old_state): + # restore properties of changed and hidden particles + for particle_info in p_properties_old_state["changed_particles"] + \ + p_properties_old_state["hidden_particles"]: + p = self.system.part.by_id(particle_info["pid"]) + p.type = particle_info["type"] + p.q = particle_info["charge"] + # destroy created particles + self.delete_created_particles(p_properties_old_state) + + def hide_particle(self, pid): + p = self.system.part.by_id(pid) + p.type = self.non_interacting_type + p.q = 0. + + def create_particle(self, ptype): + if len(self.m_empty_p_ids_smaller_than_max_seen_particle) == 0: + pid = self.system.part.highest_particle_id + 1 + else: + pid = min(self.m_empty_p_ids_smaller_than_max_seen_particle) + self.system.part.add(id=pid, type=ptype, q=self.default_charges[ptype], + pos=self.rng.random((3,)) * self.system.box_l, + v=self.rng.normal(size=3) * math.sqrt(self.kT)) + return pid + + +class ReactionEnsemble(ReactionAlgorithm): + + def calculate_acceptance_probability( + self, reaction, E_pot_diff, old_particle_numbers): + """ + Calculate the acceptance probability of a Monte Carlo move. + """ + + volume = self.system.volume() + expr = math.exp(-E_pot_diff / self.kT) + expr *= volume**reaction.nu_bar * reaction.gamma + + # factorial contribution of reactants + for i in range(len(reaction.reactant_types)): + nu_i = -reaction.reactant_coefficients[i] + N_i0 = old_particle_numbers[reaction.reactant_types[i]] + expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0) + + # factorial contribution of products + for i in range(len(reaction.product_types)): + nu_i = reaction.product_coefficients[i] + N_i0 = old_particle_numbers[reaction.product_types[i]] + expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0) + + return expr + + +class ConstantpHEnsemble(ReactionAlgorithm): + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.constant_pH = kwargs["constant_pH"] + + def add_reaction(self, **kwargs): + kwargs["reactant_coefficients"] = [1] + kwargs["product_coefficients"] = [1, 1] + super().add_reaction(**kwargs) + + def calculate_acceptance_probability( + self, reaction, E_pot_diff, old_particle_numbers): + """ + Calculate the acceptance probability of a Monte Carlo move. + """ + + ln_bf = E_pot_diff - reaction.nu_bar * self.kT * math.log(10.) * ( + self.constant_pH + reaction.nu_bar * math.log10(reaction.gamma)) + + factorial_expr = math.exp(-ln_bf / self.kT) + + # factorial contribution of reactants + nu_i = -reaction.reactant_coefficients[0] + N_i0 = old_particle_numbers[reaction.reactant_types[0]] + factorial_expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0) + + # factorial contribution of products + nu_i = reaction.product_coefficients[0] + N_i0 = old_particle_numbers[reaction.product_types[0]] + factorial_expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0) + + return factorial_expr + + +# System parameters +############################################################# +box_l = 35 +system = espressomd.System(box_l=[box_l] * 3) +SEED = 42 +N_acid = 50 +N_chain = 5 +sigma = 1. +epsilon = 1. +system.time_step = 0.01 +system.cell_system.skin = 1.0 +N_steps_MD = 1000 +N_steps_MC = 50 + +# Reaction parameters +############################################################# +pKa = 7. +pH = 7.25 +kT = 1. +friction = 1. +types = { + "HA": 0, + "A-": 1, + "H+": 2, +} +charges = { + "HA": 0., + "A-": -1., + "H+": +1., +} +params = { + "kT": kT, + "exclusion_range": sigma, + "seed": SEED, +} +if args.method == "cph": + params["constant_pH"] = pH + +# Setup +############################################################# +np.random.seed(seed=SEED) + +positions = espressomd.polymer.linear_polymer_positions( + n_polymers=N_chain, beads_per_chain=N_acid // N_chain, + bond_length=sigma + 0.1, seed=SEED) + +bond = espressomd.interactions.HarmonicBond(k=10., r_0=sigma + 0.1) +system.bonded_inter.add(bond) +for polymer_pos in positions: + bond_partner = None + for pos in polymer_pos: + p = system.part.add(pos=pos, type=types["A-"], q=charges["A-"]) + if bond_partner: + p.add_bond((bond, bond_partner)) + bond_partner = p + +for _ in range(N_acid): + pos = np.random.random(3) * system.box_l + system.part.add(pos=pos, type=types["H+"], q=charges["H+"]) + +for type_pair in itertools.combinations_with_replacement(types.values(), 2): + system.non_bonded_inter[type_pair[0], type_pair[1]].wca.set_params( + epsilon=epsilon, sigma=sigma) + +p3m = espressomd.electrostatics.P3M( + prefactor=2., accuracy=1e-2, mesh=8, cao=3, verbose=False) +dh = espressomd.electrostatics.DH( + prefactor=2., kappa=0., r_cut=0.2 * box_l) + +# energy minimize the system +system.integrator.set_steepest_descent( + f_max=0., gamma=0.1, max_displacement=0.1 * sigma) +system.integrator.run(200) +system.electrostatics.solver = p3m +system.integrator.run(1000) + +# thermalize the system +system.integrator.set_vv() +system.thermostat.set_langevin(kT=kT, gamma=friction, seed=SEED) +system.integrator.run(1000) +system.electrostatics.solver = dh + +if args.mode == "core": + if args.method == "cph": + ReactionMethod = espressomd.reaction_methods.ConstantpHEnsemble + elif args.method == "re": + ReactionMethod = espressomd.reaction_methods.ReactionEnsemble +elif args.mode == "python": + if args.method == "cph": + ReactionMethod = ConstantpHEnsemble + elif args.method == "re": + ReactionMethod = ReactionEnsemble + params["system"] = system + +# set up reaction method +RE = ReactionMethod(**params) +RE.set_non_interacting_type(type=max(types.values()) + 1) +if args.method == "cph": + RE.add_reaction( + gamma=10**-pKa, + reactant_types=[types["HA"]], + product_types=[types["A-"], types["H+"]], + default_charges={types[name]: charges[name] for name in types.keys()}) +elif args.method == "re": + RE.add_reaction( + gamma=1e-3, + reactant_types=[types["HA"]], + reactant_coefficients=[1], + product_types=[types["A-"], types["H+"]], + product_coefficients=[1, 1], + default_charges={types[name]: charges[name] for name in types.keys()}) +reaction = RE.reactions[0] +system.setup_type_map(type_list=list(types.values())) + +# equilibrate the polyelectrolyte chains +for i in range(5): + RE.reaction(steps=10 * N_steps_MC) + system.integrator.run(N_steps_MD) + + +@profile +def sample_alpha(length): + alpha_list = [] + for _ in tqdm.tqdm(range(length)): + system.integrator.run(steps=N_steps_MD) + RE.reaction(steps=N_steps_MC) + alpha = system.number_of_particles(type=types["A-"]) / N_acid + alpha_list.append(alpha) + return alpha_list + + +sample_size = 100 +tick = time.time() +alphas = sample_alpha(sample_size) +tock = time.time() + +alpha_avg = np.mean(alphas) +alpha_err = 1.96 * np.sqrt(np.var(alphas) / len(alphas)) +acceptance_rate = reaction.get_acceptance_rate() +print(f"acceptance rate = {100. * acceptance_rate:.0f}%") +print(f"alpha = {alpha_avg:.2f} +/- {alpha_err:.2f}") +print(f"runtime = {tock - tick:.2f}s") diff --git a/src/python/espressomd/analyze.py b/src/python/espressomd/analyze.py index bec8e00b49..c36e98c38c 100644 --- a/src/python/espressomd/analyze.py +++ b/src/python/espressomd/analyze.py @@ -380,6 +380,15 @@ class Analysis(ScriptInterfaceHelper): Where [0] contains the midpoints of the bins, and [1] contains the values of the rdf. + potential_energy() + Calculate the potential energy of the system, i.e. the total energy + minus the kinetic energy. + + Returns + ------- + :obj: `float` + Potential energy. + """ _so_name = "ScriptInterface::Analysis::Analysis" _so_creation_policy = "GLOBAL" @@ -387,6 +396,7 @@ class Analysis(ScriptInterfaceHelper): "linear_momentum", "center_of_mass", "nbhood", + "potential_energy", "particle_neighbor_pids", "calc_re", "calc_rg", diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py index 768b1a2072..4bea1d8d44 100644 --- a/src/python/espressomd/reaction_methods.py +++ b/src/python/espressomd/reaction_methods.py @@ -117,6 +117,8 @@ class ReactionAlgorithm(ScriptInterfaceHelper): the Reaction Ensemble algorithm and the constant pH method. Initialize the reaction algorithm by setting the standard pressure, temperature, and the exclusion range. + The exclusion range mechanism is explained in more detail + in :class:`~espressomd.reaction_methods.ExclusionRadius`. Note: When creating particles the velocities of the new particles are set according the Maxwell-Boltzmann distribution. In this step the mass of the diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index d7f9f20f53..6e1e3e282f 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -96,6 +96,10 @@ Variant Analysis::do_call_method(std::string const &name, return {}; } #endif + if (name == "potential_energy") { + auto const obs = calculate_energy(); + return obs->accumulate(-obs->kinetic[0]); + } if (name == "particle_neighbor_pids") { on_observable_calc(); std::unordered_map> dict; @@ -109,6 +113,17 @@ Variant Analysis::do_call_method(std::string const &name, }); return make_unordered_map_of_variants(dict); } + if (name == "get_pids_of_type") { + auto const type = get_value(parameters, "ptype"); + std::vector pids; + for (auto const &p : ::cell_structure.local_particles()) { + if (p.type() == type) { + pids.push_back(p.id()); + } + } + Utils::Mpi::gather_buffer(pids, context()->get_comm()); + return pids; + } #ifdef DPD if (name == "dpd_stress") { auto const result = dpd_stress(context()->get_comm()); diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 78d4a9c267..989304e166 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -77,6 +77,10 @@ sample_test(FILE test_p3m.py SUFFIX gpu LABELS "gpu") sample_test(FILE test_reaction_methods.py SUFFIX constant_pH_ensemble) sample_test(FILE test_reaction_methods.py SUFFIX reaction_ensemble) sample_test(FILE test_reaction_ensemble_complex_reaction.py) +sample_test(FILE test_monte_carlo.py SUFFIX core_cph) +sample_test(FILE test_monte_carlo.py SUFFIX core_re) +sample_test(FILE test_monte_carlo.py SUFFIX python_cph) +sample_test(FILE test_monte_carlo.py SUFFIX python_re) sample_test(FILE test_rigid_body.py) sample_test(FILE test_save_checkpoint.py) set_tests_properties(sample_save_checkpoint PROPERTIES FIXTURES_SETUP diff --git a/testsuite/scripts/samples/test_monte_carlo.py b/testsuite/scripts/samples/test_monte_carlo.py new file mode 100644 index 0000000000..1653f43e93 --- /dev/null +++ b/testsuite/scripts/samples/test_monte_carlo.py @@ -0,0 +1,47 @@ +# +# 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 . +# + +import unittest as ut +import importlib_wrapper + +mode, method = "@TEST_SUFFIX@".split("_") +assert method in ("cph", "re") + +sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@SAMPLES_DIR@/monte_carlo.py", script_suffix="@TEST_SUFFIX@", + cmd_arguments=["--mode", mode, "--method", method], sample_size=150) + + +@skipIfMissingFeatures +class Sample(ut.TestCase): + system = sample.system + + def test(self): + if method == "cph": + self.assertAlmostEqual(sample.alpha_avg, 0.29, delta=0.05) + self.assertAlmostEqual(sample.alpha_err, 0.01, delta=0.02) + self.assertAlmostEqual(sample.acceptance_rate, 0.56, delta=0.10) + else: + self.assertAlmostEqual(sample.alpha_avg, 0.33, delta=0.05) + self.assertAlmostEqual(sample.alpha_err, 0.01, delta=0.02) + self.assertAlmostEqual(sample.acceptance_rate, 0.55, delta=0.10) + + +if __name__ == "__main__": + ut.main()