Theoretical background

Below is a brief overview of the theory underpinning the correlation functions supported by dynasor as well as their analysis. For a more detailed description and derivations of the correlation functions see, for example, [BooYip80].

Dynamic structure factor

The particle density is denoted \(n(\boldsymbol{r},t)\) and defined as

\[n(\boldsymbol{r},t) = \sum _i ^N \delta (\boldsymbol{r} - \boldsymbol{r}_i(t)),\]

where \(N\) is the number of particles and \(\boldsymbol{r}_i(t)\) is the position of particle \(i\) at time \(t\). The intermediate scattering function \(F(\boldsymbol{q},t)\) is defined as the Fourier transform of the auto-correlation of the particle density

\[F(\boldsymbol{q},t)=\frac{1}{N}\left<n(\boldsymbol{q},t)n(-\boldsymbol{q},0)\right> \,= \, \frac{1}{N}\sum _i ^N \sum _j ^N \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>,\]

where \(\boldsymbol{q}\) is a wave vector. The brackets denote an ensemble average, which in the case of dynasor is replaced by a time average. (The system under study should be ergodic in order for the time average to be a suitable substitute for the ensemble average.) Here, \(F(\boldsymbol{q},t)\) is commonly referred to as the coherent part.

The incoherent intermediate scattering function, i.e., \(i=j\), is defined as

\[F_\mathrm{incoh}(\boldsymbol{q},t)=\frac{1}{N} \sum _i ^N \left < \mathrm{exp}\left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_i(0))\right] \right >\]

and describes single-particle motion, i.e., diffusive motion. \(F_\mathrm{incoh}(\boldsymbol{q},t)\) is commonly referred to as the incoherent scattering function, but sometimes also referred to as the self-part of the intermediate scattering function.

The static structure factor is given by

\[S(\boldsymbol{q}) = F(\boldsymbol{q},0),\]

whereas the dynamic structure factor \(S(\boldsymbol{q},\omega)\) is the Fourier transform of \(F(\boldsymbol{q},t)\)

\[S(\boldsymbol{q},\omega) = \int _{-\infty} ^\infty F(\boldsymbol{q},t) \, e^{-i\omega t} \mathrm{d}t,\]

where \(\omega\) is the angular frequency. \(S(\boldsymbol{q},\omega)\) exhibits peaks in the \((\boldsymbol{q},\omega)\) plane corresponding to the dispersion of lattice vibrations (phonons). The widths of these peaks are related to the phonon lifetimes. Computing \(S(\boldsymbol{q},\omega)\) via molecular dynamics simulation has the advantage of fully including anharmonic effects and allowing one to study the temperature dependence of phonon dispersions.

Note

The static and dynamic structure factors can be calculated using the compute_static_structure_factors and compute_dynamic_structure_factors functions, respectively. The procedure is illustrated for a liquid in this tutorial and for a crystalline material in this tutorial.

Velocity correlation functions

It is often convenient to also consider current correlations based on the particle velocities. The current density \(\boldsymbol{j}(\boldsymbol{r},t)\) is given by

\[\begin{split}\boldsymbol{j}(\boldsymbol{r},t) &= \sum_i^N \boldsymbol{v}_i(t) \, \delta (\boldsymbol{r} - \boldsymbol{r}_i(t)) \\ \boldsymbol{j}(\boldsymbol{q},t) &= \sum_i^N \boldsymbol{v}_i(t) \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)},\end{split}\]

where \(\boldsymbol{v}_i(t)\) is the velocity of particle \(i\) at time \(t\). This can be split into longitudinal and transverse parts as

\[\begin{split}\boldsymbol{j}_L(\boldsymbol{q},t) &= \sum_i ^N(\boldsymbol{v_i}(t) \cdot\hat{\boldsymbol{q}}) \, \hat{\boldsymbol{q}} \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)} \\ \boldsymbol{j}_T(\boldsymbol{q},t) &= \sum_i ^N \left[\boldsymbol{v_i}(t) - (\boldsymbol{v_i}(t) \cdot \hat{\boldsymbol{q}}) \, \hat{\boldsymbol{q}}\right] \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)}.\end{split}\]

Now the correlation functions can be computed as

\[\begin{split}C_L(\boldsymbol{q},t) &= \frac{1}{N}\left<\boldsymbol{j}_L(\boldsymbol{q},t)\cdot\boldsymbol{j}_L(-\boldsymbol{q},0)\right> \\ C_T(\boldsymbol{q},t) &= \frac{1}{N}\left<\boldsymbol{j}_T(\boldsymbol{q},t)\cdot\boldsymbol{j}_T(-\boldsymbol{q},0)\right>.\end{split}\]

The longitudinal current can be related to the particle density as

\[\begin{split}\frac{\partial }{\partial t}n(\boldsymbol{q}, t) &= \mathrm{i} \boldsymbol{q}\cdot \boldsymbol{j}(\boldsymbol{q}, t) \\ \omega^2 S(\boldsymbol{q},\omega) &= q^2 C_L(\boldsymbol{q},\omega),\end{split}\]

which means some features can be easier to resolve in one function than in the other. The current correlation functions can be thought of as a spatially dependent generalization of the velocity correlation function and are closely related to the phonon spectral energy density.

Note

The velocity (or current) correlation functions can be calculated using the compute_dynamic_structure_factors function by setting calculate_currents to True. An example for the application of this function can be found in the tutorial on the dynamic properties of an NiAl alloy.

Multi-component systems

In multi-component systems, one can introduce partial correlation functions, which enable separation of the contributions by groups of particles, e.g., by type or site symmetry. For example, in the case of a binary system, one can define \(F_\mathrm{AA}(\boldsymbol{q},t)\), \(F_\mathrm{BB}(\boldsymbol{q},t)\), \(F_\mathrm{AB}(\boldsymbol{q},t)\) and so on for all correlation functions introduced above. In dynasor we define the partial correlation functions as

\[F_\mathrm{AA}(\boldsymbol{q},t) = \frac{1}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in A}^{N_\mathrm{A}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>\]

where the sums run over all A atoms. Here, \(N_\mathrm{A}\) refers to the number of atoms of type A, \(N_\mathrm{B}\) refers to the number of atoms of type B, and \(N = N_\mathrm{A} + N_\mathrm{B}\) refers to the total number of atoms in the system. For the cross-term, \(F_\mathrm{AB}(\boldsymbol{q},t)\), we get

\[\begin{split}F_\mathrm{AB}(\boldsymbol{q},t) = & \frac{1}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in B}^{N_\mathrm{B}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right> + \frac{1}{N}\sum_{i\in B}^{N_\mathrm{B}} \sum_{j\in A}^{N_\mathrm{A}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right> \\=& \frac{2}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in B}^{N_\mathrm{B}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>\end{split}\]

Here, the factor of two in the final expression is due to considering both the A-B and B-A terms. These normalization choices mean the total intermediate scattering function \(F(\boldsymbol{q},t)\) (as defined above) is given by

\[F(\boldsymbol{q},t) = F_\mathrm{AA}(\boldsymbol{q},t) + F_\mathrm{AB}(\boldsymbol{q},t) + F_\mathrm{BB}(\boldsymbol{q},t)\]

In some cases, instead of analyzing the partial functions directly, it is interesting to consider linear combinations of them, for example by weighting them with relevant atomic form factors, masses, or charges.

Relation to scattering experiments

The correlation functions can be convolved with atomic weights, \(w_i\), i.e.,

\[F_\mathrm{coh}(\boldsymbol{q},t)=\frac{1}{N}\sum_i^N \sum_j^N w_i w_j \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>,\]

where the weights could, for example, be masses, charges, or even atomic form factors (scattering lengths) to allow for direct comparison to structure factors measured in various neutron and X-ray scattering experiments. For some probes, e.g., neutrons, the scattering lengths are in general constant complex numbers but differ per isotope and for the coherent and incoherent correlation functions,

\[F_\mathrm{coh}(\boldsymbol{q},t)=\frac{1}{N}\sum_i^N \sum_j^N \operatorname{Re} \left(b_i^* b_j \right) \left< \mathrm{exp}\left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>.\]

X-ray form factors and electron scattering factors, on the other hand, are \(q\)-dependent and thus take the form

\[F_\mathrm{coh}(\boldsymbol{q},t)=\frac{1}{N}\sum_i^N \sum_j^N f_i(q) f_j(q) \left< \mathrm{exp}\left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>.\]

Conveniently, since the weights only depend on the atom type, the weighting can be done as a post-processing step using partial correlation functions. For a binary (A-B) system, the total weighted correlation function is

\[F_\mathrm{coh}(\boldsymbol{q},t)=f^{AA}(q) F^{AA}_\mathrm{coh}(\boldsymbol{q},t) + f^{AB}(q) F^{AB}_\mathrm{coh} (\boldsymbol{q},t) + f^{BB}(q) F^{BB}_\mathrm{coh}(\boldsymbol{q},t).\]

Furthermore, it is frequently desirable to determine the above-mentioned quantities along specific paths between high-symmetry \(\boldsymbol{q}\)-points when working with solids. In isotropic samples, such as liquids, it is, on the other hand, usually preferable to compute these functions with respect to \(q=|\boldsymbol{q}|\), i.e., a spherical average over wave vectors.

Note

In order to convolve the dynamic structure factor with form factors, you can use the functionality in the post-processing module. The procedure is demonstrated in this tutorial.

For a one-to-one comparison with experimental measurements, instrumental effects such as broadening and kinematic constraints must also be considered. These effects can be applied conveniently to Sample objects by following the calorine tutorial on inelastic neutron scattering.

Spectral energy density method

The spectral energy density (SED) is closely related to the phonon dispersion and thus also to the current correlations above. It is a measure of how the kinetic energy is partitioned over different wavelengths and frequencies in the system. The theoretical background for SED is quite simple and can be derived in several ways. Formally, the kinetic action \(W\) of the system can be written as

\[W = \int dt \sum_n \sum_i \frac{1}{2} m_i v_i(n, t)^2.\]

Here, \(n\) is the cell index and \(i\) is the basis index, which together uniquely identify each atom in an infinite crystal. By using Parseval’s theorem for both the time integration and space summation, the above is transformed to inverse spatial \(q\)-space and inverse temporal \(\omega\)-space

\[W = \int d\omega \sum_q \sum_i \frac{1}{2} m_i v_i(q, \omega)^2.\]

This is motivated by the Wiener–Khinchin theorem, and the SED is simply defined as the quantity

\[\text{SED}(q, \omega) = \sum_i \frac{1}{2} m_i v_i(q, \omega)^2.\]

For a more in-depth discussion see, e.g., [ThoTurIut10].

Note

To compute the SED you can use the compute_spectral_energy_density function. Its application is demonstrated in the SED tutorial.

Damped harmonic oscillator model

The correlation functions can be fitted to analytical expressions for a damped harmonic oscillator (DHO) in order to extract phonon frequencies and lifetimes. The governing differential equation for the displacement (or coordinate) of a DHO, \(x(t)\), is

\[\ddot{x}(t) + \Gamma \dot{x}(t) + \omega_0^2 x(t) = f(t).\]

Here, \(\omega_0\) is the natural angular frequency and \(\Gamma\) is the damping constant in units of reciprocal time. Note that we work in units where the mass \(m=1\). Below are the analytical solutions for the equations of motion of a DHO under the assumption that the stochastic force, \(f(t)\), is white noise.

The DHO can be either underdamped or overdamped. In the underdamped case, \(\omega_0\tau > 1\), where \(\tau = \frac{2}{\Gamma}\) is commonly referred to as the lifetime of the DHO. The position autocorrelation function (ACF) in the underdamped case is given by

\[F(t) = A \mathrm{e}^{- t/\tau} \left ( \cos{ \omega_e t} + \frac{1}{\tau\omega_e}\sin{ \omega_e t}\right ),\]

where the angular eigenfrequency of the DHO is \(\omega_e = \sqrt{\omega_0^2 - 1/\tau^2}\) and \(A\) is a constant amplitude \(A=F(t=0)\). The functional form is only valid for time lag \(t\geq 0\) and should be interpreted as even, so that \(t\rightarrow -t\) for negative times. Alternatively, replace \(t\) with \(|t|\) in the expressions for the ACFs.

In the overdamped case, \(\omega_0 \tau < 1\), the autocorrelation function is given by

\[F(t) = \frac{A}{\tau_\text{L} - \tau_\text{S}} \left( \tau_\text{L} \text{e}^{-t / \tau_\text{L}} - \tau_\text{S} \text{e}^{-t / \tau_\text{S}} \right),\]

where \(\tau_\text{S}\) and \(\tau_\text{L}\) correspond to the short and long relaxation times, respectively, defined as

\[\tau_\text{S,L} = \frac{\tau}{1 \pm \sqrt{1-(\omega_0\tau)^2}}.\]

The corresponding underdamped and overdamped velocity ACFs are related to the position ACFs via \(F_v(t) = -\ddot{F}(t)\):

\[\begin{split}F_v(t) &= A \omega_0^2 \mathrm{e}^{- t/\tau} \left ( \cos \omega_e t - \frac{1}{\tau\omega_e}\sin \omega_e t\right) \quad (\text{underdamped}) \\ F_v(t) &= \frac{A}{\tau_\text{L}-\tau_\text{S}} \left( \frac{1}{\tau_\text{S}} \text{e}^{-t / \tau_\text{S}} - \frac{1}{\tau_\text{L}} \text{e}^{-t / \tau_\text{L}} \right ) \quad (\text{overdamped})\end{split}\]

The analysis can also be carried out in the frequency domain by computing the power spectral density (PSD) functions, which are Fourier transforms of the corresponding ACFs. The position PSD of a DHO is given by

\[S(\omega) = A\frac{2\Gamma \omega_0^2}{(\omega^2 - \omega_0^2)^2 + (\Gamma\omega)^2},\]

which exhibits a peak at \(\omega_\text{max} = \sqrt{\omega_0^2 - \Gamma^2/2}\).

The PSD of the velocity, \(S_v(\omega) = \omega^2 S(\omega)\), is

\[S_v(\omega) = A\frac{2\Gamma \omega_0^2 \omega^2}{(\omega^2 - \omega_0^2)^2 + (\Gamma\omega)^2}.\]

The velocity PSD has a peak at \(\omega_\text{max} = \omega_0\). The amplitude \(A\) can be related to the variance (strength) \(\sigma^2=\langle f^2 \rangle\) of the noise by \(2\Gamma\omega_0^2 A=\sigma^2\). Taking the strength to be equal to \(\sigma^2=2\Gamma k_B T\), the mean squared displacement becomes \(\langle x^2 \rangle=\frac{k_B T}{\omega_0^2}\), and the mean squared velocity is given by \(\langle v^2 \rangle=k_B T\).

The position and velocity ACF and PSD solutions for the DHO with natural angular frequency \(\omega_0 = 1\) for a few different values of the damping constant \(\Gamma\). The red line corresponds to the critical damping. For the position ACF, a critically damped oscillator goes to zero without any oscillation. When the damping is larger than \(\sqrt{2}\), the peak for the corresponding PSD remains at zero.

Note

Functions for fitting DHO parameters can be found in the DHO module. The application of these functions is illustrated in the peak fitting tutorial.

Mode projection

The primitive cell metric \(C_{ax}\) is related to the supercell metric \(S_{ax}\) via the repetition matrix \(P\) as

\[P_{ab}C_{bx}=S_{ax}.\]

where the lattice vectors are rows. The inverse lattice is given by

\[\widetilde{C}_{ay} C_{ax}=\delta_{yx}\]

i.e. \(C^{-T}\), the inverse transpose of the primitive cell metric. The commensurate reduced q-points are defined such that

\[r_a \widetilde{C}_{ax} = n_a \widetilde{S}_{ax}\]

for some integer vector \(n_a\) and \(0 \leq r_a < 1\). These properties can be accessed by mp.supercell.P, the cell metrics as mp.primitive.cell and mp.supercell.cell, the inverse cells as mp.primitive.inv_cell and mp.supercell.inv_cell, and the commensurate \(\boldsymbol{q}\)-points as mp.q_reduced. In the following, a sum over \(\boldsymbol{q}\)-points \(\sum_\boldsymbol{q}\) is understood such that \(\boldsymbol{q}\) is taken from this list indexed by \(k\). Important to note is that the \(\boldsymbol{q}\)-points are not taken from the Wigner-Seitz cell but from the interior of the inverse cell, i.e. \(-\boldsymbol{q}\) is mapped back into the inverse cell.

Each unit cell in the supercell, indexed by an integer \(l\), contains some basis atoms indexed by \(s\). The displacement of an atom is thus denoted by \(u_{ls\alpha}\) where \(\alpha\) indexes the Cartesian directions. The lattice site and basis index for each atom (indexed by \(i\)) in the supercell can be accessed by mp.supercell.offsets and mp.supercell.indices respectively, and these lists are indexed by \(i\). The supercell Cartesian equilibrium position is denoted by \(r_{ls}\). In the following, a sum over lattice sites or basis indices is taken from the unique elements in these lists.

The ModeProjector object performs the following mapping from displacements \(u\) to mode coordinates \(Q\)

\[Q_{kb} = \frac{1}{\sqrt{N_\text{unit}}} \sum_{ls\alpha} \sqrt{m_s} u_{ls\alpha} \overline{W}_{kbs\alpha} \mathrm{e}^{-iq_k \cdot r_{ls}}\]

and back

\[u_{ls\alpha} = \frac{1}{\sqrt{N_\text{unit}}} \sum_{kb} \frac{1}{\sqrt{m_s}}Q_{kb} W_{kbs\alpha} \mathrm{e}^{-iq_k \cdot r_{ls}}\]

where \(W\) are the polarization vectors and \(r_{ls}\) is the scaled position if \(\boldsymbol{q}\) is taken from the reduced set or the Cartesian position if \(\boldsymbol{q}\) is taken from the Cartesian set.

The polarization vectors are constructed from the eigenvectors and eigenvalues of the dynamical matrix

\[D_{ss'}^{\alpha\beta}(q) = \sum_l \frac{\Phi_{ss'}^{\alpha\beta}(0, l)}{\sqrt{m_s m_{s'}}} \mathrm{e}^{iq (r_s(l)-r_{s'}(0))} = \sum_b W^\alpha_{sb}(q) \omega_b^2(q) \overline{W}_{bs'}^\beta(q)\]

It is assumed that the force constants obey lattice translation symmetry. The polarization vectors \(W\) can be accessed by mp.polarizations and the frequencies as mp.omegas.

The eigenmodes \(X\) in the supercell are defined so that

\[Q_{kb} = X_{kbi\alpha} u_{i\alpha}\]

and

\[u_{i\alpha} = \overline{X}_{kbi\alpha} Q_{kb} / m_i\]

For momenta we have \(v=XP/m\) and \(P=\overline{X}v\), while for force we have \(f=Xf\) and \(F=\overline{X}f/m\).

Note

You can find more information in the comprehensive documentation of the functionality available for mode projection analysis. The mode projection tutorial illustrates how this functionality can be used in practice.