Franz J. Vesely > CompPhys Tutorial > Selected Applications > Molecular Simulation

 < >
 Part III: Ch. 6

## 6.2 Monte Carlo Method

In a canonical average such as

$\langle a \rangle_{NVT} = \frac{\textstyle 1}{\textstyle Q} \int a(1.. N) \, \exp\{-E(1..N)/kT\} d1..dN$

the denominator $Q \equiv \int \exp\{-E(1..N)/kT\} d1..d N$ is usually unknown.

But in the section on "Biased Random Walks" we learned that this is no hindrance for the calculation of averages:

Writing, for a certain $N$-particle configuration, $\Gamma_{c} \equiv \left\{1..N \right\}$ we generate a Markov chain of, say, $K$ configurations $\{ \Gamma_{c}(k), \, k \dots K\}$ such that the relative frequency of a configuration in the chain is proportional to the corresponding Boltzmann factor. An estimate for the mean value $\langle a \rangle$ is then

$\langle a \rangle = \frac{\textstyle 1}{\textstyle K} \sum \limits_{ k}^{ K} a \left[\Gamma_{c}(k) \right]$

Here is the procedure due to N. METROPOLIS:

 Metropolis Monte Carlo: Let $\Gamma_{c}(k) \equiv \left\{r_{1} \dots r_{N} \right\}$ be given; the potential energy of this configuration is $E(k) \equiv (1/2) \sum_{i} \sum_{j} u(|r_{j}-r_{i}|)$. Generate a "neighboring" configuration $\Gamma_{c}'$ by randomly moving one of the $N$ particles within a cubic region centered around its present position: $x_{j}' = x_{j}+d(\xi -0.5)$ and similarly for $y_{j},z_{j}$. Here, $d$ (= side length of the displacement cube) is a parameter to be optimized (see below), and $\xi$ is a random number from an equidistribution in $(0,1)$. The number $j$ of the particle to be moved may either be drawn among the $N$ candidates, or may run cyclically through the set of particle indices. Determine the modified total energy $E\,'$; since displacing particle $j$ affects only $N-1$ of the $N(N-1)/2$ pair distances in the system, it is not necessary to recalculate the entire double sum to get $E \, '$. If $E \, ' \leq E(k)$, we accept $\Gamma_{c}'$ as the next element of the Markov chain: $E\,' \leq E(k): \;\; \Rightarrow \; \Gamma_{c}(k+1)= \Gamma_{c}' ;\;\;{\rm go\;to\; (1)}$ If $E\,' > E(k)$, compare the quotient of the two thermodynamic probabilities, $q \equiv e^{\textstyle -[E\,'-E(k)]/kT}$ to a random number $\xi \in (0,1)$: $\begin{eqnarray} E\,' > E(k): \;\; && \\ \xi & \leq & q:\;\;\; \Rightarrow \; \Gamma_{c}(k+1)= \Gamma_{c}' ;\;\;{\rm go \; to\; (1)} \\ \xi & > & q: \;\;\; \Rightarrow \; \Gamma_{c}(k+1)= \Gamma_{c}(k) ;\;\; {\rm go \; to\; (1)} \end{eqnarray}$ This is the so-called "asymmetric rule". Remarks: - The parameter $d$ should be adjusted empirically in such a way that in step 3 approximately one out of two attempted steps $\Gamma_{c}'$ is accepted. - The random variate sampled in step 1 need not come from an equidistribution; any probability density that is symmetrical about $0$, such as a Gauss distribution, will do.
 Figure 6.2: Statistical-mechanical Monte Carlo for a model fluid with continuous potential

In the case of hard disks or spheres the 3rd step in the above recipe must be modified. Values of $E(k)$ and $E\,'$ are then restricted to $0$ or $\infty$, with Boltzmann factors $1$ or $0,$ respectively. Here is the modified part of the MC procedure:

 Let $\Gamma_{c}(k) \equiv \left\{ r_{1} \dots r_{N} \right\}$ be given. Trial move $\Gamma_{c}(k)$ $\longrightarrow \Gamma_{c}'$: $x_{j}' = x_{j}+d(\xi-0.5) \;\;\;\;\mbox{etc., for} \;\;y_{j}, z_{j}$ If particle $j$ now overlaps with any other particle, let $\Gamma_{c}(k+1)= \Gamma_{c}(k)$; otherwise let $\Gamma_{c}(k+1)= \Gamma_{c}'$.
 Figure 6.3: Monte Carlo for hard spheres

Yet another modification is needed for spin systems:

 Let $\Gamma_{\textstyle c}(k)\equiv \left\{ \sigma_{1},\dots,\sigma_{N}\right\}\;$ be given. Pick a spin $\sigma_{i}$ and tentatively invert it. The resulting energy change is $\Delta E=A \sigma_{i} \sum \limits_{j(i)}^{4} \sigma_{j}$ If $\Delta E \leq 0$, accept the inverted spin: $\sigma_{i}(k+1) = - \sigma_{i}(k)$; otherwise, draw an equidistributed $\xi \in (0,1)$ and compare it to $w \equiv \mbox{exp}[-\Delta E/kT]$; if $\xi < w$, accept $-\sigma_{i}$, else leave $\sigma_{i}$ unchanged: $\sigma_{i}(k+1) = \sigma_{i}(k)$.
 Figure 6.4: Monte Carlo simulation on an Ising lattice

PROJECT MC (FLUID): (See also here) Write a subroutine MCSTEP which performs the basic Monte Carlo step as described in Fig. 6.2: selecting at random one of the LJ particles that were placed on a lattice by STARTCONFIG, displace it slightly and apply the PBC; then compute the new potential energy (using NIC!) and check whether the modified configuration is accepted or not, given a specific temperature $T^{*}$; if accepted, the next configuration is the modified one, otherwise the old configuration is retained for another step.

Write a main routine to combine the subroutines STARTCONF, ENERGY, and MCSTEP into a working MC program.

PROJECT MC (LATTICE): (See also here) Let $N=n.n$ spins $\sigma_{i}= \pm 1; \; i, \dots N$ be situated on the vertices of a two-dimensional square lattice. The interaction energy is defined by

$E = \sum \limits_{i} E_{i} = -\frac{1}{2} \sum \limits_{i}^{N} \sum \limits_{j}^{4} \sigma_{i} \sigma_{j}$

where the sum over $j$ involves the 4 nearest neighbors of spin $i$. Periodic boundary conditions are assumed
• Write a Monte Carlo program to perform a biased random walk through configuration space.
• Determine the mean total moment $\langle M \rangle \equiv \langle \sum_{i} \sigma_{i} \rangle$ and its variance as a function of the quantity $1/kT$. Compare your results to literature data (e.g. BINDER, K.: Applications of the Monte Carlo Method in Statistical Physics. Springer, Berlin 1987).

vesely 2006

 < >