FourPhonon: An extension module to ShengBTE for computing four-phonon scattering rates and thermal conductivity

• FourPhonon calculates lattice thermal conductivity and both three-phonon and four-phonon scattering rates.Four-phonon scattering rates are computed at the RTA level using an adaptive broadening scheme.Source: PROGRAM SUMMARY · p. 2
• The package is implemented in Fortran 90 with MPI and distributed under the GNU General Public License, version 3.Source: PROGRAM SUMMARY · p. 2
• Productive calculations require HPC facilities and may take several hours to several days to complete.Source: PROGRAM SUMMARY · p. 2

• Phonon-phonon interactions govern phonon properties and lattice thermal conductivity, which can be predicted from first principles using the phonon BTE.Source: §1 · p. 2
• Four-phonon scattering can be significant beyond the conventional three-phonon treatment, including in BAs and diverse other systems.Experiments confirmed the significance of four-phonon scattering in BAs, and later studies extended its relevance across insulators, semiconductors, and semimetals.Source: §1 · pp. 2–3
• FourPhonon extends ShengBTE to compute four-phonon phase space, normal and Umklapp scattering rates, and lattice thermal conductivity.The package includes Fourthorder.py for calculating fourth-order interatomic force constants through first-principles calculations.Source: §1 · p. 3
• The paper presents theoretical background, computational details, and case studies on Si, BAs, and LiCoO2.Source: §1 · p. 3

Workflow integration

• FourPhonon adds fourth-order force constants and a CONTROL-file flag to ShengBTE, enabling four-phonon utilities within the existing workflow.An additional input file supplies fourth-order force constants, and the four_phonon flag activates the added utilities.Source: §2.1 · pp. 3–4

BTE solution

• FourPhonon treats four-phonon scattering rates at the RTA level rather than including four-phonon phase space in the iterative BTE scheme.The full formalism contains iterative terms, but the current package excludes those terms because of exceptionally high memory demands.Source: §2.1 · p. 5, §2.2 · p. 6

Scattering formalism

• The implemented formalism computes three- and four-phonon scattering rates using normal and Umklapp processes, with energy conservation enforced through a delta function.Four-phonon processes use scattering probability matrices and distinguish reciprocal-lattice-vector conditions for normal and Umklapp events.Source: §2.2 · pp. 6–7

Fourth-order force constants

• Fourthorder.py calculates fourth-order IFCs with a real-space finite-difference method based on displaced-atom forces.The displacement h is set by default to 0.04 Bohr radius and can be changed in the program header.Source: §2.3 · p. 8
• Symmetry analysis reduces the otherwise impractical fourth-order IFC workload by applying index-transposition, space-group, and translational constraints.A fourth-order IFC calculation may require tens of thousands of DFT simulations for supercells exceeding 100 atoms, motivating symmetry reduction.Source: §2.3 · p. 9

Convergence and cost

• Fourth-order IFC convergence with respect to cutoff distance must be checked manually, although second-nearest-neighbor cutoffs often suffice for many three-dimensional crystals.Increasing the cutoff substantially raises computational cost; convergence can be assessed from force constants or thermal conductivity.Source: §2.3 · p. 10

Si

• 17 × 17 × 17 q-points suffice for converged Si thermal conductivity with four-phonon scattering, compared with 25 × 25 × 25 for three-phonon scattering.This reduced mesh requirement is beneficial because four-phonon calculations are computationally expensive.Source: §3.1 · p. 11
• In Si, Umklapp processes dominate four-phonon scattering, while normal three-phonon processes do not overwhelm Umklapp processes.The authors therefore consider an RTA treatment acceptable for Si thermal-conductivity calculations.Source: §3.1 · p. 11

BAs

• Four-phonon scattering rates in BAs can be comparable to or larger than three-phonon rates even at room temperature.A single-temperature scattering-rate calculation requires about 1680 CPU hours.Source: §3.2 · p. 13
• BAs thermal conductivity calculations give 1221.9 W/mK with scalebroad=1 and 1266.9 W/mK with scalebroad=0.1, a relative difference within 4%.The temperature-dependent results show larger reductions when four-phonon scattering is included at higher temperatures and agree well with experimental data.Source: §3.2 · pp. 13, 16

LiCoO2

• LiCoO2 has four-phonon scattering rates generally comparable to three-phonon rates across the frequency range at room temperature.Unlike BAs, four-phonon rates are less significant in optical branches; redistribution contributes most, with Umklapp dominating that channel.Source: §3.3 · p. 16
• At 300 K, LiCoO2 thermal conductivity is 9.35 W/mK in-plane and 1.39 W/mK through-plane, compared with reference values of 9.7 and 1.4 W/mK.The calculations use a 3 × 3 × 3 supercell, a 10 × 10 × 10 q-point grid, and scalebroad=0.01.Source: §3.3 · p. 16

• FourPhonon provides an open-source tool for four-phonon scattering calculations while preserving ShengBTE features.Source: §4 · p. 17
• The package extends adaptive broadening to four-phonon scattering and supplies Fourthorder.py for fourth-order force constants using point-group and translational symmetries.Source: §4 · p. 17
• Examples on Si, BAs, and LiCoO2 cover common materials, systems with significant four-phonon effects, and newly reported complex systems.Source: §4 · pp. 17–18
• The examples demonstrate convergence with broadening factors and q-point grids, plus contributions from process types and scattering channels.Source: §4 · p. 18
• Current calculations use the Single Mode Relaxation Time Approximation and can iterate with three-phonon scattering rates, but not fully four-phonon iterations.Source: §4 · p. 18

• FourPhonon requires the FORCE_CONSTANTS_4TH file, generated by Fourthorder.py, in addition to ShengBTE input files.Source: Appendix A · p. 19
• Users enable four-phonon calculations by setting four_phonon=.TRUE. in the CONTROL file.Source: Appendix A · p. 19
• BTE.Numprocess_4ph reports the number of allowed four-phonon scattering processes for each mode.Source: Appendix A · p. 19
• BTE.P4 and BTE.WP4 report available and weighted phase space for four-phonon processes.Source: Appendix A · p. 19
• BTE.w_4ph, BTE.w_4ph_Umklapp, and BTE.w_4ph_normal report total, Umklapp, and normal four-phonon RTA scattering rates.Source: Appendix A · p. 19

• For a mode (q, ω), three other phonons occupy q′, q′′, and q′′′ with corresponding frequencies, where q′′′ depends on q′ and q′′.Source: Appendix B · p. 19
• The second and third phonons are independent, while the fourth phonon is determined by their choice.Source: Appendix B · p. 19
• The recombination process uses ∆ω = −ω′ −ω′′ + ω′′′ and q + q′ + q′′ = q′′′.Source: Appendix B · p. 20
• The redistribution process uses ∆ω = −ω′ + ω′′ + ω′′′ and q + q′ −q′′ = q′′′.Source: Appendix B · p. 20
• The splitting process uses ∆ω = ω′ + ω′′ + ω′′′ and q −q′ −q′′ = q′′′, after which absolute values give a common formula.Source: Appendix B · p. 20