[gmx-users] g_energy averages

[Erik Lindahl]

David L. Bostick wrote:
> .. So, if I want to find the standard error in an average value, I just
> take the fluctuation and divide it by sqrt(nsteps)?
> 
> i.e.
> If I get from the energy file
> Energy                      Average       RMSD     Fluct.
> #Surf*SurfTen               147.805    2270.69    2263.84
> 
> and the run was over 250000 steps, the standard error is
> 
> 2263.84/sqrt(250000)
> 
> even if nstenergy was say 500 steps?
> 

Yes and no.

It is only correct if each point is an independent observation. If we
have a membrane where the area fluctuates on a scale of 10ns the area
fluctuations will be very small if we only simulate it for 10ps, even
if we save 100000 points during those 10ps. We can't say anything about
the 'physical' (in contrast to purely statistical) standard error 
though, since we don't even cover a single period of the motion.

To get an accurate estimate of the standard error you will have to take 
the autocorrelation time into account. Load the curve into xmgrace, 
calculate the autocorrelation of the fluctuations, and fit e.g. a single 
exponential A*exp(-t/tau) to it - 'tau' is your correlation time.

Assuming the total length of your simulation is T, we can then estimate 
the number of independent observations as T/tau, so the standard error
is


   <rmsd fluctiations>
s=-------------------
          T/tau

In the limit where the observations are completely uncorrelated the
correlation time is the stepsize, so in that case the denominator would
be the total number of steps.



The alternative is to split the trajectory into parts that you *know* 
are longer than the correlation time, calculate an average from each of
these, and then the standard error as the fluctuations of these averages
divided by the number of parts. This will only be an upper estimate of
the standard error though, and you will still need an estimate of the
correlation time.

Cheers,

Erik