Metropolis Monte Carlo Method

The goal of this section is to introduce the Metropolis Monte Carlo Method (J. Chem. Phys. 21,1087, 1953) and to illustrate the algorithm as applied to the computation of canonical ensemble averages for the Ising model.

The Metropolis Monte Carlo method is a computational approach (i.e., an algorithm) for generating a set of N configurations of the system $\xi_1, \xi_2, \xi_3, ..., \xi_N$ such that

$$\lim_{N \rightarrow \infty} \frac{N_{\xi}}{N} = P(\xi),$$ (332)

where P($\xi$) is a given probability distribution (e.g., the Boltzmann distribution P($\xi$) = Z$^{-1}$ exp[-$\beta$ E($\xi$)]) and $N_{\xi}$ is the number of configurations $\xi$ (e.g., the number of configurations generated with a particular arrangement of spins $S_1(\xi),S_2(\xi),...,S_N(\xi)$ in the Ising model).
 

The Metropolis Monte Carlo algorithm can be described as follows:

Step (1): Pick a configuration $\xi_n$ (the initial configuration can be any configuration of the system, e.g., any arrangement of spins in the Ising model).
Step (2): Pick a trial configuration$\xi_t$ (usually a configuration similar to $\xi_n$) and compute the probability ratio $R=\frac{P(\xi_t)}{P(\xi_n)}$. Pick a random number $p$ with value between 0 and 1. Make $\xi_{n+1}=\xi_t$ if $p \leq R$. Otherwise, make $\xi_{n+1}=\xi_n$.
Step (3): Go to (2) replacing $\xi_n$ by $\xi_{n+1}$.
Step (3) is repeated N times, where N is a sufficiently large number. Note that, according to step (2), the probability of accepting a trial configuration $\xi_t$ by making $\xi_{n+1}=\xi_t$ from a configuration $\xi_n$ is

$$P_{\xi_n, \xi_t} = \begin{cases}R=\frac{P(\xi_t)}{P(\xi_n)}, \hsp.......2cm} P(\xi_t) < P(\xi_n), \\ 1, \hspace{.2cm} \text{otherwise}. \end{cases}$$ (333)

The goal of the remaining of this section is to prove that such an algorithm indeed produces an ensemble of configurations that satisfies Eq. (332).

Consider an ensemble of N configurations with N($\xi$) members of the ensemble in state $\xi$. Apply the Metropolis Monte Carlo algorithm to each member of the ensemble by setting $\xi_n=\xi$ and $\xi_t=\xi'$ in step (2), where $\xi$ and $\xi'$ are any two possible states. Note that by applying the algorithm the we generate more configurations and we therefore evolve the initial distribution. In order to show that the algorithm produces an ensemble of configurations that satisfies Eq. (332) we need to show that the any initial distribution $N(\xi)/N$ evolves towards the distribution $P(\xi)=$ and once such a distribution is reached it remains at equilibrium.

According to step (2), for any pair of states $\xi$ and $\xi'$, the number of configurations generated in state $\xi'$ by applying the algorithm to the N($\xi$) configurations in state $\xi$ is $N(\xi) P_{\xi, \xi'}$, where $P_{\xi, \xi'}$ is the probability of accepting the trial configuration $\xi'$ when $\xi_n=\xi$. In addition, the number of configurations generated in state $\xi'$ by applying the algorithm to the N($\xi'$) configurations in state $\xi'$ is (1-P$_{\xi',\xi}$) N($\xi'$). Therefore, the total number $\overline{N}(\xi')$ of configurations generated in state $\xi'$ due to any other state $\xi$ is

$$\overline{N}(\xi') = N(\xi') + \Delta N(\xi'),$$ (334)

where

$$\Delta N(\xi') = N(\xi) P_{\xi, \xi'} - N(\xi') P_{\xi',\xi},$$ (335)

is the net change in the number of configurations in state $\xi'$, relative to $N(\xi')$. According to Eqs. (333) and (335),

$$\Delta N(\xi') = N(\xi) - N(\xi') \frac{P(\xi)}{P(\xi')},$$ (336)

when $P(\xi') > P(\xi)$ and

$$\Delta N(\xi') = N(\xi) \frac{P(\xi')}{P(\xi)} - N(\xi'),$$ (337)

when $P(\xi') < P(\xi)$. Therefore, according to Eqs. (336) and (337), $\Delta N(\xi')=0$ when $N(\xi)/N=P(\xi)$ and $N(\xi')/N=P(\xi')$, i.e., the algorithm does not alter the relative population of the states when the ensemble distribution is equal to the equilibrium distribution. In addition, Eqs. (336) and (337) indicate that $\Delta N(\xi')>0$ when $N(\xi')/N<P(\xi')$ (and $\Delta N(\xi')<0$ when $N(\xi')/N>P(\xi')$), i.e., the algorithm evolves any arbitrary distribution towards the equilibrium distribution where $\frac{N_{\xi}}{N}=P(\xi)$.

Note: The most important aspect of the method is that the algorithm is able generate an ensemble of configurations with the probability distribution
P($\xi$) = Z$^{-1}$ exp[-$\beta$ E($\xi$)], simply by computing the probability ratios$P(\xi')/P(\xi)$. Therefore, the method avoids the need of computing the canonical partition function of the system$Z$, a computational task that would be computationally intractable for most real applications. This numerical technique is thus extremely useful since it allows one to compute any canonical ensembles without having to compute the canonical partition function of the system as follows,

$$\langle A \rangle \approx \bar{A} = \frac{1}{N} \sum_{\xi} N_{\xi} A(\xi),$$ (338)

where $A(\xi)$ is the value of the observable $A$ for state $\xi$ and $\bar{A}$ is the Monte Carlo estimator of $\langle A \rangle$ associated with the finite number of configurations N.
 

Exercise (due 04/24/03):

Implement the Metropolis Monte Carlo Algorithm to generate an ensemble of configurations for a 2-dimensional Ising model with (20 $\times$ 20 spins) in the absence of an external field. Compute the average value of the magnetization at various different temperatures and show that the system exhibits spontaneous magnetization when $T < 2.3 J/k$, where J is the coupling constant between spins.