Theoretical background

Below follows 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 correlations 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^{-iwt} \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 a longitudinal and transverse part 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 the other. The current correlation functions can be thought of as 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 correlations functions, which enables 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 correlaion 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 runs over all A atoms. Here, \(N_\mathrm{A}\) refers to number of atoms of type A, \(N_\mathrm{B}\) refers to 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 two in the final expression is due to considering both the A-B and B-A terms. These normalization choices means 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 weighting them with relevant atomic form factors, masses or charges etc.

Relation to scattering experiments

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

\[F_\mathrm{coh}(\boldsymbol{q},t)=\frac{1}{N}\sum _i ^N \sum _j ^Nw_iw_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 lenghts 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 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 ^Nf_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 quantites along specific paths between high-symmetry \(\boldsymbol{q}\)-points when working with solids. In isotropic samples, such as for example 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 convolute the dynamic structure factor with form factors, you can use the functionality in the post-processing module. The procedure is demonstrated in this tutorial.

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 follow 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 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 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 fixates 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 as rows. The inverse lattice is given by

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

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

\[r_a \tilde{C}_{ax} = n_a \tilde{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 and the inverse cells as mp.prim.inv_cell and mp.supercell.cell and the commensurate q-points as mp.q_reduced. In the following, a sum over q-points \(\sum_q\) it is understood that \(q\) is taken from this list indexed by \(k\). Important to note is that the q-points are not taken from the Wigner-Seitz cell but from the interior of the inverse cell, i.e. \(-q\) is mapped back into the inverse cell.

To each unit cell, indexed by an integer \(l\), in the supercell belongs some basis atoms indexed by \(s\). The displacement of an atom is thus denoted by \(u_{ls\alpha}\) where \(\alpha\) index the cartesian directions. The lattice site and basis index for each atom (indexed by \(i\)) in the supercellcan 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 are 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} \bar{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 q is taken from the reduced set or the cartesian position if 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) \bar{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 is defined so that

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

and

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

for momenta and force we have \(v=XP/m\) and \(P=\bar{X}v\) as well as \(f=Xf\) and \(F=\bar{X}f/m\), respectively.

Note

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