I am have written a Metropolis-Hastings algorithm and am currently trying to compare it to the analytical results for the 2D Ising model. The free energy seems reasonable but the heat capacity I'm getting is way off.
I am using the equation $C_{v}=\frac{\beta}{T}[\left<E^{2}\right>-\left<E\right>^{2}]$. Using $\beta=1$ would imply that $T$ is on the order of $10^{23}$ so the prefactor becomes really small. My understanding was that the variance in the energy will be humongous at such high temperatures but my variance is not that large. I think it is because I am using a $16\times16$ lattice and the variance on a finite system (I'm guessing) is bounded. But the computation time is already really slow (5 minutes for $10^{6}$ iterations of Monte Carlo).
I have tried using $C_{v}=-\beta^{2}\frac{\partial^{2} f}{\partial \beta^{2}}$ with a numerical approximation for the second derivative. But my free energy values are also numerically calculated so I sometimes get nonsensical results like negative heat capacity. I was wondering if there are other ways of estimating the heat capacity. Or if I'm stuck with using the variance of the energy, how large a lattice do I need to have?