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=1,\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
Figure 2.1:
Statistical-mechanical Monte Carlo for a model fluid
with continuous potential
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:
$
\begin{eqnarray}
x_{j}' &=& x_{j}+d(\xi -0.5)
\end{eqnarray} $
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:
If $E' > E(k)$, compare the quotient of the two thermodynamic probabilities,
$
\begin{eqnarray}
q &\equiv & e^{\textstyle -[E'-E(k)]/kT}
\end{eqnarray} $
to a random number $\xi \in (0,1)$:
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.
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:
Figure 2.2:
Monte Carlo for hard spheres
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}'\;$.
Yet another modification is needed for spin systems:
Figure 2.3:
Monte Carlo simulation on an Ising lattice
Let $\Gamma_{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_{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)$.
PROJECT MC (FLUID):
Write a subroutine MCSTEP which performs the basic Monte Carlo step
as described in Fig. 2.1: 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):
Let $N=n.n$ spins $\sigma_{i}= \pm 1; \; i=1, \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).