Pieces of iron exhibit magnetism. Magnetism can be explained at the molecular level by considering that each iron atom has a "spin" that can be "up" or "down". In a fully magnetized piece of iron, all the spins point in the same direction. Most iron is not fully magnetized, having many atoms with up and down spin. But even in this case, the iron atoms will typically be organized, so that there are large areas where all the atoms have up spin, and others where the spin is down. If the iron is heated up, and allowed to cool down, the spins will again be scattered. But gradually, over time, the arrangement of spins into common territories will return.
The idea of this project is to set up a simple two dimensional model of a piece of iron. The model should change over time, and its behavior should depend on the temperature also.
As our first task, let's just try do make a model of the material, and draw a picture of it. Let's say we have a block of material that is 1 unit wide and 1 unit tall, containing 20 atoms arranged in a simple 4 by 5 grid. We will keep track of the atoms by row I and column J. Each atom will have a spin of +1 or -1, and we can store the spin of a particular atom in spin(i,j). To start, we'll assume the spins are set randomly.
The only random number generator we have, rand, returns random numbers between 0 and 1. How can we use this to make random integers that are either -1 or 1? A useful tool for this process is the ROUND command, which rounds a number to the nearest integer.
Exercise A: Figure out how to use the ROUND command and some arithmetic to correct the spins in the following MATLAB M file ising_init.m.
function [ x, y, spin ] = ising_init ( m, n )
dx = 1 / ( m - 1 );
dy = 1 / ( n - 1 );
[ x, y ] = meshgrid ( 0:dx:1, 0:dy:1 );
spin = rand ( m, n );
Your work is correct if your final spins are random values of +1 and -1.
Now we'd like to draw a picture of our problem. Let's imagine that a negative spin is red, and a positive spin is blue. Then our picture should just be a bunch of red and blue dots at the appropriate places. We can almost do this. The command
plot ( x, y, 'ro' )
will plot red dots at every point. We're close, but now we need to
figure out how to plot red dots at some places, and blue dots elsewhere.
The command will look something like this:
plot ( x(down), y(down), 'ro', x(up), y(up), 'bo' )
where somehow down and up are selecting the points
we want.
The tool we will need is the MATLAB FIND command, which can return a list of the indices in an array for which a certain condition is true. To ask for the indices of spin where the value is equal to 1, we write:
up = find ( spin == 1 )
(Notice that we have to use two equal signs when we describe the value of spin. One equal sign is actually the assignment operator, which takes the value of the right hand side and stores it in the variable on the left hand side.)
Exercise B: use the FIND command to define the two lists up and down, and create a MATLAB M file ising_plot:
function ising_plot ( x, y, spin )
up = find ( spin == 1 )
down = ?
plot ( ? )
which makes a plot of the current spin array.
Now it's time to look at the physics in the model. The spin of an atom is affected by the spins of nearby atoms; in fact, an iron atom "prefers" to have the same spin as the average of its neighbors; that is, the same-spin-state has a lower energy. We need to be able to compute the energy associated with the spin arrangements in our simple model.
We start by assuming that only the 4 immediate neighboring atoms can affect the spin of an atom. We imagine the atom being connected to its neighbors by little lines that we call call pair_bonds. In our model, the total energy of the material will be the sum of the energy associated with each of these pair-bonds. The energy of a single pair bond will be
E = - S1 * S2
Example 1: What is the total energy of the following material, made of 20 atoms:
-1 +1 +1 +1 +1
+1 -1 -1 +1 +1
+1 +1 +1 -1 -1
-1 -1 +1 -1 -1
Answer: The total energy is: ___________________
Example 2: Now suppose that the material "wraps around". Every atom now has four neighbors. For the atoms in the last row, for instance, the "southern" neighbors are the ones in the first row, and so on. In the diagram, the wrap-around neighbors have been added, in parentheses. Be careful not to count the wrap-around bonds twice! What is the total energy?
(-1)(-1)(+1)(-1)(-1)
(+1) -1 +1 +1 +1 +1 (-1)
(+1) +1 -1 -1 +1 +1 (+1)
(-1) +1 +1 +1 -1 -1 (+1)
(-1) -1 -1 +1 -1 -1 (-1)
(-1)(+1)(+1)(+1)(+1)
Answer: The total energy is: ___________________
The MATLAB M file ising_neighbors.m returns the values of the spins to the north, south, east and west of the atom in position (I,J):
function [ north, south, east, west ] = ising_neighbors ( spin, i, j )
For the following exercises, you will need to
download a copy of ISING_NEIGHBORS.M.
Exercise C: Write a MATLAB M file ising_energy.m that can compute the total energy of an array of spins. It should look something like this:
function e = ising_energy ( spin )
[ m, n ] = size ( spin );
e = 0;
for i = 1 : m
for j = 1 : n
e = e + ???
end
end
Try out your code using the data for Example 2 and see if you get the
answer you computed by hand.
Now we need to look at the process of random fluctuations. Let's consider one atom at random. It has two "choices", to be up or down. Each choice affects the energy of the system. We really only need to look at the four neighbor atoms, and compute the local energy. So here is the energy associated with the four local pair bonds if atom (I,J) is UP:
up_energy = - +1 * ( spin(i+1,j) + spin(i-1,j) + spin(i,j-1) + spin(i,j+1) )
and, of course, the down energy is:
down_energy = - -1 * ( spin(i+1,j) + spin(i-1,j) + spin(i,j-1) + spin(i,j+1) )
Example 3: Go back to Example 2. If we look at atom (2,3),
up_energy = - +1 * ( 1 + 1 - 1 + 1 ) = 3
and, of course, the down energy is:
down_energy = - -1 * ( 1 + 1 - 1 + 1 ) = -3
In the same way that water flows downhill, physical systems tend to seek the state of lowest energy. That means that in Example 2, atom (2,3) would "want" to be down, because that reduces the total energy of the system.
Let's think about adding fluctuations to our model. Our idea will be that, at random, a single atom is picked in the array. The up and down energies for the atom are compared, and the spin of the atom is reset to whichever value has the lower energy. We will consider repeating this process until the model settles down, that is, until all the spins are equal.
So if we write this as a MATLAB script file, ising_downhill.m, we have:
function new_spin = ising_downhill ( spin, steps )
for k = 1 : steps
i = ?
j = ?
[ north, south, east, west ] = ising_neighbor ( spin, i, j );
e_up = ?
e_down = ?
if ( ? ) then
spin(i,j) = ?
else
end
end
new_spin = spin;
Exercise D: Fill in the blanks in the ising_downhill.m code. Start with a 4 by 5 grid of cells, and have ising_init.m initialize the spins. Then run ising_downhill.m until the spins "settle down". You should probably reach a point when all the spins are equal. Try this 5 times, to see the variation in the number of steps required. Repeat on an 8 by 10 grid. Does it take about the same number of steps, twice as long, or three times as long?
Now what we saw in the downhill version of our model is that the material ends up with all spins pointing the same way. But real life materials don't behave this way, in part because of the effects of heat. We can think of heat as a disruptive energy that jostles the material and keeps it from settling down. In order to see the effects of heat, we must turn this idea into a mathematical form that can become part of our program.
The Boltzmann energy distribution says, in part, that if a system is at a nonzero temperature T, and there are two states to choose from with energies E1 and E2, then the probability P1 that the system will move to state 1 instead of state 2 is
P1 = R1
R1 + R2
where
R1 = e-E1/T
with a similar formula for R2.
Briefly, if the temperature is zero, then the system will always move to the state of lower energy; at very high temperatures, the system will choose either state with equal probability. In between, the system prefers the state of lower energy, but a rising temperature increases the chance it will take a step in the "contrary" direction.
Now before we proceed, the question is, does adding the effect of temperature merely delay the time at which the system settles down to a single spin, or does something interesting happen? This is our final question for this project (but there's lots more to learn about the Ising model!)
Let us see if we can program the effects of temperature properly. We're going to write a new MATLAB script file calling ising_hot.m. It will be similar to the downhill model, except that we pass it the temperature T, and once we compute the up and down energies, we have to compute the relative probabilities, and instead of always choosing the lower energy spin, we need to "flip a coin" somehow.
function new_spin = ising_hot ( spin, steps, t )
for k = 1 : steps
i = ?
j = ?
[ north, south, east, west ] = ising_neighbor ( spin, i, j );
e_up = ?
e_down = ?
p_up = e^( - e_up / t )
p_down = e^( - e_down / t )
r_up = p_up / ( p_up + p_down )
coin_toss = rand ( 1, 1 )
if ( coin_toss < r_up ) then
spin(i,j) = +1
else
spin(i,j) = -1
end
end
new_spin = spin;
Notice that this code will not work properly if T = 0. We could fix that for now, if we need results for zero temperature, we should just use the downhill code.
Exercise E: Get the ising_hot code working. Experiment with the model, using a 10 by 10 grid and a "low" temperature of T=0.5. When you're comfortable with the program, try to find a higher value of T at which roughly half the atoms are up and half down, for quite a few steps. Now look at some temperatures in between the cold and hot states. What you should hope to see is that for the "warm" states, the material will group itself into up and down areas. These areas will expand, contract, join and split. But unlike the material in the downhill model, the material will keep both up and down spins.
If you get some nice results, you might try making an animation. Another useful tool might be a graphs of the total energy over time, or of the percent of atoms with "up" spin.
Back to the Mini Projects page.
Last revised on 08 June 2001.