VACF: Theory and Validity¶
Background for VACF and Spectra: what sets the wavenumber axis and its resolution, how apodization changes band shapes, what the normalization discards, and how the two estimators differ in the tail.
Validity¶
The frequency axis, and what sets it¶
Write \(\Delta t\) for time_step in ps,
\(W\) for window in frames and \(F\) for ftsize.
The grid spacing is set by ftsize. PQAnalysis mirrors the padded
correlation into an even extension and labels the transform with
The axis starts at \(\Delta\tilde\nu\) (there is no \(\tilde\nu = 0\)
point) and runs up to \(F\,\Delta\tilde\nu\). For time_step = 0.001
and ftsize = 5000 that is a 3.34 cm⁻¹ grid reaching 16682 cm⁻¹.
The upper limit is the Nyquist wavenumber, set by the sampling interval alone:
which is 16678 cm⁻¹ for frames written every 1 fs and 1668 cm⁻¹ for frames written every 10 fs. What matters is the interval between written frames, not the integration time step of the underlying dynamics: an X-H stretch near 3000 cm⁻¹ needs a trajectory stride below about 5.5 fs to be represented at all, and motion above the limit is aliased back into the spectrum rather than discarded. The legacy axis is one grid point longer than Nyquist, \(F/(F-1)\cdot\tilde\nu_{\max}\), and the last spectrum point duplicates its predecessor; both are deliberate legacy conventions.
Note
The axis carries a deliberate legacy calibration: the spacing is computed with a period of \(2(F-1)\) points although the underlying even extension has \(2F-1\) points. Reported wavenumbers are therefore high by a factor \((2F-1)/(2F-2)\). The absolute offset grows linearly with wavenumber and reaches about half a grid spacing at the top of the axis, so it is at most a sub-bin effect, but it is a systematic one.
Resolution comes from the correlation length, not from padding¶
The correlation extends to a maximum lag time \(T = W\,\Delta t\). Truncating a correlation at \(T\) convolves the spectrum with a kernel of width
and features narrower than that are not resolved no matter how fine the grid
is. Two lines 20 cm⁻¹ apart merge into a single peak for \(T = 0.5\) ps
(\(\delta\tilde\nu = 67\) cm⁻¹) and separate cleanly for \(T = 2.5\) ps
(\(\delta\tilde\nu = 13\) cm⁻¹). Increasing ftsize beyond the number of
correlation points only interpolates the same information onto a denser grid;
resolving finer structure requires a longer window, which in turn requires
a longer trajectory.
Warning
ftsize also truncates. The correlation is zero-padded or cut to
ftsize points, so with ftsize smaller than window + 1 everything
beyond the first ftsize lags is silently discarded before the transform.
The default ftsize is 2000 while the default window is 1000, so the
defaults are safe, but raising window without raising ftsize throws
the extra correlation away. Keep ftsize at least window + 1.
Apodization: leakage against band shape¶
A correlation that has not decayed to zero at the end of the window produces sinc-like ringing and leakage in the transform. An apodization window suppresses the discontinuity, at the cost of widening every band and lowering every peak [Harris1978]. Always report which window and which parameters were used; band widths from differently apodized spectra are not comparable.
The exponential window multiplies the correlation by
\(\exp(-a\,t)\) with \(a\) = window_param in ps⁻¹. It adds
\(a\) to the decay rate of the correlation, which broadens a Lorentzian
band by roughly
On a test band with a 2 ps dephasing time, the unapodized width is about
9 cm⁻¹, and window_param values of 2, 4 and 8 ps⁻¹ widen it to about 26, 47
and 89 cm⁻¹ while the peak height drops to 28 %, 16 % and 8 % of its
unapodized value. Apodization strong enough to tame leakage is also strong
enough to dominate the linewidth, so a width read off an apodized spectrum is a
property of the window, not of the dynamics.
Warning
hann and blackman do nothing at their default settings. Both are
built from window_start (default 0.0 ps) and window_stop
(default 1000.0 ps), and over a correlation of a few ps the resulting
factors deviate from unity by less than \(3\times10^{-5}\), so the
spectrum is indistinguishable from window_function = none. To use them,
set window_stop to the correlation length window * time_step. With
window_stop = 2.5 on a 2.5 ps correlation, hann widened the test
band from 9.2 to 15.0 cm⁻¹ and blackman to 17.5 cm⁻¹, with peak heights
falling by 38 % and 45 %. The exponential window is unaffected by this,
since it decays from window_start onwards.
The legacy hann and blackman formulas are additionally non-standard:
the Hann window is mirrored and the Blackman denominators use the stop
index rather than the window width. Both quirks are reproduced deliberately;
see PQAnalysis.analysis.vacf.spectrum.apodization_window().
The normalization discards absolute magnitude¶
The correlation is normalized to \(C(0)=1\) before anything else happens:
the direct estimator divides each time origin by its own aggregate squared
velocity norm, and the fft estimator divides by its lag-zero value. The
cosine transform is linear, so the spectrum inherits that normalization and its
amplitudes are in arbitrary units. Relative peak areas within one spectrum are
meaningful, since apodization broadens bands but conserves their area; peak
heights are only comparable between spectra that used the same apodization.
Zero-padding does not rescale amplitudes; it samples the same envelope more
finely, so a coarse grid can under-read the apex of a narrow band by a few per
cent. Absolute intensities are not available at all, for the velocity spectrum
as much as for the charge-flux spectrum.
The sums over atoms are also unweighted; no masses enter anywhere in the VACF
code. Each atom therefore contributes in proportion to its own mean squared
velocity, which at equipartition is proportional to \(1/m\), so light atoms
dominate a mixed-element selection. This is not the mass-weighted vibrational
density of states of the textbook definition. For element-resolved band
assignments, run separate analyses with a per-element target_selection.
Classical dynamics limits what a peak position means¶
The spectrum describes the classical nuclear motion actually present in the trajectory. PQAnalysis applies no quantum correction of any kind: no zero-point energy, no harmonic quantum correction factor to the intensities, no frequency scaling factor. Band positions therefore carry the classical-nuclei error and the full anharmonic and thermal shift of the underlying dynamics at the simulated temperature, and are not directly comparable to harmonic normal-mode wavenumbers (see Vibrational Analysis for the harmonic route). Comparisons with experiment must state the temperature, the level of theory and the fact that the peak positions are classical.
Estimator choice changes the statistics of the tail¶
method = direct(default) spawns an origin everygapframes while a full window still fits. Every origin covers the full window, so every lag from 0 towindowis averaged over the same number of origins.method = fftuses every frame as an origin and ignoresgap, but divides lag \(\tau\) by its own origin count \(N-\tau\). Its origin count does fall off with lag, so the far tail of the correlation, which determines the low-wavenumber structure of the spectrum, is progressively noisier. It also holds all velocities in memory.
Warning
A velocity trajectory of exactly window frames with gap = 1 takes
the legacy single-origin branch: one origin is spawned, the correlation is
an unaveraged single-origin estimate, and the final lag bin stays exactly
zero. Use a longer trajectory or a smaller window.
References¶
[Rahman1964] introduced the velocity autocorrelation function as a molecular-dynamics observable and describes its negative-lobe behavior in a liquid.
[Green1954] and [Kubo1957] establish the link between equilibrium time correlation functions and transport coefficients. The time integral of the unnormalized velocity autocorrelation function is the Green-Kubo expression for the self-diffusion coefficient.
[Wiener1930] and [Khintchine1934] prove that the autocorrelation function and the power spectrum of a stationary process are a Fourier pair, which is what the
fftestimator exploits.[Dickey1969] interprets the transformed velocity autocorrelation function as a vibrational density of states.
[Thomas2013] covers vibrational density of states, charge-flux and dipole routes to infrared spectra, and the effect of correlation depth and windowing on band shapes.
[Harris1978] tabulates the apodization windows and their trade-off between sidelobe suppression and main-lobe broadening.
[Allen2017] gives the time-origin averaging and sampling requirements for correlation functions.