304 The Monte Carlo method
10.3.1 Monte Carlo for the Ising modelThe Ising model was discussed in Chapter 7. Here we consider the Metropolis
algorithm for the Ising model in two dimensions (the extension to more than two
dimensions is straightforward). In order to formulate the Monte Carlo method, we
must make a choice for the matrixωXX′. For the two-dimensional Ising model on
anL×Lsquare lattice, we take
ωXX′= 1 /L^2 ifXandX′differ by one spin;
ωXX′= 0 otherwise.(10.23)
The realisation of the first stage of the Markov step – generating a trial configur-
ation – is then easy: we select a spin at random, and the trial configuration is the
present configuration with the selected spin turned over.
We then calculate the energy differenceE(X→X′)between the old and the
trial configuration:
E(X→X′)=E(X′)−E(X). (10.24)
This is easy as there are only nearest neighbour interactions: the energy difference
depends only on the number of neighbours which have the same spin as the selected
spin in the old configuration – this is therefore an integer number between 0 and 4.
If the energy increases from the old to the new configuration,E(X→X′)>0,
the trial state is accepted with probability exp[−βE(X→X′)]. If, however, the
energy decreases, the trial state is always accepted as the new state. Implementation
of this algorithm for the two- or three-dimensional Ising model is straightforward;
the reader is encouraged to go through this exercise. The average number of steps
between two updates of the same spin is equal toL^2. Therefore, the ‘time’ in a MC
simulation is often expressed in units ofMonte Carlo steps per spin(MCS), 1 MCS
being equal toL^2 trials.
programming exercise
Write a program for simulating the nearest neighbour two-dimensional Ising
model on a square lattice using the Metropolis Monte Carlo technique.
The following considerations should be borne in mind. It is convenient to have
a variable representing the total energy of the system. Usually,βEis recorded
rather than the energyEitself, as the behaviour of the system is fully determined
byβJ. This must be calculated at the beginning by performing a sweep over the
whole lattice. If, however, the initial state is one in which all the spins are the same,
the total energy for theL×Llattice with only nearest neighbour interactions is
simply given by− 2 βJL^2 , whereβJis the coupling constant. Every time a trial
configuration is accepted as the new one, we add the energy difference (which is