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}|)$.
  1. 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.

  2. 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 \, '$.

  3. 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".


    - 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

< >