ISING2D with Mathematica. Searching a correct way to compute the heat capacity (mean values over several iterations)

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];
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!

Unfortunately I am not a Mathematica whiz. and in general do not encourage the use of Mathematica for extensive numerical simulations like Ising model or the likes. So I will post some answers pertaining to C++ and the general equations.

First I am going to assume, that your $$E_t$$ is actually the spin-exchange energy $$-JS_iS_j$$ where $$S_i$$ are the spins and $$J$$ the interaction parameter. Heat capacity $$C$$ is then computed, as you'Ve stated via

$$C=\frac{\langle E^2\rangle - \langle E\rangle^2}{k_bT^2}$$

where the following holds

$$=-\frac{J}{N}\sum_{}S_iS_j$$

$$^2 = \frac{1}{N}\sum_{i=1}^N\left(-J\sum_{j\in n.N.} S_iS_j\right)^2$$

where in the first equation $$$$ indicates summation over nearest-neighbours, what has been split in the second one into a sum over $$i$$ and for each $$i$$ a sum over its nearest neighbours. $$N$$ is the total number of sites ($$=n^2$$ for a square grid of $$n$$ cells in each direction).

Also I'm pretty sure you are not supposed to average over multiple iterations when you are computing the observables (also:Label your axes!).

• Okay i got it now. What it did for me was the suggestion not to use Mathematica. I put everything as I thought I did it into Python and it worked. – PaladinDanse Feb 5 at 1:30