Fortran Analysis Dynamics and Errors
The aim of the program corfun is to illustrate the direct method, and the FFT method,
for calculating time correlation functions.
The program is self contained: it generates the time dependent data itself,
using a generalized Langevin equation,
for which the time correlation function is known.
The default parameters produce a damped, oscillatory, correlation function,
but these can be adjusted to give monotonic decay,
or to make the oscillations more prominent.
If the origin_interval parameter is left at its default value of 1,
then the direct and FFT methods should agree with each other to within numerical precision.
The efficiency of the direct method may be improved,
by selecting origins less frequently,
and in this case the results obtained by the two methods may differ a little.
Sample results using default program parameters are shown here. The direct method (orange line) is overlayed perfectly by the FFT result (dashed green line), as expected. The exactly known function is a blue line. There are very small discrepancies with the results of the simulation, due to the finite length of the latter.
The program diffusion reads in a sequence of configurations and calculates
the velocity autocorrelation function diffusion.out.
It is instructive to plot all three curves vs time.
The input trajectory is handled in a crude way,
by reading in successive snapshots with filenames cnf.000, cnf.001, etc.
These might be produced by a molecular dynamics program,
at the end of each block,
choosing to make the blocks fairly small (perhaps 10 steps).
As written, the program will only handle up to cnf.999.
Obviously, in a practical application,
a proper trajectory file would be used instead of these separate files.
It is up to the user to provide the time interval between successive configurations.
This will typically be a small multiple of the timestep used in the original simulation.
This value delta is only used to calculate the time, in the first column of
the output file.
A default value of 0.05 is provided as a place-holder, but
the user really should specify a physically meaningful value;
forgetting to do so could cause confusion when one attempts
to quantify the results.
To make it easier to test this program,
we have also supplied a self-contained program diffusion_test,
which generates an appropriate trajectory by numerically solving
the simple Langevin equation for N non-interacting atoms (N=250 by default) in 3D.
For this model, one specifies the temperature and friction coefficient,
which dictates the rate of exponential decay of the vacf,
and hence the diffusion coefficient.
The exact results for the above three functions are written out to diffusion_exact.out
for easy comparison with diffusion.out.
Here are some typical results using default program parameters throughout.
The results of the Langevin simulation are shown as points, plotting only every 10th value for clarity,
while the exact results are indicated as lines.
For the default program parameters, the diffusion coefficient is D=1.
The program error_calc is a self-contained illustration of the effects of
correlations on the estimation of errors for a time series.
We produce the series using a generalized Langevin equation,
in the same manner as for the correlation function program (see above).
Since the correlation time of the GLE is exactly known,
we can predict the effects, and compare with the empirical estimates
obtained by different methods.
The program contains extensive comments to explain what is being calculated at each stage.
The aim of fft3dwrap is to illustrate the way a standard Fast Fourier Transform
library routine is wrapped in a user program.
We numerically transform a 3D Gaussian function,
and compare with the analytically, exactly, known result.
User input defines the number of grid points and the box size;
sensible defaults are provided.
The library that we use for this example is FFTW.