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.
Colab Setup (optional)¶
The next cells should only be needed if you use Google Colab, so we test the environment for this. Commands may work in other environments too but are not tested. If you need to try, just replace the following by SETUP = True.
[ ]:
import sys
SETUP = "google.colab" in sys.modules
if SETUP:
print('Fetching files for Google Colab')
[ ]:
if SETUP:
# Install additional packages
! pip install h5py data_tutorials weas_widget ase meson ninja
! apt install gfortran libhdf5-dev
[ ]:
if SETUP:
from data_tutorials.data import get_data
# Get the python modules and top-level meson.build needed for this tutorial
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_4/MC_Tutorial/",
filename=["hdf5_module.py","eos_ljcs.py","meson.build"],
folder=".",
)
# Get the starting configuration, in the DATA subdirectory
get_data(
url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/Day_4/MC_Tutorial/DATA/",
filename=["mc_npt.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_4/MC_Tutorial/src",
filename=["meson.build","config_io_module.f90","hdf5_module.f90","maths_module.f90","mc_module.f90","mc_nvt.f90","mc_npt.f90","potential_module.f90"],
folder="src",
)
Preliminaries¶
Start by importing some useful Python modules and functions. Here are links to NumPy, Matplotlib, pathlib, and shutil, in case you need any details. We shall use ase, and the weas-widget tool to visualize atomic configurations. The remaining modules are supplied as part of this workshop. We also apply some Matplotlib settings.
[ ]:
# Standard and external modules
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from shutil import copy
from ase.atoms import Atoms
from ase.io import read
from ase.data import vdw_radii, atomic_numbers
from weas_widget import WeasWidget
# Modules supplied in this workshop
from hdf5_module import read_hdf5_file
from eos_ljcs import eos
# Settings for Matplotlib plots
plt.style.use(['seaborn-v0_8-talk','seaborn-v0_8-darkgrid','seaborn-v0_8-colorblind'])
plt.rc('image',cmap='viridis') # Default colour map
plt.rc('figure',figsize=(8,5)) # Default figure dimensions in inches
plt.rc('figure.constrained_layout',use=True) # Default figure layout option
plt.rc('legend',frameon=True,framealpha=1.0) # Default appearance of figure legend
plt.rc('hist',bins=100) # Default number of bins to use in histograms
# Setting for weas-widget viewer
weas_shape = {"width":"800px","height":"600px"} # Default dimensions of weas widget
The following cells should build, in the build subdirectory, both the program mc_npt needed here, and the program mc_nvt used in the earlier MC-Lennard-Jones.ipynb notebook. There may be no need to do this again, but it will do no harm if you run the cells anyway.
[ ]:
!meson setup build
[ ]:
!meson compile -C build
The main program file is mc_npt.f90, in the src subdirectory. 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. The default temperature is \(T=2\), and pressure is \(P=0.5\), but we will be changing this later. The supplied configuration file DATA/mc_npt.xyz was equilibrated using mc_nvt at \(T=2\) and \(\rho=0.5\), and this is copied into the working directory before the run.
[ ]:
this_dir = Path.cwd() # Current working directory
data_dir = this_dir / 'DATA' # Data directory
old_conf = data_dir / 'mc_npt.xyz' # Supplied configuration file
new_conf = copy(old_conf,this_dir) # Copy supplied file into the working directory
!echo '&nml nstep=1000, pressure=0.5 /' | build/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_module.f90 (in the src subdirectory). 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_module.f90 file, and have saved it, continue by building the program again.
[ ]:
!meson compile -C build
Now re-run the program.
[ ]:
new_conf = copy(old_conf,this_dir) # Copy supplied file into the working directory
!echo '&nml nstep=1000, pressure=0.5 /' | build/mc_npt
Now the volume move acceptance rate should be around 60-70%. If not, check your formulae, and the v_move routine, again.
The following cell visualizes the original configuration, which we are calling old_conf, as argon atoms, using weas-widget. Then look at the final configuration in the current working directory, by replacing old_conf with new_conf below. You should see that the final box is larger, and hence the final density is smaller, than the initial one. Therefore the simulation pressure \(P=0.5\) does not correspond to the density
\(\rho=0.5\). You can adjust weas_shape (defined in the Preliminaries section above), or try the widget’s fullscreen button ⛶ (may not work in Colab) if the viewer size does not suit your display, and you may wish to do “Open Controls/Atoms/Radius Type” to switch to van der Waals radii.
[ ]:
atoms=read(old_conf)
print('Cell dimensions (reduced units)',atoms.cell.lengths())
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
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 attributes and datasets from the HDF5 output file.
[ ]:
attr, data = read_hdf5_file('mc_npt.hdf5')
We give names to some of the important attributes, and print them out.
[ ]:
print(attr['Title'].astype(str))
print('Run steps',attr['nstep'])
N = attr['N']
P = attr['P']
T = attr['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=0.5\). 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.25-0.30. If it is very different, once more you may need to check your
formulae and the v_move routine again.
[ ]:
rho = N/data['V']
discard = 200
rhoavg = rho[discard:].mean()
print(f'Simulation average rho = {rhoavg:10.4f}')
fig, ax = plt.subplots()
ax.set_xlabel('step')
ax.set_ylabel(r'$\rho$')
ax.plot(rho);
ax.hlines(rhoavg,discard,rho.size,ls='dashed',color='C1');
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=0.5\).
[ ]:
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 which can be set in the &nml ... / namelist. 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.
[ ]:
new_conf = copy(old_conf,this_dir) # Copy supplied file into working directory
!echo '&nml pressure=1.3384/' | build/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.
[ ]:
attr, data = read_hdf5_file('mc_npt.hdf5')
[ ]:
print(attr['Title'].astype(str))
print('Run steps',attr['nstep'])
N = attr['N']
P = attr['P']
T = attr['T']
The first check is that the average density is close to the desired value \(\rho=0.5\).
[ ]:
rho = N/data['V']
rhoavg = rho.mean()
print(f'Simulation average rho = {rhoavg:10.4f}')
fig, ax = plt.subplots()
ax.set_xlabel('step')
ax.set_ylabel(r'$\rho$')
ax.plot(rho);
ax.axhline(rhoavg,ls='dashed',color='C1');
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.
[ ]:
Tavg = data['T'].mean()
print(f'Specified temperature T = {T:10.4f}')
print(f'Simulation average T = {Tavg:10.4f}')
Vavg = data['V'].mean()
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 correct expression for kappa_avg
kappa_avg = 0.0
kappa_eos = eos_fit['kappa_T']
print(f'Simulation average kappa_T = {kappa_avg:10.4f}')
print(f'Fitted EOS value kappa_T = {kappa_eos: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 \(p(U)\) and volume \(p(V)\) probability distributions are plotted next. We could use \(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 \(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 \(p(U,V)\). We will not pursue this, but for interest’s sake, we plot a 2D histogram of \(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()
ax.set_xlabel(r'$V$')
ax.set_ylabel(r'$p(V)$')
ax.hist(data['V'],density=True,label='Simulation')
ax.axvline(Veos,c='C1',label='Expected average')
ax.legend();
[ ]:
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()
ax.set_xlabel(r'$U$')
ax.set_ylabel(r'$p(U)$')
ax.hist(data['U'],density=True,label='Simulation')
ax.axvline(Ueos,c='C1',label='Expected average')
ax.legend();
[ ]:
fig, ax = plt.subplots()
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');
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.
[ ]: