Source-linked AI summary
DynaPhoPy: A code for extracting phonon quasiparticles from molecular dynamics simulations
Abel Carreras, Atsushi Togo, Isao Tanaka
TL;DR
The paper develops DynaPhoPy to address temperature-dependent anharmonic phonon properties that require analysis beyond the harmonic approximation. It uses normal-mode decomposition of molecular-dynamics trajectories to extract quasiparticle frequencies and linewidths, then tests first-principles and Tersoff simulations for silicon against Raman data. Frequency shifts agree well with experiment for both approaches, while linewidths show discrepancies between the approaches at X and L points.
Problem
The harmonic approximation does not describe temperature-dependent phonon frequencies, motivating methods that analyze anharmonic phonon behavior.
Method
DynaPhoPy projects molecular-dynamics atomic velocities onto phonon eigenvectors and fits their power spectra to extract quasiparticle frequencies and linewidths.
Results
Frequency shifts agree well with Raman experiments using both first-principles and Tersoff approaches, while linewidths disagree between approaches at X and L points.
Takeaways & Limitations
The methodology provides a systematic way to calculate microscopic anharmonic phonon properties from molecular-dynamics simulations and assess empirical potentials.
Takeaways & Limitations
Quantum lattice-dynamics effects are not considered, and first-principles molecular-dynamics simulations remain computationally expensive.
Abstract
from arXiv · showhide
We have developed a computational code, DynaPhoPy, that allow us to extract the microscopic anharmonic phonon properties from molecular dynamics (MD) simulations using the normal-mode-decomposition technique as presented by Sun et al. [T. Sun, D. Zhang, R. Wentzcovitch, 2014]. Using this code we calculated the quasiparticle phonon frequencies and linewidths of crystalline silicon at different temperatures using both of first-principles and the Tersoff empirical potential approaches. In this work we show the dependence of these properties on the temperature using both approaches and compare them with reported experimental data obtained by Raman spectroscopy [M. Balkanski, R. Wallis, E. Haro, 1983 and R. Tsu, J. G. Hernandez, 1982].
PROGRAM SUMMARY
DynaPhoPy is a Python and C code for extracting phonon quasiparticle properties from molecular-dynamics simulations. It supports Fourier-transform and maximum-entropy power-spectrum methods, but does not include quantum lattice-dynamics effects.
- The code is implemented in Python and C for UNIX/OSX systems and can run on PCs or cluster computers.
- DynaPhoPy extracts phonon quasiparticles from molecular-dynamics simulations using normal-mode decomposition.
- Quantum effects of lattice dynamics are not considered.
- Two power-spectrum methods are implemented: Fourier transformation and maximum entropy.Their stated computational scalings are O(N^2) and O(N × M), respectively.
1. Introduction
The paper addresses temperature-dependent anharmonic phonon behavior, which harmonic lattice dynamics cannot describe. It develops a molecular-dynamics analysis based on normal-mode decomposition and evaluates it for silicon using first-principles and Tersoff approaches against Raman measurements.
- The harmonic approximation does not describe the temperature dependence of phonon frequencies, motivating explicit treatment of crystal anharmonicity.
- Anharmonic phonon-phonon interactions produce frequency shifts and lifetimes represented by phonon quasiparticles observable through Raman spectroscopy and neutron scattering.
- The method fits mass-weighted velocity power spectra from molecular dynamics to spectral-function shapes to obtain quasiparticle frequencies and linewidths.Atomic velocities are projected onto phonon eigenvectors using normal-mode decomposition.
- For silicon, first-principles and Tersoff simulations reproduce experimental phonon frequency shifts well, while empirical potentials enable longer simulations and larger supercells.
2. Anharmonic model
The anharmonic model represents atomic motion through phonon quasiparticle velocities rather than purely harmonic normal-mode displacements. These quantities are defined from lattice positions, masses, eigenvectors, wave vectors, and phonon branches.
- Atomic positions are decomposed into equilibrium positions and displacements at each lattice point and atomic index.
- Within the harmonic approximation, atomic displacements are expressed as superpositions of phonon normal modes.
- The normal modes are obtained from the dynamical-matrix eigenvalue problem, with masses, eigenvectors, wave vectors, branches, and phonon coordinates defining the expansion.
- The anharmonic model defines atomic velocity using phonon quasiparticle velocities associated with the phonon modes.
- A phonon quasiparticle is considered well modeled when its velocity power spectrum has a well-defined spectral shape.
3. Methodology
The methodology projects molecular-dynamics atomic velocities onto phonon eigenvectors, constructs wave-vector- and mode-resolved power spectra, and fits them with Lorentzian functions. Peak positions provide quasiparticle frequencies, while FWHM values provide linewidths.
- Atomic velocities are projected onto phonon eigenvectors to obtain the phonon-quasiparticle velocity v_qs(t) from molecular-dynamics trajectories.
- A Lorentzian form is used to fit G_qs(ω) when phonon frequency shifts and linewidths are sufficiently small.
- The quasiparticle phonon frequency is the fitted peak position, and the linewidth is the full width at half maximum.
- The full power spectrum and wave-vector-resolved spectra are related through orthonormal phonon eigenvectors and their total power.The total power is related to the system’s vibrational kinetic energy.
- The projected-velocity normalization uses the atomic velocity and average vibrational kinetic energy, with the classical-limit expression involving the number of atoms per lattice point and Boltzmann’s constant.
- The power spectrum is computed from the quasiparticle velocity and its velocity-autocorrelation function.
4. Software overview
DynaPhoPy extracts quasiparticle phonon frequencies and linewidths from molecular-dynamics trajectories through normal-mode decomposition. Its workflow projects velocities onto wave vectors and phonon eigenvectors, calculates power spectra, and supports interpolation between commensurate and incommensurate q-points.
- DynaPhoPy extracts quasiparticle phonon frequencies and linewidths from molecular-dynamics trajectories using normal-mode decomposition.
- Mass-weighted velocities are projected first onto a wave vector and then onto a phonon eigenvector to obtain mode-resolved velocities.
- The mode-resolved velocity power spectrum is fitted to obtain quasiparticle phonon properties.
- Quasiparticle frequencies calculated at commensurate q-points are interpolated to incommensurate q-points using a crudely renormalized force-constants matrix.
- The renormalized force constants can support thermodynamic-property calculations, but this work focuses on frequency shifts and linewidths instead.
5. Computational details
The study benchmarks DynaPhoPy on crystalline silicon using first-principles and Tersoff-potential molecular dynamics. It compares computational approaches, supercell considerations, and power-spectrum methods for temperature-dependent quasiparticle properties.
- Crystalline silicon was selected as a benchmark because its cubic structure and extensive theoretical and experimental literature provide substantial reference information.
- The simulations computed silicon quasiparticle phonon frequencies and linewidths versus temperature using first-principles and Tersoff empirical-potential approaches.
- The computational setups were documented separately for first-principles VASP simulations and Tersoff-potential LAMMPS simulations.
- First-principles calculations served as the higher-prediction-level reference, while the less-demanding empirical potential enabled analysis of MD supercell-size effects.
- The study used maximum entropy for the full power spectrum and Fourier transform for frequency shifts and linewidths.
6. Results and discussion
The results use harmonic references and MD-derived power spectra to extract temperature-dependent quasiparticle frequencies and linewidths in silicon. First-principles and Tersoff calculations generally agree qualitatively, while comparisons with experiment and supercell tests establish the scope of the results.
- Power spectra: Harmonic band structures and DOS provide reference shapes without temperature effects, while Gq(ω) assigns MD spectral peaks to quasiparticle modes at commensurate wave vectors.The analysis examines projected spectra at Γ, X, and L and separates modes using phonon eigenvectors.
- Power spectra: As temperature increases, full-spectrum peaks broaden and shift to lower frequencies, reflecting anharmonicity; overlapping peaks prevent individual quasiparticle analysis in G(ω).G(ω) contains the superposition of modes represented by the 64-atom supercell, so its shape differs from the interpolated DOS.
- Power spectra: Average kinetic energies agree excellently with the classical limit at all temperatures for both VASP and LAMMPS, indicating that the ME method reproduces the full power spectrum.The degree-of-freedom count differs between codes because VASP constrains the center of mass while LAMMPS does not.
- Power spectra: Projected spectra show one LTO peak at Γ, three peaks at X, and four peaks at L, consistent with the corresponding phonon branches.The acoustic Γ mode is absent because it represents lattice translation restricted in the VASP simulations; mode-resolved spectra match the peaks in Gq(ω).
- Frequency shifts: Frequency shifts from first-principles and Tersoff calculations show similar temperature tendencies, but Tersoff systematically underestimates the experimental-contribution shift, especially for LTO at Γ and TA at X and L.Some Tersoff optical-mode shifts also show unexpected high-temperature curvature absent from first-principles results.
- Frequency shifts and linewidths: The calculated total frequency shifts generally agree with Raman data, with excellent agreement at Γ, while 64-atom supercells sufficiently describe frequency shifts and linewidths.Linewidths agree with Raman results at Γ; first-principles linewidths for TO and LA at X and L become much larger than Tersoff values at high temperatures, and 64- versus 512-atom differences are weak.
7. Conclusions
DynaPhoPy implements normal-mode decomposition to extract anharmonic phonon properties from MD simulations. Tests on silicon show agreement with Raman data for first-principles calculations, while Tersoff linewidths differ at X and L and first-principles MD remains expensive.
- The implemented normal-mode-decomposition methodology calculates microscopic anharmonic phonon properties from molecular-dynamics simulations.
- Frequency shifts and linewidths of crystalline silicon were computed at Γ, X, and L using first-principles and Tersoff approaches.
- First-principles frequency shifts and linewidths agree well with Raman spectroscopy experiments.
- Tersoff frequency shifts agree well with experiments and first-principles calculations, but linewidths show discrepancies at X and L.
- First-principles molecular-dynamics simulations remain computationally expensive.
Appendix A. The discrete Fourier transform
The discrete Fourier-transform workflow converts discretized MD velocities into power spectra with controlled frequency limits and resolution. Segment averaging improves precision, while FFT implementations reduce computation time.
- The discrete Fourier transform computes power spectra from velocities sampled at constant time intervals.
- The power spectrum shares the discrete transform’s frequency range and resolution, with resolution determined by the sampling and number of samples.
- At constant sampling interval, increasing the number of samples improves frequency resolution.
- DynaPhoPy chooses the Nyquist frequency and resolution to cover the target range and resolve quasiparticle phonon peaks.
- The full MD trajectory is divided into shorter segments whose power spectra are averaged to improve measurement precision.
- The fast Fourier transform substantially reduces computation time compared with direct discrete-transform evaluation.
Appendix B. The maximum entropy method
The maximum-entropy method estimates smoother power spectra through an autoregressive model and recursively optimized coefficients. Burg’s method minimizes forward and backward prediction discrepancies using Levinson–Durbin recursion.
- The maximum-entropy method provides a smoother power-spectrum estimate at lower computational cost than the discrete Fourier transform, but requires more samples for good estimation.
- The method assumes that atomic trajectories follow an autoregressive model whose coefficients determine the power spectrum.
- Linear prediction uses previous values for forward prediction and following values for backward prediction, producing corresponding discrepancy vectors and coefficients.
- Burg’s method minimizes the sum of forward and backward discrepancies through Levinson–Durbin recursion.
- Reflection coefficients are optimized iteratively as the autoregressive coefficient set is updated.
Appendix C. Supercell size effect on first-principles MD calculations
Larger first-principles MD supercells produce qualitatively similar temperature trends but limited-precision linewidth results. Their more linear temperature dependence suggests that curvature in smaller-supercell results may reflect finite-size effects.
- The 512-atom first-principles supercell results should be treated qualitatively because computational demands prevent the precision achieved in the main calculations.
- 512-atom linewidths at Γ are slightly overestimated relative to experiment while retaining the same temperature tendency.
- The 512-atom linewidths have a more linear temperature dependence than 64-atom results, consistent with the general linear experimental tendency.
- The curvature of linewidths at high temperatures may result from limited supercell size, an effect also seen for the Tersoff Γ-point linewidth.