Chemistry and Quantum Mechanics in 2019: Give Us Insight and Numbers.

This Perspective revisits Charles Coulson's famous statement from 1959 "give us insight not numbers" in which he pointed out that accurate computations and chemical understanding often do not go hand in hand. We argue that today, accurate wave function based first-principle calculations can be performed on large molecular systems, while tools are available to interpret the results of these calculations in chemical language. This leads us to modify Coulson's statement to "give us insight and numbers". Examples from organic, inorganic, organometallic and surface chemistry as well as molecular magnetism illustrate the points made.

Introduction

Very few informal talks given by scientists to other scientists ever reach the level of impact comparable to the speech that Prof. Charles Coulson gave at a conference (on Molecular Quantum Mechanics) banquet on June 26, 1959 in Boulder, Colorado. It should be considered as a fortunate circumstance that this speech was later published as an article of high impact.[1] In fact, with the hindsight of almost 60 years that have passed since Coulson gave his famous speech, it is truly astounding how much foresight and insight he was offering at a time where numerical theoretical calculations of molecular electronic structure were very much in their infancy. While the lecture contains a large number of memorable and important points, the sentence that probably stood out most, was Coulson’s statement of despair that he expressed with the words:

“Give us insight, not numbers”

As Coulson pointed out in his lecture, quantum chemistry was about to split into two groups of scientists:

Researchers who care for the most accurate approximate solution to the molecular Schrödinger equation, while not paying (enough) attention to chemical concepts (Coulson’s “electronic computors” or “ab-initio’ists”).

Researchers who are only interested in qualitative concepts and chemical trends while not paying (enough) attention to physical rigor or computability (Coulson’s “nonelectronic computors” or “a posteriori-ists”).

In a similar vein, Malrieu in his article on “Quantum Chemistry and its unachieved missions” stated that “...some tasks, especially the construction of models for a qualitative intelligibility of the molecular world, have been neglected to the benefit of numerism.”[2]

There are definitely situations in which reaching very high accuracy in electronic energies is of critical importance for the success of a theoretical study. For example, high numerical precision is crucial in comparing two transition states that lead to different enantiomeric products[3] or in identifying molecular species in interstellar space that have been detected by high-resolution spectroscopy (e.g., refs (4 and 5)). However, despite all accomplishments in numerical quantum chemistry, the qualitative chemical information content of the calculations should not be neglected.

The number of highly insightful contributions into the structure, bonding and reactivity of molecules that have been made on the basis of creative theoretical reasoning is far too large to do justice to in this short article (e.g., refs (6−9)). However, the central theme of this Perspective is to ask the question whether Coulson’s dichotomy is still valid in 2019. In other words, is it a necessary consequence of computing accurate approximate solutions to the Schrödinger equation that these solutions remain uninterpretable in chemical terms? This question has fascinated theoretical chemists for a very long time. In fact, Klaus Ruedenberg, whose pioneering contributions to both conceptual and numerical quantum chemistry cannot be overestimated, states in his autobiography “... the extraction of correct physical and chemical interpretations from accurate and hence necessarily complex electronic structure calculations, especially as regards bonding, has remained a challenge that has attracted my attention”.[10]

Accurate Numbers

As far as it is known, quantum mechanics describes the material world with perfect accuracy. Hence, solving the molecular (relativistic) Schrödinger equation exactly, is expected to lead to perfectly accurate chemical predictions. While exact solutions of the Schrödinger equation for many particle systems are not possible, systematic approximation methods have been developed, that can approach these solutions with the impressive accuracy of up to 99.999 999% (or 0.01 ppm relative to the total energy of a molecule). This leads to chemical predictions that are accurate to a fraction of a kcal/mol[11−14] These methods are invariably based on systematic expansions of the many particle Schrödinger equation in the framework of coupled-cluster (CC) theory or large-scale configuration interaction (CI) expansions. In these methods, the total energy coconverges with the wave function, the density and all molecular properties to the exact solution. The same, despite all undisputable successes, is clearly not true for density functional theory (DFT) based methods in their practical realizations. (e.g., refs (15−17)).

Concentrating on wave function approaches, conventional wisdom indicates that accurate approaches show a highly unfavorable scaling with system size. Hence, they are hardly applicable to “real-life” chemical problems that involve molecules with, say, 50–200 atoms. In particular, the (often referred to as “gold-standard”) CCSD(T) method[18] is well-known to scale as the seventh power of the system size. Thus, the development of low-order scaling approximations to such high-level wave function methods has been an active field of research ever since the formulation of these methods. We are not able to provide an even cursory discussion of the history and status of these approaches here. However, the past decade has witnessed significant progress in the development of low-order scaling wave function methods. In making approximations, the main difficulty lies in the very high precision that needs to be met in order to not spoil the accuracy of the parent method while, at the same time, still realizing significant computational advantages at realistic system sizes. For example, if one considers a medium sized molecule with a correlation energy of around 10 Eh (≈6275 kcal/mol), it becomes evident that “chemical accuracy” (defined as 1 kcal/mol) is only reached by approximations that recover ≈99.99% of the correlation energy given by the parent canonical method. In practice, one is interested in relative energies and given a realistic cancellation of errors, one can loosely define 99.9% as a reasonable target accuracy.

It has been shown by several groups that such high accuracy can be reached by methods that are based on the so-called Pair Natural Orbital (PNO) expansion.[19−24] PNO methods were invented in the late 1960s and have been exploited with great success in the early and mid 1970s[25−28] before they were abandoned. They were, however, revived in 2009 and have seen a rapid development since then.[20−24] In this short Perspective, we cannot review the developments in this important branch of theoretical chemistry, but rather point out that the methods have been developed to a point where they can be applied to chemistry in almost the same black box fashion that researchers have come to appreciate with DFT calculations.

There is a rapidly growing host of chemical applications of PNO based correlation methods that demonstrate that these methods work in real-life applications and deliver accurate numbers that typically are within chemical accuracy of the parent canonical approaches, while saving many orders of magnitude of calculation time.[20−24] However, the extremely high accuracy that canonical methods reach for small molecules cannot be reached with these methods, since it is difficult to push the remaining error below approximately 0.2 kcal/mol while retaining the computational advantages. As a rule of thumb, with somewhat relaxed wave function truncation thresholds (“NormalPNO”; ∼1 kcal/mol accuracy relative to CCSD(T) in standard test sets[19]), PNO based CCSD(T) calculations are typically not much more expensive than a (hybrid) DFT calculation (with, e.g., the B3LYP functional). With tight thresholds (“TightPNO”, accuracy roughly 0.25 kcal/mol relative to CCSD(T) in standard test sets[19]) they can be about 1 order of magnitude slower. However, given their linear or near linear scaling (linear scaling in linear hydrocarbon chains sets in around 30–40 carbon atoms[22]), local correlation methods are still applicable in a routine fashion to most computational chemistry applications.

A question that is closely related to the subject of the necessary and obtainable accuracy of theoretical predictions, is the overwhelmingly important question of how to falsify theoretical results.[29] Even in 2019, only very small systems can be calculated with an accuracy that leaves no room for doubt. For most contemporary chemically relevant questions, there are always uncertainties that may prevent the calculations from reaching true chemical accuracy. Such error sources include overlooked reaction pathways, unexpected conformers, errors in calculated environmental effects (e.g., solvation, protein environment, embedding treatments, ...), errors stemming from the calculation of entropic contributions or simply unexpectedly complex electronic structures that are inadequately treated by the chosen computational method (e.g., an unbalanced treatment of closed- and open-shell systems or multireference cases) among many others. Thus, a given reaction mechanism or chemical structure can, in general, not be proven through calculations. Rather, quantum chemical calculations can be instrumental in formulating a working hypothesis to be critically tested by experiment. The crucial role of spectroscopy in this context has been elaborated in ref (30).

Quantitative Electronic Structure Methods and Chemical Concepts

Local Energy Decomposition

With the advent of reliable local coupled cluster methods, fairly accurate electronic energies (within the domain of applicability of coupled cluster theory, which is assumed for the remainder of the article unless otherwise stated[31]) can be obtained for realistic systems as a matter of routine. However, returning to Coulson’s famous statement, one has to ask whether the chemical interpretation of the results is possible. Clearly, the complexity of the many particle wave function itself is far beyond what the human brain can process. However, the methods can be very fruitfully combined with a number of analysis tools that allow one to translate the results of the calculations into chemical language. While applications to well-known approaches, such as the natural bond orbital (NBO) analysis[32] or the quantum theory of atoms in molecules (AIM)[33] are straightforward, we will focus here on a recently developed approach that has been termed the “local energy decomposition” (LED)[34] and that belongs to a family of approaches that go back to Morokuma’s pioneering contribution.[35] These methods aim at decomposing the total energy (or an interaction energy) into chemically interpretable quantities such as electrostatic contributions, quantum mechanical exchange and, importantly, the dispersion energy. The latter has entered center stage in chemistry in recent years.[36]

Within the LED scheme, the total interaction energy (ΔEAB) between two fragments A and B can be expressed as[34]

Where Eelstat denotes induced and permanent electrostatic interactions,[37]Eexchange is the quantum mechanical exchange, Edispersion the London dispersion and ΔEprep summarizes a few terms that describe the deformation of the geometric (ΔEgeo-prep, i.e. the so-called “strain” energy) and electronic structure (ΔEel-prep) of the constituents upon complex formation. The ability to clearly identify electrostatic and dispersive components of the interaction energy (with results that are typically consistent with those from other popular energy decomposition techniques[38]) has been found to be particularly illuminating in a number of chemical applications, some of which are briefly discussed below.

Ab Initio Ligand Field Theory

A second class of analysis methods is based on models that have become central to chemical thinking. These models are cast in the language of effective Hamiltonians[39] and include the familiar theory of π-systems that is a cornerstone of organic chemistry or the ligand field theory that has a similar status in coordination chemistry. These theories are cast in terms of relatively simple mathematical models that involve a few adjustable parameters that are determined by fitting experiments. These parameters show characteristic chemical trends that are interpreted in terms of chemical concepts, such as the covalency of metal–ligand bonds or the effective nuclear charge at the metal center.[40] However, the parameters that enter these model Hamiltonians suffer from the lack of a rigorous physical definition. Hence, the challenge is to find a unique connection that allows one to interpret the outcome of an elaborate correlated wave function calculation in terms of the parameters that enter the model Hamiltonian. For ligand theory, it has been shown[40] that there is such a unique connection between multireference perturbation theory calculations (complete active space self-consistent field (CASSCF), followed by N-electron valence perturbation theory (NEVPT2)[41] or complete active space perturbation theory (CASPT2)[42]). The results of CASSCF/NEVPT2 or CASSCF/CASPT2 calculations for optical and magnetic properties of transition metal complexes are frequently in good to excellent agreement with experiment.

The resulting method has been termed ‘Ab Initio Ligand Field Theory’ (AILFT).[40] With the AILFT scheme, one performs a simple CASSCF calculation in which the active space consists of appropriate number n of d- or f-electrons for a given d or f system and the five metal d-based molecular orbitals (likewise seven orbitals for f-elements) and solves for all states of a given spin multiplicity or several spin multiplicities. This CASSCF calculation is followed by a NEVPT2 or CASPT2 calculation. The AILFT now constructs an effective Hamiltonian that is the closest possible to the ab initio Hamiltonian in a least-squares sense. This is an optimization problem that yields the ab initio values for the ligand field matrix V, the Racah parameters B and C and the spin–orbit coupling constant ζ. From the ligand field matrix V, one may obtain ab initio values for the ligand-field splitting 10Dq, possible low symmetry splittings or alternatively, parameters of the angular overlap model (AOM), which decompose the ligand field into σ- and π-contributions from each individual ligand.

It is important to point out that the values deduced in this way are derived in a unique way from the ab initio calculation rather than having been fitted to experiment. Thus, using AILFT, one can obtain deep insight into structure/property relationships along a series of real or hypothetical compounds and, in this way, obtain inspiration for new molecular designs, as will be illustrated below.

In the following, we will exemplify the use of the two above-mentioned analysis methods with the aim to illustrate the central point of this Perspective: Using modern wave function based correlation approaches in conjunction with chemically motivated analysis tools, it is possible to overcome Coulson’s dichotomy and obtain both: accurate numbers and chemical insights without compromising on either goal. We will try to demonstrate this point with a number of case studies from our own laboratory. However, we do not want imply in any way whatsoever that the methods described above are the only useful ones in this context. In fact, as pointed out previously, the methods that do not directly connect to actual observables are classified as “interpretation aids”.[30] Whether a given interpretation aid inspires or guides an individual researcher in their research is a highly subjective matter. Consequently, in our opinion, there is not necessarily “right” or “wrong” when it comes to interpretation aids; they simply provide different flavors of usefulness to a given individual or a community of researchers.

Reaction Mechanisms

Using the emerging local correlation methods, it becomes possible to study chemical reaction mechanisms with higher accuracy and with a higher degree of confidence in the results than was possible before. At the same time, qualitative chemical questions can be addressed that are important in order to derive design principles.

The conversion of CO2 into value-added products and fuels has attracted major attention in recent years. However, due to the thermodynamic stability and kinetics inertness of CO2, efficient CO2 functionalization requires not only high-energy input but also appropriate catalysts.[43] Homogeneous CO2 hydrogenation producing formic acid or formate represents an effective pathway for CO2 transformation.[44] In this regard, impressive reactivity has long been reported for precious-metal catalysts, whereas developments of base-metal catalysts have occurred much more recently.[45]

The CO2 hydrogenation process typically proceeds via a general mechanism shown in Figure a. In the catalytic cycle, molecular hydrogen (H2) first undergoes heterolytic bond cleavage with the assistance of a base to generate metal–hydride species (Int), which then initiates hydride transfer to CO2 to afford the final product, formate as an acid–base complex. The entire reaction thus entails two critical steps, namely, H2 splitting and hydride (H–) transfer, either of which can be the rate-determining step (RDS), as suggested by the experimental mechanistic investigations.[46] To achieve rational catalyst design, one has to first identify the key factors that control the nature of the rate-determining step (RDS) and, more importantly, its barrier height.

Figure 1

(a) Proposed catalytic cycle for CO2 hydrogenation by the representative iron complex Fe/P. (b) Newly designed Fe(II) and Co(III) complexes. (c) DLPNO-CCSD(T) free energy profile of the key reaction steps by Fe/P and Co/P complexes. (d) Correlation plots for the barrier of the key reaction steps and calculated hydricity of Int.

To address this question, we undertook a comparative mechanistic study on the CO2 hydrogenation processes mediated by Fe/P (Figure a), which has been reported to exhibit comparable catalytic activity to noble metals,[47] and its hypothetical Co(III) congener (Co/P) using the DLPNO-CCSD(T) approach. For complex Fe/P, our calculations predicted heterolytic H2-splitting to be the RDS.[48] The computed barrier (24.4 kcal/mol, Figure c) nicely reproduced the experimental value, ∼25 kcal/mol, thus lending credence to the subsequent analyses of the theoretical results, especially for the hypothetical catalysts studied subsequently. In the computations, solvent effects were treated by employing conductor-like screening solvation model at the M06L level of theory, for which methanol was chosen as the solvent, consistent with the experiment. Specifically, for the above-mentioned RDS, the solvent correction only contributes 2.8 kcal/mol to the total barrier (24.4 kcal/mol).[48] By contrast to the DLPNO-CCSD(T) results, none of a wide variety of tested density functionals gave a barrier that was compatible with experiment. With errors generally exceeding 7–8 kcal/mol, the reaction rates predicted by DFT are off by 7 orders of magnitude and more, which renders any quantitative aspect of such modeling questionable.

By contrast to the reaction with Fe/P, a facile H2-splitting process was found for the reaction catalyzed by Co/P, but the subsequent hydride transfer appears unlikely to happen (Figure c). It is clear that changing the metal center from Fe(II) to Co(III) switches the RDS from H2 splitting to hydride transfer.

Typically, as the driving force of a chemical reaction increases, the barrier decreases (Bell–Evan–Polanyi principle). In the present case, the metal–hydride bond (M–H–) is formed in the H2 splitting process, but it is broken in the hydride transfer step. Therefore, the stronger the M–H– bond, the lower the barrier for H2 splitting, whereas the weaker the M–H– bond, the lower the barrier for hydride transfer. The M–H– bonding strength is often quantified by its hydricity or hydride affinity (ΔG°H–), which measures the ability of a metal–hydride species to donate its hydride as MH → M+ + H–.[49] A more positive value of ΔG°H–(MH) (higher hydricity) means stronger M–H– bond and hence diminished hydride donating ability. Apparently, a metal center with a higher oxidation state tends to form a stronger M–H– bond and should be a poor hydride donor. The calculated hydricities for Int(Fe/P) and Int(Co/P) of 58 and 100 kcal/mol, respectively, explain the distinct activity of Fe/P and Co/P for both pivotal steps.

On the basis of the above analysis, a strategy to improve the catalytic activity of the existing catalysts can be proposed. For catalysts with low-hydricity (e.g., complex Fe/P), the RDS is likely the H2 splitting process. Strengthening the M–H– bond by pulling electron density from the metal center would enhance its hydricity and lower the RDS barrier (e.g., complex Fe/P in Figure b). In the case of catalysts with high-hydricity (e.g., complex Co/P), often hydride transfer is the RDS. Weakening the M–H– bond by pushing electron density to the metal center would reduce its hydricity and decrease the RDS barrier (e.g., complex Co/C and Co/Si in Figure b).

To verify the design strategy, the CO2 hydrogen activity of a series of Co(III) and Fe(II) complexes (Figure b) have been examined by using the same computational approach.[50] As shown in Figure d, the computed activation barriers for the H2 splitting and hydride transfer processes nicely correlate with the hydricity of Int. A useful catalyst must strike a delicate balance for both steps to be accomplished efficiently, because the two steps have a just opposite requirement for the hydricity of Int. An optimal hydricity value of 59.7 kcal/mol can be deduced from two linear-fitted lines for the two processes (Figure d). Thus, one can identify that, in addition to Fe/P, Co/C and Co/Si appear to be promising catalysts for CO2 hydrogenation.

Intermolecular Interactions

Weak intermolecular interactions, while known and conceptually understood for a long time,[51,52] have recently been realized to play an essential role in a large variety of chemical phenomena including receptor/effector binding, the relative stability of conformers, solvation phenomena among many others.[53,54] Of particular importance in this context is the London dispersion (LD) interaction, which is an attractive force that is always present in all matter and decays by the inverse sixth power of the intermolecular distance.[55] Unlike electrostatic interactions or specific, directed interactions, such as hydrogen bonds, the LD is mostly isotropic. While the largest part of the electrostatic interactions is already present in relatively simple calculations, e.g. at the Hartree–Fock level, the LD is a pure electron correlation effect that requires a high-level electron correlation treatment to be accurately predicted. Hence, it is of particular interest, from both a chemical and a methodological point of view, to be able to decompose the interaction energy into dispersive and nondispersive components. We will illustrate the insights that can be obtained from such an analysis with two chemical examples from the recent literature.

Frustrated Lewis Pairs

A bulky Lewis acid and bulky Lewis base sterically incapable of forming a Lewis adduct in solution forms a so-called “Frustrated Lewis Pair” (FLP). Intermolecular FLPs can form van der Waals adducts held together by noncovalent interactions, which have been found to catalyze a wide range of transformations involving the activation of small molecules.[56] Hence, understanding the factors that determine whether a classical Lewis adduct is stable or dissociates into a FLP is crucial to design Lewis pairs with tailored bonding features and reactivity.

In ref (57), a series of bulky Lewis pairs was studied with the aim to determine the nature and magnitude of their binding energies (BEs, Figure ). To this end, accurate DLPNO-CCSD(T) calculations were performed that were converged with respect to all technical parameters. Unfortunately, the comparison of the computed BEs with experiment is hampered by the scarce availability of accurate thermodynamic data. For only one FLP system the separate quantification of the free enthalpy (ΔH) and the free energy (ΔG) of association has been reported.[58] In this case, the measured ΔH of −17.5 kca/mol is in excellent agreement with the computed value of −17.1 kcal/mol.

Figure 2

Interaction of a series of phosphines of the type PR3 (R = H, Ph, Cy, Bu and Mes) with the bulky Lewis acid B(C6F6)3. Gas-phase free association energies (ΔG) and enthalpies (ΔH) are computed from DLPNO-CCSD(T) electronic energies and PBE-D3 harmonic frequencies. The LD (Edispersion) and the geometrical preparation (ΔEgeo-prep) contributions to the association energy are also reported.

The chosen computational protocol was then used to calculate highly accurate association energies for a wide range of Lewis Pairs and the LED scheme was used to discuss the role that LD plays in affecting the structural stability of these species. For the Lewis pairs reported in Figure , a Lewis acid B interacts with phosphines of the type PR, with R being H, Ph, Cy, Bu and Mes. In the same figure, the computed ΔH and ΔG values are compared with the LD (Edispersion) and strain energy (ΔEgeo-prep) contributions to the BEs. Remarkably enough, for PR/B complexes, Edispersion values are very similar to the final ΔH ones, thus indicating that LD strongly contributes to the association of all Lewis pairs, especially but not exclusively in the presence of bulky substituents. Its magnitude increases with the size of the substituents on the phosphorus atom (PH < PPh < PCy) for Lewis adducts. The large and repulsive ΔEgeo-prep also increases along the series. In fact, the interacting fragments rearrange significantly in the presence of bulky substituents to facilitate dative bond formation. Hence, LD and polarization effects are both crucial for the stability of these systems, consistent with chemical intuition. By contrast, LD is the only significant component of the interaction in FLPs.

Finally, our results were compared with those obtained at the DFT level for a wide range of dispersion corrected functionals. For FLPs, a good agreement between DFT and DLPNO-CCSD(T) results was obtained. Conversely, most functionals were found to underestimate BEs in classical Lewis adducts, thus limiting the confidence on DFT predictions of the relative stabilities of different species, despite the fact that FLPs and classical Lewis adducts are chemically closely related.

Agostic interactions

Agostic interactions between C–H bonds and coordinatively unsaturated transition metals (TMs) have been a key concept in organometallic chemistry for a long time[59] and their importance for catalysis has been amply documented.[60] Agostic interactions are typically understood in terms of the popular Dewar–Chatt–Duncanson (DCD) bonding model (Scheme ).[61]

Scheme 1

Schematic Representation of the Molecular Orbitals Involved in the Agostic TM···H–C interaction

In this model, an agostic interaction is thought of as arising from a donor/acceptor type orbital interaction in which electron density is transferred from the occupied C–H bonding orbital to empty orbitals on the TM (σ-donation) and from a metal d-based orbital (if occupied) into the empty, C–H σ-antibonding orbital (π-backdonation). Hence, σ-donation and π-backdonation are both expected to weaken the C–H bond and activate it for cleavage. Importantly, being based on an orbital picture, the chemical content of the DCD model should already be present at the HF level of theory, i.e. it will not require electron correlation for its qualitative explanation.

This hypothesis has recently been put to a quantitative test through the LED analysis of accurate DLPNO-CCSD(T) model calculations.[62] The classic system [EtTiCl3(dmpe)] (denoted Ti-1 hereafter; dmpe = 1,2-bis(dimethylphosphino)ethane) was the first β-agostic complex to be experimentally characterized. Its X-ray structure shows that the ethyl moiety distorts to form a close TM···HC contact.[63]Figure shows the comparison between the energy profiles associated with the rotation of the agostic methyl group in Ti-1 around the Cα–Cβ bond computed at the HF and DLPNO-CCSD(T) levels of theory. For θ = 0°(θ is the Ti–Cα–Cβ–Hβ dihedral angle), the system is in its equilibrium geometry and features a close TM···HC contact. The agostic interaction weakens upon methyl group rotation, being absent in the transition state (θ ≈ 60). Thus, the energy barrier for the rotation can be considered as a measure of the strength of the agostic interaction. Unfortunately, experimental rotational barriers are only available for very few systems. In particular, the free energy barrier for the rotation in [EtCo(C5Me5)(PMe3)]+ and [EtCo(C5H5)(PMe3)]+ complexes was measured in solution to be ∼11 and 12.5 kcal/mol, respectively.[64] The corresponding DLPNO–CCSD(T) values in the gas phase are 13.0 and 13.8 kcal/mol, respectively, thus lending credence to the quantitative accuracy of the analysis.

Figure 3

Energy profile for the rotation of the agostic methyl group around the Cα–Cβ bond in the agostic [EtTiCl3(dmpe)] complex at different levels of theory. The reference energy corresponds to the equilibrium geometry (θ = 0). HF values are denoted by gray filled circles while the DLPNO–CCSD(T) ones by black filled circles. Solid lines are spline fits. The vertical arrow represents the correlation contribution to the rotational barrier and the London dispersion component is reported in red.

Remarkably, the minimum corresponding to the agostic structure is not present in the HF potential energy surface. This result emphasizes the importance of dynamic electron correlation in these complexes. It also shows that the DCD orbital interaction model is at least incomplete for the explanation of the agostic interaction since, being a pure orbital model, the effects contained in the DCD model should at least qualitatively be contained in the HF model. However, instead of an agostic attraction, HF shows an agostic repulsion, which means that the major driving force for the formation of agostic structures has not been captured by either HF theory or the DCD model. This conclusion is of general importance and holds true for a wide range of agostic complexes.[62]

Further insight can be obtained by the LED scheme which demonstrates that at least half of the difference between the (wrong) HF and correct (DLPNO-CCSD(T)) curve can be accounted for by LD forces (Figure ). This is a surprising result as it demonstrates that quantitatively, agostic interactions cannot be explained by simply looking at orbital interactions but instead, the LD must be considered as well. This is relevant for the chemical design of agostic interactions given that LD has a different angular and distance dependence than orbital interactions that are highly directed and of very short-range. Importantly, LD does not lead to bond activation since it merely involves the interaction of fluctuating dipoles and not the transfer of charge. Hence, the results of ref (62) also explain why a large variety of agostic structures, including the textbook case just discussed, do not show any experimental evidence of significant C–H activation.[65]

Solids and Surfaces

In the field of solids and surfaces adsorption energies are central quantities for identifying key intermediate species in heterogeneous catalysis.[66−70] However, reliable and quantitative experimental values for adsorption energies on well-defined surfaces at low coverage limit are rather scarce.[71] Hence there is an increasing need to develop combined experimental and theoretical protocols with strong predictive ability. This employs theoretical methods that are reliable, robust and accurate, while they are applicable to large system sizes. This field of research is clearly dominated by DFT. A wide range of functionals have been developed for this purpose ranging from semilocal exchange and correlation functionals up to “higher-rung” functionals like the screened hybrid functionals as well as diagrammatically derived functionals based on the random phase approximation (RPA). This necessitates the importance of a systematic benchmarking against accurate reference numbers in order to access the accuracy of any DFT functional.

With the advent of local correlation methods in conjunction with embedding techniques, it recently became possible to perform “gold standard” CCSD(T) level calculations for surface systems[72] A point in case is a recent study that investigated small molecule binding to TiO2 surfaces (Figure a,b). For this model system, there exist experimental measurements of the binding energies of a variety of small molecules (H2O, NH3, CH4, CH3OH, CO2). It was shown that DLPNO-CCSD(T) calculations, if carefully done, reproduce all experimental data within the error bar of the measurements (Figure c,d).[72]

Figure 4

Graphical representation of the quantum region of the employed embedded clusters to represent the binding of (a) water (cluster Ti17O34–H2O) and (b) methane (cluster Ti17O34–CH4) over the rutile TiO2 (110) surface. The DLPNO-CCSD(T) energies are analyzed within the LED scheme to provide estimation for the respective Ti17O34–H2O and Ti17O34–CH4 electrostatic and dispersion interactions. Average experimental estimated adsorption energies (red-dot line) versus various DFT, DLPNO-MP2 and DLPNO-CCSD(T) (sticks) computed zero-point energy corrected adsorption energies for clusters (c) Ti17O34–H2O and (d) Ti17O34–CH4, respectively. Color-coding: Ti, light-gray; O, red; C, dark-gra; H, white.

For DFT functionals, the situation is more complex. As it is seen in the case of the water binding a systematic convergence toward experiment is observed in the sequence of GGA (PBE-D3, rPBE-D3), Hybrid (PBE0-D3), and Double Hybrid (B2PLYP-D3) functionals as well as DLPNO-MP2 and DLPNO-CCSD(T). Of these, DLPNO-CCSD(T) is the most accurate and deviates by only 0.1–0.2 kcal/mol from the experimentally estimated adsorption energy. In the case of methane, the computed DLPNO-CCSD(T) adsorption energy slightly overestimates the experimental value, deviating only by about 1.1–1.3 kcal/mol. By contrast, it has been shown that large variations between functionals are observed in the case of CH4 adsorption. This is indicated by the spread of ∼4 kcal/mol in the computed zero-point energy corrected adsorption energies (Ead,0) which amounts to ∼50% of the DLPNO-CCSD(T) computed Ead,0.[72] The most successfully performing functionals are presented in Figure d. In accordance to DLPNO-CCSD(T) the various DFT computed adsorption energies deviate also by about 1.4–1.1 kcal/mol. However, the Jacob’s ladder hierarchy does not hold in the case of the methane binding as the different functionals are showing either overbinding or underbinding behavior with no recognizable pattern. In fact, different minima of the computed potential energy surfaces (PES) are predicted by different functionals.[72] DLPNO-MP2 shows, as expected, a tendency for overbinding. Thus, only the DLPNO-CCSD(T) calculations were in systematic and quantitative agreement with experiment.

These results set the basis for a deeper understanding on the nature of the surface–adsorbate chemical interactions. The LED analysis reveals that in the case of Ti17O34–H2O the Ti–H2O bonding interaction is mainly electrostatic (73%) with the remaining 27% accounting for dispersion interactions (Figure a). However, the situation is reversed in the case of the Ti17O34–CH4 in which the Ti–CH4 bonding interaction is mainly due to dispersion interactions (68%) with the remaining 32% accounting for electrostatic interactions (Figure b). Thus, opposite to the case of FLPs and classical Lewis Pairs, DFT performs better for the mostly electrostatically bound molecules than for the dispersion bound systems.

Coordination Chemistry and Molecular Magnetism

Molecular magnetism is concerned with the properties of paramagnetic open-shell 3d to 5d transition metal or 4f(5f) lanthanide(actinide) coordination compounds(complexes). It has a significant impact on closely related disciplines such as molecular electronics and chemical reactivity (homogeneous and heterogeneous catalysis).[73] Depending on the number of unpaired d- or f-electrons (n), complexes in their high-spin ground states possess a total spin of S = n/2. This spin is isotropic in the sense that it aligns readily along any direction of an external magnetic field. Magnetic anisotropy arises in complexes with axially symmetric coordination geometries. It is described by zero-field splitting (ZFS) parameter D of the 2S+1 sublevels M of the spin S within a given electronic state.

The basic goal in the field of molecular magnetism is to provide molecules with magnetic moments that are highly anisotropic (large negative D); when induced by an external magnetic field, the magnetic moment persists after switching-off the field for given time (the relaxation time, single molecule magnets, SMM). The “holey grail” in the molecular magnetism is hence to increase the relaxation time to, ideally, room temperatures in order to function as a miniature switch in devices.[74]

Despite intense research efforts, molecules that show magnetic blocking at room temperature have not been found and it is unlikely that approaches that fully rely on serendipity can succeed in this endeavor. The first SMM, a poly oxo metalate with a Mn12 core and a S = 10 ground state, was discovered.[75] Using Mn12 as a lead structure, it was believed for almost two decades that better SMMs could be found by increasing the number of unpaired electrons and with it the overall ground state spin. It was later shown, partially by theoretical reasoning,[76,77] that large, oligonuclear clusters are not a necessity for SMM behavior. Subsequently, much recent progress has been obtained on the basis of much simpler and synthetically far more easily controllable mono- or dinuclear nuclear transition metal or f-element systems.

Hence, it appears obvious that it is necessary to have first principle approaches with predictive power available that, at the same time allow for chemical insights to be derived that eventually lead to new design principles. AILFT coupled to CASSCF/NEVPT2 calculations is such a tool, as will be demonstrated below with one recent example. Many others can be found in the literature (e.g.,[78−83]).

Quite surprisingly, Cobalt(II) tetra-thiolate [Co(SPh)4]2– was reported as a first example of a mononuclear SMM that shows a slow relaxation of the magnetization in the absence of an external magnetic field.[84] Interestingly, crystal structures of [Co(SPh)4]2– with a variety of different counterions have been reported. Even more surprisingly, the magnetic properties of these complexes differ drastically. While [Co(SPh)4][P(Ph)4]2 shows a powder EPR spectrum that is indicative of S = 3/2 system with a large negative ZFS, the corresponding [Co(SPh)4][N(Et)4]2 salt shows a much more “normal” S = 3/2 EPR spectrum with a small and positive ZFS. Consequently, only [Co(SPh)4][P(Ph)4]2 shows magnetic blocking and SMM behavior. Such a drastic change in the magnetic properties of a given compound following a subtle chemical variation is perhaps unprecedented. Understanding this behavior potentially opens new routes for the design of SMMs.

Subsequently, both forms of the complex were subjected to very detailed experimental studies including magnetic susceptibility, high-frequency EPR, magnetic circular dichroism as well as far-infrared spectroscopy.[78] All of these measurements can be interpreted in terms of a model, in which [Co(SPh)4][P(Ph)4]2 shows a nearly axial ZFS with a D-value of around −55 cm–1 while [Co(SPh)4][N(Et)4]2 features a more rhombic ZFS tensor with a D-value of around +9 cm–1. CASSCF/NEVPT2 calculations done on the crystal structures of both compounds lead (after reoptimization of the hydrogen positions) to near quantitative agreement with experiment (Figure ).

Figure 5

Temperature dependence of the molar magnetic susceptibility of the (Ph4P)2 [Co(SPh)4] from experiment and simulated using CASSCF and NEVPT2 results.

However, the actual cause of the highly peculiar behavior of the two compounds only became intelligible following the AILFT analysis of the CASSCF/NEVPT2 results. Using AILFT, the effective splitting of the d-orbital manifolds could be deduced. What emerges from the calculations is that [Co(SPh)4][P(Ph)4]2 and [Co(SPh)4][N(Et)4]2 differ qualitatively in terms of the low-symmetry distortions away from pure tetrahedral symmetry. [Co(SPh)4][P(Ph)4]2 shows a tetragonal elongation which leads to a splitting pattern in which the d and d are higher in energy than the d orbital. The opposite is true for [Co(SPh)4][N(Et)4]2, which shows tetragonal compression and the reverse orbital splitting pattern. These opposite geometric distortions lead to qualitatively different spin–orbit interactions between the 4A2 ground state and the first excited 4T2 term. The [Co(SPh)4][P(Ph)4]2 leads to preferential stabilization of the MS = ±3/2 magnetic sublevels thus defining a negative ZFS, while the [Co(SPh)4][N(Et)4]2 splitting pattern preferentially stabilizes MS = ±1/2 thus defining a positive ZFS.

This peculiar difference was further investigated by calculating a two-dimensional potential energy and property surface of the [Co(SPh)4]2– ion as a function of the tetrahedral angle S–Co–S and the C–S–Co–S dihedral angle describing the tilt of the phenyl moiety (Figure ). The potential energy surface clearly shows two minima that are connected by a relatively low-energy transition state and that correspond to an elongated and a flattened tetrahedron, respectively. Consistent with the qualitative analysis, the calculated D-value at the elongated tetrahedron minimum is large and negative while the D-value at the flattened tetrahedron minimum is small and positive.

Figure 6

Left: Definition of the two angles θ1 and ψ1 define two-dimensional potential energy surface. Right: Potential energy and property surface (colors represent the relative energy in kcal/mol) together with calculated ZFS isolines (black lines, numbers in cm–1) computed by SOC-CASSCF(7,5) for optimized structures of the model [Co(SCH3)4]2– (left side modified after Figure 2 of ref (86); right side modified after Figure 12 of ref (78)).

The picture that emerges from this analysis is that during crystallization the weak intermolecular interaction between the [Co(SPh)4]2– and the counterions lead the P(Ph)4-salt to lock into an elongated minimum and the N(Et)4 salt to a flattened minimum. Thus, a subtle twist of a phenyl ring is all that it takes to turn the magnetic properties of [Co(SPh)4]2– completely around and determine the difference between a normal coordination compound and a SMM.

The reason for the existence of the two minima is readily understood from the electronic structure of the Ph–S– ligand (Figure ). There are two mainly S-centered orbitals that bind to the metal ion: the in-plane (ip) lone pair forms a σ-interaction with the metal and the out-of-plane (oop) lone pair interaction in a π-fashion. However, the oop lone pair strongly interacts with the π-system of the phenyl ring and consequently, this orbital is delocalized over the sulfur and phenyl moieties. This, in turn, leads this orbital to be “rigidly” oriented perpendicular to the plane of the phenyl ring. Consequently, small rotations of the phenyl moiety immediately lead to strong changes in the π/σ-interaction of the sulfur ligand with the central metal (“misdirected valence”) thus changing the orbital splitting pattern and leading to the occurrence of two minima on the PES.

Figure 7

Ligand orbitals relevant to the metal–ligand interaction in the case of the sp2 hybridized and nonhybridized ligand orbitals (reproduced from Figure 4 of ref (86)).

There are several lessons to be learned from this study. Perhaps the most important conclusion is that subtle chemical variations in the second coordination sphere can be rationally employed to strongly influence the magnetic response of the system. While the counterion effect described above is based on serendipity, we see no reason why chemists should not be able to exploit the sensitivity of the ZFS to conformation via intelligent ligand design. In fact, alternative ligand systems have been explored and lead to even better Co(II) based SMMs.[82,86]

Conclusion

We hope that in this short Perspective article, we were able to demonstrate that quantum chemistry has come a long way since Charles Coulson’s famous speech from 1959. It is indeed now possible in many cases to compute accurate energies and accurate properties from first-principles wave function based approaches on realistic systems that contain dozens to a few hundred atoms. Moreover, the results of such calculations can be interpreted in a chemical language using a variety of tools. The examples given were selected to provide a glimpse of the many different ways accuracy can meet chemical insights in a wide variety of chemical areas. Importantly, the examples are concrete, real-life studies that are of contemporary chemical interest.

We thus believe that it is fair to modify Coulson’s statement by one word to read:

“Give us insight and numbers”[87]

This statement is meant to emphasize that accurate numbers, while nowadays achievable for much larger systems than previously possible, should not be viewed as the only and ultimate goal of a theoretical investigation. The chemical insights that arise from it are at least as important as they fulfill one of the, if not the single most important missions of theory, to inspire and guide new experiments. Chemical insights do not automatically present themselves as the results of a calculation, but require additional human effort. Appropriate methods to obtain accurate numbers and the tools to interpret them in chemical language are widely available.