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.
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
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).