phq: a Fortran code to compute phonon quasiparticle properties and dispersions

TL;DR

This paper introduces a Fortran code that calculates anharmonic phonon properties and free energy of solids from first-principles, using a phonon quasiparticle approach based on molecular dynamics simulations.

Contribution

The code uniquely combines phonon quasiparticle extraction from MD with harmonic force constant construction to compute anharmonic phonon dispersions and free energy.

Findings

Successfully applied to silicon in diamond structure

Accurately computes temperature-dependent anharmonic phonons

Provides a practical tool for first-principles anharmonic calculations

Abstract

We here introduce a Fortran code that computes anharmonic free energy of solids from first-principles based on our phonon quasiparticle approach. In this code, phonon quasiparticle properties, i.e., renormalized phonon frequencies and lifetimes, are extracted from mode-projected velocity auto-correlation functions (VAF) of modes sampled by molecular dynamics (MD) trajectories. Using renormalized frequencies as input, the code next constructs an effective harmonic force constant matrix to calculate anharmonic phonon dispersions over the whole Brillouin zone and thus the anharmonic free energy in the thermodynamic limit ($N \rightarrow \infty$). We illustrate the use of this code to compute ab initio temperature-dependent anharmonic phonons of Si in the diamond structure.

Tables

Table 1. Table 1: A description of subroutines in main.f90

Subroutine name	Description
control	General control of the program.
read_scf	Read in structural information of the primitive cell.
read_md	Read in MD information and compute atomic displacements from their static positions.
properties	Compute the heat capacity via heat fluctuation. If no information other than atomic coordinates is provided in the md.out file, comment out this subroutine in the control subroutine.
read_dyn	Read in harmonic phonon information, compute eigenmodes of the supercell and the primitive cell at different q-point and check the orthogonality.
harmonic_matrix	Compute the harmonic force constant matrix from harmonic frequencies and harmonic eigenmodes.
harmonic force	Compute the interatomic harmonic forces.
frequency	Compute the atomic velocities and perform velocity eigenmode projection, according to Eq. (2)(3).
equipartition	Compute the ionic kinetic energy from equipartition theorem.
correlation	Compute the VAFs, accourding to Eq. (1).
correlation_fit	Fit the VAFs to Eq. (5) and obtain quasiparticle properties [2].
corr_fourier	Perform FT of VAFs according to Eq. (4) and obtain quasiparticle properties.
corr_fit_fourier	Perform FT of the fitted curves. Avoid calling this subroutine if shorter running time is desired.
maximum_entropy	Perform MEM calculation and obtain the quasiparticle properties.
linear_response	Obtain MEM power spectrum via linear prediction method.
lorentzian	Obtain quasiparticle properties by performing a fitting of the MEM power spectrum to a Lorentzian function.
gamma_matrix	Compute the effective harmonic force constant matrix from renormalized frequencies and harmonic eigenmodes.
gamma_force	Compute the interatomic anharmonic forces.
dynamics_matrix	Compute the harmonic dynamical matrices from harmonic frequencies and harmonic eigenmodes. Comment out this subroutine if the effective harmonic dynamical matrices are desired.
dynamics_matrix_md	Compute the effective harmonic dynamical matrices from renormalized frequencies and harmonic eigenmodes, according to Eq. (6)(7)(8).
write_dym	Write down the effective harmonic dynamical matrices or the harmonic dynamical matrices with the same format of ph.x executable’s output, which can be read in by q2r.x executable in the Quantum ESPRESSO [16] suite to perform Eq. (9).

Table 2. Table 2: Parameters used in the input

Parameter name	Description
dt	MD time step in atomic unit. If 1 fs time step is used, simply enter dt = 20.67055273.
step_md_use	Number of MD steps needed for the calculation of mode-projected VAFs.
correlation_time	Desired decay time for VAFs in unit of dt.
pole	Parameter used in the MEM. It is used to filter high-frequency components. Usually a reasonable range of pole (200 $\sim$ 2000) can yield smooth spectrum.
supercell	Supercell used in MD related to the primitive cell. 3 integers are required to enter.
temperature	MD temperature in unit of Kelvin.
method	Enter one of the integer number 0, 1 or 2 to select renormalized frequencies obtained by which method to compute the effective harmonic force constants and dynamical matrices. 0 represents the fitting approach (Eq. (5)). 1 represents the Fourier transformation. 2 represents the maximum entropy method.

Table 3. Table 3: Parameters used in the scf.out

Parameter name	Description
ntype	Number of different elements in the primitive cell.
natom	Number of atoms in the primitive cell.
mass	Atomic mass of each element following the symbol of the element.
lattice_parameter	Scale of lattice vectors in unit of Bohr radius.
cell_parameters	Lattice vectors in cartesian cooradinate in unit of lattice_parameter.
atomic_positions	Atomic positions of atoms in the primitive cell in reduced coordinates.

Table 4. Table 4: Parameters used in the dyn.out

Parameter name	Description
q	q-point in cartesian coordinates in unit of 2$\pi$/lattice_parameter. This parameter is required.
freq	Harmonic phonon frequencies. This parameter in unit of THz is required as program input as well as the following six columns of eigenvectors. Eigenvector of each atom has x, y and z components and each component has a real part and an imaginary part. The order of atoms should be the same as that of the atoms entered in atomic_positions in scf.out.
Dielectric Tensor	Dielectric Tensor, which is a $3\times 3$ matrix of real numbers. This parameter is optional, depending on whether LO-TO splitting needs to be considered in the system.
Effective Charges	Effective Charges, which are $3\times 3$ matrices of real numbers for each of the atom. This parameter is optional, depending on whether LO-TO splitting needs to be considered in the system.

Table 5. Table 5: Parameters used in the md.out

Parameter name	Description
total_step	Total actual MD steps. In practice, recorded MD steps should be configurations after reaching thermal equilibrium and therefore less than this number.
atomic_positions	Initial atomic positions in reduced coordinates of the supercell lattice vectors.
md_step	Integer numbers of recorded MD steps, followed by the instantaneous atomic_md_positions in the MD simulation.
atomic_md_positions	Atomic positions in reduced coordinates of the supercell lattice vectors during the MD simulation.

Table 6. Table 6: Output files of the phq program

Output	Description
corr.vaf	VAF of each eigenmode of the supercell. The VAF of each eigenmode lies on each column following the correlation time on the first column.
corr_fit.vaf	Fitted exponentially decaying cosine functions (Eq. (5)) of correlation functions according to Green’s function from molecular dynamics [2]. The fitted curve of VAF of each eigenmode of the supercell lies on each column following the correlation time on the first column.
corr_fit_fourier.vaf	FT of the fitted VAF curves. The FT of fitted VAF curve of each eigenmode of the supercell lies on each column following their common frequency range lies on the first column.
corr_fourier.vaf	FT of VAFs. The FT of VAF of each eigenmode of the supercell following its frequency range lie on every two columns.
dynamical_matrix_md.mat	Effective harmonic dynamical matrices. It is an assemble of all the effective harmonic dynamical matrices from dynmatmd1, …, dynmatmdnq.
dynmatmd0	q-point list with the same format of ph.x executable’s output in the Quantum ESPRESSO suite. They are in cartesian coordinates in unit of $2\pi$/lattice_parameter.
dynmatmd1, …, dynmatmdnq	Renormalized phonon information with the same format of ph.x executable’s output in the Quantum ESPRESSO suite. nq is the total number of the q-points. They share the same set of parameters of dyn.out.
frequency.freq	Harmonic, fitted renormalized, FT renormalized and MEM renormalized phonon frequencies lie on the second, third, fourth and fifth column separately while the first column labels different eigenmodes of the supercell.
gamma_matrix.mat	Effective harmonic force constant matrix, which are $3\times 3$ matrices of real numbers between each pair of atoms in the supercell.
harmonic_matrix.mat	Harmonic force constant matrix, which are $3\times 3$ matrices of real numbers between each pair of atoms in the supercell.
msd.out	Mean square displacement of each atom in the supercell.
tau_fit.tau	Phonon quasiparticle’s lifetime of each eigenmode of the supercell obtained by fitting the VAFs to Eq. (5).
tau_fourier.tau	Phonon quasiparticle’s lifetime extracted from the FT of VAF of each eigenmode of the supercell.
tau_mem.tau	Phonon quasiparticle’s lifetime extracted from the MEM power spectrum of each eigenmode of the supercell.
vector.out	Eigenmodes of the supercell.
vector_q.out	Eigenmodes of the primitive cell at each q-point.

Equations

$$⟨V_{\mathbf{q}s}(0)\cdot V_{\mathbf{q}s}(t)⟩=\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }{V}_{\mathbf{q}s}^{\ast }(t^{\prime })V_{\mathbf{q}s}(t^{\prime }+t)d{t}^{\prime }.$$

$$V_{\mathbf{q}s}(t)=\sum _{i=1}^{N}V(t)\cdot e^{i\mathbf{q}\cdot {\mathbf{r}}_i}\cdot {\hat{\mathbf{e}}}_{\mathbf{q}s}$$

$$V(t)=V(\sqrt{M_1}{\mathbf{v}}_1(t),...,\sqrt{M_N}{\mathbf{v}}_N(t))$$

$$G_{\mathbf{q}s}(\omega )={\left|{\int }_0^{\infty }⟨V_{\mathbf{q}s}(0)\cdot V_{\mathbf{q}s}(t)⟩e^{i\omega t}dt\right|}^2.$$

$$A_{\mathbf{q}s}\cos ({\tilde{\omega }}_{\mathbf{q}s}t)e^{-t/(2{\tau }_{\mathbf{q}s})}$$

$$\tilde{D}(\mathbf{q})={\hat{\mathbf{e}}}_{\mathbf{q}}{\Omega }_{\mathbf{q}}{\hat{\mathbf{e}}}_{\mathbf{q}}^{†},$$

$${\Omega }_{\mathbf{q}}=\mathrm{diag}[{\tilde{\omega }}_{\mathbf{q}1}^2,{\tilde{\omega }}_{\mathbf{q}2}^2,...,{\tilde{\omega }}_{\mathbf{q}3N}^2].$$

$${\hat{\mathbf{e}}}_{\mathbf{q}}=[{\hat{\mathbf{e}}}_{\mathbf{q}1},{\hat{\mathbf{e}}}_{\mathbf{q}2},...,{\hat{\mathbf{e}}}_{\mathbf{q}3N}].$$

$$\tilde{D}({\mathbf{q}}^{\prime })=\sum _{\mathbf{r}}\tilde{\Phi }(\mathbf{r})\cdot e^{-i{\mathbf{q}}^{\prime }\cdot \mathbf{r}}$$

$${\tilde{\omega }}_{\mathbf{q}s}(T){|}_{P=0}={\tilde{\omega }}_{\mathbf{q}s}(T){|}_{V=V_{eq}}+{\omega }_{\mathbf{q}s}[\exp (-{\gamma }_{\mathbf{q}s}\bar{\alpha }T)-1],$$

$$\bar{\alpha }(T)=\frac{1}{T}{\int }_0^T\alpha (T^{\prime })d{T}^{\prime }.$$

$${\gamma }_{\mathbf{q}s}=-\frac{d\ln {\omega }_{\mathbf{q}s}(T)}{d\ln V(T)}.$$

$${\gamma }_{\mathbf{q}s}=-\frac{d\ln {\omega }_{\mathbf{q}s}(0)}{d\ln V(0)},$$

$$\alpha (T)=\frac{1}{V(T)}\frac{dV}{dT}.$$

$$\bar{\alpha }(T)=\frac{1}{T}\ln \frac{V(T)}{V(0)}.$$

$$S=k_B\sum _{\mathbf{q}s}[(n_{\mathbf{q}s}+1)\ln (n_{\mathbf{q}s}+1)-n_{\mathbf{q}s}\ln n_{\mathbf{q}s}],$$

$$F(T)=E_0+\frac{1}{2}\sum _{\mathbf{q}s}\hslash {\omega }_{\mathbf{q}s}-{\int }_0^{T}S(T^{\prime })d{T}^{\prime },$$

$$C(\omega )={\int }_0^{\infty }c(t)e^{i\omega t}dt$$

$$D_k=\sum _{j=0}^{N-1}{c}_j{\omega }_j{e}^{i2\pi jk/N}k=0,...,N-1,$$

$$P(f_k)=\frac{1}{N{\sum }_{j=0}^{N-1}{\omega }_j^2}{\left|D_k\right|}^2.$$

$$f_k=\frac{k}{N\Delta }=f_c\frac{k}{N}k=0,...,N-1,$$

$${\omega }_j=1j=0,...,N-1,$$

$${\omega }_j=\sin {\left[\frac{\pi (N-j)}{2N}\right]}^{2}j=0,...,N-1.$$

$$P(f)=\frac{a_0}{{\left|1+{\sum }_{j=1}^N{a}_j{e}^{i2\pi f\Delta j}\right|}^2},$$

$$y_n=\sum _{j=1}^M{a}_j{y}_{n-j}+x_n$$

$$\mathcal{H}(f)=\frac{1}{1-{\sum }_{j=1}^M{d}_j{e}^{i2\pi f\Delta j}}$$

$$a_0=\text{xms},a_j=-\text{d}(j)j=1,...,N.$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program title: phq

Licensing provisions: GPLv3

Programming language: Fortran 90

Operating system: Linux, macOS, Windows

*Nature of problem:

*Accurate free energy calculations deliver predictive thermodynamic properties of solids. Although the quasi-harmonic approximation (QHA) has been widely used for various materials, its validity at very high temperatures is not guaranteed because of intrinsic anharmonic effects. The QHA cannot be used even at low temperatures for strongly anharmonic systems, e.g., those that are stabilized by anharmonic fluctuations. When anharmonicity is non-negligible, phonon frequencies display a pronounced dependence on temperature, which impacts thermodynamic properties. Therefore, the calculation of anharmonic phonon dispersions is critical for an accurate estimation of the free energy.

*Solution method:

*We calculate phonon quasiparticle properties, i.e., renormalized phonon frequencies and lifetimes, by combing molecular dynamics (MD) trajectories and harmonic normal modes (HNM). MD velocities are projected into HNMs and the mode-projected velocity auto-correlation functions (VAF) are then calculated. Quasiparticle properties are then obtained by analyzing the VAFs, if they are well defined. Using the renormalized phonon frequencies and HNMs we build an effective (temperature-dependent) harmonic force constant matrix from which the anharmonic phonon dispersions can be obtained.

*Unusual features:

*(i) Because MD captures anharmonic effects to infinite order, this approach to compute renormalized phonon frequencies and lifetimes is expected to be more appropriate than lowest-order perturbation theory for strongly anharmonic systems or for weakly anharmonic systems at very high temperatures. (ii) By characterizing phonon quasiparticles, this approach provides a validity check of the phonon concept, a fundamental procedure missing in other approaches to compute anharmonic dispersions, e.g., the self-consistent phonon method.

1 Introduction

The phonon gas model (PGM) [1, 2, 3, 4] provides an excellent paradigm for calculating thermodynamics properties of materials. For weakly anharmonic systems, anharmonic effects reflecting the intrinsic temperature dependence of phonon frequencies are negligible. Thus, the quasiharmonic approximation (QHA) [5, 6] offers predictive free energies and thermodynamics properties. Note that the QHA accounts for the extrinsic temperature dependence of the phonon frequencies caused by volume changes (thermal expansion). However, in strongly anharmonic systems, intrinsic anharmonic effects are key to structural stabilization and must be properly accounted for. Therefore, anharmonic phonon dispersions are required for estimations of free energy and related thermodynamic properties.

Here we introduce a code of a hybrid approach which fully accounts for anharmonic effects. In this approach, we describe a phonon quasiparticle of normal mode $({\bf q},s)$ by the mode-projected velocity autocorrelation functions (VAF) [7],

[TABLE]

where,

[TABLE]

is the $s$-th branch of q-projected velocity. ${\bf\hat{e}}_{{\bf q}s}$ is the harmonic phonon polarization vector of mode (${\bf q},s$).

[TABLE]

is the weighted velocity with $3N$ components. ${\bf v}_{i}(t)(i=1,...,N)$ are atomic velocities calculated from MD trajectories of an $N$-atom supercell. $M_{i}$ is the atomic mass of the $i^{th}$ atom in the supercell.

For a well-defined phonon quasiparticle, its power spectrum,

[TABLE]

should have a Lorentzian line shape with a peak at $\widetilde{\omega}_{{\bf q}s}$ and a linewidth of 1/(2$\tau_{{\bf q}s}$), where $\tau_{{\bf q}s}$ is the lifetime of mode (${\bf q},s$) in the relaxation time approximation [2, 3, 8]. In principle, a complete decay of $\langle V_{{\bf q}s}(0)\cdot V_{{\bf q}s}(t)\rangle$ is required to obtain a reliable $G_{{\bf q}s}$. However, very long MD runs are needed to satisfy this condition for phonons with long lifetimes. This is inconvenient for ab initio MD simulations. Here, we extract $\widetilde{\omega}_{{\bf q}s}$ and $\tau_{{\bf q}s}$ from $\langle V_{{\bf q}s}(0)\cdot V_{{\bf q}s}(t)\rangle$ phenomenologically from relatively shorter MD runs. For a well-defined phonon quasiparticle of mode (${\bf q},s$), one can simply fit $\langle V_{{\bf q}s}(0)\cdot V_{{\bf q}s}(t)\rangle$ to the expression,

[TABLE]

where, $A_{{\bf q}s}$ is the vibrational amplitude. In practice, to obtain reliable $\widetilde{\omega}_{{\bf q}s}$ and $\tau_{{\bf q}s}$, the fitting needs only the numerical data of $\langle V_{{\bf q}s}(0)\cdot V_{{\bf q}s}(t)\rangle$ for the first few oscillation periods.

With $\widetilde{\omega}_{{\bf q}s}$ of those modes sampled by the MD run, an effective harmonic dynamical matrix can be constructed,

[TABLE]

where,

[TABLE]

and $\hat{{\bf e}}_{\bf q}$ is the harmonic eigenvector matrix:

[TABLE]

The effective harmonic force constant matrix $\widetilde{\Phi}({\bf r})$ can be obtained from the Fourier transformation of $\widetilde{D}({\bf q})$, from which,

[TABLE]

the effective harmonic dynamical matrix at arbitrary wave vector ${\bf q^{\prime}}$ is obtained. This way, $\widetilde{\omega}_{{\bf q^{\prime}}s}$ at any ${\bf q^{\prime}}$ in the Brillouin zone are calculated by diagonalizing $\widetilde{D}({\bf q^{\prime}})$.

This approach gives anharmonic phonon dispersions that are useful to study thermodynamic and lattice thermal transport properties of materials, especially at extreme conditions. It also offers a means of checking the existence of feasible vibrational modes even in the presence of strong anharmonicity. For example, we can witness the emergence of the cubic phase of CaSiO3 perovskite at high temperature and high pressure [9]. The concept of phonon quasiparticle also helps to easily resolve mode frequencies in the vibrational spectrum that overlap in energy even in the presence of crystal structure complexity. This advantage has been demonstrated on MgSiO3 perovskite, a typical non-trivial crystal structure with Pbnm symmetry and 20 atoms/cell. The obtained subtle and irregular frequency shifts with temperature are in good agreement with experiments [10]. Recently, we applied this approach to beryllium metal to address phase transitions at temperatures near melting. Our calculations confirmed that the long-debated and controversial hcp$\rightarrow$bcc transition occurs at very high temperatures, just before melting [11]. While the hcp structure was shown to be dynamically stable at all temperatures up to the transition temperature, the bcc phase was dynamically stabilized by anharmonic fluctuations just below the transition temperatures. At the phase boundary, both structures had phonon quasiparticles well defined and that was sufficient to enable free energy and phase boundary computations at those high temperatures. In addition, phonon lifetimes obtained this way allow us to study lattice thermal conductivity using the Boltzmann transport equation (BTE) approach. For instance, we demonstrated for the first time the breakdown of the well-known phonon minimal mean free theory in MgSiO3 perovskite [12], a result later confirmed by lowest-order perturbation theory [13].

The procedure to compute quasiparticle properties consists of four key steps. Step 1, atomic velocities are calculated by carrying out MD simulations and recording velocities [Eq. (3)]. Step 2, projecting instantaneous atomic velocities into the harmonic normal modes (HNM) [Eq. (2)]. Step 3, computing mode-projected velocity auto-correlation functions (VAF) and extracting quasiparticle frequencies and lifetimes from VAFs [Eq. (1)(4)(5)]. Step 4, constructing effective harmonic dynamical matrices using the renormalized phonon frequencies and HNMs and using them to obtain anharmonic phonon dispersions [Eq. (6)(7)(8)(9)].

In Step 3, three numerical techniques are provided to obtain quasiparticle properties: (i) Fitting the VAFs to an analytical function, Eq. (5), (ii) Fourier transforming (FT) Appendix A [14] VAFs, and (iii) using the maximum entropy method (MEM) Appendix B [14]. Also, we note that Step 4, critical to obtaining anharmonic phonon dispersions and free energies, is missing from a similar work [15].

2 Description of the program

The phq package folder is contained in the compressed phq.zip file. Inside the folder, there are LICENSE file, README.md file and three sub-folders. There is an example of diamond Si in the example sub-folder. Four makefiles are in the system sub-folder. makefile_linux is for Linux or macOS operating systems and makefile_windows is for Windows operating system. Makefiles for both gfortran and ifort compilers are also prepared. ifort compiler is recommended for faster running speed and smaller memory requirement. One can modify properly the FC flag to enable other compilers.

The source code is in the src sub-folder with three Fortran90 files: configure.f90 contains three auxiliary modules, which configure the physical quantities, parameters and text read in settings. main.f90 is the main module where: (i) the mode-projected velocities are calculated from the MD trajectory; (ii) renormalized phonon frequencies and lifetimes are extracted from the mode-projected VAFs; (iii) effective harmonic force constants and thus anharmonic phonon dispersions are obtained. phq.f90 runs the program and records the running time. A detailed description of subroutines used in the main.f90 file is listed in Table 1.

After making the file, users need to put the phq executable into the same folder of those four input files in order to run the code. Run command ./phq < input to perform the calculation. The input files will be introduced in the Input files section. The code workflow is shown in Fig. 1.

3 Input files

There are 4 input files for the phq program: input, scf.out, dyn.out and md.out, which provide the general settings, structure, harmonic phonon, and MD information separately. Users can prepare the input files according to the example provided within the code package. The instructions for preparing the input files are as following:

3.1 input

input provides general settings for the calculation. The values of the parameters should be set by the user. The details of the parameters are described in Table 2.

3.2 scf.out

scf.out provides structure information of the primitive cell from self-consistent calculations. The values of the parameters should be same as those used in first-principles calculations. The details of the parameters are described in Table 3.

3.3 dyn.out

dyn.out provides harmonic phonon information. The format follows the ph.x executable’s output of the Quantum ESPRESSO [16] suite. Quantum ESPRESSO users can simply obtain the dyn.out file by: (i) Turn off the symmetry in the self-consistent calculations by pw.x. (ii) Do the harmonic phonon calculations by ph.x. (iii) Concatenate dynmat1, …, dynmatnq together, where nq is the total number of the q-point. The details of the parameters are described in Table 4.

If other harmonic phonon calculation program is used, the dyn.out file should be prepared according to Table 4. It is further illustrated in Fig. 2, where the harmonic phonon information of the first q-point ($\Gamma$ point) of diamond silicon is shown as an example. The yellow and blue parts are required and optional respectively as the program input while the green part is not required as input. Parameters that are required or optional need to be prepared by users, while parameters that are not required can be left blank. Note that the character strings Dielectric, Effective, Diagonalizing, q and freq are arguments that will be read in by the phq program and therefore should not be modified. Dielectric Tensor and Effective Charges are only needed to be provided in the harmonic phonon information of the first q-point.

3.4 md.out

md.out provides MD information. The initial atomic coordinates of the supercell before the MD run and the instantaneous atomic coordinates during the MD simulation need to be recorded. The details of the parameters are described in Table 5. Note that for either atomic_positions or atomic_md_positions, atoms in each one primitive cell should be together instead of atoms of the same element being together. It is further illustrated in Fig. 3, which shows an example md.out file of $3\times 3\times 3$ supercell of Pbnm MgSiO3 in two formats. Fig. 3(a) takes the convention of Quantum ESPRESSO where atoms in one primitive cell are put together, while in Fig. 3(b), all Mg atoms, Si atoms and O atoms are put together. Here, format demonstrated in Fig. 3(a) should be taken for md.out. Also, the order of atoms in each of the primitive cell should be in the same order as provided in scf.out. The way to extend the primitive to generate the supercell for the MD simulation can be determined by the users.

4 Ouput files

A detailed description of all the output files of phq program are summarized in Table 6. The main output files are frequency.freq, which records the renormalized frequencies, and dynmatmd1, …, dynmatmdnq, which record the effective harmonic phonon information with the same format of ph.x executable’s output and can be read in by the q2r.x executable in the Quantum ESPRESSO suite. Anharmonic force constants, renormalized phonon dispersions, anharmonic entropy and free energy can be obtained from q2r.x executable’s output. If other Fourier interpolation program is used, either use the effective harmonic dynamical matrices file dynamical_matrix_md.mat or use the effective harmonic force constant matrix file gamma_matrix.mat to fit the format.

5 Example

Here we illustrate the performance of the phq code using diamond Si as an example. Anharmonic phonon dispersion, vibrational density of states, and anharmonic free energy are obtained. We use the density-functional theorem based Vienna ab initio simulation package (VASP) [17, 18] to carry out MD simulations on $4\times 4\times 4$ supercells (128 atoms) of Si. We employ the generalized gradient approximation of Perdew, Burke, and Ernzerhof [19] and the projector-augmented wave method [20] with an associated plane-wave basis set energy cutoff of 320 eV. Simulations are carried out at the volume with a zero static pressure. Temperatures ranging from 300 K to 1500 K are controlled through the Nosé dynamics [21]. The MD runs for 50 ps with a time step of 1 fs. Harmonic phonons are calculated using density-functional perturbation theory (DFPT) [22].

5.1 Obtaining VAFs

The mode-projected VAFs, $\langle V_{{\bf q}s}(0)V_{{\bf q}s}(t)\rangle$, are given by corr.vaf file. Fig. 4(a) displays the calculated $\langle V_{{\bf q}s}(0)V_{{\bf q}s}(t)\rangle$ at several temperatures of the LA/LO mode at ${\bf q}=W$ ([1/4, 3/4, 1/2]). At each temperature, the VAF amplitude decays exponentially. Accordingly, its power spectrum obtained via a Fourier transformation (FT) of VAFs has a Lorentzian line shape with a main single peak, Fig. 4(c). This behavior demonstrates the validity of the concept of phonon quasiparticle and allows us to identify the renormalized phonon frequency ($\widetilde{\omega}_{{\bf q}s}$) and line width ($\widetilde{\Gamma}_{{\bf q}s}$) of the phonon mode. Note that the file recording the FT is corr_fourier.vaf.

In practice, however, it is not a optimal way to find $\widetilde{\omega}_{{\bf q}s}$ and $\widetilde{\Gamma}_{{\bf q}s}$ directly from the power spectrum obtained from FT. Instead, we fit the VAFs to an exponentially decaying cosine function [2, 3, 11, 8] to extract $\widetilde{\omega}_{{\bf q}s}$ and $\widetilde{\Gamma}_{{\bf q}s}$, Eq. (5). From our previous work [10, 12, 11], $\widetilde{\omega}_{{\bf q}s}$ and $\widetilde{\Gamma}_{{\bf q}s}$ can be obtained by performing the fitting only for the first few oscillation periods of VAFs. Therefore, to do the fitting, one only needs to calculate VAFs for a relatively short correlation time. Note that the file recording the fitting is corr_fit.vaf.

The fitted curve at each temperature is shown in Fig. 4(b), matching the VAF very well. Our code also offers a way to obtain $\widetilde{\omega}_{{\bf q}s}$ by maximum entropy method (MEM) [15]. A list of $\widetilde{\omega}_{{\bf q}s}$ obtained by different methods can be found in the frequency.freq file. Users can choose renormalized frequencies obtained by fitting approach, FT or MEM to obtain the effective force constants and dynamical matrices by specifying 0, 1 or 2 for the method in the input file.

5.2 Obtaining frequency shifts

As a demonstration of the validity of the present approach, we compare the calculated thermal frequency shifts given by frequency.freq file with experimental measurements. Fig. 5(a)-5(c) display frequency shifts $\Delta\omega_{{\bf q}s}=\widetilde{\omega}_{{\bf q}s}|_{V=V_{eq}}-\omega_{{\bf q}s}$ of a TO/LO mode at the $\Gamma$ point, TA and TO modes at the $X$ point ([0, 1/2, 1/2]), and at the $L$ point ([1/2, 1/2, 1/2]) at the equilibrium volume of zero static pressure. Because experiments were conducted under ambient pressure, the renormalized phonon frequencies obtained under constant volume conditions have to be corrected. We indicate that such correction can be properly captured by considering the thermal shifts caused by the temperature induced volume expansion [10],

[TABLE]

where $\gamma_{{\bf q}s}$ is the mode Grüneisen parameter and $\bar{\alpha}$ is the temperature averaged thermal expansion coefficient:

[TABLE]

The mode Grüneisen parameter connects the phonon frequency variation with the volume expansion and can be obtained by DFPT [22]. It is defined as [23]:

[TABLE]

In practice the weak temperature dependence can be ignored [24]:

[TABLE]

where $\omega_{{\bf q}s}(0)$ is the harmonic phonon frequency at volume $V(0)$. $\alpha(T)$ is the thermal expansion coefficient which can be obtained by QHA. It is defined as:

[TABLE]

With this, according to Eq. (11), the temperature averaged thermal expansion coefficient can be onbtained:

[TABLE]

Overall, our calculated $\Delta\omega_{{\bf q}s}|_{P=0}=\widetilde{\omega}_{{\bf q}s}|_{P=0}-\omega_{{\bf q}s}$ are in good agreement with experiments [25, 26], Fig. 5(d)-5(f). We note that at the $X$ point, TA and TO modes have similar frequency shifts at constant volume. After the incorporation of the volume expansion, i.e., at constant pressure, their frequency shifts differ significantly. This is essentially caused by the different Grüneisen parameters of these two modes.

5.3 Obtaining temperature dependent phonon dispersion

Temperature dependent renormalized phonon frequencies obtained as such are only for those few q-points sampled by MD simulations, which are not sufficient to converge the calculation of thermal properties. Nevertheless, these frequencies provide a basis to obtaining phonon frequencies at q-points throughout the Brillouin zone. Here, we rely on an effective harmonic dynamic matrix, Eq. (6), and the derived effective harmonic force constants, Eq. (7), which give us the effective harmonic dynamic matrix at arbitrary q-points via a Fourier interpolation, Eq. (8) [10].

To construct the effective harmonic dynamic matrix along this line, users first find files of dynmatmd0, dynmatmd1, …, dynmatmd64 given by phq, and copy them to the phq/example/Si/postprocessing/dispersion sub-folder. Files q2r.in, dispersion.in and plotband.in files are also there. Next run:

q2r.x < q2r.in > q2r.out

matdyn.x < dispersion.in > dispersion.out

plotband.x < plotband.in > plotband.out

to calculate anharmonic phonon dispersion of diamond silicon, where matdyn.x and plotband.x are also Quantum ESPRESSO executables.

Fig. 6 (a) displays the obtained anharmonic phonon dispersions at 600 K and 1200 K. The harmonic phonon dispersions calculated by DFPT at static zero-pressure equilibrium volume is also shown for comparison. It can be seen that the temperature effect on phonon frequencies is discernible, and the overall trend is that phonon frequencies decrease as temperature increases [28]. Accordingly, the major peaks of the obtained vibrational density of states (VDoS) shift to lower frequencies with increasing temperature, Fig. 6 (b). Note that the VDoSs are calculated with a $20\times 20\times 20$ q-point grid, equivalent to a supercell of 16000 atoms in a direct MD simulation, which is well beyond the capability of ab initio methods.

5.4 Obtaining anharmonic vibrational entropy and free energy

With the renormalized phonon frequencies, the entropy is obtained by an integration over the whole Brillouin zone [6],

[TABLE]

where $n_{{\bf q}s}=[{\rm exp}(\hbar\widetilde{\omega}_{{\bf q}s}/k_{B}T)-1]^{-1}$. The free energy is obtained by numerical integration [10]:

[TABLE]

where $E_{0}$ is the static ground state energy and $\omega_{{\bf q}s}$ is the harmonic phonon frequency.

It can be seen that $\omega_{{\bf q}s}$ is key quantity to calculate $S$ and $F$. To obtain $\omega_{{\bf q}s}$ over the whole Brillouin zone, users copy the dynmatmd0, dynmatmd1, …, dynmatmd64 files obtained from phq to the phq/example/Si/postprocessing/vdos sub-folder. q2r.in and vdos.in files are also there. Then, run:

q2r.x < q2r.in > q2r.out

matdyn.x < vdos.in > vdos.out

to calculate anharmonic VDoS of diamond silicon on a $20\times 20\times 20$ q-point grid. $S$ and $F$ can be readily obtained from output VDoS file vdos according to Eq. (16)(17).

Fig. 7 compares $S$ and $F$ obtained with different q-grids: $4\times 4\times 4$ and $20\times 20\times 20$. The discernible differences $\Delta S$, Fig. 7(a), and thus $\Delta F$, Fig. 7(b) reveal that it is necessary to use a sufficiently large number of q-point in order to converge thermal properties.

6 Conclusion

In this communication we report phq, a code for extracting phonon quasiparticle properties, obtaining temperature dependent renormalized phonon dispersions, and calculating anharmonic free energies from harmonic phonon calculations and MD simulations. The program can be installed on Linux, macOS and Windows operating systems. Outputs from VASP [17, 18] and Quantum ESPRESSO [16] are both accepted as inputs here. Extension of compatibility with ABINIT [29], LAMMPS [30], Phon [31] and Phonon [32] for phonon dispersions processing are planned for the near future.

Acknowledgments

This work was supported primarily by NSF grant EAR-1503084. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number TG-DMR170058. DBZ acknowledges support from NSFC grant 11874088 and U1530401. TS is supported by NSFC 41474069. Computations are performed at the Stampede2 (the flagship supercomputer at the Texas Advanced Computing Center (TACC), University of Texas at Austin) and at Beijing Computational Science Research Center.

Appendix A Fourier transformation

The Fourier transformation (FT) of function $c(t)$ ($t\geq 0$),

[TABLE]

is done by mapping the $N$ consecutive evenly sampled points in the time domain to the frequency domain

[TABLE]

where $\omega_{j}$ is the window function. Then the discrete power spectrum is defined at $N$ nonnegative frequencies as

[TABLE]

$f_{k}$ is defined only for the zero and positive frequencies

[TABLE]

where $f_{c}$ is the upper limit of the frequency, known as Nyquist frequency. The actual frequency range considered in the calculation is limited by the window function. In phq program, by default the upper limit of the window function of each mode is set to be integral period where the amplitude of the velocity autocorrelation function (VAF) drop to 0.75 of its maximum. This value can be modified by users in subroutine correlation of the main.f90 file. By default, rectangular window function is used

[TABLE]

which is ready to obtain power spectrum of Lorentzian line shape. To avoid frequency leakage into neighboring bins, Hann window is also provided

[TABLE]

Users can enable Hann window in subroutine corr_fourier.

Appendix B Maximum entropy method

Maximum entropy method (MEM) is also known as all-poles model, autoregressive model, which is used to obtain smooth power spectrum in signal processing. Here we assume that the power spectrum of each mode follows the approximation

[TABLE]

where $a_{0}$ as well as $a_{j}$ is a set of coefficients of this approximation and $\Delta$ is the sampling interval in the time domain. $a_{0}$ and $a_{j}$ can be obtained from linear prediction method (LP) [14]

[TABLE]

where $y_{n}$ are a series of input and also predicted values by the LP, $x_{n}$ is the discrepancy between the data and the predicted value. Here we use the polarized velocities as the input of LP. A recursive procedure called Burg method [14] is adopted to minimize the discrepancy $x_{n}$, which is performed in subroutine linear_response of the main.f90 file. Regarding it to take input $x_{n}$ into output $y_{n}$, LP has a filter function

[TABLE]

Thus, the power spectrum of the $y_{n}$ equals the power spectrum of the $x_{n}$ multiplied by $\left|\mathcal{H}(f)\right|^{2}$. Therefore, the coefficients $a_{0}$ and $a_{j}$ can be related to the LP coefficients by

[TABLE]

Then the MEM power spectrum can be obtained from Eq. (24). Phonon quasiparticle properties are extracted from the power spectrum, which is done in subroutine maximum_entropy. A fitting of the MEM power spectrum of each mode to a Lorentzian line shape performed in subroutine lorentzian also give quasiparticle properties for comparison.