Computing autocorrelations of configurations in Monte Carlo simulations

Question

In the context of Monte Carlo simulations, I am trying to learn how I should ensure that the configurations of my system are not correlated for the chosen interval of measurements. I have found out that one way is to compute the autocorrelation function for a relevant observable.

Let us suppose we have a system of 1000 particles in a box. In case of a simple scalar observable $A$ that can be measured for each generated configuration (i.e. a function of particle positions), then the autocorrelation is, with $t$ denoting the Monte Carlo sweep step:

$$ C(t)=\frac{\langle A(t=0)A(t)\rangle - \langle A(t) \rangle}{\langle A^2 \rangle - \langle A\rangle^2} \propto \exp(-t/t_d) \tag{1} $$

where $A(t=0)$ is the initial value of our observable at the start of the simulation, $A(t)$ is the value corresponding to the configuration generated at MC sweep step $t,$ $\langle A^2\rangle$ is the average of the observable squared thus far, i.e. if we are at step $t$ then $\langle A^2\rangle = \frac{1}{t}\sum_{i=0}^{t} A^2_i,$ and finally if the exponential decay behaviour holds, the found time scale $t_d$ would tell us how many MC sweep steps apart should our subsequent measurements be to ensure they are not correlated.

1. Have I made any mistakes with the eq. $(1)$ and the description of the involved terms?

My last question: usually if $A$ is an order parameter of the system of particles, then once equilibrium is reached $A$ becomes a constant with small fluctuations (e.g. $A$ can be the energy of the system), and in such case, to ensure our measurements are still uncorrelated, should we compute the autocorrelation for the particle positions? That is, it would tell us how many steps it would be needed until the positions of the particles are sufficiently different from the previous configuration.

2. But if we have $N$ particles, how do we compute a position autocorrelation function for the entire system? That is, given we have $N$ particles and each have a $3$-component position vector, how will $A$ be defined in eq. $(1)$?

Answer 1

To 1.: This is not quite right, your definitions seem to be a bit messy. As you noted correctly, the averages are taken over time $t$. Consider $A(t)$ some signal of length $T$, the mean of the signal is

$$\langle A(t) \rangle = \frac{1}{T} \sum_{t=0}^{T-1} A(t) $$

Since we take the averages over time $t$ sometimes we use $\langle A \rangle$ and $\langle A(t) \rangle$ interchangeably. Note that with these definitions the denominator in your $C(t)$ is actually equal to the variance of $A(t)$ over time:

$$ \sigma_A^2 = \langle \left( A(t) - \langle A \rangle \right)^2\rangle = \langle A^2 \rangle - \langle A \rangle^2$$

So let's get back to the acf. The autocorrelation is not a function of time $t$, but rather a function of delay/time difference $\Delta t$. The proper definition would be

$$ C(\Delta t) = \frac{\langle \left( A(t + \Delta t) - \langle A \rangle \right) \left( A(t) - \langle A \rangle \right) \rangle}{\langle A^2 \rangle - \langle A \rangle^2} = \frac{\langle A(t + \Delta t) A(t)\rangle - \langle A \rangle^2}{\langle A^2 \rangle - \langle A \rangle^2}$$

So why is that? $C$ is not supposed to give us information about some point in time $t$, but what we want to know is how much an observation at a given point in time is correlated with the observation $\Delta t$ later. So the numerator of $C$ computes the covariance of the signal $A$ with itself $\Delta t$ later. Since the covariance is going to be numerically scaled with the variance, we'd much rather express autocorrelation in units of the signals variance. This is why we put the signals variance in the denominator. This way the acf for every signal is going to take values between -1 and 1, independet of the variance. Only this way we can reasonably talk about the time scale $t_d$.

To 2.: I might be wrong about this, but I would say it depends a bit on what you are trying to sample. Let's say you want to get i.i.d. samples of the particle positions. I would argue that it is sufficient to adjust $t_d$ using the systems energy, because the energy is going to be a deterministic function of the positions (assuming your phase space is only the particle positions and there is no other time-dependet/random influences on the energy). If the positions were correlated, this should manifest itself in a correlation of energy, at least in a physically realistic system.

A counter-example to this reasoning might be a pseudo-random numbers generator, where a deterministic series of numbers is processed by deterministic functions until they appear random. However, given how tricky it is to design such pseudo-RNGs, I should think that you can assume the reasoning above to work.

If you want to be slightly more serious about it anyway, you could compute the acf for each dimension of the state space independently and take the longest time-scale of all dimensions. For $N$ particles in 3D you have a $3N$ dimensional signal and would have $3N$ different time-scales. Computing this will have a huge runtime, so you might want to check here for tipps on how to do it efficiently using FFT.

If you do it this way you would still be neglecting cross-correlations between the dimension. You'd have to compute an object like

$$ C_{ij}(\Delta t) = \frac{\langle x_i(t+\Delta t) x_j(t)\rangle - \langle x_i \rangle \langle x_j \rangle}{\sigma_{x_i} \sigma_{x_j}}$$

and then choose $t_d$ such that $C_{ij}(\Delta t) < e^{-1}$ for all $i,j$ and all $\Delta t > t_d$. At this point I'm not sure if this is done in practice, as $ C_{ij}(\Delta t)$ is quite a massive thing to compute.