Computational Physics

318 The Monte Carlo method

The quotients in this expression are canonical ensemble averages corresponding to

Z 0

Z 1

=

〈M[β(U 0 −U 1 )]〉 1

〈M[β(U 1 −U 0 )]〉 0

. (10.46)

This equation is now used as follows. We perform two simulations, one with poten-
tialU 0 , and one with potentialU 1. We consider switching the potential in system
0 fromU 0 toU 1 as a trial move for which we calculate the acceptance ratio, but
which is in fact never carried out. Taking the average of the acceptance ratios of
such moves gives the numerator of(10.46). Similarly, the denominator is given as
the acceptance ratio of trial switches fromU 1 toU 0 , evaluated in the system with
potentialU 1. Again, this method is reliable only for appreciable overlap between
the two equilibrium distributions in phase space. Bennett’s method can also be
extended with a weight function similar to that of Torrie and Valleau, and Ben-
nett calculates the actual form of this function that gives the most accurate results.
Bennett furthermore describes an interpolation method to estimate the free energy
difference even if the overlap between the distributions is very small. His paper
should be consulted for details[33].

10.5.2 Chemical potential determination

Determining the chemical potential is done using relation(7.13), which enables
us to extract this thermodynamic quantity as a free energy difference between two
systems withNandN+1 particles respectively. The exponential of this free energy
difference can be written as the fraction of two partition functions:

ZN+ 1

ZN

=

N!^3 N

(N+ 1 )!^3 (N+^1 )

∫

dRN+ 1 exp[−βU(RN+ 1 )]

∫

dRNexp[−βU(RN)]

=

V

(N+ 1 )^3

〈

1

V

∫

d^3 rN+ 1 exp(−βU+)

〉

N

=e−βμ (10.47)

whereU+is the energy of a particle inserted atrN+ 1 into theN-particle system.
The prefactorV−^3 /(N+ 1 )of the expectation value on the right hand side is
exp(−βμideal), whereμidealis the chemical potential of the ideal gas. The term
within angular brackets gives the expectation value of the Boltzmann factor asso-
ciated with the energy difference for the addition of a particle anywhere in the
system. This expectation value is determined via trial additions which are regularly
performed but never accepted. After each MC step we generate a positionrat
random in the system and calculate the factor exp[−βU(rN+ 1 )]. The factor 1/V
within the brackets ensures the proper evaluation of the average. These values are
then used to calculate the required expectation value by dividing by the number
of such trial additions. This method is called ‘Widom’s particle insertion method’