Trying to compute the heat capacity for the 2D-Ising model

I'm trying compute the heat capacity $$C_v$$ out of my simulation for the 2D-Ising model which is given by $$C_v = \frac{\langle E^2 \rangle - \langle E \rangle^2}{T^2N^2}$$ ($$E$$: Energy, $$T$$: Temperature, $$N$$: Number of points). I implemented the Ising simulation using the Metropolis-method which seems to be the easiest way to simulate this:

Ising[n_, count_, B_, J_, T_] := Module[{},
e1 = 0; e2 = 0;
dE[a_, b_]= a*2*(B - b*J)/T;
Lat = 2*Table[Random[Integer], {n}, {n}] - 1;

flip =
(
L = #;
{i1, i2} = {RandomInteger[{1, n}], RandomInteger[{1, n}]};

If[i1 == n, down = 1, down = i1 + 1];
If[i2 == n, rechts = 1, rechts = i2 + 1];
If[i1 == 1, up = n, up = i1 - 1];
If[i2 == 1, links = n, links = i2 - 1];
Es = L[[down, i2]] + L[[up, i2]] + L[[i1, rechts]] + L[[i1, links]];
If[ dE[ L[[i1, i2]], Es ] < 0 ||
RandomReal[] < Exp[-dE[ L[[i1, i2]], Es ]],
L[[i1, i2]] = -L[[i1, i2]]; e1 = e1 + dE[ L[[i1, i2]], Es ];
e2 = e2 + dE[ L[[i1, i2]], Es ]*dE[ L[[i1, i2]], Es ], L;
e1 = e1; e2 = e2];
L
) &;

fliplist = NestList[flip, Lat, count];
]


This simulation is working for the plots of the lattice.

But as I try computing the terms $$\langle E \rangle$$ and $$\langle E \rangle^2$$ and using them to compute $$C_v$$ the result is not what I am expecting.

For large numbers of iterations (in this case equivalent to time units $$t$$) the mean value of the energy can be computed by:

$$\langle E \rangle = \frac{1}{count} \sum_{t=0}^{count} E_t$$

The sum over $$E_t$$ is what i thought i was doing with the e1 and e2 in the code. As i repeat the simulation for several temperatures and then use the values to calculate $$C_v$$ the result is not the heat capacity.

Do[steps = 50000;
Ising[64, steps, 0, -1, i]; Subscript[epts, i] = 1/steps e1;
Subscript[esq, i] = 1/steps e2,
{i, 1, 3, 0.1}];

Do[Subscript[c, i] =
1/(2*i*i) (Subscript[esq,
i] - (Subscript[epts, i])*(Subscript[epts, i])), {i, 1, 3, 0.1}];


The plot of those points looks like this:

This is not the expected plot for $$C_v$$. There should be a peak at around 2.3.

Could anybody help me to calculate the real $$C_v$$ because I can't find the mistake. I'm very thankful for any help!

• What are your axes? Around 2.3? So possibly you have missed this with too few of plot points? Nov 16, 2019 at 22:16

1 Answer

For this question, I can find some code of python version, Refs. here. And I rewrote this Ising simulation as Mathematica version, as follow,

1. Define a function used for generating a table (as configuration)

Initialstate[n_Integer] := 2*Table[RandomInteger[], {n}, {n}] - 1


where n as the number of points. And define a function used for getting any point's neighbor.

GetNeighbor[x_Integer, y_Integer, size_Integer] := {{If[x - 1 > 0, x - 1, size], y}, {x,
If[y - 1 > 0, y - 1, size]}, {If[x + 1 <= size, x + 1, 1], y}, {x,
If[y + 1 <= size, y + 1, 1]}} /; (0 < x <= size && 0 < y <= size)


2. Monte Carlo move using Metropolis-method

MCMove[config_, beta_] :=
Module[{a, b, s, nb, cost, size = Length@config,
configReturn = config}, Table[
a = RandomInteger[{1, size}];
b = RandomInteger[{1, size}];
s = config[[a, b]];
nb = config[[#[[1]], #[[2]]]] & /@ GetNeighbor[a, b, size] //
Total; cost = 2*s*nb;
If[cost < 0, s = s*(-1),
If[RandomReal[] < Exp[-cost*beta], s = s*(-1)]];
configReturn[[a, b]] = s, size*size]; Return@configReturn]


where config is random number table generate by Initialstate function.

3. Define some functions used to get the energy of the given configuration

Energy[config_] := Module[{energy, size = Length@config},
Total[Table[-config[[i,
j]]*(config[[#[[1]], #[[2]]]] & /@ GetNeighbor[i, j, size] //
Total), {i, 1, size}, {j, 1, size}], 2]/4
]

GenerateEnergyTable[config_, tStart_Real, tEnd_Real, nt_Integer,
eqSteps_Integer, mcSteps_Integer] :=
Module[{TList = Subdivide[tStart, tEnd, nt], config0 = config,
energyList, n = Length@config},
energyList =
Function[t, config0 = Nest[MCMove[#, 1/t] &, config0, eqSteps];
({Total@#, Total[#^2]} &@
Table[(config0 = MCMove[config0, 1/t]; Energy[config0]),
mcSteps])][#] & /@ TList;
Transpose[{TList, energyList}]
] /; (Depth[config] == 3 && (#[[1]] == #[[2]] &@Dimensions[config]))


About parameters are described below:

• tStart(tEnd) : temperature start (end) point
• nt : number of temperature points
• n : size of the lattice, N x N
• eqSteps : number MC sweeps for equillibration
• mcSteps : number of MC sweeps for calculation

In particular, the function GenerateEnergyTable can generate a table, which includes temperature, $$\sum E_i$$ and $$\sum E_i^2$$. Next we define a function used to calculate heat capacity $$C_v$$.

4. Define Heat capacity function

HeatCapacity[energyTable_, mcSteps_Integer,
n_Integer] := {#[[
1]], (#[[2]][[2]]/(mcSteps) - #[[2]][[1]]^2/(mcSteps^2))/(#[[
1]]^2*n^2)} & /@ energyTable


5. Calculation

energylist =
GenerateEnergyTable[Initialstate[16], 1.53, 3.28, 88, 500, 500];
ListPlot[HeatCapacity[energylist, 500, 16], PlotTheme -> "Scientific",
Filling -> Axis]


Runing above code, we can get one figure as below:

As you said, it has a peak at around 2.3. I hope that these help you.

Furthermore, maybe you can use some method to improve the speed of calculation, such as compiled function, parallel calculation, etc...