Ising 2d Using Montecarlo Metropolis (Markov chain method) I'm doing as a personal training the 2d Square lattice Ising model. I decided to go with metropolis Monte Carlo method using Markov chain. I'm not into this methods, but I'm just using them as a tool (and maybe that's why I'm failing). 
Im following this paper (https://arxiv.org/abs/0803.0217).
I understand the way I let my system evolve (doing the successive Monte Carlo steps) but when it comes to get some physical quantities, I just realize some algorithms sum all the spins each Monte Carlo step. 
My intuition tells me I should just finish all the Monte Carlo steps, and in the end of the program, calculate the energy per spin, the sum of all spins, the $C_v$, etc. Why I'm not seeing? Is my intuition wrong?
Thanks, and sorry if this is very poor, just tell me and I will try to elaborate more. I'm doing this with C (As I said, as a training, I'm proficient with fortran)
 A: I'm not sure I understand your question, but I think you're asking how to sample physical quantities, and this is what I'll try to give an answer to.
The general idea is that MCMC simulations give you a collection (an ensemble) of equilibrium configurations. In principle you can average over the instantaneous values of the physical quantities of interest (energy, magnetization, etc.) computed for each of these configurations. In practice, configurations that are a "few" steps away are correlated, and thus it doesn't make a lot of sense to output the observables every step. What people usually do is thus to compute these quantities every a fixed number of steps (1k, 10k, you name it) and take the average at the end of the simulation. The optimal lag value is somewhat dependent on the details of the simulations (system studied, algorithm employed, physical conditions, etc).
On the practical side, what you usually do is the following (at least for simple observables that can be straightforwardly computed while the simulation is running):


*

*You set your simulation to print the desired observables every N time steps

*You run the simulation

*You do a post-run analysis on the data you have acquired


One of the advantages of the post-processing is that you can, for instance, easily throw away the part of the signal relative to the equilibration process.
During the post-processing, you may want to also estimate the error associated to the averages you have computed. A common strategy is (again, for simple observables)


*

*Calculate the autocorrelation of the signal you want to average

*Calculate how many steps are required to decorrelate it (this is usually done by fitting the autocorrelation to an exponential and using its characteristic time $\tau$ as the decorrelation time)

*Take your average, compute the associated variance and divide it by the number of uncorrelated samples you have ($N_t/\tau$, where $N_t$ is the total number of steps). The square root of that number is the error associated to the average.

