Efficient perturbation-tracking method for directly probing the spectral phonon properties from molecular dynamics simulations

Existing methods for directly extracting the spectral phonon properties from molecular dynamics (MD) simulations, like the normal mode analysis (NMA) and spectral energy density analysis, all require a very long simulation time to produce reliable results with good convergence. So far, these methods are mainly applied in studies using small systems and with empirical potentials, as the heavy computational load has greatly hindered their further applications. Here we propose a perturbation-tracking (PT) method for directly probing the mode-wise phonon anharmonic frequencies and lifetimes. We show that results obtained from our method are in excellent agreement with those from the conventional NMA approach, using Si as the model material system. Comparing with the NMA approach, the PT method offers a greater accuracy and significant improvement of efficiency. It takes an average of two orders of magnitude and up to three orders of magnitude less simulation time to obtain the same lifetime result of a phonon mode with intermediate to high accuracy. Meanwhile, our method preserves all the dynamics of probed phonon mode from a particular state, which means it is capable of studying the transient thermal transport processes in a nonequilibrium system. Besides the exceptional efficiency, our method also comes with freedom to choose to probe only those modes of interest. This makes it ideal for use with large systems and in computationally demanding applications, such as ab initio MD simulations. Moreover, the PT method we propose here is very straightforward and easy to implement.

Figure 1

Comparison of frequency shifts of Si phonons at 600 K in the [001] direction from four different polarization branches: (a) TA, (b) LA, (c) LO, and (d) TO. The NMA results are plotted as blue squares and those from the PT method are plotted as blue crosses (these correspond to the left axis). The dispersion relations are plotted as red curves (use the right axis).

Figure 2

Convergence of probed anharmonic frequencies with simulation run length using the PT method. Two phonon modes are sampled for each different polarization branch: (a) TA, (b) LA, (c) LO, and (d) TO.

Figure 3

Examples of the decay profiles of three LA phonon modes from the PT method. The series of discrete symbols in different colors are the total probe energy of each mode normalized by the respective initial amplitude at $t=0$ and each black line represents the corresponding best exponential fit.

Figure 4

Comparison of phonon lifetimes of Si in the [100] direction at 600 K obtained from the NMA (orange squares) and the PT method (blue circles with error bars). Results from four different polarization branches are presented: (a) TA, (b) LA, (c) LO, and (d) TO.

Figure 5

Comparison of convergence of the lifetime result from the NMA approach (left) and the PT method (right) with simulation run length. Example of a Si LA mode of 12.05 THz. The blue curve with symbols in each plot is the average from multiple runs and the two black curves form an uncertainty band with its upper and lower bound being one standard deviation away from the average.

Figure 6

Variation of the RE versus simulation run length calculated using Eq. (33) for Si phonon modes. Results from lifetimes obtained by the NMA approach are plotted in the left panel and those by the PT method are plotted in the right panel. (a) RE values for sampled TA and LA modes. (b) RE values for sampled LO and TO modes.