Ab initio calculation of electron-phonon linewidths and molecular dynamics with electronic friction at metal surfaces with numeric atom-centered orbitals

Molecular motion at metallic surfaces is affected by nonadiabatic effects and electron-phonon coupling. The ensuing energy dissipation and dynamical steering effects are not captured by classical molecular dynamics simulations, but can be described with the molecular dynamics with electronic friction method and linear response calculations based on density functional theory. Herein, we present an implementation of electron-phonon response based on an all-electron numeric atomic orbital description in the electronic structure code FHI-aims. After providing details of the underlying approximations and numerical considerations, we present significant scalability and performance improvements of the new code compared to a previous implementation [Phys. Rev. B 94, 115432 (2016)]. We compare convergence behaviour and results of our simulations for exemplary systems such as CO on Cu(100), H$_2$ adsorption on Cu(111), and CO on Ru(0001) against existing plane wave implementations. We find that our all-electron calculations exhibit faster and more monotonic convergence behaviour than conventional plane-wave-based electron-phonon calculations. Our findings suggest that many electron-phonon linewidth calculations in literature may be underconverged. Finally, we showcase the capabilities of the new code by studying the contribution of interband and intraband excitations to vibrational linewidth broadening of aperiodic motion in previously unfeasibly large periodic surface models.


Introduction
Electronic and vibrational dynamics in metallic systems can be strongly coupled [1,2], due to the insufficient energy scale separation of low-lying electronic excitations and vibrational motion.This also holds true for the reactive dynamics of molecules adsorbed at metallic surfaces [3].Due to the infinitesimal bandgap present in metallic surfaces, adsorbate dynamics at metal surfaces may be affected by nonadiabatic (NA) effects with experimentally measurable consequences.This includes the ultrafast decay of adsorbate vibrational excitation, measurable as picosecond-scale vibrational lifetimes of molecular adsorbates [4,5].Ultrashort adsorbate lifetimes have been reported for systems such as carbon monoxide (CO) on Cu(100) [6][7][8][9], and on Pt(111) [10].Atomic and molecular beam scattering experiments of hydrogen atoms and diatomics such as nitric oxide exhibit highly vibrationally inelastic scattering behaviour that arises from electron-hole pair excitations due to atomic and molecular impingement [11][12][13].Equally, ultrafast time-resolved spectroscopy measurements have shown efficient activation of molecular motion and desorption driven by hot electron excitation [9,14].Beyond capturing the fundamental dynamics at surfaces, hot electrons are also relevant in the context of light-controlled chemistry and catalysis [15,16].
The accurate and efficient simulation of electron-phonon coupling (EPC) is crucial to provide insights on vibrational lifetimes and linewidths [17], critical superconducting temperatures and electron-phonon instabilities in phonon spectra [18].Several software packages have been established for the calculation of Conceptual graphic.FHI-aims software can calculate the position (R) and energy (hω) dependent electronic friction tensor, Λ q (hωqν ) using DFT.The tensor can be used in molecular dynamics with electronic friction (MDEF) (we show a molecule scattering from a metal surface) or projected along normal mode displacements uqν to calculate phonon linewidths.We mostly focus on the latter in this work, but the implementation serves both applications.
EPC effects in condensed matter based on pseudopotential plane-wave formalisms (e.g.Quantum ESPRESSO (QE) [19], EPW [20] or YAMBO [21]).Most existing approaches perform perturbative expansions starting from the harmonic phonon spectrum, i.e. they calculate the frequency renormalization or linewidth associated with EPC [18,22].In order to deliver converged and robust results, very dense Brillouin zone integration for both electrons ('k-space') and phonons ('q-space') is required, which is often achieved by transformation to a maximally localized Wannier basis for interpolation of coupling matrix elements [23,24].Most phonon linewidth calculations in literature have been calculated following this approach based on the simplified first-order perturbation theory expression introduced by Allen and coworkers (vide infra) [25].
The simulation of EPC and NA energy dissipation effects in chemical dynamics at molecule-metal interfaces requires molecular surface models that correctly capture physical surface coverage and bonding scenarios, and that explore dynamics well beyond the harmonic regime.In the case of low coverage adlayers or catalytic models, this can lead to systems featuring hundreds of atoms for which the dynamical phase space needs to be explored.This provides a significant computational challenge for existing plane-wave implementations.Only recently have atomic-orbital based EPC implementations emerged [26,27].Atomic orbital-based codes offer the benefit of a sparse and atom-centred basis, which leads to compact matrix representations of the Hamiltonian.The FHI-aims code is an all-electron full-potential-based numeric atom-centred orbitals (NAO) implementation of electronic structure theory based on density functional theory (DFT) [28].Recently within FHI-aims, Shang et al have developed a linear response implementation for the calculation of phonons and electric field response [29].FHI-aims is highly scalable on high-performance computing architectures due to the use of distributed matrix parallelism as provided by the eigenvalue solvers for petaflop applications library [30] that is part of the electronic structure infrastructure (ELSI) project [31].
EPC effects on dynamics can also be represented in an anharmonic or Cartesian representation as an alternative to describing EPC effects in the harmonic regime based on phonons.This is much less common but has the benefit of enabling an anharmonic treatment of transport and dissipation effects [32][33][34].One such approach is the MDEF method, where the dynamics of nuclei are represented via classical equations of motion (see MDEF in figure 1, top) [4, 35,36].In MDEF, forces on the nuclei derive from the potential energy landscape that arises from the ground-state electronic structure and additional force contributions that relate to dissipative interactions due to hot electron excitations in the metal.The central quantity that governs these forces is the friction tensor Λ, which is a function of atomic coordinates and can be directly calculated from EPC based on first-order perturbation theory [37].Maurer et al have presented an implementation of electronic friction based on linear response calculation of EPC in FHI-aims [38].This implementation has enabled the calculation of on-the-fly MDEF simulations based on DFT for reactive scattering of H 2 on Ag(111) and has subsequently been used to construct machine-learning-based representations of the electronic friction tensor for state-to-state scattering of H 2 on Ag(111) [32,39,40] and NO on Au(111) [41].The electronic friction tensor can also be employed to directly calculate EPC-based vibrational linewidths equivalent to the conventional approach, by simply projecting the Cartesian friction tensor into normal mode space as shown in figure 1, where ũqν refer to mass-weighted phonon displacement eigenmodes associated with phonon frequency ω qν [38].
In this work, we present a detailed introduction and assessment of ab initio electronic friction calculations based on an all-electron NAO implementation of DFT.We present a refactorization of the linear response EPC and electronic friction code in FHI-aims that delivers a substantial increase in computational efficiency and scalability.This code enables the calculation of EPCs, electronic friction tensors, and EPC-based vibrational linewidths.Details of the underlying theory, the necessary approximations and numerical choices are provided.We show that the code delivers reliable and converged predictions with improved convergence behaviour compared to existing plane wave pseudopotential EPC implementations.The improved scalability of the code enables the study of periodic systems with more than 250 atoms per unit cell, which we showcase by calculating the inter-and intraband contributions to the vibrational linewidth of aperiodic internal stretch (IS) vibration in a coherent c(2 × 2) CO monolayer adsorbed at Cu(100).

NA effects on lattice vibrations in metals
Many-body perturbation theory provides a general framework to understand how EPC and nonadiabaticity affect phonon properties [18,42].The harmonic phonon spectrum is based on adiabatic phonon frequencies ω qν that are calculated from DFT via the dynamical matrix of the system [43].This approach already considers some effects of electron-electron correlation on the frequency of the phonon.The phonon is characterized by two quantum numbers, namely ν, which refers to the band index of the phonon and q, which refers to the crystal periodic momentum vector.Phonon frequencies correspond to the poles on the frequency axis of the adiabatic one particle phonon Green's function, D A (q) [18,22,42]: where η is a positive real infinitesimal regularization constant.The phonon spectrum is given by the spectral function, Each frequency corresponds to a delta resonance of infinite lifetime.The interacting phonon Green's function D(q, ω) can be written as [42,44]: The effect of EPC is described by the NA electron-phonon self-energy Π NA ν (q, ω).EPC can lead to a shift (renormalization) of phonon frequencies, coupling of phonons, and to linewidth broadening by introducing a finite lifetime to the phonon.The real part of Π NA ν (q, ω) is associated with frequency renormalization and the imaginary part with linewidth broadening and lifetime of the phonon.The interacting Green's function D(ω) has complex poles Ωqν of D −1 (ω) (i.e.frequencies at which D −1 = 0) which have a real and a complex contribution Ω2 qν = Ω 2 qν − iγ qν .The effects of Π NA ν (q, ω) on the phonon spectrum are complex and involve mode coupling and multiple poles for the same frequency, leading to new structures in the phonon spectrum.However, the spectral function of the NA interacting phonon spectrum is typically approximated as a sum of quasiparticle peaks represented by a Lorentzian centred around the renormalized frequency Ω 2 qν with a full width at half maximum (FWHM) of γ qν [18]: Many approximations are applied in practical applications.Most practical EPC linewidth calculations only consider a subset of effects, for example, they ignore mode coupling and multiple pole structures by making the Lorentzian ansatz shown in equation (4).The EPC phonon FWHM, γ qν is most often evaluated only up to first-order in EPC with perturbation theory [18].EPC-based vibrational linewidth calculations are typically performed in the single particle picture given by independent Kohn-Sham electrons within DFT.Herein, certain approximations are implicit, e.g.ignoring vertex corrections, replacing dielectric matrix with random phase approximation plus exchange correlation ('RPA+xc') response obtained from DFT and neglecting frequency dependence of EPC matrix elements, g (more details can be found below and in [18]).
Calculated EPC linewidths can be directly compared with experiments, which often report FWHM broadening, vibrational relaxation rates, or vibrational lifetimes.The vibrational lifetime is related to the FWHM via the Heisenberg uncertainty: τ qν = h/γ qν .The relaxation rate in inverse time is defined as Γ qν = γ qν /h.

The electronic friction tensor
A common practical expression for the vibrational linewidth of a phonon, γ qν , induced by first-order EPC is defined in equation ( 5): This expression can be obtained from the imaginary part of the first-order NA phonon self-energy as shown in equation (4) [18,44] or can be derived from Fermi's golden rule [38,45].The EPC matrix elements g correspond to screened EPC.Equation ( 5) covers only interband excitations for the case of q = 0 (m > n) and interband and direct intraband excitations for the case of q ̸ = 0 (m ⩾ n).Due to momentum conservation in periodic systems, contributions with q ′ ̸ = q vanish at first-order.Higher-order electron-phonon effects (1-electron, 2-phonon scattering) can lead to additional intraband contributions, but they are not covered by this expression [22,44].We will indirectly incorporate such effects in section 4.3 by band folding in larger unit cells.In this and the following expressions, we explicitly note the sum over spin, σ.The integration of the Dirac delta in equation ( 5) requires careful convergence with respect to Brillouin zone sampling.It is common to rewrite equation ( 5) by neglecting the phonon energy in the delta and by taking the low-temperature limit as proposed by Allen [25].In most practical applications, vibrational broadening is calculated with this low temperature approximation equation ( 6): where ϵ F is the Fermi energy.The appendix of [38] shows how to arrive at equation ( 6) from equation (5).Equation ( 6) is positive definite, hence easier to converge numerically than equation (5).The disadvantage is that the temperature dependence is lost, and that one cannot resolve fine features on the scale of the phonon energy.Calculations with equation ( 5) are somewhat more demanding and have been less commonly reported [23,38,46,47].For example, Lazzeri and Mauri have used equation (5) to predict the temperature dependence of phonon line broadening in graphite and graphene [48].
The matrix elements that feature in equations ( 5) and ( 6) have units of energy and are defined as follows: with The EPC matrix elements in equation ( 8) have units of energy per length.They describe the excitation of an electron from a state nk to a state mk + q by absorption of a phonon qν.V in equation ( 8) is the self-consistent ('screened') effective potential from a Kohn-Sham DFT calculation.Its derivative with respect to atomic displacement ∂ qν V, is defined as Equation ( 9) introduces the elements of the phonon displacement eigenmode u aκ qν associated with phonon qν.Here, a is an atom index and κ indicates the 3 Cartesian directions x, y, and z.Note that we do not include mass weighting as the mass weighting is considered in equation ( 7) via M ν .By inserting equation (9) into equation (8) we find On the right-hand side of equation (10) we introduce a vector notation to denote the multiplication of the normal mode displacement vector u qν with the vector of atomic contributions to the perturbation potential gmn (k, q).
We can find an expression for the Cartesian friction tensor by inserting equations ( 7) and (10) into equation (5): The factor of 2 in this expression arises from the definition of FWHM in equation (4) and is cancelled by the factor 2 in the definition of the EPC matrix elements [49].We can bring the normal mode displacement vectors u qν outside of the sum over states, and define the mass-weighted normal mode displacement vectors ũqν = u qν / √ M ν .Herein, we have defined the Cartesian friction tensor [38] Λ A similar result for the phonon linewidth can be generated on the basis of atomic response correlation functions [35,49].Additionally, we define a related expression for the Cartesian friction tensor, with a substition that is valid for true delta functions: We can follow a similar derivation starting from the low temperature double delta expression by Allen in equation ( 6) and arrive at: These expressions of electronic friction are generalizations of the expressions defined by Maurer et al [38] for arbitrary q values and both have units of [mass/time] or [(energy×time)/length 2 ] or [action/length 2 ].
Litman et al derived an alternative expression of friction tensor based on the Laplacian transform of the spectral density that does not rely on artificial broadening during the Fermi surface integration [50]: where Ω nm = (ε n − ε m ) /h.For this expression, there is no dependence on an ad hoc broadening width parameter.In sections 4.1 and 4.2 we compare the 3 expressions of the Cartesian friction tensor above for the evaluation of vibrational linewidths of CO/Ru(0001) and the electronic friction tensor of H 2 /Cu(111).A similar analysis has been carried out previously by Calandra et al for the phonon linewidths of MgB 2 [51].
We revisit this analysis for the case of molecules adsorbed on surfaces using finite difference and localised atomic orbitals rather than DFPT and plane waves.

EPC matrix elements and expressions in local atomic orbital basis
Several previous works have discussed the implementation of EPC matrix elements in local orbital basis [24,38,52,53].We will provide a brief summary of the main expressions below.We can find expressions for the EPC matrix elements gaκ mn (k, q) by expanding the electronic eigenstates mk in a local orbital basis of crystal periodic atomic orbitals ϕ ik : Inserting equation ( 17) into equation ( 8) yields We start with h aκ ij (k, q) defined in equation ( 18) and simplify notation.We first recognize that the perturbation of the effective potential is equivalent to the Hamiltonian response: The matrix elements h aκ ij (k, q) define a first-order Hamiltonian response matrix H (1) that can be calculated in FHI-aims via density functional perturbation theory (DFPT) [29] or finite difference response [38].In this work, we will only consider finite-difference-based calculations.While h aκ ij (k, q) can be directly calculated with DFPT, in the case of finite difference evaluation of the response with consecutive self-consistent-field calculations for displaced atoms, we only have access to As previously shown, following relationship between the derivative of the expectation value and the derivative of the operator holds [35,38,52]: We can now insert equation ( 21) into equation ( 18): Using the fact that the wave functions satisfy a generalized eigenvalue equation, we can make following replacement: The second and third term on the right hand side correspond to the left and right sided derivative components of the overlap matrix S aκ,L ij,k and S aκ,R ij,k : With this, we come to following expression of the EPC matrix elements of equation ( 18): Comparing with the findings of Maurer et al [38], we can relate the potential derivative (deformation potential matrix elements) and NA coupling matrix elements as follows: As shorthand, we define an EPC element G aκ ij (k, q) (contents of brackets in equation ( 25)) in local orbital representation as In previous works, for computational convenience, approximate versions of these matrix elements have been proposed by Head-Gordon and Tully [35,38], where the dependence on individual electronic states is dropped and the linear response of the full overlap matrix is used instead: Maurer et al have assessed the impact of this approximation on vibrational relaxation rates for aperiodic systems and have found it to affect final predictions by only few percentage points [38].For lattice periodic (q = 0) perturbations in metallic systems (specifically where fractional occupations are present) the response of the Fermi level must additionally be accounted for: This applies only to intraband elements, and thus is only relevant for Allen's formula (equation ( 15)) where q = 0 intraband elements contribute and not the single delta expression of type equation ( 14) where they do not.

Delta function approximations
The numerical evaluation of the friction tensor involves either two delta functions (low temperature limit, equation ( 15)) or a single delta function (equation ( 13)).These delta functions cannot be reliably evaluated directly for ab initio calculations due to the discrete nature of the density of states represented on a finite k-point mesh.Typically the delta function(s) are replaced by a broadening function of finite width, common choices include Gaussian and Lorentzian functions.The most common choice is a Gaussian function: or as defined in QE: whereas a Lorentzian broadening function can be defined as follows: An additional normalization technique was previously proposed for single delta function calculations when relatively large smearing widths are applied [38].This is to account for the contribution of the Gaussian that falls below the integration window: We will reassess the effect of different choices of broadening in section 4.2.A slightly different approach is to replace the last term in equation ( 14) (including the reciprocal energy term) with the derivative of the smearing function [22,54], using a Gaussian as an example: The choice of smearing parameter can have a large effect on the vibrational linewidths and electronic friction tensor [22,38], by including contributions from a wider range of electron-hole pair excitations.For the calculation of vibrational linewidths, the delta function(s) are typically evaluated within the weak-coupling limit, whereby σ is chosen to be as small as possible whilst being converged with respect to computationally feasible k-meshes.However, MDEF has been most commonly applied to the simulation of atomic/molecular beam scattering experiments [11,41,55], where impinging atoms and molecules with high translational energy and possible high vibrational and rotational excitation are fired at metallic surfaces.For this case, the appropriate choice of σ is unclear, as the coupling can now be much stronger than when calculating vibrational linewidths.A common choice of σ for evaluating electronic friction is within 0.1 − 0.6 eV [38,54,55].A recent study by Spiering et al [56] justifies their use of a σ value of 0.60 eV on the basis of three main arguments: the first, originally proposed by Maurer et al [32], is that σ is less than the incidence energy of the molecule and the second, proposed by Novko et al [22], is that realistic values of broadening can be obtained from the electrical resistivity by estimating the relaxation time of the electronic-hole pairs.Novko et al [22] calculates an appropriate value of ⩽0.15 eV for CO/Cu(100), Spiering et al [56] argues that since ruthenium has a 4 times larger value of resistivity than copper, a value of 0.6 eV (4 times larger than 0.15 eV) for σ is realistic for the N 2 /Ru(0001) studied.This implies that previous studies Graphical depiction of (a) aperiodic motion of adsorbate intramolecular stretch motion within a periodic surface overlayer, (b) periodic motion with q = 0 and (c) with |q| = qmax.Crystal periodic motion is characterized by a momentum vector q and the associated perturbing potential can be represented by a phase modulated sum over displacement contributions of atom images in all periodically repeating unit cells.Aperiodic motion breaks translational invariance and can therefore only be approximated with an infinite Fourier sum over q.
might have employed a σ value that was too high for copper.Also, the resistivity of platinum is a factor of ∼10 larger than that of copper, thus should imply a σ of over 1.0 eV is realistic.The last argument is that the results do not vary significantly for σ values of 0.15, 0.25 and 0.60 eV.In this work, we will explore the effect of the value of σ for different expressions for calculating electronic friction.

Supercell approach
Our aim is to be able to describe vibrational energy dissipation during surface dynamics in the context of heterogeneous catalysis, where typically individual adsorbates react at the surface and excite hot electrons.This is also the most relevant case for modelling molecular beam scattering experiments.We therefore differentiate between aperiodic and periodic motion on the surface in the context of periodic slab calculations.Figure 2 shows a graphical depiction of the aperiodic IS vibrational motion of a single adsorbate molecule in a periodic overlayer, the coherent totally symmetric motion of the q = 0 IS phonon, and the totally antisymmetric motion of the phonon with |q| = q max .
The perturbing potential due to aperiodic motion (figure 2(a)) of an atom in the principal unit cell can be described as a Fourier sum over all q vector contributions to the potential in the principal unit cell N = 0: Let's consider coupling matrix elements that arise due to this aperiodic perturbation: Note that summation over q always implies normalization weighting with corresponding Brillouin weights w q .By moving the sum over q outside of the friction tensor description, we can see that the calculation of the vibrational relaxation rate for an aperiodic vibrational motion in the principal unit cell can be described as q-averaging over the relaxation rates of all q-modes in the irreducible Brillouin zone: where N q is the number of q points.As the code is currently restricted to q = 0 excitations, we perform this sum by increasing the size of supercells to include q > 0. This effectively folds intraband excitations into the principal Brillouin zone, which increases the number of available interband excitations.The improved memory distribution in FHI-aims allows such systems to be calculated and the improved scaling of the code allows such cells to be tackled within the walltime of typical HPC systems.Due to the lack of a MDEF framework to deal with q > 0 matrix elements and the typical expense of q > 0 calculations it is common to instead calculate friction associated with coherent motion of the overlayer (i.e q = 0 case, see figure 2(b), which should be similar to aperiodic motion in the limit of low adsorbate coverage.In previous work [41], we investigated the dependence on the q = 0 electronic friction tensor with respect to unit cell size from a p(2 × 2) coverage to a p(5 × 5) for NO/Au(111) in the context of scattering and found little change.

Implementation
The workflow for the evaluation of the friction tensor as implemented in FHI-aims is shown in figure 3. The electronic friction calculation begins with a ground state self-consistent field (SCF) DFT calculation until convergence is reached, this is required as the ground state eigenvectors and chemical potential (Fermi energy) are used in the calculation of the response.The calculation of the required response matrices (typically the first-order Hamiltonian and overlap matrices) follows.These can be read in from a prior calculation using input and output routines provided via the ELSI library, or calculated using finite difference.The subset of atoms to be included in the electronic friction calculation can be specified by introducing the keyword friction_atom after the definition of the Cartesian position of the atom in the geometry.ininput file.This means that the response for this atom in each Cartesian direction will be calculated in a new SCF cycle for both forwards and backwards directions.
If the analytical matrix elements G nm,( 1) ij (equation ( 27)) are required rather than the approximate G HGT ij (equation ( 28)), a calculation of the one-sided derivative of the overlap matrix follows.Now the coupling matrix can be obtained, the type depends on the choice of non-adiabatic coupling matrix elements.Finally, the tensor is evaluated and printed.During the tensor evaluation, several calculation parameters can be selected, including the value of σ for the broadening of the delta function, the value of the temperature that In (b) we use a larger system (with 64 k-points), which was not calculable with the original code.Calculations performed on Warwick SCRTP Avon HPC (Dell PowerEdge C6420 compute nodes each with 2 x Intel Xeon Platinum 8268 2.9 GHz 24-core processors; 48 cores per node; 180 nodes).
defines the orbital occupations, the choice of delta function, the type of friction tensor expression to be used (equation (13) or equation ( 15)), and the frequency at which the friction tensor is evaluated (if equation ( 13) is used).The computational performance and scalability of the code is described in section 3.2.

Code scalability and performance
The electronic friction code has been substantially refactored since it was originally published [38].The response matrices now support ScaLAPACK-type matrix distribution, allowing efficient distribution and scaling when N cores > N k .This allows much larger systems to be investigated, the largest system calculated in this work comprises 252 atoms, 16 k-points and 108 response coordinates (3 Cartesian coordinates per 36 response atoms).The original code would not be able to treat this system, or any system with more than approximately 10 atoms for which the response is calculated at a meaningful level of accuracy.In figure 4(a) the improved scaling behaviour of the refactored code (version 210928) is compared to the original code (version 190506).The refactored code significantly outperforms the original code even when N cores < N k , specifically for the friction tensor evaluation.The scaling of the original code saturates with cores when N cores = N k , whilst the refactored code continues to scale beyond this point.The scaling remains favourable for larger systems that were not feasible with the original code, as shown in figure 4(b).The friction tensor evaluation remains on a similar timescale to a single SCF iteration, and shows similar scaling with the number of cores.Key benefits of having an internal finite difference code (rather than calling FHI-aims externally multiple times via post-processing) is that (i) the ground-state wave function can be used as an initial guess for finite-difference displacements to reduce SCF iterations for all finite difference displacements, (ii) time consuming I/O is kept to a minimum, and (iii) evaluation of the friction tensor runs internally in parallel.Whilst the full calculation of the friction tensor is still too computationally costly for large-scale on-the-fly ab initio MDEF simulations, the code improvements will dramatically enhance our ability to gather data for the construction of machine learning models of the electronic friction tensor and the electron-phonon response [57].

FHI-aims computational details
Unless stated otherwise, all calculations are performed using the all-electron NAO code FHI-aims (2021 development version) [28] without spin polarization and with PBE exchange-correlation functional, [58] dipole correction, atomic ZORA relativistic correction, tolerances of 10 −6 eV, 10 −3 eV and 10 −5 e / a 0 3 for the total energy, the eigenvalue energies and the electronic density, respectively.Unless stated otherwise, calculations are performed with standard 'tight' basis set definitions and were evaluated using finite difference (displacement of 0.001 Å).We employ the approximate EPC matrix elements given in equation ( 28) as we find no significant difference to the matrix elements given in equation (27).
For the benchmarking of the code, we have aimed to reproduce the results given in reference [55] and reference [14] with the FHI-aims (2021 developmental version) code.Our settings are chosen to be as close as possible to the settings discussed in the two literature studies.The Gaussian smearing width employed in the Fermi surface integration is also used for all evaluations of occupations and Fermi level in the ground state SCF and finite difference calculations.
Each figure caption contains a DOI that links to a dataset in the NOMAD materials repository where all output files are available for download.

QE computational details
For the QE results shown in section 4.2, we have aimed to reproduce the results given in reference [55] Our settings are chosen to be as close as possible to the settings reported in the two literature studies.QE (version 6.8) [59] is employed using PW91 [60] exchange-correlation functional with ONCV pseudopotentials [61], plane-wave energy cutoff of 816 eV with p(2×2) Cu(111) slabs containing 4 layers of metal and a vacuum height of 10 Å.We utilized the k-mesh interpolation scheme as employed in the PHonon QE module to interpolate a 6 × 6 × 1 to an 18 × 18 × 1 k-mesh.The delta function was approximated with a standard Gaussian broadening (equation ( 30)) of 0.6 eV.
QE operates in a reduced normal mode framework.To convert the outputted vibrational mode linewidths to Cartesian friction tensor elements we project using the normal modes produced by QE.

Cartesian to internal coordinate transformation
Friction tensor elements are calculated in Cartesian coordinate framework, however two elements are shown in internal coordinates (namely Λ dd and Λ ZZ ), employing Λ int = U T Λ cart U for the transformation from Cartesian coordinates, where:

Numerical convergence of coherent phonon linewidths of CO/Ru(0001)
In this section, we compare the ability of different Fermi surface integration expressions described above to predict the q = 0 linewidth of different phonon modes of CO adsorbed on Ru(0001) in a p(1 × 1) coverage, and suggest the most appropriate expression.The linewidths predicted by different expressions as a function of smearing width are shown in figures 5(a)-(d) for the 3 types of modes shown in figure 5(e) (noting that there are 2 frustrated rotation (FR) modes).We also show results previously reported by Diesen et al using the well-established QE package as validation for our implementation [14].QE employs an Allen's formula expression (equation ( 6)), thus exhibits an asymptotic behaviour at σ → 0 and is dominated by intraband contributions.Allen's formula is the most common expression employed for phonon linewidth calculations but, as previously discussed in literature, this neglects temperature and frequency dependence whilst giving inflated linewidths at q = 0 [51].FHI-aims results employing Allen's formula (Λ II ) provide similar linewidths as Diesen et al, thereby validating the FHI-aims implementation for evaluating phonon linewidths at q = 0. Conversely, the Λ Ia expression is dominated by interband excitations and so approaches 0 for σ → 0. The Λ Ia expression quickly increases with σ, as expected as this represents a stronger coupling regime.The effect is weaker for the IS mode which is at a higher perturbing frequency, where the coupling is already stronger for q = 0 thus 'smearing in' further electron-hole pair excitations has a smaller effect.However, it is difficult to converge calculations with a small value of σ, thus typically other expressions are employed that have a weaker dependence on σ, allowing reasonable values to be obtained at a value of σ that can be converged with a smaller number of k points [38].
The Λ Ib expression excludes intraband excitations for q = 0 and is very similar to Λ Ia for σ → 0, but displays much less dependence on σ and reaches a plateau much faster than any other expression (at σ = 0.1 eV).This means that it is easier to predict converged and robust results in the weak coupling limit with this expression than with the original Λ Ia expression.At high σ, a similar linewidth is predicted for all modes and for all expressions except for Λ Ia , which generally increases drastically with increasing σ.Λ III does not depend on an ad hoc broadening parameter, and so is plotted as a horizontal line.Interestingly, for the IS and SA modes, Λ III is similar to the σ → 0 value for Λ Ia and Λ Ib .
We now compare the convergence of the different expressions for calculating vibrational linewidth with respect to the size of k-mesh in figure 6.For Λ Ia and Λ Ib with σ = 0.01 eV, reasonable convergence is only achieved for the IS mode (i.e the mode with highest frequency) for the range of k-meshes investigated.Λ Ib with σ = 0.3 eV appears well converged for all k-meshes and modes but appears to converge at a smaller linewidth for IS and a larger linewidth for FR modes compared to the Λ Ia and Λ Ib with σ = 0.01 eV.Λ II at σ = 0.3 eV remains fairly stable with respect to the number of k points, though the inflation of the friction due to intraband excitations is particularly extreme for the SA mode.Λ III achieves convergence, and interestingly, is similar in magnitude to Allen's formula expression with σ = 0.3 eV for all modes except the SA mode, where the latter is inflated by intraband contributions.
In figure 7, we compare the vibrational linewidth convergence with respect to the number of substrate layers for the different Fermi surface integration expressions.The magnitude and the convergence behaviour differ significantly depending on the expression used to evaluate the vibrational linewidth.It is clear that some expressions were not converged with respect to substrate thickness for the 6 layer cell employed in the comparisons above, including Λ Ia where a small σ value does not achieve adequate convergence for any mode.Conversely, the Λ Ib expression with 0.3 eV smearing appears to converge the earliest with respect to substrate thickness for the IS and SA modes while still predicting values close to Λ Ia for all modes.Values predicted by Λ II appear more sensitive to layer convergence than Λ Ib , which accounts for most of the discrepancy between the two expressions (including that observed in figures 5 and 6).However, Λ II differs significantly from Λ Ib for the SA mode even at 11 layers.For the FR modes, there is little difference in the substrate thickness convergence behaviour between Λ Ib and Λ II .Λ III converges fast for the high frequency IS mode, but is unstable even at 10 layers for FR modes.
Λ Ib in the low broadening limit recovers Λ Ia for all substrate thicknesses, even for the IS mode with the highest perturbing energy.For larger σ, Λ Ib converges with respect to substrate thicknesses earliest and provides similar values to Λ Ia .Thus we conclude that, of the investigated expressions, Λ Ib enables a more robust computation of linewidths compared to Λ Ia or Λ II , the latter of which appears to overestimate the SA linewidth.Numerical convergence can be achieved with σ values that are slightly larger than what can be used with Λ Ia yet without distorting the results too much.As a result, convergence can be achieved with a smaller k mesh.

Electronic friction during dissociation of H 2 /Cu(111)
Our second investigated system is the minimum energy reaction pathway (MERP) of dissociative H 2 adsorption on Cu(111) as shown in figure 8(a), where we now study two elements of the electronic friction tensor, namely the diagonal elements for the IS (Λ dd ) and the centre of mass (COM) in the Z direction (Λ ZZ ).We first compare results produced by FHI-aims (using equation (15)) and QE to prove that our implementation provides predictions that are consistent with existing software packages.We further compare to previous VASP and QE results published by Luntz and Persson [62] and Spiering and Meyer [55], respectively, in figure 8(b) for the two aforementioned elements.All results for both elements show a general qualitative trend of maximum friction at the transition state (TS) with a drop in friction for the chemisorbed product, except for Luntz's results for Λ ZZ , which show a drop at the TS.The agreement between our FHI-aims and QE results is generally excellent, except for the TS for both elements and at ≈ −0.8 Å for Λ dd .We attribute the discrepancy of the two codes at the TS to the extreme sensitivity of the friction coefficients at the saddlepoint.Some discrepancy is expected since the MERP was optimised using FHI-aims and employed for both the FHI-aims and QE results that we produce.Our results agree well with the previously published Luntz and Spiering results for Λ dd except at the TS for FHI-aims and ≈ − 0.8 Å for QE.However, for Λ ZZ , our values calculated with FHI-aims and QE are generally 4 times smaller than literature results.We can only attribute the differences in Λ ZZ to the internal coordinate transformation, as we have shown that the FHI-aims and QE results where we apply the same transformation are otherwise consistent.
Having proven the validity of our implementation, we now compare the convergence behaviour of different expressions for evaluating the electronic friction tensor to be used within the Markov MDEF scheme, calculated at zero perturbing energy.We investigate the geometry at reaction coordinate −1 from the MERP of H 2 /Cu(111) just before the TS as shown in figure 8, and present the electronic friction tensor elements in the internal coordinate framework given in the same figure.The same 2 elements from the previous section are investigated, namely the IS (dd) and z coordinate of the COM (ZZ) elements of the friction tensor.We limit our investigation to the coherent motion (q = 0), which is applicable to the molecular beam scattering case, provided the unit cell is large enough and the coverage is sparse, where significant broadening may be physical as the translational and vibrational energy can each reach several eV [12].However, it is currently unclear what magnitude of this ad hoc broadening parameter is most appropriate.For zero perturbing frequency, Λ Ia → ∞ and Λ III → 0, thus we just compare Λ Ib and Λ II .
In figure 9 we see the general behaviour of the different friction expressions with respect to changing σ.The results can be split into two groups, those that apply the low-temperature Allen's formula expression which tends to ∞ for σ → 0 and those that tend to 0 for σ → 0. It is evident that the Allen's formula expression (Λ II ) inflates the q = 0 friction values for all values of σ.This implies that previous studies employing Λ II and moderate values of σ may overestimate the effect of friction and thus the electron-hole pair dissipation channel in dynamics at surfaces [55,56].All single and double delta expressions approach a  15)) compared to results by Luntz and Persson [62] and Spiering and Meyer [55].The reaction coordinate is given by the negative difference of the height with respect to the COM height of the H2 at the TS.The Gaussian function is defined in equation (31).Use of the delta definition in equation ( 30) slightly changes the results (the largest difference is ≈15% ).σ = 0.6 eV is employed for our results.Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.07-1.

Figure
Comparison of expressions to calculate the dd and ZZ zero-frequency electronic friction tensor elements for H2/Cu(111), at the geometry preceding the transition state from MERP defined above for a range of σ values.The Fermi integration delta functions are replaced with a Gaussian function given in equation ( 30) unless an alternative is given in the legend.4 substrate layers and a k-point mesh size of 18 × 18 × 1 is employed.Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-3.plateau at a similar friction value as σ increases, though the difference remains larger for Λ dd than Λ ZZ .The specific delta function approximation employed has a smaller effect on the friction coefficient for σ > 0.2, with the derivative approximation in equation (34) inflating the friction coefficient the most.The normalisation of the Gaussian due to 'leakage' (equation ( 33)) [38] inflates the friction values with no observed benefit to stability with respect to σ.There is little difference between using Gaussian or Lorentzian broadening for this system.Thus, Λ Ib with either a Gaussian (equation ( 30)) or Lorentzian (equation ( 32)) appears to be the most appropriate expression from those investigated here as they are not inflated arbitrarily.
We now explore the behaviour of Λ Ib and Λ II with respect to changing k-mesh size and the number of layers of the substrate.Figure 10(a) shows both expressions are stable with respect to the investigated k-mesh sizes for σ = 0.3 eV, despite having different converged values.The discrepancy can be attributed to the inflation of electronic friction for q = 0 for Λ II .The friction coefficients are more sensitive to the number of substrate layers as shown in figure 10(b).Interestingly, Λ Ib dd appears most sensitive to the number of substrate layers, approaching the otherwise stable Λ II values with no obvious convergence achieved.The sensitivity of EPC on substrate thickness has been previously reported for CO/Cu(100) [22], and recorded in the experiment as a physical effect [63].Otherwise the Λ Ib ZZ , Λ II dd and Λ II ZZ values appear stable with respect to substrate thicknesses above 5 layers for the investigated σ = 0.3 eV value.
We also see that Λ II differs more between the relaxed 4 layer cell employed in figure 10(a) to the unrelaxed 4 layer cell in figure 10(b) compared to Λ Ib .Generally we find that Λ II can give very different friction values than Λ Ib for this system, though particular attention to the thickness of the substrate should be paid for this system.Neither expression is ideal due to the intrinsic limitations of the Markov approximation of MDEF, which imposes the evaluation of the relaxation rate at zero frequency.On balance, Λ Ib appears to predict the most reliable and robust results.

Interband and intraband contributions to the vibrational lifetime of aperiodic adsorbate motion in c(2 × 2) CO overlayer at Cu(100)
We now use the strong performance of the newly developed code for a showcase application, namely to calculate the vibrational rate of a single CO molecule in a c(2 × 2) overlayer of CO adsorbed at Cu(100).The lifetime of the coherent (q = 0) CO overlayer on Cu(100) has been experimentally studied by IR absorption [64], IR pump-probe [6] and time-resolved vibrational sum-frequency generation spectroscopy [9].In all cases, the coherent IS vibrational mode that is excited promptly decays via dissipation into electron-hole pair excitations and dephasing with other modes.There have been several attempts to theoretically reproduce this shortened lifetime based on calculating only the first-order electron-phonon lifetime of the q = 0 IS mode, [22,38,[65][66][67] all of which yield too large vibrational lifetimes and too small linewidths to be consistent with experiment or neglect other important contributions such as the coupling with substrate phonon modes.More recently, Novko et al provided a comprehensive characterization of IS vibrational energy dissipation for this system, [44,68] which included a description of second order electron mediated phonon-phonon coupling that accounts for the coupling of the IS mode to other modes and the dephasing contribution to dissipation of the q = 0 IS mode.This description is fully periodic and includes how intraband excitations can lead to phonon-phonon coupling.The calculations we present in this paper do not include such contributions directly as we only consider first-order interband excitations, but intraband excitations can be included nominally by band folding when large supercells are used, which our code enables.
In figure 11(a) we show a graphical depiction of various c(2 × 2) CO/Cu(100) supercells of different size, labelled by the number of CO molecules in the cell.The cells range from the primitive cell to a supercell containing 18 CO molecules, resulting in a 108 × 108 friction tensor in the case of the largest cell.We calculate the q-averaged IS lifetime (τ q IS ) as a function of σ for different k-meshes for the cell containing 9 CO molecules in figure 11(b).The 6 × 6 × 1 k-mesh agrees well with a 8 × 8 × 1 k-mesh even for σ = 0.01 eV.For converged k-meshes, τ q IS is fairly independent of σ between 0.04 eV and 0.08 eV, with a weak σ dependence below 0.04 eV.[6] (shaded region is the experimental uncertainty).In (c) a k-point mesh size of N k × N k × 1 with N k = 16, 11, 8, 8, 4 and 3 was employed for the NCO = 2, 4, 8, 9, 16 and 18 cells respectively.6 substrate layers and a perturbation energy of a single CO internal stretch (0.245 eV) in the primitive cell was employed in all cases for the evaluation of friction in equation (14).Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-1.
In figure 11(c), τ q IS is calculated for σ = 0.05 eV as a function of cell sizes (specifically the number of CO molecules present in the cell).We increase/decrease the k-mesh size when decreasing/increasing the cell size to keep the k-mesh density approximately constant and consistent with the N CO = 9 cell.With increasing cell size, τ q IS decreases until reaching convergence with supercell size for the cell containing 8 CO molecules, with around 100 Å 2 surface area.The decrease in lifetime is expected, as q > 0 modes contribute higher relaxation rates, [38] and by increasing the unit cell more such modes are considered.The calculated normal mode frequency for the symmetric IS is 2059 cm −1 for the largest cell, only 5 wavenumbers lower than the experimentally recorded frequency of 2064 cm −1 [6].The calculations for the largest cells are only possible with the improvements to the code and can act as a benchmark for future code development (i.e when including intraband excitations for smaller cells, which would be the preferred strategy when convergence is slow with respect to supercell size.) Whilst the aperiodic motion is not directly comparable to experiment in the case of the vibrational linewidth of c(2×2) CO on Cu(100), it can be seen as an 'instant dephasing approximation' , where we assume that the excited q = 0 phonon instantaneously dephases into aperiodic motion.This may be a fair approximation, as Novko et al [44] reported that q > 0 contributions dominate the IS linewidth.Our predicted lifetime of 1.85 ps is in good agreement with the 2 ps observed in experiment [6].Our supercell results show the utility of a highly efficient implementation and can guide future development of the code, particularly when extending the code towards non-momentum conserving EPC matrix elements (q > 0).

Conclusion
We have presented an all-electron NAO implementation of electronic friction and first-order electron-phonon response within FHI-aims with excellent scalability and memory efficiency that can compete with existing pseudopotential plane wave implementations of EPC.We employ the code for the calculation of vibrational linewidths and electronic friction tensor components of different molecular adsorbate overlayers on metal surfaces.The implementation is documented in detail and is part of the latest development version of FHI-aims and the next release.We showcase the scalability and utility of the code by investigating the convergence properties of vibrational relaxation rates and linewidths due to electron-phonon response.
By comparison with literature and calculations that we perform in the plane wave pseudopotential code QE, we show that our atomic orbital all-electron implementation predicts vibrational linewidths for CO/Ru(0001) and friction coefficients for H 2 /Cu(111) that are in close agreement with existing implementations.We provide a comprehensive analysis of the convergence behaviour for different expressions to calculate the vibrational linewidth of CO/Ru(0001) and the electronic friction tensor of H 2 /Cu(111) in the Markov limit.We find that, of the studied mathematical expression for EPC linewidth calculation, the expression of the friction tensor in equation ( 14) (Λ Ib ) provides the most robust and facile prediction of numerically converged results.
In our calculations, copper surfaces represent particularly difficult cases to achieve layer convergence for electron-phonon response, which likely is related to the low density-of-states of copper at the Fermi level.Finally, we showcase the capabilities of the code by describing aperiodic surface motion in large unit cells of unprecedented size by averaging over inter and intraband contributions to the EPC-induced vibrational lifetime of the IS motion of CO on Cu(100).We show that this 'instant-dephasing approximation' yields a vibrational lifetime that is in good agreement with experiments [6] and previous calculations in the literature that consider higher-order electron-phonon contributions [44].
The new implementation will enable the calculation of EPC-induced vibrational linewidths for complex systems with intrinsic length scales that require unit cells too large to be treated with existing plane-wave pseudopotential implementations.It will further enable on-the-fly ab initio MDEF simulations [32] and the large-scale generation of machine-learning training data to predict the electronic friction tensor as a function of atomic coordinates [57].Alternatively, modern EPC approaches beyond the harmonic phonon picture such as the stochastic self-consistent harmonic approximation [69] may also become viable in this context in the future.We have also shown an approach to calculate electronic friction and EPC linewidths for aperiodic motion within a periodic slab approach.In the future, we plan to extend this implementation to evaluate EPC matrix elements for q > 0 motion, which will likely require real-space interpolation approaches as previously proposed for EPC calculations with maximally localized Wannier functions [23] or atomic orbitals [24].This will be particularly useful where the integration over q-space converges too slowly to be addressed via supercell construction.

Figure 1 .
Figure1.Conceptual graphic.FHI-aims software can calculate the position (R) and energy (hω) dependent electronic friction tensor, Λ q (hωqν ) using DFT.The tensor can be used in molecular dynamics with electronic friction (MDEF) (we show a molecule scattering from a metal surface) or projected along normal mode displacements uqν to calculate phonon linewidths.We mostly focus on the latter in this work, but the implementation serves both applications.

Figure 2 .
Figure 2. Graphical depiction of (a) aperiodic motion of adsorbate intramolecular stretch motion within a periodic surface overlayer, (b) periodic motion with q = 0 and (c) with |q| = qmax.Crystal periodic motion is characterized by a momentum vector q and the associated perturbing potential can be represented by a phase modulated sum over displacement contributions of atom images in all periodically repeating unit cells.Aperiodic motion breaks translational invariance and can therefore only be approximated with an infinite Fourier sum over q.

Figure 3 .
Figure 3. Sequence of calculation operations required during the calculation of the electronic friction tensor.

Figure 4 .
Figure 4. Scaling behaviour for the electronic friction module using finite difference linear response calculations for CO/Cu(100).We compare each SCF step and the friction tensor (Λ evaluation, one per calculation).In (a) we compare the old code (190506, open symbols) with the refactored code (210928, closed symbols) (256 k-points).In (b) we use a larger system (with 64 k-points), which was not calculable with the original code.Calculations performed on Warwick SCRTP Avon HPC (Dell PowerEdge C6420 compute nodes each with 2 x Intel Xeon Platinum 8268 2.9 GHz 24-core processors; 48 cores per node; 180 nodes).
) and where X n , Y n , Z n are Cartesian coordinates of the nth H atom, r = (x 2 − x 1 ) 2 + (y 2 − y 1 ) 2 + (z 2 − z 1 ) 2 and r ′ = (x 2 − x 1 ) 2 + (y 2 − y 1 ) 2 , and M = m A + m B .The column vectors of the transformation matrix are normalised to conserve the trace of the friction tensor during projection.

Figure 5 .
Figure 5.Comparison of expressions to calculate vibrational linewidths of CO/Ru(0001) to results published by Diesen et al [14] for 4 modes at q = 0 as a function of smearing width, shown are the (a) IS, (b) surface adsorbate stretch (SA), (c) and (d) the 2 FR modes.A k-point mesh size of 96 × 96 × 1 with a 1 × 1 cell with 6 substrate layers is employed in all cases.The Fermi integration delta functions are replaced with a Gaussian function given in equation (30) unless an alternative is given in the legend.Graphical depictions of the 3 types of mode are shown in (e).Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-5.

Figure 6 .
Figure 6.Same as figure 5 but as a function of k-mesh size (N k × N k × 1).Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-5.

Figure 7 .
Figure 7.Comparison of expressions to calculate vibrational linewidths of CO/Ru(0001) for 4 modes as a function of number of substrate layers.Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-4.

Figure 8 .
Figure 8.(a) Graphical depiction of the dissociative adsorption of H2 on Cu(111) MERP studied in this work, with COM coordinates (X, Y, Z) and H2 bond length (d) labelled.(b) Two electronic friction tensor elements are plotted across MERP for QE and FHI-aims (using equation (15)) compared to results by Luntz and Persson[62] and Spiering and Meyer[55].The reaction coordinate is given by the negative difference of the height with respect to the COM height of the H2 at the TS.The Gaussian function is defined in equation(31).Use of the delta definition in equation (30) slightly changes the results (the largest difference is ≈15% ).σ = 0.6 eV is employed for our results.Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.07-1.

Figure 10 .
Figure 10.Comparison of convergence behaviour of two expressions to calculate Λ dd and Λ ZZ electronic friction elements for H2/Cu(111), at geometry before transition state in the defined MERP.The convergence behaviour with respect to k-mesh size (N k × N k × 1) is shown in (a) for a relaxed 4 layer cell.The convergence behaviour for unrelaxed cells of varying number of substrate layers (N layers ) is shown in (b), with N k = 16.In all cases σ = 0.3 eV is employed.Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-2.

Figure 11 .
Figure11.(a) Graphical depiction of the CO/Cu(100) cells investigated, labelled with the number of CO molecules present in the cell (NCO).Convergence of the q-averaged IS lifetime (τ q IS ) for (b) the NCO = 9 cell as a function of σ for different k-mesh sizes (legend) and (c) as a function of cell size for σ = 0.05 eV with comparison to the lifetime measured by Morin et al[6] (shaded region is the experimental uncertainty).In (c) a k-point mesh size of N k × N k × 1 with N k = 16, 11, 8, 8, 4 and 3 was employed for the NCO = 2, 4, 8, 9, 16 and 18 cells respectively.6 substrate layers and a perturbation energy of a single CO internal stretch (0.245 eV) in the primitive cell was employed in all cases for the evaluation of friction in equation(14).Raw data used to produce this figure is available in the NOMAD repository under: 10.17172/NOMAD/2023.09.06-1.