Marcus N. Bannermanm.campbellbannerman@abdn.ac.uk Your browser does not support SVG Access the slides here: www.marcusbannerman.co.uk/CCP5School

• Practical: Scripting
• Statistical mechanics 2
• Monte Carlo 1
• Monte Carlo 2
• Monte Carlo 3
• Extra files etc.

Scripting introduction

There are many compiled programming languages out there: C++, C, Fortran, .... These are the heavyweight languages which are great for learning how to program (i.e. type safe) and are fast. But they are designed more for the computer, not so much for you.

The interpreted languages:python, ruby, fortran, erlang, perl, gawk, bash, Java(ish)/javascript, go, rust, matlab ... are lightweight, designed to be as simple as possible for you; however, many are application specific and slow.

Python is general, powerful, and popular...

xkcd.com/353/

Hello world python script

#!/usr/bin/python

		#Above is the shebang (#!)
		#This tells the shell that this text file is a script!
		#(and what should be used to interpret it)

		#P.s. the interpreter ignores anything after a hashtag

		print "Hello, world!"

We'll create and run this now...

That should have been the quickest program you've ever seen (programming + run time).

Scripting is about fast programming, not fast execution, and this is quickest approach overall for many applications.

Python has a lot of built-in functionality and libraries for rapid implementation...

Gaussian/Normally distributed random number generator

#!/usr/bin/python

		import numpy.random

		for i in range(100):
		print numpy.random.normal(0.0, 1.0)

The examples in the practical sessions (gauss.f90 or gauss.py) use a box-muller transform to generate gaussian variables to remain comparable.

Writing data to a file

#!/usr/bin/python

		import numpy.random

		output_file = open("junk.dat", "w")
		for i in range(100):
		y = numpy.random.normal(0.0, 1.0)
		output_file.write(str(i)+' '+str(y))

Reading is just as easy, you can google how to do this.

Plotting

#!/usr/bin/python

		import numpy.random

		y_data = []
		for i in range(100):
		y_data.append(numpy.random.normal(0.0, 1.0))

		print y_data

		import pylab

		pylab.plot(y_data)
		pylab.show()

The plotting features of Python alone make it worthwhile learning.

Data manipulation and plotting

#!/usr/bin/python

		import numpy
		import pylab

		y_data = []
		for i in range(1000):
		y_data.append(numpy.random.normal(0.0, 1.0))

		n, bins, patches = pylab.hist(y_data, 50, normed=1, histtype='step')

		x_data = numpy.arange(-5.0, 5.0, 0.01)
		y_data = [numpy.exp(-0.5*x**2/(1.0*1.0))/numpy.sqrt(2.0*numpy.pi*1.0**2) for x in x_data]

		pylab.plot(x_data, y_data)
		pylab.show()

The example (gauss.py) uses its own histogramming code, rather than pylab.hist() to remain comparable to the fortran code.

Command-line arguments

#!/usr/bin/python
		import sys
		from optparse import OptionParser

		# parse command line arguments
		usage = "usage: %prog [options] arg"
		parser = OptionParser(usage)
		parser.add_option("--n", type="int", dest="nn", default=1)
		(options, args) = parser.parse_args()

		import numpy
		import pylab

		y_data = []
		for i in range(options.nn):
		y_data.append(numpy.random.normal(0.0, 1.0))

		n, bins, patches = pylab.hist(y_data, 50, normed=1, histtype='step')

		x_data = numpy.arange(-5.0, 5.0, 0.01)
		y_data = [numpy.exp(-0.5*x**2/(1.0*1.0))/numpy.sqrt(2.0*numpy.pi*1.0**2) for x in x_data]

		pylab.plot(x_data, y_data)
		pylab.show()

Command line arguments turn scripts into controllable programs, here, the last example has a controllable number of iterations.

Running other programs

#!/usr/bin/python
		import os

		#Run the ls program
		os.system("ls")

Scripts are particularly powerful when used to "control" other programs and analyse their output.

Molecular dynamics

#!/usr/bin/python
		import numpy, pylab, math
		box_length = 10
		N = 20
		dt = 0.1
		k_particle = 5
		diameter = 2.0
		k_domain = 0.1
		pos = (numpy.random.rand(N, 2) - 0.5) * box_length
		vel = numpy.zeros((N,2))

		def force(i, j):
		rij = pos[i] - pos[j]
		distance = numpy.linalg.norm(rij)
		if distance > diameter or i == j:
        return numpy.array([0,0])
	rij_norm = rij / distance
	return k_particle * (diameter - distance) * rij_norm
	
	while True:
	pylab.scatter(pos[:,0], pos[:,1], s=250)
	pylab.xlim(-box_length, +box_length)
	pylab.ylim(-box_length, +box_length)
	pylab.draw()
	pylab.pause(0.001)
	pylab.clf()
	
	for i in range(N):
        pos[i] += 0.5 * dt * vel[i]

	forces = numpy.zeros((N,2))
	for i in range(N):
        forces[i] += - k_domain * pos[i]
        for j in range(N):
        forces[i] += force(i,j)

	for i in range(N):
        vel[i] += dt * forces[i]
        
	for i in range(N):
        pos[i] += 0.5 * dt * vel[i]

This is a toy "spring" particles trapped in a spring well, using a symplectic integrator. Great example of python's power (and its lack of speed). Balance YOUR programming time against the COMPUTER run time when choosing scripting or compiled languages.

Python versions of the fortran codes are available for the first day's workshop if you want to take a look. Check the programming-exercises folder.

Good luck! Have fun! Ask if you need help!

Statistical Mechanics 2

Ludwig Boltzman, who spent much of his life studying statistical mechanics, died in 1906, by his own hand. Paul Ehrenfest, carrying on the work, died similarly in 1933. Now it is our turn to study statistical mechanics. Perhaps it will be wise to approach the subject cautiously.

(Opening lines of "States of Matter", by D.L. Goodstein).

Statistical mechanics is our bridge from the microscopic scale to the “macroscopic” thermodynamics.

We can connect microscopic details (positions, velocities) to more interesting macroscopic properties (i.e., internal energy, heat capacity).

To be able to perform Monte Carlo, we must understand the relative probabilities (and therefore feasibility) of different states or configurations of a microscopic system.

Thermodynamics observationally tells us what is possible (dS≥0), but statistical mechanics can quantify the probabilities in terms of the microscopic state.

It also proves that (dS<0) is possible, but that this is so astronomically rare it may be ignored.

Microscopics notation recap

• N spherical atoms of mass m
• position of atom α is rα
• momentum of atom α is pα
• point in phase space: Γ=(r1,…rN,p1,…,pN)
• pairwise additive potential u(|rα−rα′|)
• Hamiltonian: H(Γ)=U(r1,…rN)+K(p1,…,pN)U(r1,…rN)=12∑α≠α′u(|rα−rα′|)K(p1,…,pN)=∑αp2α2m

Phase space

Thermodynamics Recap

Basic concepts

Ergodic hypothesis: A system with fixed N, V, and E (U) is equally likely to be found in any of its Ω(N,V,E) microscopic states.

Consider two subsets of states, ΩA and ΩB. A system is more likely to be found in set A if ΩA>ΩB.

Therefore S(ΩA)>S(ΩB) as A is more likely, thus S must be a monotonically increasing function in Ω.

As states increase multiplicatively, yet entropy is additive, the relationship must be logarithmic. S(N,V,E)=kBlnΩ(N,V,E) where kB=1.3806503×10−23 J K−1 is the Boltzmann constant, present for historic reasons.

Density of states

• The density of states Ω(N,V,E) is the number of ways the system can have an energy E. In other words, it is the "volume" of phase space with an energy E.
• Volume element in phase space: dΓ=1N!dr1dp1h3⋯drNdpNh3 where h is Planck's constant.
• Mathematically, the density of states is given by Ω(N,V,E)=∫E<H(Γ)<E+δEdΓ where δE≪E.

Microcanonical ensemble (NVE)

Free energies (other ensembles)

Canonical ensemble (NVT)

P(E)∝Ω(E;Etot)Ω(E;Etot)=Ω(E)Ωsurr.(Etot−E) (Note: the N and V variables for the system and the surroundings are implicit)

Canonical ensemble

P(E,Etot)∝Ω(E)Ωsurr.(Etot−E)∝Ω(E)elnΩsurr.(Etot−E)

Performing a taylor expansion of lnΩsurr. around Etot. P(E,Etot)∝Ω(E)elnΩB(Etot)−E∂∂EtotlnΩB(Etot)+⋯

First-order is fine as E≪Etot. Note that ∂lnΩ(E)/∂E=β. At equilibrium the surroundings and system have the same temperature thus: P(E,Etot)∝Ω(E)elnΩB(Etot)−βEP(E)∝Ω(E)e−βE (If the surroundings are large, Etot is unimportant and the constant term lnΩB(Etot) cancels on normalisation)

Canonical ensemble

Canonical ensemble: Partition function

Factorization of the partition function

• The partition function can be factorized \begin{align*} Q(N,V,T) = \frac{1}{N!} \int d{\bf r}_1\cdots d{\bf r}_N e^{-\beta U({\bf r}_1,\dots{\bf r}_N)} \int \frac{d{\bf p}_1}{h^3}\cdots\frac{d{\bf p}_N}{h^3} e^{-\beta K({\bf p}_1,\dots,{\bf p}_N)} \end{align*}
• The integrals over the momenta can be performed exactly \begin{align*} \int\frac{d{\bf p}_1}{h^3}\cdots\frac{d{\bf p}_N}{h^3} e^{-\beta K({\bf p}_1,\dots,{\bf p}_N)} &= \left[\int\frac{d{\bf p}}{h^3} e^{-\frac{\beta p^2}{2m}} \right]^{N} \\ &= \left(\frac{2\pi m}{\beta h^2}\right)^{3N/2} = \Lambda^{-3N} \end{align*}
• The partition function is given in terms of a configurational integral and the de Broglie wavelength \Lambda\begin{equation*} Q(N,V,T) = \frac{1}{N!\Lambda^{3N}} \int d{\bf r}_1\cdots d{\bf r}_N e^{-\beta U({\bf r}_1,\dots{\bf r}_N)} = \frac{Z(N,V,T)}{N!\Lambda^{3N}} \end{equation*}

Average properties

• The probability of being at the phase point {\boldsymbol\Gamma} is \begin{equation*} {\mathcal P}({\boldsymbol\Gamma}) = \frac{e^{-\beta H({\boldsymbol\Gamma})}}{Q(N,V,\beta)} \end{equation*}
• The average value of a property {\mathcal A}({\boldsymbol\Gamma}) is \begin{align*} \langle{\mathcal A}\rangle &= \int d{\boldsymbol\Gamma} {\mathcal P}({\boldsymbol\Gamma}) {\mathcal A}({\boldsymbol\Gamma}) \\ &= \frac{1}{Q(N,V,\beta)}\int d{\boldsymbol\Gamma} e^{-\beta H({\boldsymbol\Gamma})} {\mathcal A}({\boldsymbol\Gamma}) \end{align*}
• If the property depends only on the position of the particles \begin{align*} \langle{\mathcal A}\rangle &= \frac{1}{Q(N,V,\beta)} \int d{\bf r}_1\cdots d{\bf r}_N e^{-\beta U({\bf r}_1,\dots,{\bf r}_N)} {\mathcal A}({\bf r}_1,\dots,{\bf r}_N) \end{align*}

Grand canonical ensemble (\mu V T)

Grand canonical ensemble

Isothermal-isobaric ensemble (NpT)

n-particle density

• The n-particle density \rho^{(n)}({\bf r}_1,\dots,{\bf r}_n) is defined as N!/(N-n)! times the probability of finding n particles in the element d{\bf r}_1\cdots d{\bf r}_n of coordinate space. \begin{align*} &\rho^{(n)}({\bf r}_1,\dots,{\bf r}_n) = \frac{N!}{(N-n)!} \frac{1}{Z(N,V,T)} \int d{\bf r}_{n+1}\cdots d{\bf r}_N\, e^{-\beta U({\bf r}_1,\dots,{\bf r}_N)} \\ &\int d{\bf r}_1\cdots d{\bf r}_n\,\rho^{(n)}({\bf r}_1,\dots,{\bf r}_n) = \frac{N!}{(N-n)!} \end{align*}
• This normalisation means that, for a homogeneous system, \rho^{(1)}({\bf r})=\rho=N/V\begin{align*} &\int d{\bf r}\, \rho^{(1)}({\bf r}) = N \qquad \longrightarrow \qquad \rho V = N \\ &\int d{\bf r}_1d{\bf r}_2\, \rho^{(2)}({\bf r}_1,{\bf r}_2) = N(N-1) \end{align*}

n-particle distribution function

• The n-particle distribution function g^{(n)}({\bf r}_1,\dots,{\bf r}_n) is defined as \begin{align*} g^{(n)}({\bf r}_1,\dots,{\bf r}_n) &=\frac{\rho^{(n)}({\bf r}_1,\dots,{\bf r}_n)} {\prod_{k=1}^n \rho^{(1)}({\bf r}_i)} \\ g^{(n)}({\bf r}_1,\dots,{\bf r}_n) &= \rho^{-n} \rho^{(n)}({\bf r}_1,\dots,{\bf r}_n) \qquad \mbox{for a homogeneous system} \end{align*}
• For an ideal gas (i.e., U=0): \begin{align*} &\rho^{(n)}({\bf r}_1,\dots,{\bf r}_n) = \frac{N!}{(N-n)!} \frac{\int d{\bf r}_{n+1}\cdots d{\bf r}_N}{\int d{\bf r}_{1}\cdots d{\bf r}_N} = \frac{N!}{(N-n)!} \frac{V^{N-n}}{V^N} \approx \rho^n \\ &g^{(n)}({\bf r}_1,\dots,{\bf r}_n) \approx 1 \\ &g^{(2)}({\bf r}_1,{\bf r}_2) = 1 - \frac{1}{N} \end{align*}

Radial distribution function

• n=2: pair density and pair distribution function in a homogeneous fluid \begin{align*} g^{(2)}({\bf r}_1,{\bf r}_2) &= \rho^{-2} \rho^{(2)}({\bf r}_1,{\bf r}_2) \\ g^{(2)}({\bf 0},{\bf r}_2-{\bf r}_1) &= \rho^{-2} \rho^{(2)}({\bf 0},{\bf r}_2-{\bf r}_1) \end{align*}
• Average density of particles at r given that a tagged particle is at the origin is \begin{equation*} \rho^{-1} \rho^{(2)}({\bf 0},{\bf r}) = \rho g^{(2)}({\bf 0},{\bf r}) \end{equation*}
• The pair distribution function in a homogeneous and isotropic fluid is the radial distribution function g(r), \begin{equation*} g(r) = g^{(2)}({\bf r}_1,{\bf r}_2) \end{equation*} where r = |{\bf r}_1-{\bf r}_2|

Radial distribution function

Radial distribution function

• Energy, virial, and compressibility equations: \begin{align*} U &= N \rho \int d{\bf r} u(r) g(r) \\ \frac{\beta p}{\rho} &= 1 - \frac{\beta\rho}{6} \int d{\bf r} w(r) g(r) \\ \rho k_BT \kappa_T &= 1 + \rho \int d{\bf r} [g(r)-1] \end{align*}
• "Pair" virial \begin{equation*} w(r) = {\bf r}\cdot \frac{\partial u}{\partial{\bf r}} \end{equation*}

Structure factor

• Fourier component \hat{\rho}({\bf k}) of the instantaneous single-particle density \begin{align*} \rho({\bf r}) &= \sum_\alpha \delta^d({\bf r}-{\bf r}_\alpha) = \int \frac{d{\bf k}}{(2\pi)^d} \hat{\rho}({\bf k}) e^{-i{\bf k}\cdot{\bf r}} \\ \hat{\rho}({\bf k}) &= \int d{\bf r} \rho({\bf r}) e^{i{\bf k}\cdot{\bf r}} \end{align*}
• Autocorrelation function is called the structure factor S(k)\begin{align*} \frac{1}{N}\langle\hat{\rho}({\bf k})\hat{\rho}(-{\bf k})\rangle &= 1 + \frac{1}{N}\left\lt \sum_{\alpha\ne\alpha'} e^{-i{\bf k}\cdot({\bf r}_\alpha-{\bf r}_{\alpha'})} \right\gt \\ S({\bf k}) &= 1 + \rho\int d{\bf r}[g({\bf r})-1] e^{-i{\bf k}\cdot{\bf r}} \end{align*}
• S({\bf k}) reflects density fluctuations and dictates scattering

Fourier decomposition

Fourier decomposition

Structure factor

• A peak at wavevector k signals density inhomogeneity on the lengthscale 2\pi/k
• k=0 limit is related to the isothermal compressibility \kappa_T\begin{equation*} S(0) = 1 + \rho \int d{\bf r} [g({\bf r})-1] = \rho k_BT \kappa_T \end{equation*}

Further reading

• D Chandler, Introduction to Modern Statistical Mechanics (1987).
• DA McQuarrie, Statistical Mechanics (2000).
• LE Reichl, A Modern Course in Statistical Physics (2009).
• JP Hansen and IR McDonald, Theory of Simple Liquids, 3ed (2006).

Monte Carlo: Part 1

God does not play dice. A. EinsteinGod does play dice with the universe. All the evidence points to him being an inveterate gambler, who throws the dice on every possible occasion. S. Hawking

You are the supreme power of your simulation, so you get to choose: \begin{align*} \langle{\mathcal A}\rangle &= \int d{\boldsymbol\Gamma}\,{\mathcal P}({\boldsymbol\Gamma})\, {\mathcal A}({\boldsymbol\Gamma}) & \text{(MC)}\\ &\text{or} &\text{(or)}\\ \langle{\mathcal A}\rangle&=\lim_{t_{obs.}\to\infty} \frac{1}{t_{obs.}}\int_0^{t_{obs.}}\mathcal A(t)\,{\rm d}t & \text{(MD)} \end{align*}

Monte Carlo Casino

Introduction

• Monte Carlo (MC) is a stochastic molecular-simulation technique. Sampling can be biased towards regions of specific interest.
• Not constrained by natural timescales (as in molecular dynamics).
• Can easily handle fluctuating numbers of particles (as in the grand canonical ensemble).
• No (physical) dynamic foundation.

Overview

• Introduction
• Practicalities
• Free energy barriers
• Advanced algorithms
• Further reading:
  • MP Allen and DJ Tildesley, Computer Simulation of Liquids (1987).
  • D Frenkel and B Smit, Understanding Molecular Simulation (2001).
  • DP Landau and K Binder, A Guide to Monte Carlo Simulations in Statistical Physics (2000).

Monte Carlo evaluation of integrals

• Uniform, probability distribution {\mathcal P}(x) along the interval [a,b]: \begin{equation*} {\mathcal P}(x) = \frac{1}{b-a} \left\{ \begin{array}{ll} 0 & \mbox{for $-\infty\lt x\lt a$} \\ 1 & \mbox{for $a \le x \le b$} \\ 0 & \mbox{for $b\lt x\lt \infty$} \\ \end{array} \right. \end{equation*}
• Randomly sample points x_1, x_2, \cdots from {\mathcal P}(x).
• The integral is related to the average of the function f: \begin{equation*} I = (b-a) \int_a^b dx\,{\mathcal P}(x) f(x) = (b-a) \langle f(x) \rangle \end{equation*}
• The average can be approximated by: \begin{equation*} I \approx \frac{(b-a)}{N} \sum_{k=1}^N f(x_k) \end{equation*}

Example

Convergence

Example: Convergence

Central limit theorem

Example: Uniform distribution

Example: Exponential distribution

Example: Galton Board

Multidimensional integrals

\begin{equation*} I = \int_{V} d{\bf x} f({\bf x}) \end{equation*}

• Uniform, probability distribution {\mathcal P}({\bf x}) within the region V:
• Randomly sample points {\bf x}_1, {\bf x}_2, \cdots from {\mathcal P}({\bf x}).
• The integral is related to the average of the function f: \begin{equation*} I = V \int_{V} d{\bf x}\, {\mathcal P}({\bf x}) f({\bf x}) = V \langle f({\bf x}) \rangle \end{equation*}
• The average can be approximated by: \begin{equation*} I \approx \frac{V}{N} \sum_{k=1}^N f({\bf x}_k) \end{equation*}

Importance sampling

Influence of importance sampling

Sampling from an arbitrary distribution

• Importance sampling requires the generation of numbers from an arbitrary probability distribution.
• One dimensional distributions can be sampled by sampling uniformly from the cumulative distribution function.
• In higher dimensions, sampling algorithms are only available for certain distributions (e.g., Gaussian distributions).
• The Metropolis method allows the sampling of arbitrary, multidimensional distributions. This relies on constructing a Markov process.

Markov process

• {\mathcal P}_\alpha^{(n)} probability that the state \alpha is occupied at step n.
• Transition probability W_{\alpha'\leftarrow\alpha} Probability that a system in state \alpha' will transition to state \alpha.
• Markov chain: \begin{equation*} {\mathcal P}_\alpha^{(n+1)} = \sum_{\alpha'} W_{\alpha\leftarrow\alpha'} {\mathcal P}_{\alpha'}^{(n)} \end{equation*}
• What guarantees a stationary, limiting distribution {\mathcal P}_\alpha?

Transition probabilities

• W_{\alpha'\leftarrow \alpha} is the probability that the system will transition to state \alpha', given that it is in state \alpha.
• The transition matrix must satisfy: \begin{equation*} \sum_{\alpha'} W_{\alpha'\leftarrow \alpha} = 1 \end{equation*}
• The choice of the transition matrix determines the properties of the Markov process.

Stationary distribution

We need to choose the transition matrix to ensure that the stationary distribution

• exists and is unique \begin{equation*} {\mathcal P}_\alpha = \sum_{\alpha'} W_{\alpha\leftarrow\alpha'} {\mathcal P}_{\alpha'} \end{equation*}
• is rapidly approached \begin{align*} \lim_{n\to\infty} {\mathcal P}_\alpha^{(n)} &= {\mathcal P}_\alpha \\ {\mathcal P}_\alpha^{(n)} &= \sum_{\alpha_1,\alpha_2,\dots,\alpha_n} W_{\alpha\leftarrow\alpha_1}W_{\alpha_1\leftarrow\alpha_2}\cdots W_{\alpha_{n-1}\leftarrow\alpha_n}{\mathcal P}_{\alpha_n}^{(0)} \\ &= \sum_{\alpha'} W_{\alpha\leftarrow\alpha'}^n{\mathcal P}_{\alpha'}^{(0)} \end{align*}

Detailed balance

• Physically, the condition of detailed balance requires that the probability of observing the system transition from state a \alpha to another state \alpha' is the same as observing it transition from \alpha' to \alpha.
• Mathematically, this is \begin{equation*} W_{\alpha\leftarrow \alpha'} {\mathcal P}_{\alpha'} = W_{\alpha'\leftarrow \alpha} {\mathcal P}_\alpha \end{equation*}
• If the transition probability satisfies the detailed balance condition, the Markov process has a unique, stationary limiting distribution.
• Detailed balance is not required for the Markov process to possess a unique, stationary limiting distribution.

Metropolis method

• The Metropolis method provides a manner to sample from any probability distribution {\mathcal P}_\alpha, using a Markov process.
• The transition probabilities are divided into two contributions: \begin{equation*} W_{\alpha'\leftarrow \alpha} = A_{\alpha'\leftarrow \alpha} g_{\alpha'\leftarrow \alpha} \end{equation*} where g_{\alpha'\leftarrow \alpha} is the probability of selecting a move to state \alpha', given the system is at state \alpha, and A_{\alpha'\leftarrow \alpha} is the probability of accepting the proposed move.
• The acceptance probability is given by \begin{equation*} A_{\alpha'\leftarrow\alpha} = \left\{ \begin{array}{ll} 1 & \mbox{if ${\mathcal P}_{\alpha'}\gt {\mathcal P}_\alpha$} \\ \frac{{\mathcal P}_{\alpha'}}{{\mathcal P}_\alpha} & \mbox{if ${\mathcal P}_{\alpha'}\lt {\mathcal P}_\alpha$} \end{array} \right. \end{equation*}
• If g_{\alpha'\leftarrow \alpha}=g_{\alpha\leftarrow \alpha'}, then the process satisfies detailed balance.

Metropolis method

Markov process: Continuous state space

• {\mathcal P}^{(n)}(x) probability that the state x is occupied at step n.
• Transition probability W(x'\leftarrow x) Probability that a system in state x will transition to state x'.
• Markov chain \begin{equation*} {\mathcal P}^{(n+1)}(x) = \int dx' W(x\leftarrow x') {\mathcal P}^{(n)}(x') \end{equation*}

Random number generators

von Neumann is simply trying to caution you against believing that random number generators (RNG) are actually random.

But GPU graphics programmers really do use a state of \sin! \begin{align*} x_i = {\rm fract}(43758.5453\sin(12.9898\,i)) \end{align*}

Pseudo-random number generators

• Typically, random deviates are available with a uniform distribution on [0,1).
• Be suspicious of intrinsic "ran" functions
• Linear congruential generators use the recursion \begin{equation*} X_{j+1} = a X_j + b \qquad \mod m \end{equation*}
• Written {\rm LCG}(m,a,b,X_0)
• Generates integers in range 0 \le X \le m
• Sequence repeats with period m AT BEST
• 32-bit: m \sim 2^{31} \sim 10^9, whereas 64-bit: m \sim 10^{18}

Tests for random number generators

• RNGs assessed by plotting d-dimensional vectors (X_{n+1},\dots,X_{n+d}) and looking for “planes” (which are bad)
• 3d plots: good (left); bad LCG(2^{31},65539,0,1) (right)
• In addition, simulate a finite-size system for which the properties are known exactly (e.g., 2D Ising model on square lattice)
• The gold-standard battery of tests are the Diehard tests.

Mersenne Twister

• Based on Mersenne prime 2^{19937}-1.
• Passes many tests for statistical randomness.
• Default random number generator for Python, MATLAB, Ruby, etc. Available in standard C++. Code is also freely available in many other languages (e.g., C, FORTRAN 77, FORTRAN 90, etc.).

Monte Carlo 2

Being secret, the work [radiation shielding for nuclear weapons] of von Neumann and Ulam required a code name. A colleague of von Neumann and Ulam, Nicholas Metropolis, suggested using the name Monte Carlo, which refers to the Monte Carlo Casino in Monaco where Ulam's uncle would borrow money from relatives to gamble.

Using lists of "truly random" random numbers was extremely slow, but von Neumann developed a way to calculate pseudorandom numbers, using the middle-square method. Though this method has been criticized as crude, von Neumann was aware of this: he justified it as being faster than any other method at his disposal, and also noted that when it went awry it did so obviously, unlike methods that could be subtly incorrect.Taken from the MC wikipedia page

Overview

• Monte Carlo in the NVT ensemble
• Practicalities
• Ergodicity and free-energy barriers
• Measuring ensemble averages
• Monte Carlo in the isobaric-isothermal ensemble
• Monte Carlo in the grand canonical ensemble
• Example: Truncated, shifted Lennard-Jones fluid
• Finite-size effects

Canonical ensemble

• Constant particle number, volume, and temperature; fluctuating energy.
• Distribution of energies: \begin{equation*} {\mathcal P}(E;N,V,\beta) \propto \Omega(N,V,E)\,e^{ - \beta E} \end{equation*}
• Distribution of states: \begin{equation*} {\mathcal P}({\boldsymbol\Gamma};N,V,\beta) \propto e^{ - \beta H({\boldsymbol\Gamma})} \end{equation*}
• Partition function: \begin{align*} Q(N,V,\beta) = \int {\rm d}E\,\Omega(N,V,E)\,e^{- \beta E} = \int {\rm d}{\boldsymbol\Gamma}\,e^{- \beta H({\boldsymbol\Gamma})} \end{align*}
• Helmholtz free energy \begin{equation*} F(N,V,T) = -k_B\,T\ln Q(N,V,\beta) \end{equation*}

Metropolis Monte Carlo in the NVT ensemble

Boltzmann distribution: {\mathcal P}({\boldsymbol\Gamma}) \propto e^{-\beta\,H({\boldsymbol\Gamma})} Monte Carlo integration \begin{align*} \langle {\mathcal A} \rangle &= \int d{\boldsymbol\Gamma}\, {\mathcal P}({\boldsymbol\Gamma}) {\mathcal A}({\boldsymbol\Gamma}) \approx \frac{1}{{\mathcal N}} \sum_{k=1}^{\mathcal N} {\mathcal A}({\boldsymbol\Gamma}_k) \end{align*} (provided {\boldsymbol\Gamma}_k is generated from {\mathcal P}({\boldsymbol\Gamma}))

Transition probability: \begin{equation*} W({\boldsymbol\Gamma}_n\leftarrow{\boldsymbol\Gamma}_o) = g({\boldsymbol\Gamma}_n\leftarrow{\boldsymbol\Gamma}_o) A({\boldsymbol\Gamma}_n\leftarrow{\boldsymbol\Gamma}_o) \end{equation*} where \begin{align*} g({\boldsymbol\Gamma}_n\leftarrow{\boldsymbol\Gamma}_o) &= \mbox{random particle move} \\ A ({\boldsymbol\Gamma}_n\leftarrow{\boldsymbol\Gamma}_o) &= \left\{ \begin{array}{ll} 1 & \mbox{if $H({\boldsymbol\Gamma}_n) \lt H({\boldsymbol\Gamma}_o)$} \\ e^{-\beta(H({\boldsymbol\Gamma}')-H({\boldsymbol\Gamma}))} & \mbox{if $H({\boldsymbol\Gamma}_n) \gt H({\boldsymbol\Gamma}_o)$} \end{array} \right. \end{align*}

Particle move

• Pick particle at random (to obey detailed balance)
• Alter x, y, and z coordinates by displacements randomly from the interval -\Delta{\bf r}_{\rm max}\lt \Delta{\bf r}\lt \Delta{\bf r}_{\rm max}
• Calculate resulting change in potential energy \Delta U
• Accept move with probability {\rm min}(1,e^{-\beta\Delta U})
  • If \Delta U\lt 0, accept.
  • If \Delta U\gt 0, generate random number R \in[0,1)
  • If R\le e^{-\beta\Delta U}, accept.
  • If R\gt e^{-\beta\Delta U}, reject and retain old configuration.

Tuning the acceptance rate

• Acceptance rate controlled by maximum displacement(s).
• Displacements result in energy change \Delta U.
• Large displacements are unlikely to be accepted.
• Small displacements are more likely to be accepted, but progress through configuration space is slow.
• This is not always a sensible choice!
• Adjust trial displacements to give maximum actual displacement per unit CPU time.

Energy autocorrelation function

Equilibration

• MC simulations converge towards equilibrium (or “thermalize”) until p(A,\tau+1) = p(A,\tau)
• Starting from a non-equilibrium configuration, the deviation from equilibrium decreases like \exp(-\tau/\tau_0), where \tau_0 is the correlation “time”
• Must only retain measurements (for ensemble averages etc.) when \tau\gg\tau_0
• See the exercises and the Molecular Dynamics simulation lectures for more information.

Diagnostics

• Useful diagnostics include binning, autocorrelations, hot and cold starts, and configurational temperature
• Binning: accumulate separate averages over blocks of (say) 1000 MC cycles, and watch for convergence of block averages
• Autocorrelations: calculate correlation function \langle A(\tau)A(0)\rangle \sim \exp(-\tau/\tau_0) for some observable A, estimate \tau_0, discard configurations for t\lt 10\,\tau_0
• “Hot/cold starts”: check that simulations starting from different configurations (e.g., disordered and crystalline) converge to the same equilibrium state

Configurational temperature

• In MD, one can use the equipartition theorem to check the (kinetic) temperature (although this can be misleading!)
• DPD potentials: MP Allen, J. Phys. Chem. B 110, 3823 (2006).
• In MC, configurational temperature provides a sensitive test of whether you are sampling the correct Boltzmann distribution for the prescribed value of T\begin{align*} k_BT_{\rm conf} &= -\frac{\langle\sum_{\alpha}{\bf F}_{\alpha}\cdot{\bf F}_{\alpha}\rangle} {\langle\sum_{\alpha}\nabla_{{\bf r}_{\alpha}}\cdot{\bf F}_{\alpha}\rangle} + O(N^{-1}) \\ {\bf F}_\alpha &= \sum_{\alpha'\ne\alpha} {\bf F}_{\alpha\alpha'} = -\sum_{\alpha'\ne\alpha} \nabla_{{\bf r}_{\alpha}} u({\bf r}_{\alpha}-{\bf r}_{\alpha'}) \end{align*}
  • H. H. Rugh, Phys. Rev. Lett. 78, 772 (1997)
  • O. G. Jepps et al., Phys. Rev. E 62, 4757 (2000)

Configurational temperature: Example

• T_{\rm conf} is sensitive to programming errors
• BD Butler et al., J. Chem. Phys. 109, 6519 (1998).
• Orientations: AA Chialvo et al., J. Chem. Phys. 114, 6514 (2004).

Ergodicity

• “All areas of configuration space are accessible from any starting point.”
• Difficult to prove a priori that an MC algorithm is ergodic.
• Many interesting systems possess configuration spaces with “traps”.
• This can render many elementary MC algorithms practically non-ergodic (i.e. it is very difficult to access certain (potentially important) regions of configuration space).
• Common examples of apparent non-ergodicity arise from free-energy barriers (e.g., in first-order phase transitions).

Ergodicity

Isobaric-isothermal ensemble

• Constant number of particles, pressure, and temperature; fluctuating volume and energy.
• Joint distribution of volume and energy: \begin{equation*} {\mathcal P}(V,E;p,\beta) \propto \Omega(N,V,E) e^{\beta p V - \beta E} \end{equation*}
• Partition function: \begin{align*} \Delta(N,p,\beta) = \int dV \int dE \Omega(N,V,E) e^{\beta pV - \beta E} \end{align*}
• Gibb's free energy \begin{equation*} G(N,p,\beta) = - k_BT \ln \Delta(N,p,\beta) \end{equation*}

MC simulations in the isobaric-isothermal ensemble

• Single-particle move (results in U\to U+\Delta U) \begin{align*} A(o\to n) = \min(1,e^{-\beta\Delta U}) \end{align*}
• Volume move V\to V+\Delta V (-\Delta V_{\rm max}\le\Delta V \le\Delta V_{\rm max} results in U\to U+\Delta U) \begin{align*} A(o\to n) &= \min(1,e^{-\beta\Delta U}) \\ \frac{p(V+\Delta V)}{p(V)} &= \frac{(V+\Delta V)^N}{V^N} e^{-\beta(\Delta U + p\Delta V)} \end{align*}
• One MC sweep consists of N_{\rm sp} attempted single-particle moves, and N_{\rm vol} attempted volume moves.
• Single-particle moves and volume moves must be attempted at random with a fixed probability (and not cycled) so that A(o\to n) = A(n\to o)

Isobaric-isothermal ensemble

Grand canonical ensemble

• Constant chemical potential, volume, and temperature; fluctuating particle number and energy.
• Joint distribution of particle number and energy: \begin{equation*} {\mathcal P}(N,E;\mu,\beta) \propto \Omega(N,V,E) e^{\beta \mu N - \beta E} \end{equation*}
• Partition function: \begin{align*} Z_{G}(\mu,V,\beta) = \sum_N \int dE \Omega(N,V,E) e^{\beta \mu N - \beta E} \end{align*}
• Free energy: pressure \begin{equation*} pV = k_BT \ln Z_G(\mu,V,\beta) + \mbox{constant} \end{equation*}

MC simulations in the grand canonical ensemble

• Single particle move (results in U\to U+\Delta U) as before.
• Particle insertion: N\to N+1 (results in U\to U+\Delta U) \begin{align*} A(n\leftarrow o) &= \min\left(1,\frac{{\mathcal P}(N+1)}{{\mathcal P}(N)}\right) \\ \frac{{\mathcal P}(N+1)}{{\mathcal P}(N)} &= \frac{(zV)^{N+1}e^{-\beta(U+\Delta U)}}{(N+1)!} \frac{N!}{(zV)^{N}e^{-\beta U}} = \frac{zV e^{-\beta \Delta U}}{N+1} \end{align*} where z=e^{\beta\,\mu-\ln\Lambda^3} (\Lambda^3 arises from the new particle's velocity integral in its unadded state (see Frenkel and Smit), N+1 is to maintain detailed balance on selecting the added particle for deletion).
• Particle deletion: N\to N-1 (results in U\to U+\Delta U) \begin{align*} A(n\leftarrow o) &= \min\left(1,\frac{{\mathcal P}(N-1)}{{\mathcal P}(N)}\right) \\ \frac{{\mathcal P}(N-1)}{{\mathcal P}(N)} &= \frac{(zV)^{N-1}e^{-\beta(U+\Delta U)}}{(N-1)!} \frac{N!}{(zV)^{N}e^{-\beta U}} = \frac{N e^{-\beta \Delta U}}{zV} \end{align*}

Ensemble averages

• Finally(!) assuming that you are sampling the equilibrium distribution, ensemble averages are calculated as simple unweighted averages over configurations \begin{equation*} \langle {\mathcal A} \rangle = \frac{1}{\mathcal N} \sum_{k} {\mathcal A}({\boldsymbol\Gamma}_k) \end{equation*}
• Statistical errors can be estimated assuming that block averages are statistically uncorrelated

Properties

• pressure \displaystyle p=\left\lt \rho k_BT - \frac{1}{3V}\sum_{j\lt k}{\bf r}_{jk}\cdot \frac{\partial u_{jk}}{\partial {\bf r}_{jk}}\right\gt
• density \displaystyle \rho=\left\lt \frac{N}{V}\right\gt
• chemical potential \displaystyle \beta\mu-\ln\Lambda^3 =-\ln\left\lt \frac{V\exp(-\beta\Delta U_+)}{N+1}\right\gt
• heat capacity \displaystyle C_V/k_B = \frac{3}{2} N + \beta^2(\langle U^2\rangle-\langle U\rangle^2) \displaystyle C_p/k_B = \frac{3}{2} N + \beta^2(\langle H^2\rangle-\langle H\rangle^2)
• isothermal compressibility \displaystyle \kappa_T = -\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_T = \frac{\langle V^2\rangle-\langle V\rangle^2} {k_BT\langle V\rangle} = \frac{V(\langle N^2\rangle-\langle N\rangle^2} {k_BT\langle N\rangle^2}

Truncated-shifted Lennard-Jones atoms

\begin{align*} u_{\rm LJ}(r) &= 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] \\ u(r) &= \left\{ \begin{array}{ll} u_{\rm LJ}(r) - u_{\rm LJ}(r_c) & \mbox{for $r\lt r_c$} \\ 0 & \mbox{for $r\gt r_c$} \\ \end{array} \right. \end{align*}

Results: Lennard-Jones fluid

Results: Lennard-Jones fluid

Results: Lennard-Jones fluid

Regardless of the ensemble (how we specify the thermodynamic state), at the same state the average properties are nearly identical.

Ensembles

• canonical ensemble (NVT)
  • single-phase properties (at fixed densities)
  • phase transition with change in volume (e.g., freezing)
  • phase transition with change in structure (e.g., solid-solid)
• isobaric-isothermal ensemble (NpT)
  • Corresponds to normal experimental conditions
  • Easy to measure equation of state
  • Density histograms {\mathcal P}(V) at phase coexistence
  • Structural phase transitions
  • Expensive volume moves (recalculation of total energy)
• grand canonical ensemble (\muVT)
  • Open systems (mixtures, confined fluids, interfaces)
  • Vapor-liquid coexistence and critical phenomena
  • Dense systems (low insertion/deletion rate)

Finite size effects

• The finite size of the system (box dimension L, number of particles N) is a limitation for all computer simulations, irrespective of periodic boundary conditions.
• Finite-size effects become apparent when the interaction range or the correlation length \xi becomes comparable to L\begin{equation*} h(r) = g(r)-1 \sim \frac{e^{-r/\xi(T)}}{r} \end{equation*}
• Systems with long-range interactions (e.g., Coulomb)
• Vapour-liquid critical point: \xi(T) diverges like |T-Tc|^{-0.63}.

Monte Carlo 3: Biased sampling

Introduction

• Histogram methods
• Quasi non-ergodicity
• Vapor-liquid transition
• Gibbs ensemble MC simulations
• Free-energy barrier in the grand-canonical ensemble
• Multicanonical simulations
• Replica exchange

Histogram extrapolation

• Consider an NVT simulation at \beta_0, where we collect an energy histogram, {\mathcal H}(E;\beta_0): \begin{align*} {\mathcal P}(E;\beta_0) &= \frac{\Omega(N,V,E)}{Q(N,V,\beta_0)} e^{-\beta_0 E} \\ {\mathcal H}(E;\beta_0) &\propto \Omega(N,V,E) e^{-\beta_0 E} \end{align*}
• Estimate for the density of states: \begin{align*} \Omega(N,V,E) \propto {\mathcal H}(E;\beta_0) e^{\beta_0 E} \end{align*}
• Using the estimate for \Omega(N,V,E), we can estimate the histogram at any other \beta\begin{align*} {\mathcal P}(E;\beta) &= \frac{\Omega(N,V,E)}{Q(N,V,\beta)} e^{-\beta E} \\ {\mathcal H}(E;\beta) &\propto \Omega(N,V,E) e^{-\beta E} \\ &\propto {\mathcal H}(E;\beta_0) e^{-(\beta-\beta_0) E} \end{align*}

Example: LJ fluid \rho\sigma^3=0.5

Solid lines are simulations, dashed lines are extrapolated histograms from k_B\,T=2\varepsilon.

Example: LJ fluid \rho\sigma^3=0.5

The extrapolated histograms allow evaluation of properties at any/all temperatures (within sampled states).

Histogram extrapolation: Other properties

• Other properties (i.e., X) can be extrapolated by collecting the joint probability distribution at \beta_0\begin{align*} {\mathcal P}(X,E;\beta_0) &= \frac{\Omega(N,V,E,X)}{Q(N,V,\beta_0)} e^{-\beta_0 E} \\ {\mathcal H}(X,E;\beta_0) &\propto \Omega(N,V,E,X) e^{-\beta_0 E} \end{align*}
• Estimate for the modified density of states: \begin{align*} \Omega(N,V,E,X) \propto {\mathcal H}(X,E;\beta_0) e^{\beta_0 E} \end{align*}
• Using the estimate for \Omega(N,V,E), we can estimate the histogram of property X at any other \beta\begin{align*} {\mathcal H}(X,E;\beta) &\propto \Omega(N,V,E,X) e^{-\beta E} \\ &\propto {\mathcal H}(X,E;\beta_0) e^{-(\beta-\beta_0) E} \end{align*}

Example: LJ fluid

Extrapolation of the virial (i.e., excess) contribution to the pressure.

Histogram combining

• Consider the case where we perform NVT simulations at several temperatures \beta_1, \beta_2,..., \beta_n where we collect the histograms: \begin{align*} {\mathcal H}(X,E;\beta_k) &\propto \Omega(N,V,E,X) e^{-\beta_k E} \end{align*}
• Estimate for the density of states as a weighted sum: \begin{align*} \Omega(N,V,E,X) \propto \sum_k \frac{w_k}{Z_k} {\mathcal H}(X,E;\beta_k) e^{\beta_k E} \end{align*} where \sum_k w_k=1 and Z_k=\sum_{E,X} \Omega(N,V,E,X) e^{-\beta_k E}.
• The uncertainty in {\mathcal H}(X,E;\beta_k) is roughly [{\mathcal H}(X,E;\beta_k)]^{1/2} (Poisson).
• Determining the weights w_k by minimizing the uncertainty of the estimate of the density of states leads to the following: \begin{align*} \Omega(N,V,E,X)&=\frac{\sum_j {\mathcal H}(E,X,\beta_j)}{\sum_k N_k\,Z_k^{-1}\,e^{-\beta_k\,E}} \end{align*} where N_k is the number of samples, now solve for Z^{-1}_k!

Quasi non-ergodicity

• An MC algorithm may be theoretically ergodic, but in some cases it can be very difficult (or impossible) to sample all of the important regions of configuration space.
• Even in “simple” systems (e.g., mono-atomic fluids) significant free-energy barriers can separate these important areas.
• A common example is the barrier between simulated coexisting phases at a first-order phase transition.
• We concentrate on the vapour-liquid phase transition.

Vapor-liquid transition

• At coexistence T_{\rm vap}=T_{\rm liq}, p_{\rm vap}=p_{\rm liq}, and \mu_{\rm vap}=\mu_{\rm liq}
• The density is the order parameter

Vapor-liquid transition: Difficulties

• In a finite-size system, the interfacial free energy (positive, unfavourable) is significant
• Free energy: F=F_{\rm int}+F_{\rm bulk}=\gamma L^2 - AL^3
• Direct simulation will not locate the coexistence densities accurately

Gibbs ensemble Monte Carlo

• Two simulation boxes, I and II
• Total number of particles N, volume V, temperature T
• Separate single-particle moves (thermal equilibrium)
• Exchange volume such that V = V_I + V_{II} (equates p)
• Exchange particles such that N = N_I + N_{II} (equates \mu)
  \longleftrightarrow
• Phase coexistence in a single simulation!
• Not good for dense phases: low probability for simultaneous particle insertion (in one box) and deletion (from the other)

• AZ Panagiotopoulos, Mol. Phys. 61, 813 (1987).
• D Frenkel and B Smit, Understanding Molecular Simulation (2001).

Gibbs ensemble Monte Carlo

Phase diagram of square-well fluids

Phase diagram of square-well fluids

MC simulations in the GC ensemble

• Single particle move (results in U\to U+\Delta U) \begin{align*} A(n\leftarrow o) = \min(1,e^{-\beta\Delta U}) \end{align*}
• Particle insertion: N\to N+1 (results in U\to U+\Delta U) \begin{align*} A(n\leftarrow o) &= \min\left(1,\frac{{\mathcal P}(N+1)}{{\mathcal P}(N)}\right) \\ \frac{{\mathcal P}(N+1)}{{\mathcal P}(N)} &= \frac{(zV)^{N+1}e^{-\beta(U+\Delta U)}}{(N+1)!} \frac{N!}{(zV)^{N}e^{-\beta U}} = \frac{zV e^{-\beta \Delta U}}{N+1} \end{align*}
• Particle deletion: N\to N-1 (results in U\to U+\Delta U) \begin{align*} A(n\leftarrow o) &= \min\left(1,\frac{{\mathcal P}(N-1)}{{\mathcal P}(N)}\right) \\ \frac{{\mathcal P}(N-1)}{{\mathcal P}(N)} &= \frac{(zV)^{N-1}e^{-\beta(U+\Delta U)}}{(N-1)!} \frac{N!}{(zV)^{N}e^{-\beta U}} = \frac{N e^{-\beta \Delta U}}{zV} \end{align*}

Free energy barrier in GC ensemble

• In the grand-canonical ensemble \begin{equation*} {\mathcal P}(N) = \frac{Q(N,V,T)}{Z_G(\mu,V,T)}e^{\beta\mu N } = \frac{Q(N,V,T)}{Z_G(\mu,V,T)} z^N \end{equation*}
• In a system at coexistence well below T_c, the particle number histogram p(N) [or p(\rho=N/V)] at constant chemical potential, temperature, and volume is bimodal, with almost no “overlap” between the peaks

Multicanonical simulations

• Apply a biasing potential \eta(N) in the Hamiltonian to cancel out the barrier. \begin{equation*} {\mathcal P}_{\rm bias}(N) \propto {\mathcal P}(N) e^{-\eta(N)} \end{equation*}
• Insertion/deletion acceptance probabilities are modifed \begin{equation*} \frac{{\mathcal P}_{\rm bias}(N_n)}{{\mathcal P}_{\rm bias}(N_o)} = \frac{{\mathcal P}(N_n)}{{\mathcal P}(N_o)} e^{-[\eta(N_n)-\eta(N_o)]} \end{equation*}
• Single-particle moves are unaffected

Multicanonical weights

• The “ideal” choice for \eta(N) would cancel out the barrier completely such that the biased probability distribution is flat, i.e., \eta(N) = k_BT\ln{\mathcal P}(N) so that {\mathcal P}_{\rm bias}(N) = {\rm constant}
• However, if we knew that in advance, there would be no need for a simulation!
• Fortunately, \eta(N) can be determined iteratively, to give successively smooth biased distributions

Determining the multicanonical weights

• [Step 1:] Choose z to be near coexistence (selected from hysteresis region)
• [Step 2:] Get good statistics in {\mathcal P}(N) for N_{\rm min}\le N\le N_{\rm max}
• [Step 3:] Generate biasing potential \begin{align*} \eta(N) &= k_BT\ln{\mathcal P}(N) \\ \eta(N\lt N_{\rm min}) &= k_BT\ln{\mathcal P}(N_{\rm min}) \\ \eta(N\gt N_{\rm max}) &= k_BT\ln{\mathcal P}(N_{\rm max}) \end{align*}
• [Step 4:] Get good statistics in {\mathcal P}_{\rm bias}(N) (now over a wider range of N)
• [Step 5:] Generate unbiased p(N)\begin{align*} {\mathcal P} \propto {\mathcal P}_{\rm bias}(N) e^{\eta(N)} \end{align*} return to Step 3

Multicanonical simulations: LJ fluid

GCMC: Histogram reweighting

• At reciprocal temperature \beta_0 and activity z_0\begin{align*} {\mathcal P}(N,E;\mu_0,V,\beta_0) = \frac{\Omega(N,V,E)}{Z_G(\mu_0,V,\beta_0)} e^{\beta_0\mu_0 N - \beta_0 U} \end{align*}
• Now consider a different activity and temperature \begin{align*} {\mathcal P}(N,E;\mu_1,V,\beta_1) &= \frac{\Omega(N,V,E)}{Z_G(\mu_1,V,\beta_1)} e^{\beta_1\mu_1 N-\beta_1 U} \\ &= {\mathcal P}(N,E;\mu_0,V,\beta_0) \left(\frac{z_1}{z_0}\right)^N e^{-(\beta_1-\beta_0)U} \\ & \qquad \times \frac{Z_{G}(\mu_1,V,\beta_1)}{Z_{G}(\mu_0,V,\beta_0)} \end{align*}
• The ratio Z_{G}(\mu_1,V,\beta_1)/Z_{G}(\mu_0,V,\beta_0) is just a normalizing factor, so that \begin{align*} \sum_{N} \int dE {\mathcal P}(N,E;\mu,V,\beta) = \sum_{N} {\mathcal P}(N;\mu,V,\beta) = \int dE {\mathcal P}(E;\mu,V,\beta) = 1 \end{align*}

Determining coexistence

• To find coexistence, tune \mu to satisfy the “equal area” rule (not equal peak height)
• These histograms were generated from a multicanonical simulation with z=0.029.

Density histogram for the LJ fluid

Lennard-Jones phase diagram

Replica exchange

• “Rough” energy landscapes are hard to sample at low temperature (get stuck in local minima)
• High-temperature simulations can glide over barriers
• Exchange complete configurations (with energies U_0 and U_1) between simulations run in parallel at different reciprocal temperatures (\beta_0 and \beta_1, respectively) \begin{equation*} \frac{{\mathcal P}(n)}{{\mathcal P}(o)} = \frac{e^{-\beta_0U_1} \times e^{-\beta_1U_0}}{e^{-\beta_0U_0} \times e^{-\beta_1U_1}} = e^{-(\beta_0-\beta_1 )(U_1-U_0 )} \end{equation*}

Replica exchange: Simple example

• \begin{align*} U(x) &= \left[\sin\frac{\pi x}{2} \sin 5\pi x\right]^2 \\ p(x,\beta) &\propto e^{-\beta U} \end{align*}

  Model one-dimensional system (after Frenkel and Smit)
• Single particle in unit box (with periodic boundary conditions)
• External potential, U(x).
• Start, x = 0.
• \Delta x = \pm0.005
• 10^6 MC moves