Gibbs Ensemble scripts

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).

samples/gibbs_ensemble/client_system.py

@@ -0,0 +1,71 @@
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+"""
+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))
+
+
+class Gibbs_Client(gibbs_ensemble.Client):
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    @property
+    def energy(self):
+        """ Energy handle (energy corrections can be added here) """
+        return super().energy
+
+    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

samples/gibbs_ensemble/create_fits.py

@@ -0,0 +1,229 @@
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+"""
+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 = 3000
+
+# 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))
+    with gzip.open(f, 'rb') as _f:
+        data = pickle.load(_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()