6.3. ANPHON: Theoretical background
Dynamical matrix
The dynamical matrix is given by
where \(M_{\kappa}\) is the atomic mass of atom \(\kappa\). By diagonalizing the dynamical matrix, one can obtain \(m\) (\(=3N_{\kappa}\)) eigenvalues \(\omega_{\boldsymbol{q}j}^{2}\) (\(j = 1, 2, \dots, m\)) and corresponding eigenvectors \(\boldsymbol{e}_{\boldsymbol{q}j}\) for each \(\boldsymbol{q}\) point. Here, \(\boldsymbol{e}_{\boldsymbol{q}j}\) is a column vector consisting of atomic polarization \(e_{\mu}(\kappa;\boldsymbol{q}j)\). Let \(D(\boldsymbol{q})\) denote a matrix form of equation (1), the eigenvalues may be written as
Next, we introduce \(m\times m\) matrices \(\Lambda\) and \(W\) which are defined as \(\Lambda(\boldsymbol{q}) = \mathrm{diag} (\omega_{\boldsymbol{q}1}^{2},\dots,\omega_{\boldsymbol{q}m}^{2})\) and \(W(\boldsymbol{q}) = (\boldsymbol{e}_{\boldsymbol{q}1},\dots,\boldsymbol{e}_{\boldsymbol{q}m})\), respectively. Then, equation (2) can be denoted as
Non-analytic correction
When one needs to capture the LO-TO splitting near the zone-center, it is necessary to account for the long-range correction to the dynamical matrix.
In the anphon code, three different approaches are implemented: the Parlinski’s method [1] (NONANALYTIC = 1), the mixed-space approach [2] (NONANALYTIC = 2), and the Ewald summation method (NONANALYTIC = 3).
The first two methods are rather simple and add the non-analytic part of the dynamical matrix with different weighting factors. The non-analytic part of the dynamical matrix is defined as
where \(V\) is the volume of the primitive cell, \(Z_{\kappa}^{*}\) is the Born effective charge tensor of atom \(\kappa\), and \(\epsilon^{\infty}\) is the electronic dielectric permittivity tensor, respectively.
Parlinski’s approach (NONANALYTIC = 1)
In the Parlinski’s approach [1] , the total dynamical matrix is given by
where \(\sigma\) is a damping factor (NA_SIGMA).
\(\sigma\) must be chosen carefully so that the non-analytic contribution becomes negligible at Brillouin zone (BZ) boundaries.
Too large \(\sigma\) can result in a discontinuity of phonon dispersion at the BZ boundaries, while too small \(\sigma\) can yield a large oscillation of the LO phonon frequencies as going away from \(\boldsymbol{q}=0\). So, an appropriate \(\sigma\) value needs to be chosen in an ad hoc way.
Mixed-space approach (NONANALYTIC = 2)
In the mixed-space approach [2] , the total dynamical matrix is given by
The second term vanishes at the \(\boldsymbol{q}\) points commensurate with the supercell except for the \(\Gamma\) point (\(\boldsymbol{q} = 0\)). No adjustable parameters exist in this approach. In many cases, the mixed-space approach give reasonable dispersion curves. However, it can still yield a discontinuity in the phonon dispersion at BZ boundaries when the boundary point is incommensurate with the supercell size.
Ewald method (NONANALYTIC = 3)
The Ewald summation technique [10] is more accurate and therefore recommended in general. In this method, we first decompose the harmonic force constants computed in a supercell into the long-range dipole-dipole (DD) components and the short-range ones as follows
After the short-range components are obtained, they are used in equation (1) to obtain the short-range dynamical matrix \(D_{\mu\nu}^{\mathrm{SR}}(\kappa\kappa^{\prime};\boldsymbol{q})\). Then, the total dynamical matrix is constructed as
The long-range terms, \(\Phi_{\mu\nu}^{\mathrm{DD}}(\ell\kappa;\ell^{\prime}\kappa^{\prime})\) and \(D^{\mathrm{DD}}_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q})\), can be calculated from \(\epsilon^{\infty}\) and \(Z_{\kappa}^{*}\) using the Ewald summation technique.
To include the non-analytic correction with NONANALYTIC > 0, one also need to give BORNINFO.
Group velocity
The group velocity of phonon mode \(\boldsymbol{q}j\) is given by
To evaluate the group velocity numerically, we employ a central difference where \(\boldsymbol{v}\) may approximately be given by
If one needs to save the group velocities, please turn on the PRINTVEL-tag.
Mode effective charge
The mode effective charge is defined as
where \(Z^{*}_{\kappa,\alpha\beta}\) is the atomic Born effective charge. To compute the mode effective charges of the zone-center modes, please set ZMODE = 1 in the &analysis field and supply BORNINFO.
Thermodynamics functions
The specific heat at constant volume \(C_{\mathrm{v}}\), the internal energy \(U\), the vibrational entropy \(S\), and the Helmholtz free energy \(F\) of individual harmonic oscillator are given as follows:
Here, \(k\) is the Boltzmann constant. These quantities are saved in the PREFIX.thermo file.
When the self-consistent phonon mode (MODE = SCPH) is selected, the anharmonic free-energy
defined by the following equation will be calculated and saved in the PREFIX.scph_thermo file:
Details of the derivation of the above expression can be found in Ref. [7].
Mean square displacement
The displacement-displacement correlation function is given by
where \(n_{\boldsymbol{q}j} = 1/(e^{\hbar\omega_{\boldsymbol{q}j}/kT}-1)\) is the Bose-Einstein distribution function. When UCORR = 1, the code prints the above correlation function in PREFIX.ucorr. The vector \(\boldsymbol{r}(\ell')\) is the lattice translation vector to the \(\ell'\) th unit cell, which can be specified by the SHIFT_UCORR tag.
When PRINTMSD is turned on, the code print the mean square displacements, which are the diagonal components of the correlation function:
Phonon DOS
When KPMODE = 2, the program anphon saves the (one) phonon density of states (DOS) to the file PREFIX.dos.
The one-phonon DOS is given by
If PDOS = 1 is given, the program also prints the atom-projected phonon DOS which is given by
In addition, TDOS-tag is available to compute the two-phonon DOS defined by
where \(\boldsymbol{G}\) is a reciprocal lattice vector. The sign \(\pm\) correspond to absorption and emission processes, respectively. Please note that the computation of the two-phonon DOS can be expensive especially when \(N_{q}\) or \(N_{\kappa}\) is large.
(Atomic) participation ratio
Participation ratio (PR) and atomic participation ratio (APR) defined in the following may be useful to analyze the localized nature of the phonon mode \(\boldsymbol{q}j\).
Participation ratio (PR)
Atomic participation ratio (APR)
For an extended eigenmode, the PR value is of order 1, whereas for a localized eigenmodes PR is of order \(1/N_{\kappa}\) [3]. APR is an atomic decomposition of PR that satisfies \(PR_{\boldsymbol{q}j} = \sum_{\kappa} (APR_{\boldsymbol{q}j,\kappa})^{2}\). To print the PR and APR, please set MODE = phonons and PRINTPR = 1 in the &analysis entry field.
Scattering phase space
When KPMODE = 2 and SPS = 1, the three-phonon scattering phase space \(P_{3}\) is calculated and saved to the file PREFIX.sps. \(P_{3}\) is defined as
where \(m\) is the number of phonon branches and
anphon also print the total scattering phase space
When SPS = 2, the three-phonon scattering phase space with the occupation factor \(W_{3}^{(\pm)}\) will be calculated and saved to the file PREFIX.sps_Bose. \(W_{3}^{(\pm)}\) is defined as
Here, \(n_{1}=n(\omega_{\boldsymbol{q}_{1}j_{1}})\) and \(n_{2}=n(\omega_{\boldsymbol{q}_{2}j_{2}})\) where \(n(\omega) = \frac{1}{e^{\hbar\omega/k_B T}-1}\) is the Bose-Einstein distribution function. Since \(n(\omega)\) is temperature dependent, \(W_{3}^{(\pm)}\) is also temperature dependent. The file PREFIX.sps_Bose contains \(W_{3}^{(\pm)}\) for all phonon modes at various temperatures specified with TMIN, TMAX, and DT tags.
Grüneisen parameter
The mode Grüneisen parameter, defined as \(\gamma_{\boldsymbol{q}j} = - \frac{\partial \log{\omega_{\boldsymbol{q}j}}}{\partial \log{V}}\), is calculated by
where \(\delta D(\boldsymbol{q})\) is a change in the dynamical matrix due to a volume change \(\delta V\), which is given by
Please set GRUNEISEN = 1 and give an appropriate FCSXML file containing cubic IFCs to print Grüneisen parameters.
Anharmonic self-energy
The anharmonic self-energy due to cubic anharmonicity to the lowest order is given by
where \(i\omega_{m}\) is the Matsubara frequency. In equation (3), we simply denoted \(\omega_{\boldsymbol{q}_{i}j_{i}}\) as \(\omega_{i}\). The matrix element \(V^{(3)}\) is given by
which becomes zero unless \(\boldsymbol{q}+\boldsymbol{q}^{\prime}+\boldsymbol{q}^{\prime\prime}\) is an integral multiple of \(\boldsymbol{G}=n_{1}\boldsymbol{b}_{1}+n_{2}\boldsymbol{b}_{2}+n_{3}\boldsymbol{b}_{3}\). Phonon linewidth \(\Gamma_{\boldsymbol{q}j}\), which is the imaginary part of the phonon self-energy, can be obtained by the analytic continuation to the real axis (\(i\omega_{m}\to \omega + i0^{+}\)) as
The computation of equation (4) is the most expensive part of the thermal conductivity calculations.
Therefore, we employ the crystal symmetry to reduce the number of triplet pairs \((\boldsymbol{q}j,\boldsymbol{q}^{\prime}j^{\prime},\boldsymbol{q}^{\prime\prime}j^{\prime\prime})\) of \(V^{(3)}\) to calculate.
To disable the reduction, please set TRISYM = 0.
Isotope scattering
The effect of isotope scatterings can be considered by the mass perturbation approach proposed by S. Tamura [4] by the ISOTOPE-tag.
The corresponding phonon linewidth is given by
where \(g_{2}\) is a dimensionless factor given by
Here, \(f_{i}\) is the fraction of the \(i\)th isotope of an element having mass \(m_i\),
and \(M_{\kappa}=\sum_{i}f_{i}m_{i}(\kappa)\) is the average mass, respectively.
The \(g_{2}\) values should be provided by the ISOFACT-tag.
The average mass \(M_{\kappa}\) is substituted by the value specified in the MASS-tag.
Lattice thermal conductivity (Peierls term)
The lattice thermal conductivity tensor \(\kappa_{\mathrm{ph}}^{\mu\nu}(T)\) is estimated within the relaxation-time approximation as
where \(V\) is the unit cell volume, \(c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T\), and \(\tau_{\boldsymbol{q}j}(T)\) is the phonon lifetime. The phonon lifetime is estimated using the Matthiessen’s rule as
The lattice thermal conductivity is saved in the file PREFIX.kl.
The spectra of the lattice thermal conductivity \(\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)\) can also be calculated by setting KAPPA_SPEC = 1 in the &analysis field. \(\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)\) is defined as
If we integrate this quantity over \(\omega\), we then obtain the bulk thermal conductivity, namely \(\kappa_{\mathrm{ph}}^{\mu\mu} = \int_{0}^{\infty} \kappa_{\mathrm{ph}}^{\mu\mu}(\omega) \; \mathrm{d}\omega\).
Cumulative thermal conductivity
The accumulative lattice thermal conductivity \(\kappa_{\mathrm{ph,acc}}^{\mu\nu}(L)\) is defined as
where \(\Theta(x)\) is the step function. This quantity can be calculated by using the script analyze_phonons.py with --calc cumulative flag.
One can also use another definition for the accumulative thermal conductivity:
In this case, the contribution to the total thermal conductivity is limited only from phonon modes whose mean-free-path along the \(\mu\)-direction is smaller than \(L\).
To calculate this, please use the --calc cumulative2 flag and specify the direction \(\mu\) by the --direction option.
Coherent component of lattice thermal conductivity
The coherent components of lattice thermal conductivity (see Ref. [8]), which are associated with the band off-diagonal components of the harmonic heat-flux operator, is calculated as
where \(c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T\) and \(\Gamma_{\boldsymbol{q}j}\) is the total phonon linewidth (half width) of phonon mode \(\boldsymbol{q}j\).
\(\boldsymbol{v}_{\boldsymbol{q}jj'}\) is a band off-diagonal generalization of the group velocity [9]. When KAPPA_COHERENT = 1 | 2 the coherent component is calculated and saved in PREFIX.kl_coherent. When KAPPA_COHERENT = 2, all components of the coherent term before summation are saved in PREFIX.kc_elem.
Delta function
To compute the phonon DOSs and the imaginary part of phonon self-energies,
it is necessary to evaluate the Brillouin-zone integration containing Dirac’s delta function.
For that purpose, we provide 3 options through the ISMEAR-tag.
When ISMEAR = 0, the delta function is replaced by the Lorentzian function as
When ISMEAR = 1, the delta function is replaced by the Gaussian function as
which decays faster than the Lorentzian function.
For both cases, \(\epsilon\) should be given by the EPSILON-tag, which must be chosen carefully
to avoid any unscientific results. \(\epsilon\) should be small enough to capture detailed phonon structures
such as phonon DOS or energy conservation surface related to three-phonon processes, but it should be large
enough to avoid unscientific oscillations. Choosing an appropriate value for \(\epsilon\) is not a trivial task
since it may depend on the phonon structure and the density of \(\boldsymbol{q}\) points.
To avoid such issues, the program anphon employs the tetrahedron method [5] by default (ISMEAR = -1)
for numerical evaluations of Brillouin zone integration containing \(\delta(\omega)\).
When the tetrahedron method is used, the EPSILON-tag is neglected.
We recommend using the tetrahedron method whenever possible.
Self-consistent phonon (SCPH) calculation
The self-consistent phonon mode (MODE = SCPH) computes temperature-dependent phonon frequencies by solving the following equation self-consistently [6]:
Here, \(\omega_{\boldsymbol{q}j}\) is the harmonic phonon frequency and \(F_{\boldsymbol{q}\boldsymbol{q}_{1},ijk\ell} = \Phi(-\boldsymbol{q}i;\boldsymbol{q}j;\boldsymbol{q}_{1}k;-\boldsymbol{q}_{1}\ell)\) is the reciprocal representation of fourth-order force constants computed using the harmonic eigenvectors. The updated phonon frequency in the \(n\)th iteration is obtained by diagonalizing the matrix \(V_{\boldsymbol{q}ij}^{[n]}\) as
where \(\omega_{\boldsymbol{q}j}^{[n]} = (\Lambda^{[n]}_{\boldsymbol{q}jj})^{\frac{1}{2}}\) and \(C^{[n]}_{\boldsymbol{q}}\) is the unitary matrix that transforms the harmonic phonon eigenvectors into anharmonic ones as \(W^{[n]}_{\boldsymbol{q}} = W_{\boldsymbol{q}}C^{[n]}_{\boldsymbol{q}}\). The matrix \(\mathcal{K}\) in Eq. (5) is defined as
\(\alpha\) is the mixing parameter, which can be changed via the MIXALPHA tag, and \(n(\omega)\) is the Bose-Einstein distribution function.
When CLASSICAL = 1 is given in input, the expectation value of the mean square displacement is computed using classical statistics, where the \(Q_{\boldsymbol{q}}\) matrix is replaced by
The SCPH equation is solved on the irreducible \(\boldsymbol{q}\) grid defined by the KMESH_INTERPOLATE tag.
The \(\boldsymbol{q}_{1}\) grid in Eq. (5), given by the KMESH_SCPH tag,
can be finer than the \(\boldsymbol{q}\) grid. After the SCPH iteration converges, the code computes the anharmonic correction to the harmonic force constant \(\Delta D(\boldsymbol{r}(\ell))\) as follows:
\(\Delta D(\boldsymbol{r}(\ell))\) is saved in PREFIX.scph_dfc2.
The most computationally expensive part is the calculation of matrix elements of \(F_{\boldsymbol{q}\boldsymbol{q}_{1},ijk\ell}\).
When SELF_OFFDIAG = 0 (default), the code only computes the elements of \(F_{\boldsymbol{q}\boldsymbol{q}_{1},iikk}\).
Therefore, the computational complexity is \(\mathcal{O}(N_{q}^{\mathrm{irred.}}N_{q_{1}}m^{2})\).
When SELF_OFFDIAG = 1, the off-diagonal elements are also calculated, and the computational complexity is \(\mathcal{O}(N_{q}^{\mathrm{irred.}}N_{q_{1}}m^{4})\).