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.
Setup (optional)¶
This step is needed only if you run on google colab. Instructions may work on other systems but are not tested.
[ ]:
! pip install h5py data_tutorials weas_widget ase
! apt install gfortran libhdf5-dev
Now bring the data needed for the exercises.
[ ]:
from data_tutorials.data import get_data
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/",
filename=["dat_to_ase.py", "eos_lj.py", "hdf5_module.py"],
folder=".",
)
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/",
filename=["Makefile","config_io_module.f90","hdf5_module.f90","maths_module.f90","mc_gibbs.f90","mc_gibbs_module.f90","mc_module.f90","potential_module.f90"],
folder=".",
)
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/",
filename=["config_one_old.dat", "config_two_old.dat"],
folder=".",
)
Preliminaries¶
Start by importing some useful Python modules.
[ ]:
import matplotlib.pyplot as plt
import numpy as np
from hdf5_module import read_file
from dat_to_ase import argon
from weas_widget import WeasWidget
from eos_lj import eos
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis')
plt.rc('legend',frameon=True,framealpha=1.0)
plt.rc('hist',bins=100) # Default number of bins to use in histograms
The system simulated here is the same as that in the accompanying MC-MUCA.ipynb notebook: the cut-and-shifted Lennard-Jones potential with cut-off \(r_{\mathrm{cut}}= 2.5\sigma\). Again, reduced units are employed, so the Lennard-Jones parameters are \(\varepsilon=1\) and \(\sigma=1\), and in addition Boltzmann’s constant is taken to be unity \(k_{\mathrm{B}}=1\). The critical point is at \(T_{\mathrm{c}}=1.0779\), \(\rho_{\mathrm{c}}=0.3190\).
Just as in the accompanying notebook, the temperature of interest here is \(T=0.95\), at which liquid and vapour phases may coexist with densities \(\rho_{\mathrm{liq}}\approx 0.622\) and \(\rho_{\mathrm{vap}}\approx 0.0665\) (see J Vrabec, GK Kedia, G Fuchs, H Hasse, Molec Phys, 104, 1509 (2006)). At coexistence, the pressure is \(P_{\text{coex}}\approx 0.045\) and the chemical potential \(\mu_{\text{coex}}\approx-3.14\).
[ ]:
# Store the literature coexistence values for this temperature, for later use
rho_vap = 0.0665
rho_liq = 0.622
P_coex = 0.045
mu_coex = -3.14
Gibbs ensemble Monte Carlo simulations¶
Run the following cells to make the program and start the simulation. While it is running, read through the program description below, which refers to the important program files.
[ ]:
!make mc_gibbs
[ ]:
!echo '&nml /' | ./mc_gibbs
The simulations keep the total number of atoms \(N_1+N_2\) and the total volume \(V_1+V_2\) fixed.
The program files of most interest here are mc_gibbs.f90 and mc_gibbs_module.f90. Open both of them and take a look.
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 in mc_gibbs_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.
The program mc_gibbs simulates two systems (one liquid, one vapour) at the same time. Typical starting configurations are provided in the files config_one_old.dat and config_two_old.dat. The format of the configuration files used here is as follows:
n
xbox ybox zbox
x1 y1 z1
x2 y2 z2
x3 y3 z3
: : :
xn yn zn
where the first line gives the number of atoms, the second line gives the box dimensions (in this case cubic boxes are employed) and the subsequent lines give the coordinates of each atom. The following cell should let you visualize config_one_old.dat. Feel free to edit the filename, to look at config_two_old.dat and (if the program has finished) the final configurations config_one.dat and config_two.dat.
[ ]:
atoms=argon('config_one_old.dat')
v=WeasWidget()
v.from_ase(atoms)
v
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, or alternatively just re-do the run.
The next cells open the HDF5 file, and read the simulation parameters (attributes) and datasets, in a way that is hopefully familiar from earlier workshops.
[ ]:
params, data = read_file('mc_gibbs.hdf5')
[ ]:
print(params['Title'].astype(str))
print('Number of steps',params['nstep'])
Ntot = params['N1+N2']
Vtot = params['V1+V2']
T = params['T']
print(f'Total N1+N2 = {Ntot:10d}')
print(f'Total V1+V2 = {Vtot:10.4f}')
print(f'Temperature T = {T:10.4f}')
print('Shape of N dataset ',data['N'].shape)
print('Shape of V dataset ',data['V'].shape)
Notice that the program produces separate step-by-step values of \(N\) and \(V\) for each system. (The same is true for the other datasets as well). It will be instructive to look at histograms for each system separately, starting with the density: \(\mathcal{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 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 = data['N'] / data['V']
rhoavg = np.mean(rho,axis=0) # Two simulation averages, one in each system
print(f'Average density in system 1 = {rhoavg[0]:10.4f}')
print(f'Expected vapour density = {rho_vap:10.4f}')
print(f'Average density in system 2 = {rhoavg[1]:10.4f}')
print(f'Expected liquid density = {rho_liq:10.4f}')
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$\rho$')
ax.set_ylabel(r'$\mathcal{P}(\rho)$')
ax.hist(rho[:,0],density=True,label='Box 1')
ax.hist(rho[:,1],density=True,label='Box 2')
ax.axvline(x=rhoavg[0],c='C2',label='Box 1 average')
ax.axvline(x=rhoavg[1],c='C3',label='Box 2 average')
ax.axvline(x=rho_vap,ls='dashed',c='C4',label='vapour')
ax.axvline(x=rho_liq,ls='dashed',c='C5',label='liquid')
ax.legend()
plt.tight_layout()
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 data['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 literature value (given at the start of this notebook)? And, out of interest, is the fitted EOS approximately correct?
[ ]:
P = rho*T + data['W']/(3*data['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}')
P_liq = eos(density=rho_liq,temperature=0.95)['P']
P_vap = eos(density=rho_vap,temperature=0.95)['P']
print(f'EOS liquid state pressure P = {P_liq:10.4f}')
print(f'EOS vapour state pressure P = {P_vap:10.4f}')
Let us construct histograms as we did before.
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$P$')
ax.set_ylabel(r'$\mathcal{P}(P)$')
ax.hist(P[:,0],density=True,label='Box 1')
ax.hist(P[:,1],density=True,label='Box 2',alpha=0.8)
ax.axvline(x=P_coex,ls='dashed',c='C2',label='coexistence')
ax.legend()
plt.tight_layout()
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 data['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)\).
[ ]:
Z = data['Z']
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}')
z_liq = eos(density=rho_liq,temperature=0.95)['z']
z_vap = eos(density=rho_vap,temperature=0.95)['z']
print(f'EOS liquid state activity z = {z_liq:10.4f}')
print(f'EOS vapour state activity z = {z_vap: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}')
mu_liq = eos(density=rho_liq,temperature=0.95)['mu']
mu_vap = eos(density=rho_vap,temperature=0.95)['mu']
print(f'EOS liquid state chem pot mu = {mu_liq:10.4f}')
print(f'EOS vapour state chem pot mu = {mu_vap: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(figsize=(8,5))
ax.set_xlabel(r'$Z$')
ax.set_ylabel(r'$\mathcal{P}(Z)$')
ax.hist(Z[:,0],density=True,label='Box 1')
ax.hist(Z[:,1],density=True,label='Box 2',alpha=0.5)
ax.axvline(x=1/z_coex,ls='dashed',c='C2',label='coexistence')
ax.legend()
plt.tight_layout()
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.
[ ]: