Python Analysis Dynamics and Errors

Correlation function program

The aim of the program corfun.py 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.

In the current implementation, the direct method makes little or no use of NumPy's efficient array manipulation functions, and so is very slow. For this reason, the default run length is much shorter than for the Fortran example. An alternative, much faster, direct method is also provided at the end of the program. This uses the NumPy correlate library function. The selected mode 'full' corresponds to aligning the data array v with itself, at all possible offsets that result in some overlap of values, before computing the products (for each offset) and summing. In this case, the resulting array is symmetric in time (offset), and only the values from the mid-point onwards are required; these must be normalized in the standard way and identical results to the slow direct method (with origin_interval=1) and FFT method, are obtained.

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 library function result is not shown, as it is indistinguishable from the first two. The exactly known function is a blue line. There are clear discrepancies with the results of the simulation, as expected, due to the rather short duration of the latter. The very slow direct method could easily be edited out of the program, and the run length dramatically increased, to obtain more accurate results.

Diffusion program

The program diffusion.py reads in a sequence of configurations and calculates the velocity autocorrelation function $\langle\mathbf{v}(0)\cdot\mathbf{v}(t)\rangle$ (vacf), the mean square displacement $\langle |\Delta\mathbf{r}(t)|^2\rangle$ (msd), where $\Delta\mathbf{r}(t)=\mathbf{r}(t)-\mathbf{r}(0)$, and the cross-correlation between velocity and displacement $\langle\mathbf{v}(0)\cdot\Delta\mathbf{r}(t)\rangle$. Any of these may be used to estimate the diffusion coefficient, as described in the text. The output appears in 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.py, 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 results are very similar to those obtained from the analogous Fortran example.

Error calculation

The program error_calc.py 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.

FFT program

The aim of fft3dwrap.py 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 the built-in NumPy one.