Correlations between Density-Based Bond Orders and Orbital-Based Bond Energies for Chemical Bonding Analysis.

Quantum chemistry-based codes and methods provide valuable computational tools to estimate reaction energetics and elucidate reaction mechanisms. Electronic structure methods allow directly studying the chemical transformations in molecular systems involving breaking and making of chemical bonds and the associated changes in the electronic structure. The link between the electronic structure and chemical bonding can be provided through the crystal orbital Hamilton population (COHP) analysis that allows quantifying the bond strength by computing Hamilton-weighted populations of localized atomic orbitals. Another important parameter reflecting the nature and strength of a chemical bond is the bond order that can be assessed by the density derived electrostatic and chemical (DDEC6) method which relies on an electron and spin density-partitioning scheme. Herein, we describe a linear correlation that can be established between the DDEC6-derived bond orders and the bond strengths computed with the COHP formalism. We demonstrate that within defined boundaries, the COHP-derived bond strengths can be consistently compared among each other and linked to the DDEC6-derived bond orders independent of the used model. The validity of these correlations and the effective model independence of the electronic descriptors are demonstrated for a variety of gas-phase chemical systems, featuring different types of chemical bonds. Furthermore, the applicability of the derived correlations to the description of complex reaction paths in periodic systems is demonstrated by considering the zeolite-catalyzed Diels-Alder cycloaddition reaction between 2,5-dimethylfuran and ethylene.

Introduction

Modern computational chemistry provides a powerful toolbox for studying the fundamental aspects of chemical bonding and chemical reactivity.[1−4] A wide range of practical methodologies has been developed so far to investigate the electronic aspects of chemical bonding, making use of electron density partitioning schemes or based on the direct analysis of the electronic wavefunctions. For example, the quantum theory of atoms in molecules (QTAIM) provides a framework for the topological analysis of the electron density.[5] Local descriptors such as the electron, Laplacian, and energy densities can be computed in the framework of QTAIM and utilized for the qualitative and quantitative analysis of chemical bonds.[6,7] The potential electron density has been shown recently to be particularly useful for estimating the effective force constants of a chemical bond as a measure of the bond strength.[8,9] An alternative partitioning scheme is the density derived electrostatic and chemical (DDEC6) method[10−12] which involves spherical averaging of the atomic electron densities. In this method, the dressed exchange hole approach is employed to compute the DDEC6-based bond orders (BOs) that can be regarded as quantitative descriptors, reflecting the strength of the chemical bonds. The energy decomposition analysis (EDA) partitions the interaction between a pair of atoms into energy components, namely, the electrostatic, polarization, charge transfer, exchange, and correlation contributions by referencing the wavefunction and electron density of the chemical system to those of the isolated reference ions.[13,14] The crystal orbital overlap population (COOP) analysis put forward by Hofmann[15] and the crystal orbital Hamilton population (COHP) analysis introduced by Dronskowski and co-workers[16] utilize the electronic wavefunction to derive the bonding information. The COOP and COHP schemes enable direct quantification of the (anti)bonding orbital overlap and the strength of interatomic bonds, respectively.

Among the different available bonding analysis formalisms, the COHP and DDEC6 methods are particularly attractive for practical applications in computational chemistry, given the chemically intuitive nature of the respective bond quantifiers. Recent studies demonstrate the power of these approaches for the theoretical analysis of catalytic reactions and scaling laws on transition-metal surfaces,[17,18] transition-metal oxides,[19] and zeolites.[20]

The DDEC6 charge partitioning yields consistently accurate results for a wide range of materials and bonding types. The DDEC6 methodology assigns atomic electron and spin distributions to each atom in a chemical system.[12] This approach provides a number of important advantages over other available methods because (1) it avoids the assumption of a constant BO to electron density overlap,[21,22] (2) does not require the use of the method-dependent first-order density matrix,[23−26] (3) does not use the bonding/antibonding orbital occupancies which fail for longer bonds,[27] and (4) avoids the computationally expensive exchange–correlation hole partitioning approach.[28] The COHP yields complementary information to the density of states (DOS). Although the DOS provides insight into the probability of finding an electron in a particular atomic orbital as a function of electron energy, the COHP enables one to identify whether the respective electron contributes to a bonding, antibonding, or nonbonding interaction.[16] The COHP formalism enables a direct quantification of these energy-partitioned contributions by using the Hamiltonian.

Recently, we thoroughly investigated the Diels–Alder cycloaddition (DAC) between 2,5-dimethylfuran (DMF) and ethylene over third-row d-block-[29] and alkali-exchanged[30] faujasite models using both the DDEC6 and COHP methods. In the former study,[29] we found a qualitative correlation between the integrated COHP (ICOHP) and BO values computed for the interaction between the various d-block cations and the carbon atoms of either DMF or ethylene. In the latter study,[30] we found a qualitative trend between the ICOHP-computed interaction strength between the active sites and the carbon atoms of DMF in line with the Lewis acidity of the alkali cations, albeit these interactions were ionic for which the COHP analysis is less well suited.[16] Furthermore, it has been demonstrated that the variation of the ICOHP values and the DDEC6-derived BOs with interatomic distances generally follows the trend of the heuristic Pauling BOs.[12,31,32] Therefore, the observed trends and the notion that one could potentially correlate chemical intuitive BOs with chemical bond strength, make an exploration of the correlation between the ICOHP-computed bond strength and DDEC6-derived BO (ICOHP–BO correlation) appealing.

However, such an effort is hampered by the fundamental problem in which periodic models share no absolute energy reference point when using the COHP formalism.[16,33,34] The lack of an absolute energy reference makes the COHP method model dependent and thus prohibits a comparison of bond strengths computed for pairs of atoms in different structural models. On the other hand, the DDEC6-derived BOs are relatively insensitive to the choice of basis sets and exchange–correlation functionals and can directly be compared between the different systems.[12] Considering these notions, it is thus important to investigate methodological and chemical limitations of a possible ICOHP–BO correlation.

Herein, we report on a theoretical and quantum chemical investigation which aimed at establishing the ICOHP–BO correlation and exploring the framework in which such a correlation is possible. The ultimate goal is to use this correlation to describe and compare the changes to interatomic bonds in catalytic reactions. To this end, we will first introduce the COHP- and DDEC6-derived BO formalisms. Then, gas-phase molecules will be studied to explore the boundaries within which the correlation remains valid. Finally, the validity of such correlation is evaluated for the case of the zeolite-catalyzed DAC between DMF and ethylene in periodic zeolite models of different chemical compositions.

Introduction to the Methods

COHP Analysis

The COOP and COHP analyses allow the description of bonding in molecules and solids. They allow for the deconvolution of the band structure into atomic orbitals, quantify the degree of net orbital overlap, and they can also be used to determine the bond strength.

In order to apply methods[35−39] using localized basis sets to analyze material properties for plane wave-based calculations, projection schemes were introduced.[40,41] Also, within the COHP method employed in this work, the plane-wave (PW) wavefunction is transformed into a wavefunction based on a combination of localized atomic orbitals (LCAOs):[16,36,37,42,43]where in eq , n is the nth basis function in a set with m total basis functions that all make up the wavefunction, J is the band number, and c is the the orbital coefficient matrix of the orbital φ. Consequently, the energy-partitioned band structure can be transformed into orbital pair contributions to obtain a localized DOS for atoms i and j, Figure a. This is also called the projected DOS (pDOS):

Figure 1

The pDOS of N2 split into pDOSs of every individual N-atom (red = 2s and black = 2p) in (a). The total pCOOP and −pCOHP of the nitrogen molecule are shown in (b,c), respectively.

Using the overlap population (P):one obtains an overlap population-weighted pDOS. Here, S = ⟨φ|φ⟩ is the overlap of atomic orbitals φ and φ.

Because the DOS is an energy-partitioned property, the overlap population-weighted pDOS shares this property. Using this and P one can define regions where the atomic orbital overlap is bonding, antibonding, and nonbonding in nature. The resulting function is called the COOP (COOP(E)) as introduced by Hoffmann[15] (Figure b): in which f is the occupation number of each band J. Integration of the COOP(E)-function up to the Fermi level will yield the net orbital overlap between atoms i and j. As P is used, the COOP(E) function is basis-set dependent.[34]

The COOP function can be rewritten by replacing the overlap matrix with the Hamiltonian matrix. By convention, the COHP function is defined as COHP(E) and can be computed according to:where H represents the Hamiltonian matrix element between atomic orbitals φ and φ, and c and c are the coefficients of these atomic orbitals in the molecular orbital Ψ (Ψ = ∑cφ). A positive value for −COHP(E) symbolizes a bonding electronic interaction between the atomic orbitals i and j, whereas a negative value describes an antibonding interaction. A value of zero is associated with a nonbonding interaction (see for instance Figure c). The integrated value of −COHP(E), ICOHP, is a measure for the bond strength. This formulation provides a good approximation of the bond energy as long as the repulsive energy of the nuclei is canceled by the double-counted electrostatic interactions.[44] Note, though, that the energy computed with eq , accounts for pair-wise interactions but does not account for many-particle interactions, which may also influence the strength of the interatomic bond under study.

Within the COHP formalism, the lack of an absolute zero energy reference prohibits one to compare ICOHP values obtained in different structural models directly with each other. This can be appreciated by examining eq , which is an expression for the crystal’s total cohesive energy (Ecoh), obtained by subtracting the total energy of the atoms from the total energy of the crystal. The equation is adopted from ref (16), which also contains the full derivation:

In the above equation, ε refers to occupied one-electron eigenvalues. In combination with f, summation over all J (third term) yields the band structure energy. The first term is the off-site COHP. The second term is the on-site COHP and the third term is the band structure energy of the atom. The fourth, fifth, and sixth relate to the charge density difference, exchange–correlation difference, and Madelung term, respectively. Note that the third term arises from the summation of the total energy of a single reference atom whose energy is determined by the occupation of J bands each with energy ε. The fourth and fifth terms are the differences in Coulomb and exchange–correlation energies between the separate reference atoms and that of the actual crystal under study, respectively. The sixth term is the nucleus–nucleus repulsion term arising from the Schrödinger equation used for the crystal.

The importance of eq is that it tells us that the bond energies (off-site COHP) are only the real bond energies if the second and third lines in eq cancel out exactly. Additionally, the second, third, and fourth terms carry undetermined constants. We can, therefore, rewrite eq essentially as:where C is the total sum of the errors that not exactly cancel out each other and the undetermined constants.

Density-Derived Electrostatic and Chemical Method

Molecular bonding can be studied in terms of the chemically intuitive BOs. The density-derived electrostatic and chemical method was introduced by Manz and Limas in 2016 (DDEC6)[10] and is a revised version of its predecessors[45,46] (e.g. DDEC3). The DDEC6-based BOs can provide one ab initio BOs without the assumption of constant BOs to atomic charge-density overlap ratios.[12] The exact derivation of the equations necessary to both compute the BOs and execute the underlying density derived electrostatic and chemical (DDEC6) charge partitioning is explained elsewhere.[10−12] Fundamentally, formation of a bond is assumed to arise from electron exchange between two atoms close enough to exhibit overlapping electron densities.

Manz defined the BO of an atom pair A (in the unit cell) and j (atoms in both unit cell and periodic images):[12]where B is the BO between atom A and j, CE is the contact exchange, and Λ is the dressed exchange hole delocalization term. The term CE describes the electron exchange between atoms A and j in a material, formulated in:where any ρ⃗avg is the average spherical electron density of atom i as a function of the atomic electron distribution and atomic spin magnetization density vector obtained through DDEC6-based partitioning of the electron density. The term ρ⃗avg is the sum of all ρ⃗avg found in the material (unit cell + periodic images). Note that, this equation deals with the dressed exchange hole, which is an adjusted (either more contracted or more diffuse) exchange hole to obtain more accurate BOs. The second term in eq is the dressed exchange hole delocalization term, defined according to eq :where Χcoord.nr. accounts for coordination number effects, Χpairwise for pair-wise interactions, and Χcon. is a constraint on the density-derived localization index, B. The latter is a matrix that equals the total number of the dressed exchange electrons in the material (unit cell + periodic images). These terms are constraints and scaling relationships to keep the BOs well behaved.

Computational Details

Models

The first model was a box with cell edges of 20 × 20 × 20 Å3. The molecules were located in the center of the periodic box. An N2 molecule was always located on one of the cell edges as a reference. The studied molecules were: ethane (1), ethylene (2), acetylene (3), propylene (4), cyclooctyne (5), 1,3,5-hexatriene (6), hexa-1,5-dien-3-yne (7), benzene (8), toluene (9), para-xylene (10), benzoic acid (11), terephthalic acid (12), benzodithioic acid (13), benzene-1,4-bis(carbodithioic) acid (14), furan (15), 2-methylfuran (16), DMF (17), furan-2-carboxylic acid (18), furan-2,5-dicarboxylicacid (19), furan-2-carbodithioic acid (20), furan-2,5-bis(carbodithioic) acid (21), acetonitrile (22), hydrogen cyanide (23), cyanogen fluoride (24), cyanogen chloride (25), cyanogen bromide (26), cyanogen iodide (27), cyanogen (28), acrylonitrile (29), acetic acid (30), oxaclic acid (internal H-bond) (31), ethanedithioic acid (32), ethanebis(dithioic) acid (33), and pyridine-2,5-dicarboxylic acid (34), pyridine (35), pyrazine (36), piperidine (37), piperazine (38), benzamidine (39), terephthalamidine (40), aniline (41), and dihydropyrazine (42). Atom indices can be found in Figure .

Figure 2

Molecular library including 42 molecules (part I) and the Diels–Alder reactive system (part II) represented by the intermediates and TSs involved in the DAC reaction between DMF and C2H4 analyzed in this work. Black Arabic numbers between brackets are the molecular indices. Gray, blue, red, and pink Arabic numbers are the carbon, nitrogen, oxygen, and sulfur atom indices, respectively. For obvious cases, the atom indices are omitted. Red or Blue Roman numerals indicate the various groups as used in this work. Some groups are used twice to analyze carbon heteroelement/carbon–carbon bonds, that is, IX/VII and XI/XII.

Different groups of bonds were defined on the basis of the class of molecules (e.g. furanics vs aromatics) and the presence of heteroatoms (e.g. oxygen vs sulfur). The former accounts for different stoichiometries, for example, the carbon-to-oxygen ratio changes the number of valence electrons on a per element basis. The second accounts for a variation in the valence principle quantum number of the atoms involved in the ICOHP analysis. The groups were (I) CC bonds in hydrocarbon molecules including aromatic compounds containing only carbon and hydrogen atom groups (1–10). This group is referred to as the H,C-only hydrocarbon group. (II) CC bonds in molecules containing carboxylic acid functionalities (11, 12, 30, and 31), (III) CC bonds in molecules (13, 14, 32, and 33) containing dithioic acid functions, (IV) CO bonds in furanic compounds (15–21), (V) CO bonds in aromatic compounds containing carboxylic acid functionalities (11, 12, and 34), (VI) CC bonds in furanic compounds containing carboxylic acid functionalities (15–19), (VII) CC bonds in furanic compounds having dithioic acid functions (20 and 21), (VIII) CS bonds in aromatic compounds containing dithioic acid functionalities (13 and 14), (IX) CS bonds in furanic compounds containing dithioic acid functionalities (20 and 21), (X) CN bonds in halogen cyanides, acrylonitrile, and acetonitrile (22–29), (XI) CC bonds in N-heterocyclic cycles, cyanogen, and acrylonitrile (28, 29, 34–42), and (XII) CN bonds in N-heterocyclic cycles, cyanogen, and acrylonitrile (28, 29, 34–42). The members of each group can be looked up in Figure and the Supporting Information, Tables S1–S12.

The NN bond in the dinitrogen molecule placed at a large distance from the investigated molecule in the same periodic box served as a reference for the COHP-computed bond strengths. We hypothesized that for an adiabatic system of noninteracting molecules, the strength of the NN bond in a noninteracting N2 molecule should always be the same, irrespective of the chemical composition of the system.

The second and third types of models were adopted from previous work.[29,30] These were periodic rhombohedral low-silica alkali-exchanged faujasite[30] and the high-silica third-row d-block cation-exchanged faujasite[29] models, respectively. Briefly, the low-silica alkali-exchanged zeolites (Si/Al = 2.4, Si34Al14O96M14, M = Li+, Na+, K+, Rb+, Cs+) are characterized by a high accessible active site density in the faujasite supercage and are referred to as MY. The high-silica third-row d-block cation-exchanged faujasites hold a single active site and an appropriate amount of framework aluminum substitutions to compensate for the charge of the d-block cation. These models are referred to as TMFAU (TM = Cu(I), Cu(II), Zn(II), Ni(II), Cr(III), Sc(III), and V(V)).

Diels–Alder Cycloaddition

The DAC reaction between DMF and ethylene in both TMFAU and MY has been studied in previous studies using periodic DFT calculations.[29,30,47,48] We directly used these models in our current work. Briefly, the initial state (IS) consists of DMF and ethylene both coadsorbed in either TMFAU or MY. The IS is referred to as 1DAC. The DAC transition state (TS) is referred to as TS and the adsorbed cycloadduct (FS) as 2DAC. During the DAC reaction, three π-bonds are converted into two σ-bonds and one π-bond. The two σ-bonds are completely new bonds, which do not exist in 1DAC. As this reaction has been analyzed in different chemical environments and involves several bonds undergoing significant changes, it provides ample opportunity to investigate both the scaling and reproducibility of the ICOHP and BO analyses. Examples of some of the evaluated structures are given in Figure . The IS (1DAC/Cu(I)FAU) and TS (TS/Cu(I)FAU) of the synchronous concerted DAC reaction over Cu(I)FAU are shown in panels (a) and (b), respectively. Panel (c) displays the first TS of the two-step DAC reaction over Cu(II)FAU (TS1/Cu(II)FAU). 1DAC/LiY and 1DAC/KY are shown in panels (d,e), respectively and the DAC TS in KY in panel (f).

Figure 3

Example geometries for 1DAC/Cu(I)FAU and TS/Cu(I)FAU shown in (a,b), respectively. The first TS of the two-step DAC pathway over Cu(II)FAU is shown in panel (c). Selected cuts from the periodic unit cell from 1DAC/LiY and 1DAC/KY in (d,e), respectively. The DAC TS in KY is shown in (f).

Electronic Structure Calculations

Periodic DFT calculations were performed with the Vienna Ab initio Simulation Package (VASP).[49−51] For all the systems, the k-point mesh was set to the Gamma point. The cut-off energy was 500 eV, employing a plane-wave basis set. To approximate the exchange and correlation energy, the PBE-functional was used.[52] This was complemented by the projected augmented wave scheme to describe the electron–ion interactions.[53] The DFT-D3 method with Becke–Johnson damping was used to account for long-range dispersive interactions.[54,55] The gas-phase models were optimized from scratch. The root-mean-square force convergence criterion was set to 0.015 eV/Å. The TMFAU and MY zeolite models were adopted from previous work[29,30] and subjected to a single-point calculation only to obtain the wavefunction and the electron density.

COHP-Analysis

The COHP analysis was performed with the Lobster 2.2.1 code, upon a transformation of the (plane) wave functions from VASP into a localized basis set (STO).[16,36,37,42,43] In addition, an automatic rotation of the basis set was applied. The pair-wise interatomic interaction strength was computed by integrating the COHP up to the Fermi level (ICOHP). For a proper COHP analysis, the number of bands was set to the total number of orbitals present in the model in each calculation.

BO Analysis

BOs were analyzed using the Chargemol code.[56] To obtain accurate electron densities, the VASP calculations were performed using a 2.5 times increased fast Fourier transform (FFT) grid density. This grid is more than sufficient for accurate BO analysis. The effects of plane-wave energy cutoff, k-point mesh, and FFT grid spacing on computed DDEC6 properties have been studied in detail by Limas and Manz.[57]

Results

Bond Strength Quantification in Gas-Phase Molecules

The first part of the study focused solely on the molecular library including 42 gas-phase molecules. Here, we specifically targeted the possible ICOHP–BO correlation. The dependency of the ICOHP–BO correlation on the chemical composition of the evaluated molecules was also investigated. The results are shown in Figure . The fitted parameters are reported in Table .

Figure 4

ICOHP–BO correlation for groups I to XII. The inset shows the residuals of each fit. The lines represent the linear fits.

Table 1

Linear Fits a × BO + b of Each Group

id. no.	a (eV/BO)	aerror (eV/BO)	b (eV)	berror (eV)	R2-adj
I	–5.04	0.13	–1.12	0.21	0.98
II	–2.6	0.16	–5.39	0.22	0.94
III	–2.81	0.40	–5.08	0.57	0.75
IV	–5.89	0.3	–2.95	0.43	0.95
V	–6.67	0.12	–1.72	0.22	0.99
VI	–3.24	0.91	–4.52	1.29	0.37
VII	–2.5	0.78	–6.03	1.1	0.54
VIII	–5.01	0.04	–1.87	0.05	0.99
IX	–4.93	0.13	–1.81	0.19	0.99
X	–5.43	0.53	–3.10	1.42	0.86
XI	–4.72	0.39	–2.09	0.56	0.79
XII	–7.27	0.14	1.26	0.25	0.99

The residuals from Figure imply that a linear fit is appropriate to describe the data. The standard errors of the fits of groups I, II, V, VIII, IX, and XII are found to be smallest with values of 0.13, 0.16, 0.12, 0.04, 0.13, and 0.12 eV/BO, respectively (Table ). In contrast, groups III, IV, VI, VII, X, and XI exhibit lower R2-adjusted values accompanied by the larger errors, being 0.40, 0.3, 0.91, 0.78, 0.53, and 0.39 eV/BO, respectively.

The good fit of group I is attributed to the fact that the molecules of this group only contain carbon and hydrogen without the presence of functional groups that could potentially affect the CC bond strength via electron delocalization or ionic interactions. Interestingly, the large range of the H/C ratio (3 in 1 and 1 in 8) has no significant influence on the ICOHP–BO correlation. This is possibly because hydrogen atoms cannot induce too significant changes to the electronic structure of CC bonds. The poor fits for groups VI and VII indicate that their ICOHP–BO correlations are sensitive to changes in the stoichiometry. For instance, in group VI, the O/C ratio in 15 is 1/4 and increases to 0.83 in 19. Additionally, group VI contains molecules with and without functional side groups, that is, methyl and carboxylic acid site groups. The inductive effect of these functional groups is different and will affect the electronic structure of the CC bonds significantly, thus reducing the quality of the ICOHP–BO correlation. Another example is group VII. Similar S/C ratios are found for group VII as compared to those in group VI, but a third-row element is also present within the molecules. The addition of a carbodithioic acid side group (20 → 21) increases the sulfur content significantly. This adds extra orbitals with a principal quantum number of 3 (3s and 3p), potentially shifting the reference energy and/or the unknown constant C. Because of this, the CC bond ICOHP–BO correlation does not hold well anymore. Therefore, although the aforementioned O/C ratios in groups VI and VII are smaller than the H/C ratio in group I, it is hypothesized that the introduction of a third element and the presence of functional side groups affect the electronic structure of the molecules and thus reduces the quality of the ICOHP–BO correlation. The above observations might explain why group II exhibits such a surprisingly good fit. All its members contain carboxylic acid groups and only the elements C and O.

Furthermore, fits of the ICOHP–BO correlations for groups IV, V, VIII, and IX are good. We note that these groups only consist of a limited set of CS and CO bonds such that each of these fits essentially interpolates two points. However, the data in these groups illustrate the reproducibility of the BO and ICOHP analyses and with all other groups indicating linear correlations; we are confident that a linear interpolation is acceptable.

The poor fit for group X, as illustrated by the large error value, is caused by the presence of group 17 elements (denoted as Hal). Although the various CHal bonds in the cyanogen halides exhibit a trend by themselves as a function of their valence shell principle quantum number (data not shown), the cyanogen halides cannot be used to fit CN bonds. Rather, a good correlation could have been obtained when using acrylonitrile or cyanide derivatives as members of group X. This statement is supported by the good fit of group XII for which only small standard deviations are found for the slope and intercept (0.12 and 0.25 eV/BO, respectively, with an R2-adjusted of 0.99). We note that the poor fit for group XI is explained along the same lines as those for groups VI and VII. The nitrogen atoms are present at different positions within the cyclic molecules and are sometimes part of a functional group. Therefore, although the ICOHP–BO correlation for CN bonds in XII holds well, that for CC bonds in XI breaks because of inductive effects. From these data, we infer that inductive effects of functional side groups reduce the ICOHP–BO correlations for CC bonds involving β and γ carbon atoms. The ICOHP values for CN, CO, and CS bonds can consistently be correlated to BO values.

We also note that none of the ICOHP–BO slopes (Figure ) intersect with the point (0, 0). This is an unexpected observation. We expected the BO and ICOHP values to reach zero when the interatomic distance is infinitely large. More research is necessary to elucidate the meaning of this nonzero intercept.

In an attempt to increase the quality of the fits, the N2 reference molecule was used. We took one ICOHP value of the NN bond in N2 within each group against which we referenced all other members of the respective group. The reference compounds are marked with an asterisk in Tables S1–S12 in the Supporting Information. The results are shown in Table S13. Scaling within each group was performed by computing the deviation (in percentage) of each NN bond with respect to the NN reference bond. Subsequently, the CC/CN/CS/CO bonds were scaled with the same percentage. This procedure yielded hardly any changes to the fitting parameters.

The established correlations allow either the ICOHP or BO value to be estimated when the other quantity is known. For instance, ethane (1), ethylene (2), and acetylene (3) have BOs (ICOHP values) of 1.224 (−6.64 eV), 1.706 (−8.86 eV), and 2.993 (−15.45 eV) respectively. Propylene (4) has a BO of 1.968 (−10.6 eV) for the C1–C2 bond. In benzene (8), each CC bond has a BO of ca. 1.552 (−9.1 eV). Other H,C-only molecules with conjugated bonds involve hexatriene (6) and hexa-1,5-dien-3-yne (7). The C1–C2 and C3–C4 bonds in 6 are characterized by BO values of 2.594 and 2.364 and ICOHP values of −14.59 and −14.19 eV, respectively. The C1–C2 bond in 7 has a BO of 1.945 and an ICOHP of −10.54 eV. On the basis of these findings, we can define regimes within the ICOHP–BO correlation of H,C-only hydrocarbons. Namely, a BO of 1–1.5 is characterized by an ICOHP of ca. −6 to −9 eV; CC bonds with a strength of −9 to −12 eV have BOs of ca. 1.5–2; and the CC bonds with BOs of 2.5 and 3 are characterized by ICOHP values of −14 and −16 eV, respectively. Such regimes can also be defined for all of the other studied chemical bonds.

Minor changes to carbon–carbon BOs via the addition of substituents are more challenging to probe with the ICOHP method (or vice versa). Adding a methyl side group to benzene (8) changes it into toluene (9). The result is that the C1–C2/C1–C6 bonds in 9 have a BO 0.14 lower than those in 8 (BO = 1.55). Addition of the second methyl group at the para position yields para-xylene (10) and results in the same effect for the C3–C4/C4–C5 bonds. Lowering of the BO values results in a consistent and concomitant bond weakening; the absolute ICOHP bond strength becomes 0.2–0.3 eV smaller. Additionally, the BOs of the C2–C3/C5–C6 bonds increase marginally but consistently with ca. 0.01 per added methyl substituents. However, the changes to these BO values are too small to provide reliable changes in ICOHP values.

As the fits in Figure illustrated, the ICOHP–BO correlation changes upon varying the stoichiometry. For instance, replacing the methyl substituents by carboxylic acid functionalities shifts the ICOHPCC values to more negative values. The C2–C3/C5–C6 bonds have BO values of approximately 1.569/1.565 in benzoic acid (11) and 1.571 each in terephthalic acid (12), ca. 0.01 higher than in benzene. However, ICOHP values shift to −9.38/–9.36 and −9.55/–9.54 eV. The respective ICOHP values change with ca. 0.34 and 0.53 eV with respect to the ICOHP values found in benzene. Replacement of the carboxylic acid groups by dithioic acid functionalities hardly shifts the bond energies for the same BO. Replacement of oxygen by sulfur results in C2–C3/C5–C6 BOs of 1.578/1.577 with the corresponding ICOHP values of −9.42/–9.42 eV in benzodithioic acid (13). These BOs are 0.01 higher than in 11 and equal to those in 12. ICOHP values have changed with only −0.06 and +0.1 eV. Yet, in benzene(bis)carbodithoic acid (14), the C2–C3/C5–C6 BOs increase with only approximately 0.02 whilst the ICOHP values increase with ca. 0.16/0.22 to −9.58/–9.66 eV. Furthermore, the CC bonds in ethanedithioic acid (35) and ethane(bis)dithioic acid (36) are 1.233 (−7.78 eV) and 1.169 (−8.66 eV), respectively. These seem to be best related to the C1–C7 and C4–C8 bonds in 13 and 14 of which the BOs are 1.164 (−8.19) and 1.167 (−8.39 eV). However, 32 and 33 differ clearly from the correlation found in 13 and 14. The relative amount of sulfur increased from C/S = 0.5 in 33 to C/S = 3.5 and 2 in 13 and 14.

Summarizing, the change in the ICOHP value as a function of increasing BO can be studied between different models directly, provided the stoichiometry remains relatively similar. CC and CO bonds can be compared when originating from the same compound as they come from the same model with the same unknown total constant C. Inductive effects of functional side groups reduce the ICOHP–BO correlations for CC bonds involving β and γ carbon atoms. The ICOHP values for CN, CO, and CS bonds can consistently be correlated to BO values. Importantly, the trends displayed in Figure allow one to distinguish between bonds with different BOs, once the ICOHP value is known. In addition to these findings, we define four prerequisites:

An ICOHP–BO correlation of good quality can only be established when the stoichiometry of the evaluated molecules does not vary significantly. For instance, CC bonds in H,C-only hydrocarbons and furanic compounds should not be compared directly. The presence of the oxygen atom changes the ICOHP–BO correlation.

ICOHP–BO correlation of CC bonds involving β and γ carbon atoms can only be constructed when the number of functional groups or the stoichiometry within the functional group remains similar. The functional side groups cause inductive effects which negatively affect the quality of the ICOHP–BO correlation for the CC bonds involving β and γ carbon atoms.

ICOHP–BO trends are approximately linear for a species containing CHet bonds with Het being a heteroelement whose principle quantum number changes (i.e., going down a group). However, there is not necessarily a strong ICOHP–BO correlation for the other bond types in these molecules. For instance, there is an ICOHP–BO correlation for carbon–halide bonds in cyanogen halides, but the CN bonds in these cyanogen halides do not exhibit a strong ICOHP–BO correlation.

The changes in ICOHP /BO values have to be sufficiently large (ΔBO ≈ 0.25). In the evaluated trends presented here, CC bonds in H,C-only hydrocarbons can be studied with greater accuracy and with smaller BO-margins (in the order of ΔBO ≈ 0.2) than most other evaluated trends.

Studying Bond Evolution in Chemical Reactions

The second part of the study was dedicated to the investigation of the reproducibility and BO–ICOHP correlation on a more practical example of the DAC reaction between DMF and C2H4 in periodic models of cation-exchanged faujasite zeolites.

The mechanism of the DAC/D reaction between 2,5-DMF and ethylene in alkali (Li, Na, K, Rb, and Cs)-exchanged faujasite catalysts was investigated in previous work.[47] Two models were used, namely, the isolated-site high-silica and low-silica faujasite model containing a high density of accessible active sites representative of the as-synthesized catalyst. The results indicated that the DAC reactivity trend was inverted in the second more realistic model as compared to the single-site model.[47] These reactivity differences were rationalized in a follow-up work based on an in-depth electronic structure analysis.[30] The results indicated that there are only ionic interactions because of the absence of effective alkali cation···reactant orbital overlap. The effects associated with the substrate–catalyst orbital interaction become important when the reaction is carried out using first-row d-block (Cu(I), Cu(II), Zn(II), Ni(II), Cr(III), Sc(III), and V(V)) cation-exchanged faujasites.[29] The overlap between the d-orbitals of the active site and the MOs of the substrate strongly affect the reaction energetics and can even alter the mechanism of the chemical transformation.

In this work, we studied the variations in ICOHP and BO as a function of the CC bond lengths in 1DAC, TS, and 2DAC per catalyst. This results in ICOHP–BO correlations for every evaluated catalyst model. Selected results of the ICOHP and BO analyses focused on the CC bonds that are plotted versus the interatomic distance in Figure a,b. Note that these plots also include the C1–C6 and C4–C5 interactions present in 1DAC, TS, and 2DAC. Thus, the resulting ICOHP and BO values for the C1–C6 and C4–C5 bonds are correlated to distances longer than the equilibrium CC bond length.

Figure 5

Plots of ICOHP (a) and BO values (b) for CC pair-wise interactions relevant to 1DAC, TS, and 2DAC in various cation-exchanged faujasites and the molecular library. CC pair-wise interactions exceeding 1.6 Å relate to the C1···C6 and C4···C5 interactions in 1DAC and TS. The changes to the ICOHP values for the C2–C3 and C5–C6 bonds during the DAC reaction in MY models are shown in (c).

The results in Figure a,b show that we reproduced the same asymptotic ICOHP and BO curves for the CC bonds in all evaluated models. The trends for the CC bond are thus independent of the chemical composition of the periodic faujasite model. On the one hand, the qualitatively similar CC bond length—ICOHP correlation seems logical as we only measured the CC bonds originating from DMF and C2H4. On the other hand, the chemical surrounding is drastically changed. TMFAU catalysts hold d-block cations with different d-shell fillings and exhibit zeolite matrices with different Al contents. Similarly, MY catalysts hold different alkali cations with different principle quantum numbers for the valence shells and exhibit different degrees of framework basicity.[58] The origin of this behavior requires further investigation.

Additionally, the gas-phase models yield slightly higher ICOHP values for the same bond length as compared to the TMFAU and MY models. Meanwhile, the BO−bond length correlation does not vary upon changing from the gas phase to periodic faujasite models. This is a feature that has to be studied in greater depth also. Nevertheless, we anticipate that the change in the ICOHP versus CC bond distance trend results in a change in ICOHP–BO correlation.

With the ICOHP versus CC bond length trend established to be similar for all evaluated zeolite models in this work, two of the CC bonds in the DMF + C2H4 DAC reactions can reliably be compared. Inspection of Figure c shows that the ICOHP values in 1DAC for the ethylene (C5–C6) bond are ordered as CsY < NaY ≈ LiY < KY < RbY. This is well in line with the ethylene adsorption geometries and the properties of the MY models, discussed elsewhere.[30,47] Briefly, because of the large size of the Cs+ cations, ethylene establishes multiple noncovalent interactions with the accessible sites in the faujasite supercage. This results in polarization of the C5–C6 bond and a reduction of electron density in between the carbon atoms. In LiY and NaY, ethylene is adsorbed on one active site only. The Li and Na cations are both relatively strong Lewis acids as compared to the other alkali cations. In contrast, potassium and rubidium cations are relatively weaker Lewis acids, which polarize the C5–C6 bond less, which results in the highest C5–C6 bond energies found for these systems. However, while 1DAC can be analyzed in a chemical meaningful way, investigations aimed at the other states in Figure c become uncertain because of the relatively small differences between the ICOHP values and the inherent numerical inaccuracies in the employed ab initio method.

The difficulty of discussing relatively small variations of the ICOHP values per state in depth is illustrated by the weak ICOHP–BO correlation. For instance, the C5–C6 bond in 1DAC,RbY has a BO of 2.13 and an ICOHP of −12.87 eV, whereas 1DAC,CsY has a BO of 2.18 and an ICOHP of −12.22 eV. Still, these BOs are significantly higher than those in the TS for which the computed BOs (ICOHP values) are ca. 1.8 (about −10.7 to −11 eV). Continuing along the reaction coordinate toward 2DAC, the C5–C6 ICOHP and BO values decrease to values between −8.12 to −8.38 and 1.08–1.09, respectively. For C2–C3, BOs range from 1.414 to 1.432 in 1DAC and from 1.79 to 1.81 in 2DAC. The respective ICOHPs range from −9.55 to −9.78 eV in 1DAC and from −12.01 to −12.22 in 2DAC. The C2–C3 bond reaches approximately similar values as C5–C6 in state TS.

In summary, the zeolite-based ICOHP interatomic distance trends show significant reproducibility. The formation and cleavage of the CC bonds during the DAC reaction between DMF and ethylene can be studied which involves different models (e.g. 1DAC, TS, 2DAC). Additionally, a reactivity trend on the basis of the catalyst stoichiometry can be established.

Statistical Analysis

The third part of the study was devoted to the investigation of the ICOHP–BO correlation for CC bonds in more chemically complex zeolite-based systems. This additional analysis was required as we observed a change in the gradient of the ICOHP bond length correlation. Plots of the ICOHP–BO correlation in Cu(II)FAU and KY can be found in Figure a,b, respectively. These panels also show the 95% confidence interval and the prediction limit (green and gray lines). The ICOHP–BO correlations of all TMFAU are plotted together in Figure c. The correlations of all MY models can be found in Figure d. Plots showing the linear ICOHP–BO correlation and the associated 95% confidence intervals and prediction limits for every cation individually can be found in the Supporting Information, Figures S1 and S2, respectively. We have also attempted to fit the ICOHP–BO correlations with a polynomial fit. These results can be found in Figures S3 and S4. The resulting studentized residuals obtained by fitting the ICOHP–BO correlation with linear and polynomial plots, are shown in Figure e. The resulting linear fitted parameters are displayed in Table . Parameters for a polynomial fit can be found in Table S14.

Figure 6

ICOHP–BO correlation for CC bonds during the DAC reaction between DMF and C2H4 over cation-exchanged faujasite MY and TMFAU periodic models. In panels (a,b) the green and gray lines demarcate the 95% confidence interval and prediction limit, respectively. Panels (c,d) show all ICOHP–BO correlations obtained for all evaluated cation-exchanged faujasites, without 95% confidence intervals and prediction limits. Panel (e) shows the studentized residuals of the linear and polynomial fits using the data of all TMFAU models.

Table 2

Fitted Parameters for the Linear Fits

ICOHP = a2 × x (BO < 2)
system	a2 (eV/BO)
Cu(I)FAU	–7.3 ± 0.10
Cu(II)FAU	–7.79 ± 0.07
Zn(II)FAU	–7.44 ± 0.08
Ni(II)FAU	–7.53 ± 0.09
Cr(III)FAU	–7.63 ± 0.09
Sc(III)FAU	–7.39 ± 0.09
V(V)FAU	–8.13 ± 0.15
LiY	–7.44 ± 0.16
NaY	–7.29 ± 0.14
KY	–7.28 ± 0.15
RbY	–7.34 ± 0.14
CsY	–7.35 ± 0.15

Both the linear and polynomial fits exhibit R2-adjusted values of ca. 0.98. In Figure e, the large residuals (>3) for data points above a BO of 2 indicate that such points are outliers within the framework of a linear fit. Although the polynomial fit seems to be the best as indicated by a distribution of residuals around a value of zero in panel (e) of Figure , we opted for the linear fit. Such a linear fit allows for a chemical intuitive and chemically relevant interpretation of the ICOHP–BO correlations. The data points related to BO > 2 were, therefore, removed from the datasets to obtain the linear fit. The linear fit was deemed reasonably because (1) the DAC reaction (and many other reactions) does not exceed CC double bonds and (2) the standard errors of the linear fits are small. Furthermore, all parameters of the polynomial fits are characterized by large standard errors. Apart from such correlation being difficult to interpret chemically, these large standard errors make predictions dubious.

The obtained nonlinearity of the ICOHP–BO correlations in cation-exchanged faujasites when compared to the correlations obtained with the gas-phase library is rather surprising. The exact reason for the change in ICOHP–BO correlation is unknown and further research into this matter is required. Yet, we note that the COHP analysis treats chemical bonds as pair-wise interactions with no consideration for the effects of the chemical surrounding on the bond being evaluated.[16] Therefore, we hypothesize that the increased chemical complexity of the cation-exchanged faujasites might be the cause of the nonlinear ICOHP–BO correlation as compared to the molecular gas-phase models involving isolated entities. For instance, donor–acceptor interactions of the exchangeable cations with the substrates will affect the bonds in the substrates and are thus suspected to be at the origin of the nonlinear ICOHP–BO correlation. This hypothesis is further strengthened by the second prerequisite defined after evaluating the ICOHP–BO correlations using the molecular library. This prerequisite involved the inductive effect of functional groups on the CC bonds involving β and γ carbon atoms.

Limiting our discussion to the linear fits, the obtained results indicate that the values for a2 are very similar for the different models. The smallest slope is found for KY with a2 = −7.28 ± 0.15 eV/BO. The largest slope is found for V(V)FAU with a2 = −8.13 ± 0.15 eV/BO. In spirit of our hypothesis, the ICOHP values in V(V)FAU might be significantly affected by the polarizing power of the pentavalent V(V) cation. If this system is omitted from the series, one obtains a maximum slope for Cu(II)FAU with an a2 of −7.79 ± 0.07 eV/BO. The ICOHP–BO correlations in TMFAU models are characterized by a standard error of less than 0.1 eV/BO (V(V)FAU excluded). The trends in MY have standard errors below 0.15 eV/BO. The result is thus an error of 0.2–0.3 eV maximum upon changing the CC BO from 1 to 2. The total ICOHP changes from 7.3 to 7.8 eV. Note that although we constraint the linear fit to intersect the y-axis at the point (0, 0), unconstraint fitting yielded only minor values for the intercept (intercept < 0.5 eV, data not shown). This is markedly different from the gas phase correlations for which significant values for the intercepts were found. The underlying reason for this difference is in need of further research.

Summarizing, the changes to the ICOHP trends with respect to the interatomic distance in TMFAU and MY models as compared to gas-phase systems resulted in a polynomial ICOHP–BO trend instead of a linear one. As a first-order approximation, the ICOHP–BO correlation can be described with a linear fit up to a BO of 2. The resulting fits are practically equivalent. For the models evaluated herein, the fits allow for a description of CC and CO bond strength changes in the DAC reaction along the reaction coordinate both within the same system and across different systems. The practical relevance of this finding is in the fact that it provides us with a new descriptor that can be used directly to compare the reactivity of different catalysts for the same type of reaction as a function of their chemical properties. Such property–activity relationships are of great interest to the catalysis community.

Conclusions

In this work, we investigated the possibility and limitations of correlating DDEC6-derived BOs and COHP-computed bond strengths to quantify chemical bonding in gas-phase molecules and periodic zeolite models. We show that the ICOHP analysis allows obtaining reproducible results when limiting the analysis to one reaction class, albeit in chemical systems with substantially different chemical compositions. Our study implies that the strengths of chemical bonds estimated using the COHP approach can be successfully employed to quantitatively analyze the changes in bonding patterns along chemical conversion routes.

When applied to gas-phase molecules, the ICOHP–BO correlations can be established for different types of chemical bonds and these correlations exhibit a pronounced sensitivity to the stoichiometry of the chemical system and valence shell principle quantum number of the involved atoms. Furthermore, the ICOHP–BO correlations for the same class of substances are affected by the presence of functional groups, reflecting the inductive effects and short-range electrostatic interactions. We identify four key prerequisites for establishing the ICOHP–BO correlations, namely: (1) the individual correlations may be constructed for molecules with similar chemical composition; (2) ICOHP–BO correlations of CC bonds involving β and γ carbon atoms can only be constructed when the number of functional groups or the stoichiometry within the functional group remains similar; (3) for molecules containing CHet bonds with Het being an heteroelement down a group (i.e., different valence principle quantum numbers), an ICOHP–BO correlation may only be established for the CHet bonds, but not for the other bonds in the molecules; (4) only for substantially large variations in the BO/ICOHP values, their direct comparison among the different molecules is possible. The exact threshold depends on the particular system investigated and the accuracy of the correlation (ΔBO ≈ 0.2–0.25).

Despite these promising findings, our study reveals a number of phenomena that require additional theoretical analysis. In particular, it is not clear why for the ICOHP–BO correlations of chemical bonds in gaseous molecules, the extrapolation of BO to zero gives rise to finite values of ICOHP, whereas the ICOHP value vanishes in periodic models. Furthermore, the fundamental origin for the different ICOHP bond distance trends obtained for the periodic and gaseous systems requires further in-depth theoretical analysis beyond the scope of this initial work.

Significantly, the presented results on the ICOHP–BO correlations demonstrate the applicability of the ICOHP analysis method beyond a single structural model. The transferability of the ICOHP parameters is illustrated by considering the DAC of DMF and ethylene catalyzed by faujasite-type zeolite catalysts as a model chemical process. When applied to intermediates belonging to the same reaction class, the ICOHP analysis yields consistent results for different periodic zeolites with varied chemical compositions. Importantly, our study demonstrates the possibility of the direct bond strength quantification using the ICOHP analysis and its applicability to study the formation and cleavage of chemical bonds during the catalytic reactions. These data together with the computed BOs can provide detailed quantifiable bonding information on the reacting chemical systems necessary for constructing quantitative structure–activity relations in complex chemical systems.