## Liquid-State Physical Chemistry: Fundamentals, Modeling, and Applications (2013)

### 9. Modeling the Structure of Liquids: The Simulation Approach

### 9.3. The Monte Carlo Method

Whereas, in MD the structure and dynamics are calculated, in MC calculations attention is focused on the structure alone.

The method consists essentially of distributing molecules in space, calculating the associated energy, modifying the molecular distribution, and then repeating the calculation so as to obtain eventually a sequence of, say *n*, molecular distributions in configuration space that samples the desired ensemble (*NVE*, *NVT*, or otherwise) properly. If *f* represents any observable dependent only on the configuration, then each element of the sequence will generate a sample *f _{j}* of the observable and the mean

*f*

_{m}=

*n*

^{−}^{1}Σ

*is an estimator of the ensemble average ⟨*

_{j}f_{j}*f*⟩. The essence of MC simulations is to construct the sequence in a smart fashion. The procedure is best explained by using ideas from calculating multidimensional integrals in numerical mathematics.

Suppose one has the integral

(9.16)

To evaluate this integral, the following trick is used. One writes

(9.17)

where *f*(*x*) > 0 over the interval [*a*,*b*] and *f*(*x*) is integrable. We treat *f*(*x*) as a probability density function so that . Then we can write

(9.18)

where <*h*(*x*)> is the average value over the interval [*a*,*b*] assuming that *x* is distributed as *f*(*x*). Hence, it follows that if we choose points in the interval [*a*,*b*] at random__ ^{5)}__ with probability

*f*(

*x*) we can evaluate <

*h*(

*x*)> and by multiplying with the normalization factor

*F*, calculate

*G*.

In statistical thermodynamics we may wish to evaluate integrals such as

__(9.19)__

where the configurational partition function *Q* acts as a normalizing factor. Consequently, when we choose position vectors according to the distribution function

(9.20)

we can evaluate integrals of the type of __Eq. (9.19)__. Unfortunately, *Q*, when acting as normalization factor, generally cannot be evaluated, but it is still possible to determine thermodynamic properties and structure.

The most frequently used method is that of Metropolis *et al*. [10], based on Markov processes. In this method one starts with an initial configuration of molecules and changes the configuration over a long sequence. The transition probabilities *T _{ij}* of the transition from state

^{6)}*i*to state

*j*are taken in such a way that in the infinite limit they tend to be distributed according to exp[−

*Φ*(

**r**

_{1},…,

**r**

*)/*

_{N}*kT*]. Note that the total probability to end up in any state starting in state

*i*is Σ

*, so that we should have Σ*

_{j}T_{ij}*= 1. If is the probability to go from state*

_{j}T_{ij}*i*to state

*j*in

*m*steps, we have . It can be shown for Markov processes that, provided all states are accessible from any other state in a finite number of steps, the limit of for

*m*→ ∞ becomes

*P*, independent of

_{j}*i*. This implies that, after a sufficient number of steps, the state of the system is independent of the initial configuration. In the limit we have

__(9.21)__

stating that the probabilities are positive, that they are normalized, and that each state eventually is reached from any other state. We choose

__(9.22)__

From __Eq. (9.22)__ it follows that we have the detailed balance condition

__(9.23)__

If we substitute this result in __Eq. (9.21)__–__(9.23)__, we obtain the *P _{i}* = Σ

*= Σ*

_{j}T_{ji}P_{j}*=*

_{j}T_{ij}P_{i}*P*Σ

_{i}*or Σ*

_{j}T_{ij}*= 1, satisfying our earlier requirement. Hence, by choosing the relationship between*

_{j}T_{ij}*T*and

_{ij}*T*according to

_{ji}__Eq. (9.22)__, all requirements for a Markov process are fulfilled. What remains to be done is to choose the appropriate values for

*T*.

_{ij}Metropolis *et al*. used for *T _{ij}* the following recipe

(9.24)

The constants *A _{ij}* are chosen as

*A*=

_{ij}*A*, and determined as follows:

_{ji}· A molecule is chosen at random.

· A change in position (δ*x*, δ*y*, δ*z*) is selected.

· *A _{ij}* = 0 if any of |δ

*x*|, |δ

*y*| and |δ

*z*| > Δ, and

*A*= (2Δ)

_{ij}^{−3}otherwise.

With these rules a large number of configurations is generated with the constant Δ chosen in such a way that about half of the configurations are accepted.

Calculations are initiated by providing the positions of *N* molecules in the cell chosen; this starting configuration can be selected as indicated in Section 9.1. Normally the first – say 200 000 – configurations are discarded in order to limit the influence of the starting configuration, and thereafter the properties are evaluated using, say 500 000, configurations. The above-outlined procedure is for the canonical ensemble, and enables one to calculate the ensemble average by direct sampling of phase space:

(9.25)

where *G _{j}*(

**r**

_{1},

*…*,

**r**

*) represents the value of*

_{N}*G*in configuration

*j*. For example, the internal energy

*U*is given by

(9.26)

where *U*_{kin} is the kinetic energy, given for spherical molecules by *U*_{kin} = 3*NkT*/2 and is due to the term *Λ*^{−3N}. This can also be written as

(9.27)

so that the internal energy can be estimated by averaging the potential energy of all configurations in the MC sequence. The heat capacity is then

(9.28)

with *C _{V}*

_{,kin}= 3

*Nk*/2 the contribution of the kinetic energy to

*C*. The pressure can be calculated from the virial theorem expression (see

_{V}__Chapter 3__) reading

(9.29)

with ∇* _{i}* the gradient operator for molecule

*i*. Finally, using the number of particles

*n*(

*r*) with distance between

*r*and

*r*+ Δ

*r*from a reference particle and taking subsequently each particle as reference, from the expression

(9.30)

one can obtain the correlation function *g*(*r*).

Today, the MC method is used much like the MD method, more or less routinely, to calculate the properties of thermodynamic systems. Both, the MC and MD methods have the property of providing, in principle, exact calculations for a given interaction model, and the typical accuracy is about 1–2%. Provided that the same potentials are used and also the same conditions, the equilibrium results of the MD method agree to a large degree with those of the MC method (see __Figure 9.3__). Whereas, the MD method yields the full dynamics, it is somewhat more complex to eliminate time correlations when compared to the MC method.

Problem 9.3*

Is it possible to calculate *Q _{N}* in principle from a MC simulation?