Open In Colab

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.

Setup (optional)

This step is needed only if you run on google colab, be aware the kernel will restart, this is normal. Instructions may work on other systems but are not tested.

[ ]:

! pip install h5py data_tutorials ipympl weas_widget ase
! apt install gfortran libhdf5-dev
get_ipython().kernel.do_shutdown(restart=True)

Allow widgets to be shown in google colab

[ ]:

from google.colab import output
output.enable_custom_widget_manager()

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/DATA/",
    filename=["mc_muca_0.hdf5","mc_muca_1.hdf5","mc_muca_2.hdf5","mc_muca_3.hdf5","mc_muca_4.hdf5","mc_muca_5.hdf5","mc_muca_long.hdf5"],
    folder="DATA/",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/DATA/",
    filename=["config_0.dat","config_1.dat","config_2.dat","config_3.dat","config_4.dat","config_5.dat"],
    folder="DATA/",
)
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_module.f90","mc_muca.f90","mc_muca_module.f90","potential_module.f90"],
    folder=".",
)
get_data(
    url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_5/Phase_Equilibria/",
    filename=["config_liq.dat", "config_vap.dat"],
    folder=".",
)

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. In the following, reduced units are employed, so \(\varepsilon=1\) and \(\sigma=1\), and in addition Boltzmann’s constant is taken to be unity \(k_{\text{B}}=1\). The simulation programs use the cut-and-shifted Lennard-Jones potential 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\).

Preliminaries

Start by importing all the required Python modules and functions.

[ ]:

%matplotlib ipympl
import matplotlib.pyplot as plt
import numpy as np
import os
import shutil
from scipy.optimize import brentq
from ipywidgets import FloatSlider, VBox
from glob import glob
from dat_to_ase import argon
from weas_widget import WeasWidget
from hdf5_module import read_file
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

For the chosen potential, in reduced units, the critical point is at \(T_{\text{c}}=1.0779\), \(\rho_{\text{c}}=0.3190\). The state point of interest here, however, is at a somewhat lower temperature, \(T=0.95\). At this temperature, the liquid and vapour phases may coexist, with densities \(\rho_{\text{liq}}\approx 0.622\) and \(\rho_{\text{vap}}\approx 0.0665\) (see J Vrabec et al, Molec Phys, 104, 1509 (2006)). The pressure of these two phases at coexistence is \(P_{\text{coex}}\approx 0.045\), and the chemical potential \(\mu_{\text{coex}}\approx -3.14\).

Locating these coexistence densities by simulation is part of the aim of this exercise.

[ ]:

# 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

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 \(\mathcal{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 \(\mathcal{P}_{\mu}(N)\), one for each phase, and adjust \(\mu\) to give equal peak areas, which is the condition for two-phase coexistence.

The following bash instruction should build the program used in this notebook: mc_muca.

[ ]:

!make mc_muca

The two most important program files are mc_muca.f90 and mc_muca_module.f90. Open these, and look through them to figure out what is going on.

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, 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 config_liq.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 and the subsequent lines give the coordinates of each atom. The following cell should allow you to visualize this configuration.

[ ]:

atoms=argon('config_liq.dat')
v=WeasWidget()
v.from_ase(atoms)
v

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 above configuration file is copied to the input file config_old.dat. Also, the file mc_muca_old.hdf5, if it exists, is removed. While the program is running, feel free to take a look at the file mc_muca_module.f90 which contains the appropriate move routines.

[ ]:

shutil.copy('config_liq.dat','config_old.dat')
try:
    os.remove('mc_muca_old.hdf5')
except FileNotFoundError:
    pass
!echo '&nml /' | ./mc_muca
shutil.copy('mc_muca.hdf5','mc_muca_liq.hdf5')

You’ll need to wait until the run has finished before proceeding. The main object of interest is the probability distribution \(\mathcal{P}_{\mu}(N)\). An un-normalized version of this, \(h(N)\), as well as other quantities accumulated as functions of \(N\), are stored as datasets in mc_muca.hdf5, which we have copied to mc_muca_liq.hdf5 for safe keeping. Notice that, for this program, we have not output any step-by-step datasets.

[ ]:

params, data = read_file('mc_muca_liq.hdf5')

The following cells read in run parameters, and the \(h(N)\) histogram. The range of \(N\) sampled during the run is printed, then \(h(N)\) is normalized and plotted.

[ ]:

print(params['Title'].astype(str))
print('Number of steps',params['nstep'])
Nmax  = params['Nmax']
V     = params['V']
T     = params['T']
mu    = params['mu']
L     = params['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}')
Nsampled,=np.nonzero(data['h']>0.5)
print(f'N range sampled = {Nsampled[0]:10d}{Nsampled[-1]:10d}')

[ ]:

p = data['h']/data['h'].sum()
N = np.arange(Nmax+1)
fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_ylim([0, 0.06])
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\mu}(N)$')
ax.plot(N,p)
ax.axvline(x=np.sum(p*N),c='C1') # Average N in this simulation
ax.axvline(x=rho_vap*V,ls='dotted',c='k')
ax.axvline(x=rho_liq*V,ls='dotted',c='k');

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 a value of \(N\) that matches the average density multiplied by the volume of the box. This is plotted as a vertical line in the above plot, as well as 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 configuration file config_vap.dat containing no atoms. Start the run in the following cell. It should take less time than the previous one. Note that, once again, just to be sure, any file mc_muca_old.hdf5 containing weights is removed, so the run should be a straightforward grand canonical simulation.

[ ]:

shutil.copy('config_vap.dat','config_old.dat')
try:
    os.remove('mc_muca_old.hdf5')
except FileNotFoundError:
    pass
!echo '&nml /' | ./mc_muca
shutil.copy('mc_muca.hdf5','mc_muca_vap.hdf5')

Once more, wait until the run has finished before proceeding. Use the following cells to plot the \(h(N)\) histogram file, once more normalized.

[ ]:

params, data = read_file('mc_muca_vap.hdf5')

[ ]:

print(params['Title'].astype(str))
print('Number of steps',params['nstep'])
Nmax = params['Nmax']
V    = params['V']
T    = params['T']
mu   = params['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}')
Nsampled,=np.nonzero(data['h']>0.5)
print(f'N range sampled = {Nsampled[0]:10d}{Nsampled[-1]:10d}')

[ ]:

p = data['h']/data['h'].sum()
N = np.arange(Nmax+1)
fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_ylim([0, 0.06])
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\mu}(N)$')
ax.plot(N,p)
ax.axvline(x=np.sum(p*N),c='C1') # Average N in this simulation
ax.axvline(x=rho_vap*V,ls='dotted',color='k')
ax.axvline(x=rho_liq*V,ls='dotted',color='k');

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*} \mathcal{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\mathcal{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 \(\mathcal{P}_{\mu}(N)\); but the aim here is to find weights that will make \(\mathcal{P}_{\Phi}(N)\approx\) constant. This will be true if \begin{equation*} \Phi(N) = F(N) . \end{equation*} This can then be used to construct an (improved) estimate for the multicanonical weights. Roughly, this is based on the formula above, i.e. having measured \(\mathcal{P}_{\Phi}(N)\) from a simulation using weights \(\Phi(N)\), we estimate the free energy \begin{equation*} F(N) = -k_{\text{B}}T\ln\mathcal{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 \(\mathcal{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 \(\mathcal{P}_{\Phi}(N)/\mathcal{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:

  1. h, the number histogram \(h(N) \propto \mathcal{P}_{\Phi}(N)\) which we have seen already;

  2. mu, an estimate of \(\mu(N)\) from which we can calculate \(F(N)\) and hence \(\Phi(N)\);

  3. wt, the statistical weight of that estimate.

The following cell performs such a sequence of runs, \(1, 2, \ldots, 5\), refining the weights after each run. Run number \(0\) is the one just conducted (with no weights). The new configuration, and new weights, produced by that run are copied into the starting files for run \(1\). Copies of the relevant files are also stored for later analysis, each labelled by the run number. This will take a few minutes, and it is possible that you will not see any output until all the runs are finished.

[ ]:

shutil.copy('mc_muca.hdf5','mc_muca_0.hdf5')
shutil.copy('config.dat','config_0.dat')
shutil.move('mc_muca.hdf5','mc_muca_old.hdf5')
shutil.move('config.dat','config_old.dat')
for run in range(1,6):
    ! echo '&nml /' | ./mc_muca
    shutil.copy('mc_muca.hdf5','mc_muca_'+str(run)+'.hdf5')
    shutil.copy('config.dat','config_'+str(run)+'.dat')
    shutil.move('mc_muca.hdf5','mc_muca_old.hdf5')
    shutil.move('config.dat','config_old.dat')
    print('Completed run ',str(run))

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 whose names match the given pattern, and plot the \(h(N)\) histograms, normalized to give \(\mathcal{P}_{\Phi}(N)\).

[ ]:

files=glob('mc_muca_?.hdf5')
files.sort()
fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\Phi}(N)$')
for file in files:
    params, data = read_file(file)
    p = data['h']/data['h'].sum()
    N = np.arange(Nmax+1)
    ax.plot(N,p,label=file[-6:-5])
ax.axvline(x=rho_vap*V,ls='dotted',color='k')
ax.axvline(x=rho_liq*V,ls='dotted',color='k')
ax.legend(loc='upper center',ncol=2);

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 the first statement in the cell above with files=glob('DATA/mc_muca_?.hdf5').

You may also like to visualize the final configuration from each of your runs, config_1.dat, config_2.dat etc., 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 the DATA subdirectory. The supplied file DATA/config_4.dat contains a mixed configuration.

[ ]:

atoms=argon('config_1.dat')
v=WeasWidget()
v.from_ase(atoms)
v

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/mc_muca_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.

[ ]:

params, data = read_file('DATA/mc_muca_long.hdf5')

[ ]:

Nmax = params['Nmax']
V    = params['V']
T    = params['T']
print(f'Maximum N  Nmax = {Nmax:10d}')
print(f'Volume        V = {V:10.4f}')
print(f'Temperature   T = {T:10.4f}')
Nsampled,=np.nonzero(data['h']>0.5)
print(f'N range sampled = {Nsampled[0]:10d}{Nsampled[-1]:10d}')

[ ]:

N = np.arange(Nmax+1)
fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\Phi}(N)$')
p = data['h']/data['h'].sum()
ax.plot(N,p)
ax.axvline(x=rho_vap*V,ls='dotted',c='k')
ax.axvline(x=rho_liq*V,ls='dotted',c='k');

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, especially the later ones.

The data['mu'] array contains values of \(\mu(N)\) for \(N=1\ldots N_{\text{max}}\). 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)\) (a similar formula was used in the Fortran code to compute the weights), for this same range of \(N\). We take \(F(0)=0\) as a reference point, although we don’t store this value. Note that the array elements F[0], F[1] etc. store the values of \(F\) for \(N=1,2,\ldots\), so we redefine the N array accordingly, to cover the range \(1\ldots N_{\text{max}}\).

[ ]:

F = np.cumsum(data['mu'])
N = np.arange(1,Nmax+1)

The following function calculates the unweighted probability distribution \(\mathcal{P}_{\mu}(N)\) for any given \(\mu\), using the formula \(\mathcal{P}_{\mu}(N) \propto Q(N) \exp(\beta\mu N)\) where \(Q(N)=\exp(-\beta F(N))\). The arrays F and N are 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\). The function also returns the total probabilities in the first half and the second half of this distribution, assuming that these will correspond to vapour and liquid, respectively. Hopefully, \(\mathcal{P}_{\mu}(N)\) will be very small in the middle of the range, so the result should not depend sensitively on this choice.

[ ]:

def prob ( T, mu, F, N ):
    """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 = 1 .. Nmax inclusive
    N : float, NumPy array
        sequence of N = 1 .. Nmax inclusive
    Returns
    -------
    float, Numpy array
        normalized probability distribution for N = 1 .. Nmax inclusive
    float, scalar
        total probability weight in first half of array
    float, scalar
        total probability weight in second half of array
    """

    P = np.exp ( ( -F + mu*N ) / T ) # Un-normalized probabilities
    p = P/np.sum(P)                  # Normalized probabilities
    m = p.size//2                    # Mid-point
    p1 = np.sum(p[:m])               # Sum of first half
    p2 = np.sum(p[m:])               # Sum of second half
    return p, p1, p2

The next few cells set up an interactive plot for \(\mathcal{P}_{\mu}(N)\), with a slider to adjust \(\mu\). You should find that the plot is very sensitive to the value of \(\mu\). The values of the areas under the two peaks are displayed along with the plot. The aim is to get two peaks with equal areas.

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(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)

[ ]:

# 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
ax.set_ylim([0, 0.08])
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\mu}(N)$')
mu = slider.value
p, p_vap, p_liq = prob ( T, mu, F, N )
ann_vap=ax.annotate('{:8.5f}'.format(p_vap),xy=(50,0.06),ha='center')
ann_liq=ax.annotate('{:8.5f}'.format(p_liq),xy=(208,0.06),ha='center')
lines=ax.plot(N, p)
ax.axvline(x=rho_vap*V,ls='dotted',c='k')
ax.axvline(x=rho_liq*V,ls='dotted',c='k')
def update(change):
    mu = change.new
    p, p_vap, p_liq = prob ( T, mu, F, N )
    lines[0].set_data(N,p)
    ann_vap.set_text('{:8.5f}'.format(p_vap))
    ann_liq.set_text('{:8.5f}'.format(p_liq))
    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 left and right areas 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 left and right areas); the numerical task is to solve the equation f(mu)=0 for mu. This can be done by a method such as Newton-Raphson, if the derivative of f(mu) with respect to mu is also computed. However, 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 just needs f, and a pair of limits mu_min and mu_max which should bracket the root. brentq returns the root, i.e. the value of mu satisfying f(mu)=0.

[ ]:

def f(mu):
    p, p_vap, p_liq = prob ( T, mu, F, N )
    return p_vap-p_liq

[ ]:

# Insert your code here

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 rough estimate given in the introduction to this notebook.

To finish this exercise, we shall plot \(\mathcal{P}_{\mu}(N)\), and the effective, or Landau, free energy which we shall call \(\mathcal{F}_{\mu}(N)=-k_{\text{B}}T \ln \mathcal{P}_{\mu}(N)\), at coexistence.

At the exact coexistence point, the plot of \(\mathcal{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\mathcal{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(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_ylim([0, 0.03])
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{P}_{\mu}(N)$')
p, p_vap, p_liq = prob ( T, mu, F, N )
ax.plot(N,p)
ax.annotate('{:8.5f}'.format(p_vap),xy=(50,0.025),ha='center')
ax.annotate('{:8.5f}'.format(p_liq),xy=(208,0.025),ha='center')
ax.axvline(x=rho_vap*V,ls='dotted',c='k')
ax.axvline(x=rho_liq*V,ls='dotted',c='k');

As explained above, the free energy associated with two interfaces is \(\Delta\mathcal{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\mathcal{F}_{\mu} \approx 2\times25\times 0.15 \approx 7.5\). Recall that \(\Delta\mathcal{F}_{\mu}=-k_{\mathrm{B}}T\ln\mathcal{P}_{\mu}(N)+C\) and in our units \(k_{\mathrm{B}}=1\). An equivalent formula is \(\Delta\mathcal{F}_{\mu}(N)=F(N)-\mu N\) (look at the definition of prob(T,mu,F,N) above). A plot of \(\Delta\mathcal{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\mathcal{F}_{\mu}\). Let’s see!

[ ]:

A     = params['L'][0]*params['L'][1] # Cross-sectional area of simulation box
gamma = 0.15                          # Approximate literature value of surface tension

[ ]:

fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_ylim([0,8])
ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\mathcal{F}_{\mu}(N)$')
p, p_vap, p_liq = prob ( T, mu, F, N )
Free = -T*np.log ( p ) # This formula is equivalent to the following one
# Free = F - mu*N      # This formula is equivalent to the preceding one
Free = Free - np.min(Free)
ax.plot(N,Free)
ax.axvline(x=rho_vap*V,ls='dotted',c='k')
ax.axvline(x=rho_liq*V,ls='dotted',c='k')
ax.axhline(y=2*A*gamma,ls='dotted',c='k');

The result should be reasonably close to the literature value of \(\Delta\mathcal{F}_{\mu} \approx 7.5\). Remember, a serious study would use a larger system.

Further Work

As a by-product of the flat-histogram approach, we have effectively done \(NVT\) simulations across a wide range of densities, and can therefore obtain some useful thermodynamic results. For instance, the pressure \(P(\rho)\) may be calculated using the virial expression, if we collate results according to the value of \(N\). However, we can also calculate this 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, along with the essential parameters, and recalculate F and N. The following cells combine these together, and plot the resulting equation of state \(P(\rho)\) against 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 replace the predictions of that equation by the known coexistence pressure.

[ ]:

params, data = read_file('DATA/mc_muca_long.hdf5')

[ ]:

Nmax = params['Nmax']
V    = params['V']
T    = params['T']
F   = np.cumsum(data['mu'])
N   = np.arange(1,Nmax+1)
G   = data['mu'] * N
P   = (G - F) / V
rho = N / V

# Fill Peos array
Peos = [ eos(density=den,temperature=T)['P'] for den in rho ]
Peos = np.array(Peos)

# Replace values in coexistence region by literature value
Peos[(rho>rho_vap) & (rho<rho_liq)] = P_coex

[ ]:

fig, ax = plt.subplots(figsize=(8,4))
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.canvas.footer_visible = False
fig.set_tight_layout(True)
ax.set_xlabel(r'$\rho$')
ax.set_ylabel(r'$P(\rho)$')
ax.plot(rho,P,label='simulated')
ax.plot(rho,Peos,label='fitted')
ax.axvline(x=rho_liq,ls='dotted',color='k')
ax.axvline(x=rho_vap,ls='dotted',color='k')
ax.legend();

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.

[ ]: