Self-diffusion coefficient

vacf.py calculates self-diffusion coefficient from particle velocities. The self-diffusion coefficient, \(D\), is computed from an exponential fit to the running integral of velocity auto-correlation function (VACF) using the following Green-Kubo relation

\[D = \frac{1}{3} \int_0^\infty \left\langle \boldsymbol{v}_i(t) \cdot \boldsymbol{v}_i(t_0) \right\rangle dt\]

where \(\boldsymbol{v}_i(t)\) denotes the velocity of particle \(i\) at any specific time \(t\).

Usage:

python vacf.py -h

Tip

An example data file, velocity.data, is available in the example directory. The required values to run the example can be found in the md.param file.

vacf.acf(velocities, max_lag)

Computes the velocity autocorrelation function (VACF) for a set of particles, which is a measure of how the velocity of a particle at a given time is correlated with its velocity at a later time.

Parameters:

  • velocities (np.ndarray) – A 3D NumPy array of shape (N, 3, M) where N is the number of particles, 3 represents the x, y, z dimensions, and M is the number of time steps.

  • max_lag (float) – Portion of the total simulation length to be used as maximum lag time.

Returns:

vacf – A 1D NumPy array representing the average VACF over all particles.

Return type:

np.ndarray

Notes

This implementation leverages the Wiener-Khinchin theorem, which states that the autocorrelation of a signal is the inverse Fourier transform of its power spectral density. This allows for a fast calculation of the autocorrelation. The maximum lag time is set to 30% of the total number of steps, for better statistics, as correlations often decay to zero long before this point.

vacf.diffusion(vacf, time, timestep)

Calculates the diffusion coefficient from the velocity autocorrelation function (VACF).

Parameters:

  • vacf (np.ndarray) – The velocity autocorrelation function. This is a 2D NumPy array representing the correlation over time.

  • time (np.ndarray) – The time array corresponding to the VACF. This should be a 1D NumPy array of the same length as vacf.

  • timestep (float) – The time step between consecutive points in the time array.

Returns:

  • integral (np.ndarray) – The running integral of the VACF, divided by 3, representing a quantity proportional to the mean squared displacement over time.

  • func (callable) – The exponential function used for fitting. It has the form \(a + b e^{-c x}\).

  • opt (np.ndarray) – The optimal parameters (\(a\), \(b\), \(c\)) found by the curve fitting.

  • Rsqrd (float) – The coefficient of determination (\(R^2\)) for the exponential fit, indicating how well the model fits the data.

Notes

The diffusion coefficient (\(D\)) is given by the long-time limit of the integral, which corresponds to the parameter \(a\) in the fitted function. Specifically, \(D = a\). The VACF is divided by 3, assuming a 3-dimensional system, to relate it to the mean squared displacement in one dimension.