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

Question

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!

Answer 1

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

$<E>=-\frac{J}{N}\sum_{<i,j>}S_iS_j$

$<E>^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 $<i,j>$ 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!).