Computational Physics Lectures: Statistical physics and the Ising Model

Ensembles

In statistical physics the concept of an ensemble is one of the cornerstones in the definition of thermodynamical quantities. An ensemble is a collection of microphysics systems from which we derive expectations values and thermodynamical properties related to experiment. As an example, the specific heat (which is a measurable quantity in the laboratory) of a system of infinitely many particles, can be derived from the basic interactions between the microscopic constituents. The latter can span from electrons to atoms and molecules or a system of classical spins. All these microscopic constituents interact via a well-defined interaction. We say therefore that statistical physics bridges the gap between the microscopic world and the macroscopic world. Thermodynamical quantities such as the specific heat or net magnetization of a system can all be derived from a microscopic theory.

Famous Ensembles

The table lists the most used ensembles in statistical physics together with frequently arising extensive (depend on the size of the systems such as the number of particles) and intensive variables (apply to all components of a system), in addition to associated potentials.

Canonical Ensemble

One of the most used ensembles is the canonical one, which is related to the microcanonical ensemble via a Legendre transformation. The temperature is an intensive variable in this ensemble whereas the energy follows as an expectation value. In order to calculate expectation values such as the mean energy $\langle E \rangle$ at a given temperature, we need a probability distribution. It is given by the Boltzmann distribution

 
$$\begin{equation*} P_i(\beta) = \frac{e^{-\beta E_i}}{Z} \end{equation*}$$

 
with $\beta=1/k_BT$ being the inverse temperature, $k_B$ is the Boltzmann constant, $E_i$ is the energy of a microstate $i$ while $Z$ is the partition function for the canonical ensemble defined as

The partition function is a normalization constant

In the canonical ensemble the partition function is

 
$$\begin{equation*} Z=\sum_{i=1}^{M}e^{-\beta E_i}, \end{equation*}$$

 
where the sum extends over all microstates $M$.

Helmoltz free energy, what does it mean?

The potential of interest in this case is Helmholtz' free energy. It relates the expectation value of the energy at a given temperatur $T$ to the entropy at the same temperature via

 
$$\begin{equation*} F=-k_{B}TlnZ=\langle E \rangle-TS. \end{equation*}$$

Helmholtz' free energy expresses the struggle between two important principles in physics, namely the strive towards an energy minimum and the drive towards higher entropy as the temperature increases. A higher entropy may be interpreted as a larger degree of disorder. When equilibrium is reached at a given temperature, we have a balance between these two principles. The numerical expression is Helmholtz' free energy.

Thermodynamical quantities

In the canonical ensemble the entropy is given by

 
$$\begin{equation*} S =k_{B}lnZ +k_{B}T\left(\frac{\partial lnZ}{\partial T}\right)_{N, V}, \end{equation*}$$

 
and the pressure by

 
$$\begin{equation*} p=k_{B}T\left(\frac{\partial lnZ}{\partial V}\right)_{N, T}. \end{equation*}$$

 
Similarly we can compute the chemical potential as

 
$$\begin{equation*} \mu =-k_{B}T\left(\frac{\partial lnZ}{\partial N}\right)_{V, T}. \end{equation*}$$

Thermodynamical quantities, the energy in the canonical ensemble

For a system described by the canonical ensemble, the energy is an expectation value since we allow energy to be exchanged with the surroundings (a heat bath with temperature $T$).

This expectation value, the mean energy, can be calculated using

 
$$\begin{equation*} \langle E\rangle =k_{B}T^{2}\left(\frac{\partial lnZ}{\partial T}\right)_{V, N} \end{equation*}$$

 
or using the probability distribution $P_i$ as

 
$$\begin{equation*} \langle E \rangle = \sum_{i=1}^M E_i P_i(\beta)= \frac{1}{Z}\sum_{i=1}^M E_ie^{-\beta E_i}. \end{equation*}$$

Energy and specific heat in the canonical ensemble

The energy is proportional to the first derivative of the potential, Helmholtz' free energy. The corresponding variance is defined as

 
$$\begin{equation*} \sigma_E^2=\langle E^2 \rangle-\langle E \rangle^2= \frac{1}{Z}\sum_{i=1}^M E_i^2e^{-\beta E_i}- \left(\frac{1}{Z}\sum_{i=1}^M E_ie^{-\beta E_i}\right)^2. \end{equation*}$$

 
If we divide the latter quantity with $kT^2$ we obtain the specific heat at constant volume

 
$$\begin{equation*} C_V= \frac{1}{k_BT^2}\left(\langle E^2 \rangle-\langle E \rangle^2\right), \end{equation*}$$

 
which again can be related to the second derivative of Helmholtz' free energy.

Magnetic moments and susceptibility in the canonical ensemble

Using the same prescription, we can also evaluate the mean magnetization through

 
$$\begin{equation*} \langle {\cal M} \rangle = \sum_i^M {\cal M}_i P_i(\beta)= \frac{1}{Z}\sum_i^M {\cal M}_ie^{-\beta E_i}, \end{equation*}$$

 
and the corresponding variance

 
$$\begin{equation*} \sigma_{{\cal M}}^2=\langle {\cal M}^2 \rangle-\langle {\cal M} \rangle^2= \frac{1}{Z}\sum_{i=1}^M {\cal M}_i^2e^{-\beta E_i}- \left(\frac{1}{Z}\sum_{i=1}^M {\cal M}_ie^{-\beta E_i}\right)^2. \end{equation*}$$

 
This quantity defines also the susceptibility $\chi$

 
$$\begin{equation*} \chi=\frac{1}{k_BT}\left(\langle {\cal M}^2 \rangle-\langle {\cal M} \rangle^2\right). \end{equation*}$$

Our model, the Ising model in one and two dimensions

The model we will employ in our studies of phase transitions at finite temperature for magnetic systems is the so-called Ising model. In its simplest form the energy is expressed as

 
$$\begin{equation*} E=-J\sum_{< kl>}^{N}s_ks_l-{\cal B}\sum_k^Ns_k, \end{equation*}$$

 
with $s_k=\pm 1$, $N$ is the total number of spins, $J$ is a coupling constant expressing the strength of the interaction between neighboring spins and ${\cal B}$ is an external magnetic field interacting with the magnetic moment set up by the spins.

The symbol $< kl>$ indicates that we sum over nearest neighbors only. Notice that for $J>0$ it is energetically favorable for neighboring spins to be aligned. This feature leads to, at low enough temperatures, a cooperative phenomenon called spontaneous magnetization. That is, through interactions between nearest neighbors, a given magnetic moment can influence the alignment of spins that are separated from the given spin by a macroscopic distance. These long range correlations between spins are associated with a long-range order in which the lattice has a net magnetization in the absence of a magnetic field.

Boltzmann distribution

In order to calculate expectation values such as the mean energy $\langle E \rangle$ or magnetization $\langle {\cal M} \rangle$ in statistical physics at a given temperature, we need a probability distribution

 
$$\begin{equation*} P_i(\beta) = \frac{e^{-\beta E_i}}{Z} \end{equation*}$$

 
with $\beta=1/kT$ being the inverse temperature, $k$ the Boltzmann constant, $E_i$ is the energy of a state $i$ while $Z$ is the partition function for the canonical ensemble defined as

 
$$\begin{equation*} Z=\sum_{i=1}^{M}e^{-\beta E_i}, \end{equation*}$$

 
where the sum extends over all microstates $M$. $P_i$ expresses the probability of finding the system in a given configuration $i$.

Energy for a specific configuration

The energy for a specific configuration $i$ is given by

 
$$\begin{equation*} E_i =-J\sum_{< kl>}^{N}s_ks_l. \end{equation*}$$

Configurations

 
$$\begin{equation*} \begin{array}{cccccccccc} \uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow&\uparrow&\dots&\uparrow&\downarrow\\ 1&2&3&\dots& i-1&i&i+1&\dots&N-1&N\end{array} \end{equation*}$$

 
In order to illustrate these features let us further specialize to just two spins.

With two spins, since each spin takes two values only, we have $2^2=4$ possible arrangements of the two spins. These four possibilities are

 
$$\begin{equation*} 1= \uparrow\uparrow\hspace{1cm} 2= \uparrow\downarrow\hspace{1cm} 3= \downarrow\uparrow\hspace{1cm} 4=\downarrow\downarrow \end{equation*}$$

Boundary conditions, free ends

What is the energy of each of these configurations?

For small systems, the way we treat the ends matters. Two cases are often used.

In the first case we employ what is called free ends. This means that there is no contribution from points to the right or left of the endpoints. For the one-dimensional case, the energy is then written as a sum over a single index

 
$$\begin{equation*} E_i =-J\sum_{j=1}^{N-1}s_js_{j+1}, \end{equation*}$$

Free ends and the energy

If we label the first spin as $s_1$ and the second as $s_2$ we obtain the following expression for the energy

 
$$\begin{equation*} E=-Js_1s_2. \end{equation*}$$

 
The calculation of the energy for the one-dimensional lattice with free ends for one specific spin-configuration can easily be implemented in the following lines

for ( j=1; j < N; j++) {
     energy += spin[j]*spin[j+1]; 
}

where the vector $spin[]$ contains the spin value $s_k=\pm 1$.

Free ends and energy

For the specific state $E_1$, we have chosen all spins up. The energy of this configuration becomes then

 
$$\begin{equation*} E_1=E_{\uparrow\uparrow}=-J. \end{equation*}$$

 
The other configurations give

 
$$\begin{equation*} E_2=E_{\uparrow\downarrow}=+J, \end{equation*}$$

 
$$\begin{equation*} E_3=E_{\downarrow\uparrow}=+J, \end{equation*}$$

 
and

 
$$\begin{equation*} E_4=E_{\downarrow\downarrow}=-J. \end{equation*}$$

Periodic boundary conditions

We can also choose so-called periodic boundary conditions. This means that the neighbour to the right of $s_N$ is assumed to take the value of $s_1$. Similarly, the neighbour to the left of $s_1$ takes the value $s_N$. In this case the energy for the one-dimensional lattice reads

 
$$\begin{equation*} E_i =-J\sum_{j=1}^{N}s_js_{j+1}, \end{equation*}$$

 
and we obtain the following expression for the two-spin case

 
$$\begin{equation*} E=-J(s_1s_2+s_2s_1). \end{equation*}$$

Energy with PBC

In this case the energy for $E_1$ is different, we obtain namely

 
$$\begin{equation*} E_1=E_{\uparrow\uparrow}=-2J. \end{equation*}$$

 
The other cases do also differ and we have

 
$$\begin{equation*} E_2=E_{\uparrow\downarrow}=+2J, \end{equation*}$$

 
$$\begin{equation*} E_3=E_{\downarrow\uparrow}=+2J, \end{equation*}$$

 
and

 
$$\begin{equation*} E_4=E_{\downarrow\downarrow}=-2J. \end{equation*}$$

Simple code for PBC

If we choose to use periodic boundary conditions we can code the above expression as

    jm=N;
    for ( j=1; j <=N ; j++) {
        energy += spin[j]*spin[jm]; 
        jm = j ;
    }

The magnetization is however the same, defined as

 
$$\begin{equation*} {\cal M}_i=\sum_{j=1}^N s_j, \end{equation*}$$

 
where we sum over all spins for a given configuration $i$.

Summing up

The table lists the energy and magnetization for both free ends and periodic boundary conditions.

Reorganizing

We can reorganize according to the number of spins pointing up, as shown in the table here

Our model, the Ising model in one and two dimensions

It is worth noting that for small dimensions of the lattice, the energy differs depending on whether we use periodic boundary conditions or free ends. This means also that the partition functions will be different, as discussed below. In the thermodynamic limit we have $N\rightarrow \infty$, and the final results do not depend on the kind of boundary conditions we choose.

For a one-dimensional lattice with periodic boundary conditions, each spin sees two neighbors. For a two-dimensional lattice each spin sees four neighboring spins. How many neighbors does a spin see in three dimensions?

Ising model in one and two dimensions

In a similar way, we could enumerate the number of states for a two-dimensional system consisting of two spins, i.e., a $2\times 2$ Ising model on a square lattice with {\em periodic boundary conditions}. In this case we have a total of $2^4=16$ states. Some examples of configurations with their respective energies are listed here

 
$$\begin{equation*} E=-8J\hspace{1cm}\begin{array}{cc}\uparrow & \uparrow \\ \uparrow & \uparrow\end{array} \hspace{0.5cm} E=0\hspace{1cm}\begin{array}{cc}\uparrow & \uparrow \\ \uparrow & \downarrow\end{array} \hspace{0.5cm} E=0\hspace{1cm}\begin{array}{cc}\downarrow & \downarrow \\ \uparrow & \downarrow\end{array} \hspace{0.5cm} E=-8J\hspace{1cm}\begin{array}{cc}\downarrow & \downarrow \\ \downarrow & \downarrow\end{array} \end{equation*}$$

List of configurations with energies and magnetic moment

In the table here we group these configurations according to their total energy and magnetization.

Phase Transitions and Critical Phenomena

A phase transition is marked by abrupt macroscopic changes as external parameters are changed, such as an increase of temperature. The point where a phase transition takes place is called a critical point.

We distinguish normally between two types of phase transitions; first-order transitions and second-order transitions. An important quantity in studies of phase transitions is the so-called correlation length $\xi$ and various correlations functions like spin-spin correlations. For the Ising model we shall show below that the correlation length is related to the spin-correlation function, which again defines the magnetic susceptibility. The spin-correlation function is nothing but the covariance and expresses the degree of correlation between spins.

Phase Transitions and Critical Phenomena, correlation length

The correlation length defines the length scale at which the overall properties of a material start to differ from its bulk properties. It is the distance over which the fluctuations of the microscopic degrees of freedom (for example the position of atoms) are significantly correlated with each other. Usually it is of the order of few interatomic spacings for a solid. The correlation length $\xi$ depends however on external conditions such as pressure and temperature.

Classification of phase transitions

First order/discontinuous phase transitions are characterized by two or more states on either side of the critical point that can coexist at the critical point. As we pass through the critical point we observe a discontinuous behavior of thermodynamical functions. The correlation length is normally finite at the critical point. Phenomena such as hysteris occur, viz. there is a continuation of state below the critical point into one above the critical point. This continuation is metastable so that the system may take a macroscopically long time to readjust. A classical example is the melting of ice. It takes a specific amount of time before all the ice has melted. The temperature remains constant and water and ice can coexist for a macroscopic time. The energy shows a discontinuity at the critical point, reflecting the fact that a certain amount of heat is needed in order to melt all the ice

Second-order phase Transitions

Second order or continuous transitions are different and in general much difficult to understand and model. The correlation length diverges at the critical point, fluctuations are correlated over all distance scales, which forces the system to be in a unique critical phase. The two phases on either side of the critical point become identical. The disappearance of a spontaneous magnetization is a classical example of a second-order phase transitions. Structural transitions in solids are other types of second-order phase transitions.

Phase Transitions and Critical Phenomena

Eherenfest definition of phase Transitions

Using Ehrenfest's definition of the order of a phase transition we can relate the behavior around the critical point to various derivatives of the thermodynamical potential. In the canonical ensemble we are using, the thermodynamical potential is Helmholtz' free energy

 
$$\begin{equation*} F= \langle E\rangle -TS = -kTln Z \end{equation*}$$

 
meaning $ lnZ = -F/kT = -F\beta$. The energy is given as the first derivative of $F$

 
$$\begin{equation*} \langle E \rangle=-\frac{\partial lnZ}{\partial \beta} =\frac{\partial (\beta F)}{\partial \beta}. \end{equation*}$$

 
and the specific heat is defined via the second derivative of $F$

 
$$\begin{equation*} C_V=-\frac{1}{kT^2}\frac{\partial^2 (\beta F)}{\partial\beta^2}. \end{equation*}$$

Phase Transitions and Critical Phenomena

We can relate observables to various derivatives of the partition function and the free energy. When a given derivative of the free energy or the partition function is discontinuous or diverges (logarithmic divergence for the heat capacity from the Ising model) we talk of a phase transition of order of the derivative. A first-order phase transition is recognized in a discontinuity of the energy, or the first derivative of $F$. The Ising model exhibits a second-order phase transition since the heat capacity diverges. The susceptibility is given by the second derivative of $F$ with respect to external magnetic field. Both these quantities diverge.

The Ising Model and Phase Transitions

The Ising model in two dimensions with ${\cal B} = 0$ undergoes a phase transition of second order. What it actually means is that below a given critical temperature $T_C$, the Ising model exhibits a spontaneous magnetization with $\langle {\cal M} \rangle\ne 0$. Above $T_C$ the average magnetization is zero. The mean magnetization approaches zero at $T_C$ with an infinite slope. Such a behavior is an example of what are called critical phenomena. A critical phenomenon is normally marked by one or more thermodynamical variables which vanish above a critical point. In our case this is the mean magnetization $\langle {\cal M} \rangle\ne 0$. Such a parameter is normally called the order parameter.

The Ising Model and Phase Transitions, mean magnetization

It is possible to show that the mean magnetization is given by (for temperature below $T_C$)

 
$$\begin{equation*} \langle {\cal M}(T) \rangle \sim \left(T-T_C\right)^{\beta}, \end{equation*}$$

 
where $\beta=1/8$ is a so-called critical exponent. A similar relation applies to the heat capacity

 
$$\begin{equation*} C_V(T) \sim \left|T_C-T\right|^{-\alpha}, \end{equation*}$$

 
and the susceptibility

 
$$\begin{equation*} \chi(T) \sim \left|T_C-T\right|^{-\gamma}, \end{equation*}$$

 
with $\alpha = 0$ and $\gamma = -7/4$.

The Ising Model and Phase Transitions, correlation length

Another important quantity is the correlation length, which is expected to be of the order of the lattice spacing for $T$ is close to $T_C$. Because the spins become more and more correlated as $T$ approaches $T_C$, the correlation length increases as we get closer to the critical temperature. The discontinuous behavior of the correlation $\xi$ near $T_C$ is

 
$$\begin{equation} \xi(T) \sim \left|T_C-T\right|^{-\nu}. \tag{1} \end{equation}$$

The Ising Model and Phase Transitions, critical behavior

A second-order phase transition is characterized by a correlation length which spans the whole system. The correlation length is typically of the order of some few interatomic distances. The fact that a system like the Ising model, whose energy is described by the interaction between neighboring spins only, can yield correlation lengths of macroscopic size at a critical point is still a feature which is not properly understood.

The Ising Model and Phase Transitions, critical temperature

In our actual calculations of the two-dimensional Ising model, we are however always limited to a finite lattice and $\xi$ will be proportional with the size of the lattice at the critical point. Through finite size scaling relations it is possible to relate the behavior at finite lattices with the results for an infinitely large lattice. The critical temperature scales then as

 
$$\begin{equation} T_C(L)-T_C(L=\infty) \propto aL^{-1/\nu}, \tag{2} \end{equation}$$

 
with $a$ a constant and $\nu$ defined in Eq. (1).

The Ising Model and Phase Transitions, correlation length

The correlation length for a finite lattice size can then be shown to be proportional to

 
$$\begin{equation*} \xi(T) \propto L\sim \left|T_C-T\right|^{-\nu}. \end{equation*}$$

 
and if we set $T=T_C$ one can obtain the following relations for the magnetization, energy and susceptibility for $T \le T_C$

 
$$\begin{equation*} \langle {\cal M}(T) \rangle \sim \left(T-T_C\right)^{\beta} \propto L^{-\beta/\nu}, \end{equation*}$$

 
$$\begin{equation*} C_V(T) \sim \left|T_C-T\right|^{-\gamma} \propto L^{\alpha/\nu}, \end{equation*}$$

 
and

 
$$\begin{equation*} \chi(T) \sim \left|T_C-T\right|^{-\alpha} \propto L^{\gamma/\nu}. \end{equation*}$$

The Metropolis Algorithm and the Two-dimensional Ising Model

In our case we have as the Monte Carlo sampling function the probability for finding the system in a state $s$ given by

 
$$\begin{equation*} P_s=\frac{e^{-(\beta E_s)}}{Z}, \end{equation*}$$

 
with energy $E_s$, $\beta=1/kT$ and $Z$ is a normalization constant which defines the partition function in the canonical ensemble. As discussed above

 
$$\begin{equation*} Z(\beta)=\sum_se^{-(\beta E_s)} \end{equation*}$$

 
is difficult to compute since we need all states.

The Metropolis Algorithm and the Two-dimensional Ising Model

In a calculation of the Ising model in two dimensions, the number of configurations is given by $2^N$ with $N=L\times L$ the number of spins for a lattice of length $L$. Fortunately, the Metropolis algorithm considers only ratios between probabilities and we do not need to compute the partition function at all. The algorithm goes as follows

1. Establish an initial state with energy $E_b$ by positioning yourself at a random configuration in the lattice
2. Change the initial configuration by flipping e.g., one spin only. Compute the energy of this trial state $E_t$.
3. Calculate $\Delta E=E_t-E_b$. The number of values $\Delta E$ is limited to five for the Ising model in two dimensions, see the discussion below.
4. If $\Delta E \le 0$ we accept the new configuration, meaning that the energy is lowered and we are hopefully moving towards the energy minimum at a given temperature. Go to step 7.
5. If $\Delta E > 0$, calculate $w=e^{-(\beta \Delta E)}$.
6. Compare $w$ with a random number $r$. If $r \le w$, then accept the new configuration, else we keep the old configuration.
7. The next step is to update various expectations values.
8. The steps (2)-(7) are then repeated in order to obtain a sufficently good representation of states.
9. Each time you sweep through the lattice, i.e., when you have summed over all spins, constitutes what is called a Monte Carlo cycle. You could think of one such cycle as a measurement. At the end, you should divide the various expectation values with the total number of cycles. You can choose whether you wish to divide by the number of spins or not. If you divide with the number of spins as well, your result for e.g., the energy is now the energy per spin.

The Metropolis Algorithm and the Two-dimensional Ising Model, practical issues

The crucial step is the calculation of the energy difference and the change in magnetization. This part needs to be coded in an as efficient as possible way since the change in energy is computed many times. In the calculation of the energy difference from one spin configuration to the other, we will limit the change to the flipping of one spin only. For the Ising model in two dimensions it means that there will only be a limited set of values for $\Delta E$. Actually, there are only five possible values.

Five possible energy differences

To see this, select first a random spin position $x,y$ and assume that this spin and its nearest neighbors are all pointing up. The energy for this configuration is $E=-4J$. Now we flip this spin as shown below. The energy of the new configuration is $E=4J$, yielding $\Delta E=8J$.

 
$$\begin{equation*} E=-4J\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \uparrow & \uparrow & \uparrow\\ & \uparrow & \end{array} \hspace{1cm}\Longrightarrow\hspace{1cm} E=4J\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \uparrow & \downarrow & \uparrow\\ & \uparrow & \end{array} \end{equation*}$$

 
The four other possibilities are as follows

 
$$\begin{equation*} E=-2J\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \downarrow & \uparrow & \uparrow\\ & \uparrow & \end{array} \hspace{1cm}\Longrightarrow\hspace{1cm} E=2J\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \downarrow & \downarrow & \uparrow\\ & \uparrow & \end{array} \end{equation*}$$

 
with $\Delta E=4J$,

 
$$\begin{equation*} E=0\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \downarrow & \uparrow & \uparrow\\ & \downarrow & \end{array} \hspace{1cm}\Longrightarrow\hspace{1cm} E=0\hspace{1cm}\begin{array}{ccc} & \uparrow & \\ \downarrow & \downarrow & \uparrow\\ & \downarrow & \end{array} \end{equation*}$$

 
with $\Delta E=0$,

 
$$\begin{equation*} E=2J\hspace{1cm}\begin{array}{ccc} & \downarrow & \\ \downarrow & \uparrow & \uparrow\\ & \downarrow & \end{array} \hspace{1cm}\Longrightarrow\hspace{1cm} E=-2J\hspace{1cm}\begin{array}{ccc} & \downarrow & \\ \downarrow & \downarrow & \uparrow\\ & \downarrow & \end{array} \end{equation*}$$

 
with $\Delta E=-4J$ and finally

 
$$\begin{equation*} E=4J\hspace{1cm}\begin{array}{ccc} & \downarrow & \\ \downarrow & \uparrow & \downarrow\\ & \downarrow & \end{array} \hspace{1cm}\Longrightarrow\hspace{1cm} E=-4J\hspace{1cm}\begin{array}{ccc} & \downarrow & \\ \downarrow & \downarrow & \downarrow\\ & \downarrow & \end{array} \end{equation*}$$

 
with $\Delta E=-8J$.

The Metropolis Algorithm and the Two-dimensional Ising Model, elements of program

This means in turn that we could construct an array which contains all values of $e^{\beta \Delta E}$ before doing the Metropolis sampling. Else, we would have to evaluate the exponential at each Monte Carlo sampling. For the two-dimensional Ising model there are only five possible values. It is rather easy to convice oneself that for the one-dimensional Ising model we have only three possible values. The main part of the Ising model program is shown here

/* 
   Program to solve the two-dimensional Ising model 
   The coupling constant J = 1
   Boltzmann's constant = 1, temperature has thus dimension energy
   Metropolis sampling is used. Periodic boundary conditions.
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include "lib.h"
using namespace  std;
ofstream ofile;
// inline function for periodic boundary conditions
inline int periodic(int i, int limit, int add) { 
  return (i+limit+add) % (limit);
}
// Function to read in data from screen  
void read_input(int&, int&, double&, double&, double&);
// Function to initialise energy and magnetization
void initialize(int, double, int **, double&, double&);
// The metropolis algorithm 
void Metropolis(int, long&, int **, double&, double&, double *);
// prints to file the results of the calculations  
void output(int, int, double, double *);

//  main program
int main(int argc, char* argv[])
{
  char *outfilename;
  long idum;
  int **spin_matrix, n_spins, mcs;
  double w[17], average[5], initial_temp, final_temp, E, M, temp_step;

  // Read in output file, abort if there are too few command-line arguments
  if( argc <= 1 ){
    cout << "Bad Usage: " << argv[0] << 
      " read also output file on same line" << endl;
    exit(1);
  }
  else{
    outfilename=argv[1];
  }
  ofile.open(outfilename);
  //    Read in initial values such as size of lattice, temp and cycles
  read_input(n_spins, mcs, initial_temp, final_temp, temp_step);
  spin_matrix = (int**) matrix(n_spins, n_spins, sizeof(int));
  idum = -1; // random starting point
  for ( double temp = initial_temp; temp <= final_temp; temp+=temp_step){
    //    initialise energy and magnetization 
    E = M = 0.;
    // setup array for possible energy changes
    for( int de =-8; de <= 8; de++) w[de+8] = 0;
    for( int de =-8; de <= 8; de+=4) w[de+8] = exp(-de/temp);
    // initialise array for expectation values
    for( int i = 0; i < 5; i++) average[i] = 0.;
    initialize(n_spins, double temp, spin_matrix, E, M);
    // start Monte Carlo computation
    for (int cycles = 1; cycles <= mcs; cycles++){
      Metropolis(n_spins, idum, spin_matrix, E, M, w);
      // update expectation values
      average[0] += E;    average[1] += E*E;
      average[2] += M;    average[3] += M*M; average[4] += fabs(M);
    }
    // print results
    output(n_spins, mcs, temp, average);
  }
  free_matrix((void **) spin_matrix); // free memory
  ofile.close();  // close output file
  return 0;
}

Coding energy differences

The array $w[17]$ contains values of $\Delta E$ spanning from $-8J$ to $8J$ and it is precalculated in the main part for every new temperature. The program takes as input the initial temperature, final temperature, a temperature step, the number of spins in one direction (we force the lattice to be a square lattice, meaning that we have the same number of spins in the $x$ and the $y$ directions) and the number of Monte Carlo cycles.

Efficient expression for energy differences

For every Monte Carlo cycle we run through all spins in the lattice in the function metropolis and flip one spin at the time and perform the Metropolis test. However, every time we flip a spin we need to compute the actual energy difference $\Delta E$ in order to access the right element of the array which stores $e^{\beta \Delta E}$. This is easily done in the Ising model since we can exploit the fact that only one spin is flipped, meaning in turn that all the remaining spins keep their values fixed. The energy difference between a state $E_1$ and a state $E_2$ with zero external magnetic field is

 
$$\begin{equation*} \Delta E = E_2-E_1 =J\sum_{< kl>}^{N}s_k^1s_{l}^1-J\sum_{< kl>}^{N}s_k^2s_{l}^2, \end{equation*}$$

 
which we can rewrite as

 
$$\begin{equation*} \Delta E = -J \sum_{< kl>}^{N}s_k^2(s_l^2-s_{l}^1), \end{equation*}$$

 
where the sum now runs only over the nearest neighbors $k$.

Final energy difference

Since the spin to be flipped takes only two values, $s_l^1=\pm 1$ and $s_l^2=\pm 1$, it means that if $s_l^1= 1$, then $s_l^2=-1$ and if $s_l^1= -1$, then $s_l^2=1$. The other spins keep their values, meaning that $s_k^1=s_k^2$. If $s_l^1= 1$ we must have $s_l^1-s_{l}^2=2$, and if $s_l^1= -1$ we must have $s_l^1-s_{l}^2=-2$. From these results we see that the energy difference can be coded efficiently as

 
$$\begin{equation} \Delta E = 2Js_l^1\sum_{< k>}^{N}s_k, \tag{3} \end{equation}$$

 
where the sum runs only over the nearest neighbors $k$ of spin $l$. We can compute the change in magnetisation by flipping one spin as well. Since only spin $l$ is flipped, all the surrounding spins remain unchanged.

Change in magnetization

The difference in magnetisation is therefore only given by the difference $s_l^1-s_{l}^2=\pm 2$, or in a more compact way as

 
$$\begin{equation} M_2 = M_1+2s_l^2, \tag{4} \end{equation}$$

 
where $M_1$ and $M_2$ are the magnetizations before and after the spin flip, respectively. Eqs. (3) and (4) are implemented in the function metropolis shown here

void Metropolis(int n_spins, long& idum, int **spin_matrix, double& E, double&M, double *w)
{
  // loop over all spins
  for(int y =0; y < n_spins; y++) {
    for (int x= 0; x < n_spins; x++){
      // Find random position
      int ix = (int) (ran1(&idum)*(double)n_spins);
      int iy = (int) (ran1(&idum)*(double)n_spins);
      int deltaE =  2*spin_matrix[iy][ix]*
	(spin_matrix[iy][periodic(ix,n_spins,-1)]+
	 spin_matrix[periodic(iy,n_spins,-1)][ix] +
	 spin_matrix[iy][periodic(ix,n_spins,1)] +
	 spin_matrix[periodic(iy,n_spins,1)][ix]);
      // Here we perform the Metropolis test
      if ( ran1(&idum) <= w[deltaE+8] ) {
	spin_matrix[iy][ix] *= -1;  // flip one spin and accept new spin config
        // update energy and magnetization
        M += (double) 2*spin_matrix[iy][ix];
        E += (double) deltaE;
      }
    }
  }
} // end of Metropolis sampling over spins

A small note

Note that we loop over all spins but that we choose the lattice positions $x$ and $y$ randomly. If the move is accepted after performing the Metropolis test, we update the energy and the magnetisation. The new values are used to update the averages computed in the main function.

Initialization

We need also to initialize various variables. This is done in the function here.

// function to initialise energy, spin matrix and magnetization
void initialize(int n_spins, double temp, int **spin_matrix, 
		double& E, double& M)
{
  // setup spin matrix and intial magnetization
  for(int y =0; y < n_spins; y++) {
    for (int x= 0; x < n_spins; x++){
      if (temp < 1.5) spin_matrix[y][x] = 1; // spin orientation for the ground state
      M +=  (double) spin_matrix[y][x];
    }
  }
  // setup initial energy
  for(int y =0; y < n_spins; y++) {
    for (int x= 0; x < n_spins; x++){
      E -=  (double) spin_matrix[y][x]*
	(spin_matrix[periodic(y,n_spins,-1)][x] +
	 spin_matrix[y][periodic(x,n_spins,-1)]);
    }
  }
}// end function initialise

The Metropolis Algorithm and the Two-dimensional Ising Model, elements of program

Here follows an alternative Ising model code using the Mersenne twister engine as described in the c++ "random class":" ".

/* 
   Program to solve the two-dimensional Ising model 
   with zero external field and no parallelization using the Mersenne twister engine for generating random 
   numbers.
   The coupling constant J = 1
   Boltzmann's constant = 1, temperature has thus dimension energy
   Metropolis sampling is used. Periodic boundary conditions.
   The code needs an output file on the command line and the variables mcs, nspins,
   initial temp, final temp and temp step.
   Run as
   ./executable Outputfile numberof spins number of MC cycles initial temp final temp tempstep
   ./test.x Lattice 100 10000000 2.1 2.4 0.01
   Compile and link as 
   c++ -O3 -std=c++11 -Rpass=loop-vectorize -o Ising.x IsingModel.cpp -larmadillo
*/

#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <random>
#include <armadillo>
#include <string>
using namespace  std;
using namespace arma;
// output file
ofstream ofile;

// inline function for periodic boundary conditions
inline int periodic(int i, int limit, int add) { 
  return (i+limit+add) % (limit);
}
// Function to initialise energy and magnetization
void InitializeLattice(int, mat &, double&, double&);
// The metropolis algorithm including the loop over Monte Carlo cycles
void MetropolisSampling(int, int, double, vec &);
// prints to file the results of the calculations  
void output(int, int, double, vec);

// Main program begins here

int main(int argc, char* argv[])
{
  string filename;
  int NSpins, MCcycles;
  double InitialTemp, FinalTemp, TempStep;
  if (argc <= 5) {
    cout << "Bad Usage: " << argv[0] << 
      " read output file, Number of spins, MC cycles, initial and final temperature and tempurate step" << endl;
    exit(1);
  }
  if (argc > 1) {
    filename=argv[1];
    NSpins = atoi(argv[2]);
    MCcycles = atoi(argv[3]);    
    InitialTemp = atof(argv[4]);
    FinalTemp = atof(argv[5]);
    TempStep = atof(argv[6]);
  }
  // Declare new file name and add lattice size to file name
  string fileout = filename;
  string argument = to_string(NSpins);
  fileout.append(argument);
  ofile.open(fileout);
  // Start Monte Carlo sampling by looping over T first
  for (double Temperature = InitialTemp; Temperature <= FinalTemp; Temperature+=TempStep){
    vec ExpectationValues = zeros<mat>(5);
    // start Monte Carlo computation
    MetropolisSampling(NSpins, MCcycles, Temperature, ExpectationValues);
    output(NSpins, MCcycles, Temperature, ExpectationValues);
  }
  ofile.close();  // close output file
  return 0;
}

// function to initialise energy, spin matrix and magnetization
void InitializeLattice(int NSpins, mat &SpinMatrix,  double& Energy, double& MagneticMoment)
{
  // setup spin matrix and initial magnetization
  for(int x =0; x < NSpins; x++) {
    for (int y= 0; y < NSpins; y++){
      SpinMatrix(x,y) = 1.0; // spin orientation for the ground state
      MagneticMoment +=  (double) SpinMatrix(x,y);
    }
  }
  // setup initial energy
  for(int x =0; x < NSpins; x++) {
    for (int y= 0; y < NSpins; y++){
      Energy -=  (double) SpinMatrix(x,y)*
	(SpinMatrix(periodic(x,NSpins,-1),y) +
	 SpinMatrix(x,periodic(y,NSpins,-1)));
    }
  }
}// end function initialise


// The Monte Carlo part with the Metropolis algo with sweeps over the lattice
void MetropolisSampling(int NSpins, int MCcycles, double Temperature, vec &ExpectationValues)
{
  // Initialize the seed and call the Mersenne algo
  std::random_device rd;
  std::mt19937_64 gen(rd());
  // Then set up the uniform distribution for x \in [[0, 1]
  std::uniform_real_distribution<double> distribution(0.0,1.0);
  //  Allocate memory for spin matrix
  mat SpinMatrix = zeros<mat>(NSpins,NSpins);
  //    initialise energy and magnetization 
  double Energy = 0.;     double MagneticMoment = 0.;
  // initialize array for expectation values
  InitializeLattice(NSpins, SpinMatrix, Energy, MagneticMoment);
  // setup array for possible energy changes
  vec EnergyDifference = zeros<mat>(17); 
  for( int de =-8; de <= 8; de+=4) EnergyDifference(de+8) = exp(-de/Temperature);
  for (int cycles = 1; cycles <= MCcycles; cycles++){
    // The sweep over the lattice, looping over all spin sites
    for(int x =0; x < NSpins; x++) {
      for (int y= 0; y < NSpins; y++){
	int ix = (int) (distribution(gen)*(double)NSpins);
	int iy = (int) (distribution(gen)*(double)NSpins);
	int deltaE =  2*SpinMatrix(ix,iy)*
	  (SpinMatrix(ix,periodic(iy,NSpins,-1))+
	   SpinMatrix(periodic(ix,NSpins,-1),iy) +
	   SpinMatrix(ix,periodic(iy,NSpins,1)) +
	   SpinMatrix(periodic(ix,NSpins,1),iy));
	if ( distribution(gen) <= EnergyDifference(deltaE+8) ) {
	  SpinMatrix(ix,iy) *= -1.0;  // flip one spin and accept new spin config
	  MagneticMoment += (double) 2*SpinMatrix(ix,iy);
	  Energy += (double) deltaE;
	}
      }
    }
    // update expectation values  for local node
    ExpectationValues(0) += Energy;    ExpectationValues(1) += Energy*Energy;
    ExpectationValues(2) += MagneticMoment;    
    ExpectationValues(3) += MagneticMoment*MagneticMoment; 
    ExpectationValues(4) += fabs(MagneticMoment);
  }

} // end of Metropolis sampling over spins

void output(int NSpins, int MCcycles, double temperature, vec ExpectationValues)
{
  double norm = 1.0/((double) (MCcycles));  // divided by  number of cycles 
  double E_ExpectationValues = ExpectationValues(0)*norm;
  double E2_ExpectationValues = ExpectationValues(1)*norm;
  double M_ExpectationValues = ExpectationValues(2)*norm;
  double M2_ExpectationValues = ExpectationValues(3)*norm;
  double Mabs_ExpectationValues = ExpectationValues(4)*norm;
  // all expectation values are per spin, divide by 1/NSpins/NSpins
  double Evariance = (E2_ExpectationValues- E_ExpectationValues*E_ExpectationValues)/NSpins/NSpins;
  double Mvariance = (M2_ExpectationValues - Mabs_ExpectationValues*Mabs_ExpectationValues)/NSpins/NSpins;
  ofile << setiosflags(ios::showpoint | ios::uppercase);
  ofile << setw(15) << setprecision(8) << temperature;
  ofile << setw(15) << setprecision(8) << E_ExpectationValues/NSpins/NSpins;
  ofile << setw(15) << setprecision(8) << Evariance/temperature/temperature;
  ofile << setw(15) << setprecision(8) << M_ExpectationValues/NSpins/NSpins;
  ofile << setw(15) << setprecision(8) << Mvariance/temperature;
  ofile << setw(15) << setprecision(8) << Mabs_ExpectationValues/NSpins/NSpins << endl;
} // end output function



    

Two-dimensional Ising Model, energy per spin and specific heat

The following Python program, based on the above C++ codes, plots the expectation value of the energy and its fluctuation, that is the specific heat. Both quantities are plotted per spin and genererated for a $20\times 20$ lattice.

                 (spin_matrix.item(periodic(i,size,-1),j) + spin_matrix.item(i,periodic(j,size,1)))

Two-dimensional Ising Model and analysis of spin values

The following python code displays the values of the spins as function of temperature. The blue color corresponds to spin up states while red represents spin down states. Increasing the temperature as input parameter, see the parameters below, results in a a net magnetization which becomes zero. At low temperatures, the system is highly ordered with essentially only one specific spin value.

print "%15.8E %15.8E %15.8E %15.8E %15.8E %15.8E\n" % (temp, E_av, E_variance, M_av, M_variance, Mabs_av)