Cooperative small molecule activation by apolar and weakly polar bonds through the lens of a suitable computational protocol

Abstract

Small molecule activation processes are central in chemical research and cooperativity is a valuable tool for the fine-tuning of the efficiency of these reactions. In this contribution, we discuss recent and remarkable examples in which activation processes are mediated by bimetallic compounds featuring apolar or weakly polar metal–metal bonds. Relevant experimental breakthroughs are thoroughly analyzed from a computational perspective. We highlight how the rational and non-trivial application of selected computational approaches not only allows rationalization of the observed reactivities but also inferring of general principles applicable to activation processes, such as the breakdown of the structure–reactivity relationship in carbon dioxide activation in a cooperative framework. We finally provide a simple yet unbiased computational protocol to study these reactions, which can support experimental advances aimed at expanding the range of applications of apolar and weakly polar bonds as catalysts for small molecule activation.

---

Diego Sorbelli

Diego Sorbelli received a PhD degree in Theoretical Chemistry and Computational Modeling under the supervision of Prof. Paola Belanzoni and Dr Leonardo Belpassi from the University of Perugia in 2023. His PhD research focused on the modeling of transition metal complexes for applications in catalysis and small molecule activation. After his PhD, he joined the University of Chicago, where he currently works as a postdoc in the group of Prof. Giulia Galli. His research now focuses on the computational modeling and characterization of organic materials for applications in electronics.

Leonardo Belpassi

Leonardo Belpassi is a Senior researcher at National Research Council CNR, Institute of Chemical Science and Technologies [SCITEC-CNR] c/o University of Perugia (Italy). His primary research activities are the study of relativistic effects in the chemistry of systems containing heavy elements; invention and derivation of theoretical ideas in software; and application to a wide range of chemical problems, including weak interactions and gold based catalysis. He is an author of the relativistic Dirac–Kohn–Sham density functional theory module of the BERTHA code.

Paola Belanzoni

Paola Belanzoni studied Chemistry at the University of Perugia (Italy), where she obtained her PhD. She spent several periods at the Vrije Universiteit in Amsterdam, working in collaboration with Prof. Evert Jan Baerends, within the Theoretical Inorganic Chemistry field. She is an Associate Professor at the Department of Chemistry, Biology and Biotechnology of the University of Perugia since 2007. Her current research interests mainly concern homogeneous gold catalysis, organometallic compounds and reactivity, bio-inspired inorganic catalysts, spin states in catalysis, heavy-element containing compound properties, and small molecule activation processes from the computational perspective.

---

Introduction

Activation of inert and abundant small molecules is a central topic in chemical research. These processes are critical for a variety of pressing issues, such as production of fuels, energy storage and conversion, removal of environmentally hazardous and toxic gases from the atmosphere and their use as synthetic building blocks.¹

Notably, the design and implementation of efficient and sustainable strategies to target small molecule activation is challenging, particularly because these substrates can be relatively inert. This has engaged researchers for decades: a multitude of successful approaches have been developed over the years, resulting in a myriad of strategies to activate small molecules both in heterogeneous and homogeneous phase. The latter processes, especially, are highly desirable since they usually enable high selectivity towards these substrates, as well as allowing the implementation of flexible design rules for improving the activity of the species. In this context, well-established and successful approaches encompass (among others) frustrated Lewis pairs,^2,3 main group ambiphiles,⁴ low-valent main group species,^5–7 transition metal complexes with redox-active ligands⁸ and homo- or heterobimetallic compounds.^9–11 Particularly, for CO₂ capture and conversion, attempts by photochemical activation in catalytic transformations deserve to be mentioned.^12–15

Bimetallic species have received a lot of attention in this scenario since the introduction of a second metal moiety implements an additional tunable parameter that makes the design of these compounds even more flexible. In fact, cooperative reactivity in metal–metal bonds can be also finely modulated by tuning the degree of cooperativity between the two centres (single-site vs. multi-site reactivity) and the metal–metal bond order and polarity.^9,10

When dealing with such activation processes, there is an intricate relationship between the structure and the electronic properties of the bimetallic compound and their matching with the electronic properties of the small molecule substrate. In this context, a synergic effort between experiments and theoretical calculations is essential to sketch such a complex picture. Computational endeavours aim at providing not only an accurate description of the possibly unusual and intricate reaction mechanisms that occur in this environment, but also an unbiased description of the electronic properties of the species involved – both the small molecule and the active species – so that they can reinforce and drive experimental discovery.

In this Contribution we report relevant recent breakthroughs in the cooperative activation processes mediated by either weakly polar or apolar bonds in electron rich bimetallic compounds. By discussing the activation of relevant small molecules separately (carbon dioxide, dihydrogen, nitrous oxide, alkynes, water, ammonia), we illustrate the main experimental advances together with detailed theoretical insights, highlighting findings that subvert previously established paradigms both in small molecule activation and in traditional catalytic processes mediated by a single metal moiety. Ultimately, we provide, based on our recent findings, a combined computational protocol that should help the readership navigate through the intricate realm of cooperative activation processes by relying on the ad hoc application of popular and available electronic structure analysis tools.

Carbon dioxide activation

The activation of carbon dioxide is a central and urgent matter in chemical research. CO₂ is, indeed, on one hand a hazardous greenhouse gas (its atmospheric concentration increased dramatically in the last few years over 400 ppm).¹⁶ On the other hand, it is, from a thermodynamical standpoint, an extremely stable molecule ( ¹⁷), thus its activation, capture and eventual use as a C1 building block are mostly hindered by its inertness.

Several cooperative CO₂ activation processes have been reported. Most of these reactions feature the presence of homobimetallic compounds in homogeneous phase, where the apolar M–M bond (M = block s, p or d metal) acts to capture and/or reduce carbon dioxide. A couple of relevant and representative examples are given in Scheme 1, consisting in the cooperative reactivity towards CO₂ of two Mg(I) units in a Mg–Mg complex (1, Scheme 1), that results in the high yield formation of the Mg(II) carbonate complex 2,¹⁸ and the reductive oxygenation of CO₂ by a Ti–Ti double bond in the dititanium complex 3.¹⁹

Scheme 1 Relevant examples of cooperative CO₂ activation with s- (1), d- (3) and p-block (6) homobinuclear compounds featuring apolar E–E bonds. (R = SiⁱPr₃; Ar* = C₆H₂{C(H)Ph₂}₂Me-2,6,4).

These examples are instructive since, in both cases, computational studies have followed experimental findings to elucidate the central role of the M–M (M = Mg, Ti) σ non-polar bonds serving as an electron reservoir populating the LUMO of carbon dioxide and driving its activation.^20,21

Further evidence in this direction comes from the extensive experimental and computational insights into the activation of carbon dioxide mediated by apolar bonds that has been obtained upon the study of CO₂ capture and reduction by p-block cooperation in Group 14 dimetallenes and dimetallynes featuring apolar E–E (E = Si, Ge, Sn) bonds.²² One of the first examples is the reactivity of a sterically hindered digermyne 6 (Scheme 1), which has been reported to insert CO₂ into the Ge–Ge bond (intermediate 7) and to further reduce it to CO at temperatures as low as −40 °C (8).²³

In the pivotal work by Frenking and coworkers, the mechanism of this reduction has been thoroughly explored using spectroscopic and computational techniques. Subsequently, various examples of the reactivity of CO₂ with apolar bonds in the field of homodinuclear main-group element species have been reported. For instance, pincer-ligand supported dimetallynes [ArNEt₂E]₂ (E = Ge, Sn, ArNEt₂ = C₆H₃(CH₂NEt₂)₂-2,6) have been shown to reversibly react with CO₂ under mild conditions to generate insertion compounds of the type [ArNEt₂E(CO₂)EArNEt₂].²⁴ Analogous to digermynes, the products are not stable at ambient temperature and, after uptaking of a second CO₂ molecule, they convert into the corresponding carbonate compounds with accompanying generation of CO.

Only very recently, a few examples of unconventional heterobimetallic compounds featuring weakly polar M–M′ bonds, i.e. metal–aluminyls, have been reported to insert CO₂ into their M–Al bonds (see 9–14 in Scheme 2).

Scheme 2 Experimentally reported examples of complexes derived from CO₂ insertion into Al–M (M = Cu, Ag, Au, Zn, Fe) bonds.

In 2019, Hicks and coworkers reported the CO₂ insertion into an Al–Au bond in the [^tBu₃PAuAl(NON)] (NON = 4,5-bis(2,6-diisopropylanilido)-2,7-di-tert-butyl-9,9-dimethylxanthene) complex (11 in Scheme 2).²⁵ Subsequently, in 2021, Aldridge and coworkers disclosed analogous reactivity of the Al–Cu and Al–Ag bonds in the [^tBu₃PMAl(NON)] (M = Cu, Ag) analogues (9 and 10 in Scheme 2).²⁶ A closely related set of Al–Cu bond reactions with CO₂ were reported by Liu and coworkers (12 in Scheme 2).²⁷ In 2021, Roy and coworkers reported CO₂ insertion into an Al–Zn bond (13 in Scheme 2).²⁸ Finally, in 2022, Mankad and coworkers reported CO₂ insertion into an Al–Fe bond by a [L^dipp(Me)AlFeCp(CO)₂] (L^dipp = HC{(CMe)(2,6-i-Pr₂C₆H₃N)}₂, Cp = η⁵-C₅H₅) complex (14 in Scheme 2).²⁹

In these examples, the CO₂ insertion into M–Al bonds occurs in mild conditions, always yielding an insertion product where the CO₂ carbon atom is coordinated to M and the two oxygen atoms are bonded to the aluminyl moiety. For coinage metal aluminyl complexes (9, 10, 11 in Scheme 2), depending on the Group 11 metal nature, the product of carbon dioxide insertion is stable to further reaction for the [^tBu₃PAuAl(NON)] complex, whereas for the silver analogue [^tBu₃PAgAl(NON)] the insertion product leads to the corresponding carbonate complex (and CO), and for the copper system [^tBu₃PCuAl(NON)] it proceeds rapidly to the carbonate and CO even at low temperatures.

Assessing the polarity of these M–Al bonds is not as trivial as in the previously discussed E–E homobinuclear compounds for obvious reasons. Notably, however, some experimental evidence here suggests that we are dealing with low-polar M–Al bonds. For instance, compounds 9–11 exhibit a very similar reactivity despite the different electronegativities of Cu, Ag and Au, while the formation of 13via a covalent Zn–Al bond was discussed in close relationship with its ionic Li–Al counterpart. Furthermore, a thorough mechanistic investigation via EPR spectroscopy for the formation of 14 highlighted that CO₂ is activated by initial formation of Al^II and Fe^I metalloradicals, as a result of the homolytic dissociation of the weak low-polar Al–Fe bond.

Further characterization of the features of these peculiar M–Al bonds came after our recent computational exploration of the mechanistic pathway for CO₂ insertion into the Au–Al bond in [^tBu₃PAuAl(NON)], where we found striking analogies between the Au–Al and Ge–Ge reactivity for the reduction of CO₂ to CO by [LGe-GeL] (see Fig. 1).³⁰

Fig. 1 Sketched structures of the stationary points along the reaction profile calculated for the CO₂ insertion into the Ge–Ge bond leading to 7 (left)²³ and into the Au–Al bond leading to 11 (right).³⁰ TSI = transition state for step I, INT = intermediate, TSII = transition state for step II, PC = product complex.

For both [LGe–GeL] and [^tBu₃PAuAl(NON)], the reaction mechanism is a two-step process, where the nucleophilic attack to the CO₂ carbon atom (rate-determining step) has a relatively low free energy activation barrier (well below 20 kcal mol⁻¹), a stable intermediate is formed and a second transition state is located with a lower activation barrier associated with the CO₂ rearrangement leading to the thermodynamically stable insertion product.

Such analogies hold even stronger in terms of electronic structure of these compounds. In the case of digermynes, it is often accepted that their high reactivity can be attributed to their nonnegligible biradical character.^31–33 Furthermore, computational studies on their reactivity suggest a central role of their electron-sharing Ge–Ge bond acting as a nucleophilic site, supported by the presence of empty p orbitals on the Ge atoms acting as electrophiles.³⁴

To shed light into the nature of the [^tBu₃PAuAl(NON)]–CO₂ interaction, we used the energy decomposition analysis (EDA)³⁵ and natural orbitals for chemical valence (NOCV)^36,37 methods, revealing striking analogies with digermynes, as shown in Fig. 2.

Fig. 2 Left: Isodensity surfaces associated with the Δρ₁ and Δρ₂ NOCV deformation densities for the transition state structure TSI. The charge flux is red-to-blue. The isodensity value is 5 m_e a₀⁻³ for both surfaces. Right: Most significant contributions from the [^tBu₃PAuAl(NON)] and CO₂ fragments molecular orbitals (MOs) to the interactions described by the NOCV deformation densities. Isodensity values for all MO isosurfaces is 30 m_e a₀⁻³. All isosurfaces are taken and adapted with permission from ref. 30 Copyright 2021, American Chemical Society.

For [^tBu₃PAuAl(NON)] complex, we have in fact demonstrated that the nucleophilic behaviour of the electron-sharing, weakly polarized Au–Al bond combined with Al acting as electrophile is the driving force for the CO₂ insertion reactivity occurring via a cooperative radical-like mechanism. To help the readership to understand the concept of “radical-like mechanism”: with this term we refer to a singlet closed-shell pathway, along which, upon homolitic dissociation of the electron-sharing bond, the substrate insertion (and/or dissociation) is stabilized by the two resulting fragments behaving as radicals would. The main interaction in the CO₂ insertion process has been shown by us to be electron donation from the Au–Al σ bond towards the CO₂ LUMO, assisted by a secondary interaction where electron donation occurs from the CO₂ HOMO towards the Al centre (mostly an empty 3p_z orbital). Notably, a later study employing high-level multireference ab initio calculations further confirmed the low polarity and high covalency of the Au–Al bond as a driving force for the reactivity.³⁸

In the same context, we have also analysed the Cu–Al bond in [^tBu₃PCuAl(NON)] and Ag–Al bond in [^tBu₃PAgAl(NON)], which have been found to have a similar electron-sharing nature, explaining the qualitatively analogous mechanism, in agreement with the experimental features for the CO₂ activation with M–Al (M = Cu, Ag, Au) compounds.

These metal–aluminyl heterobimetallic complexes, in analogy with main-group homonuclear E–E compounds, do not present intrinsic polarization at the M–Al bond, although exhibiting cooperative reactivity towards CO₂ (and other small molecules) in a way that may resemble heterobimetallic polarized M–M′ bonds. These compounds thus represent a unique and particularly interesting playground to learn more about how to assess bond polarization and its effects on the reactivity of heterobimetallic compounds. The simplest way to evaluate bond polarity in such compounds would be resorting in descriptors as simple and effective as atomic electronegativity and atomic charges which, however, may be too simple for some systems. If we think of the gold–aluminyl complex, Au has a much higher electronegativity than Al (2.54 vs. 1.61 on the Pauling scale), thus suggesting a strong Au(δ−)–Al(δ+) polarization. However, some inconsistencies become immediately clear from Fig. 3a: Cu and Ag, for instance, are much less electronegative than Au but a practically identical reactivity is observed for the corresponding Cu/Ag–Al complexes; the comparison between Zn and Al electronegativities (1.65 vs. 1.61, respectively) suggests no bond polarization, yet an analogous reactivity of Zn–Al compounds has been reported experimentally (and similarly for Fe–Al). Likewise, the computation of QTAIM charges on the gold–aluminyl complex reveals an exceptional −0.83 e charge on Au, which would fully support the Au(δ−)–Al(δ+) bond polarization. Yet, as shown in Fig. 3b, the charge on Au spans an extremely wide range of values depending on the approach used for their computation (−0.83/+ 0.22 e), because of, as we demonstrated, the diffuse nature of the aluminyl anion that makes the atomic charge assessment particularly challenging in this context.

Fig. 3 (a) Pauling electronegativities for relevant elements of the periodic table. (b) Atomic charge on the Au atom for the [^tBu₃PAuAl(NON)] complex calculated with different approaches (data taken from ref. 30).

The computational design of metal–aluminyl compounds for carbon dioxide insertion further assesses and confirms the central role of the weak bond polarity in this type of reactivity. For instance, switching from Au–Al to Au–Ga and Au–In bonds, we expect, according to Fig. 3a, no relevant changes in the reactivity: there is no reason not to expect Au(δ−)–Ga(δ+) and Au(δ−)–In(δ+) polarization based on atomic electronegativities (1.81 for Ga and 1.78 for In). We demonstrated, however, that a kinetically and thermodynamically favourable reactivity with CO₂ should only be expected for the gold–aluminyl complex.³⁹ Upon further analysis with several computational approaches, we demonstrated that the Au–Al bond features the most (non-polar) electron-sharing character among the Group 13 anions, which is, in turn, responsible for its highest efficiency and analogous conclusions can be drawn for boryl compounds as well.⁴⁰

Another evidence for the uniqueness of metal–aluminyl compounds (and, again, for the importance of critical bond polarity assessment) comes from a comparison between aluminyl and carbene coordination bonds. It has been extensively reported that there is a strong isolobal analogy between aluminyls and carbenes, both featuring (i) a lone pair on the Al/C atom, respectively, associated with their donating ability, which is responsible for their nucleophilic behaviour and (ii) an empty p orbital (3p for Al and 2p for C), associated with their acceptor capability, which is, in turn, responsible for their electrophilicity (see Scheme 3a).⁴¹

Scheme 3 Schematic representation of the (a) isolobal analogy between singlet carbenes and aluminyls: lone pair and vacant 2p/3p orbitals on the X = C and X = Al⁻ sites of carbenes and aluminyls, respectively; (b) electron-sharing bond in a gold–aluminyl complex; (c) Dewar–Chatt–Duncanson bonding model in a gold–carbene analogue (dative bond).

We quantitatively compared the two ligands in terms of bonding with gold and reactivity towards CO₂, resulting in two clear different trends. While all aluminyls feature non-tunable Au–Al bonds and a similarly exergonic and kinetically accessible CO₂ insertion, Au–C bonds are of dative nature and tunable via carbene structural modification and, therefore, no insertion reactivity can be expected thermodynamically.⁴²

These facts highlight the first important take home message of this Feature Article. In these systems the weak (or absent) bond polarity is the key ingredient for an efficient and feasible CO₂ activation process and thus the unbiased assessment of bond polarization is critical for their understanding, as will be more extensively discussed from a computational perspective in the last section of this article.

In this framework, the shortcomings of other paradigms related to carbon dioxide activation emerge. For instance, the feasibility and efficiency of CO₂ activation is conventionally described in terms of a merely structural and/or electrostatic perspective i.e. (i) bending of the OCO angle from 180°, (ii) elongation of at least one of the two CO bonds, and (iii) polarization of the charge on C and O, including the charge transfer to or from the outside. The first type of CO₂ activation through the bending of the OCO angle is what is more commonly considered as a good indicator of the nature and extent of the metal–CO₂ interaction.⁴³

Scheme 4 displays a partial Walsh diagram for CO₂, which clearly shows how the molecule can be activated by its geometrical (bending) modification. This schematic diagram depicts the energy levels of the most important CO₂ frontier molecular orbitals as a function of the OCO angle. The most stable state is the linear shape. As the angle decreases, the frontier orbitals change their energy levels, followed by the orbital topological mixture of each other. It is known that when the angle becomes close to 90°, the HOMO and LUMO orbital energies become almost the same value, i.e. the molecule is so activated that theoretically electrons can be accepted and released without barrier. The degree of activation of CO₂ is thus often monitored by following its bending distortion along the reaction path. The interpretative framework lies exactly on Walsh diagram, which shows that, upon reduction, gas-phase CO₂ accepts electron charge in its LUMO, which is of antibonding (π*) character and becomes energetically more favoured in a bent structure.

Scheme 4 Qualitative Walsh diagram for the OCO angle bending from 180° to 120°. Most important frontier molecular orbitals of CO₂ are depicted.

This simplified interpretative scheme unfortunately neglects one very important factor involved in the activation process, that is the substrate mode of activation. In the presence of apolar or weakly polar M–M′ bonds, the cooperative activation equally involves two centres (M and M′) through the electron-sharing M–M′ bond that serves as a nucleophilic site for the reaction, possibly leading to non-trivial topologies of the potential energy surfaces (PESs). This is exactly the case for [^tBu₃PAuAl(NON)], as shown in Fig. 4.

Fig. 4 Potential energy surface (PES) in the region surrounding transition state TSI for complex [^tBu₃PAuAl(NON)]. Insets: Position on the PES and schematic structure of TSI and TSI′. Energy has been shifted in each case according to the minimum energy structure. Reproduced with permission from ref. 39 Copyright 2022, American Chemical Society.

In a work by us, we explored the potential energy surface around TSI (see Fig. 1), revealing that the Au–Al bimetallic CO₂ activation occurs in a PES which is extremely flat in a wide range of interfragment distances, indicating that similar energies correspond to very different structures.³⁹ Remarkable changes in the structural parameters characterizing the coordination of carbon dioxide to the complex (i.e. Au–C and Al–O distances, as well as the CO₂ bending angle) correspond to almost negligible changes in their energies (all structures in Fig. 4 are located within a 3 kcal mol⁻¹ range). The PES is so flat (Fig. 4) that two different transition states (TSI and TSI′) could be located along the same reaction coordinate. The two TSs have very different structures: TSI features shorter Au–C (2.403 Å) and longer Al–O (2.569 Å) contacts with respect to TSI′ (2.818 and 2.314 Å, respectively). Most importantly, TSI’ features a much larger OCO bending angle (159.5°) with respect to TSI (146.3°), which, according to discussion of Scheme 4, would suggest a more efficient activation pathway through TSI′. Instead, the two transition states lie extremely close in energy (TSI’ is 0.9 kcal mol⁻¹ higher than TSI in terms of Gibbs’ free energy), meaning that CO₂ activation is equally efficient from a kinetic standpoint along both pathways, despite the remarkable structural differences. Thorough exploration reveals that, along the PES, the cooperative nature of the process induced by the weakly polar nucleophilic Au–Al acts so that the orbital interactions between [^tBu₃PAuAl(NON)] and CO₂ always efficiently counterbalance the variable and increasing distortion penalty of CO₂, leading to an efficient activation along multiple directions.

These findings reveal how, generally, in the framework of cooperative activation processes, using simple structural descriptors to monitor and/or design new activation processes may lead to an erroneous interpretation of the results (notably similar conclusions on the OCO bonding angle have been drawn recently for activation processes mediated by semiconductor oxides⁴⁴).

Dihydrogen activation

Dihydrogen is a promising alternative fuel source due to its potential to provide energy with minimal environmental impact.⁴⁵ One of the drawbacks for its large-scale use as a fuel is its difficult storage and transportation. Thus, the possibility of relying on cooperative activation processes to reversibly capture H₂ is surely attractive, but its high kinetic and thermodynamic inertness make such processes particularly complicated to implement.⁴⁶

Several examples of cooperative activation mediated by apolar E–E bonds have been reported and, notably, very early examples of such reactivity involve aryl- and amido-digermynes analogous to 6. In the pivotal work by Power and coworkers, Ar′GeGeAr′ (Ar′ = [C₆H₃-2,6(C₆H₃-2,6ⁱPr₂)₂]) has been reported to easily activate H₂ under mild experimental conditions (1 atm, 25 °C) in solution (hexane) delivering H–H bond dissociation on the apolar Ge–Ge bond.⁴⁷ Following this pioneering contribution, Jones and coworkers reported on a dihydrogen facile activation by a RGeGeR compound (R = [N(SiMe₃)(C₆H₂Me{C(H)Ph₂}₂-4,2,6)]), which occurs both in solution (toluene or d₆-benzene, 1 atm, 20 °C) and the solid state at temperatures as low as −10 °C.⁴⁸ Later on, the amido digermyne RGeGeR (R = [N(SiⁱPr₃)(2,6-[C(H)Ph₂]2-4-ⁱPrC₆H₂)]) has also been shown to easily activate H₂ in solution in mild conditions (toluene, 20 °C).⁴⁹ For all compounds, experiments and theoretical investigations outline a reaction mechanism that proceeds via a singly-bridged intermediate [RGe(μ-H)GeHR] species (see 15 in Scheme 5), which subsequently, upon isomerization, yields different hydrogenation products, experimentally revealed depending on the steric hindrance of the substituents.³⁴ Since then, several examples of analogous reactivity have been reported for Group 14 dimetallenes and dimetallynes, performing H₂ activation mediated by apolar Si–Si and Sn–Sn bonds (see 16 and 17 in Scheme 5 for two relevant examples).¹⁸

Scheme 5 (a) Experimentally reported examples of complexes formed from H₂ insertion into E–E (E = Ge, Sn, Si) bonds; (b) computationally predicted complexes formed from H₂ insertion into Al–M (M = Cu, Ag, Au) bonds.

Computational studies on these compounds and their reactivity have shed light on the features of the dihydrogen activation. EDA-NOCV-based analysis by Frenking and coworkers highlighted that the reactivity towards H₂ is induced by the population of the antibonding σ* LUMO of H₂ by the Ge–Ge centred HOMO of the digermyne, supported by population of the LUMO of the digermyne by the bonding σ HOMO of H₂. The reaction mechanism has been described as a two-step H atom abstraction from H₂ followed by radical pair recombination induced by the substantial diradical character of digermynes.³⁴

Inspired by this experimental evidence and by the strong analogy between the apolar Ge–Ge and weakly polar M–Al (M = Cu, Ag, Au) bonds, very recently we reported that analogous H₂ activation by [^tBu₃PMAl(NON)] complexes may be in principle feasible and ideal for the use of [^tBu₃PAuAl(NON)] as hydrogenation catalyst.⁵⁰ In our computational work, we predicted that [^tBu₃PMAl(NON)] (M = Cu, Ag, Au) complexes should indeed react with dihydrogen, with a very similar mechanism to that computed by Frenking and coworkers for digermynes. The first hydrogenation was reported to occur in a single step which leads to the formation of the singly bridged [^tBu₃PM(μ-H)AlH(NON)] species (18–20, Scheme 5), with a relatively low energy barrier. Calculations have predicted the experimentally accessible formation of a singly-bridged [^tBu₃PM(μ-H)Al(H)(NON)] species for all the metal–aluminyl complexes and an additional and preferred doubly-bridged [^tBu₃PM(μ-H)₂Al(NON)] product for copper and silver (21–22, Scheme 5).

The striking resemblance with digermynes appears even more evident in terms of electronic structure. We analyzed the nature of the complex–H₂ interaction along the reaction path, revealing that it consists in (i) the depletion from the Au–Al σ bonding orbital (HOMO) to the empty σ* MO of H₂ and (ii) the population of the LUMO of the complex (mainly of Al 3p_z character) from the HOMO of H₂, consistent with digermyne–H₂ interactions as described earlier. Notably, interactions between [^tBu₃PMAl(NON)] (M = Cu, Ag) at the corresponding transition states are of analogous nature, with NOCV results showing analogous depletion/accumulation patterns.

These analogies lead to two notable implications. The first is of a more practical nature and concerns possible usage of these compounds as catalysts for a combined H₂ activation and subsequent substrate hydrogenation. The role of gold molecular complexes in catalysis is mainly known in the π-activation of unsaturated substrates,⁵¹ and therefore this catalytic application would be clearly novel. In particular, the H₂ cleavage with gold–aluminyl complex has been calculated to be thermoneutral (ΔG = 3.3 kcal mol⁻¹): this is ideal for its use as catalyst in hydrogenation reactions. In our work, by selecting ethylene as a substrate, we have modelled the alkene hydrogenation process catalysed by [^tBu₃PAuAl(NON)], which has been found to occur in a single, exergonic step (ΔG = −25.6 kcal mol⁻¹), via a concerted transition state, where both hydrogens are simultaneously transferred to the substrate, forming ethane and regenerating the gold–aluminyl catalyst, with a kinetically accessible barrier (ΔG^# = 21.7 kcal mol⁻¹).

Secondly, this computational exploration allows to further characterize this type of reactivity mediated by weakly polar bonds compared to polarized dative bonds. Bourissou and coworkers experimentally and computationally characterized the reactivity of a Pt(0)–Al(III) compound featuring a dative Pt → Al bond towards H₂ (Fig. 5a).⁵²

Fig. 5 (a) Experimentally reported reactivity of a Pt → Al complex with H₂.⁵² (b) and (c) Isosurfaces of the HOMO of the gold–aluminyl complex (b) and platinum–aluminium complex (c). The isovalue is 3 m_e a₀⁻³. The isosurfaces have been taken and adapted with permission from ref. 50 Copyright 2023, Royal Society of Chemistry.

In their work, a two-step mechanism for the formation of the product, which strongly resembles complex 20 (Fig. 5), with a bridging hydrogen between Al and Pt and a hydride on Pt, is reported. These analogies prompted us to compare the complex–H₂ interaction in the two complexes, revealing that the apparent similarities arise from substantial differences in the electronic structure of the two complexes.

At a first sight, the complex–H₂ interaction is analogous: we find it to be mainly driven by the population of the σ* antibonding LUMO of H₂ by the bonding HOMO of the Pt–Al complex. However, at a qualitative level, the nature of the MO is different, being mainly located on the Pt moiety of the complex in Al–Pt and almost equally shared between the two metals in Al–Au (Fig. 5b and c). The atomic orbital contributions to the MO corroborate this picture, showing a mostly Pt-centered MO (72.2%) with negligible contributions from Al (3.0%), consistent with a dative bond and a Pt(δ−) polarization. In the case of the gold–aluminyl complex, instead, the HOMO is centred on both Au and Al, contributing almost equally (19.5% and 15.9%, respectively), consistent with the electron-sharing bond described above. Consequently, while the former reactivity can be formally described as an oxidative activation at the Pt(0) centre, the latter is most likely a radical pair stabilization reactivity, similar to digermynes. An unusual picture of the reactivity with heterobimetallic compounds clearly emerges from these results, that can be exploited to discover new routes for cooperative homogeneous activation processes.

Nitrous oxide activation

Nitrous oxide (N₂O) is a greenhouse gas with a global warming potential of more than 300 CO₂ equivalents,⁵³ with extensive anthropogenic emissions mainly due to agricultural activities.^53,54 It is currently the main responsible for ozone depletion,^55,56 so clearly efficient strategies for its activation are highly desirable, also because it is a valuable synthetic tool behaving as an oxygen-atom transfer agent, with release of non-hazardous N₂ as by-product.^57,58 The main drawback regarding its activation and reduction to N₂ (which is highly exergonic) is that this process is kinetically hindered, with an associated activation barrier of more than 60 kcal mol⁻¹.^59–61

Notably, kinetically favourable N₂O reduction ubiquitously occurs in nature mediated by metalloenzymes equipped with multiple metal centres acting cooperatively, such as the four copper sites in nitrous oxide reductase (N₂OR).^59,62 Inspired by them, several biomimetic and multicopper compounds, as well as homo- and heterobimetallic (and -trimetallic) species, have been reported to efficiently reduce N₂O cooperatively¹¹ and notably, efficient reduction mediated by Group 14 dimetallenes and dimetallynes has been extensively explored,²² as shown in Scheme 6.

Scheme 6 (a) Experimentally reported examples of complexes derived from N₂O activation with E–E (E = Ge, Si) bonds; (b) computationally predicted complex derived from N₂O activation with the Al–Au bond.

In a very recent computational work by us, we reported that analogous N₂O activation with weakly polar, electron-sharing Au–Al bond in [^tBu₃PAuAl(NON)] complex may be feasible, resulting in the formation of the [^tBu₃PAuOAl(NON)] oxo-species (see 25 in Scheme 6) and N₂ extrusion.⁶³ This reaction has been predicted to be accessible from a kinetic point of view in mild conditions (Gibbs’ free energy activation barrier of 15.2 kcal mol⁻¹). The reaction mechanism is qualitatively different from those illustrated so far for coinage metal–aluminyl reactions with CO₂ and H₂. In this case, the reaction is reported to occur through a one-step mechanism with a concerted transition state, where N₂O coordinates to the two active Au and Al sites with a μ-1,2-O,N binding mode.

A thorough analysis of the electronic structure and interactions occurring along this reaction pathway has been carried out within a similar approach as previously discussed for the other substrates. The analysis revealed, once again, that the main driving force of such favourable predicted activation lies in the electron-sharing Au–Al bond serving as a nucleophilic site for the reaction and populating the π* LUMO of N₂O, assisted by the population of the 3p_z orbital of Al by the terminal 1,O-coordination of N₂O. Furthermore, computational catalytic modelling indicates that the resulting Au–O–Al complex may serve in this context as a catalyst for several oxidation processes, including the challenging oxidation of phosphines.

Interestingly, the μ-1,2-O,N coordination mode we reported for the activation is highly unusual. Generally, there has been an effort in the literature to try to understand the various N₂O binding modes occurring at bi- or multimetallic centres.¹¹ For instance, the mechanism of action of nitrous oxide reductase has been studied thoroughly, suggesting that the reduction proceeds with an η²-1,3-N₂O binding at the two active copper centres (Scheme 7a).⁶⁴

Scheme 7 (a) Proposed N₂O binding mode at the two copper centres of the nitrous oxide reductase enzyme. (b) and (c) Observed N₂O binding modes over copper centres in biomimetic dicopper compounds.

The same reactivity with a biomimetic multicopper assembly has been proposed to proceed via a μ-1,O binding (Scheme 7b),⁶⁵ while the reactivity of a biomimetic mixed-valence dicopper complex has been proposed to proceed via an η¹-1,O coordination (Scheme 7c).⁶⁶ The μ-1,2-O,N binding mode at the transition state unveiled by our calculations appears to be really uncommon in this framework. For this reason, in analogy with carbon dioxide activation, we explored the PES around the concerted transition state. As shown in Fig. 6 (and in analogy with carbon dioxide) the topology of the PES is particularly flat around this transition state.

Fig. 6 Potential energy surface (PES) in the region neighbouring the transition state TS for the reaction of complex [^tBu₃PAuAl(NON)] with N₂O. Insets: Position on the PES and schematic structure of TS and TSI*. Energy has been shifted in each case according to the minimum energy structure. Reproduced with permission from ref. 63 Copyright 2023, Wiley-VCH GmbH.

Our calculations show indeed a flat topology of the PES in a wide range of Au–N and Al–O distances around the transition state, leading eventually to locate an alternative almost isoenergetic transition state along the same reaction coordinate where, however, nitrous oxide is coordinated to the complex with a more common μ-1,O binding mode. Notably, electronic structure analyses unravel that, despite the structural difference, the main interaction between nitrous oxide and the complex remains the charge transfer from the Au–Al bond towards the N₂O LUMO, supported by the population of the vacant 3p_z orbital of Al by N₂O HOMO.

These findings further highlight that, while inspection of N₂O binding modes has been often subject of study since considered to be revealing of the nature of the complex-N₂O interaction, in the case of such bimetallic cooperative reactivity, the topology of the PES allows different binding modes to coexist towards the same reaction coordinate.

Other small molecule activation processes

As mentioned in the Introduction, one of the advantages of using bimetallic compounds in the framework of small molecule activation processes is that the addition of a second metal moiety modifies the properties of the first, paving the way for unconventional reactivity patterns. In the case of gold, for instance, examples of reactivity of mononuclear gold complexes towards carbon dioxide are scarce,⁶⁷ in sharp contrast with the facile CO₂ capture by the bimetallic Au/Al complex.

An even more striking example of this paradigm breakdown involves the reactivity of coinage metal compounds with alkynes. Coinage metal (and particularly gold) complexes are widely used in organometallic chemistry as powerful π-acid catalysts, with their soft Lewis acid behaviour enabling the facile activation of alkynes (and generally unsaturated substrates) towards nucleophilic attack.⁵¹

A recent work by Aldridge and coworkers showed how the [^tBu₃PCuAl(NON)] complex reacts with internal alkynes leading to the formation of syn alkyne insertion products under kinetic control, with following formation of the thermodynamic anti insertion product (Scheme 8).⁶⁸

Scheme 8 Experimental insertion of internal 3-hexyne in the Cu–Al bond of the [^tBu₃PCuAl(NON)] complex.⁶⁸

In a following computational work by us, we have both rationalized the experimental results for the copper–aluminyl complex, which more favourably forms the syn insertion product than the anti, and also predicted a similar, although less efficient, reactivity with gold–aluminyl.⁶⁹ The reaction mechanism for 3-hexyne insertion into Al–M (M = Cu, Au) bonds has been reported to be qualitatively very similar to those illustrated for coinage metal–aluminyl reactions with CO₂, H₂ and N₂O. Despite a different geometrical rearrangement of both the complex and alkyne moieties, the analysis of the substrate-complex interaction at relevant transition states along the reaction path has revealed that analogous alkyne-complex interaction occurs. Namely, (i) the electron sharing M–Al bond acts as a nucleophile by donating electron charge to the π* LUMO of the alkyne in a radical-like mechanism and (ii) a second significant component of the interaction is represented by a reverse alkyne HOMO-to-complex (mainly empty Al 3p_z orbital) charge transfer.

While perfectly consistent with the previously reported reactivity of metal–aluminyl complexes, this type of reactivity is in striking contrast with the “classical” reactivity of mononuclear coinage metal complexes towards alkynes. In fact, while the more electrophilic and soft gold centers are more efficient than copper ones in the π-activation of alkynes, in this context we observed an inverted trend, where the copper–aluminyl complex is more efficient towards alkyne insertion, consistent with experimental results. As shown in Scheme 9, this can be rationalized by the remarkable difference in the electronic structure of the bimetallic complexes induced by the aluminyl moiety.

Scheme 9 (a) Electrophilic behaviour of [ML]⁺ (M = Cu, Au) complexes for alkyne activation (Au more efficient than Cu); (b) nucleophilic behaviour of the M–Al (M = Cu, Au) bond in coinage metal–aluminyl complexes towards alkynes (Cu more efficient than Au).

The electronic structure analysis we carried out has highlighted that the difference has a common origin, namely the relativistically stabilized 6s atomic orbital of gold. In fact, it has been extensively discussed – and widely accepted – that the exceptional reactivity of gold towards alkynes is related to the stabilization of the 6s orbital of gold induced by relativistic effects (Scheme 9a).⁵¹ The atomic ns (n = 4, 6 for Cu and Au, respectively) orbital is involved in the formation of the empty [LM]⁺ fragment orbital, and the exceptionally stable 6s orbital of Au makes it more suitable for accepting electron density from a substrate, thus enhancing the electrophilic power of gold [LAu]⁺-type complexes.

This becomes, however, a detrimental feature in the case of the nucleophilic reactivity of M–Al compounds (Scheme 9b). In these compounds, the ns atomic orbital is involved in the M–Al bonding, but in the case of Au, the 6s orbital is less engaged in the bonding with respect to the 4s orbital of Cu, leading to a lowered contribution to the bond from the gold fragment with respect to copper. As a result, the Au–Al bond is less nucleophilic than the Cu–Al bond and, consequently, a less efficient alkyne insertion should be expected.

Such trend inversion in standard reactivity paradigms can be also expected beyond small molecule insertion into M–M′ bonds. One interesting example is the activation of the N–H bond in ammonia. A few examples of N–H bond activation processes mediated by apolar bonds, such as Group 14 dimetallenes and dimetallynes are available.²² On the other hand, examples of such activation mediated by mononuclear transition metal complexes are scarce, due to the strong coordinating ability of NH₃ that prevents the activation process to occur.⁷⁰

Our computational studies on gold–aluminyl compounds revealed that the introduction of the aluminyl moiety and the subsequent low polarity of the resulting Au–Al bond could create an ideal scenario for using transition metal compounds for this scope. We demonstrated that the [^tBu₃PAuAl(NON)] complex could, in principle, activate the N–H bond in ammonia in mild conditions, through a kinetically and thermodynamically favourable one-step process (Scheme 10).⁶³

Scheme 10 Sketched structures of the stationary points along the reaction profile calculated for the O–H/N–H bond activation of H₂O/NH₃ with the Au–Al bond. RC = reactant complex, TS = transition state, PC = product complex.

The N–H bond activation has been found to be exergonic by 9.3 kcal mol⁻¹, with an activation barrier of 13.3 kcal mol⁻¹. We also computed the O–H bond activation process of water finding an even more favourable process, that is exergonic by 20.8 kcal mol⁻¹, with an even smaller activation barrier of 5.5 kcal mol⁻¹. The H₂O- and NH₃-complex interactions have been shown to arise from the same components of the CO₂-, H₂-, N₂O- and alkyne-complex ones, i.e. the electron donation from the Au–Al covalent bond towards the O–H/N–H antibonding molecular orbital is the driving force of the reaction, combined with a sizable electron donation from the lone pair on the O/N atom towards the LUMO of the complex.

All these examples provide interesting evidence for a scenario in which the introduction of a second metal moiety and subsequent formation of a low-polar M–M′ bond can induce changes in established paradigms for the reactivity of mononuclear transition metal compounds. This means that previously inaccessible reactivity of coinage metal mononuclear compound can become feasible, paving a new direction for the molecular design of these compounds towards small molecule activation.

Theoretical tools for cooperative small molecule activation

At this point, we have hopefully convinced the reader that computational insights on cooperative small molecule activation processes are essential to provide a proper rationalization of the observed reactivity, as well as to give valuable insights into the molecular design of multicentre compounds with improved reactivity. However, a black-box use of computational tools in this framework may lead to pitfalls and, eventually, to an erroneous picture of the interactions at play. For this reason, in this last section, we would like to illustrate a composite computational setup to use in these situations that may help the readership in rationally characterizing these reactions.

The steps of the proposed computational protocol are reported in Fig. 7. We divide the scheme in three major blocks, being (1) bond nature/polarity assessment, (2) reaction mechanism and (3) complex–substrate interaction analysis. We discuss below the computational tools which are of invaluable use across each block in more detail, focusing on the reaction of the gold–aluminyl complex [^tBu₃PAuAl(NON)] with carbon dioxide as a case study.

Fig. 7 Schematic workflow for the application of the composite protocol for the theoretical analysis of cooperative small molecule activation processes.

(1) Bond nature/polarity assessment

(a) Comparative energy decomposition analysis (EDA). One crucial step for providing the least biased description of a chemical bond is based on an empirical screening using the widely popular energy decomposition analysis (EDA)³⁵ scheme (1a in Fig. 7).

In general, the EDA scheme aims to provide a simple yet effective description of chemical bonds. For the bond between two generic E and E′ fragments, their interaction energy (ΔE_int) at their in-adduct geometry can be decomposed within the EDA framework as follows:

ΔE_int = ΔE^Pauli + ΔV_elst + ΔE_oi + ΔE_disp (1)

where ΔE^Pauli corresponds to the Pauli repulsion interaction between occupied orbitals on the two fragments, ΔV_elst represents the quasiclassical electrostatic interaction between the unperturbed charge distribution of the fragments at their final positions, ΔE_disp takes into account the dispersion contribution and ΔE_oi is the orbital interaction, which arises from the orbital relaxation and the orbital mixing between the fragments, and accounts for electron pair bonding, charge transfer, and polarization. EDA calculations can be easily carried out with commonly used and available quantum chemistry softwares, such as ADF,⁷¹ ORCA⁷² and BERTHA,⁷³ where the recent relativistic 4-component EDA implementation including spin–orbit coupling is available.⁷⁴ For less experienced users, we note that that ADF also provides a step-by-step tutorial for the setup of an EDA calculation.⁷⁵

One crucial step for the analysis of a generic E–E′ bond is the definition of the separate E and E′ fragments (and the overall complex) in terms of charge and spin multiplicity. While for the overall complex this is straightforward, in the case of the fragments, their definition requires caution. As previously discussed by Frenking, Schwerdtfeger and coworkers, the least biased way to describe the fragments in an EDA calculation is to empirically test the EDA results with different possible fragmentation schemes, i.e. starting either from closed shell [E]^+/− and [E′]^−/+ or open shell doublet [E]˙ and [E′]˙ fragments.⁷⁶ They showed that the most suitable E–E′ bond fragmentation is the one with an associated less stabilizing ΔE_oi value, corresponding to the fragments with an electronic structure which more closely resembles that they acquire upon bond formation.

Taking the previously discussed [^tBu₃PAuAl(NON)] complex as an example, the most obvious choice for describing the formation of the gold-aluminyl bond would be choosing singlet closed shell cationic gold ([^tBu₃PAu]⁺) and anionic aluminyl ([Al(NON)]⁻) fragments. As shown in Table 1, this fragmentation results in a ΔE_oi value of −105.3 kcal mol⁻¹. If we choose the [^tBu₃PAu]⁻ and [Al(NON)]⁺ fragments, as expected, we get a more unlikely fragmentation scheme, with an associated orbital interaction energy of −225.6 kcal mol⁻¹.

Table 1 Comparative EDA results for the [^tBu₃PAuAl(NON)] complex from ref. 30

Fragmentation scheme	ΔEoi (kcal mol^−1)
[^tBu3PAu]^+ + [Al(NON)]^−	−105.3
[^tBu3PAu]^− + [Al(NON)]^+	−225.6
[^tBu3PAu]˙ + [Al(NON)]˙	−71.5

---

Instead, starting from doublet open-shell gold and aluminyl fragments, [^tBu₃PAu]˙ and [Al(NON)]˙, respectively, we get a much less stabilizing ΔE_oi value (−71.5 kcal mol⁻¹), clearly indicating that the Au–Al bond in this complex is best described by open-shell fragments, i.e. as an electron-sharing bond.

(b) Natural orbitals for chemical valence (NOCV) analysis. Once selected the most suitable fragmentation scheme for describing the bond of interest, a qualitative and quantitative assessment of the bond nature and polarity can be performed by coupling the EDA results with the natural orbitals for chemical valence (NOCV)^36,37 and charge displacement (CD)⁷⁷ approaches.

Briefly, the NOCV approach is based on the rearrangement of the electron density occurring when a chemical bond is formed (Δρ′). This deformation density can be brought into diagonal contributions in terms of NOCVs. In the NOCV scheme, the charge rearrangement taking place upon bond formation is obtained from the occupied orbitals of the two fragments suitably orthogonalized to each other and renormalized (promolecule). The resulting electron density rearrangement can be expressed in terms of NOCV pairs which are defined as the eigenfunctions of the so-called “valence operator”^78–80 as follows:

where φ_+k and φ_−k are the NOCV pairs orbitals and ν_±k are the corresponding eigenvalues. When the adduct is formed from the promolecule, a fraction ν_k of electrons is transferred from the φ_−k to the φ_+k orbital, which are envisaged as donor and acceptor orbitals, respectively.

The NOCV approach can be coupled with the EDA decomposition within the ETS-NOCV framework.⁸¹ With this tool, the EDA orbital interaction term ΔE_oi can be further decomposed into NOCV pairwise orbital contributions which associate an energy contribution (ΔE^k_oi) to each NOCV deformation density . NOCV calculations can be easily carried out using ORCA, BERTHA or ADF softwares and we refer the reader to a detailed tutorial on how to setup NOCV calculations.⁷⁵

In practice, following the gold–aluminyl example, after applying the EDA for the Au–Al bond formation starting from [^tBu₃PAu]˙ and [Al(NON)]˙ fragments, we compute the NOCV deformation densities and employ the ETS-NOCV decomposition of the orbital interaction energy term (Fig. 8).

Fig. 8 Isodensity surfaces (1 m_e a₀⁻³) of the two main NOCV deformation densities (top) and (bottom) (charge flux is red → blue) for the Au–Al bond in [^tBu₃PAuAl(NON)] together with the corresponding eigenvalue (|ν_k|) and associated orbital interaction energy contribution (ΔE^k_oi) adapted and reprinted with permission from ref. 30 Copyright 2021, American Chemical Society.

The application of the NOCV approach yields one main component (k = 1) of the Au–Al interaction in the complex (which in this case, since we are starting from open shell fragments, consists of two components, i.e. α and β). At first, a qualitative inspection of the associated isosurfaces is needed, revealing that the two components are of analogous nature (accumulation/depletion regions are mainly localized on Au and Al and have a σ cylindrical symmetry) but describe charge fluxes in opposite directions. is an aluminium-to-gold charge flux, while represents a gold-to-aluminium charge transfer. A first semi-quantitative perspective on the bond nature and polarity can be obtained by inspecting the eigenvalues (|ν_k|, see eqn (2)) and the orbital interaction energy contribution (ΔE^k_oi) associated to each. In this case, similar values are associated to both and (|ν_k| = 0.45 and 0.42 e, respectively), while ΔE^k_oi is −32.7 and −24.5 kcal mol⁻¹, respectively, (see Fig. 8) suggesting that the two opposite charge fluxes may result in a weak bond polarization.

(c) Charge displacement (CD) analysis. A quantitative measure of the bond polarization in these systems can be obtained by coupling the NOCV scheme with the charge displacement (CD) analysis, resulting in the CD-NOCV scheme.⁸² The CD analysis allows to quantify the amount of electronic charge that is transferred between the two fragments upon formation of the E–E′ bond. The charge displacement function (Δq) is defined as the partial progressive integration on a suitable z-axis of the deformation density Δρ′:The CD function, Δq(z), quantifies at each point of the bond axis the exact amount of electron charge that, upon formation of the bond, is transferred from the right to the left across a plane perpendicular to the bond axis through z. In the CD-NOCV scheme, the density rearrangement due to the bond formation between two fragments, (Δρ′), can be partitioned into different NOCV deformation densities and therefore one is able to quantify the charge transfer (CT) associated to each different component. It must be noted that only few of the NOCV pairs contributes to the chemical bond. Therefore, when the CD-NOCV analysis is carried out, usually only the first components are investigated in order to understand which significant chemical contribution to the bond they represent. In eqn (3), the integration axis is usually conveniently chosen as the bond axis between the two fragments constituting the adduct and usually, we choose to evaluate the charge transfer between E and E′ by taking the CD value at the “isodensity boundary”, i.e. the z-point where equally valued isodensity surfaces of the isolated fragments become tangent.^77,83 The application of CD-NOCV scheme is quite straightforward, consisting in the orientation of the bond of interest along the integration axis and subsequent integration of the NOCV components calculated in the previous step. All can be done by printing the NOCV deformation densities on cube files and then proceeding with integration (and/or manipulation) with a set of publicly available tools.⁸⁴

As a result, one obtains CD-NOCV curves, as shown for the gold–aluminyl complex in Fig. 9. The CD-NOCV curves reflect what could be qualitatively inferred from the NOCV isosurfaces. The curve associated to (red curve in Fig. 9) is positive over the whole molecular region, indicating an aluminium-to-gold charge transfer, while the curve (blue curve in Fig. 9) is entirely negative, consistent with an opposite gold-to-aluminium charge transfer. In this case, however, we can associate a charge transfer (CT) value to each component by evaluating its extent at the isodensity boundary. As a result, we find a CT_1α value of 0.30 e and CT_1β value of −0.27 e, confirming that the two charge fluxes are quantitatively close. To definitively assess bond polarization, we can compute the CD curve associated to the overall deformation density (grey curve in Fig. 9), which indicates the overall amount of charge transferred between the two fragments. The curve approaches zero in the bond region and the associated CT value is 0.05 e, definitively confirming that the net polarization of the Au–Al bond is extremely weak, consistent with a poorly polarized electron-sharing bond.

Fig. 9 CD-NOCV curves associated to the first two NOCV deformation densities (, in red and , in blue) for the analysis of the Au–Al bond in the [^tBu₃PAuAl(NON)] complex. The grey curve is associated to the overall deformation density Δρ′. The red dots indicate the position of the nuclei along the z axis. The vertical solid line marks the isodensity boundary between the fragments. Positive (negative) values of the curve indicate right-to-left (left-to-right) charge transfer. Adapted and reprinted with permission from ref. 30 Copyright 2021, American Chemical Society.

2 Reaction mechanism

(a) Structures computation. After bond characterization, the second step of the protocol consists in computationally characterizing the reaction mechanism of the process under study. There are several ways to account for experimental conditions (e.g. implicit solvation, temperature, etc…), allowing to simulate activation processes as much accurately as possible. This step allows to also computationally screen the feasibility of such processes in terms of thermodynamics and kinetics based on Gibbs’ free energies. Computing reaction mechanisms can be non-trivial, since a lot of input from the scientist is required in choosing appropriate reaction coordinates for exploring the reaction pathway. This is especially true in the case of cooperative activation processes, since the not necessarily known interplay of the features of the involved active centres can make the identification of an appropriate reaction coordinate complicated.

Note that there are some automated tools that can be helpful to assess the proper reaction coordinate, but we advise that they should be used with extreme caution in the case of cooperative activation processes mediated by apolar or weekly polar bonds. One of these tools is the intrinsic reaction coordinate (IRC) approach.⁸⁵ The minimum energy reaction path can be in general defined as the steepest-descent path from a transition state to the neighbouring local minima. The IRC approach uses a similar definition, but it relies on mass-weighted coordinates, following the maximum instantaneous acceleration to determine the features of the reaction path.⁸⁶ In practice, the IRC approach uses as input a transition state (or a structure that is reasonably close to it) and determines the energy profile and its length, shape and curvature automatically.

While particularly attractive, the IRC approach can fail in circumstances where the potential energy surface has a complex topology around the transition state. It has been discussed in the literature that a single IRC cannot define a particular preferred reaction pathway when the potential energy surface is flat⁸⁷ and we already encountered the same issue before when studying computationally CO₂ activation processes.⁸⁸

The case of the carbon dioxide activation by a gold–aluminyl complex is a textbook example of this type of situation. We extensively discussed the particularly complex topology of the flat PES in the case of this reaction, and this prevented the possibility for us to use the IRC to help providing insights into the connectivity between transition states and local minima. As shown in Fig. 10, when we used the IRC approach on TSI (see Fig. 1), it resulted in a non-physical shape of the reaction path, basically with a transition state connecting two identical local minima, which is obviously an unphysical result due to the flat topology of the PES.

Fig. 10 One-dimensional IRC minimum energy reaction path for the first step of the reaction of the [^tBu₃PAuAl(NON)] complex with CO₂. Insets: IRC-optimized geometry for the two minima adjacent to TSI. Adapted and reprinted with permission from ref. 40 Copyright 2022, American Chemical Society.

While particularly helpful in the computational exploration of reaction pathways, this example shows that tools like the IRC approach can lead to unphysical results for cooperative activation processes featuring apolar bonds, where a non-black-box PES exploration should always be preferred.

(b) Activation strain model (ASM) analysis. Once the reaction mechanism has been computed, it is useful to characterize the pathway in more detail to grasp more information on the computed energetics. A straightforward yet effective way to do so is by decomposing the relative energy of the calculated stationary points by using the activation strain model (ASM).^89–91

The model has a practically analogous starting point to that of the EDA approach. It allows to decompose the relative energy of a stationary point in two contributions: a penalty arising from the distortion of the reactants from their relaxed geometry to their in-adduct one (ΔE_dist) and a (usually) stabilizing interaction contribution arising from the interaction between reacting fragments (ΔE_int). Thus, the relative energy of a stationary point (ΔE) can be expressed as follows:

ΔE(ζ) = ΔE_dist(ζ) + ΔE_int(ζ) (4)

where ζ represents the reaction coordinate. For a detailed protocol describing the necessary calculations to be carried out within the ASM scheme, we refer the reader to ref. 91.

In the case of the gold-aluminyl complex reacting with CO₂, the ASM approach was particularly useful to understand the factors behind such favourable reaction kinetics. As shown in Fig. 11, we applied ASM to decompose the first part of the reaction pathway, i.e. the first activation barrier, by using eqn (4) to decompose the relative energies of “RC” and TSI” stationary points.

Fig. 11 Activation strain model (ASM) decomposition of the electronic energy activation barrier between RC and TSI for the reaction of [^tBu₃PAuAl(NON)] complex and carbon dioxide. Adapted with permission from ref. 26 Copyright 2021, American Chemical Society.

The results of the ASM analysis revealed a few important features that keep this barrier acceptably low to enable the reaction to proceed as fast at room temperature. ASM results unveil that, while the overall distortion penalty ΔE_dist for this step is high (21.8 kcal mol⁻¹), this is efficiently counterbalanced by the stabilizing interaction between the complex and carbon dioxide (ΔE_int = −12.9 kcal mol⁻¹), resulting in a particularly low activation barrier. Furthermore, decomposition of the distortion penalty revealed that most of it arises from the geometrical distortion (bending and anti-symmetric C–O bond stretching) of CO₂ (20.3 kcal mol⁻¹). This means that the gold–aluminyl complex, despite undergoing a notable geometrical distortion upon approaching TSI, has a very small associated energy penalty (1.5 kcal mol⁻¹), revealing that the complex is particularly flexible, thus favouring the kinetically accessible CO₂ activation process.

3 Complex–substrate interaction analysis

(a) Energy decomposition analysis (EDA). In the previous step, we discussed the positive effect of the complex–CO₂ interactions in stabilizing the activation barrier at TSI. The interaction energy term (ΔE_int) is defined as the interaction energy between the two fragments (CO₂ and complex) at their in-adduct (RC and TSI) geometry, i.e. it has the same definition as the EDA approach illustrated in step 1. This means that we can characterize these interactions by relying on EDA approach.

The application of the EDA analysis is totally analogous to that discussed in step 1a. In this case, however, we have a much simpler situation. At both RC and TSI the two fragments are interacting, but no actual bonds are formed, so we can confidently assign neutral charges to both the complex and carbon dioxide, avoiding additional fragmentation screening in step 1a (note that this may vary from case to case and caution should always be used).

(b) Natural orbitals for chemical valence (NOCV) analysis. Based on the same considerations made in step 3a, we can obtain the relevant chemically intuitive components of the complex–substrate interactions by decomposing the overall deformation density with the NOCV approach.

Similar to the EDA approach, the resulting NOCV calculation is analogous to that for the M–M′ bond shown in step 1b, since we can decompose deformation densities in terms of NOCV according to eqn (2) and use the ETS-NOCV approach to decompose the orbital interaction energy term accordingly. Note, however, that it is often useful to bring the analysis one step further in this context. As a reminder, the NOCV φ_−k and φ_+k are envisaged as donor and acceptor orbitals, respectively. For the sake of interpretation, a population analysis can also be performed in order to single out, for φ_−k and φ_+k orbitals, which molecular orbitals (MOs) of the two constituting fragments contribute to the interaction (with a resulting associated coefficient accounting for the magnitude of the contribution). In the context of small molecule activation this is particularly useful since it allows to bring the discussion and interpretation at the MO level, which is particularly common in this framework (see Walsh diagram in Scheme 4) by, however, relying on an unbiased approach.

(c) Charge displacement (CD) analysis. The use of the CD analysis (eqn (3)) in this case is less trivial, since multiple substrate–complex interactions occur, and it is often not possible to unambiguously define an integration axis (bond axis). It is possible, however, to use the CD approach to calculate the complex-to-substrate (and vice versa) charge transfer by relying on a slightly different formalism.⁹²

Within this approach, the electron density rearrangement (Δρ′), which typically shows charge accumulation regions (positive values) and charge depletion regions (negative values), defines two different positive functions, Δρ⁺ and Δρ⁻, each equal to the magnitude of the appropriate portion, i.e.:

Δρ^+/− = max[±Δρ(r),0] (5)

so that

Δρ(r) = Δρ⁺(r) = Δρ⁻(r) (6)

By defining two arbitrary regions that are associated with the interacting fragments, we can evaluate the charge transfer as follows:By combining eqn (6) and (7), CT can also be expressed as:Ultimately, this approach can also be applied in the CD-NOCV framework. By combining eqn (2) and (7), we can use this approach for calculating the charge transfer associated to each NOCV deformation density as follows:Despite the spatial regions associated to the two interacting fragments being arbitrarily defined, this approach is particularly suitable for the analysis of the interaction between the complex and the substrate when the two fragments are well-separated in space.

In the case of the gold–aluminyl complex, the combination of the approaches illustrated in steps 3a–c allowed to characterize the interaction between the gold–aluminyl complex and carbon dioxide at TSI in fine detail. We chose TSI since it is the first stationary point along the PES in which significant interactions between the two fragments occur, allowing, however, to unambiguously define them in space.

Firstly, the application of the EDA approach (step 3a) reveals that the overall interaction energy (−12.9 kcal mol⁻¹, Fig. 11) is made favourable (negative) by a strong orbital interaction between the two fragments (−53.3 kcal mol⁻¹). Then, we move to step 3b, where we use the NOCV and ETS-NOCV approaches to characterize qualitatively and quantitatively the nature of such strong orbital interaction.

On a qualitative ground, the NOCV approach is in perfect agreement with the nature of the Au–Al bond we outlined in step 1.

From the NOCV analysis, we find that only two NOCV pairs (k = 1, 2) contribute significantly to the complex–CO₂ interaction. We already showed the isosurfaces for these two deformation densities and discussed the molecular orbital contributions in Fig. 2. We were able to achieve the results by going through step 3b. We visualized at first the donor/acceptor NOCV pairs (φ_−k/φ_+k, k = 1, 2 see eqn (2)) associated to the related deformation densities. By inspection of φ₋₁ and φ₊₁ we confirmed that Δρ₁ is a complex-to-CO₂ charge transfer since φ₋₁ (donor) is mainly localized on Au and Al, while φ₊₁ (acceptor) is almost entirely localized on CO₂. Furthermore, we break down the contributions to both from the MOs of the fragments by population analysis, revealing that the donor φ₋₁ shows mainly contributions from high-lying occupied MOs (HOMO and HOMO−2) of the complex centred on both Au and Al, while the acceptor has contributions mostly from the LUMO of carbon dioxide, as one expects.

Similarly, inspection of φ₋₂ and φ₊₂ confirmed a CO₂-to-complex charge flux for Δρ₂. In particular, the donor φ₋₂ (centred on CO₂) has mainly contributions from the HOMO of the substrate, while the acceptor φ₊₂ is centred on the complex and has contributions from both the low-lying LUMO of the complex and a high-lying unoccupied orbital (LUMO+15), both displaying mainly contributions from the empty 3p_z atomic orbital of Al.

Both step 3b (ETS-NOCV analysis) and step 3c (CD) analysis allow to provide a quantitative perspective on the extent of these interactions, which becomes crucial when comparing reactivity of slightly different complexes (different Au–M complexes, for instance, as we studied in ref. 39). The analyses report that Δρ₁ has an associated ΔE¹_oi of −41.2 kcal mol⁻¹ and an associated CT¹ of 0.33 e, while the second component, smaller in magnitude yet not negligible, has a ΔE²_oi of −4.0 kcal mol⁻¹ and an associated CT¹ of 0.05 e.

Note that these approaches are not new, yet they have proven to be effective in characterizing a wide variety of reactivities (see, for instance, ref. 93 for a recent Feature Article on their use in the framework of frustrated Lewis pairs by Fernandez). Furthermore, some additional tools may be useful to give an even more insightful picture of the forces at play. For instance, in previous works, we also relied on using tools from conceptual density functional theory (CDFT) to describe the molecular electronegativity and nucleophilic/electrophilic areas of the complex,⁹⁴ or intrinsic bond orbitals (IBO)⁹⁵ for a complementary characterization of the metal–metal bond nature. However, while the plethora of approaches within DFT and post-Hartree–Fock methods can surely be useful and insightful, we believe that the composite protocol presented here is a versatile, unbiased, accessible and easy-to-use way to quantitatively characterize homogeneous cooperative small molecule activation processes.

Conclusions and outlook

The efficient activation of small molecules is nowadays a pressing environmental and technological challenge due to their atmospheric abundance and their advantageous chemical reutilization. Cooperativity in small molecule activation has emerged as a particularly flexible and valuable tool to improve these processes. In this framework, an accurate and unbiased computational description of these reactions is essential to further drive the experimental exploration and discovery of novel increasingly efficient processes. In this Feature Article, we present recent experimental advances on cooperative small molecule activation processes mediated by weakly polar or apolar bonds together with our recent extensive work on the theoretical characterization of these reactions. We show how a thorough and unbiased theoretical characterization of the electronic structure of the involved species is essential for rationalizing the experimentally observed phenomena. Furthermore, we also demonstrate how theoretical approaches allow to inductively extract general principles and knowledge on an entire class of processes from single cases, sometimes resulting in unveiling unexpected shifts in previously established paradigms. The experimental work in this field has been intense due to the promising and tunable features of cooperative processes and we have reason to believe that the attention towards these systems will continue to increase. In this scenario, this Feature Article hopefully provides the readership with a broad overview of these processes, as well as with a simple and effective modus operandi to offer solid computational support to experimental discovery. We aim to inspire chemists to explore species containing apolar o weakly polar bonds and expand the range of their applications as catalysts for small molecules activation.

Author contributions

All authors contributed equally to the conceptualization, writing and editing of this Feature Article.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been funded by the European Union – NextGenerationEU under the Italian Ministry of University and Research (MUR) National Innovation Ecosystem grant ECS00000041 – VITALITY – CUP B43C22000470005. Funding by the European Union – NextGenerationEU, AdP POR H2 project L. A. 1.1.35 CUP B93C22000630006 is also gratefully acknowledged. LB acknowledges support from the Ministero dell’Università e della Ricerca (MUR) through PRIN 2022 project “REVOLUTION” (2022HRZH7P).

Notes and references

1. W. B. Tolman, Activation of Small Molecules, Wiley, 2006.
2. D. W. Stephan and G. Erker, Angew. Chem., Int. Ed., 2010, 49, 46–76.
3. D. W. Stephan, Dalton Trans., 2009, 3129–3136.
4. R. D. Dewhurst, M. A. Légaré and H. Braunschweig, Commun. Chem., 2020, 3, 1–4.
5. P. P. Power, Nature, 2010, 463, 171–177.
6. C. Weetman and S. Inoue, ChemCatChem, 2018, 10, 4213–4228.
7. S. Fujimori and S. Inoue, Eur. J. Inorg. Chem., 2020, 3131–3142.
8. C. Camp and J. Arnold, Dalton Trans., 2016, 45, 14462–14498.
9. M. Navarro, J. J. Moreno, M. Pérez-Jiménez and J. Campos, Chem. Commun., 2022, 58, 11220–11235.
10. J. Campos, Nat. Rev. Chem., 2020, 4, 696–702.
11. S. Sinhababu, Y. Lakliang and N. P. Mankad, Dalton Trans., 2022, 51, 6129–6147.
12. Y. Wang, E. Chen and J. Tang, ACS Catal., 2022, 12, 7300–7316.
13. E. Fujita, Coord. Chem. Rev., 1999, 185–186, 373–384.
14. H. Kumagai, Y. Tamaki and O. Ishitani, Acc. Chem. Res., 2022, 55, 978–990.
15. P. R. Yaashikaa, P. Senthil Kumar, S. J. Varjani and A. Saravanan, J. CO2 Util., 2019, 33, 131–147.
16. Webpage: https://gml.noaa.gov/ccgg/trends/monthly.html, accessed 19 October 2023.
17. Webpage: https://webbook.nist.gov/cgi/cbook.cgi?ID=C124389&Mask=1, accessed 18 October 2023.
18. R. Lalrempuia, A. Stasch, A. Stasch and C. Jones, Chem. Sci., 2013, 4, 4383–4388.
19. A. F. R. Kilpatrick and F. G. N. Cloke, Chem. Commun., 2014, 50, 2769–2771.
20. C. E. Kefalidis, C. Jones and L. Maron, Dalton Trans., 2016, 45, 14789–14800.
21. A. F. R. Kilpatrick, J. C. Green and F. G. N. Cloke, Organometallics, 2015, 34, 4830–4843.
22. F. Hanusch, L. Groll and S. Inoue, Chem. Sci., 2021, 12, 2001–2015.
23. J. Li, M. Hermann, G. Frenking and C. Jones, Angew. Chem., Int. Ed., 2012, 51, 8611–8614.
24. A. Caise, L. P. Griffin, C. McManus, A. Heilmann and S. Aldridge, Angew. Chem., Int. Ed., 2022, 61, e202117496.
25. J. Hicks, A. Mansikkamäki, P. Vasko, J. M. Goicoechea and S. Aldridge, Nat. Chem., 2019, 11, 237–241.
26. C. McManus, J. Hicks, X. Cui, L. Zhao, G. Frenking, J. M. Goicoechea and S. Aldridge, Chem. Sci., 2021, 12, 13458–13468.
27. H. Y. Liu, R. J. Schwamm, M. S. Hill, M. F. Mahon, C. L. McMullin and N. A. Rajabi, Angew. Chem., Int. Ed., 2021, 60, 14390–14393.
28. M. M. D. Roy, J. Hicks, P. Vasko, A. Heilmann, A. M. Baston, J. M. Goicoechea and S. Aldridge, Angew. Chem., Int. Ed., 2021, 60, 22301–22306.
29. S. Sinhababu, M. R. Radzhabov, J. Telser and N. P. Mankad, J. Am. Chem. Soc., 2022, 144, 3210–3221.
30. D. Sorbelli, L. Belpassi and P. Belanzoni, J. Am. Chem. Soc., 2021, 143, 14433–14437.
31. A. Hinz, A. Schulz and A. Villinger, Angew. Chem., Int. Ed., 2016, 55, 12214–12218.
32. L. Zhao, F. Huang, G. Lu, Z. X. Wang and P. V. R. Schleyer, J. Am. Chem. Soc., 2012, 134, 8856–8868.
33. Y. Jung, M. Brynda, P. P. Power and M. Head-Gordon, J. Am. Chem. Soc., 2006, 128, 7185–7192.
34. M. Hermann, C. Goedecke, C. Jones and G. Frenking, Organometallics, 2013, 32, 6666–6673.
35. L. Zhao, M. von Hopffgarten, D. M. Andrada and G. Frenking, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, 1345.
36. M. Mitoraj and A. Michalak, J. Mol. Model., 2007, 13, 347–355.
37. A. Michalak, M. Mitoraj and T. Ziegler, J. Phys. Chem. A, 2008, 112, 1933–1939.
38. I. F. Leach, D. Sorbelli, L. Belpassi, P. Belanzoni, R. W. A. Havenith and J. E. M. N. Klein, Dalton Trans., 2023, 52, 11–15.
39. D. Sorbelli, L. Belpassi and P. Belanzoni, Inorg. Chem., 2022, 61, 1704–1716.
40. D. Sorbelli, E. Rossi, R. W. A. Havenith, J. E. M. N. Klein, L. Belpassi and P. Belanzoni, Inorg. Chem., 2022, 61, 7327–7337.
41. J. Hicks, P. Vasko, J. M. Goicoechea and S. Aldridge, Angew. Chem., Int. Ed., 2021, 60, 1702–1713.
42. D. Sorbelli, L. Belpassi and P. Belanzoni, Chem. Sci., 2022, 13, 4623–4634.
43. H. J. Freund and M. W. Roberts, Surf. Sci. Rep., 1996, 25, 225–273.
44. A. Mazheika, Y. G. Wang, R. Valero, F. Viñes, F. Illas, L. M. Ghiringhelli, S. V. Levchenko and M. Scheffler, Nat. Commun., 2022, 13, 1–13.
45. M. K. Singla, P. Nijhawan and A. S. Oberoi, Environ. Sci. Pollut. Res., 2021, 28, 15607–15626.
46. G. J. Kubas, Chem. Rev., 2007, 107, 4152–4205.
47. G. H. Spikes, J. C. Fettinger and P. P. Power, J. Am. Chem. Soc., 2005, 127, 12232–12233.
48. J. Li, C. Schenk, C. Goedecke, G. Frenking and C. Jones, J. Am. Chem. Soc., 2011, 133, 18622–18625.
49. T. J. Hadlington, M. Hermann, J. Li, G. Frenking and C. Jones, Angew. Chem., Int. Ed., 2013, 52, 10199–10203.
50. D. Sorbelli, L. Belpassi and P. Belanzoni, Chem. Sci., 2023, 14, 889–896.
51. D. J. Gorin and F. D. Toste, Nature, 2007, 446, 395–403.
52. M. Devillard, R. Declercq, E. Nicolas, A. W. Ehlers, J. Backs, N. Saffon-Merceron, G. Bouhadir, J. C. Slootweg, W. Uhl and D. Bourissou, J. Am. Chem. Soc., 2016, 138, 4917–4926.
53. E. A. Davidson and D. Kanter, Environ. Res. Lett., 2014, 9, 105012.
54. D. J. Wuebbles, Science, 2009, 326, 56–57.
55. A. R. Ravishankara, J. S. Daniel and R. W. Portmann, Science, 2009, 326, 123–125.
56. M. J. Prather, Science, 1998, 279, 1339–1341.
57. K. Severin, Chem. Soc. Rev., 2015, 44, 6375–6386.
58. F. Le Vaillant, A. Mateos Calbet, S. González-Pelayo, E. J. Reijerse, S. Ni, J. Busch and J. Cornella, Nature, 2022, 604, 677–683.
59. A. Pomowski, W. G. Zumft, P. M. H. Kroneck and O. Einsle, Nature, 2011, 477, 234–237.
60. W. G. Zumft and P. M. H. Kroneck, Adv. Microb. Physiol., 2006, 52, 107–227.
61. L. Zhang, A. Wüst, B. Prasser, C. Müller and O. Einsle, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 12822–12827.
62. S. I. Gorelsky, S. Ghosh and E. I. Solomon, J. Am. Chem. Soc., 2006, 128, 278–290.
63. D. Sorbelli, L. Belpassi and P. Belanzoni, Chem. – Eur. J., 2023, 29, e202203584.
64. E. M. Johnston, C. Carreira, S. Dell’Acqua, S. G. Dey, S. R. Pauleta, I. Moura and E. I. Solomon, J. Am. Chem. Soc., 2017, 139, 4462–4476.
65. I. Bar-Nahum, A. K. Gupta, S. M. Huber, M. Z. Ertem, C. J. Cramer and W. B. Tolman, J. Am. Chem. Soc., 2009, 131, 2812–2814.
66. E. Otten, R. C. Neu and D. W. Stephan, J. Am. Chem. Soc., 2009, 131, 9918–9919.
67. K. Sekine, Gold Bull., 2017, 50, 203–209.
68. C. McManus, A. E. Crumpton and S. Aldridge, Chem. Commun., 2022, 58, 8274–8277.
69. D. Sorbelli, L. Belpassi and P. Belanzoni, Inorg. Chem., 2022, 61, 21095–21106.
70. J. I. Van der Vlugt, Chem. Soc. Rev., 2010, 39, 2302–2322.
71. G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders and T. Ziegler, J. Comput. Chem., 2001, 22, 931–967.
72. F. Neese, WIREs Comput. Mol. Sci., 2012, 2, 73–78.
73. S. Rampino, L. Belpassi, F. Tarantelli and L. Storchi, J. Chem. Theory Comput., 2014, 10, 3766–3776.
74. D. Sorbelli, P. Belanzoni, L. Storchi, O. Bizzarri, B. Bizzarri, E. Mosconi and L. Belpassi, Mol. Phys., 2023, e2245061.
75. Webpage: Energy Decomposition Analysis (EDA)—Tutorials, 2023.1 documentation, https://www.scm.com/doc/Tutorials/Analysis/EDA.html, accessed 18 October 2023.
76. P. Jerabek, P. Schwerdtfeger and G. Frenking, J. Comput. Chem., 2019, 40, 247–264.
77. L. Belpassi, I. Infante, F. Tarantelli and L. Visscher, J. Am. Chem. Soc., 2008, 130, 1048–1060.
78. R. F. Nalewajski and J. Mrozek, Int. J. Quantum Chem., 1994, 51, 187–200.
79. R. F. Nalewajski, J. Mrozek and A. Michalak, Int. J. Quantum Chem., 1997, 61, 589–601.
80. T. Lu and F. Chen, J. Phys. Chem. A, 2013, 117, 3100–3108.
81. M. P. Mitoraj, A. Michalak and T. Ziegler, J. Chem. Theory Comput., 2009, 5, 962–975.
82. G. Bistoni, S. Rampino, F. Tarantelli and L. Belpassi, J. Chem. Phys., 2015, 142, 084112.
83. N. Salvi, L. Belpassi and F. Tarantelli, Chem. – Eur. J., 2010, 16, 7231–7240.
84. Webpage: GitHub – BERTHA-4c-DKS/pycubescd, https://github.com/BERTHA-4c-DKS/pycubescd, accessed 19 October 2023.
85. L. Deng and T. Ziegler, Int. J. Quantum Chem., 1994, 52, 731–765.
86. C. Gonzalez and H. B. Schlegel, J. Phys. Chem., 1990, 94, 5523–5527.
87. D. H. Ess, S. E. Wheeler, R. G. Iafe, L. Xu, N. Çelebi-Ölçüm and K. N. Houk, Angew. Chem., Int. Ed., 2008, 47, 7592–7601.
88. D. Sorbelli, P. Belanzoni, L. Belpassi, J. W. Lee and G. Ciancaleoni, J. Comput. Chem., 2022, 43, 717–727.
89. F. M. Bickelhaupt and K. N. Houk, Angew. Chem., Int. Ed., 2017, 56, 10070–10086.
90. I. Fernández and F. M. Bickelhaupt, Chem. Soc. Rev., 2014, 43, 4953–4967.
91. P. Vermeeren, S. C. C. van der Lubbe, C. Fonseca Guerra, F. M. Bickelhaupt and T. A. Hamlin, Nat. Protoc., 2020, 15, 649–667.
92. G. Bistoni, L. Belpassi and F. Tarantelli, J. Chem. Theory Comput., 2016, 12, 1236–1244.
93. I. Fernández, Chem. Commun., 2022, 58, 4931–4940.
94. P. Geerlings, F. De Proft and W. Langenaeker, Chem. Rev., 2003, 103, 1793–1873.
95. G. Knizia, J. Chem. Theory Comput., 2013, 9, 4834–4843.

---