MATE60002/70002: Theory and Simulation of Materials
Assignment 3: Metropolis Monte Carlo
February 2022
1 Introduction
In the lecture accompanying this workshop and assignment, I discuss the thermodynamics and
statistical physics (aka statistical mechanics) of materials. You will need to understand
some of the fundamentals of these topics to complete the assessed exercise, which involves
investigating ferroic order in crystals using a well-known toy model, known as the Ising
model, and a computational tool known as Metropolis Monte Carlo.
On Blackboard you can find a set of notes written by Andrew Horsfield, which provides a
brief recap of some key concepts in statistical mechanics and thermodynamics and explains
the Metropolic Monte Carlo method.
Monte Carlo methods
Metropolis Monte Carlo is one of many varieties of so-called
Monte Carlo simulation. ‘Monte Carlo’ is a reference
to the famous casinos in Monte Carlo and Monte Carlo
methods are methods that use random numbers to calculate
quantities of interest. For example, on the right there is
a picture of a circle of area pir2 inscribed in a square of
area 4r2. The dots mark points that have been chosen
at random. If there were a total of N these randomly-
chosen points, of which M of them were within the pink
circular region, we could calculate pi from the relation∗
M/N =
(
pir2
)
/
(
4r2
)
= pi/4.
∗You can try this yourself with a short program that calls a random number generator to choose the green
points. Or go to a pub with a dartboard, draw a square around it, and imbibe to randomise your throws.
2 MATE60002/70002: TSM: Paul Tangney
(Anti)Ferroic order
A ferroic crystal is a crystal in which each primitive unit cell
possesses some energetically-important property or structural feature that can be quantified
by a vector, and in which the crystal’s potential energy depends on the degrees to which these
vectors are aligned (↑↑) or antialigned (↑↓). Alignment of the vectors is known as ferroic
order and antialignment of the vectors is known as antiferroic order. In a ferromagnet
(e.g., Fe, Ni, or Co) or an antiferromagnet (e.g., α−Fe2O3) the vectors are magnetic
moments; in a ferroelectric (e.g., BaTiO3) or antiferroelectric (e.g., PbZrO3) they are electric
dipole moments; and in a ferroelastic material, they are strains of the primitive unit cells
along particular crystal axes. Multiferroics are materials that possess more than one form
of ferroic order, which tend to be coupled. For example, BaTiO3 is both ferroelectric and
ferroelastic and the dipole moment and strain of each primitive cell are intimately coupled to
one another.
Ising model
When we don’t fully understand a phenomenon that occurs in some class of materials, and we
want to use theory or simulation to understand its gross features better, we sometimes use so-
called toy models. A toy model is a simplified and idealized representation of one aspect of
the material’s structure and energetics. Its purpose is to provide a crude basic understanding
of how that one aspect of the material depends on initial conditions, boundary conditions,
and/or external stimuli. For example, to study ferroic order we will use a well-known toy model
called the Ising model, in which the complicated interdependent spatiotemporal distributions
of nuclei and electrons are replaced with a single vector at each lattice point of a crystal. The
vector at the i th lattice point is denoted by ~si and we assume that the energy of its interaction
with the vector ~sj at the j th lattice point is U intij ≡ −J ~si ·~sj , where J > 0 if there is a tendency
for the arrows to align (ferroic order), and J < 0 if there is a tendency for them to antialign
(antiferroic order). In addition, there might be an externally-applied field ~h that couples to
each vector ~si with a potential energy of Uexti = −~h ·~si , where ~h could be an electric field,
a magnetic field, and/or a strain, depending on what the set of vectors {~si} represents and
which kind of ferroic order is being investigated. The negative sign of −~h ·~si simply means
that the direction of ~h has been defined as the direction of each ~si for which the energy Uexti of
its coupling to the external field is minimized. The total potential energy of the Ising model
when the set of vectors defining the state of the system is s ≡ (~s1,~s2, · · · ) can be expressed as
H(s) ≡ −J
∑
〈i ,j〉
~si ·~sj − ~h ·
∑
i
~si , (1)
where the first sum is restricted to pairs 〈i , j〉 that are close enough to interact with one
another.
Replacing all nuclei and electrons in a unit cell with a single arrow is an enormous simplification
of the crystal’s structure and energetics, but it can be very useful: it can provide us with
valuable ideas about how ferroic materials behave. We can then use these ideas to develop a
theory, or to build a more accurate or more appropriate model for use in computer simulations.
MATE60002/70002: TSM: Paul Tangney 3
2 Metropolic Monte Carlo simulations of the 2-d Ising Model
We will study a two dimensional square lattice of spins
under periodic boundary conditions. The sides of the
simulation cell, and the vectors joining each spin to its
four nearest neighbours, are all parallel to the x or y
axes. Each spin can be in one of four directions, which
are along the two axes that are rotated by 45◦ from
the x and y axes.
The energy is calculated from Eq. 1, where ~h is a vector
defining a uniform applied field, and J is a scalar whose
magnitude is set to one for most purposes, because all
energies (e.g., kBT = kT) are expressed in units of |J|.
When J = 1 the system is ferroic and when J = −1
it is antiferroic.
The program that you will use is called Ising2dMMC.m
and you can find it on Blackboard. It generates initial
structures randomly, unless told to restart from a
previous run, and then it performs nSweep sweeps,
where each sweep consists of nStepsPerSweep steps, and each step is an attempt to change
the sign of one randomly-chosen spin component. The structure evolves as more spins are
flipped and, even if the initial microstate is one that is highly improbable under the physical
conditions being simulated, as the simulation progresses it will move into microstates that
are more probable and therefore more physical. The Metropolis algorithm is described briefly
below, and much more detail can be found online. There is also a great deal of information
available online about the Ising model.
Metropolis algorithm
Step 1: Generate an initial microstate (i.e., assign directions to all the spins) and set isweep = 0.
Step 2: If isweep=nSweeps, then skip forward to Step 8. Otherwise, increase the value of
isweep by 1 and set istep= 0.
Step 3: If istep=nStepsPerSweep, then return to Step 2. Otherwise, increase the value of
istep by 1.
Step 4: Select one component of one of the spins at random.
Step 5: Calculate the change in energy, ∆E , when the spin is flipped.
Step 6: If either ∆E ≤ 0 or ∆E > 0 and 〈−∆E/kBT 〉 < R, where R is a random number
between 0 and 1, then flip the spin. Otherwise do not flip the spin.
Step 7: Return to Step 3.
Step 8: End the simulation.
4 MATE60002/70002: TSM: Paul Tangney
Assignment 3
This will be marked out of 50. Task 1 is worth 30 marks and Task 2 is worth 20 marks.
As in previous assignments, you will be marked on presentation, on clarity, on how well
you understand the theory behind the simulations, on how well you perform the simulations
(equilibrate!), on how appropriately you interpret the results, and on how well you identify
and quantify sources of significant uncertainty.
Task 1
For this task you should submit nine graphs and three sentences.
Simulate each of the following three systems over an interesting or appropriate range of
temperatures:
1. A 10× 10 lattice with J = 1 and ~h = 0 (variable hmag= 0).
2. A 10× 10 lattice with J = −1 and ~h = 0.
3. A 10× 10 lattice with J = 1 and |~h|/J = 0.1, with ~h parallel to the horizontal axis.
In each case, create the following three plots:
A. Magnetization magnitude | ~M| (variable MagMag) versus kBT/J (internal variable kT,
input variable InputkT).
B. Energy H (variable Energy) versus kBT .
C. Scatter plot (function scatter) of H versus | ~M|.
Each set of three plots (e.g., 2A, 2B, and 2C) should be accompanied by one sentence that
identifies, and suggests an explanation for, one difference between one of the plots from that
set (e.g., 2B) and the same plot from simulations of the other two systems (e.g., 1B and 3B).
You may choose to identify and explain any difference that you find interesting, striking, or
important, but sentences will be awarded more marks if they are clear, if they identify an
interesting and original/nontrivial difference, and if they demonstrate understanding.
I recommend automating this task. For example, you could first run a long equilibration, as
follows:
>> [t E Mmag Mdir f S] = Ising2dMMC(0.3, 10, LongEquil, 0., 0., 100, ‘new’);
and then run the remaining simulations automatically as follows:
MATE60002/70002: TSM: Paul Tangney 5
>> for ktemp = 0.3:-0.01:0.05
[t E Mmag Mdir f S] = Ising2dMMC(ktemp, 10, ShortEquil,0., 0., 0, ‘restart’);
[t E Mmag Mdir f S] = Ising2dMMC(ktemp, 10, Production,0., 0., 0, ‘restart’);
end
where LongEquil, ShortEquil, and Production are the numbers of sweeps in the simulations.
They are chosen large enough to equilibrate initially, to equilibrate starting from a microstate
at the previous temperature, and to calculate averages, respectively. I have deliberately chosen
an uninteresting temperature range for this example because I want you to choose your own.
Task 2
Perform a series of simulations to investigate the dependence of any quantity or property
on any other. Quantities that you might find interesting to investigate include, but are not
limited to,
- Energy, H
- kBT
- ~M (magnitude or direction)
- ~h (magnitude or direction)
- Lattice size (LatticeSize) and whether it is even or odd (this is relevant when J = −1).
- Number of steps and sweeps taken to reach equilibrium.
- Specific heat capacity (per spin).
- | ~M| in the initial microstate (use internal variables probSpinRight and probSpinUp).
- log W (M), where W (M) is the number of microstates with magnetization of magnitude
| ~M| = M (see function MagnetizationEntropy for more explanation).
- The average sizes of ordered domains.
- The probability of an ordered domain growing as a function of its size.
You can choose J = 1 or J = −1, and you can choose any values for the other parameters
of the simulation, but choose them wisely! For example, there is no point in studying the
dependence of | ~M| on temperature for a range of values of kBT that is much smaller than
the value of |J| because if you simulate at very low temperatures, the magnetization changes
very little and you are unlikely to be able to detect small changes unless you run very long
simulations. Play with the program before deciding on something interesting to investigate
and on the range of values of each quantity that you will study.
You do not need to investigate any of the quantities listed. You could choose a different
quantity, and change the code to calculate it.
6 MATE60002/70002: TSM: Paul Tangney
DO NOT simply repeat part of Task 1. Do something different.
Submit the following:
1. A list of all technical details of your simulations, such as the parameters used and any
significant changes you have made to the code.
2. One plot.
3. One sentence describing the what you have investigated.
4. One sentence identifying, and suggesting an explanation for, an interesting feature of
the plot.
You may also provide one or two images of the spin structure, like the one on Page 3, if they
help to illustrate the point made in your sentence(s), but only do so if they serve an important
purpose.