Monte Carlo at Constant Pressure¶
This workshop is intended to demonstrate the operation of a constant-pressure Monte Carlo program.
This exercise builds on some of the ideas already covered in the MC-Lennard-Jones.ipynb notebook.
Introduction¶
In a system with pairwise-additive interactions \(u(r)\) the virial expression for the pressure is given (see lecture notes) in terms of the pair virial function \(r (du/dr)\). This function is calculated and accumulated in the same way as the pair potential. In a constant-pressure simulation the value of \(P\) is input as a parameter, and there is no need to use the expression just mentioned, although it is always good to check that the measured average matches the input value.
In such a simulation, the system size varies: the box volume is changed randomly \(V\rightarrow V+\Delta V\), and the move accepted or rejected using a formula which depends on the specified value of \(P\). Such a volume-changing move is typically attempted once per step, after \(N\) attempted molecule displacements, where \(N\) is the number of molecules. We can obtain the equation of state in the form of the average density \(\rho\) at given \(P\), rather than averaging \(P\) at fixed \(\rho\).
The system being studied here is the same one that was investigated in the MC-Lennard-Jones.ipynb notebook, namely the cut-and-shifted Lennard-Jones potential with \(r_{\text{cut}}=2.5\) (in reduced units). Further details are given in that notebook. Moreover, the state point of interest is the same, \(T=2.0\), \(\rho=0.5\), which lies in the supercritical region of the phase diagram. This state point was discussed in lectures.
Setup (optional)¶
The next cell needs to be run only if you use google colab. Commands may work in other environments too but are not tested.
[ ]:
! pip install h5py data_tutorials
! apt install gfortran libhdf5-dev
now grab the data needed for tutorial
[ ]:
from data_tutorials.data import get_data
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_4/MC_Tutorial/",
filename=["hdf5_module.py","eos_lj.py"],
folder=".",
)
# get the files for MC programme
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_4/MC_Tutorial/",
filename=["Makefile", "config_old.bak","config_old.dat","config_io_module.f90","hdf5_module.f90","maths_module.f90","mc_module.f90","mc_npt.f90","mc_npt_module.f90","potential_module.f90"],
folder=".",
)
Preliminaries¶
Start by importing some useful Python modules and functions.
[ ]:
import numpy as np
import matplotlib.pyplot as plt
import shutil
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
The following cell should build the program mc_npt.
[ ]:
!make mc_npt
The main program file is mc_npt.f90. Open this file and take a look at it. The program carries out a Monte Carlo simulation for a system of \(N\) atoms at specified temperature \(T\) and pressure \(P\).
The basic program structure in the file mc_npt.f90 is very similar to mc_nvt.f90. Simulation parameters are read from standard input using a namelist. A run consists of nstep steps. In every step, \(N\) attempts are made to move individual atoms, selected randomly; in addition, here, there is a call to the volume-move routine. There are a few more variables and parameters associated with the volume moves. The volume \(V\) is added to the datasets of variables saved
step-by-step to the HDF5 file, which is here called mc_npt.hdf5. In this program, for simplicity, we do not save configurations at intervals, and will not be calculating \(g(r)\).
Start by doing a short run. Just in case the originally supplied configuration file config_old.dat has been overwritten, it is copied from the backup.
[ ]:
shutil.copy('config_old.bak','config_old.dat')
!echo '&nml nstep=1000 /' | ./mc_npt
The program does not appear to be working properly: the volume move acceptance ratio is zero. The first exercise will be to fix the program so that it samples the constant-pressure ensemble correctly.
Fixing the program¶
Examine the v_move subroutine in the file mc_npt_module.f90. The trial move is defined by choosing \(\Delta V\) uniformly within a specified range \(\pm\Delta V_\text{max}\). The volume change is stored in the variable dv. (There are various different, non-equivalent, ways of choosing \(\Delta V\), which may need a different acceptance formula from the one given below). From the new volume, we calculate the volume scaling factor v_scale, and we deduce the scaling
factor r_scale which applies to the box lengths, and to the atomic coordinates. Hopefully, that all looks correct in the program.
- It is then necessary to calculate the total potential energy of the scaled system. The potential energy change \(\Delta U\) associated with the trial move is a key element of the acceptance probability. Referring to the lecture notes, this may be written \(\min [1,e^{-\delta}]\) where :nbsphinx-math:`begin{equation*}
e^{-delta} equiv left(frac{V_text{new}}{V_text{old}}right)^N , e^{-beta P(V_text{new}-V_text{old})}, e^{-beta(U_text{new}-U_text{old})} = left(frac{V+Delta V}{V}right)^N , e^{-beta PDelta V} , e^{-beta Delta U} .
end{equation*}` Here we have written \(V_\text{old}=V\), \(V_\text{new}=V+\Delta V\), with similar expressions for \(U\), and \(\beta=1/k_{\text{B}}T\) as usual. The program defines a variable delta, to represent the value \(\delta\) here. This is passed to the metropolis function, which handles the \(\min [1,e^{-\delta}]\) part of the calculation and returns the accept/reject decision. Rearrange the above equation to give an explicit expression for
\(\delta\); you will need to take logs of both sides. Compare your expression for \(\delta\), i.e. delta, with the program. You should see that the potential energy part is handled correctly in the statement
delta = beta * ( u_new - u )
and for constant-volume MC, this is all we would need. However, for constant pressure, two extra terms need to be added to delta, as the above expression indicates. One or both of these two extra terms may be wrong in the supplied code. They may be easily expressed in terms of the variables v_scale and dv, as well as the number of atoms N and variables beta and pressure. Make the necessary correction(s).
Once you’ve finished editing the mc_npt_module.f90 file, and have saved it, continue by making the program again, and re-running in the following cells.
[ ]:
!make mc_npt
[ ]:
!echo '&nml nstep=1000 /' | ./mc_npt
Now the volume move acceptance rate should be around 60-70%. If not, check your formulae, and the v_move routine, again.
Even for such a short run, it will be worth taking a look at the output. Just as for the \(NVT\) program, the next cell reads in the simulation parameters and datasets from the HDF5 output file.
[ ]:
params, data = read_file('mc_npt.hdf5')
We give names to some of the important parameters, and print them out.
[ ]:
print(params['Title'].astype(str))
print('Run steps',params['nstep'])
N = params['N']
P = params['P']
T = params['T']
print(f'Number of atoms N = {N:10d}')
print(f'Pressure P = {P:10.4f}')
print(f'Temperature T = {T:10.4f}')
First, look at the evolution of the density. We expect to see some early variation of this quantity, since the initial configuration (with density \(\rho=0.5\)) was not equilibrated at \(P=1.0\). A (very) rough estimate of the average density at this pressure may be obtained by discarding the early part of the simulation. This is done in the following cell. As a guide, the value of rhoavg should be around 0.43. If it is very different, once more you may need to check your formulae
and the v_move routine again.
[ ]:
rho = N/data['V']
rhoavg = rho[200:].mean()
print(f'Simulation average rho = {rhoavg:10.4f}')
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel('step')
ax.set_ylabel(r'$\rho$')
ax.plot(rho);
ax.axhline(rhoavg,ls='dashed',color='C1')
plt.tight_layout()
We can double check by calculating the value of \(P\) that would be expected at this average density according to the fitted EOS. This should be reasonably close to the simulation input value \(P=1.0\).
[ ]:
eos_fit = eos(temperature=T,density=rhoavg)
Peos = eos_fit['P']
print(f'Fitted EOS Pressure P = {Peos:10.4f}')
print(f'Specified value of P = {P:10.4f}')
As a last check, compare the average configurational temperature for this very short run with the simulation input value.
[ ]:
Tavg=data['T'][200:].mean()
print(f'Simulation average T = {Tavg:10.4f}')
print(f'Specified value of T = {T:10.4f}')
If you wish, feel free to experiment with the simulation parameters. You should be able to adjust dv_max to increase or decrease the volume move acceptance ratio, and you can choose different pressures and temperatures, comparing with the expected equation of state. Remember that the runs we have done so far are very short indeed; for more accuracy, you could use the default number of steps, but this will take several minutes.
The state point of most interest to us is the one discussed in the Monte Carlo lecture: \(\rho=0.5\) and \(T=2.0\). In the lecture, you were given the corresponding value of \(P\). Also, you measured it in the last workshop, by running mc_nvt. If neither of those values is to hand, the fitted equation of state gives an estimate.
[ ]:
eos_fit=eos(temperature=2.0,density=0.5)
Peos = eos_fit['P']
print(f'Fitted EOS Pressure P = {Peos:10.4f}')
Conduct a run with \(P\) estimated from one of these sources (you may wish to change the numerical value in the next cell). The run should take a few minutes.
[ ]:
!echo '&nml pressure=1.3384/' | ./mc_npt
You should wait until the program finishes before proceeding.
Once more, we read in the simulation results and calculate some averages, to confirm that all is well. There should be no need to discard the early part of the run, since the initial configuration was already equilibrated at the desired state point. The first check is that the average density is close to the desired value \(\rho=0.5\). Then we check the average (configurational) temperature and pressure against the set values, and the activity \(z\), estimated by Widom test particle insertion (which gives \(1/z\)), compared with the fitted EOS value.
[ ]:
params, data = read_file('mc_npt.hdf5')
[ ]:
print(params['Title'].astype(str))
print('Run steps',params['nstep'])
N = params['N']
P = params['P']
T = params['T']
Vavg = data['V'].mean()
print(f'Number of atoms N = {N:10d}')
print(f'Simulation average rho = {N/Vavg:10.4f}')
Tavg = data['T'].mean()
print(f'Specified temperature T = {T:10.4f}')
print(f'Simulation average T = {Tavg:10.4f}')
Wavg = data['W'].mean()
Pavg = N*T/Vavg + Wavg/(3*Vavg)
Peos = eos_fit['P']
print(f'Specified pressure P = {P:10.4f}')
print(f'Simulation average P = {Pavg:10.4f}')
print(f'Fitted EOS pressure P = {Peos:10.4f}')
zavg = 1/data['Z'].mean()
zeos = eos_fit['z']
print(f'Simulation estimate z = {zavg:10.4f}')
print(f'Fitted EOS value z = {zeos:10.4f}')
Hopefully the simulation averages agree reasonably well with the predictions of the fitted EOS and the input parameters of the simulation. As always, if we wanted to properly judge whether our results agree with expectations, we would need to estimate statistical errors on the simulation averages, which is not the topic for today.
Isothermal compressibility¶
- The volume fluctuations may be used to calculate the isothermal compressibility, which was defined in lectures :nbsphinx-math:`begin{equation*}
kappa_T = -frac{1}{V}left(frac{partial V}{partial P}right)_T = frac{1}{rho}left(frac{partialrho}{partial P}right)_T = frac{langle V^2rangle - langle V rangle^2}{k_{text{B}}T langle Vrangle} .
end{equation*}` (The symbols \(\beta_T\) and \(\chi_T\) are also sometimes used for this same quantity, but \(\kappa_T\) will be used here).
Referring to the last formula above, which was also given in lectures, calculate, in the following cell, the isothermal compressibility from the data['V'] array, and other quantities output by the program.
Compare the isothermal compressibility, calculated by this route, with the value given in lectures, and the value expected from the fitted EOS, which is already calculated in the cell.
[ ]:
# Insert your code here
kavg = 0.0
keos = eos_fit['kappa_T']
print(f'Simulation average kappa_T = {kavg:10.4f}')
print(f'Fitted EOS value kappa_T = {keos:10.4f}')
It is also possible to calculate this thermodynamic derivative numerically, by measuring \(\langle \rho\rangle\) at two different pressures (one slightly higher, one slightly lower, than the desired pressure) in two separate runs of mc_npt. If you decide to try this, don’t forget that an equilibration period will be needed in both cases. It is an interesting question to ask: Which method is likely more efficient, in terms of CPU time required to achieve the same statistical error: volume
fluctuations or numerical differentiation?
Volume and energy distributions¶
Histograms of the potential energy \(\mathcal{P}(U)\) and volume \(\mathcal{P}(V)\) probability distributions are plotted next. We could use \(\mathcal{P}(U)\) in a reweighting scheme to estimate results at a different \(T\), but the same \(P\), as mentioned in the earlier workshop (and discussed in the 3rd MC lecture). Similarly, we could use \(\mathcal{P}(V)\) to estimate results at a different \(P\), but the same \(T\). If we wished to reweight our simulation results to obtain estimates at a different \(T\) and \(P\), we would need to use the joint distribution \(\mathcal{P}(U,V)\). We will not pursue this, but for interest’s sake, we plot a 2D histogram of \(\mathcal{P}(U,V)\) below. This shows that \(U\) and \(V\) are closely correlated with each other, which is not really surprising.
[ ]:
Veos = N/eos_fit['density'] # desired average density = 0.5
Vavg = data['V'].mean()
Vstd = data['V'].std()
print(f'Fitted EOS value V = {Veos:10.4f}')
print(f'Simulation average V = {Vavg:10.4f}')
print(f'Standard deviation V = {Vstd:10.4f}')
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$V$')
ax.set_ylabel(r'$\mathcal{P}(V)$')
ax.hist(data['V'],density=True,label='Simulation')
ax.axvline(Veos,c='C1',label='Expected average')
ax.legend()
plt.tight_layout()
[ ]:
Ueos = N*eos_fit['u']
Uavg = data['U'].mean()
Ustd = data['U'].std()
print(f'Fitted EOS value U = {Ueos:10.4f}')
print(f'Simulation average U = {Uavg:10.4f}')
print(f'Standard deviation U = {Ustd:10.4f}')
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel(r'$U$')
ax.set_ylabel(r'$\mathcal{P}(U)$')
ax.hist(data['U'],density=True,label='Simulation')
ax.axvline(Ueos,c='C1',label='Expected average')
ax.legend()
plt.tight_layout()
[ ]:
fig, ax = plt.subplots(figsize=(8,5))
ax.set_xlabel('$U$')
ax.set_ylabel('$V$')
ax.hist2d(data['U'],data['V'],bins=50,density=True)
ax.axhline(Veos,c='gray')
ax.axvline(Ueos,c='gray')
plt.tight_layout()
Constant-pressure heat capacity¶
It should be possible to estimate the constant-pressure heat capacity per atom \(c_P = C_P/N\). The value at this state point was given in the lecture, and the estimate from the fitted equation of state is included in the results from the eos function, see below.
To calculate this from a constant-\(NPT\) simulation, a formula was given in the lecture notes. It involves “enthalpy” fluctuations. \begin{equation*} c_P/k_{\text{B}} = C_P/N k_{\text{B}} = \frac{3}{2} + \frac{\langle H^2\rangle - \langle H\rangle^2}{N(k_{\text{B}} T)^2} \end{equation*} Here the “enthalpy” is defined as \(H=U+PV\) where \(U\) is the total potential energy.
Use this formula to calculate \(c_P\) Remember that in our reduced units, \(k_{\text{B}}=1\). The numerical answer should be similar to the one given in the table in your lecture notes, and also similar to the value produced by eos_lj.py, which already appears in the cell below.
[ ]:
# Insert your code here
cavg = 0.0
ceos = eos_fit['c_P']
print(f'Simulation average c_P = {cavg:10.4f}')
print(f'Fitted EOS value c_P = {ceos:10.4f}')
There is a bit of subtlety in the above formula. \(H\), as defined just now, is not really the enthalpy. It is missing the kinetic energy contribution. It is not the non-ideal part of the enthalpy either, because the \(PV\) term is not zero for an ideal gas. In the expression for \(c_P/k_B\) given in lectures, and above, a contribution \(\frac{3}{2}\) appears explicitly. The ideal gas value, of course, is \(c_P/k_B=\frac{5}{2}\), so the value \(\frac{3}{2}\) only represents part of the ideal gas contribution, due to the kinetic energy fluctuations. The remaining ideal gas contribution is not missing: it is included in the \(PV\) term.
This concludes the notebook.
[ ]: