Time Correlation Functions

Many physical observables depend on the (Fourier transform of the) Time Correlation Functions (TCF) of some microscopic quantities. Here are some examples:

Infrared, Raman and Sum Frequency Generation spectra

frequency dependent dielectric constant and dielectric susceptibility

Vibrational Density Of States

The way to compute these quantities depends on their definition of course, but a general way to evaluate the TCF from the time series of microscopic quantities wil lbe provided in this following.

Cross Correlation

The Time Correlation Functions (TCF) between two observables \(A,B\) is defined through the cross correlation (similar, but not the same of the convolutuon) of their time series:

\[\begin{split}(A \star B)(\tau) \triangleq & \int_{-\infty}^{\infty} \overline{A(t)} B(t+\tau) \, dt \\ \triangleq & \int_{-\infty}^{\infty} \overline{A(t-\tau)} B(t)\,dt\end{split}\]

A naive implementation of the previous formula would have a computation cost scaling as \(\mathcal{O}\left(N\right)\) where \(N\) is the size of the arrays (assumed to be of the same length for simplicity).

However, we can exploit an analogous of the convolution theorem which allows to express the TCF as:

\[A \star B = \mathcal{F}^{-1} \left[ \overline{\mathcal{F}}\left[ A \right] \cdot \mathcal{F}\left[ B \right] \right]\]

where \(\mathcal{F}\left[ \cdot \right]\) is the Fourier transform of a function defined as:

\[\mathcal{F}[ f ](\omega) \triangleq \int_{-\infty}^{+\infty} e^{-i \omega t} f(t) \, dt\]

The Fourier transfrom has a computational cost of \(\mathcal{O}\left(N\log N\right)\) thanks to the Fast Foruier Transform (FFT) implementations.
For this reason, we will provide a way to evaluate the TCF using this “Fourier transform trick”.
Here is a simple function that returns the TCF of an array with itself:

Function Documentation

tcf.tcf.autocorrelate(A: ndarray, axis: int | None = 0) → ndarray[source]

Compute the autocorrelation function of an array along the specified axis.

Parameters:

A (np.ndarray) – array containing the time series
axis (Optional[int], optional) – axis corresponding to time, defaults to 0

Returns:

autocorrelation function of the array input

Return type:

np.ndarray

import numpy as np
from typing import Optional

def autocorrelate(A: np.ndarray, axis: Optional[int] = 0) -> np.ndarray:
    """
    Compute the autocorrelation function of an array along the specified axis.
    
    :param A: array containing the time series
    :type A: np.ndarray
    :param axis: axis corresponding to time, defaults to 0
    :type axis: Optional[int], optional
    :return: autocorrelation function of the array input
    :rtype: np.ndarray
    """
    len_a = A.shape[axis]
    len_fft = 2 * len_a
    
    # Normalization factor (decreasing window size)
    norm_tcf = np.arange(len_a, 0, -1, dtype=int)
    
    # Expand normalization dimensions    
    norm_tcf = np.expand_dims(norm_tcf, tuple(i for i in range(A.ndim) if i != axis))
    
    # Compute FFT of A
    ftA = np.fft.rfft(A, axis=axis, n=len_fft)
    
    # Power spectrum (multiply by complex conjugate)
    ftA *= np.conj(ftA)
    
    # Inverse FFT to get the auto-correlation
    autocorr = np.fft.irfft(ftA, axis=axis, n=len_fft)
    
    # Slice the auto-correlation to original signal length
    autocorr = np.take(autocorr, np.arange(len_a), axis=axis)
    
    return autocorr / norm_tcf

Usage Examples

import numpy as np
from tcf import autocorrelate

# Create an array with 1000 random numbers
arr = np.random.rand(1000)

# Compute the autocorrelation function
autocorr = autocorrelate(arr)

# Plot the autocorrelation
import matplotlib.pyplot as plt
plt.plot(autocorr)
plt.xlabel('Lag')
plt.ylabel('Autocorrelation')
plt.show()

Attention

In practical calculations, we do not have infinite time series with the time \(t\) spanning in \((-\infty,+\infty)\), but the time ranges from \(t \in [0,T]\) only.

Moreover, the numerical routines actually implement a monolateral Fourier transform \(\tilde{\mathcal{F}}[ \cdot ]\) of the following form:

\[\tilde{\mathcal{F}}[ f ](\omega) \triangleq \int_{0}^{T} e^{-i \omega t} f(t) \, dt\]

Infrared Spectrum

The Vibrational Infrared Absorption spectrum \(S\left(\omega\right)\), at thermal equilibrium, can be evaluated using the time series of the dipole using the following expression:

\[\begin{split}S\left(\omega\right) & = \alpha\left(\omega\right) n\left(\omega\right) \\ & = \frac{\pi \omega}{3 \hbar c \Omega c \varepsilon_0} \left( 1 - e^{-\beta \hbar \omega} \right) \, \mathcal{F}\left[\boldsymbol{\mu}\star\boldsymbol{\mu}\right]\left(\omega\right) \\ & \stackrel{\beta \ll 1}{ \approx} \, \frac{\beta\omega^2}{6c\Omega\varepsilon_0} \, \mathcal{F}\left[\boldsymbol{\mu}\star\boldsymbol{\mu}\right]\left(\omega\right)\end{split}\]

where

\(\alpha\left(\omega\right)\) is the Beer-Lambert absorption coefficient,

\(n\left(\omega\right)\) is the refractive index of the material,

\(\beta = \frac{1}{k_B T}\) is the thermodynamic beta,

\(c\) is the speed of light,

\(\Omega\) is the system volume

\(\varepsilon_0\) is the vacuum permittivity,

and clearly

\[\mathcal{F}\left[\boldsymbol{\mu}\star\boldsymbol{\mu}\right]\left(\omega\right) \triangleq \int_{-\infty}^{+\infty} e^{-i \omega t} \boldsymbol{\mu}\star\boldsymbol{\mu}(t) \, dt\]

It is important to stress a couple of things:

the TCF should actually be calculated for the fluctuation of the dipole w.r.t. to its average value:

\[\boldsymbol{\mu}\star\boldsymbol{\mu}(t) \longrightarrow C_{\delta\boldsymbol{\mu},\delta\boldsymbol{\mu}}(t) \quad \text{with} \quad \delta\boldsymbol{\mu} = \boldsymbol{\mu} - \braket{\boldsymbol{\mu}}\]

and this means that, before computing any Fourier transform, one should remove the mean value of the time series. However, we will keep using the same notation as before for simplicity.

the definition of the infrared spectrum requires an integration over time which actually range in \((-\infty,+\infty)\).

The TCF of the dipole with itself is also time-reversible, i.e.

\[\left(\boldsymbol{\mu}\star\boldsymbol{\mu}\right)(t) = \left(\boldsymbol{\mu}\star\boldsymbol{\mu}\right)(-t)\]

and this implies that its Fourier transform in purely real.

However, in numerical implementation we can only evaluate the monolateral Fourier transform \(\tilde{\mathcal{F}}[ \cdot ]\), whose output is complex since no integration over time for \(t\lt 0\) occurs.

For this reason, to guarantee the correct calculation of the infrared spectrum, what is actually implemented numerically is

\[\begin{split}S\left(\omega\right) = & \frac{\beta\omega^2}{6c\Omega\varepsilon_0} \, 2 \, {\rm Re} \tilde{\mathcal{F}}\left[\boldsymbol{\mu}\star\boldsymbol{\mu}\right]\left(\omega\right) \\ = & \frac{\beta\omega^2}{3c\Omega\varepsilon_0} {\rm Re} \int_{0}^{+\infty} e^{-i \omega t} \left(\boldsymbol{\mu}\star\boldsymbol{\mu}\right)\left(t\right) \, dt\end{split}\]

A common trick

Numerical simulations providing the dynamics over time of the dipole (as well as any other quantity) can be used to evaluate the previous formulas only if the simulat time is long enough and the dynamics is properly captured, i.e. if the integration time step is small enough.
In certain systems, like liquid water, properly sampling the long time behavior of \(\left(\boldsymbol{\mu}\star\boldsymbol{\mu}\right)\left(t\right)\), i.e. the low frequency behavior of its Fourier transform, can be challenging due to its intrisically “low dynamics”.
This means that to properly sample the low frequrncy region of \(S\left(\omega\right)\), extremely long simulations would be needed.
There is however a possible, or partial, solution that can be applied to the infrared spectrum, which consists in relating \(\left(\dot{\boldsymbol{\mu}}\star\dot{\boldsymbol{\mu}}\right)\left(t\right)\) with \(\boldsymbol{\mu}\star\boldsymbol{\mu}(t)\), i.e. the TCF of the dipole time derivative.
It is not hard to show that:

\[\omega^2 \mathcal{F}\left[\boldsymbol{\mu}\star\boldsymbol{\mu}\right]\left(\omega\right) = \mathcal{F}\left[\dot{\boldsymbol{\mu}}\star\dot{\boldsymbol{\mu}}\right]\left(\omega\right)\]

This relation shows that the TCF of the dipole time derivative decays faster, and we can also use this TCF to evaluate the infrared spectrum:

\[S\left(\omega\right) = \frac{\beta}{3c\Omega\varepsilon_0} {\rm Re} \int_{0}^{+\infty} e^{-i \omega t} \left(\dot{\boldsymbol{\mu}}\star\dot{\boldsymbol{\mu}}\right)\left(t\right) \, dt\]