Simulation of the 2D ising ferromagnet, which has Hamiltonian (from 10.40 Set 9)

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

where $s_i \in {-1, 1}$ is the spin on lattice site $i$, the first summation runs over nearest neighbor pairs $\langle i, j \rangle$, $J$ is the coupling constant between neighboring spins, and $h$ is the external magnetic field.