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.
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.
| Microcanonical | Canonical | Grand Canonical | Pressure canonical | |
|---|---|---|---|---|
| Exchange of heat | no | yes | yes | yes |
| with the environment | ||||
| Exchange of particles | no | no | yes | no |
| with the environemt | ||||
| Thermodynamical | \( V, \cal M, \cal D \) | \( V, \cal M, \cal D \) | \( V, \cal M, \cal D \) | \( P, \cal H, \cal E \) |
| parameters | \( E \) | \( T \) | \( T \) | \( T \) |
| \( N \) | \( N \) | \( \mu \) | \( N \) | |
| Potential | Entropy | Helmholtz | \( PV \) | Gibbs |
| \( N \) | \( N \) | \( \mu \) | \( N \) | |
| Energy | Internal | Internal | Internal | Enthalpy |
| \( N \) | \( N \) | \( \mu \) | \( N \) | |
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
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 \).
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.
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*} $$
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*} $$
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.
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*} $$
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.
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 \).
The energy for a specific configuration \( i \) is given by $$ \begin{equation*} E_i =-J\sum_{< kl>}^{N}s_ks_l. \end{equation*} $$
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*} $$
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*} $$
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 \).
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*} $$
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*} $$
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*} $$
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 \).
The table lists the energy and magnetization for both free ends and periodic boundary conditions.
| State | Energy (FE) | Energy (PBC) | Magnetization |
|---|---|---|---|
| \( 1= \uparrow\uparrow \) | \( -J \) | \( -2J \) | 2 |
| \( 2=\uparrow\downarrow \) | \( J \) | \( 2J \) | 0 |
| $ 3=\downarrow\uparrow$ | \( J \) | \( 2J \) | 0 |
| $ 4=\downarrow\downarrow$ | \( -J \) | \( -2J \) | -2 |
We can reorganize according to the number of spins pointing up, as shown in the table here
| Number spins up | Degeneracy | Energy (FE) | Energy (PBC) | Magnetization |
|---|---|---|---|---|
| 2 | 1 | \( -J \) | \( -2J \) | 2 |
| 1 | 2 | \( J \) | \( 2J \) | 0 |
| 0 | 1 | \( -J \) | \( -2J \) | -2 |
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?
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*} $$
In the table here we group these configurations according to their total energy and magnetization.
| Number spins up | Degeneracy | Energy | Magnetization |
|---|---|---|---|
| 4 | 1 | \( -8J \) | 4 |
| 3 | 4 | \( 0 \) | 2 |
| 2 | 4 | \( 0 \) | 0 |
| 2 | 2 | \( 8J \) | 0 |
| 1 | 4 | \( 0 \) | -2 |
| 0 | 1 | \( -8J \) | -4 |
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.
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.
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 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.
| System | Transition | Order Parameter |
| Liquid-gas | Condensation/evaporation | Density difference \( \Delta\rho=\rho_{liquid}-\rho_{gas} \) |
| Binary liquid | mixture/Unmixing | Composition difference |
| Quantum liquid | Normal fluid/superfluid | \( < \phi> \), \( \psi \) = wavefunction |
| Liquid-solid | Melting/crystallisation | Reciprocal lattice vector |
| Magnetic solid | Ferromagnetic | Spontaneous magnetisation \( M \) |
| Antiferromagnetic | Sublattice magnetisation \( M \) | |
| Dielectric solid | Ferroelectric | Polarization \( P \) |
| Antiferroelectric | Sublattice polarisation \( P \) |
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*} $$
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 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.
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 \).
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}. \label{eq:xi} \end{equation} $$
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.
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}, \label{eq:tc} \end{equation} $$ with \( a \) a constant and \( \nu \) defined in Eq. \eqref{eq:xi}.
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*} $$
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.
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
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.
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 \).
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;
}
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.
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 \).
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, \label{eq:deltaeising} \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.
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, \label{eq:deltamising} \end{equation} $$ where \( M_1 \) and \( M_2 \) are the magnetizations before and after the spin flip, respectively. Eqs. \eqref{eq:deltaeising} and \eqref{eq:deltamising} 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
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.
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
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
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.
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.