Multicanonical 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 this 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. Be aware that the kernel will restart. 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 ipympl weas_widget ase meson ninja
! apt install gfortran libhdf5-dev
get_ipython().kernel.do_shutdown(restart=True)
[ ]:
if SETUP:
# Allow widgets to be shown in Google Colab
from google.colab import output
output.enable_custom_widget_manager()
[ ]:
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 specimen HDF5 files in DATA subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/DATA/",
filename=["0.hdf5","1.hdf5","2.hdf5","3.hdf5","4.hdf5","5.hdf5","long.hdf5"],
folder="DATA/",
)
# Get starting configurations and specimen extxyz files in DATA subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/DATA/",
filename=["0.xyz","1.xyz","2.xyz","3.xyz","4.xyz","5.xyz","mc_muca_vap.xyz","mc_muca_liq.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_module.f90","mc_muca.f90","mc_gibbs.f90","potential_module.f90"],
folder="src/",
)
Preliminaries¶
Start by importing some useful Python modules and functions. Here are links to information about NumPy, SciPy, Matplotlib, ipywidgets, 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.
[ ]:
%matplotlib ipympl
# Standard and external modules
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
from pathlib import Path
from shutil import copy
from ipywidgets import FloatSlider, VBox
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,4)) # 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-Gibbs.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_{\text{c}}=1.0779\), \(\rho_{\text{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 default value of \(\mu\) for the preliminary grand canonical MC simulations in the mc_muca program is \(\mu= -3.14\). This exercise
aims to confirm (and possibly refine) some of these values.
[ ]:
# 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']
In this notebook you will conduct grand canonical simulations, at given values of chemical potential \(\mu\), volume \(V\) and temperature \(T\), followed by a series of weighted simulations at the same \(V\) and \(T\) aimed at sampling all values of \(N\) with equal probability: a so-called flat-histogram simulation.
Then you will process the flat-histogram results, to calculate the probability \(p_{\mu}(N)\) of particle number \(N\) in the constant-\(\mu VT\) ensemble, for any chosen value of \(\mu\). The aim is to see two peaks in \(p_{\mu}(N)\), one for each phase, and adjust \(\mu\) to give equal peak areas, i.e. \(50\%\) probability of observing each phase, which is the condition for two-phase coexistence.
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
Grand canonical Monte Carlo simulations¶
The program mc_muca carries out a simulation of a single system, with fixed box dimensions, at a specified temperature \(T\), in which the number of particles \(N\) is allowed to vary. Optionally, it reads in data from a file, mc_muca.hdf5, and calculates a set of weights \(\Phi(N)\), which will be useful later. To begin with, however, this file is not supplied, and in these circumstances all the weights are set to the values appropriate for the grand canonical ensemble,
namely \(\Phi(N)=\mu N\), where \(\mu\) is the chemical potential. Default choices for \(T\) and \(\mu\) are in the file mc_muca.f90; \(\mu\) should be reasonably close to the coexistence value at the chosen \(T\).
The first run will start from a configuration file at liquid density supplied in the file DATA/mc_muca.xyz. The file is in extended XYZ format, and looks 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, 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 this file, converted into units appropriate for argon atoms, using weas-widget. As usual, you may wish to adjust weas_shape (defined in Preliminaries above) or use the widget’s ⛶ button (may not work in Colab) to change the viewer size. You can use the controls Atoms/Radius Type/VDW to give a space-filling view. Rotate it around, to see that we are using an elongated simulation box.
[ ]:
this_dir = Path.cwd() # Current working directory
data_dir = this_dir / 'DATA' # Data directory
old_conf = data_dir / 'mc_muca.xyz' # Supplied configuration file
new_conf = copy(old_conf,this_dir) # Copy supplied file into current directory
atoms=read(new_conf)
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
The box dimensions are \(5\times5\times15\). The reason for using such an elongated box will become clear later; note that it is on the margins of acceptability, being only twice as wide as \(r_{\mathrm{cut}}= 2.5\).
Start the run in the following cell. Before running, the file mc_muca.hdf5, if it exists, is removed.
The main program file (in the src subdirectory) is mc_muca.f90, and the particle insertion and deletion is handled by the subroutine n_move in mc_module.f90. While the program is running, feel free to open these files, and look through them to figure out what is going on.
[ ]:
hdf5_file = this_dir / 'mc_muca.hdf5' # HDF5 file
hdf5_file.unlink(missing_ok=True) # Remove it, if it exists
!echo '&nml /' | build/mc_muca
You’ll need to wait until the run has finished before proceeding. Then, read in the simulation results from the HDF5 file.
[ ]:
attr, data = read_hdf5_file(hdf5_file)
The main object of interest is the probability distribution \(p_{\mu}(N)\). An un-normalized version of this, \(h(N)\), is stored as the dataset data['h'], the histogram containing the number of occurrences of each value of \(N\) during the simulation. This is an integer NumPy array of Nmax+1 elements, with indices starting at N=0 and ending at N=Nmax.
Notice that, for this program, we have not output any step-by-step datasets. The other ones are data['wt'], containing statistical weights used and refined by the simulation program, and data['mu'], containing estimates of \(\mu(N)\), for \(N=1,2,\ldots,N_\text{max}\), which we will turn to later on.
First, we print some simulation parameters (attributes).
[ ]:
print(attr['Title'].astype(str))
print('Number of steps',attr['nstep'])
Nmax = attr['Nmax']
V = attr['V']
T = attr['T']
mu = attr['mu']
L = attr['L']
print(f'Box Lx,Ly,Lz = {L[0]:10.4f}{L[1]:10.4f}{L[2]:10.4f}')
print(f'Maximum N Nmax = {Nmax:10d}')
print(f'Volume V = {V:10.4f}')
print(f'Temperature T = {T:10.4f}')
print(f'Chem pot mu = {mu:10.4f}')
Now we normalize the \(h(N)\) histogram, print the average \(N\) for this simulation, and also the range of \(N\) sampled during the run. The expected coexistence values of \(N\), for this volume, are given for comparison. Then the distribution is plotted.
[ ]:
h = data['h']
p = h/h.sum()
N_avg = np.average(np.arange(Nmax+1),weights=p) # Average N in this simulation
N_vap, N_liq = rho_vap*V, rho_liq*V # Literature coexistence values of N for this volume
sampled, = h.nonzero() # Values of N sampled (will be one contiguous block)
lo, hi = sampled[0], sampled[-1]
print(f'Simulation sampled N {lo:3d} - {hi:3d}')
print(f'Simulation average N {N_avg:6.1f}')
print(f'EOS coexistence N {N_vap:6.1f}{ N_liq:6.1f}')
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim=(0,0.06)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\mu}(N)$')
ax.plot(p)
ax.axvline(N_avg)
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
plt.show()
It is most likely that the system will stay in its initial phase, liquid (but it doesn’t matter if a switch to vapour happens). In most cases, this will show a single peak, around the simulation average value of \(N\), which is plotted as a solid vertical line in the above plot. We also plot dotted lines corresponding to the expected values of the coexisting liquid and vapour densities from the literature (as mentioned at the top of the worksheet). We could calculate similar averages, such as the pressure and the potential energy, from our simulation, which is just a grand-canonical ensemble simulation of the liquid phase. However, in this workshop, we want to focus on the histogram of \(N\).
The second run will start from a supplied configuration file DATA/mc_muca_empty.xyz containing no atoms. Start the run in the following cell. It should take less time than the previous one. Note that to ensure that the program carries out a straightforward grand canonical simulation, the file mc_muca.hdf5 containing weights is removed.
[ ]:
old_conf = data_dir / 'mc_muca_empty.xyz' # Supplied, empty, configuration
copy(old_conf,new_conf) # Copy to mc_muca.xyz in current directory
hdf5_file.unlink(missing_ok=True) # Remove HDF5 file, if it exists
!echo '&nml /' | build/mc_muca
Once more, wait until the run has finished before proceeding.
Use the following cells to plot the \(h(N)\) histogram file, once more normalized.
[ ]:
attr, data = read_hdf5_file(hdf5_file)
[ ]:
print(attr['Title'].astype(str))
print('Number of steps',attr['nstep'])
Nmax = attr['Nmax']
V = attr['V']
T = attr['T']
mu = attr['mu']
print(f'Box Lx,Ly,Lz = {L[0]:10.4f}{L[1]:10.4f}{L[2]:10.4f}')
print(f'Maximum N Nmax = {Nmax:10d}')
print(f'Volume V = {V:10.4f}')
print(f'Temperature T = {T:10.4f}')
print(f'Chem pot mu = {mu:10.4f}')
[ ]:
h = data['h']
p = h/h.sum()
N_avg = np.average(np.arange(Nmax+1),weights=p) # Average N in this simulation
N_vap, N_liq = rho_vap*V, rho_liq*V # Literature coexistence values of N for this volume
sampled, = h.nonzero() # Values of N sampled (will be one contiguous block)
lo, hi = sampled[0], sampled[-1]
print(f'Simulation sampled N {lo:3d} - {hi:3d}')
print(f'Simulation average N {N_avg:6.1f}')
print(f'EOS coexistence N {N_vap:6.1f}{ N_liq:6.1f}')
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim = (0,0.06)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\mu}(N)$')
ax.plot(p)
ax.axvline(N_avg)
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
plt.show()
It is most likely that the system will stay in its initial phase, vapour (but it doesn’t matter if a switch to liquid happens).
If these runs were extended, it would become more likely (provided \(\mu\) is very close to the coexistence value) that the system would eventually sample both phases, but to do this means crossing a free energy barrier at intermediate values of \(N\). At best, the sampling in this region will be poor. The simulation needs some help!
Multicanonical Monte Carlo simulations¶
Recall from the lecture that in a multicanonical MC simulation with weights \(\Phi(N)\), the distribution of the number of atoms \(N\) is \begin{equation*} p_{\Phi}(N) \propto Q(N) \, e^{\beta\Phi(N)} \end{equation*} where \(\beta=1/k_{\text{B}}T\), or \begin{equation*} -k_{\text{B}}T\ln p_{\Phi}(N) = F(N) - \Phi(N) + \text{constant} \end{equation*} where \(Q(N)\) is the canonical partition function, and \(F(N)=-k_{\text{B}}T\ln Q(N)\) the Helmholtz free energy. For simplicity, the dependence on \(VT\) is not written explicitly here. If \(\Phi(N)=\mu N\) this is just the grand canonical distribution \(p_{\mu}(N)\); but the aim here is to find weights that will make \(p_{\Phi}(N)\approx\) constant. This will be true if \begin{equation*} \Phi(N) = F(N) . \end{equation*} Our first guess at the multicanonical weights \(\Phi(N)\) will not be correct, but this equation can be used to iteratively improve them. Having measured \(p_{\Phi}(N)\) from a simulation using weights \(\Phi(N)\), we estimate the free energy \begin{equation*} F(N) = -k_{\text{B}}T\ln p_{\Phi}(N) + \Phi(N) , \end{equation*} and then set \(\Phi(N) := F(N)\) for the next run. The hope is that these new \(\Phi(N)\) will generate a flatter \(p_{\Phi}(N)\) in the next run. In practice the program works with differences in free energies, i.e. estimates of the chemical potential \(\mu(N)\equiv F(N)-F(N-1)\), and ratios of probability histograms \(p_{\Phi}(N)/p_{\Phi}(N-1)\).
At the completion of each run, mc_muca outputs the necessary data to the file mc_muca.hdf5. Specifically, it contains the following datasets:
h, the number histogram \(h(N) \propto p_{\Phi}(N)\) which we have seen already;mu, an estimate of \(\mu(N)\) from which we can calculate \(F(N)\) and hence \(\Phi(N)\);wt, the statistical weight of that estimate.
We shall perform a series of runs, each one reading in the final configuration and weights from the previous one. The initial run, run \(0\), is the one just conducted (with no weights). We save the results of this run, and all subsequent ones, in a subdirectory named MUCA, which we refer to as result_dir.
[ ]:
result_dir = this_dir / 'MUCA' # Subdirectory to contain results
result_dir.mkdir(exist_ok=True) # Create subdirectory, unless it exists already
name = '0.hdf5'
copy(hdf5_file,result_dir/name) # Copy HDF5 file to result directory, renaming it
name = '0.xyz'
copy(new_conf,result_dir/name); # Copy configuration to result directory, renaming it
The following cell performs the sequence of runs, \(1, 2, \ldots, 5\), refining the weights after each run, and saving the results in result_dir. This will take a few minutes.
[ ]:
for run in range(1,6):
! echo '&nml /' | build/mc_muca
name = str(run)+'.hdf5'
copy(hdf5_file,result_dir/name) # Copy HDF5 file into result directory, renaming it
name = str(run)+'.xyz'
copy(new_conf,result_dir/name) # Copy configuration file into result directory, renaming it
print('Completed run ',str(run))
print('All done')
If time is pressing, it is OK to interrupt the kernel, and carry out the next step with whatever results have been generated so far. Otherwise, of course, feel free to take a break!
The idea of these runs is to show a steady improvement in the sampling, generating probability histograms that are increasingly flat. Because these runs are not very long, for this system, things may not work out perfectly, but hopefully they will give an idea of what to expect. Let’s see!
The next cell will look for all the HDF5 files in the result directory, and plot the \(h(N)\) histograms, normalized to give \(p_{\Phi}(N)\).
[ ]:
files = sorted(result_dir.glob('*.hdf5'))
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim = (0,0.06)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\Phi}(N)$')
for file in files:
attr, data = read_hdf5_file(file)
h = data['h']
p = h/h.sum()
ax.plot(p,label=file.stem)
N_vap, N_liq = rho_vap*attr['V'], rho_liq*attr['V']
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
ax.legend(loc='upper center',ncol=2)
plt.show()
In case things didn’t go well, or just for comparison, we have provided a similar set of files in the DATA subdirectory. Just replace result_dir by data_dir at the start of the cell above.
You may also like to visualize the final configuration from each of your runs, 1.xyz, 2.xyz etc., stored in result_dir, as done in the next cell. These might be liquid-like, vapour-like, or a mixture of the two. Once more, we have provided some example configurations in data_dir; file 4.xyz in that directory contains a mixed configuration.
[ ]:
file = result_dir / '1.xyz'
atoms = read(file)
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
Exercise: locating coexistence¶
Ideally, after refining the weights enough times to generate an essentially flat sampled distribution, a very long run should be carried out, perhaps \(10\times\) or \(100\times\) longer. You may like to try that after the workshop is finished, but for now the results of such a run have been supplied in the file DATA/long.hdf5. The next cell loads the essential data from that file, which has the same format as the HDF5 files described above. This time, we will be interested in two
datasets: data['h'] and data['mu']. We start by plotting the first of these.
[ ]:
file = data_dir / 'long.hdf5'
attr, data = read_hdf5_file(file)
[ ]:
Nmax = attr['Nmax']
V = attr['V']
T = attr['T']
print(f'Maximum N Nmax = {Nmax:10d}')
print(f'Volume V = {V:10.4f}')
print(f'Temperature T = {T:10.4f}')
[ ]:
h = data['h']
p = h/h.sum()
N_vap, N_liq = rho_vap*V, rho_liq*V # Literature coexistence values of N for this volume
sampled, = h.nonzero() # Values of N sampled (will be one contiguous block)
lo, hi = sampled[0], sampled[-1]
print(f'Simulation sampled N {lo:3d} - {hi:3d}')
print(f'EOS coexistence N {N_vap:6.1f}{ N_liq:6.1f}')
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim = (0.003,0.004)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\Phi}(N)$')
ax.plot(p)
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
plt.show()
The long run generates a reasonably flat histogram; the high-\(N\) end is slightly less well sampled. Compare the vertical scale of this plot with the results of your sequence of multicanonical runs.
The data['mu'] array contains values of \(\mu(N)\) for \(N=1\ldots N_{\text{max}}\). (Note: due to Python/NumPy indexing, \(\mu(1)\) is stored in element 0, \(\mu(2)\) in element 1, etc). This will be used to compute \(F(N)\), where \(F(N)\) is the corresponding Helmholtz free energy. We use the cumulative sum function, based on \(F(N)=F(N-1)+\mu(N)\) for \(N=1\ldots N_{\text{max}}\). A similar formula was used in the Fortran code to compute the
weights; we define \(F(0)=0\) as a reference point.
[ ]:
F = np.empty(Nmax+1,dtype=np.float64)
F[0] = 0.0
F[1:] = np.cumsum(data['mu'])
The following function calculates the unweighted probability distribution \(p_{\mu}(N)\) for any given \(\mu\), using the formula \(p_{\mu}(N) \propto Q(N) \exp(\beta\mu N)\) where \(Q(N)=\exp(-\beta F(N))\). The array F is supplied to this function. In the following cells, mu is the chosen value of \(\mu\) (don’t confuse with the data['mu'] array that we used in the previous cell) and recall that \(\beta=1/k_{\text{B}}T\).
[ ]:
def prob ( T, mu, F ):
"""Probability distribution for number of atoms at given chemical potential.
Arguments
---------
T : float, scalar
temperature of simulation
mu : float, scalar
chosen chemical potential
F : float, NumPy array
estimates of F(N) for N = 0 .. Nmax inclusive
Returns
-------
float, Numpy array
normalized probability distribution for N = 0 .. Nmax inclusive
"""
muN = mu*np.arange(F.size) # mu*N for every N
P = np.exp ( ( -F + muN ) / T ) # Un-normalized probabilities
p = P/np.sum(P) # Normalized probabilities
return p
The next few cells set up an interactive plot for \(p_{\mu}(N)\), with a slider to adjust \(\mu\). You should find that the plot is very sensitive to the value of \(\mu\).
The total probabilities in the first half (vapour) and the second half (liquid) of this distribution are also calculated and displayed. Hopefully, \(p_{\mu}(N)\) will be very small in the middle of the range, so the result should not depend sensitively on this choice of splitting. The aim is to get two peaks with equal probability.
After your first try, you should adjust the values of mu_min and mu_max, bringing them closer together, and re-run the cell, so as to zero in more accurately on the coexistence value of \(\mu\). It should be possible to determine \(\mu\) to 3 or 4 decimal places.
[ ]:
# Set up the figure. Just do this once.
with plt.ioff():
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
[ ]:
# Return to this cell when you want to adjust mu_min and mu_max values to bracket the coexistence value
# Remember to keep mu_min < mu_max (they are both negative)
mu_min=-3.3
mu_max=-3.0
mu_mid=0.5*(mu_max+mu_min) # Ensure that starting value is in range
mu_del=(mu_max-mu_min)/50 # Give a reasonable number of intervals in this range
# Set up the slider widget
slider=FloatSlider(min=mu_min,max=mu_max,step=mu_del,value=mu_mid,
readout_format='.4f',description='mu')
slider.layout.width = '80%'
[ ]:
plt.cla() # Clear the axes
ylim = (0.0,0.08)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\mu}(N)$')
mu = slider.value
p = prob ( T, mu, F )
vap,liq=tuple(np.array_split(p,2)) # fractions of vapour and liquid
ann_vap=ax.annotate('{:8.5f}'.format(vap.sum()),xy=(50,0.06),ha='center')
ann_liq=ax.annotate('{:8.5f}'.format(liq.sum()),xy=(208,0.06),ha='center')
lines=ax.plot(p)
N_vap, N_liq = rho_vap*V, rho_liq*V
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
def update(change):
mu = change.new
p = prob ( T, mu, F )
lines[0].set_ydata(p)
vap,liq=tuple(np.array_split(p,2)) # fractions of vapour and liquid
ann_vap.set_text('{:8.5f}'.format(vap.sum()))
ann_liq.set_text('{:8.5f}'.format(liq.sum()))
fig.canvas.draw()
fig.canvas.flush_events()
slider.observe(update, names='value')
VBox([slider,fig.canvas])
The value of \(\mu\) may be extracted from the slider, for further refinement.
[ ]:
mu=slider.value
print(mu)
Optional exercise: automatic optimization¶
Adjusting \(\mu\) by hand like this is rather fiddly. It is more satisfactory to optimize this parameter automatically, to give equal vapour and liquid fractions within a small numerical tolerance. If time permits, try to do this; otherwise skip to the next section, on Coexistence.
To help, below we define a Python function f(mu) (the difference between vapour and liquid fractions); the numerical task is to solve the equation f(mu)=0 for mu. Here, we suggest that you use the function brentq(f,mu_min,mu_max) which has been imported from the scipy.optimize sub-package. This function just needs f, and a pair of limits mu_min and mu_max which should bracket the root, and it returns the root, i.e. the value of mu satisfying f(mu)=0.
[ ]:
def f(mu):
p = prob ( T, mu, F )
vap,liq=tuple(np.array_split(p,2))
return vap.sum()-liq.sum()
[ ]:
# Insert your code here to determine an accurate value of mu
[ ]:
print(f'Best value of mu = {mu:15.10f}')
Coexistence¶
After you have found your best estimate of the coexistence value of \(\mu\), compare it with the value estimated from the Gibbs simulation programs (covered in the accompanying notebook), and with the fitted EOS estimate mu_coex given in the introduction to this notebook.
To finish this exercise, we shall plot \(p_{\mu}(N)\), and the effective, or Landau, free energy which we shall call \(F_{\mu}(N)=-k_{\text{B}}T \ln p_{\mu}(N)\), at \(\mu=\mu_{\text{coex}}\).
At the exact coexistence point, the plot of \(F_{\mu}(N)\) vs \(N\) will have a horizontal region in the middle, corresponding to a range of densities over which the system consists of two slab-like regions (one liquid, one vapour) in the simulation box. Using an elongated \(5\times 5\times 15\) box helps favour these configurations. It is possible to estimate the surface tension \(\gamma\) of the liquid-vapour interface by comparing the height of this plateau with the free energy minima. The difference is \(\Delta F_{\mu}=2A\gamma\) where \(A=5\times5\) is the cross-sectional area. The factor \(2\) is because there are two interfaces. At this temperature, \(T=0.95\), the value of the surface tension is \(\gamma\approx0.15\) (see S Stephan et al., J Phys Chem C, 122, 24705 (2018)).
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim = (0.0,0.03)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$p_{\mu}(N)$')
p = prob ( T, mu, F )
ax.plot(p)
vap,liq=tuple(np.array_split(p,2))
ax.annotate('{:8.5f}'.format(vap.sum()),xy=(50,0.025),ha='center')
ax.annotate('{:8.5f}'.format(liq.sum()),xy=(208,0.025),ha='center')
N_vap, N_liq = rho_vap*V, rho_liq*V
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
plt.show()
As explained above, the free energy associated with two interfaces is \(\Delta F_{\mu}=2A\gamma\). In the simulated system, \(A=5\times5=25\), and \(\gamma\approx0.15\) at this temperature, \(T=0.95\) (see S Stephan et al., J Phys Chem C, 122, 24705 (2018)), so \(\Delta F_{\mu} \approx 2\times25\times 0.15 \approx 7.5\). Recall that \(F_{\mu}=-k_{\mathrm{B}}T\ln p_{\mu}(N)+C\) and in our units \(k_{\mathrm{B}}=1\). An
equivalent formula is \(F_{\mu}(N)=F(N)-\mu N+C\) (look at the definition of prob(T,mu,F) above). A plot of \(F_{\mu}(N)\) against \(N\) should have two free energy minima, one for each phase, separated by a flat plateau corresponding to the two-phase system containing two interfaces. It should be possible to read off the value of \(\Delta F_{\mu}\). Let’s see!
[ ]:
L = attr['L'] # Simulation box dimensions
A = L[0]*L[1] # Cross-sectional area
gamma = 0.15 # Approximate literature value of surface tension
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ylim=(0.0,8.0)
ax.set_ylim(ylim)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$F_{\mu}(N)$')
p = prob ( T, mu, F )
Free = -T*np.log ( p ) # You may use either this formula ....
#Free = F - mu*np.arange(F.size) # .... or this equivalent one
Free = Free - np.min(Free) # Define zero at Fmin, for plot
ax.plot(Free)
N_vap, N_liq = rho_vap*V, rho_liq*V
ax.vlines([N_vap,N_liq],*ylim,linestyles='dotted',colors='k')
ax.axhline(y=2*A*gamma,linestyle='dotted',color='k')
plt.show()
The result should be reasonably close to the literature value of \(\Delta F_{\mu} \approx 7.5\). Remember, a serious study would use a larger system, and longer runs.
Further Work¶
As a by-product of the flat-histogram approach, we have effectively done \(NVT\) simulations across a wide range of densities, and could, in principle, obtain some useful thermodynamic results. For instance, the pressure \(P(\rho)\) could be calculated using the virial expression, if we had averaged and tabulated this function according to the value of \(N\). We didn’t do this in the simulations. However, we can also calculate \(P(\rho)\) from the datasets at hand. Recall that
\(G=F+PV\) and that \(G=\mu N\) for a single-component system. We read in the data['mu'] array again, which contains \(\mu(N)\) for \(N=1,\ldots, N_\text{max}\). Then we can calculate \(F(N)\) and \(G(N)\), and hence obtain \(P(N)\) for the given \(V\), hence \(P(\rho)\).
The following cells perform this calculation, plot the resulting equation of state \(P(\rho)\), and compare with the fitted EOS described in earlier workshops. Remember, though, that the fitted EOS is only valid in bulk phases. Therefore, in the coexistence region, we remove those values just before plotting.
[ ]:
file = data_dir / 'long.hdf5'
attr, data = read_hdf5_file(file)
[ ]:
Nmax = attr['Nmax']
V = attr['V']
T = attr['T']
# data['mu'] contains values of mu(N) for N=1,2,3,...
# For convenience we formally insert mu[0]=0
# This ensures that the indices of the array mu match the values of N
mu = np.empty(Nmax+1,dtype=np.float64)
mu[0] = 0.0
mu[1:] = data['mu']
# All the arrays are now indexed correctly by N=0,1,2,3,...
F = np.cumsum(mu)
G = mu * np.arange(Nmax+1)
P = (G - F) / V
rho = np.arange(Nmax+1) / V
# Fill Peos array, taking care to set Peos[0]=0 explicitly
Peos = np.empty(Nmax+1,dtype=np.float64)
Peos[0] = 0.0
Peos[1:] = np.array([ eos(density=den,temperature=T)['P'] for den in rho[1:] ], dtype=np.float64)
[ ]:
fig, ax = plt.subplots()
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
ax.set_xlabel(r'$\rho$')
ax.set_ylabel(r'$P(\rho)$')
ax.plot(rho,P,label='simulated')
keep = (rho<rho_vap) | (rho>rho_liq) # Keep only values in single phase regions
Peos=Peos[keep]
rho =rho [keep]
ax.plot(rho,Peos,label='fitted EOS')
ylim=ax.get_ylim()
ax.vlines([rho_vap,rho_liq],*ylim,linestyles='dotted',colors='k')
ax.legend(loc='upper center')
plt.show()
In the coexistence region, the simulation results contain contributions from the interfaces. Perfect agreement with the fitted EOS should not be expected, but hopefully it is close.
This concludes the notebook. If you have completed the Gibbs simulation notebook as well, it is the end of this workshop.
[ ]: