7
$\begingroup$

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: 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!

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

1 Answer 1

9
$\begingroup$

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:

The heat capacity for the 2D-Ising model

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...

$\endgroup$

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.