(* You will need the computer to generate random numbers, which are not truely random but are rather computed based on some seeded value. If you set the seed to be the same thing every time, you'll always produce the same "random" numbers, which can sometimes be useful when debugging errors. *)
SeedRandom[1234]

Get comfortable using Mathematica to generate and display configurations of a grid.

First we make a grid and initialize it with either all ones or with random ones and minus ones.

(* sl is the side length *)
sl = 20;
UpConfiguration = Table[1, {i, 1, sl}, {j, 1, sl}];
RandomConfiguration = Table[RandomChoice[{-1, 1}], {i, 1, sl}, {j, 1, sl}];

(* This grid could be viewed as a matrix *)
MatrixForm[RandomConfiguration]

It’s much prettier if we visualize this with shaded boxes, which is pretty trivial with Mathematica:

ArrayPlot[
Transpose[UpConfiguration],
Mesh->All,
ColorRules->{1->Black,-1->White},
ImageSize->400,
(* This DataReversed line is there so the spin at position {i, j} will be plotted i units to the right and j units up. *)
DataReversed->{True, False} 
]

ArrayPlot[
Transpose[RandomConfiguration],
Mesh->All,
ColorRules->{1->Black,-1->White},
ImageSize->400,
(* This DataReversed line is there so the spin at position {i, j} will be plotted i units to the right and j units up. *)
DataReversed->{True, False}
]

Rather than rewrite the same piece of code over and over again, it’s useful to just define a function:

PlotConfig[config_]:=ArrayPlot[Transpose[config],Mesh->All,ColorRules->{1->Black,-1->White},ImageSize->400, DataReversed->{True, False}]
(* So the random configuration can now be plotted with: *)
PlotConfig[RandomConfiguration]

Compute the energy

We are going to assume that the sides of the box wrap around so that the spins along the left column are neighbors with the spins along the right column. Similarly, the top and bottom spins touch. Put another way, our spins actually live on the surface of a sphere and we’re choosing to stretch them out onto a plane when we plot them. The (slightly) tricky thing we have to figure out is which spins are neighbors with which other spins since the energy only adds contributions from the closest neighbors.

Neighbors[spinlocation_, sl_]:=Module[{spinx, spiny},
{spinx, spiny} = spinlocation;

(* Naively you might think that the following will give you the neighboring spin locations, but these neighbors could go outside of the grid, which runs only from 1 to sl. Adjust the next four lines of code to make sure that the edges are handled appropriately. *)
rightneighbor = {spinx + 1, spiny};
leftneighbor = {spinx - 1, spiny};
upneighbor = {spinx, spiny + 1};
downneighbor = {spinx, spiny - 1};

(* This final line is what the function will return, the locations of the four neighboring spins. *)
{rightneighbor, leftneighbor, upneighbor, downneighbor}
]

(* Test your Neighbors function to make sure it gives you what you'd expect *)
Neighbors[{1, 1}, sl]
Neighbors[{20, 1}, sl]
Neighbors[{4,5}, sl]

(* Compute the energy coming from the spin at position spinlocation's interactions with its four neighbors *)
SingleSpinEnergy[spinlocation_, sl_, config_]:= (* Write a function to compute the energy coming from interactions between the spin at spinlocation = {spinlocationx, spinlocationy} and its four neighbors. You will probably want to make use of your Neighbors function! *)

(* Test out your energy function by calculating the energy from several spins and looking at the grid *)
SingleSpinEnergy[{4, 5}, 20, RandomConfiguration]
SingleSpinEnergy[{8, 8}, 20, RandomConfiguration]
SingleSpinEnergy[{20, 18}, 20, RandomConfiguration]

PlotConfig[RandomConfiguration]

(* The total energy is the sum of all the single spin energies (divided by two to prevent double counting) *)
Energy[config_, sl_] := -(1/2) Total[Table[SingleSpinEnergy[{i, j}, sl, config], {i, 1, sl}, {j, 1, sl}],2];

Here is a function that accepts a configuration, tries to flip a spin at a random position, computes the change in energy, then randomly returns either the new trial configuration or the original.

The function is called SlowMCStep because I’m computing the change in energy in a stupid way, which makes the function really slow.

SlowMCStep[config_, sl_, T_]:=Module[{spinflipx, spinflipy, EnergyBeforeSpinFlip, EnergyAfterSpinFlip, spinFlippedConfig, ΔE, AcceptanceProbability},
spinflipx = RandomInteger[{1, sl}];
spinflipy = RandomInteger[{1, sl}];

spinFlippedConfig = config;
spinFlippedConfig[[spinflipx, spinflipy]] *= -1;

EnergyBeforeSpinFlip = Energy[config, sl];
EnergyAfterSpinFlip = Energy[spinFlippedConfig, sl];
ΔE = EnergyAfterSpinFlip - EnergyBeforeSpinFlip;

AcceptanceProbability = Min[1, Exp[-ΔE / T]];
(* Return either the flipped configuration or the original configuration randomly with probability given by AcceptanceProbability *)
If[RandomReal[]≤AcceptanceProbability, spinFlippedConfig, config]
]

Configuration= RandomConfiguration;
(* Take 1000 Monte Carlo steps, storing each configuration as a list called MCTrajectory. *)
MCTrajectory = Table[Configuration = SlowMCStep[Configuration, sl, 1], {trials, 1, 1000}];
(* Create an animation of those configurations *)
ListAnimate[
Table[
PlotConfig[MCTrajectory[[i]]],
{i, 1, 1000, 10}
]
]

The smarter way to compute the difference in energy would use the function SingleSpinEnergy rather than Energy. Write a new function to take a much faster Monte Carlo step. Use this faster Monte Carlo algorithm for all of the calculations in the problem set.

FastMCStep[config_, sl_, T_]:=Module[{spinflipx, spinflipy, EnergyBeforeSpinFlip, EnergyAfterSpinFlip, spinFlippedConfig, ΔE, AcceptanceProbability},
Write this!
]

Configuration= RandomConfiguration;
(* Take 1000 Monte Carlo steps, storing each configuration as a list called MCTrajectory. *)
MCTrajectory = Table[Configuration = FastMCStep[Configuration, sl, 1], {trials, 1, 10000}];
(* Create an animation of those configurations *)
ListAnimate[
Table[
PlotConfig[MCTrajectory[[i]]],
{i, 1, 10000, 10} (* I am only viewing every 10 frames *) 
]
]

Note that the magnetization is easy to compute:

(* Here is one way *)
M = Total[Total[RandomConfiguration]]

(* Here is another way *)
M = Total[RandomConfiguration, 2]

(* Here is a third way *)
M = Total[Flatten[RandomConfiguration]]

(* This is not a way *)
NotM = Total[RandomConfiguration]

MTimeSeries = Table[
(* Loop over the MC timesteps with variable i*)
(* Record the MC timestep and the Magnetization for that timestep *)
{i, Total[MCTrajectory[[i]], 2] / (sl * sl)},
 {i, 1, Length[MCTrajectory]}
];

ListPlot[
MTimeSeries,
Frame->True,
FrameLabel->{Style["MC Steps",FontSize->18], Style["Magnetization per site",FontSize->18]},
ImageSize->300,
AspectRatio->1
]