Contributions of IQA electron correlation in understanding the chemical bond and non-covalent interactions

Abstract

The quantum topological energy partitioning method Interacting Quantum Atoms (IQA) has been applied for over a decade resulting in an enlightening analysis of a variety of systems. In the last three years we have enriched this analysis by incorporating into IQA the two-particle density matrix obtained from Møller–Plesset (MP) perturbation theory. This work led to a new computational and interpretational tool to generate atomistic electron correlation and thus topologically based dispersion energies. Such an analysis determines the effects of electron correlation within atoms and between atoms, which covers both bonded and non-bonded “through -space” atom–atom interactions within a molecule or molecular complex. A series of papers published by us and other groups shows that the behavior of electron correlation is deeply ingrained in structural chemistry. Some concepts that were shown to be connected to bond correlation are bond order, multiplicity, aromaticity, and hydrogen bonding. Moreover, the concepts of covalency and ionicity were shown not to be mutually excluding but to both contribute to the stability of polar bonds. The correlation energy is considerably easier to predict by machine learning (kriging) than other IQA terms. Regarding the nature of the hydrogen bond, correlation energy presents itself in an almost contradicting way: there is much localized correlation energy in a hydrogen bond system, but its overall effect is null due to internal cancelation. Furthermore, the QTAIM delocalization index has a connection with correlation energy. We also explore the role of electron correlation in protobranching, which provides an explanation for the extra stabilization present in branched alkanes compared to their linear counterparts. We hope to show the importance of understanding the true nature of the correlation energy as the foundation of a modern representation of dispersion forces for ab initio, DFT, and force field calculations.

Introduction

Electron correlation and energy partitioning

The London dispersion force [1] (also called dispersion force, London force, instantaneous dipole-induced dipole force, or simply dispersion) is one of the forces that act between atoms in molecules or molecular complexes, and which is important in the description of supramolecular systems. The London force is one of the components of van der Waals force and is typically known to be small compared to the electrostatic force. The dispersion force can keep growing because it is dominated [2] by small additive pairwise interactions. Furthermore, the additive nature of dispersion manifests itself macroscopically, as it is believed to explain the capability of a gecko to effortlessly adhere to smooth surfaces [3] such as glass.

London’s dispersion forces are covered by electron correlation, which is defined as the difference between the exact non-relativistic solution of the Schrödinger equation, and the Hartree-Fock limit, within the Born-Oppenheimer approximation. Electron correlation is typically divided into dynamical and static correlation. While the latter is important for molecules where the ground state is described well only with more than one (nearly-) degenerate determinant, the former is responsible for London dispersion. In order to evaluate the correlation energy of a system one must evaluate the two-particle density matrix, which can be quite a daunting task. The concepts of both dispersion and electron correlation are simple and well understood but there is a gap between the theoretical basis (stabilized many years ago) and the implementation (computationally costly) of these quantities in modern molecular modeling. The purpose of this review is to summarize the work that our group has carried out in an attempt to bridge the gap between theory and implementation of dispersion in molecular modeling. Our strategy involves a better understanding of the real-space partitioning of molecules. By the end of the review, we hope to make clear how we intend to bridge this gap and give examples of how it affects the reliability of predictions.

One of the most glaring examples of the theory-implementation dichotomy is how often quantum chemical methods neglect or poorly represent dispersion, as is the case for Density functional theory (DFT) [4,5,6,7]. Fortunately, DFT methodologies can be amended to include dispersive forces. Two widely used dispersion corrections to DFT are the DFT-Dn method [8,9,10] and the variants of the Tkatchenko–Scheffler (TS) method [11,12,13]. Both schemes basically incorporate an extra energy contribution, E_disp, to an underlying density functional E_DFT, as shown in Eq. 1,

$$ {E}_{total}={E}_{DFT}+{E}_{disp}={E}_{DFT}+\sum \limits_{P<Q}\frac{C_{6, PQ}}{r_{PQ}^6}{f}_{damp}\left({r}_{PQ}\right) $$

(1)

where r_PQ and C₆ are, respectively, the interatomic distance and the sixth-order dispersion coefficient between two atoms (P and Q) and f_damp a damping function to be discussed later.

Although these schemes are well regarded they suffer from a few concerning shortcomings. For example, imagine a situation where the underlying [14] (not corrected) density functional is already able to describe a given system by itself properly, and a dispersion correction is used. In that case, both the $\frac{C_6}{r_{PQ}^6}$ and the damping function f_damp will now be applied to an already well-represented system. That will cause distortions in the natural (and physically meaningful) electronic profile of a system, overestimating the attraction between non-bonded atoms, which is known as overbinding. As a second example of a shortcoming, the expressions used to add dispersion to a functional are not physically robust enough to capture the complexity of dispersion forces in their entirety. The introduction of damping functions is a testament of the shortcomings of these corrections, because the correct potential should dampen itself as a consequence of the natural behavior of the system. Even more concerning, as a third example, is the fact that the inclusion of higher order dispersion coefficients (C₈ and C₁₀, which should account for quadrupole interactions, improving the description of the potential) can in fact reduce the accuracy of the functional [15].

To include dispersion by means of quantum mechanics is a computationally expensive task that requires a detailed and intricate implementation. However, a number of methods now exist to generate electron correlation energies. These methods range from empirical potentials in simulation methods [16,17,18] to the abovementioned a posteriori corrections in DFT, or even ab initio wave function methods [19,20,21]. These methods have been applied to a wide variety of important chemical systems including ones occurring in solid state chemistry [22], biochemistry [23] and, more recently, nanoscience [24]. Recent work has seen a great improvement in methods [25] that account for electron correlation energy. Although electron correlation is proportionally small, especially compared to other terms such as the classical electrostatic interaction, it was shown to be responsible for macroscopic phenomena such as the adhesion of some gecko species to smooth surfaces. While electron correlation energy can be obtained for the whole system, an energy partition schemes can also provide localized atomistic information. The analysis of the correlation energy by itself, which is enabled by partition schemes is particularly interesting to obtain extra insight on dispersion forces.

Partitioning methods are widely reported in the literature. They serve the purpose of extracting chemical insight from systems and to attribute physical meaning to different energy terms. One of the first proposed partitioning methods, the interference and valence partitioning [26] consisted of separating the total energy of a system into (i) an unperturbed reference (quasi-classical or valence part) and (ii) a quantum mechanical interference (two-particle effects). These two terms are still the foundation of many other partitioning schemes. Another example of partitioning, namely, the distributed multipole analysis [27, 28], is limited to Gaussian functions and leaves the actual partitioning decision to an arbitrary criterion of proximity to a nucleus.

Energy partitioning is often used to obtain dispersive forces between two separate molecules (fragments) by symmetry adapted perturbation theory (SAPT) [29]. The limitation here is that dispersion within SAPT is defined only between these fragments as opposed to additionally within a fragment. As a result, the dispersive interaction between the two polarizable termini of a U-shaped chain-like molecule, for example, cannot be defined. However, this problem does not present itself in the Quantum Theory of Atoms in Molecules [30] (QTAIM), nor in the approach of Interacting Quantum Atoms [31] (IQA). Within QTAIM, atomic boundaries are sharp (as opposed to fuzzy) and lead to a space-filling (i.e., no overlap, no gaps) partitioning of the molecular electron density. The electron density can be represented by contour lines (lines whose density is constant), and the gradient paths, which make up the gradient vector field. The points where the gradient vanishes are called critical points, which define topological features (such as bonds, nuclei, rings, etc.). The outer limits of a QTAIM atom are typically defined by a surface for which the electron density is 0.001 a.u.. This criterion ensures that 99% of the electronic charge is enclosed within the atomic boundaries. In contrast, inner boundaries consist of bundles of gradient paths forming the so-called zero-flux or interatomic surfaces.

Finally, we briefly mention that through the IQA approach, the total energy of a molecular system is recovered by the sum of each atomic energy for all atoms in the molecules. The detailed formalism and the equations involved in IQA will be introduced in the next section. As an example of how naturally IQA handles dispersion wherever it occurs (intramolecular or intermolecular) we have recently [32] shown that a non-covalent C...H interaction in monomeric glycine has a dispersive interaction of −0.3 kJ mol⁻¹. This value is very similar to the C...H electron correlation energy associated with a carbon and a hydrogen interacting through-space from different molecules in the ethylene dimer (−0.4 kJ mol⁻¹).

Interacting quantum atoms (IQA)

The Interacting Quantum Atoms [31] (IQA) is an increasingly popular energy decomposition scheme [33] based on QTAIM [30] developed by the Bader group. The IQA approach was proposed in 2005 and was inspired by earlier work [34, 35] published in 2001, which invoked six-dimensional integration to obtain potential energies. IQA splits the molecule’s total energy, ${E}_{IQA}^{Total}$, into a sum of intra-atomic and interatomic energy components. The total molecule energy is then recovered by summing each individual atomic contribution energy, ${E}_{IQA}^A$, according to

$$ {E}_{IQA}^{Total}=\sum \limits_{A=1}^N{E}_{IQA}^A $$

(2)

where N is the total number of atoms in the system and A labels the atoms. Each atomic term can be expanded as a sum of intra- and interatomic contributions such that Eq.2 becomes

$$ {E}_{IQA}^{Total}=\sum \limits_{A=1}^N{E}_{Intra}^A+\frac{1}{2}\sum \limits_{A=1}^N\sum \limits_{B\ne A}^N{V}_{Inter}^{AB} $$

(3)

where ${V}_{Inter}^{AB}$ corresponds to the potential energy between atoms A and B. The intra-atomic encompasses the kinetic, T, as well as the electron-electron, V_ee, and electron-nucleus potential energies V_en

$$ {E}_{Intra}^A={T}^A+{V}_{ee}^{AA}+{V}_{en}^{AA} $$

(4)

${E}_{intra}^A$ can be seen as a measure of the intrinsic stability of an atom in the molecule and turns out [36, 37] to behave like classic steric repulsion in van der Waals complexes, as proven by successful fits to the Buckingham potential. The interatomic contribution can also be split into several terms

$$ {V}_{Inter}^{AB}={V}_{nn}^{AB}+{V}_{en}^{AB}+{V}_{ne}^{AB}+{V}_{ee}^{AB} $$

(5)

where the subscripts e and n stand for electron and nucleus, respectively. The term ${V}_{en}^{AB}$ refers to the potential energy between the electrons from atom A and the nucleus of atom B, whereas ${V}_{ne}^{AB}$ is the potential energy between the nucleus of atom A and electrons from atom B.

The electron correlation term [38] is hidden inside both intra-atomic and interatomic electron-electron potential, as stated by Eqs. 6 and 7,

$$ {V}_{ee}^{AB}={V}_{coul}^{AB}+{V}_{xc}^{AB} $$

(6)

$$ {V}_{ee}^{AA}={V}_{coul}^{AA}+{V}_{xc}^{AA} $$

(7)

Note that ${V}_{xc}^{AB}$ and ${V}_{xc}^{AA}$ both combine the exchange and electron correlation energy. All interatomic terms are obtained direct from a six-dimensional integration of the electron density over the volumes of both topological atoms A and B.

The IQA terms from density matrices [31]

Once the desired multi-electron wavefunction, Ψ, is found, the first- and second-order density matrices can be respectively obtained as follows:

$$ {\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{{\mathbf{1}}^{\prime }}\right)={N}_e\int \varPsi \left({\boldsymbol{x}}_{\mathbf{1}},{\boldsymbol{x}}_{\mathbf{2}},\dots, {\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}}\right){\varPsi}^{\ast}\left({\boldsymbol{x}}_{\mathbf{1}}^{\prime },{\boldsymbol{x}}_{\mathbf{2}},\dots, {\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}}\right)d{\boldsymbol{x}}_{\mathbf{2}}d{\boldsymbol{x}}_{\mathbf{3}}\dots d{\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}} $$

(8)

and

$$ {\rho}_2\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right)=\frac{N_e\left({N}_e-1\right)}{2}\int \varPsi \left({\boldsymbol{x}}_{\mathbf{1}},\dots, {\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}}\right){\varPsi}^{\ast}\left({\boldsymbol{x}}_{\mathbf{1}},\dots {\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}}\right)d{\boldsymbol{x}}_{\mathbf{3}}\dots d{\boldsymbol{x}}_{{\boldsymbol{N}}_{\boldsymbol{e}}} $$

(9)

where x_n is the vector encompassing both the spatial coordinates, denoted by r_n, and the spin coordinates for electron n while N_e is the total number of electrons.

According to the Born-Oppenheimer approximation, the multi-electron Hamiltonian is given by (in atomic units)

$$ {\hat{H}}_{N_e}=\hat{T}+{\hat{V}}_{en}+{\hat{V}}_{ee}=-\sum \limits_{n=1}^{N_e}\frac{1}{2}{\nabla}_n^2-\sum \limits_{n=1}^{N_e}\sum \limits_{A=1}^N\frac{Z_A}{r_{nA}}+\sum \limits_{n=1}^{N_e}\sum \limits_{l>n}^{N_e}\frac{1}{r_{nl}} $$

(10)

where $\hat{T}$ is the one-electron kinetic energy operator, while ${\hat{V}}_{en}$ and ${\hat{V}}_{ee}$ correspond respectively to the attractive electron-nucleus potential and electron-electron repulsive potential. The quantities r_nA and r_nl respectively refer to the distance between electron n and nucleus A, and that between electron nand electron l while Z_A is the nuclear charge associated with atom A. The total energy is given by

$$ {E}_{total}={E}_{electronic}+{V}_{nn} $$

(11)

where the potential energy between two nuclei is simply defined as

$$ {V}_{nn}=\frac{1}{2}\sum \limits_{A=1}^N\sum \limits_{B\ne A}^N\underset{nn}{\overset{AB}{V}}=\frac{1}{2}\sum \limits_{A=1}^N\sum \limits_{B\ne A}^N\frac{Z_A{Z}_B}{r_{AB}} $$

(12)

The electronic energy is obtained solving the Schrödinger equation

$$ {\hat{H}}_{N_e}\varPsi \left({\mathbf{x}}_{\mathbf{1}},...,{\mathbf{x}}_{{\mathbf{N}}_{\mathbf{e}}}\right)={E}_{electronic}\varPsi \left({\mathbf{x}}_{\mathbf{1}},...,{\mathbf{x}}_{{\mathbf{N}}_{\mathbf{e}}}\right) $$

(13)

which results in

$$ {E}_{eletronic}={\int}_{\infty}\hat{\mathrm{T}}{\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{{\mathbf{1}}^{\prime }}\right)d{\boldsymbol{r}}_{\mathbf{1}}+{\int}_{\infty }{\hat{\mathrm{V}}}_{en}{\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{{\mathbf{1}}^{\prime }}\right)d{\boldsymbol{r}}_1+{\int}_{\infty }{\int}_{\infty }{\hat{\mathrm{V}}}_{ee}{\rho}_2\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right)d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(14)

In order to compute individual terms for each atom the topological partitioning method is invoked.

Using only the gradient of the electron density,∇ρ(r), one obtains the so-called interatomic surfaces (denoted S), which are boundaries of an atomic basin (denoted Ω) and are defined by

$$ \nabla \rho \left(\mathbf{r}\right)\cdotp \mathbf{n}\left(\mathbf{r}\right)=0\kern1em \forall \mathbf{r}\in \mathrm{S} $$

(15)

where n(r) is the vector normal to the interatomic surface. With all the atomic basins in the system properly identified, the atomic properties are extracted by integrate over the atomic space. The kinetic energy contribution for atom A can be written as:

$$ {T}^A=\underset{\varOmega_A}{\int}\hat{T}{\rho}_1\left(\mathbf{r}\right)d\mathbf{r} $$

(16)

The monoelectronic interatomic terms are obtained as follows:

$$ {V}_{en}^{AB}={\int}_{\varOmega_A}{\hat{\mathrm{V}}}_{en}^B{\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{{\mathbf{1}}^{\prime }}\right)d{\boldsymbol{x}}_{\mathbf{1}}=-{\int}_{\varOmega_A}\frac{\rho_1\left({\boldsymbol{r}}_{\mathbf{1}}\right){Z}_B}{r_{1B}}d{\boldsymbol{r}}_{\mathbf{1}} $$

(17)

The two-electron terms are obtained from the second-order density matrix, which describes how electrons interact with each other

$$ {V}_{ee}^{AB}={\int}_{\varOmega_A}{\int}_{\varOmega_B}{\rho}_2\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(18)

If A = B, then

$$ {V}_{ee}^{AA}=\frac{1}{2}{\int}_{\varOmega_A}{\int}_{\varOmega_A}{\rho}_2\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(19)

Eq. 18 contains the Coulomb, ${V}_{Coul}^{AB}$, exchange, ${V}_x^{AB}$, and correlation, ${V}_{corr}^{AB}$, energy terms, which are made visible by rearranging,

$$ {V}_{ee}^{AB}={\int}_{\varOmega_A}{\int}_{\varOmega_B}{\rho}_2\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}}={\int}_{\varOmega_A}{\int}_{\varOmega_B}\rho \left({\boldsymbol{r}}_{\mathbf{1}}\right)\rho \left({\boldsymbol{r}}_{\mathbf{2}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}}-{\int}_{\varOmega_A}{\int}_{\varOmega_B}{\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){\rho}_1\left({\boldsymbol{r}}_{\mathbf{2}},{\boldsymbol{r}}_{\mathbf{1}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}}+{\int}_{\varOmega_A}{\int}_{\varOmega_B}{\rho}_2^{corr}{r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}}={V}_{coul}^{AB}+{V}_x^{AB}+{V}_{corr}^{AB} $$

(20)

Exchange and correlation terms are typically lumped together in the exchange-correlation term or ${V}_x^{AB}+{V}_{corr}^{AB}={V}_{xc}^{AB}$. Note that ${V}_x^{AB}$ is related to the covalency degree between atoms A and B, as well the bond order [40]. The classical electrostatic term,${V}_{coul}^{AB}$, on the other hand is related to the bond polarity and the degree of ionicity. The final term, ${V}_{corr}^{AB}$ is responsible to increase the magnitude of the nucleus-electron potential energy and decrease the electron-electron repulsion. All these terms constitute the so-called fine structure [39] of ρ₂. Figure 1 summarizes the energy terms arising in Eq. 20 and provides keywords for each, which helps to explain their respective chemical meaning.

Fig. 1

Representation of the fine structure of the second-order density matrix ρ₂. Each term is presented with its physical meaning. The “+” and “-” signs refer to the mathematical relationship between terms as shown in Eq.20

The exchange-correlation term can be obtained directly from the wavefunction

$$ {V}_{xc}^{AB}=\frac{1}{2}{\int}_{\varOmega_A}{\int}_{\varOmega_B}{\rho}_{xc}\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){r}_{12}^{-1}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(21)

where${\rho}_{xc}\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right)=-{\rho}_1\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right){\rho}_1\left({\boldsymbol{r}}_{\mathbf{2}},{\boldsymbol{r}}_{\mathbf{1}}\right)+{\rho}_2^{corr}\left({\boldsymbol{r}}_{\mathbf{1}},{\boldsymbol{r}}_{\mathbf{2}}\right)$.

IQA dynamic electron correlation energy

The details of our approach have been exhaustively explained in previous publications [32, 41,42,43,44,45]. The IQA method has been developed for over two decades, and was initially applied mainly at Hartree-Fock level [31], but was more recently made compatible with density functional theory [38, 46, 47]. So far there are two methods to account for electron correlation within the Interacting Quantum Atoms [31, 38, 46,47,48,49] (IQA) paradigm: the 6D and the 3D Electrostatic Potential (ESP) integration models. The 6D integration approach was the original partition scheme proposed by our group [50] and generates intra-atomic and interatomic terms

$$ {V}_{corr}^{AB}=\sum \limits_{j=1}^{N_G}\sum \limits_{k=1}^j{k}_{jk}\sum \limits_{l=1}^{N_G}\sum \limits_{m=1}^l{k}_{lm}{d}_{jk lm}^{corr}{\int}_{\varOmega_A}{\int}_{\varOmega_B}\frac{G_{jk}\left({\boldsymbol{r}}_{\mathbf{1}}-{\boldsymbol{R}}_{jk}\right){G}_{lm}\left({\boldsymbol{r}}_{\mathbf{2}}-{\boldsymbol{R}}_{lm}\right)}{r_{12}}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(22)

where N_G is the number of primitive Gaussian basis functions G_jk, where this function is obtained by applying the Gaussian product theorem to the primitive Gaussians G_j and G_k, and similarly for G_lm. The elements$\kern0.50em {d}_{jklm}^{corr}$ correspond to the two-particle density matrix (2PDM) while k_pqis defined as follows:

$$ {k}_{pq}{e}^{-{\alpha}_{pq}{\left|\mathbf{r}-{\mathbf{R}}_{\mathbf{pq}}\right|}^2}={G}_{pq}\left(\mathbf{r}-{\mathbf{R}}_{\mathrm{pq}}\right) $$

(23)

where α_pqcorresponds to the sum of the exponents α_pand α_q of the Gaussian primitives centered at R_pand R_qwhile ${\mathbf{R}}_{\mathrm{p}\mathrm{q}}=\left({\alpha}_p{\mathbf{R}}_{\mathrm{p}}+{\alpha}_q{\mathbf{R}}_{\mathrm{p}}\right){\alpha}_{pq}^{-1}$.

The 2PDM is the result of removing the Hartree-Fock contributions from the second-order Møller-Plesset perturbation theory (MP2) MP2-2PDM. The two 3D volume integrals in Eq. 22 are coupled (in a single 6D integral) because the inter-electronic distance r₁₂ depends both on r₁ andr₂. The interaction energy,${V}_{corr}^{AB}$, of each atom with itself (A = B) or with one of the other atoms (A ≠ B) is obtained from the 2PDM via a 6D quadrature integration. The 3D ESP [51] approach was more recently developed with the intent of ultimate implementation in our polarizable multipolar [52] topological force field FFLUX [53,54,55,56,57], since 3D ESP is faster and more accurate when compared to the 6D approach. The master equation for the 3D ESP integration is given just below, where any derivation details are not repeated here but can be found in the self-contained account of ref. [51],

$$ {V}_{corr}^{A{A}^{\prime }}=\sum \limits_{B\ne A}^N\underset{corr}{\overset{AB}{V}}=\sum \limits_{B\ne A}^N\sum \limits_{j=1}^{N_{basis}}\sum \limits_{k=1}^j{k}_{jk}\sum \limits_{l=1}^{N_{basis}}\sum \limits_{m=1}^l{k}_{lm}{d}_{jk lm}^{corr}{\int}_{\varOmega_A}{\int}_{\varOmega_B}\frac{G_{jk}\left({\boldsymbol{r}}_{\mathbf{1}}-{\boldsymbol{R}}_{\boldsymbol{jk}}\right){G}_{lm}\left({\boldsymbol{r}}_{\mathbf{2}}-{\boldsymbol{R}}_{\boldsymbol{lm}}\right)}{r_{12}}d{\boldsymbol{r}}_{\mathbf{1}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(24)

which is equivalent to

$$ \sum \limits_{j=1}^{N_{basis}}\sum \limits_{k=1}^j{k}_{jk}\sum \limits_{l=1}^{N_{basis}}\sum \limits_{m=1}^l{k}_{lm}{d}_{jk lm}^{corr}{\int}_{\varOmega_A}{G}_{jk}\left({\boldsymbol{r}}_{\mathbf{1}}-{\boldsymbol{R}}_{\boldsymbol{jk}}\right)d{\boldsymbol{r}}_{\mathbf{1}}\sum \limits_{B\ne A}^N{\int}_{\varOmega_B}\frac{G_{lm}\left({\boldsymbol{r}}_{\mathbf{2}}-{\boldsymbol{R}}_{\boldsymbol{lm}}\right)}{r_{12}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(25)

from which follows

$$ \sum \limits_{j=1}^{N_{basis}}\sum \limits_{k=1}^j{k}_{jk}\sum \limits_{l=1}^{N_{basis}}\sum \limits_{m=1}^l{k}_{lm}{d}_{jk lm}^{corr}{\int}_{\varOmega_A}{G}_{jk}\left({\boldsymbol{r}}_{\mathbf{1}}-{\boldsymbol{R}}_{\boldsymbol{jk}}\right)d{\boldsymbol{r}}_{\mathbf{1}}{\int}_{Full\ space}\frac{G_{lm}\left({\boldsymbol{r}}_{\mathbf{2}}-{\boldsymbol{R}}_{\boldsymbol{lm}}\right)}{r_{12}}d{\boldsymbol{r}}_{\mathbf{2}} $$

(26)

resulting in

$$ {V}_{corr}^{A{A}^{\prime }}=\sum \limits_{j=1}^{N_{basis}}\sum \limits_{k=1}^j{k}_{jk}\sum \limits_{l=1}^{N_{basis}}\sum \limits_{m=1}^l{d}_{jk lm}^{corr}{\int}_{\varOmega_A}{G}_{jk}\left({\boldsymbol{r}}_{\mathbf{1}}-{\boldsymbol{R}}_{\boldsymbol{jk}}\right){V}_{lm}^{ESP}\left({\boldsymbol{r}}_{\mathbf{1}}\right)d{\boldsymbol{r}}_{\mathbf{1}} $$

(27)

where k_lm has been absorbed by the ${V}_{lm}^{ESP}$term.

The full space integral in Eq. 26 is solved analytically.

The AA’ (or 3D ESP) approach is a very recent achievement [51]. Atoms can now be calculated 50 times faster with greater accuracy compared to the AB (or 6D) version. However, the price of using the AA’ algorithm is the loss of pairwise inter-atomic electron energies, and thus the insight these bring along. Most of our discussion in this review regards the AB version because we make use of the insight from inter-atomic electron energies. The latter are not required for the force field FFLUX, which only needs the correlation energy of a given atom interacting with all other atoms (i.e., AA’).

Discussion

Quantifying electron bond correlation: Where does it lie?

One of the first and most surprising observations due to the quantification of the electron correlation of the chemical bond is the appearance of a positive bond correlation energy. The positive nature of the bond electron correlation is unexpected because this energy encompasses dispersion effects, which are generally stabilizing and hence come with a negative energy. However, there is a physical explanation for this observation: while only employing the Hartree-Fock field, only Fermi correlation is considered (the Fermi hole [58] follows the electron and exclude same-spin electrons), which means that there is no correlation between electrons with opposite spin states. When Møller–Plesset perturbation is applied, Coulomb correlation becomes present in the system, and opposite-spin electrons are no longer independent [59]. This fact will cause repulsion between electrons, especially in covalent multiple bonds of low polarity. The effect of positive energies is naturally enhanced in multiple bonds, which contain a larger electron density. This observation is supported by the ascending order of the bond correlation energies in the series of F₂ (single bond, 86.9 kJ mol⁻¹), O₂ (double bond, 250.6 kJ mol⁻¹) and N₂ (triple bond, 400.1 kJ mol⁻¹). In fact, the two factors responsible for the bond correlation to be at its most positive are apolarity and multiple bonding. Figure 2 shows the interatomic electron correlation in a number of small representative systems as a function of interatomic distance. This figure does not show F₂, O₂ and N₂ in order to keep the linear energy scale and to present all other systems with sufficient clarity. There is a clear clustering pattern that separates covalent bonds (dark blue), hydrides (light blue), intramolecular hydrogen-hydrogen interactions (pink), hydrogen bonds (red) and weak intermolecular interactions (green). The notation “dim” refer to dimeric species, and “hal bond” to halogen bond.

Fig. 2

Interatomic electron correlation as a function of interatomic distance for a number of representative systems including van der Waals complexes. All values above the dashed red line are positive

In summary, the electron correlation energies shown here do not necessarily represent dispersion in the sense of perturbation theory. The positive values observed for most bonds do not correspond to what is described as a dispersive interaction, i.e., small and stabilizing. Only for larger atomic distances can electron correlation energies be interpreted as dispersion. However, at short interatomic distances, correlation is a measurement of electron repulsion, in response to the electronically populated interatomic region. In most molecules, bond correlation is actually positive and therefore repulsive. The few exceptions are all hydrides and highly ionic molecules, such as NaF and NaCl, for which the bond correlation energies range from 1.5 kJ mol⁻¹ (NaCl) to −9.1 kJ mol⁻¹ (BH₄⁻). It is only in dispersion-dominated van der Waals complexes, such as H₂...H₂ and the helium dimer, that the well-known, weak, long-range and stabilizing nature of dispersion can be seen, with interatomic correlation ranging from a hundredth to tens of kilojoules per mol.

Delocalization index

The nature of covalent bond correlation is reminiscent of the QTAIM delocalization indices [60] (DI), as can be seen from Table 1. This QTAIM delocalization index is a parameter-free, intuitive and consistent way of describing a bond order. The delocalization index is also able to detect the more subtle bonding effects that underpin most practical organic and inorganic chemistry. The difference between the Hartree-Fock and post-Hartree-Fock delocalization indices and bond correlation display the same trend: correlation “moves” electron density from the interatomic region towards the basin for highly populated bonds (covalent), and vice versa for sparsely populated bonds (ionic). This effect is loosely reminiscent of “Le Chatelier Principle” where electron correlation responds as if to desire to ease the burden on heavily populated covalent bonds by taking away some of the electron density of that region, while putting more electron density in the depleted interatomic region, characteristic of ionic bonds. Figure 3 shows this correlation.

Table 1 Bond correlation behavior compared with the QTAIM delocalization index. Both properties are heavily influenced by bond type and polarity

Fig. 3

Bond correlation energies (at fourth-order Møller-Plesset with Single, Double, and Quadruple substitutions MP4(SDQ,FULL)/6–31++G(d,p) level) as a function of the differences between the Hartree-Fock and Configuration Interaction with Singles and Doubles (CISD) delocalization [61] indices

The example of the O₂ molecule illustrates how to interpret Fig. 3 and Table 1. We see a large positive bond correlation energy for O₂’s double bond (y-axis) and a value of 0.83 for the DI_HF-DI_CISD (x-axis) meaning that the bond order is reduced from 2 to 1.17 when electron correlation is added to the molecule through the configuration interaction method.

Aromaticity and anti-aromaticity: Electron delocalization [32]

Now that we have discussed how multiple bonds affect bond correlation we should also address how resonance, aromaticity, and anti-aromaticity manifest themselves in bond correlation. IUPAC defines aromatic and anti-aromatic systems as follows [62]: “a cyclically conjugated molecular entity with stability (due to delocalization) significantly greater than that of a hypothetical localized structure (e.g. Kekulé structure) is said to possess aromatic character. If the structure is of higher energy (less stable) than such a hypothetical classical structure then the molecular entity is ‘anti-aromatic’.”

We have analyzed³¹ the bond correlation energies of both benzene and cyclobutadiene in an attempt to identify how electron delocalization manifests itself through an IQA analysis. For benzene, 62.2 kJ mol⁻¹ of energy was retrieved for all the identical the C-C bonds, while for cyclobutadiene 29.8 kJ mol⁻¹ was recovered for the single C-C bond, and 159.9 kJ mol⁻¹ for the double bond. On average, the C-C bonds in benzene are 32.7 kJ mol⁻¹ more stable than the average of the bonds (two types) in cyclobutadiene. As we stated before, the positive value for bond correlation indicates electron repulsion, which is aggravated by multiple bonds. Due to the delocalized nature of the electrons in the C-C bonds of benzene, the burden of Coulomb electron correlation is distributed over a larger interatomic volume, diminishing the repulsion effect. This effect does not occur in cyclobutadiene because it is a localized system.

Ionicity and covalency or both? The electron correlation answer [45]

The bond correlation energy, as defined by the MPn-IQA framework, deepens insight in chemical bonding method, especially when contrasted with a measure of charge transfer. A simple proposal is to calculate the charge transfer between atoms according to QTAIM net atomic charges. In diatomics, this measure unambiguously gages ionic character. However, a general polyatomic molecule poses the problem that one cannot determine to what extent electronic charge lost by one given atom moved to another atom. In other words, the charge transfer we need in order to characterize a bond is a pairwise property and a set of atomic charges does not reveal pairwise relationships (unless in a diatomic). Nevertheless, polyatomics of the type AB_n (where A and B are elements) allow the pairwise charge transfer we need, by their molecular symmetry. For example, in MgF₂, the more electropositive magnesium has lost electronic charge to the two fluorine atoms in equal measure. Hence, one can unambiguously deduce the degree of charge transfer (and hence ionicity) from the atomic charges: a simple division by two of the Mg charge suffices in this case.

Fig. 4 plots the QTAIM charge transfer against absolute values (a) of bond correlation energy and (b) bond exchange energy. The QTAIM charge transfers were obtained from the absolute charge on the central atom being divided by the number of ligands attached to it. The similarities between panels a and b in Fig. 4 are expected given the reasonable correlation between bond correlation and exchange observed in Fig. 2 of reference [45]. It is important to realize that the bond exchange energy and charge transfer are, in principle, two independent measures. Indeed, charge transfer measures ionicity, the opposite of which should actually be called non-ionicity rather than covalency, as valence bond theory would one guide to do. Covalency is measured by bond exchange, which again presents a sliding scale from (very) covalent (high value of bond exchange energy) to non-covalent. Again, in this two-dimensional world of ionicity and covalency, the opposite of covalency is not ionicity. Indeed, it is perfectly possible to have a bond that is both covalent and ionic: a polar bond. A burning question is then why simple explanations of valence bond theory, and indeed of wider undergraduate education, leave one with the impression that ionicity and covalency are opposites in a one-dimensional world.

Fig. 4

Charge transfer against absolute values of (a) bond exchange and (b) bond correlation. Negative values for bond correlation are marked in purple. Note that charge transfer is defined to be always positive and that there are two distinct bands regarding the number of electrons shared/donated when the molecule is formed

The Protobranching model: Relative stabilities of alkanes [63]

Undergraduate textbooks comment on the relative stability of branched alkanes compared to those of normal (i.e., linear) isomers, and discuss “the alkane branching effect”. McKee and Schleyer [63] noticed that the stabilization factor is related to the 1,3-alkyl-alkyl interactions, as shown in Fig. 5. Each occurrence of that interaction is called a protobranch (pb).

Fig. 5

1,3-alkyl-alkyl interactions, or protobranching, in n-propane (left) and n-pentane (right)

As made clear by McKee and Schleyer, the enhanced stability of branched alkanes can only be properly reproduced by calculation if correlation energies are considered. This is why an analysis of protobranching should be looked at in the context of MPn-IQA.

Fig. 6 was constructed using data from Table 1 of reference [63] and shows the energy changes, ΔE, in kJ mol⁻¹, for the isodesmic reaction where the product corresponds to the branched alkane formed from the equivalent amount of ethane (Eq. 28),

$$ n{C}_2{H}_6\overset{\Delta E}{\to }{C}_{\left(n+1\right)}{H}_{\left(2n+4\right)}+\left(n-1\right){CH}_4 $$

(28)

Fig. 6

Importance of including correlation effects to reproduce experimental values for the energy difference ΔE governing Eq. 28

Setting n = 2 describes how two ethane molecules form propane (which shows protobranching) and a methane molecule. Figure 6 shows the energy difference associated with this reaction, ΔE, calculated using four different methods: HF, B3LYP, B3LYP-D3, and MP2. Hartree-Fock theory does not account for correlation energy while B3LYP only accounts for short-range correlation, and therefore these two methods do not reproduce the energy changes properly as demonstrated by their respective Root Mean Squared Error (RMSE) values of 45.7 and 37.6 kJ mol⁻¹. However, the version of B3LYP as modified with the D3 empirical correction for dispersion and, prominently, MP2 calculations reproduce the experimental energies with low error: RMSE = 20.1 and 4.1 kJ mol⁻¹, respectively. The experimental energies were calculated using heats of formation data at 0 K. McKee and Schleyer calculated the correlation energies as E_MP2 − E_HF and showed that the more branches the alkane possesses, the greater the correlation contribution (see Fig. 2 in reference [63]).

Very recently we conducted a comprehensive, as of yet unpublished, MP4-IQA analysis of protobranching. The results shows that, on average, almost 20% of the stabilizing contribution of the protobranching corresponds to interatomic correlation energy, another 20% corresponds to interatomic exchange energy, while the remaining 60% corresponds to intra-atomic Coulomb energy. However, the stabilizing intra-atomic Coulombic energy is countered by the increase of intra-atomic kinetic energy, which contributes to steric repulsion [36, 37]. Similarly, the stabilizing interatomic exchange energy is also countered by the destabilizing intra-atomic exchange, which makes the interatomic correlation responsible for counterbalancing the Coulombic repulsion between the atoms and providing the extra stabilization property of the protobranching. Figure 7 summarizes this result.

Fig. 7

Average absolute values for the IQA terms of the protobranching in isobutane, isopentane and neopentane. Red bars represent destabilizing factors, while green bars represent stabilizing factors. The double arrows indicate pairs that largely cancel each other

Non-covalent interactions

Hydrogen bond

We have also investigated how electron correlation acts in the context of its most commonly application: non-covalent interactions. The hydrogen bond seems to be a good starting point for this investigation, also inspired by its relevance in biochemical processes. The ubiquitous importance of the hydrogen bond encompasses the peculiar solvent capability of water, DNA binding, and the three-dimensional structures of folded proteins, enzymes and antibodies.

Contrary to chemical common sense, the hydrogen bond can be considered a three-atom system (as opposed to a two-atom system) since there is a considerable amount of exchange between the two electronegative atoms [64] involved in the bond. An electron correlation analysis leads to the same conclusions obtained as for the X...X’ interaction (O…O′, O…N, F…F′ and F…O), because an appreciable amount of positive correlation (an average of +1.9 kJ mol⁻¹) exists between these two atoms. The net effect of these three-atom-system interactions is that there is, to an extent, a cancelation in the dispersive forces in the hydrogen bond, making the “hydration correlation energy” nearly negligible. For the Gly…H₂O complex the intermolecular dispersion energies are actually rather large for hydrogen bond interactions but these specific pairwise interactions tend to cancel one another. Hence, even in a protein encased by a solvation shell, these numerous dispersion interactions will add up to an insignificant dispersion contribution. Table 2 shows the overall cancelation effect of the dispersion energy in hydrogen bonded systems, because the dispersion contributions to the hydrogen bond itself (H…X) are partially canceled by the interaction between the electronegative atoms.

Table  2 Hydrogen bond correlation (H…X) and electronegative-electronegative correlation (X…X’) energies (kJ mol⁻¹) for several systems. The last column shows the high degree of cancelation

Based on this cancelation, most of the cohesion provided by electron correlation in water clusters, comes from O − H bonds and H…H intra-atomic interactions that become more negative. This observation is rather surprising since molecular complexes are commonly thought of being stabilized by non-covalent interactions. In fact, the only non-covalent O…H and O…O interactions involved in the hydrogen bond are larger than 1 kJ mol⁻¹ and tend to largely cancel each other out. Furthermore, IQA electron correlation energies are quite transferable to the water pentamer, especially from the water dimer, which has been energy-topologically (prior to IQA) studied [65] in terms of charge-transfer and polarizability.

Another interesting and important property observed when investigating the electron correlation in water clusters [44] is its high transferability. In particular, the intra-atomic correlation value for the oxygen atom in a single water molecule is −487 kJ mol⁻¹. When acting as hydrogen-bond acceptor in the water dimer, the correlation energy of that same oxygen becomes −495 kJ mol⁻¹, representing a decrease of 8 kJ mol⁻¹. Another decrease of the very same amount in energy is observed when the oxygen accepts another hydrogen-bond, forming the water trimer. In this system the intra-atomic correlation value for the oxygen is equal to −503 kJ mol⁻¹.

The transferability is not restricted to the intra-atomic terms. Indeed, interatomic intramolecular electron correlation energy between the water oxygen and hydrogen, ${V}_{ee, corr}^{OH}$, and between both hydrogen atoms, ${V}_{ee, corr}^{HH}$, show similar transferability. For the water monomer, dimer, trimer, tetramer, and pentamer, the ${V}_{ee, corr}^{OH}$ values range from 12 to 20 kJ mol⁻¹, while ${V}_{ee, corr}^{HH}$ values range from −4 to −2 kJ mol⁻¹. The narrowness of these ranges is remarkable given the variation in the environment of these atom pairs.

The last, but not less remarkable, phenomenon is the (near) cancelation of the interatomic intermolecular electron correlation terms. In the water monomer, the intermolecular ${V}_{ee, corr}^{O{O}^{\ast }}$ term is equal to destabilizing amount of 2 kJ mol⁻¹, whereas the ${V}_{ee, corr}^{O{H}^{\ast }}$term has the stabilizing value of −1.7 kJ mol⁻¹. The terms almost perfect cancel each other. Note that the superscript * is used to indicate that the atoms belong to different molecules. Similar behavior is observed in larger clusters such as trimer, tetramer, and pentamer.

The surprising cohesion provided by hydrogen-hydrogen dispersion

The second largest non-covalent contribution worth mentioning is the intramolecular H...H interaction. We noticed that there is substantial dispersion energy between non-bonded hydrogen atoms, and that interaction is always negative, ranging from −0.3 kJ mol⁻¹ (PH₄⁺) to −4.9 kJ mol⁻¹. This observation highlights the importance of the often neglected role of intra-molecular interatomic dispersion in the stabilization of molecules. In fact, these hydrogen-hydrogen dispersive interactions have been reported before for hydrocarbons [63] in connection with protobranching. Table 3 shows how important these intramolecular interactions can be, and how these interactions are dependent on the internuclear distance between interacting hydrogen atoms. Not only can these interactions be as large as 5 kJ mol⁻¹ in absolute value (as in the case of NH₃) but they can also occur more than once (e.g., in methane where there are 6 H…H interactions resulting in a total energy of 6 x − 4.1 = −24.6 kJ mol⁻¹), thereby much stabilizing these systems.

Table 3 The intramolecular dispersion interactions between hydrogen atoms

Machine learning: The Kriging method [42]

Machine Learning (ML) techniques have become very popular in many fields of science. This popularity stems from its high capacity to generate, analyze, and predict data at a low cost. Although the calculation of electron correlation energies often requires an inspired scientist with enough computational power, the use of ML tools can make such a task easier.

A ML predictive model is built with a training set, which in our case corresponds to a set of atomic coordinates (input) and corresponding energies (output), and a test set that contains the same features, but different points from the training set. The test set is used to evaluate errors and optimize the model.

Three different types of ML techniques, that is, Random Forest (RF), Support Vector Regression (SVR), and Kriging [39, 56, 66,67,68,69], were applied to predict correlation energies of three systems: the water monomer, the H₂-He complex and the water dimer. In the first system, intramolecular-interatomic correlation energies were investigated, while in the last two, the intermolecular correlation energies are in the spotlight, since they are correlated with intermolecular dispersion forces.

The following summary shows how small training sets can be while still producing excellent results.

• The water monomer: for this system, 100 distorted geometries were generated around the equilibrium geometry, with a O-H bond range of 0.86 Å - 1.06 Å and an H-O-H angle range of 65.1^o to 145.9^o. For each geometry the correlation energy was obtained using the program MORFI, an in-house modified version of the program MORPHY [70], and altered to obtain atomic and bond correlation energies, at the MP2/uncon-6-31++G(d,p) level of theory. Forty percent of the data set was utilized as the training set, and 60% as the test set. The optimized Kriging model was able to predict the correlation energies with a RMSE of 0.12 kJ mol⁻¹ in a range of 33.7 kJ mol⁻¹. In 90% of the predictions, the absolute error was lower than 0.01 kJ mol⁻¹. The SVR and RF prediction errors are equal to 2.9 and 3.7 kJ mol⁻¹, respectively.

• The H₂-He complex: a similar procedure was applied here. The only difference is that the training set contains only 15 geometries and the test set 35 geometries. The Kriging method was able to predict the energies with RMSE values bellow 0.01 kJ mol⁻¹ in a range of 3.6 kJ mol⁻¹. The RMSE of the SVR and RF models were respectively equal to 0.4 and 0.3 kJ mol⁻¹.

• The water dimer: A data set containing 100 distorted geometries around the equilibrium was generated and the correlation energies were, once again, calculated at MP2/uncon-6-31++G(d,p) level of theory. In this case, the training set contains 70 geometries, and the test set 30. Kriging predictions were obtained with a RMSE bellow 0.2 kJ mol⁻¹ while SVR and RF delivered to 3.9 and 4.4 kJ mol⁻¹, respectively.

The results presented here highlight the advantages of using ML techniques to construct and investigate potential surface of molecular and supramolecular systems. The correlation energies can be predicted with low error by the Kriging method, requesting only a small training set.

Conclusions

The observations made throughout this review point to two perhaps conflicting yet forthright conclusions: (i) there is much chemical insight contained in the correlated part of the two-particle density matrix (correlation energy), and (ii) relatively little is known about its nature since only a few groups are investigating correlation energy by evaluating the two-particle density matrix. The presence and influence of electron correlation is ubiquitous in a myriad of chemical systems and situations but its importance is often overshadowed by larger electronic effects (i.e., exchange and Coulomb energies) and by its intricate nature, as specific software, and complicated expressions are required to accurately gage atomic electron correlation energies. This situation is disconcerting given the fact that electron correlation contains the coveted dispersion forces. Many of the methods developed to include dispersion to systems, despite experiencing success in varying levels, continue to require extra corrections. Research on the nature of electron correlation has a major impact on the development of methods that will eventually become the foundation of a realistic and robust representation of dispersion forces. Electron correlation energy is proven to be an important contribution inside the IQA framework. The chemical insight carried by this contribution sheds light over the most fundamental aspects of chemistry, including the nature of the chemical bond and aromaticity. Furthermore, the formation of supramolecular systems, such as the water cluster, is highly dependent of how the electron correlation terms between pair of atoms interact with each other.

The high transferability discovered in topologically partitioned correlation energy terms is an important asset, not only for chemistry but also for force field design. Since the beginning of chemical science, chemists were always concerned about grouping atoms, molecules or functional groups by its periodic properties or characteristic behaviors. The transferability of electron correlation energy brings possibilities for refining our knowledge of functional groups, as demonstrated by the protobranching effect. The electron correlation obtained through MPn (n = 2, 3 and 4) densities was shown to have great potential and relevance for all aspects of chemistry. More work is currently being done in order to use our software and capabilities to explain the stabilization of branched alkanes (protobranching) and also to check if the well-known r⁻⁶ behavior for dispersion can be recovered from IQA-formulated correlation. Another running development is the extension of our method to Coupled Cluster densities.