Morten Hjorth-Jensen Email morten.hjorth-jensen@fys.uio.no, Department of Physics and Center of Mathematics for Applications, University of Oslo and National Superconducting Cyclotron Laboratory, Michigan State University
Date: Fall 2015
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 |
| Energy | Internal | Internal | Internal | Enthalpy |
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
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
and the pressure by
Similarly we can compute the chemical potential as
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
or using the probability distribution $P_i$ as
If we divide the latter quantity with $kT^2$ we obtain the specific heat at constant volume
and the corresponding variance
This quantity defines also the susceptibility $\chi$
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
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
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
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
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];
}The other configurations give
and
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
and we obtain the following expression for the two-spin case
The other cases do also differ and we have
and
jm=N;
for ( j=1; j <=N ; j++) {
energy += spin[j]*spin[jm];
jm = j ;
}The magnetization is however the same, defined as
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
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. Strong correlations make a perturbative treatment impossible. From a theoretical point of view, the way out is renormalization group theory. The table lists some typical system with their pertinent order parameters.
| 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
meaning $ lnZ = -F/kT = -F\beta$. The energy is given as the first derivative of $F$
and the specific heat is defined via the second derivative of $F$
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.
Critical phenomena have been extensively studied in physics. One major reason is that we still do not have a satisfactory understanding of the properties of a system close to a critical point. Even for the simplest three-dimensional systems we cannot predict exactly the values of various thermodynamical variables. Simplified theoretical approaches like mean-field models discussed below, can even predict the wrong physics. Mean-field theory results in a second-order phase transition for the one-dimensional Ising model, whereas we saw in the previous section that the one-dimensional Ising model does not predict any spontaneous magnetization at any finite temperature. The physical reason for this can be understood from the following simple consideration. Assume that the ground state for an $N$-spin system in one dimension is characterized by the following configuration
It is possibleto show that the mean magnetization is given by (for temperature below $T_C$)
where $\beta=1/8$ is a so-called critical exponent. A similar relation applies to the heat capacity
and the susceptibility
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
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. Stated differently, how can the spins propagate their correlations so extensively when we approach the critical point, in particular since the interaction acts only between nearest spins? Below we will compute the correlation length via the spin-sin correlation function for the one-dimensional Ising model.
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
with $a$ a constant and $\nu$ defined in Eq. (1).
The correlation length for a finite lattice size can then be shown to be proportional to
and if we set $T=T_C$ one can obtain the following relations for the magnetization, energy and susceptibility for $T \le T_C$
and
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
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
Establish an initial state with energy $E_b$ by positioning yourself at a random configuration in the lattice
Change the initial configuration by flipping e.g., one spin only. Compute the energy of this trial state $E_t$.
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.
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.
If $\Delta E > 0$, calculate $w=e^{-(\beta \Delta E)}$.
Compare $w$ with a random number $r$. If
then accept the new configuration, else we keep the old configuration.
The next step is to update various expectations values.
The steps (2)-(7) are then repeated in order to obtain a sufficently good representation of states.
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 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$.
The four other possibilities are as follows
with $\Delta E=4J$,
with $\Delta E=0$,
with $\Delta E=-4J$ and finally
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
which we can rewrite as
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
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
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 spinsNote 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
In [ ]:
# coding=utf-8
#2-dimensional ising model with visualization
import numpy, sys, math
import pygame
#Needed for visualize when using SDL
screen = None;
font = None;
BLOCKSIZE = 10
def periodic (i, limit, add):
"""
Choose correct matrix index with periodic
boundary conditions
Input:
- i: Base index
- limit: Highest \"legal\" index
- add: Number to add or subtract from i
"""
return (i+limit+add) % limit
def visualize(spin_matrix, temp, E, M, method):
"""
Visualize the spin matrix
Methods:
method = -1:No visualization (testing)
method = 0: Just print it to the terminal
method = 1: Pretty-print to terminal
method = 2: SDL/pygame single-pixel
method = 3: SDL/pygame rectangle
"""
#Simple terminal dump
if method == 0:
print "temp:", temp, "E:", E, "M:", M
print spin_matrix
#Pretty-print to terminal
elif method == 1:
out = ""
size = len(spin_matrix)
for y in xrange(size):
for x in xrange(size):
if spin_matrix.item(x,y) == 1:
out += "X"
else:
out += " "
out += "\n"
print "temp:", temp, "E:", E, "M:", M
print out + "\n"
#SDL single-pixel (useful for large arrays)
elif method == 2:
size = len(spin_matrix)
screen.lock()
for y in xrange(size):
for x in xrange(size):
if spin_matrix.item(x,y) == 1:
screen.set_at((x,y),(0,0,255))
else:
screen.set_at((x,y),(255,0,0))
screen.unlock()
pygame.display.flip()
#SDL block (usefull for smaller arrays)
elif method == 3:
size = len(spin_matrix)
screen.lock()
for y in xrange(size):
for x in xrange(size):
if spin_matrix.item(x,y) == 1:
rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)
pygame.draw.rect(screen,(0,0,255),rect)
else:
rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)
pygame.draw.rect(screen,(255,0,0),rect)
screen.unlock()
pygame.display.flip()
#SDL block w/ data-display
elif method == 4:
size = len(spin_matrix)
screen.lock()
for y in xrange(size):
for x in xrange(size):
if spin_matrix.item(x,y) == 1:
rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)
pygame.draw.rect(screen,(255,255,255),rect)
else:
rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)
pygame.draw.rect(screen,(0,0,0),rect)
s = font.render("<E> = %5.3E; <M> = %5.3E" % E,M,False,(255,0,0))
screen.blit(s,(0,0))
screen.unlock()
pygame.display.flip()
def monteCarlo(temp, size, trials, visual_method):
"""
Calculate the energy and magnetization
(\"straight\" and squared) for a given temperature
Input:
- temp: Temperature to calculate for
- size: dimension of square matrix
- trials: Monte-carlo trials (how many times do we
flip the matrix?)
- visual_method: What method should we use to visualize?
Output:
- E_av: Energy of matrix averaged over trials, normalized to spins**2
- E_variance: Variance of energy, same normalization * temp**2
- M_av: Magnetic field of matrix, averaged over trials, normalized to spins**2
- M_variance: Variance of magnetic field, same normalization * temp
- Mabs: Absolute value of magnetic field, averaged over trials
"""
#Setup spin matrix, initialize to ground state
spin_matrix = numpy.zeros( (size,size), numpy.int8) + 1
#Create and initialize variables
E = M = 0
E_av = E2_av = M_av = M2_av = Mabs_av = 0
#Setup array for possible energy changes
w = numpy.zeros(17,numpy.float64)
for de in xrange(-8,9,4): #include +8
w[de+8] = math.exp(-de/temp)
#Calculate initial magnetization:
M = spin_matrix.sum()
#Calculate initial energy
for j in xrange(size):
for i in xrange(size):
E -= spin_matrix.item(i,j)*\
(spin_matrix.item(periodic(i,size,-1),j) + spin_matrix.item(i,periodic(j,size,1)))
#Start metropolis MonteCarlo computation
for i in xrange(trials):
#Metropolis
#Loop over all spins, pick a random spin each time
for s in xrange(size**2):
x = int(numpy.random.random()*size)
y = int(numpy.random.random()*size)
deltaE = 2*spin_matrix.item(x,y)*\
(spin_matrix.item(periodic(x,size,-1), y) +\
spin_matrix.item(periodic(x,size,1), y) +\
spin_matrix.item(x, periodic(y,size,-1)) +\
spin_matrix.item(x, periodic(y,size,1)))
if numpy.random.random() <= w[deltaE+8]:
#Accept!
spin_matrix[x,y] *= -1
M += 2*spin_matrix[x,y]
E += deltaE
#Update expectation values
E_av += E
E2_av += E**2
M_av += M
M2_av += M**2
Mabs_av += int(math.fabs(M))
visualize(spin_matrix, temp,E/float(size**2),M/float(size**2), method);
#Normalize average values
E_av /= float(trials);
E2_av /= float(trials);
M_av /= float(trials);
M2_av /= float(trials);
Mabs_av /= float(trials);
#Calculate variance and normalize to per-point and temp
E_variance = (E2_av-E_av*E_av)/float(size*size*temp*temp);
M_variance = (M2_av-M_av*M_av)/float(size*size*temp);
#Normalize returned averages to per-point
E_av /= float(size*size);
M_av /= float(size*size);
Mabs_av /= float(size*size);
return (E_av, E_variance, M_av, M_variance, Mabs_av)
# Main program
size = 100
trials = 100000
temp = 2.5
method = 3
#Initialize pygame
if method == 2 or method == 3 or method == 4:
pygame.init()
if method == 2:
screen = pygame.display.set_mode((size,size))
elif method == 3:
screen = pygame.display.set_mode((size*10,size*10))
elif method == 4:
screen = pygame.display.set_mode((size*10,size*10))
font = pygame.font.Font(None,12)
(E_av, E_variance, M_av, M_variance, Mabs_av) = monteCarlo(temp,size,trials, method)
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)
pygame.quit();
In [ ]: