In the introductory lecture, I introduced state-counting using a simple, one-dimensional Ising system. To be precise, this was a particular case known as a Ising system of noninteracting spins, because each spin made its own independent contribution to the Hamiltonian.
Now, we will consider a more interesting Ising system; namely, that of interacting spins on a 2d square lattice. What do we mean by interacting? We mean that the Hamiltonian depends on pairs of spins:
where the notation \(\left \langle ij\right \rangle \) denotes that we are summing over unique pairs of nearest neighbors, and, as before a spin \(s_i\) is either +1 or -1. How many unique pairs of nearest neighbors are there on a lattice of \(N\) spins (assuming periodic boundary conditions)? To answer this question, we must know the coordination, \(z\), of the lattice; that is, how many nearest neigbhors does one lattice position have? For a square lattice, \(z\) = 4. Each spin therefore contributes two unique spin pairs to the system, so there are \(Nz/2\) unique nearest neighbor pairs.
Think about this Hamiltonian. It says that energy is minimal when all \(N\) spins have the same alignment, either all up or all down. Imagine a microscopic observable called the magnetization, or average spin orientation, \(M\): \begin{equation} M \equiv \frac {1}{N} \sum _i^N s_i \end{equation} Then, \(M\) = +1 and -1 are two energetically equivalent microstates. We can then expect that if there is any thermal energy in the system which randomly flips spins, the amount of this thermal energy (that is, the temerature) will somehow control the observable magnetization. We “observe” magnetization using an ensemble average: \begin{equation} \label {eq:magnetization} \left \langle M\right \rangle = \frac {\sum _{s_1=-1}^{+1} \cdots \sum _{s_N=-1}^{+1} \left [\sum _i^N s_i\right ]\exp \left [-\beta \mathscr {H}\left (s_1,\dots ,s_N\right )\right ]}{\sum _{s_1=-1}^{+1} \cdots \sum _{s_N=-1}^{+1} \exp \left [-\beta \mathscr {H}\left (s_1,\dots ,s_N\right )\right ]} \end{equation}
Our physical intution tells us that as \(T\rightarrow \) 0, \(\left \langle M\right \rangle \rightarrow \) 1, and as \(T\rightarrow \infty \), \(\left \langle M\right \rangle \rightarrow \) 0. The fascinating thing about an Ising magnet is that there is a finite temperature called the critical temperature, \(T_c\). If we start out with a “hot” system, and cool it to just below \(T_c\), the absolute value of the magnetization spontaneously jumps from 0 to some finite positive value. In other words, the system undergoes a phase transition from a disordered phase to a partially ordered phase. A Metropolis Monte Carlo simulation can allow us to probe the behavior of an Ising system and learn how the system behaves near criticality. The rest of this case study will be devoted to showing you the inner workings of a C program which simulates the Ising lattice using Metropolis MC, as a first implementation of this technique. In the suggested exercises appearing at the end of this case study, you will modify this code slightly to compute averages values of certain observables.
But first, I recommend a visit to the website of Peter Young at UC Santa Cruz. You will see a Java implentation of a Monte Carlo simulation of a 2-dimensional Ising magnet (One of many on the web; google “ising simulation” and you’ll get a nice sample.) This is a fun little Java applet that lets you play with an Ising system. You can change the temperature of the simulation: making it cold will “freeze” the system, and making it hot “melts” it. Near the critical temperature, \(T_c\), relatively large regions of mostly-up spins compete with regions of mostly down spin. In one of the exercises, we’ll learn how to measure an observable called the correlation length, which characterizes the size of these domains and is a useful signature of criticality.
In this section, we will dissect piece-by-piece a small program (written in C) which implements an NVT Metropolis Monte Carlo simulation of a 2D Ising lattice. Click here to download the code. You can compile the code using the command
cfa@abrams01:/home/cfa> gcc -O3 -o ising ising.c -lm -lgsl
and you can then run it at the command line as ./ising. The code will conduct a canoncial Metropolis Monte Carlo simulation of an Ising lattice of size \(L\times L\) at temperature \(T\) (both specified by the user at run time on the command line), and it computes both the average energy per spin \(\left \langle \epsilon \right \rangle \) and the average spin value, \(\left \langle s_1\right \rangle \). Periodic boundaries are employed in calculating the nearest-neighbor interactions.
An abbreviated listing of the code follows. Some comments in the full, downloadable code have been omitted for space, and I have instead explained each code fragment in accompanying text.
_________________________________________________________________________________________
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <gsl/gsl_rng.h> Standard headers for a C program, and the header for the GSL random number generator. |
|
int E ( int ** F, int L, int i, int j ) {
return -2*(F[i][j])*
(F[i?(i-1):(L-1)][j]+F[(i+1)%L][j]+
F[i][j?(j-1):(L-1)]+F[i][(j+1)%L]);
}
A function that accepts as arguments the 2D array of spins, F; the side-length, L, and a position (i,j), and returns the change in energy upon flipping spin (i,j), without actually flipping it. All variables are of integer type, int. The syntax ** F means that F is a pointer to a pointer to an int. It is a way of signifying that F is a 2D array. To access the i,j element of F, the syntax is F[i][j]. The syntax i?(i-1):(L-1) expands as, “If i is non-zero, return i-1; otherwise (if i is zero), return L-1.” This implements periodic boundaries when looking to the west of spin i,j. The syntax (i+1)%L returns (i+1) mod L, which also implements periodic boundaries when looking to the east of spin i,j. |
|
void samp ( int ** F, int L, double * s, double * e ) {
int i,j;
*s=0.0;
*e=0.0;
for (i=0;i<L;i++) {
for (j=0;j<L;j++) {
*s+=(double)F[i][j];
*e-=(double)(F[i][j])*
(F[i][(j+1)%L]+F[(i+1)%L][j]);
}
}
*s/=(L*L);
*e/=(L*L);
}
A function that samples the current 2D array of spins and computes the average spin, \(\left \langle s_1\right \rangle \), and the average energy per spin, \(\left \langle \epsilon \right \rangle \). The syntax x += y expands as x = x + y; the same is true for other “incremental” operators, -=, *=, and /=. Because s and e are “passed by reference” (so that their values can be changed by the function), we have to dereference them with the * operator to access their contents; that is *s means “the contents of s”. |
|
void init ( int ** F, int L, gsl_rng * r ) {
int i,j;
for (i=0;i<L;i++) {
for (j=0;j<L;j++) {
F[i][j]=2*(int)gsl_rng_uniform_int(r,2)-1;
}
}
}
A function that randomly initializes the 2D array of spins. The function gsl_rng_uniform_int(r,2) returns either 0 or 1 randomly, and we want each spin to be either 1 or -1. Note that we have to pass in the random number generator we created (as you will see) in the main program (below). |
|
int main (int argc, char * argv[]) {
Main program begins here. |
|
int ** F; int L = 20; int N; double T = 1.0; F is the the 2D array of spins; L is the sidelength of the array (default value is 20); N is the total number of spins = L\(^2\); T is the dimensionless temperature = \(k_BT/J\) (default is 1.0), where \(J\) is the unit energy of the Hamiltonian. |
|
int nCycles = 1000000; int fSamp = 1000; The number of cycles to run and the sample interval. A cycle is a set of \(N\) attempted spin flips. |
|
int nSamp; int de; double b; double x; int i,j,a,c; Computed variables: number of samples taken, change in energy upon spin flip, Boltzmann factor, random number, and loop counters. |
|
double s=0.0, ssum=0.0; double e=0.0, esum=0.0; Observables, average spin, \(s_1\), and average energy per spin, \(\epsilon \), and their respective accumulators for their ensemble averages, all initialized to 0. |
|
gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937); unsigned long int Seed = 23410981; This line creates a random number generator of the ”Mersenne Twister” type, which is much better than the default random number generator. This is why we need the GNU Scientific Library. |
|
for (i=1;i<argc;i++) {
if (!strcmp(argv[i],"-L"))
L=atoi(argv[++i]);
else if (!strcmp(argv[i],"-T"))
T=atof(argv[++i]);
else if (!strcmp(argv[i],"-nc"))
nCycles = atoi(argv[++i]);
else if (!strcmp(argv[i],"-fs"))
fSamp = atoi(argv[++i]);
else if (!strcmp(argv[i],"-s"))
Seed = (unsigned long)atoi(argv[++i]);
}
Here we parse the command line arguments. |
|
gsl_rng_set(r,Seed); Seed the random number generator. |
|
N=L*L; F=(int**)malloc(L*sizeof(int*)); for (i=0;i<L;i++) F[i]=(int*)malloc(L*sizeof(int)); Compute the number of spins, \(N\), and allocate memory for the 2D array of spins. |
|
init(F,L,r); T=1.0/T; Initialize the 2D array of spins, and convert the temperature to reciprocal temperature for computational convenience. |
|
for (c=0;c<nCycles;c++) {
Let \(c\) loop from 0 to nCycles-1. |
|
for (a=0;a<N;a++) {
Let \(a\) loop from 0 to N-1. |
|
i=(int)gsl_rng_uniform_int(r,L); j=(int)gsl_rng_uniform_int(r,L); Randomly select spin \((i,j)\). The function gsl_rng_uniform_int() is a from the GSL. |
|
de = E(F,L,i,j); b = exp(de*T); x = gsl_rng_uniform(r); Compute \(\Delta \mathscr {U}\), \(e^{-\beta \Delta \mathscr {U}}\), and select a uniform random variate, \(x\). The function gsl_rng_uniform() is from the GSL. |
|
if (x<b) {
F[i][j]*=-1;
}
Evaluate the acceptance test, and if passed, actually flip the spin. The last 6 lines of code can be equivalently stated with just one line: if (gsl_rng_uniform(r)<exp(E(F,L,i,j)*T)) F[i][j]*=-1;
|
|
} Close the loop over \(N\) spin flip trials. |
|
if (!(c%fSamp)) {
samp(F,L,&s,&e);
ssum+=s;
esum+=e;
nSamp++;
}
If the cycle number is divisible by the requested sample frequency, (as determined by asking whether c mod fSamp is zero), take a sample, and update the accumulators. |
|
} Close the outer loop over nCycle cycles. |
|
fprintf(stdout, "# The average magnetization is %.5lf\n", ssum/nSamp); fprintf(stdout, "# The average energy per spin is %.5lf\n", esum/nSamp); fprintf(stdout,"# Program ends.\n"); } Output final information. |