Monte Carlo of Lennard-Jones atoms¶
In this notebook you will use the Monte Carlo method to compute the properties of a system of Lennard-Jones atoms.
There is an accompanying notebook MC-integration.ipynb which discusses Monte Carlo sampling and Monte Carlo integration. The two notebooks are independent, so you can study them in either order, although the other one introduces more basic ideas, so we recommend you do that one first.
A third notebook, MC-Pressure.ipynb, focused on constant-pressure Monte Carlo simulation of Lennard-Jones atoms, will be the subject of a later workshop. It follows on from some of the material covered here.
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_nvt.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
Introduction¶
The Lennard-Jones potential is \begin{equation*} u_{\mathrm{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 adopted, so \(\varepsilon=1\) and \(\sigma=1\), and in addition Boltzmann’s constant is taken to be unity \(k_{\text{B}}=1\).
The system of interest here consists of atoms interacting through the cut-and-shifted Lennard-Jones potential defined by \begin{equation*} u(r) = \begin{cases} u_{\mathrm{LJ}}(r) - u_{\mathrm{LJ}}(r_{\text{cut}}) & r \leq r_{\text{cut}} \\ 0 & r> r_{\text{cut}} \end{cases} \end{equation*} where \(r_{\text{cut}}\) is the cutoff distance. The choice made here is \(r_{\text{cut}}= 2.5\sigma\). For this potential, the critical point is at \(T_{\text{c}}=1.0779\), \(\rho_{\text{c}}=0.3190\) (again in reduced units). This exercise concentrates on the supercritical state point, \(T=2.0\), \(\rho=0.5\).
Thermodynamic quantities, and the pair distribution function \(g(r)\), for this state point, were presented in lectures; in this exercise you will be comparing with those results. More generally, an accurate fitted equation of state for the fluid region has been developed by M Thol et al., Int J Thermophys, 36, 25 (2015). A Python program implementing their formulae is supplied in the file eos_lj.py. This can be run interactively from the
command line, in which case it will ask for the values of \(T\) and \(\rho\), and will print out a range of thermodynamic data for that state point. Instead, this worksheet has imported a function eos from eos_lj, which takes \(T\) and \(\rho\) as arguments, and returns various thermodynamic quantities in a dictionary. The next cell prints out values for this state point, which should match the values given in the lecture reasonably well. We shall compare several of these
fitted values with our simulations in this workshop.
[ ]:
eos_fit = eos(temperature=2.0,density=0.5)
for key, value in eos_fit.items():
print(f'{key:20s} {value:10.4f}')
Monte Carlo program¶
The instruction in the following cell should prepare to build the programs of interest to us here, in a directory build. We are using the Meson build system to do this. Normally, you should only need to do this setup once. However, if at some point you want a completely clean restart, you can delete the build directory and run the command again.
[ ]:
!meson setup build
The following cell should compile both the code of interest here, mc_nvt, and the constant-pressure code, mc_npt, which we will use in a later workshop. If you ever change the Fortran source files, you can compile again. If you want to see more details of the compilation, add the flag --verbose at the end of the line. If you want to clean the build directory, add the flag --clean at the end of the line.
[ ]:
!meson compile -C build
The next cell copies the supplied atomic configuration DATA/mc_nvt.xyz to the input file mc_nvt.xyz required by the mc_nvt program. Then it executes the program with default parameters. While it is running, read through the following description of some of the features of the program, looking at some of the program source files.
[ ]:
this_dir = Path.cwd() # Current working directory
data_dir = this_dir / 'DATA' # Data directory
old_conf = data_dir / 'mc_nvt.xyz' # Supplied configuration file
new_conf = copy(old_conf,this_dir) # Copy supplied file into the working directory
!echo '&nml /' | build/mc_nvt
The Fortran source files may be found in the src subdirectory.
Take a look at the main program file mc_nvt.f90. This carries out a Monte Carlo simulation of \(N\) atoms, in a fixed volume \(V\), at specified temperature \(T\). Simulation parameters are provided through standard input in the form of a namelist, a feature of Fortran that allows default values to be specified within the program, while letting you change them, using keywords, if you wish. So, a run using default parameters is initiated with echo '&nml /' | build/mc_nvt, but
you could choose a different maximum displacement by running with the command echo '&nml dr_max=0.3 /' | build/mc_nvt.
Have a look at the overall structure of the code in mc_nvt.f90: there is a loop over steps. Each step consists of an attempt to move \(N\) atoms; the appropriate routine is in the file mc_module.f90. In this case, the atoms are chosen randomly; you might like to consider whether this is the only valid approach. Also in that module is the routine for estimating the chemical potential by Widom test particle insertion.
The LJ potential details, including the value of \(r_{\text{cut}}\), are specified in the file potential_module.f90. The potential energy and virial functions are calculated here. You may also be interested in the routine which calculates the configurational temperature from the Laplacian and squared forces. These key properties, just mentioned, are stored at each step, and output at the end of the run to a file mc_nvt.hdf5, which you will read in shortly, for analysis. (We are using
HDF5 format for this, but there will be no need to look closely at the details). Standard output is just used for the crucial information needed to confirm that the program is running. At (increasing) intervals, the step number, CPU time consumed so far, and the cumulative move acceptance ratio are printed, as you will see in the preceding cell. At regular intervals gap steps, the program stores all the atomic positions, and these are output to the mc_nvt.hdf5 file at the end. These will
be used shortly to calculate the pair distribution function.
Think carefully about all the calculations done in the program, referring back to the lecture notes. If there is anything that seems unclear, feel free to ask! When you have finished examining the program source files, you can close them.
By now the run in the cell above should have finished; the statement Program ends will be printed, along with the CPU time taken by the run. Check this, and move on to the following cells. (Don’t start executing the program a second time!)
The program reads an input configuration, and writes an output configuration. This is handled by routines within config_io_module.f90, but there is no need to look at this file. An initial configuration of atoms mc_nvt.xyz was copied from the DATA subdirectory into the working directory, and overwritten at the end of the run. 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 (actually cubic here), 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 the supplied file, which we’re referring to as old_conf. To help weas-widget do this, we convert into units appropriate for argon atoms. If the viewer size does not suit your display, you can adjust weas_shape (defined in the Preliminaries section above), or try using the widget’s fullscreen button ⛶ (may not work in Colab). The initial view uses small spheres (based on the covalent
radius), but if you open the controls (top left of the widget) you can select “Atoms/Radius Type”, and switch to “VDW”, which gives a more realistic impression of the density.
[ ]:
atoms=read(old_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
Feel free to also look at the final configuration in the current working directory, by replacing old_conf by new_conf in the above cell.
Reading the HDF5 file¶
We are using a Hierarchical Data Format (HDF5) file mc_nvt.hdf5 to store output from the simulation. HDF5 is a very flexible format for storing data, but we are only using some simple features: specifically mc_nvt.hdf5 is a completely flat file containing just two kinds of object: attributes and datasets. The attributes are used to store a few simulation parameters, such as the number of particles, box lengths, and temperature. Step-by-step values are stored in the datasets. These
objects are key-value pairs, like Python dictionaries. The function read_hdf5_file is provided to read these in. It prints a list of all the keys, returns the attributes as a dictionary of values, which we store in attr, and returns the datasets as a dictionary of NumPy arrays, which we store in data.
[ ]:
attr, data = read_hdf5_file('mc_nvt.hdf5')
We store some of the attributes in named variables: the number of atoms N, volume V, temperature T, and array of box lengths L.
[ ]:
print(attr['Title'].astype(str))
print('Run steps',attr['nstep'])
N = attr['N']
V = attr['V']
T = attr['T']
L = attr['L']
print(f'Number of atoms N = {N:10d}')
print(f'Volume V = {V:10.4f}')
print(f'Temperature T = {T:10.4f}')
print(f'Box lengths L = {L[0]:10.4f}{L[1]:10.4f}{L[2]:10.4f}')
The dataset data['T'] contains the step-by-step values of the configurational temperature. It will be interesting to see if the average of this quantity agrees with the input temperature.
[ ]:
Tavg=data['T'].mean()
print(f'Simulation average T = {Tavg:10.4f}')
print(f'Specified value of T = {T:10.4f}')
Simulation Results¶
In the following we will look at some of the other quantities.
The activity is \(z=\exp(\mu/k_{\text{B}}T)\). However, the Widom test-particle insertion method gives us, after averaging, \(\exp(-\mu/k_{\text{B}}T)\), and this is what data['Z'] contains. So the simulation estimate of the activity \(z\) is the inverse of the average of data['Z']. Let’s see if it agrees with the value obtained from the fitted eos function. You can also compare with the value tabulated in the lecture notes for this state point.
[ ]:
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}')
Now: something for you to do! The average of the virial \(W\) can be used to give the simulation pressure. Referring to the Statistical Mechanics lecture notes, here is the formula. \begin{equation*} P = \frac{ N k_{\text{B}} T}{V} + \frac{\langle W\rangle}{3V} \end{equation*} Calculate this in the cell below. Bear in mind that all the variables you need have already been read from the file, and that \(k_{\text{B}}=1\) in our reduced simulation units. Compare with the value quoted in the lecture notes for this state point, and with the value given by the approximate equation of state, which is already in the cell.
[ ]:
# Insert the correct formula for Pavg here
Pavg = 0.0
Peos = eos_fit['P']
print(f'Simulation average P = {Pavg:10.4f}')
print(f'Fitted EOS value P = {Peos:10.4f}')
Hopefully the results are quite close. If we wanted, we could further analyze the step-by-step data for all these quantities to estimate the statistical error on the simulation averages, but this is not the topic of the current workshop.
Instead, we shall take a closer look at the potential energy \(U\) stored in data['U']. The next cell plots \(U\) as a function of step. You might like to plot a subset of the data, over fewer steps, to get an idea of how correlated successive values are.
[ ]:
fig, ax = plt.subplots()
ax.set_xlabel('step')
ax.set_ylabel(r'$U$')
ax.plot(data['U']);
The next cell plots a probability histogram \(p(U)\). Take a close look at this, making sure that it looks sensible, especially compared with the mean value and standard deviation of the data. We also compare with the fitted EOS values of \(U\) and \(u=U/N\).
[ ]:
Uavg = data['U'].mean()
Ustd = data['U'].std()
uavg = Uavg/N
ueos = eos_fit['u']
print(f'Simulation average U = {Uavg:10.4f}')
print(f'Standard deviation U = {Ustd:10.4f}')
print(f'Simulation average u = {uavg:10.4f}')
print(f'Fitted EOS value u = {ueos: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(N*eos_fit['u'],c='C1',label='Fitted EOS')
ax.legend();
Heat Capacity¶
Now: something more for you to do!
It should be possible to use this to estimate the heat capacity at constant volume \(C_V\), or \(c_V=C_V/N\) per atom. Referring to the lecture notes, here is the formula: \begin{equation*} c_V/k_{\text{B}} = C_V/N k_{\text{B}} = \frac{3}{2} + \frac{\langle U^2\rangle - \langle U\rangle^2}{N(k_{\text{B}} T)^2} . \end{equation*} Do this calculation in the cell below, to give \(c_V\) (heat capacity per atom) as calculated in your simulation. Remember, that \(k_{\text{B}}=1\) in our reduced units, and that the other required values were already read in above. Compare this \(c_V\) with the value quoted in the lecture for this state point, and also the value returned by the fitted EOS, which appears in the cell below.
[ ]:
# Insert the correct formula for cavg here
cavg = 0.0
ceos = eos_fit['c_V']
print(f'Simulation average c_V = {cavg:10.4f}')
print(f'Fitted EOS value c_V = {ceos:10.4f}')
Pair Distribution Function¶
Don’t worry if you run out of time before tackling the remaining topics: the pair distribution function and ensemble reweighting. You can always return to this notebook in later workshops (you should not need to re-run the simulation, and after you import the necessary Python modules and functions at the top, you may skip to this point).
The next cell re-reads the mc_nvt.hdf5 file (in case you are returning here afresh).
[ ]:
attr, data = read_hdf5_file('mc_nvt.hdf5')
Once more, we assign names to a few of the important attributes. Having checked that the box is cubic, L is redefined to be a scalar rather than an array.
[ ]:
print(attr['Title'].astype(str))
print('Run steps',attr['nstep'])
N = attr['N']
V = attr['V']
T = attr['T']
L = attr['L']
print(f'Number of atoms N = {N:10d}')
print(f'Volume V = {V:10.4f}')
print(f'Temperature T = {T:10.4f}')
assert np.allclose ( L, L[0] ), print('Error: we are assuming a cubic box')
L = L[0]
print(f'Box length L = {L:10.4f}')
The data['r'] dataset contains a set of configurations (atom positions). Interestingly, but not unexpectedly, the order of indices of the data['r'] array is reversed, compared with the Fortran order. Here, the step index comes first.
[ ]:
print(data['r'].shape) # Should be (nr,N,3)
nr = data['r'].shape[0]
We are going to calculate \(g(r)\) out to half the box length (it is a cubic box). This will involve accumulating a histogram of pair distances. For this simple example we’ll specify the number of bins, and this determines the bin width. We do the counting in a crude way, including both \(ji\) and \(ij\) for each pair. This is not necessarily the fastest approach, but the aim here is to show clearly what we are calculating. The loop over configurations should take a few seconds.
[ ]:
D = L/2 # Half the box length
b = 200 # Number of bins
h = np.zeros(b) # Histogram initialized to zero
for r in data['r']:
d = r[:,np.newaxis,:] - r[np.newaxis,:,:] # Set of all separation vectors (N,N,3)
d = np.fabs(d) # Absolute values of vector components
d = np.where(d<D,d,L-d) # Simple PBC in this case
d = np.sqrt(np.sum(d**2,axis=-1)) # Set of all separation distances (N,N)
h1,rr = np.histogram(d,bins=b,range=(0.0,D)) # Separation histogram & bin edges
h = h + h1 # Accumulate histogram
h[0] = 0 # Remove the counts arising from the diagonal of d
h = h / nr # Normalise by number of configurations
h = h / N # Normalise by number of atoms
Now stop to consider what we have calculated.
For each configuration, h1 counts the number of pair separations lying in the range corresponding to each histogram bin. We have counted each distinct pair twice (the \(N\times N\) matrix d is symmetric). This is equivalent to considering each atom \(i\) in turn: each bin of h1 counts the number of neighbours \(j\) whose distances from \(i\) fall into that bin, and the results are summed over all \(i\). The double counting corresponds to \(i\) being a
neighbour of \(j\) as well as \(j\) being a neighbour of \(i\). There are unwanted counts corresponding to the diagonal elements \(i=j\), but these can be removed when the loop is finished.
The h array simply sums these results for all nr stored configurations, and after normalizing by nr it contains the average number of pair separations lying in the range corresponding to each histogram bin. Again, think of this by considering each atom \(i\) in turn: each bin of h contains the average number of neighbours \(j\) whose distances from \(i\) fall into that bin, summed over all such atoms \(i\).
After further normalizing by N, each bin of h contains the average number of neighbours whose distance from any given atom lies in the range covered by that bin.
The rr array gives the bin edges, i.e. the separation values delimiting each bin. It has b+1 elements. We can use these values to compute the volume of the spherical shell corresponding to each bin. This in turn allows us to calculate what h would be in an ideal gas of the same density as our system. We proceed to do this, calling the result h_id: the ratio h/h_id is \(g(r)\). Before plotting, we use the rr array again to compute the mid-point (\(r\)-value) of
each bin. Have a look at the results. Does the \(g(r)\) plot look sensible? More importantly, have you understood the way we calculate \(g(r)\)? Feel free to ask if anything is unclear!
[ ]:
V_shell = np.diff((4.0*np.pi/3.0)*rr**3)
rho = N/V
h_id = rho * V_shell
g = h / h_id
r = ( rr[1:] + rr[:-1] ) / 2 # Mid-points of all the bins
fig, ax = plt.subplots()
ax.set_xlabel(r'$r$')
ax.set_ylabel(r'$g(r)$')
ax.plot(r,g)
ax.axhline(1,ls='dashed',c='C1');
Using the pair distribution function¶
This part is entirely optional, and can be skipped. An almost identical question appears in the MD tutorials: there is no need to answer both! In principle, the average potential energy may be calculated from \(g(r)\) using the formula \begin{equation*} u_{\text{avg}} = U/N = 2\pi\rho\int_0^{\infty} u(r) \, g(r) \, r^2 \, dr , \end{equation*} and a similar formula applies to the non-ideal contribution to the pressure. You can calculate this using a simple SciPy/NumPy numerical integration routine. Don’t forget that the interactions here are given by the cut-and-shifted LJ potential \(u(r)\); the next cell gives a suitable function for this. The integral ranges over \(0 \leq r \leq r_{\text{cut}}\) for this potential. How does the result compare with the average from the MC simulation (calculated in an earlier cell), and with the fitted EOS?
[ ]:
def u(r):
"""Lennard-Jones cut-and-shifted potential, r may be scalar or NumPy array."""
rc = 2.5 # Assumed cutoff distance
rc2 = 1.0 / rc**2
rc6 = rc2**3
uc = 4.0*(rc6-1.0)*rc6
r2 = 1.0 / r**2
r6 = r2**3
ulj = 4.0*(r6-1.0)*r6
return np.where ( r<rc, ulj-uc, 0.0 )
Use the above function, and the arrays g and r, in the next cell to calculate the desired integral by quadrature. For example, the np.trapezoid(f,r) function is provided by NumPy. This uses the trapezoidal rule, where r is an array containing the sample points, and f is an array containing values of the integrand evaluated at those points.
Again, we can compare with the fitted EOS.
[ ]:
# Insert code to calculate uavg
uavg = 0.0
eos_fit = eos(temperature=T,density=N/V)
ueos = eos_fit['u']
print(f'g(r) integral for u = {uavg:10.4f}')
print(f'Fitted EOS value u = {ueos:10.4f}')
Ensemble Reweighting¶
You may prefer to skip this, and come back to it after the third MC lecture, where we go into the topic in more detail.
For convenience, the next cell re-reads the HDF5 file. Then the probability histogram for the potential energy, \(p(U)\), is re-calculated.
[ ]:
attr, data = read_hdf5_file('mc_nvt.hdf5')
[ ]:
N = attr['N']
V = attr['V']
T = attr['T']
[ ]:
p, U = np.histogram(data['U'],bins=100,density=True)
dU = np.diff(U) # Get differences in bin edges, i.e. bin widths, for later use
U = (U[:-1]+U[1:])/2 # Convert bin edges into midpoint values
This is the distribution for the simulation run at \(T=2.0\). It is possible to use this data to estimate \(p(U)\), and hence calculate \(\langle U\rangle\) and other functions of potential energy, at nearby temperatures.
The distribution function may be expressed as \(p(U)\propto \Omega(U) \exp(-\beta U)\) where \(\Omega(U)\) is the density of states and \(\beta=1/k_{\text{B}}T\). From the same formula at a nearby temperature \(T_1\), it follows that the distribution \(p_1(U)\) at this nearby temperature is \begin{equation*} p_1(U) \propto p(U) \times \exp\bigl[(\beta-\beta_1) U\bigr] \end{equation*} where \(\beta_1=1/k_BT_1\). The function \(p_1(U)\) needs to be normalized, after calculating it this way, so that \(\int p_1(U) \, dU = 1\). The following cell attempts to do this for \(T_1=1.8\), using the just-calculated \(p(U)\) curve. It also compares the results with the fitted equation-of-state function at both temperatures, adding the expected values of \(\langle U\rangle\) as vertical dashed lines in the plot.
One subtlety, in general, is that the exponential function may produce very large or very small values, possibly leading to underflow and overflow issues. We should be dealing with float64 datasets, and hence the variables derived from them should also be float64. Provided \(\beta\) and \(\beta_1\) are not too different from one another, we expect this to be sufficient, but we print the normalization factors just to emphasize this danger.
[ ]:
beta = 1/T
T1 = 1.8
beta1 = 1/T1
# Confirm that p(U) is already normalized
norm = np.sum(p*dU) # This is the integral of p
print(f'p(U) norm = {norm:.5g}')
# Compute distribution at T1, and normalize
p1 = p * np.exp((beta-beta1)*U)
norm = np.sum(p1*dU) # This is the integral of p1
print(f'p1(U) norm = {norm:.5g}')
p1 = p1/norm
# Compute average values from original and scaled P(U) distributions
# Compare with expected mean values from fitted EOS
Uavg = np.average(U,weights=p)
U1avg = np.average(U,weights=p1)
Ueos = N * eos(temperature=T, density=N/V)['u']
U1eos = N * eos(temperature=T1,density=N/V)['u']
print(f'Uavg at T = {Uavg:10.1f}')
print(f'Ueos at T = {Ueos:10.1f}')
print(f'Uavg at T1 = {U1avg:10.1f}')
print(f'Ueos at T1 = {U1eos:10.1f}')
# Plot distributions
fig, ax = plt.subplots()
ax.set_xlabel(r'$U$')
ax.set_ylabel(r'$p(U)$')
line, = ax.plot(U,p,label='T={:5.2f}'.format(T))
ax.axvline(Ueos,ls='--',c=line.get_color())
line, = ax.plot(U,p1,label='T={:5.2f}'.format(T1))
ax.axvline(U1eos,ls='--',c=line.get_color())
ax.set_ylim(bottom=0)
ax.legend();
Now, experiment with \(T_1\). Try, say, \(T_1=2.2\), \(1.5\) or \(2.5\). What is the effect on \(p_1(U)\) as \(T_1\) gets further from \(T\)?
You may also like to compare with the results of a fresh MC run at temperature \(T_1\) (taking care to equilibrate at the new temperature first).
Further work¶
There’s plenty of scope to experiment with the mc_nvt program. For example, you might tinker with the parameter dr_max, observing its effect on the move acceptance ratio. Consider the question: how should the “optimal” value of dr_max be determined? What do we want to optimize?
This concludes the notebook.
[ ]: