Computational protocol to evaluate electron-phonon interactions within density matrix perturbation theory

We present a computational protocol, based on density matrix perturbation theory, to obtain non-adiabatic, frequency-dependent electron-phonon self-energies for molecules and solids. Our approach enables the evaluation of electron-phonon interaction using hybrid functionals, for spin-polarized systems, and the computational overhead to include dynamical and non-adiabatic terms in the evaluation of electron-phonon self-energies is negligible. We discuss results for molecules, as well as pristine and defective solids.


INTRODUCTION
The study of electron−phonon interaction in solids can be traced back to the early days of quantum mechanics, 1 and it has been instrumental in explaining fundamental properties of solids, including conventional superconductivity. 2 However, it was not until recent years that electron−phonon interaction was computed from first-principles, 3−7 leading to nonphenomenological predictions of transport properties of solids 8 and of electron−phonon renormalizations of band structures. 7,9−13 Although early studies relied on semi-empirical models 14−16 to study electron−phonon interaction, modern investigations typically employ the frozen phonon (FPH) approach, 7,17 density functional perturbation theory (DFPT), 4,5,10,11,18 or molecular dynamics (MD) simulations, 12,19,20 with electron−electron and electron−ion interactions described at the level of density functional theory (DFT). 21 In two recent papers, 10,11 we combined first-principles calculations of electron−electron and electron−phonon selfenergies in molecules and solids, within the framework of DFPT. We developed an approach that enables the evaluation of electron−phonon coupling at the G 0 W 0 level of theory 22−25 for systems with hundreds of atoms and the inclusion of nonadiabatic and temperature effects at no additional computational cost. Recent developments have also extended the DFPT method to study electron−phonon interactions using DFT + U 26 to obtain an improved description of systems with strong electronic correlation. We also computed 12 electron− phonon renormalizations of energy gaps by using the pathintegral MD (PIMD) methods to investigate anharmonic effects in crystalline and amorphous solids. The DFPT-, FPH-, and MD-based methods are addressing different regimes and different problems; the use of DFPT and FPH is appropriate for systems whose atomic constituents all vibrate close to their equilibrium positions, although anharmonic effects have been included in some FPH calculations. 7 The assumption of close to equilibrium vibrations is however not required when applying PIMD, which thus has a wider applicability; for example, it can be used to study amorphous materials and molecules and solids exhibiting prominent anharmonic effects, for example, molecular crystals 27,28 and several perovskites. 29 However, the calculation of electron−phonon renormalizations using FPH-and MD-based methods is carried out within the Allen−Heine−Cardona (AHC) 30,31 formalism, which neglects dynamical and non-adiabatic terms of the electron−phonon self-energies. These effects have been shown to be essential to describe electron−phonon interactions in numerous polar materials, 32,33 for example, SiC. Perturbation-based methods, on the other hand, can accurately compute electron−phonon self-energies within and beyond the AHC formalism, thus including non-adiabatic and/or frequency-dependent effects into the self-energy. Another benefit of DFPT-based methods is the ability to explicitly evaluate the electron−phonon coupling matrices, which are useful quantities, for example, in the study of mobilities 34−40 and polaron formation. 41−46 Here, we generalize the perturbation-based approach of refs 10 and 11 to enable efficient calculations of electron−phonon interaction with hybrid functionals, by using density matrix perturbation theory (DMPT) 47,48 to compute phonons and electron−phonon coupling matrices. Our implementation takes advantage of the Lanczos algorithm, 49 which enables the calculations of electron−phonon self-energies beyond the AHC approximation, at no extra cost.
DMPT has been used in the literature to compute excitation energies and absorption spectra in molecules and solids in conjunction with time-dependent DFT (TDDFT) 50,51 and to solve the Bethe−Salpeter equation (BSE). 48,52−56 In the latter case, DMPT has been applied to obtain the variation of singleparticle wavefunctions due to the perturbation of an electric field. However, DMPT is a general formalism that can be used to compute the response of a system to perturbations of any form, including perturbations caused by atomic displacements.
In this paper, we first derive a formalism for phonon calculations within DMPT, starting from the quantum Liouville equation in Section 2; we then verify our results by comparing them with those of FPH and PIMD calculations in Section 3. We then present calculations of electron−phonon interactions in small molecules (Section 4) and pristine (Section 5) and defective diamonds (Section 6) using hybrid functionals, and we conclude the paper in Section 7 with a summary of our findings.

METHODOLOGY
Using Hartree atomic units (ℏ = e = m e = 1), we describe the electronic structure of a solid or molecule within Kohn−Sham (KS) DFT, and we consider the quantum Liouville's equation to describe perturbations acting on the system where [·, ·] denotes a commutator and H KS (t) is the KS Hamiltonian with K as the kinetic operator, V H as the Hartree potential, V ext as the external potential, and V xc as the exchange−correlation potential. The KS Hamiltonian does not depend explicitly on time and depends implicitly on time through the timedependent density matrix γ that can be written in terms of KS single-particle orbitals ψ nσ where σ is the spin polarization, n is the band index, and N occ σ is the number of occupied bands in the spin channel σ. Below we present calculations performed by sampling the Brillouin zone with only the Γ point and hence omit labeling of eigenstates with k-points.
Given a time-dependent perturbation ∂V ext (t) acting on the Hamiltonian, the first-order change in the density matrix ∂γ(t) satisfies the following equation where is the Liouville super-operator Here, we use the notation ∂ to represent a change in potentials (∂V), wavefunctions ∂ψ, charge densities ∂ρ(r), and density matrices ∂γ(r, r′); γ 0 and H 0 KS are the density matrix and the KS Hamiltonian of the unperturbed system, respectively.
Taking the Fourier transform of eq 4, we rewrite it in the frequency domain In phonon calculations, we adopt the Born-Oppenheimer approximation, 57 and no retardation effects are included. Hence, we only need to solve eq 6 at ω = 0 The equation above can be cast in the following form Ä  (8) where c is the projection operator onto the virtual bands; 1e , 2e , 1d , and 2d , defined below in eqs 12 and 13 and 15−18, are related to the variation of exchange−correlation In phonon calculations, the external perturbation is static ∂V ext (ω = 0), and eq 8 can be further simplified since a nσ = b nσ = ∂ψ nσ for static perturbations, and eq 8 becomes  (13) where f Hxc = v c + f xc is the sum of the bare Coulomb potential v c and the exchange−correlation kernel (14) with ρ(r) being the electron density. When using hybrid functionals, the operators are  (18) where f Hxc loc = v c + f xc loc is the sum of the bare Coulomb potential and the local part of the exchange−correlation kernel  (19) and the parameter α is the fraction of the Hartree−Fock exchange included in the definition of the hybrid functional. Note that the 1d and 2d operators are zero for LDA/GGA functionals. Once we have the solutions a nσ of the Liouville equation (eq 8 or eq 10), that is, the change in wavefunction ∂ψ nσ , we can compute the change in the density matrix with eq 9; the change in density is then given using the following expression  (20) and force constants are obtained as follows By diagonalizing the dynamical matrix where M I and M J are atomic masses, we obtain the frequency ω ν of mode ν and its polarization ξ Iα,ν .
To compute the electron−phonon coupling matrices in the Cartesian basis (23) or in the phonon-mode basis where ξ Iα,ν is the νth vibrational mode, we need to evaluate the change in the self-consistent (scf) potential ∂V scf . The scf potential is given by the sum of the Hartree potential V H , the local part of the exchange−correlation potential V xc loc , and the non-local Hartree-Fock exchange V xc nl . Thus, the change in the scf potential ∂V scf |ψ nσ ⟩ is the sum of the following three terms The Fan−Migdal and Debye−Waller self-energies can then be computed as occupation number of the KS eigenvalues ε mσ obeying the Fermi-Dirac distribution. The Debye−Waller self-energy is derived within the rigid-ion approximation (RIA), 30,59,60 which approximates second-order electron−phonon coupling matrices with first-order ones.
Using the frequency-dependent Fan−Migdal self-energy, the renormalized energy levels can be evaluated self-consistently with initial guess ω 0 = ε nσ and using the Lanczos 49 algorithm to evaluate the frequency-dependent Fan−Migdal self-energy (for a detailed derivation, see refs 10 and 11). Calculations using full frequency-dependent self-energies with DFPT have been recently reported in the literature. 10,32,43 We refer to the FM self-energy in eq 28 as the non-adiabatic fully frequency-dependent (NA-FF) self-energy. If the frequency dependence is considered within the adiabatic approximation, the self-energy is We refer to eq 31 as the adiabatic fully frequency-dependent (A-FF) self-energy.
In our formulation the evaluation of self-energies can be carried out simultaneously at multiple frequencies using the Lanczos algorithm; however, we introduce below approximations leading to frequency-independent self-energies for comparison with results present in the literature, obtained, for example, with the AHC formalism. 30,31 In particular, we evaluate the FM self-energy by applying the so-called on-themass-shell (OMS) approximation, that is, by setting ω = ε nσ in the expressions of the A-FF and NA-FF self-energies. In the former case, we obtain the adiabatic AHC (A-AHC) 30,31 approximation, and in the latter case, we obtain the nonadiabatic AHC (NA-AHC) approximation We summarize the various levels of approximations applied to evaluate the FM self-energy in Table 1; the corresponding DW self-energies are the same for all levels of approximation. Thus, we also use the acronyms A-AHC, NA-AHC, A-FF, and NA-FF to denote the level of theory adopted for the total selfenergy (FM + DW) and for the electron−phonon renormalization of fundamental gaps.

VERIFICATION
To verify the implementation of the method described above in the WEST 24 package, we first computed the phonon frequencies of selected solids (diamond, silicon, and silicon carbide) and the vibrational modes of selected molecules (H 2 , N 2 , H 2 O, and CO 2 ) and compared our results with those of the FPH approach. The displacements used for FPH calculations are 0.001 Å for all molecules and solids. In Tables  2 and 3, we summarize our results obtained at the PBE0 61 level of theory and obtained by solving either the Liouville equation or by using the FPH approach. The lattice constants used for diamond, silicon, and silicon carbide are 3.635, 5.464, and 4.372 Å, respectively, and the cell used for molecules is a cube of edge 10.583 Å. For verification purposes, we only computed the phonon modes at the Γ point in the Brillouin zone of the solids. We used an energy cutoff of 60 Ry for the solids and 50 Ry for the molecules and the SG15 62 ONCV 63 pseudopotentials for all solids and molecules. Table 2 shows that the absolute difference in the phonon frequencies computed with the FPH approach and the method implemented here is small for silicon and silicon carbide, 0.19 and 0.07 cm −1 , respectively. The corresponding difference for diamond is larger but still acceptable being below 5 cm −1 . In Table 3, we compare the vibrational frequencies of H 2 , N 2 , H 2 O, and CO 2 molecules computed by solving the Liouville equation and applying the frozen-phonon (FPH) approach. We found again that the differences are small for H 2 , N 2 , and CO 2 (below 1 cm −1 ), albeit slightly larger for H 2 O. The largest difference is found in the case of H 2 O (15.29 cm −1 ), and this is most likely due to the numerical inaccuracy of the FPH approach.
a In the main text, we use AHC and A-AHC interchangeably.  To verify our calculations of electron−phonon interactions, we carried out a detailed study of the renormalization of the HOMO−LUMO gaps (E g ) of the CO 2 , Si 2 H 6 , HCN, HF, and N 2 molecules, with the results for CO 2 being summarized in Table 4 and the rest in Table 5. Table 4 summarizes the renormalizations to the E g of the CO 2 molecule obtained within the A-AHC formalism and using DFPT, FPH, and PIMD 12 at the LDA, 64 PBE, 65 PBE0, 61 and B3LYP 66−68 levels of theory, respectively. With the LDA and PBE/GGA functionals, the solution of the Liouville equation yields the same results as the method proposed in refs 10 and 11, as expected. When solving the Liouville equation with the DFPT method, the RIA is adopted; however, the latter approximation is not used in the FPH approach, leading to a slight difference between the FPH and Liouville results. In addition, we carried out calculations with the hybrid functionals, PBE0 and B3LYP, and compared our results with those of the FPH and PIMD approaches. 12 The PIMD approach circumvents the RIA and also includes ionic anharmonic effects. Since the RIA is adopted and anharmonicity is not included in the Liouville approach, differences on the order of ∼30 meV, relative to PIMD, are considered as acceptable. We note that the computed renormalizations of the gap of CO 2 reported in the literature, 69 −680.7 and −716.2 meV with LDA 64 and PBE + TS 70 functionals, respectively, are significantly different from those obtained here. We also note that ref 69 reports the result at the B3LYP level of theory, − 4091.6 meV, which is one order of magnitude larger than the corresponding LDA and PBE + TS results, hence calling into question the numerical accuracy of the data. Such significant differences between our and the results of ref 69 probably stems from the different choices of basis functions, localized basis functions in ref 69 and plane waves in this work.
In addition to CO 2 , we also computed the energy gap renormalizations of Si 2 H 6 , HCN, HF, and N 2 molecules with the B3LYP functional; these are shown in Table 5 In summary, we have verified our implementation of phonon and electron−phonon interaction by comparing results computed with Liouville's equation and those obtained with DFPT, FPH, and PIMD methods. At the LDA/PBE level of theories, we obtain exactly the same results as with DFPT, as expected; at the hybrid functional level of theory, the results obtained with Liouville's equation are comparable with those of the FPH and PIMD methods, with reasonable differences compatible with the different approximations employed in the three different approaches.

ELECTRON−PHONON RENORMALIZATION OF ENERGY GAPS IN SMALL MOLECULES
Having verified our implementation, we carried out a study of the renormalization of the HOMO−LUMO gap of molecules in the G2/97 test set 71 with LDA, PBE, PBE0, and B3LYP functionals, 50 Ry plane wave energy cutoff, and SG15 62 ONCV 63 pseudopotentials. The results are summarized in Tables 6 and 7 and are illustrated in Figure 1. Table 6 summarizes the renormalizations computed with the A-AHC formalism. For most of the molecules, using hybrid functionals does not significantly change the gap renormalization relative to LDA or PBE results. For example, the energy gap renormalizations of the H 2 molecule computed with LDA, PBE, PBE0, and B3LYP functionals are 58, 61, 63, and 63 meV, respectively. However, hybrid functionals do reduce the magnitude of gap renormalization in several systems, and CO 2 and CH 3 Cl are representative examples. In CO 2 , the renormalization is reduced from the −415 meV (PBE) to −137 meV (PBE0) level of theory; in CH 3 Cl, it is reduced from −149 meV (PBE) to −59 meV (PBE0).
We report in Table 7 our results within the non-adiabatic AHC (NA-AHC) framework. The removal of the adiabatic approximation significantly influences the computed magnitude of the gap renormalization in most of the molecules, with some exceptions, for example, CO 2 . For example, the H 2 gap renormalization computed using PBE0 varies from 63 meV (AHC) to −377 meV (NA-AHC). The most significant differences are found for the F 2 and H 2 O 2 molecules, where the gap renormalizations computed at the PBE0 level of theory are 25 and −72 meV, respectively, within the AHC approach and −2914 and −891 meV, when using the non-adiabatic AHC method.  We emphasize that neither the AHC nor the non-adiabatic AHC formalism correctly describes the self-energies in the full energy range, and thus, we suggest that the frequencydependent self-energies should always be computed.

ELECTRON−PHONON RENORMALIZATION OF THE ENERGY GAP OF DIAMOND
We computed the electron−phonon renormalization of the energy gap in diamond within the AHC formalism and by computing the NA-FF self-energies self-consistently (see Table  1 and eq 30). The calculations for diamond were carried out in a 3 × 3 × 3 supercell with 60 Ry plane wave energy cutoff and SG15 62 ONCV 63 pseudopotentials.
In Figure 2, we present the temperature-dependent indirect gap renormalization computed with the PBE, PBE0, and dielectric-dependent hybrid (DDH) functionals, 72,73 where the fraction of exact exchange (0.18) in DDH is chosen to be the inverse of the dielectric constant of diamond (5.61). 72 Within the same level of approximation, for example, the AHC formalism (circles in the plot), the PBE, PBE0, and DDH results are almost the same for temperatures lower than 400 K, but their difference increases at higher temperatures. With the same functional, for example, the PBE0 functional (orange lines in the plot), the results obtained with the fully frequencydependent non-adiabatic self-energies are lower than those obtained with the AHC formalism. In general, the use of the hybrid functional does not significantly modify the trend of the ZPRs computed at the PBE level, as a function of temperature.
In Figure 2, we also report the renormalization of the indirect gap of diamond obtained with the FPH approach and the PBE0 functional. The results obtained with the FPH (purple line) approach and the Liouville equation (orange lines in the plot) are essentially the same below 300 K, but they differ as T is increased. The difference between the AHC/NA-FF and FPH approaches is always smaller than 10 meV at all temperatures, and it is reasonable considering that the FPH approach does not adopt the RIA, which is instead used within the AHC and NA-FF approaches.
A comparison of the computed and measured renormalized energy gap of diamond is given in Figure 3 and Table 8.    Table 1 for the definition of the approximations). We used different functionals specified in the inset.   Table 1. The results obtained with the FPH approach and the PBE0 functional are also reported for comparison. The renormalization energy at zero temperature has been shifted to zero.  Table 1.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article 100 K is 5.45 eV. 74 By including electron−phonon renormalization, we can see that the results computed at the PBE0 level of theory agree relatively well with the experimental measurements (see Figure 3 and Table 8). The renormalized indirect gap computed with the PBE0 functional at 100 K is 5.899 eV when the AHC formalism is used, and it is 5.735 eV when the NA-FF self-energies are used. The renormalized indirect gaps computed with the DDH functional at 100 K are 5.308 (AHC) and 5.148 eV (NA-FF). As expected, the DDH results are closer to experimental measurements compared with those of the PBE0 functional since the fraction of exact exchange is chosen according to the system-specific dielectric constant. Overall, we find that computing electron−phonon interactions at the hybrid level of theory is a promising protocol to obtain quantitative results, comparable to those of experiments. We note that ref 19 reported −262 and −278 meV for the indirect gap renormalization of diamond obtained with the PBE functional, a 3 × 3 × 3 supercell, and the projector augmented-wave (PAW) method, 75 using one-shot 76 and standard 20,77 Monte-Carlo (MC) approaches, respectively. These results are in good agreement with our result of −281 meV obtained with the AHC method at 0 K. Reference 19 also reported −320 and −315 meV with a 5 × 5 × 5 supercell, using one-shot and standard MC approaches, respectively, indicating that the full converged result is about 10% larger in absolute value relative to what obtained with a 3 × 3 × 3 supercell.

APPLICATION TO SPIN DEFECTS IN DIAMOND
Spin defects have been extensively studied due to their potential applications in quantum technologies. 78−81 To accurately predict the electronic structures of spin defects, we computed their electronic properties using electron− phonon renormalizations, and we considered a single-boron defect and the NV − center shown in Figure 4. The calculations were carried out in a 2 × 2 × 2 cubic cell with a 60 Ry plane wave energy cutoff and SG15 62 ONCV 63 pseudopotentials (63 atoms for the NV − center and 64 atoms for the single-boron defect).
In Tables 9 and 10, we report the electronic energy levels and ZPR for both defects obtained with the PBE, PBE0, and DDH functionals. We found that electron−phonon inter-  Energy levels are referred to the HOMO energy level, and the labels of energy levels are given in Figure 4c. Energy levels are referred to the HOMO energy level, and the labels of energy levels are given in Figure 4d.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article actions weakly affect the energy levels of the NV − center, which exhibit localized wavefunctions; they are instead more significant for the single-Boron defect with delocalized wavefunctions. In the NV − center, the ZPR of the LUMO computed with the PBE functional is only −35 meV and that of the HOMO is negligible. In addition, the hybrid functionals PBE0 and DDH yield results similar to PBE results. For the boron defect, with the PBE (DDH) functional, the ZPRs of HOMO and LUMO are 111 meV (121) and 126 meV (241), respectively.

CONCLUSIONS
In this paper, we computed phonon frequencies and electron− phonon interaction at the level of hybrid DFT by using DMPT and by solving the Liouville equation. Using this approach, we obtained phonon frequencies and energy gap renormalizations for molecules and solids by evaluating the non-adiabatic full frequency-dependent electron−phonon self-energies, thus circumventing the static and adiabatic approximations adopted in the AHC formalism, at no extra computational cost. We investigated the electronic properties of small molecules using LDA, PBE, B3LYP, and PBE0 functionals. We also carried out calculations of the electronic structure of diamond with the PBE, PBE0, and DDH functionals and found that the hybrid functionals PBE0/DDH noticeably improve the renormalized energy gap compared to experimental measurements. In addition, we studied the electron−phonon renormalizations of defects in diamond, and we concluded that electron− phonon effects are essential to fully understand the electronic structures of defects, especially those with relatively delocalized states. We note that the use of hybrid functionals affects both the perturbing potentials and wavefunctions entering the definition of the electron−phonon coupling matrices. In general, we expect the wavefunctions computed with hybrid functionals to be of better quality compared to those obtained with LDA/GGA. The perturbing potentials computed with hybrid functionals are non-local compared to the ones obtained with LDA/GGA. Thus, we expect both differences in wavefunctions and potentials, relative to LDA/GGA, to impact calculations with hybrid functionals.
In conclusion, computing electron−phonon interactions at the hybrid functional level of theory is a promising protocol to accurately describe the electronic structure of molecules and solids, and DMPT is a general technique that allows one to do so in an efficient and accurate manner by evaluating non adiabatic and full frequency-dependent electron−phonon selfenergies. orcid.org/0000-0002-8001-5290; Email: gagalli@ uchicago.edu