A periodic energy decomposition analysis method for the investigation of chemical bonding in extended systems

The understanding in chemistry is most often based on heuristic concepts. The most powerful and most widely used concept is chemical bonding. Classifications such as covalent, ionic, or metallic bonding are central in discussing trends in different compounds and predicting new reactivity. An in-depth understanding of the chemical bond in a system therefore paves the way to predict trends and explore new reactivity. Computational methods are thereby especially helpful to complement experimental determinations of bonding parameters like structure and electron density distribution with quantitative analysis of observable and non-observable quantities.

While the discussion of chemical bonding is a major pillar for molecular chemistry,¹ the efforts to extend the concepts toward periodic structures are scarcer. Still, quantitative analysis of chemical bonding in extended structures leads to greater fundamental understanding of the systems investigated and can help with predicting new trends and reactivity for ab initio material’s design and functionalization of surfaces and interfaces.² The advent of quantum chemistry introduced the wave function and an intrinsically delocalized picture of the electronic structure in molecules and solids. Thus, the need to regain descriptions of localized phenomena like chemical bonding was immediately evident. Subsequently, a wide range of methods for analysing the electronic structure of compounds was developed. These can be broadly classified as follows: (i) electron-density- (real space) based, (ii) orbital- (Hilbert space) based methods, (iii) methods that model experimentally observable quantities, and (iv) energy-based methods. For molecular quantum chemistry, all these methods have been thoroughly explored in the past.³ For extended structures, the available approaches and experiences are more limited⁴ and shall be briefly summarized.

Direct space methods including the Quantum Theory of Atoms in Molecules and Crystals (QTAIMAC),⁵ the electron localization function (ELF),⁶ and similar approaches have been utilized to solve many questions of solid-state chemistry, also owing to the fact that they can be applied to experimentally derived electron densities as well.⁷ These methods have the advantage of a well-defined quantum mechanical framework but they do not provide rigorous answers, e.g., about the existence of a chemical bond in a system.⁸ This information can only be gained by interpretation of the results based on chemical experience and knowledge.

Hilbert space methods do not analyze the electron density as a whole, but rely on the interpretation of atomic contributions to molecular properties. To this end, they are usually based on atomic orbitals resulting from the basis set approximation in many quantum chemical methods. The most heavily used Hilbert space methods are the approaches to determine partial charges of atoms in molecules and solids⁹ despite their known shortcomings.¹⁰ They are commonly used to determine bond polarity and ionic/covalent character of bonding although charge is a scalar quantity which does not carry information about the charge distribution; this can significantly alter interpretation.¹¹ 

A popular method in surface science is the analysis of density of states (DOS) and partial DOS (pDOS) as more easily accessible representation of the band structure.¹² Closely related, an analysis of bonding and antibonding character of orbitals in an energy interval of a DOS/pDOS can be carried out with the crystal orbital overlap population (COOP)¹³ and the related crystal orbital Hamilton population (COHP) method.¹⁴ Other important methods are Wannier-Type Atomic Orbitals (WTAO),¹⁰ the d-band model by Hammer and Nørskov,¹⁵ reactivity indices¹⁶ and concepts,¹⁷ or computation of spectroscopic properties.¹⁸ Furthermore, natural bond orbital (NBO) analysis was recently extended to periodic systems by Dunnington and Schmidt¹⁹ and Galeev et al.²⁰ Bond order methods date back to Pauling²¹ and the “mobile order” defined by Coulson.²² These methods were subsequently reformulated many times.²³ 

None of the methods described above addresses the question of energetic contributions to chemical bonds, although energy changes are the ultimate driving forces in bond formation. The quantitative description of these contributions is far less developed. One exception is the energy density analysis by Nakai et al., which relies on partitioning of the total energy into atomic contributions and has recently been extended to periodic systems.²⁴ 

For molecular calculations, the energy-based analysis methods are well established and have been used successfully in the past.²⁵ The symmetry-adapted perturbation theory (SAPT)²⁶ method is the most popular method following a perturbational approach. Others include the method by Hayes and Stone.²⁷ More prominent are the variational methods employing a block-partitioning of the Fock matrix. A non-exhaustive list of approaches in this category is the Constrained Space Orbital Variation (CSOV) scheme developed by Bagus et al.,²⁸ the natural energy decomposition analysis (NEDA) by Glendening and Streitwieser,²⁹ the reduced variational self-consistent field (RVS-SCF) approach by Stevens and Fink,³⁰ the AIM-based decomposition scheme by Francisco et al.,³¹ a steric analysis by Liu,³² and the schemes by Mayer³³ and Korchowiec.³⁴ A slightly different method employing absolutely localized molecular orbitals (ALMO) is the recently developed ALMO-EDA by Head-Gordon and co-workers³⁵ which can also be applied to correlated wave functions.³⁶ The block-localized wave function (BLW) method (BLW-ED) based on valence-bond theory has been put forward by Mo.³⁷ 

The most heavily used variational method is the EDA³⁸ based on developments by Morokuma³⁹ and the related extended transition state (ETS) method by Ziegler and Rauk⁴⁰ (jointly called EDA in the following). Variations of this method are found in the recent literature.⁴¹ A valuable extension is the recently developed EDA-NOCV (Natural Orbitals for Chemical Valence) method, which offers additional insight by providing charge and energy analysis in a combined fashion.⁴² Charge transfer and charge redistribution have even been used in the past to define the formation of a surface chemical bond upon adsorption but without the ability for quantitative analysis.^4(c) The strength of these methods is the ability to quantify different bonding contributions in terms of energy. Only the CSOV and the EDA were extended to periodic systems but were only used to a limited extent in cluster-based⁴³ and pilot studies.⁴⁴ Thus, until today, quantitative energy-based analysis of chemical bonding in periodic systems has not been documented for surface-adsorbate interactions in a broader study.

We now introduce the periodic Energy Decomposition Analysis (pEDA) method which decomposes the interaction energy in periodic systems into chemically intuitive terms. Theory and implementation of the method are outlined, the validity and accuracy are tested, and the method is applied to prototype surface science questions spanning different bonding scenarios.

First, we outline the working principle of the molecular EDA method as developed by Ziegler/Rauk⁴⁰ and Morokuma.³⁹ The approach used in the EDA is the investigation of the intrinsic bond energy for the interaction of two fragments A and B forming a molecule AB by separating the bond formation process into several sub-steps. The bond dissociation energy ΔE_bond (which is the negative of the dissociation energy without zero-point vibrational corrections D_e) is thereby composed of the energy changes from a promotion step and an interaction step (Scheme 1(a)) with the respective energy terms (ΔE_prep + ΔE_int),

In the initial step, the ground state (GS) configuration of the fragments A^GS and B^GS is distorted to the particular reference states A and B they have in the molecule AB. This includes geometric distortion and electronic excitation and requires a preparation energy ΔE_prep which is given by the energies of relaxed and distorted fragments,

The promoted fragments can be represented by Slater determinants Ψ_A and Ψ_B built from fragment orbitals at A and B, respectively. They are now interacting to form AB. In this step, the intrinsic bond energy ΔE_int results from the energy difference between the molecule (E_AB) and the prepared fragments (E_A, E_B). The basic idea of the EDA scheme now is the partitioning of the interaction energy ΔE_int into well-defined terms (Scheme 1(a)),

First, the distorted fragments are brought from infinity into the position they occupy in the molecule AB without optimisation of the resulting wave function. The associated energy change for this step is the quasiclassical electrostatic interaction (ΔE_elstat) between the two charge distributions ρ_A and ρ_B of both fragments. This term is usually attractive and gives a quantitative estimate of the electrostatic bonding contributions, which are neglected in a purely orbital-interaction based analysis.⁴⁵ The energy after this step corresponds to a simple product wave function {Ψ_AΨ_B} which does not fulfil the Pauli exclusion principle. In the next step, this product wave function is antisymmetrized (⁠$A ̂$⁠) and normalized (N). This leads to the intermediate wave function Ψ⁰ with the energy E⁰,

The orbitals making up Ψ⁰ are obtained via separate Löwdin orthogonalization of occupied and virtual fragment orbitals with Gram-Schmidt orthogonalization of occupied on virtual orbitals in between. While the Löwdin scheme leads to orthogonal functions, which resemble the non-orthogonal set, the Gram-Schmidt step is more suitable to orthogonalize a set of functions with respect to a set of other functions which shall not be changed.⁴⁶ The corresponding term ΔE_Pauli is destabilizing since the orthogonalization procedure adds additional nodes to the orbitals and thus leads to an increase of the kinetic energy.⁴⁷ In the final step, the frozen molecular orbitals of the intermediate wave function are allowed to relax and the optimal wave function Ψ_AB for the molecule AB with the energy E_AB is found. This energy (ΔE_orb) can be expressed in terms of the density ρ⁰ of the intermediate wave function Ψ⁰ and the density ρ_AB of the final wave function Ψ_AB,

Following the discussion in recent years about the importance of dispersion interactions in density functional theory (DFT)-based periodic and non-periodic computations,⁴⁸ we chose to add a dispersion term to the results of the pEDA scheme. The validity of treating this term separately from the EDA has been shown before⁴⁹ and an improved performance was found.⁵⁰ Therefore, we chose to obtain via the semi-empirical correction scheme DFT-D3 put forward by Grimme et al.⁵⁷ via the difference of fragment and complex dispersion energies,

The alternative approach to use dispersion-corrected density functionals (e.g., VV10⁵¹) would disable the separate discussion of this bonding contribution although dispersion would then be treated at the density level.

For the simulation of extended systems, we chose an approach based on atom centered, atomic orbital (AO) basis functions, {ϕ_μ}, to generate the corresponding Bloch functions, {Φ_μ(k)}. Linear combination of these basis functions leads to a set of crystal orbitals (CO) {φ_i(k)} for every point in reciprocal space k, where the coefficients c_i,μ(k) are the transformation matrix elements of T(k),

As a consequence of the Bloch theorem, our wave functions are now dependent on the reciprocal wave vecto k. The first Brillouin zone in reciprocal space contains an infinite number of k-points but it is enough to evaluate only a sub-set of k-points to derive at a good approximation of the wave function and the density of the extended system (k-space sampling).

For every k-point, the Kohn-Sham equations have to be solved to derive Ψ_A(k) and Ψ_B(k) which are described in the basis of the Bloch functions by matrices T_A(k) and T_B(k). The total charge density of the fragments, ρ_A and ρ_B, is then the weighted sum of all densities at the k-points sampled,

For the calculation of the energy terms, ΔE_Pauli and ΔE_orb, the definition of the orthogonalized wave function, Ψ⁰(k), is necessary. Based on work by Philipsen and Baerends (PB),^44(a) the algorithm to construct the intermediate wave function can be separated in three steps. The initial step is the combination of the separated eigenstates, Ψ_A(k) and Ψ_B(k), to form the normalized eigenstate N{Ψ_A(k)Ψ_B(k)}, which is described by expansion of the Bloch basis {Φ_μ(k)} with the transformation matrix T₁(k) for every k-point. This transformation matrix can be further separated in two parts which characterize the occupied COs, denoted as T_o,1(k), and the virtual COs, denoted as T_v,1(k). Now, the occupied COs of A and B are orthogonalized via the Löwdin formalism. Only, S_o/o,1(k), the overlap matrix for the occupied COs of the initial eigenstate is needed to calculate the transformation matrix, T_o,2(k), for the orthogonalized, occupied fragment COs,

Then, the virtual COs of A and B are orthogonalized onto the newly formed occupied fragment COs via the Gram-Schmidt formalism. Here, the overlap matrix for the orthogonalized occupied COs, S_o/o,2(k), is the identity matrix and the correction matrix T^GS(k) is described solely by the overlap matrix between occupied and virtual crystal orbitals S_o/v,2(k),

with

The last step is the Löwdin orthogonalization of the newly formed virtual CO set T_v,2(k). Here, S_v/v,2(k) is the overlap matrix for these virtual COs,

So, the occupied and the virtual space of the transformation matrix T⁰(k) is constructed by T_o,2(k) and T_v,3(k) and expands the Bloch basis {Φ_μ(k)} in terms of the orthogonalized fragment COs {λ_i(k)},

Now, the intermediate wave function, Ψ⁰(k), is described by orthogonalized fragment COs and forms the intermediate crystal density ρ⁰, which can be assigned the energy value E⁰.

The energy difference between the separated fragments and this intermediate eigenstate is calculated as follows:

The energy difference ΔE⁰ can be expressed as the differences of the kinetic energy T of the electrons, the potential energy (Coulomb energy V_Coul), and the exchange-correlation energy (V_XC) between the two fragment states A and B and the orthogonalized, intermediate state Ψ⁰,

Since the resulting wave function Ψ⁰(k) incorporates the quasiclassical, electrostatic interaction and the Pauli exclusion principle, the term $Δ V Coul 0$ comprises the Coulomb interaction due to the overlapping fragment densities, ΔE_elstat, and the change of the Coulomb interaction due to the orthogonalization of the fragment wave functions, ΔV_P,Coul,

Now, the energy term, ΔE_Pauli, corresponding to the Pauli exclusion principle, is calculated only by terms arising due to the orthogonalization scheme,

In the last step, the intermediate wave function is allowed to relax to the final wave function Ψ_AB with the electron density ρ_AB. The corresponding energy change is called the orbital relaxation energy term

The overall interaction and bond dissociation energies can now be calculated according to Eqs. (1) and (3) (see Scheme 1(b)). If the fragments need to be chosen representing an electronic state which is not the ground state, the occupation numbers of fragment orbitals can be specified for the purpose of analysis.⁴⁷ At the current state, this is possible in the pEDA if k-space sampling is restricted to the Γ-point due to the ambiguity of setting consistent occupation numbers for bands, which might change their energetic ordering at different points in k-space. The same restriction currently applies to the decomposition into spin polarized fragments.

Unconstrained structural optimizations for the molecular systems investigated here were carried out in the framework of DFT by employing the BP86⁵² exchange-correlation functional together with the triple zeta (TZ2P) basis set⁵³ applying the frozen core approximation, an accuracy setting of 5, and incorporating relativistic effects within the zeroth order regular approximation (ZORA)⁵⁴ without periodic boundary conditions (PBC). All calculations utilizing the frozen core approximation employed a “small” frozen core and are indexed by a suffix (fc) at the basis set level (e.g., TZ2P(fc)). These structures were analyzed with the EDA method as outlined above. Other basis sets (double zeta (DZ), triple zeta (TZP), and quadruple zeta (QZ4P))⁵³ were used for convergence studies. All non-periodic calculations (optimizations and EDA) were carried out with the Amsterdam Density Functional (ADF) molecular modeling suite (version 2012).⁵⁵ Subsequently, the resulting structures were re-optimized under one-dimensional PBC with a unit cell length of 50 Å. Based on these structures, pEDA calculations were performed with k-space sampling restricted to the Γ-point.

The surface-adsorbate complexes investigated here were calculated applying two-dimensional PBC in the surface plane. For the adsorption of H₂ on Cu(001) and Pd(001) surfaces, the metal was approximated with two layers of atoms and a (2 × 3) super cell with lattice constants a(Cu) = 3.61 Å and a(Pd) = 3.95 Å. The H₂ molecule was placed parallel to the surface in bridging position. This corresponds to the setup used in the previous EDA study on this system.^44(a) 

For the adsorption of CO on the MgO(001) surface, the insulator surface was approximated with three layers of Mg atoms and the corresponding O atom layers employing a lattice parameter of 4.25 Å. Surface rippling was approximated by displacing the positions of the Mg/O atoms inward/outward by 0.001 Å, respectively. The CO molecule was placed perpendicular to the surface above a Mg atom with a Mg—C distance of 2.61 Å. Three different super cells were chosen, resulting in coverages of θ = 0.5, 0.25, and 0.125. This corresponds to the setup used in the previous CSOV study on this system.^44(b) 

The optimal adsorption structures of CO and C₂H₂ molecules on the reconstructed Si(001)c(4 × 2) surface were found by carrying out spin-polarized, constrained structural optimizations (applying three-dimensional PBC) with the two bottom layers of the six-layer silicon slab kept fixed. These calculations were performed using the GGA functional developed by Perdew, Burke and Ernzhofer (PBE)⁵⁶ considering dispersion effects via the DFT-D3 scheme with an improved damping function⁵⁷ (in the following PBE-D3). The projector-augmented wave (PAW) method was employed allowing for a kinetic energy cutoff of 350 eV for the plane wave (PW) basis set.⁵⁸ The Brillouin zone was sampled by a 4 × 2 × 1 Γ-centered Monkhorst-Pack k-mesh.⁵⁹ Based on the resulting structures, a Si₁₅H₁₆ cluster was derived, approximating the infinite surface in a finite model. While the Si positions were kept frozen, the dangling bonds were saturated by H atoms. These capping atoms were placed along the broken Si—Si bond and set to a typical Si-H bond length of 1.480 Å.⁶⁰ 

Bonding analysis for these extended structures was carried out via pEDA calculations. If not otherwise noted, the pEDA calculations were done by employing the BP86 or PBE exchange-correlation functional with a TZ2P(fc) basis set, an accuracy setting of 5, applying ZORA, and taking into account the Γ-point only.

Structural optimizations and PBE-D3 calculations of extended systems were carried out with the VASP code (version 5.2.12).⁶¹ Bond analysis calculations employing PBC (pEDA) were carried out with a development version of the BAND code.⁶² 

The pEDA was tested for a wide range of bonding scenarios by choosing three molecular and four extended test systems. For the molecular systems, three different prototypical bonding types were investigated: (i) main group donor-acceptor bonding (H₃N —BH₃, Scheme 2(a)), (ii) transition metal donor-acceptor bonding (CO —Cr(CO)₅, Scheme 2(b)), and (iii) shared electron bonding (H₃C —CH₃, Scheme 2(c)). The fragmentation applied in the EDA and pEDA decompositions is indicated in Scheme 2. The inclusion of non-periodic systems in this convergence study enables us to validate the numerical results against the established EDA procedure and quantify the outcome of using spin polarization for unrestricted fragments (Scheme 3).

Nevertheless, the main aim of the new method is thus to investigate the chemical bonding in periodic systems. Thus, we chose four different bonding scenarios: (i) adsorption of H₂ on the Cu(001) and Pd(001) surfaces (Figure 1), (ii) adsorption of CO on the MgO(001) surface (Figure 2), and (iii) adsorption of CO (Figure 3) and C₂H₂ (Figure 4) on the Si(001)c(4 × 2) surface exhibiting the well-known buckled dimer reconstruction. On the one hand, these examples cover the regimes of conducting (i), insulating (ii), and semiconducting (iii) surfaces. On the other hand, the two examples for the semiconducting surface are typical examples for donor-acceptor (CO) and shared-electron (C₂H₂) bonding to the surface as we will see later on (Scheme 4). The metallic system was chosen since it enables us to compare the results to the previous implementation of a periodic EDA by Philipsen and Baerends.^44(a) The insulating system was investigated by periodic CSOV analysis before and thus lends the opportunity to compare these two methods.^44(b) 

A robust analysis method needs to provide reliable results with respect to the computational parameters in the calculation. Therefore, we checked the convergence behavior of the pEDA results for the following settings: (i) basis set and frozen core usage, (ii) accuracy setting, and (iii) k-space sampling.

We tested basis sets of DZ, TZP and TZ2P, and QZ4P quality with and without frozen core approximation (large and small frozen core). The general precision parameter (accuracy) determines the generation of integration points and many other parameters related to the accuracy of the results.⁶³ The parameter for the k-space sampling can be set to include the Γ-point only (setting 1), use a linear tetrahedron method (setting 2, 4, …) or a quadratic tetrahedron method (setting 3, 5, …).⁶⁴ The influence of these parameters will be checked in Secs. IV B 1 and IV B 2 for molecular and extended systems.

For the three molecular test systems, we derive an accuracy parameter setting of 5 for well-converged results. The TZ2P(fc) basis set is a good trade-off between a too inflexible basis (DZ) and a large basis with linear dependency issues (QZ4P). The frozen core approximation (small fc) leads to small errors but considerable savings in computer time. Further discussion is found in the supplementary material.⁸³ 

In a further step, we compared the results from pEDA analysis with the established EDA analysis for the setup derived above (TZ2P(fc) with accuracy 5, see Table I).

The donor-acceptor bonds are discussed first. In both cases, the EDA and pEDA results agree very well with an absolute deviation of <1% for the energy terms in the upper part of the table. The intermediate terms, ΔT⁰ and ΔV⁰, exhibit larger differences. This is probably due to the differences in the basis set definitions between the non-periodic and periodic programs used, whereby the latter code neglects some diffuse functions to improve the numerical stability of the program.⁵³ 

The deviations are more significant for the shared-electron bond in ethane. This can be traced back to the different descriptions of the fragments in both approaches. The EDA was extended to unrestricted fragments previously, enabling a successful application of the method to shared-electron bonding in many molecular systems like the CN-dimer^65(a) or an explanation for the rotational barrier in ethane.^65(b) Nevertheless, the original EDA implementation neglects spin polarization in the fragments leading to a “spin-flip” error (ΔE_spin-flip) for the bond energies.⁶⁶ Although this error is often small, it amounts to −121.3 kJ mol⁻¹ for this fragmentation. The pEDA includes this spin polarization of the fragments. The improper description of the fragment density has an indirect effect on all analysis terms, leading to higher values for ΔE_Pauli, ΔE_elstat and ΔE_orbital. Note that the bonding description in ethane does not change qualitatively compared to a previous, detailed study of this prototype molecule.^65(b) Thus, it becomes clear that the usage of spin-polarized fragments in the pEDA leads to a refined bonding description for shared-electron bonding.⁶⁷ 

For the extended test systems, the results for the convergence with respect to basis set size and frozen core approximation are found in Table II. The convergence with respect to k-space sampling is summarized in Figure 5 for the examples H₂ on Cu (Fig. 5, left) and C₂H₂ on Si (Fig. 5, right). The underlying numbers and results for the other systems can be found in the supplementary material.⁸³ 

For the basis set study, the findings are comparable to the molecular systems. The frozen core approximation leads to rather small changes in the pEDA terms with few exceptions: the DZ(fc) result for CO on MgO(001) as well as the QZ4P(fc) results for the same system and the C₂H₂ on Si(001) case differ significantly from the all-electron calculations. For the most relevant TZP and TZ2P basis sets, the saving in computing time (approximately 20%) is well worth the slight errors introduced which are below chemical accuracy (4.2 kJ mol⁻¹) for the main pEDA terms.

Regarding the basis set, the TZP and TZ2P results are again very similar for the systems investigated. Surprisingly, even the rather small DZ basis set delivers a very good qualitative picture with the exception of the DZ(fc) results discussed above. While the absolute numbers can differ by up to 15%, most results are much closer to TZP. Care has to be taken when trends are discussed with the DZ basis set, since description of the relative contributions of ΔE_elstat and ΔE_orb can differ in comparison to the larger basis sets. The results for QZ4P are again noteworthy. Initially, it seems that the triple zeta results are not converged yet with respect to the larger basis set. A closer look reveals that the large deviations observed in going from TZ2P to QZ4P (up to 44 kJ mol⁻¹ as in the ΔE_orb term for C₂H₂ on Si(001)) can be explained again by linear dependencies which show up in the QZ4P calculations and which are indicated in Table II. As another check for the convergence of the TZ2P results, the ΔE_int values were compared to results from a calculation based on a plane wave basis set shown in the most-right column of Table II. This agreement points toward a balanced description of fragments and complex and thus a control of the basis set superposition error (BSSE).⁶⁸ It becomes clear that for the metallic and the semiconducting system, the TZ2P basis set delivers well-converged dissociation energies.

Next, we want to discuss the convergence with respect to the sampling of reciprocal space as shown in Figure 5. For the metallic system (Fig. 5, left), a tight sampling of more than 15 k-points per Å⁻¹ is needed to converge the pEDA terms <1%. It is a general feature of the computation of metallic systems with PBC that a tight k-mesh is needed and the influence on the EDA terms was noted before.^44(a) Specific for the method here is the different convergence behaviors for the pEDA terms. While ΔE_elstat converges very quickly (<5 k-points per Å⁻¹), ΔE_orb and ΔE_Pauli are much more strongly dependent on a good sampling of the Brillouin zone. This has to be taken into account when bonding trends (relative contributions of these terms) are discussed for a limited k-point sampling. A much quicker convergence is observed for the semiconducting example C₂H₂ on Si(001) (Fig. 5, right) where even sampling at the Γ-point is a very good approximation.⁶⁹ For the insulating system, CO on MgO(001), it is sufficient to include the Γ-point for converged numbers (except for the high coverage regime, see Table S6 in the supplementary material⁸³). The accuracy parameter convergence is much less system dependent (Figure S1⁸³) and is converged to <1% for an accuracy setting of 4 or 5 (except for the small energy terms arising for CO on MgO(001)).

The convergence tests for several extended systems thus lead to the conclusion that a triple zeta basis set with frozen core approximation together with an accuracy setting of 4-5 should be used and the k-space sampling has to be adjusted to the electronic structure of the surface investigated. This standard will be employed in Subsection IV C to discuss the chemical bonding in the test systems.

The dissociative adsorption of H₂ is a non-activated process on Pd, while an activated process is observed on Cu. This leads to the aforementioned investigation by PB focusing on the role of Pauli repulsion in this adsorption process.^44(a) It was concluded that the relief of Pauli repulsion is the decisive factor for the different behaviors on Pd and Cu. The difference in our implementation is now that we are able to explicitly calculate the ΔE_elstat and ΔE_Pauli terms, while the ΔE⁰(=ΔE_elstat + ΔE_Pauli) and ΔE_orb terms are directly comparable to the analysis by PB. Are our results in agreement with PB and do we gain additional insight by the further decomposition of the steric interaction term?

A series of pEDA analyses was conducted for H₂ approaching the surfaces (Figure 1). The resulting energy terms along the reaction coordinate (d(M—H₂)) are shown in Figure 6 and the respective numbers can be found in the supplementary material.⁸³ The same data in the format used by PB can be found in the supplementary material for direct comparison (Figure S2).⁸³ 

We confirm the finding by PB that the interaction energy is lower, the Pauli repulsion higher, and the orbital interaction weaker on Cu compared to Pd along the reaction coordinate (Figure 6(a)). Additionally, we can now state that the electrostatic stabilization term ΔE_elstat is also more favorable on the Pd surface (Figure 6(b)). All differences get larger along the reaction coordinate. Thus, we obtain an additional stabilization term on the Pd surface of electrostatic nature which adds to the picture of a non-activated process on this surface for the dissociative adsorption of H₂. Nevertheless, the Pauli repulsion stays the leading term as found by PB and thus dominates the pEDA differences between the two surfaces especially at larger distances. The different trends for the surfaces can be understood in terms of electron density at the surface as indicated via the local density of states (LDOS) shown in Figure 7.

The electron density at the Cu surface (Figure 7(a)) is thereby found to be significantly reduced in comparison to the Pd surface (dark blue areas in Figure 7(b)). Thus, the approaching H₂ molecule interacts with more surface electrons in the latter case which is reflected in higher values for the attractive terms, ΔE_orb and ΔE_elstat, as well as a higher value for ΔE_Pauli on Pd along the reaction coordinate (Figure 6).

Carbon monoxide is a typical probe molecule for surface science, while the MgO(001) surface is a prototype for an insulating surface. Therefore, this system has been widely investigated in the past and among other methods, it has been analyzed by the CSOV method in a periodic implementation by Hernandez, Zicovich-Wilson, and Sanz (HWS).^44(b) This analysis confirmed previous findings⁷⁰ indicating electrostatic forces and Pauli repulsion as the main driving forces in this weakly bound system. Recently, the adsorption energy in this system was derived by highly accurate ab initio methods and the extrapolated adsorption energy found 21.0 ± 1.0 kJ mol⁻¹ is in very good agreement with the most reliable experimental values.⁷¹ 

To test our pEDA implementation, we reconsidered the setup by HWS and analyzed CO adsorption on MgO for three different coverages (see Figure 2). The pEDA results are found in Table III. First of all, the adsorption energy derived is in very good agreement with the above-mentioned high level value of 21.0 kJ mol⁻¹ (ΔE_int = − 22.2 kJ mol⁻¹ for θ = 0.125). There is certainly a degree of error cancellation for our computational setup (possible error sources being, e.g., relaxation, zero-point vibrational energy, thermal corrections, and slab setup) but it lends confidence to the pEDA results besides the rather small energy terms. As seen in Subsection IV C 1, care has to be taken to use well converged computational parameters for such a weakly bound system.

The bonding picture given by the pEDA shows domination by a strong Pauli repulsion and electrostatic attraction, while the ΔE_orb term is only approximately a third of the attractive interactions. This agrees well with the CSOV results by HWS which can be quantified by summing up selected pEDA terms to deliver their CSOV counterparts: the sum of ΔE_Pauli and ΔE_elstat can be compared to the frozen orbital (FO) term, while the ΔE_orb term corresponds to a sum of the polarisation (Pol) and charge transfer (CT) terms. We see the same trend in the orbital term in both methods—with an absolute difference of approximately 5 kJ mol⁻¹—but a qualitative difference in the FO term. While the CSOV gives a rather constant value of approximately 8 kJ mol⁻¹ across the three coverage regimes, the pEDA leads to a decrease of the respective sum in going from θ = 0.5 (17.5 kJ mol⁻¹) to θ = 0.125 (4.2 kJ mol⁻¹). In consequence, the bonding energy (BE) in the CSOV analysis is constant across the three coverages investigated, while the pEDA leads to an increase of the interaction energy ΔE_int. While a quantitative agreement should not be expected due to the different density functionals used and the neglected contributions outlined above, it is still noteworthy that the trend in our data is consistent with the experimentally observed linear increase of the adsorption energy towards lower coverage.⁷² Thus, we conclude that the bonding in CO on MgO(001) is indeed dominated by Pauli repulsion and electrostatics which nearly cancel leading to a weak bonding as pointed out before,^44(b) but the relative importance of these terms shifts upon decreasing coverages and results in a higher bonding energy in the low coverage regime.

Carbon monoxide can also be taken as a probe molecule for semiconductor surfaces. The most investigated and technologically most relevant semiconductor surface is Si(001). This surface exhibits a well-known c(4 × 2) reconstruction revealing buckled dimers with one Si atom lying above the dimer plane (Si_up) and one Si atom below the plane (Si_down).⁷³ Analysis of the band structure reveals one occupied and one unoccupied surface states.⁷⁴ In a localized view on the bonding, this can be understood as the Si_up atom exhibiting an electron lone-pair and the Si_down atom an empty orbital. Therefore, this surface dimer is prone to electrophilic (Si_up) or nucleophilic (Si_down) attack.⁷⁵ Thus, the interaction of the CO molecule with the Si_down atom can be used as a typical example for a donor acceptor bond on a semiconductor surface (Scheme 4(a)).

This bonding mode was confirmed by experimental and computational studies.⁷⁶ We therefore do not aim at reinvestigating the adsorption behavior of this system. From analysis of vibrational signatures, a covalent bonding character with typical donation/back-donation characteristics of CO was concluded.^76(c) This is in agreement with the Dewar-Chatt-Duncanson model for transition metal complexes⁷⁷ which is known as the Blyholder model in surface science.⁷⁸ How can we quantify this view on the bonding of CO on Si(001)?

The structural optimization of one CO molecule on the Si(001)c(4 × 2) unit cell (resulting in a coverage of θ = 0.25 with respect to surface Si atoms) as outlined in Sec. III leads to a structure which is in quantitative agreement with previous PBE results regarding the bond lengths (Figure 3(a)).^76(c) The results for the bonding analysis on this system are shown in Table IV. The bond dissociation energy (ΔE_bond = − 102.3 kJ mol⁻¹) is higher compared to the previous finding (E_ads = − ΔE_bond = 0.92 eV = 88.8 kJ mol⁻¹).^76(c) This difference can be understood by the inclusion of dispersion effects in our study (ΔE_disp = 10.7 kJ mol⁻¹) and the dissimilarities in the computational setup (e.g., slab size and basis set). The pEDA results now give a rather high Pauli repulsion term (ΔE_Pauli = 855.1 kJ mol⁻¹) and near equal contribution of attractive electrostatic (46.0%) and orbital (54.0%) terms. This ratio was also found in a recent EDA study on CO bonding to a molecular Lewis acid.⁷⁹ According to the models discussed above, the main contribution to the orbital term stems from the interaction of the highest occupied molecular orbital (HOMO) of the CO molecule with the unoccupied surface state which is mainly localized at the Si_down atom (Figure 3(a)). The strong overlap between the two orbitals explains the high value for the orbital interaction term.

The rather localized nature of this interaction enables us to analyze this bond additionally in a cluster model. A Si₁₅H₁₆ cluster cut out from the slab setup bearing one CO molecule is shown in Figure 3(b). On the one hand, we can now check the numerical accuracy of the new implementation by comparing pEDA data on the cluster with EDA results. On the other hand, we can discuss the influence of including the full environment of the surface in the periodic slab model against the truncated nature of the cluster model.

The EDA and pEDA data for the cluster model are in nearly quantitative agreement (<1%) confirming the accuracy of the new method. The comparison of pEDA results for cluster and slab approach shows that the bond of CO to Si(001) can be described rather well with a small cluster model. Nevertheless, a closer look reveals that the dispersion energy term (ΔE_disp) and the preparation energy of the surface (ΔE_prep (Si)) differ more strongly than the other terms since an extended system is replaced by a zero-dimensional cluster. Another difference can be seen in the orbital/band energies shown in Figure 3. While the HOMO energy for the CO does not depend on the application of PBC, the energy for the lowest unoccupied molecular orbital (LUMO) of the cluster differs significantly from the lowest unoccupied crystal orbital (LUCO) of the surface. This is a relevant difference if one is interested in effects like energy level matching for solar cell optimization.⁸⁰ 

As a last example, we analyzed the bonding of acetylene to a Si(001)c(4 × 2) substrate as an example for shared-electron bonding between surface and adsorbate (Scheme 4(b)). This enabled us to test the capabilities of the pEDA to include spin polarization in unrestricted fragments. This system has been much investigated before experimentally⁸¹ as well as computationally⁸² and is a prototype system for surface functionalization of semiconductors. Similar to the previous system, we compared the pEDA results for the periodic system (Figure 4(a)) to EDA results on a cluster model (Si₁₅H₁₆, Figure 4(b)) and pEDA results for the same cluster (Table V).

The adsorption structure shown in Figure 4(a) reveals two unequal Si—C bond lengths (1.896 and 1.913 Å) which are due to the interactions with the surrounding buckled dimers. The C—C bond of acetylene increases upon adsorption from 1.205 Å to 1.358 Å. These bonding parameters are in good agreement with typical values for a C=C double and two Si—C single bonds in line with the Lewis picture sketched in Scheme 4(b).

The bond dissociation energy (ΔE_bond) for the three approaches chosen is very similar and in good agreement with the literature values.⁸² In contrast, the bonding analysis terms show notable differences. As was seen for the molecular case of ethane above, the pEDA exhibits numerically smaller values for all analysis terms compared to the EDA. As in the case of ethane discussed above, this can be traced back to the inclusion of fragment spin polarization for the unrestricted fragments in the new method. The differences are significant with deviations of up to 132.5 kJ mol⁻¹ (ΔE_Pauli). The largest relative deviations between the two methods are as expected seen in the preparation energy (ΔE_prep) of both fragments where the spin-flip error has to be taken into account for the EDA. A recalculation of the cluster approach with the pEDA delivers very similar terms and reveals the intrinsic differences between cluster and slab approach: higher ΔE_elstat and ΔE_prep (Si) terms in the finite model. The relative magnitude of the pEDA terms points toward a higher orbital contribution compared to CO on this substrate (Table IV) and larger terms in all cases due to two bonds being formed at the same time. We see that the preparation energy of the acetylene fragment is very significant; thus, the linear adsorbate requires a strong bending and an additional spin-flip to enable bonding to the Si surface.

The frontier orbitals for the fragments after the preparation step are shown in Figure 4. For finite and infinite approaches, the singly occupied orbitals are set up for two covalent σ-bonds between adsorbate and surface atoms. As seen for the previous system, the orbital energies (and also the energy differences) differ for both approaches which might be relevant for charge-transfer investigations.

We present an implementation of the energy decomposition analysis method for periodic systems (pEDA) on the density functional level. This enables the analysis of energy contributions for surface-adsorbate bonding in diverse systems. The implementation enables the computation of all energy terms (specifically, separate electrostatic and Pauli contributions compared to previous work by Philipsen and Baerends) and thus leads to quantification of preparation, electrostatic, electron-repulsion, and orbital contributions to the chemical bond in periodic systems together with a consideration of dispersion energy. We tested the pEDA against the established non-periodic EDA method for molecular systems with prototypical bonding scenarios (donor-acceptor bonding for main group and transition metals and shared-electron bonding). Furthermore, we tested the pEDA for four surface-adsorbate complexes representing insulating (CO on MgO(001)), metallic (H₂ on M(001), M = Pd, Cu), and semiconducting (CO and C₂H₂ on Si(001)c(4 × 2)) substrates. We confirmed the numerical stability of the pEDA and derived a reliable set of computational parameters (TZ2P(fc) basis set, accuracy 5, k-sampling at the Γ-point for semiconductors). The comparison to previous results showed new valuable insight being gained by the pEDA and the analysis of shared-electron bonding reveals the necessity for including spin polarisation of fragments in an unrestricted Kohn-Sham formalism. The need for a periodic description of the extended systems in contrast to a cluster approach for detailed bonding analysis was shown. The results presented lend confidence that the pEDA will be a powerful tool for the analysis of chemical bonding in surface and materials science in the future.

Funding by the Deutsche Forschungsgemeinschaft (DFG) through Project No. TO 715/1-1 and the Research Training Group 1782 is gratefully acknowledged. We thank SCM/Amsterdam for kindly providing a developer’s license for the BAND code and support. We are indebted to the late Professor Tom Ziegler (Canada) for helpful discussions and support. Computational resources were provided by HRZ Marburg, LOEWE-CSC Frankfurt, and HLR Stuttgart. We are very grateful to Professor Gernot Frenking (Marburg) for sparking our interest in the world of chemical bonding unicorns with his endless enthusiasm for the topic.