Breaking covalent bonds in the context of the many-body expansion (MBE). I. The purported “first row anomaly” in XH_n (X = C, Si, Ge, Sn; n = 1–4)

The Many-Body Expansion (MBE) is a concept based on combinatorial mathematics that was first introduced over 300 years ago and is usually employed to count the number of elements in the union of finite sets.¹ Its first application to chemical physics problems considered individual water molecules as “bodies” connected via hydrogen bonds as it was used to estimate the non-additive three-body term by partitioning the energy of a water trimer.² Since then the MBE analysis based on the energies of distinct, non-overlapping sub-fragments has been applied to hydrogen bonded clusters by some of us^3–9 and others^10–27 to quantify the importance of non-additive terms in the binding energies of aqueous clusters. Recently, the details of the MBE based on high level electronic structure calculations related to the size of the orbital basis set and the level of electron correlation used in the expansion^28–30 as well as a molecular dynamics protocol based on the MBE (MBE-MD)³¹ were reported. This type of MBE for hydrogen bonded molecular systems (including ions), where the definition of a “body” is straightforward and the system is partitioned in non-overlapping sub-fragments (monomers, dimers, trimers, etc.) by just breaking hydrogen bonds, has laid the foundation for the development of accurate, ab initio based, many-body interaction potentials for water.^32–53 This is to be contrasted to different partitions that are, for instance, based on overlapping sub-fragments, such as the Molecular Tailoring Approach (MTA),^54–59 or other fragmentation based approaches.^60–64 Note that the MBE has also been applied to incorporate molecular orbitals as “bodies” in order to extrapolate the total correlation energy of molecules^65–72 and solids.^73–75 

In this paper, we extend the general idea of the MBE to the breaking of covalent bonds, that is, considering a polyatomic molecule as consisting of a collection of atoms, diatomics, triatomics, etc. A recent study⁷⁶ of the carbon, silicon, and germanium hydrides, XH_n (n = 1–4) based on the spin-coupled generalized valence bond theory has examined the qualitative changes between carbon and silicon/germanium in the H_n−1X–H bond energies (D_e) on the context of the “first raw anomaly”⁷⁷ by expanding on Kutzelnigg’s argument⁷⁸ based on the hybridization of the X atom bond orbitals due to the increase in the spatial separation of the ns and np orbitals between atoms in the first row and the following rows of the Periodic Table. The results presented in our study offer (vide infra) an alternative explanation for this result.

Our analysis is based on casting the atomization energy (ΔE_atomiz⋅) of a molecule of N atoms computed with respect to its constituent atoms i (i = 1, …, N) in their respective ground states (⁠$Ei0$⁠) in the usual way⁴ as

where

and $Eip$⁠, $Eijpq$⁠, $Eijkpqr$⁠, $Eijklpqrs$⁠, etc., are the energies of the in situ (ground or excited) states of atom (i), dimer (ij), trimer (ijk), tetramer (ijkl), etc., in the molecule. It is obvious that if the in situ state of an atom in the molecule (⁠$Eip$⁠) is the atom’s ground state (⁠$Ei0$⁠), the 1-B term for that atom, given by Eq. (2), is zero.

Since the concept of the in situ state of an atom in the molecule is central to our analysis, we will further discuss it using Fig. 1, as discussed by Heitler.⁷⁹ Consider a carbon atom in its ground (³P) atomic state. The pairing of the four (2s²2p²) valence electrons in this electronic state is not appropriate to accommodate bonding with four hydrogen atoms to form CH₄. They should rather be promoted to the ⁵S excited electronic state lying 4.18 eV above⁸⁰ the ³P ground state in order to form the 4 equivalent C–H bonds. In order words, the 2s and the three 2p orbitals of the carbon atom should be hybridized into four equivalent sp³ orbitals forming the atomic ⁵S state, which is the in situ electronic state of the carbon atom in CH₄. Examples of the concept of in situ atomic states have been previously discussed for molecules^81,82 and metal aqueous clusters.^83,84 As it is evident from Eqs. (1)–(4), the choice of the in situ electronic state of an atom will affect the MBE.

The calculations for the ground states of the XH_n hydrides, X = C, Si, Ge, Sn and n = 1–4, were performed at the Coupled Cluster Singles, Doubles and perturbative Triples [CCSD(T) and RCCSD(T)]^85,86 levels of theory. Additionally, two excited states for CH and one for the remaining diatomic and all triatomic molecules were computed. For the CH_n series, we employed Dunning’s augmented correlation consistent basis sets, aug-cc-pVxZ, x = D, T, Q.^87,88 For the SiH_n, GeH_n, and SnH_n species, we employed the aug-cc-pVQZ basis sets for H,^87,88 Si,⁸⁹ and Ge,⁹⁰ and for Sn Peterson’s⁹¹ aug-cc-pVQZ-PP augmented correlation consistent basis sets, which employ accurate small-core (1s²2s²2p⁶3s²3p⁶3d¹⁰; 28 electrons) relativistic pseudopotentials.

The atomization energies, the XH_n−1–H dissociation energies, and the individual many body terms were corrected for basis set superposition error (BSSE)^92,93 as described in Ref. 94. Note that the deformation (or relaxation) energy in that reference that arises from geometrical distortions due to the interaction is zero in the present case since the atomization energy is computed with respect to the atoms. The BSSE-corrected interaction energies (ΔE′_atomiz.) and the BSSE-corrected many-body terms, ΔE′(n-B), e given in Eqs. (5)–(10),⁴ 

etc., where $EGs(M)$ refers to the total energy of the molecule M computed at the geometry G with basis set s.

In order to evaluate the appropriateness of RCCSD(T), which is a single-reference method, we checked the single (t₁) and the double (t₂) amplitudes⁹⁵ as well as the T₁ and D₁ diagnostics.^96,97 For all calculations of the present study, the t₁ and t₂ amplitudes were small and the T₁ and D₁ diagnostics were T₁ < 0.02 and D₁ < 0.04, expect for the ⁴Σ⁻ state of XH, for which the diagnostics were T₁ ∼ 0.03 and D₁ ∼ 0.07. We can, therefore, conclude that the single reference RCCSD(T) method is an appropriate methodology to be employed for the systems in this study. All calculations were carried out with the MOLPRO suite of codes.⁹⁸ 

An important detail of the above analysis is that, for each electronic state of the molecule, the two- and higher-body terms are not necessarily computed at the diatomic, triatomic, etc., ground electronic states but at the respective electronic states of these fragments that are formed from the in situ electronic states of the constituent atoms in the full molecule. We will further elaborate on this important detail in Sec. III when we present the case of XH₂.

The results in this subsection are presented in Figs. 2–13 and Tables I–V. Below we will discuss the individual members of the XH_n series for each n as well as the trends with the atom identity.

XH (X = C, Si, Ge, Sn): The analysis for the ground (²Π), first (⁴Σ⁻), and second (²Δ) excited states of CH is shown in Fig. 2. Energy differences are taken with respect to the ground state of the two atoms, viz., C(³P) + H(²S). The ²Π ground state of CH correlates with the corresponding ground state atomic states, i.e., the in situ atomic state of C in the ground state of CH is the atomic ground ³P state. In other words, the 1-B term is zero since there is no need to promote the ground ³P state of the C atom to form the C–H bond in its ground ²Π state. In this case, the 2-B term is identical with the atomization energy ΔE_atomiz⋅ However, the situation is different for the first and second excited states of the CH molecule. The in situ electronic state of the C atom in the first ⁴Σ⁻ excited state of CH is the atomic ⁵S state, so the (positive) 1-B term is the energy required to promote the C atom from ³P → ⁵S. Accordingly, the 2-B term is the energy difference between the C(⁵S) + H(²S) and CH(⁴Σ⁻) energy levels, which is also the difference between the total atomization energy (ΔE_atomiz⋅) and the 1-B term. Similarly, the in situ electronic state of the C atom in the second ²Δ state if CH is the atomic ²D state and the 1-B and 2-B terms are computed accordingly. The results with the various basis sets (AVDZ, AVTZ, and AVQZ) are denoted with different colors for the different energy levels obtained with the three basis sets in Fig. 2. The individual 2-B terms for the ground (²Π), first (⁴Σ⁻), and second (²Δ) excited states of CH are decreasing in that order (cf. Table I), and the presence of the repulsive 1-B term for the last two determines the overall order in the atomization energies with respect to the atomic ground states, as shown in Fig. 2. The σ bond length of the ⁴Σ⁻ state is shorter than the corresponding value of the X²Π state by 0.03 Å, showing that the bonding of the ⁴Σ⁻ state is stronger than in the X²Π state. The 2-B term of the ⁴Σ⁻ state is significantly larger (double) than the 2-B term of X²Π. This happens because the formation of a σ bond with a highly open shell system, such as C(⁵S), significantly stabilizes the carbon atom and consequently the whole C–H system. On the contrary, the C(³P) is lower in energy than C(⁵S), it has only two unpaired electrons, and, thus, even though the σ bond stabilize the C–H molecule, this stabilization is not dramatically large.

The situation for the ground (²Π) and the first (⁴Σ⁻) excited states of SiH, GeH, and SnH is similar to that for CH; the results for the XH series (X = C, Si, Ge, Sn) are graphically summarized in Fig. 3, and the individual numbers (including the ones corrected for BSSE) are listed in Table I. The decrease of the individual 2-B terms for the respective ground (²Π) and excited (A⁴Σ⁻) states in the XH series monotonically follows the trends of increasing X–H separations (cf. Table V) in an almost linear fashion.

XH₂ (X = C, Si, Ge, Sn): The MBE analysis for the ground (X³B₁) and the first excited (a¹A₁) states of CH₂ is schematically shown in Fig. 4, where the zero in the energy scale is taken with respect to the C(³P) + 2 × H(²S) asymptote. The ground (X³B₁) state of CH₂ correlates with the C(⁵S) + H(²S) + H(²S) atomic states so the 1-B term is positive and corresponds to the ³P → ⁵S promotion energy. The 2-B term is negative, whereas there is a smaller positive 3-B term. The MBE based on the various fragments allows for the further attribution of the individual terms to specific interactions between the constituent atoms. As can be seen from Table II, the main contribution to the total 2-B term (being the sum of the individual interactions between C–H and H–H) comes from the attractive (stabilizing) 2-B_CH term with the remaining 2-B_HH term being quite small and repulsive (destabilizing), amounting to just ∼4% of the total 2-B interaction. In contrast, for the first (a¹A₁) excited state of CH₂, also shown in Fig. 4, the 1-B term is zero (since the in situ electronic state of C in the (a¹A₁) state of CH₂ is the atomic ³P state), and, in this case, both the 2- and 3-B terms are negative. This is expected for the 3-B term for the ¹A₁ state of CH₂ and for the rest of the XH₂ molecules (see below) because this state corresponds to a closed shell system, and, thus, the simultaneous existence of the three atoms further stabilize it. On the contrary, this not the case for the X³B₁ state in which the system is not closed shell and, thus, it follows the usual trend, i.e., the many-body terms oscillate between the positive and negative values.

The analysis for the combined results for the XH₂ series, X = C, Si, Ge, Sn, are schematically shown in Fig. 5 and listed in Table II. A notable difference between CH₂ and the rest of the XH₂ series (X = Si, Ge, Sn) is that the order of the ground/first excited state is reversed. Indeed, as can be seen from Figs. 4 and 5 and the results of Table I, for CH₂, the ground state is the ³B₁ and the first excited the ¹A₁ state, whereas, for X = Si, Ge, Sn, the ground state of XH₂ is the ¹A₁ and the first excited the ³B₁ state. Therefore, the in situ electronic state of the carbon atom in the ground ³B₁ state of CH₂ is C(⁵S), whereas the in situ electronic state of the Si/Ge/Sn in the ground state of SiH₂/GeH₂/SnH₂ is Si/Ge/Sn (³P) (see Figs. 4 and 5). The difference on the type of the ground state results from the fact that the Si, Ge, and Sn atoms are larger than C, forming significant larger X–H bond lengths and this favors the formation of two X–H bonds from the X(³P) rather than from the X(⁵S) states for X = Si, Ge, Sn. Note that the bond arising from the atomic ⁵S state is stabilized via shorter X–H bond lengths than the one arising from the X(³P) atomic state. Additionally, the 2-B_XH and 2-B_HH terms are decreasing in the series for both the ground and the excited states, in line with the corresponding monotonic increase of the X–H bond distance (cf. Table V). Note that the 3-B term for the ground (X³B₁) state is positive (destabilizing), whereas the one for the first excited (a¹A₁) state is negative (stabilizing) in the series. For both states, the 3-B is much smaller than the 2-B term, as seen from Fig. 5 and Table II.

In the following, we elaborate on an important detail of the MBE analysis related to the electronic states used to compute the energies of the fragments. As it can be seen from Figs. 2 and 3, the ²Π ground state of XH (X = C, Si, Ge, Sn) is formed from the respective X(³P) + H(²S) states, whereas the first excited (⁴Σ⁻) state from the X(⁵S) + H(²S) states. The ground ³B₁ state of CH₂ is formed (cf. Fig. 4) from the C(⁵S) + 2H(²S) states and, when evaluating the MBE for that state, the two-body “CH” term is computed with CH at its excited ⁴Σ⁻ state [also formed from C(⁵S) + H(²S)]. In contrast, the first excited a¹A₁ state of CH₂ is formed (cf. Fig. 4) from C(³P) + 2H(²S) and the two-body “CH” term for that state is computed with CH at its ground ²Π state [also formed from C(³P) + H(²S)]. The calculation of the two-body “XH” term in XH₂ for X = Si, Ge, Sn is evaluated in a similar manner while also noting that for these molecules the order of the ground and first excited states flip with respect to CH₂ (see Fig. 5 and Table II). This protocol is also followed in the evaluation of the MBE for the larger molecules.

XH₃ (X = C, Si, Ge, Sn): The MBE for the ground (X²A₁) state of CH₃ is schematically shown in Fig. 6. The ground state of CH₃ correlates with the C(⁵S) + H(²S) + H(²S) + H(²S) atomic states, viz., there is a positive 1-B term corresponding to the promotion ³P → ⁵S for the carbon atom. The 3-B term is smaller than the 2-B and positive, whereas the 4-B is even smaller than the 3-B and negative. The trends for the MBE of the XH₃ series, X = C, Si, Ge, Sn, schematically shown in Fig. 7 and listed in Table III, are similar to the ones for CH₃ with the terms oscillating between the positive (destabilizing) and negative (stabilizing) values, whereas they overall decrease in magnitude with rank k after the 2-B term, as can also be seen from Table III. The MBE, therefore, seems to be “converging” albeit in a slow, oscillating manner. Again, the total 2-B_XH (sum of three identical components) is the major contributor to the 2-B term with the 2-B_HH term being quite small and repulsive as in the XH₂ series; they both decrease in magnitude in the series. With regard to the 3-B term, the signs of the individual components are reversed when compared to the 2-B term: the 3-BX_HH is positive, whereas the 3-B_HHH is very small and negative.

XH₄ (X = C, Si, Ge, Sn): The in situ atomic state of carbon in the ground (X¹A₁) state of CH₄ is ⁵S, giving rise to a positive 1-B energy. The negative (stabilizing) 2-B term is the largest one among the higher order terms, which are oscillating between the positive and negative values and are decreasing in size after the 2-B term as can be seen from Fig. 8. The same trends are observed for the rest of the XH₄ series (X = C, Si, Ge, Sn) as shown in Fig. 9 and Table IV. The 3-B_XHH and 3-B_HHH terms are similar in sign, magnitude, and trend across the XH₄ series as in the XH₃ series. The same behavior is seen for the 4-B_XHHH and 4-B_HHHH terms (cf. Table IV). A notable difference is that the total (positive) 3-B term almost cancels the total (negative) 4-B term. The 5-B term is small (∼10% of the atomization energy).

Trends in the MBE of XH_n (X = C, Si, Ge, Sn): The MBE for the XH_n series is summarized in Figs. 10–13 for X = C, Si, Ge, and Sn, respectively. The trends in both the sign and magnitude of the individual MBE terms are similar across the series, namely, that the 1-B term is zero for n = 1 and positive (energy for the ³P → ⁵S promotion of the carbon atom) for n = 2–4 with the 2-B term being the largest one in the MBE and the rest of the terms oscillating in sign and diminishing after that. The behavior of the MBE is similar across the series. The correlation between the X–H distances and the magnitude of the individual 2-B_XH terms, discussed earlier, is shown in Fig. 14 for the XH_n series. The linear correlation across the XH_n series for each n confirms that the X–H distance is a descriptor of the magnitude of the individual 2-B_XH terms. Note the difference in the magnitude of the 2-B terms for the XH series due to the absence of a 1-B term for this case. This is because the individual 2-B terms are determined with respect to the C(³P) + H(²S) for XH and with respect to the C(⁵S) + n × H(²S) for the XH_n, n = 2, 3, 4 series. It is interesting that this descriptor is independent of n, i.e., the size of the XH_n molecules. The linear fits corresponding to the same in situ atoms are almost identical as shown in Fig. 14 

The variation of the individual 2-B_HH and 2-B_XH terms for the XH_n series, X = C, Si, Ge, Se and n = 2–4 is shown in Fig. 15. In the left panel, which shows the variation of the 2-B_HH across the series, the filled symbols trace states in which X is at the in situ ⁵S atomic state, whereas the open symbols denote states in which X is at the in situ ³P atomic state. The larger relative increase in the 2-B_HH term across the CH_n series is due to the fact that the hydrogen atoms are closer to one another with increasing n as the result of the shorter C–H distances being in the 1.08–1.12 Å range (cf. Table V). In contrast, the 2-B_HH terms for XH_n, where X = Si, Ge, and Sn, are almost the same for n = 2–4, a result of the longer X–H distances (1.50–1.80 Å range, cf. Table V) that bring the hydrogen atoms further apart. The variation of the 2-B_XH term across the series is shown on the right panel of Fig. 15, where the two different curves refer to the calculation when the heavy atom is at the in situ ⁵S and ³P atomic states in the molecule. As a general trend, when X is at the atomic ³P state in the molecule, we observe a monotonic decrease of the 2-B_XH term with the size of the X atom. This is consistent with the fact that the X–H bond distance, already being identified as a descriptor of the 2-B_XH term (cf. Fig. 14), increases across the series. When X is at the in situ atomic ⁵S state, the 2-B_XH is also decreasing across the series, except for Ge. This “anomaly” is attributed to the fact that the ³P → ⁵S energy is the largest one for Ge (cf. Tables I–IV). It is 6.1% larger than C, 12.5% larger than Si, and 8.2% larger than Sn, resulting in the 2-B_Ge-H term of Ge(⁵S) to be further stabilized.

The incremental H_n−1X–H bond energies for X = C, Si, Ge, Sn are shown in Fig. 16. In the left panel of that figure, the energies are computed with respect to the lowest in energy (ground states) of the H_n−1X molecule. Following our previous discussion in Sec. III A, these are ³B₁ for CH₂ and ¹A₁ for SiH₂, GeH₂, and SnH₂ (indicated on the left panel of Fig. 16 for n = 2). The left panel of Fig. 16 is the same as the right panel of Fig. 9 of Ref. 76, but extended for X = Sn. As discussed in this earlier publication and reproduced by our results, the trend in the incremental bond energies with respect to the lowest in energy states is qualitatively different between the CH_n and the rest of XH_n series: indeed, there exists a “peak” for CH₃ (the bond energy for CH₃ is larger than the ones for CH₂ and CH₅), whereas there is a “dip” for the XH_n, X = Si, Ge, Sn, series (the bond energy for XH3 is lower than that for XH₂ and XH₅ when X = Si, Ge, Sn).

However, when the bond energies are computed with respect to the molecular states having the same in situ atomic state (i.e., not necessarily the lowest energy ones but the ones correlating with the same atomic states), the situation is different as shown in the right panel of Fig. 16. In that panel, the various states, with respect to which the individual incremental energies were computed, are noted. The H_n−1X–H incremental bond energies are obtained with respect to the H_n−1X molecule, where X is at either the ⁵S (solid line) or at the ³P (dotted line) atomic states. When this protocol is followed, the variation of the bond energies with n is qualitatively similar across the series and the concept of the “first row anomaly” is not justified. One can nevertheless argue that the (first row) carbon atom can be still considered different than the (second and subsequent rows) Si, Ge, and Sn atoms in lieu of the “flipping” of the ground/excited states for n = 2 in the series (CH₂ vs SiH₂/GeH₂/SnH₂). However, this qualitative change is not a property of the individual atoms per se but rather their different behavior in just the XH₂, but not the rest of the hydrides.

We have described a novel decomposition scheme for the MBE in molecules and applied it to the XH_n series (X = C, Si, Ge, Sn; n = 1–4). The protocol allows for the decomposition of the atomization energy in terms of atoms, dimers, trimers, etc. In the present implementation of the MBE that is based on the breaking of covalent bonds to define subsystems, the 1-B term is both qualitative and quantitatively different than the analogous 1-B term that has been popularized in the MBE of aqueous ionic systems. When breaking hydrogen bonds to define the various subsystems, the 1-B term represents the geometrical distortion⁹⁴ of the individual fragments from their gas phase geometries due to the interaction with other neighboring molecules or ions. It amounts to a few kcal/mol and it can be indirectly probed experimentally by infrared (IR) spectroscopy,^99,100 which records the change in the position of the individual vibrational bands from the gas phase as a result of the change in the fragment’s geometry from the isolated species. In contrast, when breaking covalent bonds, the 1-B represents the electronic excitation of an individual atom¹⁰¹ to bring it to its in situ electronic state in the molecule, and it amounts to a few eV (i.e., it is about two orders of magnitude larger) and can be probed by Ultraviolet–Visible (UV) spectroscopy.

We have found that the MBE in the XH_n series is oscillating between positive and negative values, and it is converging with increasing rank of the expansion. Among the individual terms, the 2-B is by far the largest one in the MBE. The X–H distance is a good descriptor of the strength of the 2-B_XH term, and its variation across the series can be attributed to the respective geometrical changes. The analysis can offer an alternative explanation for the purported “first row anomaly” that is based upon the different variation of the incremental bond energies for CH_n compared to XH_n, X = Si, Ge, Sn, when these bond energies are evaluated with respect to the lowest energies (ground states) of the members in the series. However, there is a “flipping” between the ground and the first excited state in CH₂ (ground is the ³B₁ and excited is the ¹A₁ states) compared to XH₂, where X = Si, Ge, and Sn (ground is the ¹A₁ and excited is the ³B₁ state). When these incremental binding energies are evaluated with respect to the molecular states that have the same in situ atomic state, their variation is the same with n for all X (C, Si, Ge, Sn) in the XH_n series and the concept of the “first row anomaly” is, thus, not justified.

Note that the present analysis is simplified due to the presence of just one heavy atom (C, Si, Ge, Sn) in each of the systems studied and can be much more complex when more than one heavy atom is involved. Nevertheless, it provides a straightforward extension of the popular MBE for hydrogen bonded systems to incorporate the breaking of covalent bonds and offers valuable insights into the chemical bonding of chemical systems.

S.S.X. was supported from the Center for Scalable Predictive methods for Excitations and Correlated phenomena (SPEC), which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division as part of the Computational Chemical Sciences (CCS) program at Pacific Northwest National Laboratory. Battelle operates the Pacific Northwest National Laboratory for the U.S. Department of Energy. This research used computer resources provided by the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.