diff --git a/samples/gibbs_ensemble/Readme.md b/samples/gibbs_ensemble/Readme.md
new file mode 100644
index 00000000000..abd77c7619a
--- /dev/null
+++ b/samples/gibbs_ensemble/Readme.md
@@ -0,0 +1,17 @@
+# Gibbs ensemble simulations using ESPResSo
+
+This sample code shows how to implement the Gibbs ensemble into ESPResSo using
+the `multiprocessing` module to create two separate systems and communicate to
+them through the main script.
+
+The program consists of 3 files: In the `run_sim.py` script, the client boxes
+are started and given instructions. The clients use the `Gibbs_Client` subclass
+in `client_system.py` which inherits the necessary methods for communication and
+for the trial moves from the `Client` class in `gibbs_ensemble.py`. This subclass
+allows for energy correction terms as is seen in this sample. Because ESPResSo
+only allows one instance of a system in one process, the system is created using
+the `set_up_system()` function and initialized once the `run()` method in the
+`Client` class is called.
+
+All equations are taken from Frenkel, Smit: *Understanding Molecular Simulation* 2002,
+doi:[10.1016/B978-0-12-267351-1.X5000-7](https://doi.org/10.1016/B978-0-12-267351-1.X5000-7).
diff --git a/samples/gibbs_ensemble/client_system.py b/samples/gibbs_ensemble/client_system.py
new file mode 100644
index 00000000000..1b7ef3abbcd
--- /dev/null
+++ b/samples/gibbs_ensemble/client_system.py
@@ -0,0 +1,87 @@
+#
+# Copyright (C) 2013-2021 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 .
+#
+"""
+This script uses the System class from the :file:`gibbs_ensemble.py`
+script and creates a subclass in which you can create your system, add
+energy corrections etc.
+"""
+
+import espressomd
+import gibbs_ensemble
+import numpy as np
+
+espressomd.assert_features("LENNARD_JONES")
+
+# Lennard Jones parameters
+LJ_EPSILON = 1.0
+LJ_SIGMA = 1.0
+LJ_CUTOFF = 2.5 * LJ_SIGMA
+LJ_SHIFT = - (np.power(LJ_SIGMA / LJ_CUTOFF, 12) -
+ np.power(LJ_SIGMA / LJ_CUTOFF, 6))
+
+
+def calc_shift_correction(density):
+ """ Shift correction by integrating the constant part from 0 to the cutoff
+ radius """
+ return 4. / 3. * np.pi * LJ_CUTOFF**3 * LJ_SHIFT * 4 * LJ_EPSILON * density
+
+
+def calc_tail_correction(density):
+ """ Lennard-Jones truncation correction (See Frenkel, Smit, Understanding
+ molecular simulation, eq. (3.2.5)) """
+ return 8. / 3. * np.pi * density * LJ_EPSILON * LJ_SIGMA**3 * \
+ (1. / 3. * np.power(LJ_SIGMA / LJ_EPSILON, 9) -
+ np.power(LJ_SIGMA / LJ_EPSILON, 3))
+
+
+class Gibbs_Client(gibbs_ensemble.Client):
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+
+ @property
+ def energy(self):
+ """ Additional correction terms """
+ return super().energy # + \
+ # calc_tail_correction(self.density) + \
+ # calc_shift_correction(self.density)
+
+ def init_system(self):
+ """ Initialize the system (this is executed only when the run
+ method of the process is executed) """
+ self.system = set_up_system(self.init_box_length, self.init_num_part)
+
+
+def set_up_system(box_length, num_part):
+ """ Use this function to create the system """
+
+ system = espressomd.System(box_l=3 * [box_length])
+
+ system.cell_system.set_n_square(use_verlet_lists=False)
+ system.cell_system.skin = 0
+ system.time_step = 0.01
+
+ system.part.add(pos=np.random.random((num_part, 3)) * box_length)
+
+ system.non_bonded_inter[0, 0].lennard_jones.set_params(
+ epsilon=LJ_EPSILON,
+ sigma=LJ_SIGMA,
+ cutoff=LJ_CUTOFF,
+ shift=LJ_SHIFT)
+
+ return system
diff --git a/samples/gibbs_ensemble/create_fits.py b/samples/gibbs_ensemble/create_fits.py
new file mode 100644
index 00000000000..74e1dc922e8
--- /dev/null
+++ b/samples/gibbs_ensemble/create_fits.py
@@ -0,0 +1,232 @@
+#
+# Copyright (C) 2021 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 .
+#
+"""
+Take in all simulation data, compute equilibrium values for every
+datafile and take means for simulations with same temperature.
+Plot the phase diagram and the fitted critical point.
+"""
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy import signal
+import pickle
+import gzip
+import matplotlib.pyplot as plt
+import argparse
+import logging
+
+
+parser = argparse.ArgumentParser(
+ description="""Computes mean values for every file given as input,
+ computes the mean values for same temperatures and different seeds,
+ tries to compute the critical point using rectilinear and scaling
+ law and plots the phase diagram.""")
+parser.add_argument(
+ 'files',
+ nargs='+',
+ help="(All) files of the simulated system")
+args = parser.parse_args()
+
+# Set the logging level
+logging.basicConfig(
+ format='%(asctime)s %(levelname)s: %(message)s',
+ datefmt='%H:%M:%S',
+ level=logging.INFO
+)
+
+# Warmup length to skip at beginning
+WARMUP_LENGTH = int(3e3)
+
+# Lists for (unsorted) densities and temperature
+dens_liquid = []
+dens_gas = []
+errors_liquid = []
+errors_gas = []
+kTs = []
+
+# 3D critical exponent (See Frenkel, Smit: Understanding Molecular
+# Simulation 2002, p.217)
+BETA = 0.32
+
+
+def autocorr(x):
+ """ Compute the normalized autocorrelation function """
+ _x = x - np.mean(x)
+ return (signal.convolve(
+ _x, _x[::-1], mode='full', method='auto')[(np.size(_x) - 1):]) / (np.sum(_x * _x))
+
+
+def relaxation_time(x):
+ """ Compute the relaxation time """
+ corr = autocorr(x)
+ popt, _ = curve_fit(lambda x, a: np.exp(-a * x),
+ range(np.size(corr)), corr, p0=(1. / 15000))
+ return popt[0]
+
+
+def std_error_mean(x):
+ """ Compute the standard error of the mean with correlated samples """
+ autocorr_time = relaxation_time(x)
+ N_eff = np.size(x) / (2 * autocorr_time)
+ return np.sqrt(np.var(x) / N_eff)
+
+
+def scaling_law(kT, B, kT_c):
+ """ Scaling law, see Frenkel, Smit, eq. (8.3.6) """
+ return B * np.abs(kT - kT_c)**BETA
+
+
+def rectilinear_law(kT, p_c, A, kT_c):
+ """ Rectilinear law, see Frenkel, Smit, eq. (8.3.7) """
+ return p_c + A * np.abs(kT - kT_c)
+
+
+# Load data
+for f in args.files:
+
+ logging.info("Loading {}...".format(f))
+ try:
+ with gzip.open(f, 'rb') as _f:
+ data = pickle.load(_f)
+ except FileNotFoundError:
+ raise FileNotFoundError("File {} not found.".format(f))
+
+ kT = data["temperature"]
+
+ # Calculate densities of both phases
+ density1 = np.mean(data["densities_box01"][WARMUP_LENGTH:])
+ density2 = np.mean(data["densities_box02"][WARMUP_LENGTH:])
+
+ dens_liquid.append(np.max((density1, density2)))
+ dens_gas.append(np.min((density1, density2)))
+
+ # Calculate errors
+ errors_liquid.append(
+ std_error_mean(data["densities_box01"]) if density1 > density2 else
+ std_error_mean(data["densities_box02"]))
+ errors_gas.append(
+ std_error_mean(data["densities_box01"]) if density1 < density2 else
+ std_error_mean(data["densities_box02"]))
+ kTs.append(kT)
+
+# Sort temperatures and densities respectively
+order = np.argsort(kTs)
+kTs = np.array(kTs)[order]
+dens_liquid = np.array(dens_liquid)[order]
+dens_gas = np.array(dens_gas)[order]
+errors_liquid = np.array(errors_liquid)[order]
+errors_gas = np.array(errors_gas)[order]
+
+# Average over simulation results with different seeds
+dens_liquid = [np.mean([dens_liquid[i] for (i, T) in enumerate(
+ kTs) if T == kT]) for kT in np.unique(kTs)]
+dens_gas = [np.mean([dens_gas[i] for (i, T) in enumerate(kTs) if T == kT])
+ for kT in np.unique(kTs)]
+
+# Compute Errors as maximum of each simulation error
+errors_liquid = [np.max([errors_liquid[i] for (i, T) in enumerate(
+ kTs) if T == kT]) for kT in np.unique(kTs)]
+errors_gas = [np.max([errors_gas[i] for (i, T) in enumerate(kTs) if T == kT])
+ for kT in np.unique(kTs)]
+kTs = np.unique(kTs)
+
+# Convert to arrays to do math
+dens_liquid = np.array(dens_liquid)
+dens_gas = np.array(dens_gas)
+
+# Needed for scaling law and rectilinear law
+y_scaling = dens_liquid - dens_gas
+y_rectilinear = 0.5 * (dens_liquid + dens_gas)
+
+# Fit using scaling law
+fit_scaling, p_cov_scaling = curve_fit(
+ scaling_law, kTs, y_scaling, p0=(1., 1.), maxfev=6000)
+
+# Print fit values
+logging.info("Fits using scaling law: B, T_c: {}".format(fit_scaling))
+logging.info("Fit uncertainty: {}".format(np.diag(p_cov_scaling)))
+
+# Critical temperature obtained via scaling law
+kT_c_scaling = fit_scaling[1]
+
+# Fit using rectilinear law
+fit_rectilinear, p_cov_rectilinear = curve_fit(
+ rectilinear_law, kTs, y_rectilinear, p0=(
+ 0.3, 0.2, kT_c_scaling), maxfev=6000)
+
+# Print fit values
+logging.info(
+ "Fits using rectilinear law: p_c, A, T_c: {}".format(fit_rectilinear))
+logging.info("Fit uncertainty: {}".format(np.diag(p_cov_rectilinear)))
+
+# Critical point obtained via rectilinear law
+p_c = fit_rectilinear[0]
+kT_c = fit_rectilinear[2]
+
+# Save results
+with gzip.open('result.dat.gz', 'wb') as f:
+ pickle.dump(
+ {"dens_liquid": dens_liquid,
+ "dens_gas": dens_gas,
+ "errors_liquid": errors_liquid,
+ "errors_gas": errors_gas,
+ "kTs": kTs,
+ "kT_c_scaling": kT_c,
+ "p_c": p_c},
+ f)
+
+# Plot results
+ms = 40
+
+plt.errorbar(dens_liquid,
+ kTs,
+ xerr=errors_liquid,
+ fmt='o',
+ c='red')
+
+plt.errorbar(dens_gas,
+ kTs,
+ xerr=errors_gas,
+ fmt='o',
+ c='red',
+ label='simulation data')
+
+plt.scatter(p_c,
+ kT_c,
+ s=ms, c='black',
+ marker='o',
+ label='crit. point')
+
+plt.scatter(0.5 * (dens_liquid + dens_gas),
+ kTs,
+ s=ms - 10, c='black',
+ marker='x')
+
+plt.plot(rectilinear_law(kTs,
+ p_c,
+ fit_rectilinear[1],
+ fit_rectilinear[2]),
+ kTs)
+
+plt.xlabel(r"Particle density $\rho*$")
+plt.ylabel(r"Temperature $T*$")
+plt.legend()
+
+plt.savefig("phase_diagram.pdf")
+plt.show()
diff --git a/samples/gibbs_ensemble/gibbs_ensemble.py b/samples/gibbs_ensemble/gibbs_ensemble.py
new file mode 100644
index 00000000000..8a66ee39878
--- /dev/null
+++ b/samples/gibbs_ensemble/gibbs_ensemble.py
@@ -0,0 +1,157 @@
+#
+# Copyright (C) 2013-2021 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 .
+#
+"""
+Client class containing the methods used for communicating and manipulating
+the system MessageID class for communication identifiers.
+"""
+
+from enum import unique, Enum, auto
+from multiprocessing import Process
+import os
+import numpy as np
+
+
+@unique
+class MessageID(Enum):
+ # Message identifiers
+ ENERGY = auto()
+ END = auto()
+ MOVE_PART = auto()
+ MOVE_PART_REVERT = auto()
+ CHANGE_VOLUME = auto()
+ CHANGE_VOLUME_REVERT = auto()
+ PART_ADD = auto()
+ PART_ADD_REVERT = auto()
+ PART_REMOVE = auto()
+ PART_REMOVE_REVERT = auto()
+ INFO = auto()
+
+
+class Client(Process):
+ def __init__(self,
+ box_name,
+ pipe,
+ seed,
+ box_length,
+ num_part):
+ super().__init__(name=box_name)
+ self.init_box_length = box_length
+ self.init_num_part = num_part
+ self.system = None
+ self.old_box_length = 0.0
+ # Pipe for communicating
+ self.pipe = pipe
+ self.seed = seed
+
+ @property
+ def density(self):
+ return len(self.system.part[:].id) / float(np.prod(self.system.box_l))
+
+ @property
+ def energy(self):
+ return self.system.analysis.energy()["total"]
+
+ def send_energy(self):
+ self.pipe.send([MessageID.ENERGY, self.energy])
+
+ def end_simulation(self):
+ # Additional simulation checkpointing etc.
+ os._exit(os.EX_OK)
+
+ def send_info(self):
+ self.pipe.send([MessageID.INFO, self.system.box_l[0],
+ len(self.system.part[:].id)])
+
+ def move_particle(self, DX_MAX):
+ """ Move random particle inside the box """
+ self.old_particle = np.random.choice(self.system.part)
+ self.old_pos = self.old_particle.pos
+ self.old_particle.pos = self.old_pos + \
+ (0.5 - np.random.random(size=3)) * DX_MAX
+ self.send_energy()
+
+ def revert_move_particle(self):
+ """ Revert particle move """
+ self.old_particle.pos = self.old_pos
+ self.old_particle = None
+
+ def add_particle(self):
+ """ Add a particle at random inside box """
+ self.last_added_particle = self.system.part.add(
+ pos=np.random.random(3) * self.system.box_l)
+ self.send_energy()
+
+ def revert_add_particle(self):
+ """ Revert last particle add """
+ self.last_added_particle.remove()
+ self.last_added_particle = None
+
+ def remove_particle(self):
+ """ Remove a random particle """
+ self.old_particle = np.random.choice(self.system.part).to_dict()
+ self.system.part[self.old_particle["id"]].remove()
+ self.send_energy()
+
+ def revert_remove_particle(self):
+ """ Revert last particle remove """
+ self.system.part.add(self.old_particle)
+ self.old_particle = None
+
+ def change_volume(self, new_box_l):
+ """ Change volume to new box length """
+ self.old_box_length = self.system.box_l[0]
+ self.system.change_volume_and_rescale_particles(new_box_l)
+ self.send_energy()
+
+ def revert_change_volume(self):
+ """ Revert volume change """
+ self.system.change_volume_and_rescale_particles(self.old_box_length)
+
+ def run(self):
+ """ Start while loop, wait for tasks and handle them until MessageID.END is sent """
+ # Set random seed
+ np.random.seed(self.seed)
+
+ # Start the system (needs to be done here because
+ # ESPResSo can only handle one system per process)
+ self.init_system()
+
+ # Send initial energy
+ self.send_energy()
+
+ messagehandler = {
+ MessageID.MOVE_PART: self.move_particle,
+ MessageID.MOVE_PART_REVERT: self.revert_move_particle,
+ MessageID.CHANGE_VOLUME: self.change_volume,
+ MessageID.CHANGE_VOLUME_REVERT: self.revert_change_volume,
+ MessageID.PART_ADD: self.add_particle,
+ MessageID.PART_ADD_REVERT: self.revert_add_particle,
+ MessageID.PART_REMOVE: self.remove_particle,
+ MessageID.PART_REMOVE_REVERT: self.revert_remove_particle,
+ MessageID.INFO: self.send_info,
+ MessageID.END: self.end_simulation
+ }
+
+ while True:
+ """ Can be left at while True, we can simply terminate the process via process.terminate() in the parent process """
+ # Catch errors during handling
+ msg = self.pipe.recv()
+
+ # Handle message
+ messagehandler[msg[0]](*msg[1:])
diff --git a/samples/gibbs_ensemble/run_for_several_kT.sh b/samples/gibbs_ensemble/run_for_several_kT.sh
new file mode 100755
index 00000000000..5b42ae9e73f
--- /dev/null
+++ b/samples/gibbs_ensemble/run_for_several_kT.sh
@@ -0,0 +1,43 @@
+#!/bin/sh
+#
+# Copyright (C) 2021 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 .
+#
+
+python=$(which python3)
+pypresso=../../build/pypresso
+
+if [ ! -e "${pypresso}" ]
+then
+ echo Invalid path to pypresso script: $pypresso
+ exit 1
+fi
+
+# Run the simulations for the temperatures given as parameters
+for kT in $*
+do
+ for seed in 1 2 3 4 5
+ do
+
+ test -s "temp_${kT}_seed_${seed}.dat.gz" ||
+ "${pypresso}" run_sim.py ${kT} --log --steps 1000000 --seed ${seed} > "temp_${kT}_seed_${seed}.out" 2>&1 &
+ done
+done
+wait
+
+# Plot the results, create fits
+"${python}" create_fits.py temp_*.dat.gz
diff --git a/samples/gibbs_ensemble/run_sim.py b/samples/gibbs_ensemble/run_sim.py
new file mode 100644
index 00000000000..c004d09e692
--- /dev/null
+++ b/samples/gibbs_ensemble/run_sim.py
@@ -0,0 +1,524 @@
+#
+# Copyright (C) 2013-2021 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 .
+#
+"""
+Simulate a Gibbs-ensemble of a Lennard-Jones fluid at a fixed temperature.
+This script does the Monte-Carlo part of the Gibbs-ensemble, however for the
+energy calculation of the systems, two instances of the Client class in
+:file:`gibbs_ensemble.py` initialized (and possibly updated) in
+:file:`gibbs_client_system.py` are started which each run a different instance
+of ESPResSo.
+
+The Gibbs-ensemble implemented in these scripts closely refers to chapter 8 of
+Daan Fenkel and Berend Smit 'Understanding Molecular Simulation, Second edition'.
+All equation references noted in this script can be found there.
+
+Due to the cutoff and shifting of the LJ-potential the simulated points in the
+phase diagram differ from the long range uncut LJ-potential. The phase diagram
+of the used potential can be found in 'B. Smit, J. Chem. Phys. 96 (11), 1992,
+Phase diagrams of Lennard-Jones fluids'.
+"""
+
+from multiprocessing import Pipe
+import numpy as np
+from tqdm import tqdm
+import gzip
+import pickle
+import argparse
+import logging
+from enum import Enum, unique
+
+from gibbs_ensemble import MessageID
+from client_system import Gibbs_Client
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('T', type=float, help="Temperature of the system")
+parser.add_argument(
+ '-N',
+ type=int,
+ default=256,
+ help="Number of particles, default 256")
+parser.add_argument(
+ '--density',
+ type=float,
+ default=0.3,
+ help="Init density, default 0.3")
+parser.add_argument(
+ '--steps',
+ '-s',
+ type=int,
+ default=100000,
+ help="Number of steps, default is 1e5")
+parser.add_argument(
+ '--warmup',
+ '-w',
+ type=int,
+ default=5000,
+ help="Warmup steps, default 5000")
+parser.add_argument(
+ '--seed',
+ type=int,
+ default=17,
+ help="Random seed")
+parser.add_argument(
+ '--log',
+ action='store_true',
+ help="Pipe information to log file instead of stdout")
+parser.add_argument(
+ '--debug',
+ action='store_true',
+ help="Set for debugging info (use with small step number)")
+args = parser.parse_args()
+
+
+@unique
+class Moves(Enum):
+ # Particle moves
+ MOVE_PARTICLE = 1
+ EXCHANGE_VOLUME = 2
+ EXCHANGE_PARTICLE = 3
+
+
+# Number of client processes
+NUM_OF_CLIENTS = 2
+
+# Global number of particles
+GLOBAL_NUM_PART = args.N
+try:
+ assert((GLOBAL_NUM_PART % 2) == 0)
+except AssertionError:
+ raise TypeError("Wrong number of particles, must be even")
+
+# Init particle density
+init_density = args.density
+
+# Global volume
+GLOBAL_VOLUME = GLOBAL_NUM_PART / init_density
+
+# Init box _length
+init_box_length = np.cbrt(GLOBAL_VOLUME / NUM_OF_CLIENTS)
+
+# Temperature
+kT = args.T
+
+# Displacement
+DX_MAX = 0.4
+DV_MAX = 0.2
+
+# Number of steps, warmup steps
+steps = args.steps
+warmup = args.warmup
+# Perform about 100 system checks during the simulation
+check = int(steps / 1000)
+
+# Trial move probabilities
+MOVE_CHANCE = 0.5
+VOLUME_CHANCE = 0.01
+EXCHANGE_CHANCE = 1 - MOVE_CHANCE - VOLUME_CHANCE
+
+# Random seed
+seed = args.seed
+np.random.seed(seed)
+box_seeds = np.random.randint(0, np.iinfo(np.int32).max, size=NUM_OF_CLIENTS)
+
+# Filename
+filename = "temp_{}_seed_{}".format(kT, seed)
+
+# Set the logging level
+logging.basicConfig(
+ filename=(filename + ".log" if args.log else ""),
+ format='%(asctime)s %(levelname)s: %(message)s',
+ datefmt='%H:%M:%S',
+ level=logging.DEBUG if args.debug else logging.INFO
+)
+
+# Number of moves, number of accepted moves
+num_moves = {
+ Moves.MOVE_PARTICLE: 0,
+ Moves.EXCHANGE_VOLUME: 0,
+ Moves.EXCHANGE_PARTICLE: 0}
+
+num_accepted_moves = {
+ Moves.MOVE_PARTICLE: 0,
+ Moves.EXCHANGE_VOLUME: 0,
+ Moves.EXCHANGE_PARTICLE: 0}
+
+
+class Box:
+ def __init__(self, box_name, pipe, seed, box_length, num_part):
+ self.box_name = box_name
+ self.pipe = pipe[0]
+ self.energy = 0.0
+ self.old_energy = 1e10
+ self.box_length = box_length
+ self.old_box_length = self.box_length
+ self.num_part = num_part
+ self.client = Gibbs_Client(
+ box_name,
+ pipe[1],
+ seed,
+ box_length,
+ num_part)
+
+ @property
+ def volume(self):
+ return self.box_length**3
+
+ @property
+ def density(self):
+ return self.num_part / self.volume
+
+ @property
+ def diff_energy(self):
+ return np.clip((self.energy - self.old_energy), -100 * kT, 100 * kT)
+
+ def send_data(self, data):
+ # For debugging
+ #assert(type(data) == list)
+ logging.debug("Send to {}: {}".format(self.box_name, data))
+ self.pipe.send(data)
+
+ def recv_data(self):
+ msg = self.pipe.recv()
+ logging.debug("Receive from {}: {}".format(self.box_name, msg))
+ return msg
+
+ def recv_energy(self):
+ msg = self.recv_data()
+ # Debugging
+ try:
+ assert(msg[0] == MessageID.ENERGY)
+ except BaseException:
+ raise ConnectionError(
+ "Error during energy return of {}, got this instead: \n{}".format(
+ self.box_name, msg))
+ self.old_energy = self.energy
+ self.energy = msg[1]
+
+
+def validate_info(boxes):
+ """ Check if info returned from boxes is the same as on the parent side """
+ [box.send_data([MessageID.INFO]) for box in boxes]
+ for box in boxes:
+ msg = box.pipe.recv()
+ try:
+ assert(msg[0] == MessageID.INFO)
+ except BaseException:
+ raise ConnectionError(
+ "Connection to {} seems to be broken.".format(
+ box.box_name))
+ np.testing.assert_equal(
+ box.box_length,
+ msg[1],
+ err_msg="Server side box length (actual) differs from client side (desired).")
+ np.testing.assert_equal(
+ box.num_part,
+ msg[2],
+ err_msg="Server side box length (actual) differs from client side (desired)")
+ logging.debug(
+ "Validation correct. Values of {}:\nBox length:\t{}\nNum Part:\t{}.".format(
+ box.box_name, box.box_length, box.num_part))
+
+
+def choose_random_box(boxes):
+ """ Choose a box at random (Assumes 2 Clients) """
+ random_int = np.random.randint(2)
+ try:
+ assert(boxes[random_int].num_part > 0)
+ return random_int
+ except BaseException:
+ return (random_int + 1) % 2
+
+
+def perform_move_particle(boxes):
+ """ Perform a particle move, check and revert if necessary """
+
+ # Choose random box
+ box = boxes[choose_random_box(boxes)]
+
+ # Send message to Client
+ box.send_data([MessageID.MOVE_PART, DX_MAX])
+ box.recv_energy()
+
+ # Debug
+ logging.debug("Performing particle move in {}.".format(box.box_name))
+
+ # ---- Check move ----
+
+ # Check move probability (see (8.3.1) in Frenkel, Smit)
+ acc = np.exp(- 1. / kT * box.diff_energy)
+
+ # Draw random number
+ rand = np.random.random()
+
+ # Check if move was valid
+ if rand > acc:
+ # Reject move, restore old configuration
+ box.send_data([MessageID.MOVE_PART_REVERT])
+ box.energy = box.old_energy
+ return False
+
+ # Debug
+ logging.debug("Particle move accepted.")
+
+ # Only if rand < acc
+ return True
+
+
+def perform_volume_change(boxes):
+ """ Perform a random move in log(V_1/V_2) """
+
+ # Store old box length, old volumes
+ for box in boxes:
+ box.old_box_length = box.box_length
+ old_volume_1 = boxes[0].volume
+ old_volume_2 = boxes[1].volume
+
+ # Debugging purposes
+ if args.debug:
+ np.testing.assert_almost_equal(
+ GLOBAL_VOLUME, old_volume_1 + old_volume_2)
+
+ # Random move in log(V_1/V_2) (See Algorithm 18 in Frenkel, Smit)
+ lnvn = np.log(old_volume_1 / old_volume_2) + \
+ (np.random.random() - 0.5) * DV_MAX
+
+ # Calculate new box lengths
+ new_volume_1 = GLOBAL_VOLUME * np.exp(lnvn) / (1 + np.exp(lnvn))
+ new_volume_2 = GLOBAL_VOLUME - new_volume_1
+ boxes[0].box_length = np.cbrt(new_volume_1)
+ boxes[1].box_length = np.cbrt(new_volume_2)
+
+ # Debug
+ logging.debug(
+ "Perform volume change: {}: {} -> {}; {}: {} -> {}.".format(
+ boxes[0].box_name,
+ old_volume_1,
+ new_volume_1,
+ boxes[1].box_name,
+ old_volume_2,
+ new_volume_2))
+
+ # Send new box lengths, recv new energies
+ [box.send_data([MessageID.CHANGE_VOLUME, box.box_length]) for box in boxes]
+ [box.recv_energy() for box in boxes]
+
+ # ---- Check move ----
+
+ # Calculate acceptance probability (See (8.3.3) in Frenkel, Smit)
+ acc = np.power(new_volume_1 / old_volume_1, boxes[0].num_part + 1) * \
+ np.power(new_volume_2 / old_volume_2, boxes[1].num_part + 1) * \
+ np.exp(- 1. / kT * (boxes[0].diff_energy + boxes[1].diff_energy))
+
+ # Draw random number
+ rand = np.random.random()
+
+ # Check if move was valid
+ if rand > acc:
+ # Reject move, restore old configuration
+ for box in boxes:
+ box.send_data([MessageID.CHANGE_VOLUME_REVERT])
+ box.box_length = box.old_box_length
+ box.energy = box.old_energy
+ return False
+
+ # Debug
+ logging.debug("Volume change accepted.")
+
+ # Only if rand < acc:
+ return True
+
+
+def perform_particle_exchange(boxes):
+ """ Remove a particle of a box chosen via choose_random_box() and add it to the other """
+
+ # Choose donor_box
+ box_num = choose_random_box(boxes)
+ donor_box = boxes[box_num]
+ recipient_box = boxes[(box_num + 1) % 2]
+
+ # Send moves, update num_part
+ donor_box.send_data([MessageID.PART_REMOVE])
+ recipient_box.send_data([MessageID.PART_ADD])
+ donor_box.num_part -= 1
+ recipient_box.num_part += 1
+ [box.recv_energy() for box in boxes]
+
+ # Debugging purposes
+ if logging.debug:
+ np.testing.assert_equal(
+ donor_box.num_part +
+ recipient_box.num_part,
+ GLOBAL_NUM_PART)
+
+ logging.debug(
+ "Exchange particle of {} to {}.".format(
+ donor_box.box_name,
+ recipient_box.box_name))
+
+ # ---- Check move ----
+
+ # This is the acceptance probability if the donor_box is chosen with equal
+ # probability
+ acc = np.exp(
+ np.log((donor_box.num_part + 1) * recipient_box.volume /
+ (recipient_box.num_part * donor_box.volume))
+ - 1. / kT * donor_box.diff_energy
+ - 1. / kT * recipient_box.diff_energy)
+
+ # Draw random number
+ rand = np.random.random()
+
+ # Check if move was valid
+ if rand > acc:
+ # Reject move, restore old configuration
+ donor_box.send_data([MessageID.PART_REMOVE_REVERT])
+ recipient_box.send_data([MessageID.PART_ADD_REVERT])
+ donor_box.num_part += 1
+ recipient_box.num_part -= 1
+ donor_box.energy = donor_box.old_energy
+ recipient_box.energy = recipient_box.old_energy
+ return False
+
+ # Debug
+ logging.debug("Particle exchange accepted.")
+
+ return True
+
+
+def choose_move():
+ """ Choose one of the trial moves with the probabilities stated at the beginning of the script """
+
+ rand = np.random.random()
+
+ # Choose move
+ if rand < MOVE_CHANCE:
+ current_move = Moves.MOVE_PARTICLE
+ elif rand < MOVE_CHANCE + VOLUME_CHANCE:
+ current_move = Moves.EXCHANGE_VOLUME
+ else:
+ current_move = Moves.EXCHANGE_PARTICLE
+
+ return current_move
+
+
+# -------------------- Start of program -------------------- #
+
+# Lists for boxes and pipes
+boxes = []
+pipes = []
+
+# Start processes
+for i in range(NUM_OF_CLIENTS):
+ pipes.append(Pipe())
+ boxes.append(Box("Box0{}".format(i),
+ pipes[i],
+ box_seeds[i],
+ init_box_length,
+ int(GLOBAL_NUM_PART / 2)))
+
+ # Start the client
+ boxes[i].client.start()
+
+ # Receive initial energy
+ boxes[i].recv_energy()
+
+logging.info("-------------------- Start of program --------------------")
+
+logging.info(
+ "Simulation parameters:\nRandom seed: {},\tWarmup steps: {},\tSteps: {},\tFilename: {}".format(
+ args.seed,
+ warmup,
+ steps,
+ filename))
+
+logging.info("Warming up.")
+
+# Warmup
+for _ in (tqdm(range(warmup)) if not args.log else range(warmup)):
+ acceptance_flag = perform_move_particle(boxes)
+ if acceptance_flag:
+ num_accepted_moves[Moves.MOVE_PARTICLE] += 1
+
+logging.info(
+ "Particle move acceptance rate during warmup: {:.2f}%".format(
+ num_accepted_moves[Moves.MOVE_PARTICLE] /
+ warmup *
+ 100))
+
+# Reset the counter for particle moves
+num_accepted_moves[Moves.MOVE_PARTICLE] = 0
+
+# List of observables we want to measure during the simulation
+densities_box01 = []
+densities_box02 = []
+
+logging.info("Starting simulation.")
+
+# Run the simulation
+for count in (tqdm(range(steps)) if not args.log else range(steps)):
+
+ current_move = choose_move()
+
+ if current_move == Moves.MOVE_PARTICLE:
+ accepted = perform_move_particle(boxes)
+ elif current_move == Moves.EXCHANGE_VOLUME:
+ accepted = perform_volume_change(boxes)
+ else:
+ accepted = perform_particle_exchange(boxes)
+
+ # Add current move to counter
+ num_moves[current_move] += 1
+ if accepted:
+ num_accepted_moves[current_move] += 1
+
+ # Append simulation observables
+ densities_box01.append(boxes[0].density)
+ densities_box02.append(boxes[1].density)
+
+ if count % check == 0:
+ validate_info(boxes)
+
+
+logging.info("Percentage of moves done:\n \
+ Particle moves, Volume changes, Particle exchanges:\n \
+ {:.2f}%, {:.2f}%, {:.2f}% \n \
+ Acceptance rates for each:\n \
+ {:.2f}%, {:.2f}%, {:.2f}%".format(
+ num_moves[Moves.MOVE_PARTICLE] / steps * 100, num_moves[Moves.EXCHANGE_VOLUME] /
+ steps * 100, num_moves[Moves.EXCHANGE_PARTICLE] / steps * 100,
+ num_accepted_moves[Moves.MOVE_PARTICLE] / num_moves[Moves.MOVE_PARTICLE] * 100, num_accepted_moves[Moves.EXCHANGE_VOLUME] /
+ num_moves[Moves.EXCHANGE_VOLUME] * 100, num_accepted_moves[Moves.EXCHANGE_PARTICLE] /
+ num_moves[Moves.EXCHANGE_PARTICLE] * 100
+))
+
+logging.info("Sending close signal to boxes.")
+[box.send_data([MessageID.END]) for box in boxes]
+
+logging.info("Saving data.")
+with gzip.open(filename + ".dat.gz", 'wb') as f:
+ pickle.dump({
+ 'args': args,
+ 'steps': steps,
+ 'temperature': kT,
+ 'densities_box01': [float(i) for i in densities_box01],
+ 'densities_box02': [float(i) for i in densities_box02],
+ }, f)
+
+logging.info("-------------------- End of program --------------------")
diff --git a/samples/gibbs_ensemble_client.py b/samples/gibbs_ensemble_client.py
deleted file mode 100644
index 8e5223cccbb..00000000000
--- a/samples/gibbs_ensemble_client.py
+++ /dev/null
@@ -1,197 +0,0 @@
-#
-# Copyright (C) 2013-2019 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 .
-#
-"""
-Client part of the Gibbs ensemble simulation. This script handles the
-simulation boxes and communicates the energies to the host. The Monte-Carlo
-part of the simulation is done by the :file:`gibbs_ensemble_socket.py` script.
-"""
-
-import espressomd
-import socket
-import numpy as np
-import pickle
-import struct
-import argparse
-
-espressomd.assert_features("LENNARD_JONES")
-
-seed = None
-init_box_l = None
-particle_number = None
-
-lj_epsilon = None
-lj_sigma = None
-lj_cutoff = None
-lj_shift = None
-
-HOST = 'localhost'
-PORT = 31415
-
-# Message identifiers
-MSG_START = 0
-MSG_END = 1
-MSG_MOVE_PART = 2
-MSG_MOVE_PART_REVERT = 21
-MSG_CHANGE_VOLUME = 3
-MSG_EXCHANGE_PART_ADD = 4
-MSG_EXCHANGE_PART_ADD_REVERT = 41
-MSG_EXCHANGE_PART_REMOVE = 5
-MSG_EXCHANGE_PART_REMOVE_REVERT = 51
-MSG_ENERGY = 6
-
-parser = argparse.ArgumentParser()
-parser.add_argument('-n', '--number-particles', type=int, nargs=1)
-parser.add_argument('-lj', '--lj-params', type=float, nargs=4)
-parser.add_argument('-s', '--seed', type=int, nargs=1)
-parser.add_argument('-bl', '--box-length', type=float, nargs=1)
-
-args = parser.parse_args()
-
-if args.number_particles:
- particle_number = args.number_particles[0]
-if args.lj_params:
- lj_epsilon = args.lj_params[0]
- lj_sigma = args.lj_params[1]
- lj_cutoff = args.lj_params[2]
- lj_shift = args.lj_params[3]
-if args.seed:
- seed = args.seed[0]
-if args.box_length:
- init_box_l = args.box_length[0]
-
-
-class Manipulator():
- '''
- Manipulator class, which has all the system changing functions.
- '''
-
- def __init__(self, system):
- self._system = system
- self._old_part_position = []
- self._old_part_idx = None
-
- def move_particle(self):
- '''Moves a particle to a new random position. The old position is saved in _old_part_position.'''
- sel = self._system.part.select(lambda p: p.id > -1).id
- self._old_part = self._system.part[np.random.choice(sel)]
- self._old_part_position = self._old_part.pos
- self._old_part.pos = np.random.rand(3) * box_l
-
- def move_particle_revert(self):
- '''Revert the last movement step.'''
- self._old_part.pos = self._old_part_position
-
- def remove_particle(self):
- '''Removes a random particle. The old position and index are saved in _old_part_position and _old_part_idx.'''
- sel = self._system.part.select(lambda p: p.id > -1).id
- self._old_part = self._system.part[np.random.choice(sel)]
- self._old_part_position = self._old_part.pos
- self._old_part.remove()
-
- def remove_particle_revert(self):
- '''Revert the last remove particle step.'''
- self._system.part.add(pos=self._old_part_position)
-
- def change_volume_and_rescale_particles(self, new_box_l):
- '''Changes the volume and rescales the particle positions.'''
- self._system.change_volume_and_rescale_particles(new_box_l)
-
- def insert_particle(self):
- '''Inserts a particle at a random position.'''
- self._system.part.add(pos=np.random.rand(3) * self._system.box_l)
-
- def remove_last_added_particle(self):
- '''Removes the last added particle (reverts insert particle).'''
- self._system.part[self._system.part.highest_particle_id].remove()
-
-
-def recv_data(socket):
- '''Receives data and return it.'''
- # The first 4 bytes encode the length of the messages received.
- buf = b''
- while len(buf) < 4:
- buf += socket.recv(4 - len(buf))
- length = struct.unpack('!I', buf)[0]
- msg = pickle.loads(socket.recv(length))
- return msg
-
-
-def send_data(data, socket):
- '''Send the data packet.'''
- # The first 4 bytes encode the length of the messages sent.
- length = struct.pack('>I', len(data))
- packet = length + data
- socket.send(packet)
-
-
-# set random seed
-np.random.seed(seed)
-
-# init socket
-socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
-socket.connect((HOST, PORT))
-msg = recv_data(socket)
-
-# init system
-box_l = init_box_l
-system = espressomd.System(box_l=[box_l, box_l, box_l])
-system.cell_system.set_n_square()
-
-system.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon=lj_epsilon,
- sigma=lj_sigma,
- cutoff=lj_cutoff,
- shift=lj_shift)
-system.time_step = 0.01
-manipulator = Manipulator(system)
-
-# places the particles
-for i in range(particle_number):
- system.part.add(pos=np.random.rand(3) * box_l, type=0)
-
-# send the initial energy
-energy = system.analysis.energy()['total']
-send_data(pickle.dumps([MSG_ENERGY, energy]), socket)
-
-while msg[0] != MSG_END:
- # receive command to execute next step
- msg = recv_data(socket)
- if msg[0] == MSG_END:
- break
- elif msg[0] == MSG_MOVE_PART:
- manipulator.move_particle()
- elif msg[0] == MSG_CHANGE_VOLUME:
- box_l = msg[1]
- manipulator.change_volume_and_rescale_particles(box_l)
- elif msg[0] == MSG_EXCHANGE_PART_ADD:
- manipulator.insert_particle()
- elif msg[0] == MSG_EXCHANGE_PART_REMOVE:
- manipulator.remove_particle()
- elif msg[0] == MSG_MOVE_PART_REVERT:
- manipulator.move_particle_revert()
- elif msg[0] == MSG_EXCHANGE_PART_ADD_REVERT:
- manipulator.remove_last_added_particle()
- elif msg[0] == MSG_EXCHANGE_PART_REMOVE_REVERT:
- manipulator.remove_particle_revert()
-
- # calculation energy and send it to the host
- energy = system.analysis.energy()['total']
- send_data(pickle.dumps([MSG_ENERGY, energy]), socket)
-
-# closing the socket
-socket.close()
diff --git a/samples/gibbs_ensemble_socket.py b/samples/gibbs_ensemble_socket.py
deleted file mode 100644
index 56a676aa3fe..00000000000
--- a/samples/gibbs_ensemble_socket.py
+++ /dev/null
@@ -1,449 +0,0 @@
-#
-# Copyright (C) 2013-2019 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 .
-#
-"""
-Simulate a Gibbs-ensemble of a Lennard-Jones fluid at a fixed temperature.
-This script does the Monte-Carlo part of the Gibbs-ensemble, however for the
-energy calculation of the systems, two instances of the :file:`gibbs_ensemble_client.py`
-script are executed, which each run a different instance of ESPResSo.
-
-The Gibbs-ensemble implemented in these scripts closely refers to chapter 8 of
-Daan Frenkel and Berend Smit 'Understanding Molecular Simulation, Second edition'.
-All equation references noted in this script can be found there.
-
-Due to the cutoff and shifting of the LJ-potential the simulated points in the
-phase diagram differ from the long range uncut LJ-potential. The phase diagram
-of the used potential can be found in 'B. Smit, J. Chem. Phys. 96 (11), 1992,
-Phase diagrams of Lennard-Jones fluids'.
-"""
-
-import socket
-import numpy as np
-import pickle
-import subprocess
-import struct
-import random
-import matplotlib.pyplot as plt
-import argparse
-
-from espressomd import assert_features
-assert_features("LENNARD_JONES")
-
-parser = argparse.ArgumentParser()
-parser.add_argument('-s', '--seed', type=int, nargs=1)
-parser.add_argument('-S', '--steps', type=int, nargs=1)
-parser.add_argument('-w', '--warm_up_steps', type=int, nargs=1)
-parser.add_argument('-E', '--espresso-executable', nargs=1)
-parser.add_argument('-C', '--client-script', nargs=1)
-parser.add_argument('-n', '--number-of-particles', type=int, nargs=1)
-parser.add_argument('-l', '--box-length', type=float, nargs=1)
-parser.add_argument('-T', '--temperature', type=float, nargs=1)
-
-args = parser.parse_args()
-
-# system parameters
-seed = None
-steps = 50000
-warm_up_steps = 1000
-# number of particles in both boxes combined
-global_num_particles = 256
-# starting box length
-init_box_l = 7.55
-# temperature
-kT = 1.15
-# maximum of the volume exchanged in the logarithmic space
-DV_MAX = 0.5
-PARTICLE_RADIUS = 0.5
-
-# LJ-parameters
-LJ_EPSILON = 1.0
-LJ_SIGMA = 2.0 * PARTICLE_RADIUS
-LJ_CUTOFF = 2.5
-LJ_SHIFT = 0.0
-
-# Monte-Carlo parameters
-INIT_MOVE_CHANCE = 0.16
-EXCHANGE_CHANCE = 0.8
-VOLUME_CHANCE = 1.0 - INIT_MOVE_CHANCE - EXCHANGE_CHANCE
-
-# socket parameters
-HOST = 'localhost'
-PORT = 31415
-NUMBER_OF_CLIENTS = 2
-
-# Message identifiers
-MSG_START = 0
-MSG_END = 1
-MSG_MOVE_PART = 2
-MSG_MOVE_PART_REVERT = 21
-MSG_CHANGE_VOLUME = 3
-MSG_EXCHANGE_PART_ADD = 4
-MSG_EXCHANGE_PART_ADD_REVERT = 41
-MSG_EXCHANGE_PART_REMOVE = 5
-MSG_EXCHANGE_PART_REMOVE_REVERT = 51
-MSG_ENERGY = 6
-
-# script locations
-espresso_executable = "../pypresso"
-client_script = "./gibbs_ensemble_client.py"
-
-if args.seed:
- seed = args.seed[0]
-if args.steps:
- steps = args.steps[0]
-if args.warm_up_steps:
- warm_up_steps = args.warm_up_steps[0]
-if args.espresso_executable:
- espresso_executable = args.espresso_executable[0]
-if args.client_script:
- client_script = args.client_script[0]
-if args.number_of_particles:
- global_num_particles = args.number_of_particles[0]
-if args.temperature:
- kT = args.temperature[0]
-if args.box_length:
- init_box_l = args.box_length[0]
-
-global_volume = 2.0 * init_box_l**3
-
-
-class Box:
- '''
- Box class, which contains the data of one system and controls the
- connection to the respective client.
- '''
- box_l = init_box_l
- box_l_old = init_box_l
- n_particles = int(global_num_particles / 2)
- energy = 0.0
- energy_corrected = 0.0
- energy_old = 1.0e100
- conn = 0
- adr = 0
-
- def recv_energy(self):
- '''Received the energy data from the client.'''
- msg = self.recv_data()
- if msg[0] == MSG_ENERGY:
- self.energy = msg[1]
- return 0
- else:
- print("ERROR during energy recv")
- return 1
-
- def recv_data(self):
- '''Received the data send from the client.'''
- # The first 4 bytes encode the length of the messages received.
- buf = b''
- while len(buf) < 4:
- buf += self.conn.recv(4 - len(buf))
- length = struct.unpack('!I', buf)[0]
- d = self.conn.recv(length)
- msg = pickle.loads(d)
- return(msg)
-
- def send_data(self, data):
- '''Send the data packet to the client.'''
- # The first 4 bytes encode the length of the messages sent.
- length = struct.pack('>I', len(data))
- packet = length + data
- self.conn.send(packet)
-
-
-def calc_tail_correction(box, lj_epsilon, lj_sigma, lj_cutoff):
- '''
- Calculates the tail correction to the energies of the box.
- '''
- # eq 3.2.5
- return 8.0 / 3.0 * np.pi * box.n_particles / box.box_l**3 * lj_epsilon * \
- lj_sigma**3 * (1.0 / 3.0 * np.power(lj_cutoff / lj_sigma, -9) -
- np.power(lj_cutoff / lj_sigma, -3))
-
-
-def calc_shift_correction(box, lj_epsilon, lj_cutoff, lj_shift):
- '''
- Calculates the shift correction to the energies of the box.
- '''
- # difference in the potential integrated from 0 to cutoff distance
- return -8.0 / 3.0 * np.pi * box.n_particles / box.box_l**3 * \
- lj_epsilon * np.power(lj_cutoff, 3) * 4.0 * lj_shift
-
-
-def move_particle(boxes, global_num_particles):
- '''
- Tries a displacement move and stores the new energy in the corresponding box
- '''
- if random.randint(0, global_num_particles - 1) < boxes[0].n_particles:
- rand_box = 0
- else:
- rand_box = 1
- boxes[rand_box].send_data(pickle.dumps([MSG_MOVE_PART, 0]))
- boxes[rand_box].recv_energy()
- return rand_box
-
-
-def exchange_volume(boxes, global_volume):
- '''
- Tries a volume exchange move and stores the new energy in the boxes
- '''
- boxes[0].box_l_old = boxes[0].box_l
- boxes[1].box_l_old = boxes[1].box_l
-
- # calculate the exchanged volume
- rand_box = random.randint(0, NUMBER_OF_CLIENTS - 1)
- vol2 = global_volume - boxes[rand_box].box_l**3
- lnvn = np.log(boxes[rand_box].box_l**3 / vol2) + \
- (random.random() - 0.5) * DV_MAX
- vol1 = global_volume * np.exp(lnvn) / (1 + np.exp(lnvn))
- vol2 = global_volume - vol1
- boxes[rand_box].box_l = np.cbrt(vol1)
- boxes[rand_box].send_data(pickle.dumps(
- [MSG_CHANGE_VOLUME, boxes[rand_box].box_l]))
-
- boxes[(rand_box + 1) % 2].box_l = np.cbrt(vol2)
- boxes[(rand_box + 1) % 2].send_data(pickle.dumps([MSG_CHANGE_VOLUME,
- boxes[(rand_box + 1) % 2].box_l]))
-
- boxes[rand_box].recv_energy()
- boxes[(rand_box + 1) % 2].recv_energy()
- return rand_box
-
-
-def exchange_particle(boxes):
- '''
- Tries a particle exchange move and stores the new energy in the boxes
- '''
- rand_box = random.randint(0, 1)
- if boxes[rand_box].n_particles == 0:
- rand_box = (rand_box + 1) % 2
- boxes[rand_box].n_particles -= 1
- boxes[(rand_box + 1) % 2].n_particles += 1
- boxes[rand_box].send_data(pickle.dumps([MSG_EXCHANGE_PART_REMOVE, 0]))
- boxes[(rand_box + 1) %
- 2].send_data(pickle.dumps([MSG_EXCHANGE_PART_ADD, 0]))
-
- boxes[rand_box].recv_energy()
- boxes[(rand_box + 1) % 2].recv_energy()
- return rand_box
-
-
-def check_make_move(boxes, inner_potential, rand_box):
- '''
- Returns whether the last displacement move was valid or not and reverts the changes
- if it was invalid.
- '''
- if random.random() > inner_potential:
- boxes[rand_box].send_data(pickle.dumps([MSG_MOVE_PART_REVERT, 0]))
- boxes[rand_box].recv_energy()
- return False
- return True
-
-
-def check_exchange_volume(boxes, inner_potential):
- '''
- Returns whether the last volume exchange move was valid or not and reverts the changes
- if it was invalid.
- '''
- volume_factor = \
- (boxes[0].box_l**3 / boxes[0].box_l_old**3)**(boxes[0].n_particles + 1) * \
- (boxes[1].box_l**3 / boxes[1].box_l_old **
- 3)**(boxes[1].n_particles + 1)
- if random.random() > volume_factor * inner_potential:
- boxes[0].send_data(pickle.dumps(
- [MSG_CHANGE_VOLUME, boxes[0].box_l_old]))
- boxes[0].box_l = boxes[0].box_l_old
- boxes[1].send_data(pickle.dumps(
- [MSG_CHANGE_VOLUME, boxes[1].box_l_old]))
- boxes[1].box_l = boxes[1].box_l_old
-
- boxes[0].recv_energy()
- boxes[1].recv_energy()
- return False
- return True
-
-
-def check_exchange_particle(boxes, inner_potential, rand_box):
- '''
- Returns whether the last particle exchange move was valid or not and reverts the changes
- if it was invalid.
- '''
- exchange_factor = (boxes[rand_box].n_particles) / \
- (boxes[(rand_box + 1) % 2].n_particles + 1.0) * \
- boxes[(rand_box + 1) % 2].box_l**3 / boxes[rand_box].box_l**3
- if random.random() > exchange_factor * inner_potential:
- boxes[rand_box].n_particles += 1
- boxes[(rand_box + 1) % 2].n_particles -= 1
- boxes[rand_box].send_data(pickle.dumps(
- [MSG_EXCHANGE_PART_REMOVE_REVERT, 0]))
- boxes[(rand_box + 1) %
- 2].send_data(pickle.dumps([MSG_EXCHANGE_PART_ADD_REVERT, 0]))
-
- boxes[rand_box].recv_energy()
- boxes[(rand_box + 1) % 2].recv_energy()
- return False
- return True
-
-
-# init socket
-s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
-s.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
-s.bind((HOST, PORT))
-s.listen(NUMBER_OF_CLIENTS)
-
-boxes = []
-
-random.seed(seed)
-
-# start clients and connections
-for i in range(NUMBER_OF_CLIENTS):
- boxes.append(Box())
- lj_arguments = ["-lj", str(LJ_EPSILON), str(LJ_SIGMA), str(LJ_CUTOFF),
- str(LJ_SHIFT)]
- arguments = ["-n", str(boxes[i].n_particles), "-s",
- str(random.randint(0, np.iinfo(np.int32).max)), "-bl",
- str(boxes[i].box_l)]
- subprocess.Popen([espresso_executable] + [client_script] +
- arguments + lj_arguments)
-
- boxes[i].conn, boxes[i].adr = s.accept()
- boxes[i].send_data(pickle.dumps([MSG_START, 0]))
-
- boxes[i].recv_energy()
- boxes[i].energy_old = boxes[i].energy + \
- calc_tail_correction(boxes[i], LJ_EPSILON, LJ_SIGMA, LJ_CUTOFF) + \
- calc_shift_correction(boxes[i], LJ_EPSILON, LJ_CUTOFF, LJ_SHIFT)
-
-# observables
-densities = [[], []]
-list_indices = []
-accepted_steps = 0
-accepted_steps_move = 1
-accepted_steps_volume = 1
-accepted_steps_exchange = 1
-move_steps_tried = 1
-volume_steps_tried = 1
-exchange_steps_tried = 1
-
-# rand_box defines on which box the move acts on.
-rand_box = 0
-
-# only do displacements during the warm up
-move_chance = 1.0
-# Monte-Carlo loop
-for i in range(steps):
- print("\rIntegrating: {0:3.0f}%".format(
- 100 * i / steps), end='', flush=True)
-
- # warm up ends after 1000 steps
- if i == warm_up_steps:
- move_chance = INIT_MOVE_CHANCE
-
- # rand defines which move to make.
- rand = random.random()
-
- # choose step and send the command to execute to the clients, then receive
- # the energies.
- if rand <= move_chance:
- rand_box = move_particle(boxes, global_num_particles)
- move_steps_tried += 1
-
- elif rand <= move_chance + VOLUME_CHANCE:
- rand_box = exchange_volume(boxes, global_volume)
- volume_steps_tried += 1
-
- else:
- rand_box = exchange_particle(boxes)
- exchange_steps_tried += 1
-
- # calculate the correction energies of the lj tail and shift.
- if rand <= move_chance:
- boxes[rand_box].energy_corrected = boxes[rand_box].energy + \
- calc_tail_correction(boxes[rand_box], LJ_EPSILON, LJ_SIGMA, LJ_CUTOFF) + \
- calc_shift_correction(
- boxes[rand_box],
- LJ_EPSILON,
- LJ_CUTOFF,
- LJ_SHIFT)
- boxes[(rand_box + 1) % 2].energy_corrected = \
- boxes[(rand_box + 1) % 2].energy_old
-
- else:
- boxes[0].energy_corrected = boxes[0].energy + \
- calc_tail_correction(boxes[0], LJ_EPSILON, LJ_SIGMA, LJ_CUTOFF) + \
- calc_shift_correction(boxes[0], LJ_EPSILON, LJ_CUTOFF, LJ_SHIFT)
- boxes[1].energy_corrected = boxes[1].energy + \
- calc_tail_correction(boxes[1], LJ_EPSILON, LJ_SIGMA, LJ_CUTOFF) + \
- calc_shift_correction(boxes[1], LJ_EPSILON, LJ_CUTOFF, LJ_SHIFT)
-
- # test if the move will be accepted and undo the last step if it was not
- # accepted.
- delta_energy = boxes[0].energy_corrected + boxes[1].energy_corrected - \
- boxes[0].energy_old - boxes[1].energy_old
- # limitation to delta_energy to circumvent calculating the exponential of
- # too large inner potentials, which could cause some error messages.
- if delta_energy < -10.0 * kT:
- delta_energy = -10.0 * kT
-
- inner_potential = np.exp(-1.0 / kT * delta_energy)
- accepted = True
-
- if rand <= move_chance:
- accepted = check_make_move(boxes, inner_potential, rand_box)
- elif rand <= move_chance + VOLUME_CHANCE:
- accepted = check_exchange_volume(boxes, inner_potential)
- else:
- accepted = check_exchange_particle(boxes, inner_potential, rand_box)
-
- if accepted:
- # keep the changes.
- boxes[0].energy_old = boxes[0].energy_corrected
- boxes[1].energy_old = boxes[1].energy_corrected
- accepted_steps += 1
- if rand <= move_chance:
- accepted_steps_move += 1
- elif rand <= move_chance + VOLUME_CHANCE:
- accepted_steps_volume += 1
- else:
- accepted_steps_exchange += 1
- densities[0].append(boxes[0].n_particles / boxes[0].box_l**3)
- densities[1].append(boxes[1].n_particles / boxes[1].box_l**3)
-
- list_indices.append(i)
-
-
-print("Acceptance rate global:\t {}".format(accepted_steps / float(steps)))
-print("Acceptance rate moving:\t {}".format(
- accepted_steps_move / float(move_steps_tried)))
-print("Acceptance rate volume exchange:\t {}".format(
- accepted_steps_volume / float(volume_steps_tried)))
-print("Acceptance rate particle exchange:\t {}".format(
- accepted_steps_exchange / float(exchange_steps_tried)))
-
-
-plt.figure()
-plt.ylabel('density')
-plt.xlabel('number of steps')
-plt.plot(list_indices[100:], densities[0][100:], label="box 0")
-plt.plot(list_indices[100:], densities[1][100:], label="box 1")
-plt.legend()
-plt.show()
-
-# closing the socket
-for i in range(NUMBER_OF_CLIENTS):
- boxes[i].send_data(pickle.dumps([MSG_END, 0]))
-s.close()
diff --git a/testsuite/scripts/samples/test_gibbs_ensemble.py b/testsuite/scripts/samples/test_gibbs_ensemble.py
index 423e22b3dbb..665bbd7d7ab 100644
--- a/testsuite/scripts/samples/test_gibbs_ensemble.py
+++ b/testsuite/scripts/samples/test_gibbs_ensemble.py
@@ -19,9 +19,8 @@
import importlib_wrapper
sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
- "@SAMPLES_DIR@/gibbs_ensemble_socket.py",
- cmd_arguments=["-S", "1000", "-C", "@SAMPLES_DIR@/gibbs_ensemble_client.py",
- "-E", "@CMAKE_BINARY_DIR@/pypresso"])
+ "@SAMPLES_DIR@/gibbs_ensemble/run_sim.py",
+ cmd_arguments=["0.7", "-N", "256", "--warmup", "1000", "--steps", "1000", "--seed", "10"])
@skipIfMissingFeatures