A periodic Energy Decomposition Analysis (pEDA) method for the Investigation of Chemical Bonding in Extended Systems

The development and first applications of a new periodic energy decomposition analysis (pEDA) scheme for extended systems based on the Kohn-Sham approach to density functional theory are described. The pEDA decomposes the binding energy between two fragments (e.g. the adsorption energy of a molecule on a surface) into several well-defined terms: preparation, electrostatic and dispersion interaction, Pauli repulsion and orbital relaxation energies. The pEDA presented here for an AO-based implementation can handle restricted and unrestricted fragments for 0D to 3D systems considering periodic boundary conditions with and without the determination of fragment occupations. For the latter case, reciprocal space sampling is enabled. The new method gives comparable results to established schemes for molecular systems and shows good convergence with respect to the basis set (TZ2P), the integration accuracy and k-space sampling. Four typical bonding scenarios for surface adsorbate complexes were chosen to highlight the performance of the method representing insulating (CO on MgO(001)), metallic (H$_2$ on M(001), M = Pd, Cu) and semiconducting (CO and C$_2$H$_2$ on Si(001)c(4x2)) substrates. These examples cover the regimes of metallic, semiconducting and insulating substrates as well as bonding scenarios ranging from weakly interacting to covalent (shared electron and donor acceptor) bonding. The results presented lend confidence, that the pEDA will be a powerful tool for the analysis of surface-adsorbate binding in the future, enabling the transfer of concepts like ionic and covalent binding, donor-acceptor interaction, steric repulsion and others to extended systems.


Introduction
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, 1 the efforts to extend the concepts toward extended systems are scarcer. Still, quantitative analysis of chemical bonding in extended systems leads to greater fundamental understanding of the systems investigated and can lead to predicting new trends and reactivity for ab initio material's design and functionalization of surfaces and interfaces. 2 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 and (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. 3 For extended structures, the available approaches and experiences are more limited 4 and shall be briefly summarized.
Direct space methods including the Quantum Theory of Atoms in Molecules and Crystals (QTAIMAC) 5 , the electron localization function (ELF) 6 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. 7 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. 8 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 9 despite their known shortcomings. 10 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, which can significantly alter interpretation. 11 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. 12 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) 13 and the related COHP method. 14 Other important methods are Wannier-Type Atomic Orbitals (WTAO) 10 , the d-band model by Hammer and Nørskov 15 , reactivity indices 16 and concepts 17 or computation of spectroscopic properties. 18 Furthermore, natural bond orbital (NBO) analysis was recently extended to periodic systems by Dunnington et al. 19 and Galeev et al. 20 Bond order methods date back to Pauling 21 and the "mobile order" defined by Coulson 22 and were reformulated many times. 23 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. 24 For molecular calculations, the energy-based analysis methods are well established and have been used successfully in the past. 25 The symmetry-adapted perturbation theory (SAPT) 26 method is the most popular method following a perturbational approach. Other include the method by Hayes and Stone. 27 More prominent are the variational methods employing a blockpartitioning of the Fock matrix. A not exhaustive list of approaches in this category is the Constrained Space Orbital Variation (CSOV) scheme developed by Bagus and Illas, 28 the natural energy decomposition analysis (NEDA) by Glendening and Weinhold,29 the reduced variational self-consistent field (RVS-SCF) approach by Stevens and Fink, 30 the AIM-based decomposition scheme by Francisco et al. 31 , a steric analysis by Liu 32 and the schemes by Mayer 33 and Korchoviec. 34 A slightly different method employing absolutely localized molecular orbitals is the recently developed ALMO-EDA by Head-Gordon and co-workers 35 which can also be applied to correlated wave functions. 36 The block-localized wave function method (BLW-ED) based on valence-bond theory has been put forward by Mo. 37 The most heavily used variational method is the Energy Decomposition Analysis (EDA) 38 based on developments by Morokuma 39 and the related extended transition state (ETS) method by Ziegler and Rauk 40 (jointly called EDA in the following). Variations of this method are found in the recent literature. 41 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. 42 Charge transfer and charge redistribution has even been used in the past to define the formation of a surface chemical bond upon adsorption but without the ability for quantitative analysis. 4c 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 43 and pilot studies. 44 Thus until today, quantitative energy-based analysis of chemical bonding in periodic systems has neither 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 is tested and the method is applied to prototype surface science questions spanning different bonding scenarios.

2.1
The EDA method First, we outline the working principle of the molecular EDA method as developed by Ziegler/Rauk 40 and Morokuma 39 . 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 (which is the negative of the dissociation energy without zero-point vibrational corrections De) is thereby composed of the energy changes from a promotion step and an interaction step (Scheme 1a) with the respective energy terms ( ). (1) Scheme 1 In the initial step, the ground state (GS) configuration of the fragments A GS and B GS are distorted to the particular reference state A and B they have in the molecule AB. This includes geometric distortion and electronic excitation and requires a preparation energy which is given by the energies of relaxed and distorted fragments. ( The promoted fragments can be represented by Slater determinants and built from fragment orbitals at A and B respectively. They are now interacting to form AB. In this step, the intrinsic bond energy results from the energy difference between the molecule ( ) and the prepared fragments ( , ). The basic idea of the EDA scheme now is the partitioning of the interaction energy into well-defined terms (Scheme 1a). (3) 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 ( ) between the two charge distributions and 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. 45 The energy after this step corresponds to a simple product wave function which does not fulfil the Pauli Exclusion Principle. In the next step, this product wave function is antisymmetrized ( ) and normalized (  48 and an improved performance was found. 49 Therefore we chose to obtain via the semi-empirical correction scheme DFT-D3 put forward by Grimme et al. 56 via the difference of fragment and complex dispersion energies: The alternative approach to use dispersion-corrected density functionals (e.g. VV10 50 ) would disable the separate discussion of this bonding contribution.
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, . Linear combination of these basis functions leads to a set of crystal orbitals (CO) for every point in reciprocal space k, where the coefficients are the transformation matrix elements of .
As a consequence of the Bloch theorem, our wave functions are now dependent on the reciprocal wave vector 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 and which are described in the basis of the Bloch functions by matrices and .
The total charge density of the fragments, and , is then the weighted sum of all densities at the k-points sampled: These charge densities are now used to calculate the electrostatic energy term, : .
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, , is the identity matrix and the correction matrix is described solely by the overlap matrix between occupied and virtual crystal orbitals The last step is the Löwdin orthogonalization of the newly formed virtual CO set . Here is the overlap matrix for these virtual COs.
The energy difference can be expressed as the differences of the kinetic energy of the electrons, the potential energy (Coulomb energy ) and the exchange-correlation energy ( ) between the two fragment states A and B and the orthogonalized, intermediate state .
Since the resulting wave functions incorporates the quasiclassical, electrostatic interaction and the Pauli exclusion principle, the term comprises the Coulomb interaction due to the overlapping fragment densities, , and the change of the Coulomb interaction due to the orthogonalization of the fragment wave functions, .
Now the energy term, , 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 with the electron density . The corresponding energy change is called the orbital relaxation energy term: The overall interaction and bond dissociation energies can now be calculated according to eq.
(1) and eq. 3 (see Scheme 1b). Occupation numbers can be set for the fragments chosen if kspace sampling is restricted to the -point. The same restriction currently applies to the decomposition into spin polarized fragemtns. were performed with k-space sampling restricted to the Γ-point.

Extended systems
The surface-adsorbate complexes investigated here were calculated applying two-dimensional The optimal adsorption structures of CO and C2H2 molecules on the reconstructed Si(001)c(4x2) 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 PBE 55 functional considering dispersion effects via the DFT-D3 scheme with an improved damping function 56 (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 basis set. 57 The Brillouin zone was sampled by a 4x2x1 Γ-centered Monkhorst-Pack k-mesh. 58 Based on the resulting structures, a Si15H16 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 lengths of 1.480 Å. 59 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

Convergence behavior of pEDA terms
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 double zeta (DZ), triple zeta (TZP, TZ2P) and quadruple zeta (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. 62 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, …). 63 The influence of these parameters will be checked in the following sections for molecular and extended systems.

Molecular systems
For the three molecular test systems the results for the convergence w.r.t. basis set size and frozen core approximation are found in Table 1. The influence of the accuracy parameter on the pEDA terms is shown in Table 2. The latter parameter turns out to be well converged with a setting of 5 in comparison to precise computations (accuracy 7) for all systems. The smaller setting (accuracy 3) leads to rather large deviations of the terms especially for the transition metal complex. Note, that the default setting is quite low (accuracy 3.5) and should thus be adjusted for pEDA analyses.
Regarding the basis set, we need to discuss two effects: The influence of the frozen core approximation (small cores in all cases, large core calculations lead to too large errors) and the size of the basis set. As was found for molecular systems in non-periodic calculations before 52 the frozen core approximation introduces only small errors in comparison to the all-electron results. For the TZ2P basis set, these errors are ≤ 0.5 %. At the same time SCF convergence is improved and the computing time is reduced by approximately 30%. For the double-ζ basis set without polarization functions (DZ), the deviations are much larger -notably for the transition metal complex. Here, the frozen core region is too large for the small and inflexible DZ basis set. Table 1   Table 2 The results for the basis set convergence require a closer look. The DZ basis set is in all systems not sufficient to deliver converged pEDA energy terms and large errors remain also due to the missing polarization functions. The TZP and TZ2P basis sets result in differences of < 1.25 % in all cases. The largest basis set in a convergence study is usually taken as the reference -QZ4P in this case. From first glance it looks like the triple zeta results are not converged since deviations of up to 3.6 % persist. But a closer look at the QZ4P calculations reveals several issues with linear dependencies leading to a reduced numerical accuracy. Thus, the TZP and TZ2P basis sets with frozen core approximation are considered the most reliable setup.
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 3). Table 3 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

Extended Systems
For the extended test systems the results for the convergence w.r.t. basis set size and frozen core approximation are found in Table 4. The convergence w.r.t. k-space sampling and accuracy parameter are summarized in Figure 5 for the examples H2 on Cu ( Figure 5a) and C2H2 on Si (Figure 5b). The underlying numbers and results for the other systems can be found in the Supporting Information. Next, we want to discuss the convergence w.r.t the sampling of reciprocal space as shown in Figure 5. It was recognized before, that the convergence of the pEDA terms for a metallic system can depend strongly on the k-space sampling. This is due to the separate orthogonalization steps for occupied and virtual space and discontinuities can be generated due to change of occupations from one k-point to the next. 35a This is reproduced in our case, and for the system H2 on Cu (Figure 5a) a tight sampling of more than 15 k-points per Å -1 is needed to converge the pEDA terms below an accuracy of 1% relative to the finest k-mesh.
Nevertheless, this is not unique to the pED, since it is a general feature of the computation of metallic systems with PBC. Specific for the method here is the different convergence behavior for the pEDA terms. While ΔEelstat converges very quickly (< 5 k-points per Å -1 ), ΔEorb and ΔEPauli are much more strongly dependent on a good sampling of the Brioullin 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 C2H2 on Si(001) (Figure 5b). Around 10 k-points per Å -1 deliver well converged results and even the sampling at the -point is a very good approximation. A closer look at the underlying numbers (see Table S10) reveals an alternating behavior of the results with the usage of a linear (even k-space settings) and quadratic (uneven) tetrahedron method for the sampling (see Method section for details). As suggested, only results from the same method should be used and the quadratic method delivers more accurate results. 65 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 S4).
The accuracy parameter convergence is much less system dependent ( Figure 5). For all test systems, the pEDA terms are converged to <1 % for an accuracy setting of 4 or 5. This is confirmed by the numerical results in the Supporting Information with one exception. The results for CO on MgO(001) show a non-monotonous convergence w.r.t the accuracy parameter mostly due to changes in the ΔEelstat term. It might be concluded that the weak interaction in this system leads to small absolute energy terms which pose a difficulty for deriving well converged analysis results.
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 the next sub-sections to discuss the chemical bonding in the test systems.

4.3
Applications of the pEDA 4.3.1 Dissociative adsorption of H2 on Pd(001) and Cu (001) The dissociative adsorption of H2 is a non-activated process on Pd while an activated process is observed on Cu. This lead to the aforementioned investigation by Philipsen and Baerends (PB) focusing on the role of Pauli repulsion in this adsorption process. 35a It was concluded that the relief of Pauli repulsion is the decisive factor for the different behavior on Pd and Cu. The difference in our implementation is now that we are able to explicitly calculate the ΔEelstat and ΔEPauli term while the ΔE 0 (= ΔEelstat + ΔEPauli) and ΔEorb term 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 H2 approaching the surfaces (Figure 1).
The resulting energy terms along the reaction coordinate (d(M-H2)) are shown in Figure 6 and the respective numbers can be found in the Supporting Information. The same data in the format used by PB can be found in the SI for direct comparison ( Figure S1).

Figure 6
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 6a). Additionally, we can now state that the electrostatic stabilization term ΔEelstat is also more favorable on the Pd surface ( Figure 6b). 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 H2. 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.

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

Adsorption of CO on MgO(001)
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 has been analyzed by the CSOV method in a periodic implementation by Hernandez, Zicovich-Wilson and Sanz (HWS). 44b This analysis confirmed previous findings 66 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 -1 is in very good agreement with the most reliable experimental values. 67 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 5. First of all, the adsorption energy derived is in very good agreement with the above-mentioned high level value of 21.0 kJ mol -1 (Eint = -22.2 kJ mol -1 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, slab setup) but it lends confidence to the pEDA results besides the rather small energy terms. As seen in the previous sub-section, care has to be taken to use well converged computational parameters for such a weakly bound system. Table 5 The bonding picture given by the pEDA shows domination by a strong Pauli repulsion and electrostatic attraction, while the Eorb 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 ΔEPauli and ΔEelstat can be compared to the frozen orbital (FO) term while the ΔEorb 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 approx. 5 kJ mol -1 -but a qualitative difference in the FO term. While the CSOV gives a rather constant value of approx.
8 kJ mol -1 across the three coverage regimes, the pEDA leads to an decrease of the respective sum in going from  = 0.5 (17.5 kJ mol -1 ) to  = 0.125 (4.2 kJ mol -1 ). In consequence, the binding energy (BE) in the CSOV analysis is constant across the three coverages investigated while the pEDA leads to an increase of the interaction energy ΔEint. 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. 68 Thus, we conclude that the bonding in CO on MgO(001) is indeed dominated by Pauli repulsion and electrostatics but the relative importance of these terms shifts upon decreasing coverages and results in a higher binding energy in the low coverage regime.

Adsorption of CO on Si(001)
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(4x2) reconstruction revealing buckled dimers with one Si atom lying above the dimer plane (Siup) and one Si atom below the plane (Sidown). 69 Analysis of the band structure reveals one occupied and one unoccupied surface state. 70 In a localized view on the bonding, this can be understood as the Siup atom exhibiting an electron lone-pair and the Sidown atom an empty orbital. Therefore, this surface dimer is prone to electrophilic (Siup) or nucleophilic (Sidown) attack. 71 Thus, the interaction of the CO molecule with the Sidown atom can be used as a typical example for a donor acceptor bond on a semiconductor surface (Scheme 4a).
This bonding mode was confirmed by experimental and computational studies. 72 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. 72c This is in agreement with the Dewar-Chatt-Duncanson model for transition metal complexes 73 which is known as the Blyholder model in surface science. 74 How can we quantify this view on the bonding of CO on Si(001)? Table 6 The structural optimization of one CO molecule on the Si(001)c(4x2) unit cell (resulting in a coverage of  = 0.25 w.r.t. surface Si atoms) as outlined in the method section leads to a structure which is in quantitative agreement with previous PBE results regarding the bond lengths (Figure 3a). 72c The results for the bonding analysis on this system are shown in Table 6. The bond dissociation energy (Ebond = -102.3 kJ mol -1 ) is higher compared to the previous finding (Eads = −Ebond = 0.92 eV = 88.8 kJ mol -1 ). 72c This difference can be understood by the inclusion of dispersion effects in our study (Edisp = 10.7 kJ mol -1 ) and the dissimilarities in the computational setup (e.g. slab size, basis set). The pEDA results now give a rather high Pauli repulsion term (EPauli = 855.1 kJ mol -1 ) 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 binding to a molecular Lewis acid. 75 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 Sidown atom ( Figure 3a). 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 Si15H16 cluster cut out from the slab setup bearing one CO molecule is shown in Figure 3b. On the one hand, we can now check the numerical accuracy 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 (Edisp) and the preparation energy of the surface (Eprep(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. 76

Adsorption of Acetylene on Si(001)
As a last example, we analyzed the bonding of acetylene to a Si(001)c(4x2) substrate as an example for shared-electron binding between surface and adsorbate (Scheme 4b). This enabled us to test the capabilities of the pEDA to work with unrestricted fragments. This system has been much investigated before experimentally 77 as well as computationally 78 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 4a) to EDA results on a cluster model (Si15H16, Figure 4b) and pEDA results for the same cluster (Table 7). Table 7 The adsorption structure shown in Figure 4a (Table 6) 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 binding 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.