Thermodynamics of the 2D Ising model

Question

I want to study the thermodynamics of the 2D Nearest Neighbour Ising model (calculate the average energy, susceptibility, etc.). I have the Hamiltonian

$$\mathcal{H} = J\sum_{\langle i j \rangle} s_i s_j + h\sum_i s_i$$

where $s_i \in \{-1, 1\}$. I put these spins $s_i$ on a lattice of size $n\times n$. To calculate the partition function, I need to sum over all possible spin configurations (permutations) of the spins $s_i$.

However, I am not sure how I should go about summing over all possible permutations in Mathematica and implementing the 'nearest neighbour' interaction for a 2D grid.

Answer 1

The energy can be calculated exactly as written in the formula using Sum and correct indices:

energy[s_] := J*Sum[(s[[i + 1, j]] + s[[i, j + 1]] + 
  s[[i - 1, j]] + s[[i, j - 1]])*s[[i,j]], 
  {i, 2, Length[s] - 1}, {j, 2, Length[s] - 1}] - h*Total[s, 2]

To calculate all possible spin configurations, you can use Tuples:

n = 3;
configs = Partition[#, n] & /@ Tuples[{0, 1}, {n^2}];

Note that the number of configurations is $2^{n^2}$, which grows extremely fast. You will not be able to enumerate them for $n\gtrsim 5$.

You can now calculate the energy of all configurations:

energy /@ configs
(* {0, -h, -h, -2 h, -h, -2 h, -2 h, ..., -8 h + 4 J, -9 h + 4 J} *)

---

Instead of pre-generating all possible configurations, you can generate them on-fly, as suggested by @Roman, for example with @Szabolcs' answer. This will not throw a memory allocation failure, but it will still take a lot of time to calculate:

Sum[energy[Partition[IntegerDigits[k, 2, n^2], n]], {k, 0, 2^(n^2) - 1}]