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];
 ) &;

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: C_v vs Temperature

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!


1 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>^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!).

  • $\begingroup$ 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. $\endgroup$ Feb 5, 2019 at 1:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.