The Model and Method

Here we will consider the simple case of a 2 dimensional square lattice with interactions only between nearest neighbours. In this case

$$E = -J\sum_{i,j_i}\bi{S}_i\cdot\bi{S}_{j_i}$$ (4.28)

where $j_i$ is only summed over the 4 nearest neighbours of $i$.

This model can be studied using the Metropolis method as described in the notes, where the state can be changed by flipping a single spin. Note that the change in energy due to flipping the $k$th spin from $\downarrow$ to $\uparrow$ is given by

$$\Delta E_k = -J\sum_{j_k}\bi{S}_{j_k} .$$ (4.29)

The only quantity which actually occurs in the calculation is

$$Z_k = \exp\left(-\Delta E_k/k_{\rm B} T\right)\quad,$$ (4.30)

and this can only take one of five different values given by the number of neighbouring $\uparrow$ spins. Hence it is sensible to store these in a short array before starting the calculation. Note also that there is really only 1 parameter in the model, $J/k_{\rm B} T$, so that it would make sense to write your program in terms of this single parameter rather than $J$ and $T$ separately.

The calculation should use periodic boundary conditions, in order to avoid spurious effects due to boundaries. There are several different ways to achieve this. One of the most efficient is to think of the system as a single line of spins wrapped round a torus. This way it is possible to avoid a lot of checking for the boundary. For an $N\times N$ system of spins define an array of $2 N^2$ elements using the shortest sensible variable type: char in C(++). It is easier to use $1$ for spin $\uparrow$ and $0$ for spin $\downarrow$, as this makes the calculation of the number of neighbouring $\uparrow$ spins easier. In order to map between spins in a 2d space $S_{\bi{r}}$ and in the 1d array $S_k$ the following mapping can be used.

$$\bi{S}_{\bi{r} + \delta\bi{x}} \mapsto \bi{S}_{k+1},\quad \b... ...,\quad \bi{S}_{\bi{r} - \delta\bi{y}} \mapsto \bi{S}_{k+N^2-N}$$ (4.31)

where the 2nd $N^2$ elements of the array are always maintained equal to the 1st $N^2$. This way it is never necessary to check whether one of the neighbours is over the edge. It is important to remember to change $\bi{S}_{k + N^2}$ whenever $\bi{S}_k$ is changed.

The calculation proceeds as follows:

1. Initialise the spins, either randomly or aligned.
2. Choose a spin to flip. It is better to choose a spin at random rather than systematically as systematic choices can lead to spurious temperature gradients across the system.
3. Decide whether to flip the spin by using the Metropolis condition (see notes).
4. If the spin is to be flipped, do so but remember to flip its mirror in the array.
5. Update the energy and magnetisation.
6. Add the contributions to the required averages.
7. Return to step 2 and repeat.