Gibbs ensemble simulations¶
Introduction¶
In this workshop, you will use grand canonical, multicanonical, and Gibbs ensemble Monte Carlo simulations to investigate liquid-vapour coexistence in the Lennard-Jones system.
The instructions for the workshop are contained in two notebooks: MC-MUCA.ipynb and MC-Gibbs.ipynb. They do not depend on each other: you can tackle them in either order, or both at once.
Colab Setup (optional)¶
The next cells should only be needed if you use Google Colab, so we test the environment for this. Commands may work in other environments too but are not tested. If you need to try, just replace the following by SETUP = True.
[ ]:
import sys
SETUP = "google.colab" in sys.modules
if SETUP:
print('Fetching files for Google Colab')
[ ]:
if SETUP:
! pip install h5py data_tutorials weas_widget ase meson ninja
! apt install gfortran libhdf5-dev
[ ]:
if SETUP:
from data_tutorials.data import get_data
# Get the python modules and top-level meson.build needed for this tutorial
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/",
filename=["eos_ljcs.py", "hdf5_module.py","meson.build"],
folder=".",
)
# Get starting configurations in DATA subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/DATA/",
filename=["mc_gibbs_liq.xyz","mc_gibbs_vap.xyz"],
folder="DATA/",
)
# Get the Fortran source files, and meson.build file, in the src subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/src/",
filename=["meson.build","config_io_module.f90","hdf5_module.f90","maths_module.f90","mc_gibbs.f90","mc_muca.f90","mc_module.f90","potential_module.f90"],
folder="src/",
)
Preliminaries¶
Start by importing some useful Python modules and functions. Here are links to information about NumPy, Matplotlib, pathlib, and shutil, in case you need any details. We shall use ase, and the weas-widget tool to visualize atomic configurations. The remaining modules are supplied as part of this workshop. We also apply some Matplotlib settings.
[ ]:
# Standard and external modules
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
from shutil import copy
from ase.atoms import Atoms
from ase.io import read
from ase.data import vdw_radii, atomic_numbers
from weas_widget import WeasWidget
# Modules supplied in this workshop
from hdf5_module import read_hdf5_file
from eos_ljcs import eos, rho_coex
# Settings for Matplotlib plots
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis') # Default colour map
plt.rc('figure',figsize=(8,5)) # Default figure dimensions in inches
plt.rc('figure.constrained_layout',use=True) # Default figure layout option
plt.rc('legend',frameon=True,framealpha=1.0) # Default appearance of figure legend
plt.rc('hist',bins=100) # Default number of bins to use in histograms
# Setting for weas-widget viewer
weas_shape = {"width":"800px","height":"600px"} # Default shape of weas widget
The system simulated here is the same as that in the accompanying MC-MUCA.ipynb notebook: the cut-and-shifted Lennard-Jones potential. A reminder: the Lennard-Jones potential is \begin{equation*}
u_{\text{LJ}}(r) = 4\varepsilon \left[\left( \frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6 \right]
\end{equation*} where \(r\) is the distance between the atoms, \(\varepsilon\) is an energy characterizing the strength of the interaction, and \(\sigma\) is a length scale that characterizes the size of the atoms. Reduced units are defined so that \(\varepsilon=1\) and \(\sigma=1\), and in addition Boltzmann’s constant is taken to be unity \(k_{\text{B}}=1\). The cut-and-shifted Lennard-Jones potential used in the simulation program is defined by
:nbsphinx-math:`begin{equation*}
- u(r) =
begin{cases} u_{text{LJ}}(r) - u_{text{LJ}}(r_{text{cut}}) & r leq r_{text{cut}} \ 0 & r> r_{text{cut}} end{cases}
end{equation*}` where the cut-off distance is taken to be \(r_{\text{cut}}= 2.5\sigma\).
For the chosen potential, in reduced units, the critical point is at \(T_{\mathrm{c}}=1.0779\), \(\rho_{\mathrm{c}}=0.3190\). The temperature of interest is \(T=0.95\), and this is the default value in the simulation codes. At this temperature, liquid and vapour phases may coexist. The coexistence pressure \(P_{\text{coex}}\approx 0.045\) and chemical potential \(\mu_{\text{coex}}\approx -3.14\), as well as the coexisting densities \(\rho_{\text{liq}}\approx 0.622\) and \(\rho_{\text{vap}}\approx 0.0665\), have been studied by simulation (see e.g. J Vrabec et al., Molec Phys, 104, 1509 (2006)). Here we use the values predicted by the fitted equation of state (EOS) of M Thol et al., Int J Thermophys, 36, 25 (2015). The aim of this exercise is to simulate coexistence without prior knowledge of these quantities, and obtain estimates of them.
[ ]:
# Store the fitted EOS coexistence values for this temperature, for later use
T_coex = 0.95
rho_vap, rho_liq = rho_coex(T_coex)
# Reassure ourselves that they actually satisfy the coexistence conditions
eos_vap = eos ( T_coex, rho_vap )
eos_liq = eos ( T_coex, rho_liq )
print(' vapour liquid')
print(f"rho = {rho_vap :10.7f} {rho_liq :10.7f}")
print(f"P = {eos_vap['P'] :10.7f} {eos_liq['P'] :10.7f}")
print(f"mu = {eos_vap['mu']:10.7f} {eos_liq['mu']:10.7f}")
P_coex = eos_liq['P']
mu_coex = eos_liq['mu']
Gibbs ensemble Monte Carlo simulations¶
The following two cells will build both the programs for this workshop, mc_muca and mc_gibbs, in the build directory. Even if you already did this, it will do no harm to run the cells again.
[ ]:
!meson setup build
[ ]:
!meson compile -C build
The next cell runs the Gibbs ensemble program, after copying the initial configurations from the supplied files, into the current working directory. While it is running, read through the program description below, which refers to the important program files.
[ ]:
this_dir = Path.cwd() # Current working directory
data_dir = this_dir / 'DATA' # Data directory
old_conf_1 = data_dir / 'mc_gibbs_1.xyz' # Supplied configuration file 1
new_conf_1 = copy(old_conf_1,this_dir) # Copy supplied file into working directory
old_conf_2 = data_dir / 'mc_gibbs_2.xyz' # Supplied configuration file 2
new_conf_2 = copy(old_conf_2,this_dir) # Copy supplied file into working directory
!echo '&nml /' | build/mc_gibbs
The program mc_gibbs simulates two systems (one vapour, one liquid) at the same time. The simulations keep the total number of atoms \(N_1+N_2\) and the total volume \(V_1+V_2\) fixed.
All the source files are in the src subdirectory. The main program mc_gibbs.f90 takes its run parameters from standard input using a Fortran namelist, making it easy to specify them through a key=value mechanism, while allowing the unspecified parameters to take default values which are built into the program. You should see that the default run length is 20000 steps. This, and the other default values, should be enough for our purposes. Each step consists of
A number of attempted single-atom moves, equal to the number of particles in each system.
A number,
nswap=20by default, of attempted particle exchanges (either way) between the systemsAn attempted volume exchange between the systems.
The routines that actually perform these last two types of move are n_swap and v_swap in mc_module.f90. If you have any questions about the code, by all means ask! The program simply outputs the cumulative move acceptance rates at (increasing) intervals, to confirm that the program is running. Values of all the quantities of interest are stored at each step, and output to a file mc_gibbs.hdf5 at the end of the run, for analysis in the following cells.
Once you have finished looking through these files you may close them.
Typical starting configurations are provided in the files DATA/mc_gibbs_1.xyz and DATA/mc_gibbs_2.xyz, which we’re referring to as old_conf_1 and old_conf_2. The densities are approximately equal to the coexistence values, but it is not essential that they be exactly right. The files are in extended XYZ format, and look like this:
n
Lattice=[box(1) 0.0 0.0 0.0 box(2) 0.0 0.0 0.0 box(3)] Properties=species:S:1:pos:R:3
X r(1,1) r(2,1) r(3,1)
X r(1,2) r(2,2) r(3,2)
X r(1,3) r(2,3) r(3,3)
: : :
X r(1,n) r(2,n) r(3,n)
The first line gives the number of atoms, the second line gives the box dimensions and details of the columns in the rest of the file, and the subsequent lines give the species and (x,y,z)-coordinates of each atom. In this format, the Lattice entry on the second line gives the three vectors \([\mathbf{a}\;\mathbf{b}\;\mathbf{c}]\) of a general triclinic periodic cell. Our cells are orthorhombic (actually cubic here), so each vector points along one of the Cartesian axes, and the nonzero
elements are the corresponding box lengths. Everything is in LJ reduced units. We use X as the atomic species for each atom: in realistic simulations, a chemical element symbol would appear instead, and lengths would most commonly be given in Angstroms.
The following cell should allow you to visualize old_conf_1, converted into units appropriate for argon atoms, using weas-widget. As usual, you may use the widget controls Atoms/Radius Type/VDW to give a space-filling view. If the viewer size doesn’t suit your display, you can adjust weas_shape (defined in Preliminaries above), or try the widget’s fullscreen button ⛶ (may not work in Colab). Feel free to change old_conf_1 to
old_conf_2 in the cell below, and (if the program has finished) also look at the final configurations new_conf_1 and new_conf_2, which are in the current working directory.
[ ]:
atoms=read(old_conf_1)
print('Cell dimensions (reduced units)',atoms.cell.lengths())
print('Cell angles (degrees) ',atoms.cell.angles())
atoms.center()
atoms.numbers[:] = atomic_numbers['Ar'] # Visualize as argon atoms
sigma = vdw_radii[atomic_numbers['Ar']]*2 # Argon diameter in Angstroms
atoms.set_cell(atoms.get_cell() * sigma, scale_atoms=True) # Scale all positions and box sides
viewer=WeasWidget(from_ase=atoms,viewerStyle=weas_shape)
viewer
Wait until the run has completed before proceeding. At this temperature, it is most likely that the two systems have stayed in their original phase (liquid or vapour) throughout. If this is the case, you should see that the acceptance ratios for single-particle moves are significantly different (for simplicity, the maximum particle displacement dr_max is the same in both systems). The acceptance ratio for volume displacements should be moderately high, but for particle swaps it is rather
low. For this reason we attempt more than one swap per step.
In the following, we must bear in mind that any averages calculated in a single system may be invalid, because of the possibility that the phases might have swapped. This is more likely at higher temperatures, closer to the critical temperature. We shall do our analysis in terms of histograms of the calculated properties, computed over each system separately, assuming that no swaps have happened. If this were not the case, it would be possible to combine the data into a single set and analyse it.
The next cells open the HDF5 file, and read the attributes (simulation parameters) and datasets, in a way that is hopefully familiar from earlier workshops.
[ ]:
attr, data = read_hdf5_file('mc_gibbs.hdf5')
[ ]:
print(attr['Title'].astype(str))
print('Number of steps',attr['nstep'])
Ntot = attr['N1+N2']
Vtot = attr['V1+V2']
T = attr['T']
print(f'Total N1+N2 = {Ntot:10d}')
print(f'Total V1+V2 = {Vtot:10.4f}')
print(f'Temperature T = {T:10.4f}')
Now we consider the step-by-step datasets. Notice that the program produces separate values for each system, so these arrays are of dimension (nstep,2).
[ ]:
N = data['N']
V = data['V']
W = data['W']
Z = data['Z']
print('Shape of N dataset ',N.shape)
It will be instructive to look at histograms for each system separately, starting with the density: \(p(\rho)\). This should give a clue as to whether the phases swapped during the run or not. We also compute averages, for each system, which may match the expected literature values given at the top of this worksheet. We could also combine the data from the two systems before histogramming, which would be a better approach if some phase swapping occurred.
[ ]:
rho = N / V
rhoavg = np.mean(rho,axis=0) # Two simulation averages, one in each system
print(' System 1 System 2')
print('Average density {:10.4f}{:10.4f}'.format(*rhoavg))
print('Expected density {:10.4f}{:10.4f}'.format(rho_vap,rho_liq))
fig, ax = plt.subplots()
ax.set_xlabel(r'$\rho$')
ax.set_ylabel(r'$p(\rho)$')
ax.hist(rho[:,0],density=True,label='System 1')
ax.hist(rho[:,1],density=True,label='System 2')
ylim=ax.get_ylim()
ax.vlines([rho_vap,rho_liq],*ylim,linestyles='dotted',colors='k',label='Expected')
ax.legend(loc='upper center');
All being well, there should be two distinct peaks around the expected values of density. The Gibbs ensemble should automatically adjust both systems to give coexistence.
Optional: pressure calculation¶
The coexistence pressure is not specified. In principle the run should calculate it, and we have virial datasets W available from both boxes.
The following cell computes the run-averaged pressures in both systems, through the usual virial expression. Are they approximately equal to each other? Do they match the expected fitted EOS value (given at the start of this notebook)?
[ ]:
P = rho*T + W/(3*V)
Pavg = P.mean(axis=0) # Two simulation averages, one in each system
print(f'Average pressure in system 1 P = {Pavg[0]:10.4f}')
print(f'Average pressure in system 2 P = {Pavg[1]:10.4f}')
print(f'Expected coexistence pressure P = {P_coex:10.4f}')
Let us construct histograms as we did before.
[ ]:
fig, ax = plt.subplots()
ax.set_xlabel(r'$P$')
ax.set_ylabel(r'$p(P)$')
ax.hist(P[:,0],density=True,label='System 1')
ax.hist(P[:,1],density=True,label='System 2',alpha=0.8)
ax.axvline(x=P_coex,linestyle='dashed',color='k',label='Expected')
ax.legend();
The distributions are different in the vapour (narrower) and liquid (broader), but both should be centred on (approximately) the same value.
Optional: chemical potential calculation¶
Chemical potentials may be estimated through Widom test particle insertion in both systems. The Z dataset actually gives an estimate of \(\exp(-\beta\mu)\) where \(\beta=1/k_{\mathrm{B}}T\) and the chemical potential \(\mu\) is defined in a convention where the thermal de Broglie wavelength \(\Lambda=1\). So, this is the inverse of the activity \(z=\exp(\beta\mu)\).
[ ]:
zavg = 1/Z.mean(axis=0) # Two simulation averages, one in each system
print(f'System 1 average activity z = {zavg[0]:10.4f}')
print(f'System 2 average activity z = {zavg[1]:10.4f}')
z_coex = np.exp(mu_coex/T)
print(f'Expected coexistence activity z = {z_coex:10.4f}')
muavg = T*np.log(zavg) # two simulation estimates of chemical potential, one in each system
print(f'System 1 estimated chem pot mu = {muavg[0]:10.4f}')
print(f'System 2 estimated chem pot mu = {muavg[1]:10.4f}')
print(f'Expected coexistence chem pot mu = {mu_coex:10.4f}')
Hopefully, some measure of agreement should be seen.
However, the underlying histograms are quite different from what we have seen before. It is simplest to look directly at the distributions of \(Z\) values returned by the simulation, one in each box.
[ ]:
fig, ax = plt.subplots()
ax.set_xlabel(r'$Z$')
ax.set_ylabel(r'$p(Z)$')
ax.hist(Z[:,0],density=True,label='System 1')
ax.hist(Z[:,1],density=True,label='System 2',alpha=0.5)
ax.axvline(x=1/z_coex,linestyle='dashed',color='k',label='Expected')
ax.legend();
You might like to plot histograms for the two boxes separately (assuming that one box has remained liquid, and the other vapour, throughout). You could also experiment with the axis limits using ax.set_xlim(...) etc.
Although the average values are (or should be) close to the expected value, this histogram illustrates the different character of sampling for test particle insertion, compared to the usual variables such as \(\rho\), \(P\), as seen earlier in the worksheet.
In the liquid most insertions involve large positive potential energies (overlaps) and hence values close to zero of the quantity being averaged, which is roughly \(\exp(-\beta\Delta U)/\rho\). However, a lucky insertion in a “hole” with, say, 11 or 12 neighbours, might generate \(\beta\Delta U\approx -11\) or \(-12\) (in reduced units), and values of \(\exp(-\beta\Delta U)\) of order \(10^5\).
In the vapour, most insertions will have \(\Delta U\approx 0\), some will involve overlaps, and a small fraction will generate negative values of \(\Delta U\) due to attractive interactions with one or more atoms. The distribution of values of \(\exp(-\beta\Delta U)/\rho\) will be less dramatic than in the liquid, but still quite extreme compared with those that we are used to for \(\rho\), \(P\), etc.
This concludes the notebook, i.e. the Gibbs simulation part of this workshop.
[ ]: