Searching for a shell

From velocity autocorrelation function to phonon DOS

MD  Phonon 

The autocorrelation function is a cross-correlation function that shows the correlation of one signal (function) with a delayed version of it self. Essentially, It’s the phonoic analogue of the diagonal part of the Green’s function (for more details of quantum Green’s function, see :link: this post and Eq. 11 of :link: this post).

To derive the relation between velocity autocorrelation function and the phonon density of states, let us first start with the Fourier transforming the velocity of an atom $i$, $v_i(t)$.

F[vi(t)](ω)=∫−∞∞vi(t)eiωtdt,(1)\mathcal{F}[v_i(t)](\omega) = \int_{-\infty}^{\infty}v_i(t)e^{i\omega t} dt, \tag{1}

The power spectrum (which tells us how the energy is distributed over the frequency) of this velocity function is:

∣vi(ω)∣2=vi∗(ω)vi(ω)=∫−∞∞∫−∞∞vi∗(t)vi(t′)eiω(t−t′)dtdt′.(2)\begin{aligned} |v_i(\omega)|^2 &= v^*_i(\omega) v_i(\omega) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} v^*_i(t) v_i(t') e^{i\omega(t-t')} dt dt'. \end{aligned} \tag{2}

Substitute $t-t’$ with $t’’$, the power spectrum can be re-written as:

∣vi(ω)∣2=∫−∞∞∫−∞∞vi∗(t′′+t′)vi(t′)eiω(t′′)dt′′dt′.(3)|v_i(\omega)|^2 = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} v^*_i(t''+t') v_i(t') e^{i\omega(t'')} dt'' dt'. \tag{3}

Note that the exponential $e^{i\omega t}$ in Eq. 1 doesn’t contain a minus sign. This is just a difference in convention, if we put in the minus sign, then vi∗(t′′+t′)v^*_i(t''+t') in Eq. 3 would become vi∗(−t′′+t′)v^*_i(-t''+t'). However, this does not affect any of the following analysis.

Before moving on, let’s assume that harmonic approximation is valid (which, in theory, we should have zero-width phonon branches), the atomic displacement with respect to time can be written as:

rj(t)=∑sQs,je−iωst,r_j(t) = \sum_s Q_{s,j} e^{-i\omega_s t},

where $s$ labels the phonon branch and $j$ is a composite index for the vibrations of atom $i$ on three spacial directions. Q is the amplitude of the vibration (with $\sqrt{\frac{1}{2} m_j}$ absorbed in) and $\omega$ is the frequency of the phonon mode.

With the time-dependent displacement derived, the time-dependent velocity can be expressed as:

vj(t)=ddtrj(t)=∑sQs,j(−iωs)e−iωst.(4)v_j(t) = \frac{d}{dt} r_j(t) = \sum_s Q_{s,j} (-i\omega_s) e^{-i\omega_s t}. \tag{4}

Inserting Eq. 4 to Eq. 3, and notice we are now using $j$ instead of $i$ to incorporate three spacial dimensions, we get:

∣vj(ω)∣2=∫−∞∞∫−∞∞∑s,s′Qs,j∗(iωs)eiωs(t′′+t′)Qs′,j(−iωs′)e−iωs′t′eiωt′′dt′′dt′=∫−∞∞∫−∞∞∑s,s′Qs,j∗Qs′,j(ωsωs′)ei(ω+ωs)t′′ei(ωs−ωs′)t′dt′′dt′=∫−∞∞∑s,s′Qs,j∗Qs′,j(ωsωs′)ei(ω+ωs)t′′δ(ωs−ωs′)dt′′=∑s∫−∞∞∣Qs,j∣2ωs2ei(ω+ωs)t′′dt′′(5)\begin{aligned} |v_j(\omega)|^2 &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \sum_{s,s'} Q^*_{s,j}(i\omega_s)e^{i\omega_s(t''+t')} Q_{s',j} (-i\omega_s') e^{-i\omega_{s'} t'} e^{i\omega t''} dt'' dt'\\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \sum_{s,s'} Q^*_{s,j}Q_{s',j}(\omega_s \omega_{s'})e^{i(\omega+\omega_s)t''} e^{i(\omega_s - \omega_{s'})t'} dt'' dt'\\ &= \int_{-\infty}^{\infty} \sum_{s,s'} Q^*_{s,j}Q_{s',j}(\omega_s \omega_{s'})e^{i(\omega+\omega_s)t''} \delta(\omega_s - \omega_{s'}) dt''\\ &= \sum_s \int_{-\infty}^{\infty} |Q_{s,j}|^2 \omega_s^2 e^{i(\omega+\omega_s)t''} dt''\\ \end{aligned} \tag{5}

noting that, due to :link: equipartition theorem, the total energy of an oscillator is $k_BT$ where $k_B$ is the Boltzmann constant and $T$ is temperature (there might be a factor of 2 missing here, but again, this won’t actually affect the final result as it merely acts as a scaling factor). Since we can express the total energy as pure kinetic energy at the equilibrium position, we get:

Etot=Ekin=12mv2=Q2ω2=kBT.E_{\text{tot}} = E_{\text{kin}} = \frac{1}{2} m v^2 = Q^2 \omega^2 = k_BT.

Substituting ∣Qs,j∣2ωs2|Q_{s,j}|^2 \omega_s^2 with $k_BT$, Eq. 5 becomes:

∣vj(ω)∣2=kBT∑s∫−∞∞ei(ω+ωs)t′′dt′′=kBT∑sδ(ω+ωs).\begin{aligned} |v_j(\omega)|^2 &= k_B T \sum_s \int_{-\infty}^{\infty} e^{i(\omega+\omega_s)t''} dt''\\ &= k_B T \sum_s \delta(\omega+\omega_s). \end{aligned}

If we add contributions from all atoms (assuming there are $N$ atoms) in all three directions:

∑j∣vj(ω)∣2=3NkBT∑sδ(ω+ωs)\sum_j |v_j(\omega)|^2 = 3N k_B T \sum_s \delta(\omega+\omega_s)

remember that the phonon DOS is expressed as ∑sδ(ωs−ω)\sum_s \delta(\omega_s-\omega), we can write phonon DOS as:

ρ(−ω)=∑sδ(ωs−ω)=13NkBT∑j∣vj(−ω)∣2\rho(-\omega) = \sum_s \delta(\omega_s-\omega) = \frac{1}{3Nk_BT} \sum_j |v_j(-\omega)|^2

by the knowledge of Fourier transform of the time-dependent velocity function in Eq. 3, we get:

ρ(−ω)=13NkBT ∑j∫−∞∞∫−∞∞vj∗(t′′+t′)vj(t′)e−iω(t′′)dt′′dt′.(6)\rho(-\omega) = \frac{1}{3Nk_BT}\ \sum_j \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} v^*_j(t''+t') v_j(t') e^{-i\omega(t'')} dt'' dt'. \tag{6}

Remembering that the system is at its equilibrium state and the velocity signal is more or less periodic. For this reason, the following integral diverges:

∫−∞∞v⃗(t+t′)⋅v⃗(t′)dt′\int_{-\infty}^{\infty} \vec v(t+t')\cdot \vec v(t')dt'

Noting that the DOS function is a Dirac delta function, at $\omega = \omega_s$, the the DOS should indeed diverge. However, since we only care about the difference of the $\rho$ with different $\omega$s in our final plots, instead of using the integral above, we can re-write the integral as:

∑t′=0Nt′v⃗(t+t′)⋅v⃗(t′)dt′,\sum_{t'=0}^{N_{t'}} \vec v(t+t')\cdot \vec v(t') dt',

This expression can also be seen as a scaled velocity autocorrelation function:

Cv(t)dt′=dt′∑t′=0Nv⃗(t+t′)⋅v⃗(t′)=⟨v⃗(t)⋅v⃗(0)⟩dt′(7)\begin{aligned} C_{v}(t) dt' &= dt' \sum_{t'=0}^N \vec v(t+t')\cdot \vec v(t')\\ &= \braket{\vec v(t)\cdot \vec v(0)} dt' \end{aligned} \tag{7}

Alternatively, we can re-normalize this autocorrelation function by doing something like:

Cv(t)=⟨v⃗(t)⋅v⃗(0)⟩⟨v⃗(0)⋅v⃗(0)⟩,C_{v}(t) = \frac{\braket{\vec v(t)\cdot \vec v(0)}}{\braket{\vec v(0)\cdot \vec v(0)}},

Which, normalizes the maximum correlation to $1$ and we can get rid of the prefactor $dt’$. This is what you see in, e.g. :link: 10.1016/j.cpc.2011.04.019 or, :link: 10.1039/c8nr07373b.

Inserting Eq. 7 into Eq. 6, we get:

ρ(ω)=13NkBT ∑j∫−∞∞⟨vj(t′′)⋅vj(0)⟩eiω(t′′)dt′′.(8)\rho(\omega) = \frac{1}{3Nk_BT}\ \sum_j \int_{-\infty}^{\infty} \braket{v_j(t'')\cdot v_j(0)} e^{i\omega(t'')} dt''. \tag{8}

With Eq. 8, noting the decreasing nature of the velocity autocorrelation function so that the integral always converges, we see that if we have the knowledge of the velocity of each atom across a period of time (since we cannot do infinity, a we need to have a cut-off time for the integral), we can calculate the phonon density of states.

In practice, Eq. 8 should be rewrite in a discrete form as:

ρ(ω)=13NkBT ∑j1Nt′′∑t′′=0Nt′′⟨vj(t′′)⋅vj(0)⟩eiω(t′′).(9)\rho(\omega) = \frac{1}{3Nk_BT}\ \sum_j \frac{1}{N_{t''}} \sum_{t''=0}^{N_{t''}} \braket{v_j(t'')\cdot v_j(0)} e^{i\omega(t'')}. \tag{9}

In summary, the Phonon density of states can be expressed as the Fourier transformed velocity autocorrelation function.

NOTE: I’ve ignored the $2\pi$ prefactor in all Fourier transformations.

Extension to periodic system

For periodic systems, the phonon bandstructure can also be interpreted as phonon DOS at each $k$ point. To obtain this “spectrum” function. we need to Fourier transform the velocity function with respect to the unit cell vectors. For example, lets say we have calculated a $\sqrt[3]{N}\times\sqrt[3]{N}\times\sqrt[3]{N}$ supercell (so that we have a total of $N$ unit cells):

Here, we label the atoms in each cell with $j$ (again a composite label with three direction added) and the unit cell position is expressed as $\vec R$ so that the velocity of atom $j$ at cell $\vec R$ (in three direction) is:

vj(R⃗,t).(10)v_{j}(\vec R, t). \tag{10}

Fourier (discrete) transform Eq. 9 with respect to $\vec R$, we get:

F[vj(R⃗,t)](k⃗,t)=1N∑R⃗vj(R⃗,t)eik⃗R⃗=vj(k⃗,t),(11)\begin{aligned} \mathcal{F}[v_{j}(\vec R, t)](\vec k,t) = \frac{1}{\sqrt{N}} \sum_{\vec R} v_j(\vec R,t) e^{i\vec k \vec R} = v_{j}(\vec k, t), \end{aligned} \tag{11}

where we have $N$ unit cell in the supercell.

Using Eq. 10, we can express the phonon DOS at each $\vec k$ point as:

ρ(k⃗,−ω)=13NkBT ∑j∫−∞∞⟨vj(k⃗,t′′)⋅vj(k⃗,0)⟩e−iω(t′′)dt′′.(12)\rho(\vec k, -\omega) = \frac{1}{3Nk_BT}\ \sum_j \int_{-\infty}^{\infty} \braket{v_{j}(\vec k, t'')\cdot v_{j}(\vec k, 0)} e^{-i\omega(t'')} dt''. \tag{12}

Here, the commensurate k-point can depends on supercell size. In other words, the larger your supercell is, the finner the k-points you can sample, and a better graph you’ll be able to get.

Physical explanation

Okay, the derivation is nice and easy, but what’s the physics behind?

Let’s consider a single atom that’s vibrating at its equilibrium position. Its velocity vs time can be plotted as:

The autocorrelation function $\int_{-\infty}^{\infty} v(t+t^{\prime})v(t^{\prime}) dt’$ of this velocity signal at different $t$ is:

We see that the autocorrelation function retains the periodicity of the orignal function. Such behavior will be reflected as a sharp peak in the Fourier spectrum.

For a more comprehensive understanding of the autocorrelation function, see :link: this and :link: this post.

See it in action

Download: :file_folder: code

Now it’s time to try it out! Using water molecule (so that I don’t have to deal with the k-depenedence of periodic systems, which requries a huge super cell to achieve sufficient sampling) as an example, we are going to compare the phonon frequencies calculated with DFT (DFPT) to the ones that we have obtained using our MD mtehod.

The three Raman active mode from our DFT calculation have frequencies: 47 THz($a_1$), 111 THz($b_2$), 114 THz($a_1$), comparable to :link: the experiemental values.

With a molecular dynamic calculation of a single water molecule thermalized at 300K, we get the following plot:

Okay, it works!!

Author | Chengcheng Xiao

Currently a postdoctoral research associate at Imperial College London. Predicting electron behaviour since 2016.