Molecular Interactions From the Density Functional Theory for Chemical Reactivity: The Interaction Energy Between Two-Reagents.

Reactivity descriptors indicate where a reagent is most reactive and how it is most likely to react. However, a reaction will only occur when the reagent encounters a suitable reaction partner. Determining whether a pair of reagents is well-matched requires developing reactivity rules that depend on both reagents. This can be achieved using the expression for the minimum-interaction-energy obtained from the density functional reactivity theory. Different terms in this expression will be dominant in different circumstances; depending on which terms control the reactivity, different reactivity indicators will be preferred.

1 Introduction

When a reagent approaches a reactant molecule, the reactant is perturbed. If the perturbing reagent stabilizes the molecule (lowering its energy) or destabilizes the molecule comparatively little (raising its energy by a relatively small amount), then a reaction may occur. For example, regioselectivity can be predicted by finding the location where the reagent imparts the greatest additional stability to the reaction substrate. By comparing results for a series of reagents and a given molecule, or for a series of molecules and a single reagent, absolute chemical reactivity trends can be understood, and even predicted.

This basic viewpoint on chemical reactivity is the fundamental basis of the density functional theory for chemical reactivity (DFT-CR), often called chemical density functional theory or conceptual DFT (Parr and Yang, 1989; Chermette, 1999; Geerlings et al., 2003; Ayers et al., 2005a; Chattaraj et al., 2006; Gazquez, 2008; Liu, 2009a; Johnson et al., 2012b; De Proft et al., 2014; Fuentealba and Cardenas, 2015; Miranda-Quintana, 2018). The fundamental mathematics of DFT-CR is functional perturbation theory (Ayers et al., 2005a) using whichever representation of the chemical system (or “ensemble”) is deemed to be most convenient (Nalewajski and Parr, 1982; DeProft et al., 1997; Liu and Parr, 1997; Perez et al., 2008; Cardenas et al., 2009a). In this paper we will only consider the most common, “closed-system picture” (canonical ensemble). However, the same approach could be used for any of its Legendre (Johnson et al., 2012a) transforms (Ayers and Fuentealba, 2009).

In the closed-system picture, the state of a molecule is specified by its number of electrons, N, and the external potential that binds those electrons, v(r). When a reagent approaches, the number of electrons changes (from electron transfer between reagents) and the external potential changes (because the electrons in the molecule are attracted to the nuclei and repelled by the electrons of the approaching reagent). The change in the energy of a molecule with non-degenerate ground state (Cardenas et al., 2011a; Bultinck et al., 2013, 2014; Bultinck et al., 2015) is then,

Here ΔV is the change in nuclear-nuclear repulsion energy, ΔN is the number of electrons transferred from the reagent to the molecule under scrutiny and Δv(r) is the change in the molecular external potential due to the presence of the reagent. These quantities depend on the particular identity of the reagent. The coefficients in this expansion—the chemical potential μ, (Parr et al., 1978) the chemical hardness η, (Parr and Pearson, 1983; Pearson, 1997; Ayers, 2007b) the Fukui function f(r), (Parr and Yang, 1984; Yang et al., 1984; Ayers and Levy, 2000; Ayers P. W. et al., 2009) the dual descriptor Δf(r), (Fuentealba and Parr, 1991; Morell et al., 2005, 2006; Ayers et al., 2007; Cardenas et al., 2009b; Geerlings et al., 2012) the electron density ρ(r), and the linear-response (or polarizability) kernel χ(r,r’) (Berkowitz and Parr, 1988; Senet, 1996; Ayers, 2001; Sablon et al., 2010; Geerlings et al., 2014; Franco-Perez et al., 2015a)—are often taken to depend only on which molecule we are studying, and not on what reagent is being considered. The coefficients in the expansion determine how much (or how little) the reactant of interest is stabilized by a given reagent. Since energetically favorable interactions are associated with high reactivity, the coefficients in this expansion are reactivity indicators for the molecule being studied (Ayers et al., 2005a).

The assumption that the fine details of the reagent do not control the qualitative aspects of the reactivity is implicitly invoked when the coefficients of the Taylor expansion are used as reactivity indicators. This is often an adequate assumption. The mere idea that we can define scales of acidity and basicity (proclaiming, for example, that hydrochloric acid is a stronger acid than acetic acid) presumes that the properties of reactants (in this case, the acidity) do not depend very strongly on the identity of the attacking reagents (the base) or the solvent in which the reaction occurs. Similarly, when we say that nitrobenzene is a deactivating meta-director for electrophilic aromatic substitution we imply that the regioselectivity properties of nitrobenzene are largely independent of the choice of electrophile. Extending these ideas: most molecules are reactive at just one (and certainly no more than a few) molecular site(s); these reactivity preferences recur across a wide range of reagents, solvents, and experimental conditions. In cases like this, just a few coarse details of the reagent (e.g., the reagent’s electronegativity and hardness, which determine whether ΔN is small or large) and the reagent’s charge will suffice to describe the full range of possible reactions. In these cases, applying DFT-CR is as simple, and as complicated, as finding an appropriate coefficient (or linear combination of coefficients) from Eq. 1. The “perturbative perspective” on chemical reactivity says by choosing simple “model reagents” (often just a point charge with some capacity for accepting/donating electrons), one can decipher the reactivity preferences of the molecule of interest. There are many examples where this approach has been successfully employed to define new reactivity indicators (Parr et al., 1999; Ayers and Parr, 2000; Ayers and Parr, 2001; Ayers et al., 2005b; Anderson et al., 2007a; Anderson et al., 2007b; Cedillo et al., 2007; Gazquez et al., 2007; Cardenas et al., 2011b; Cardenas, 2011; Alain Miranda-Quintana R. and Ayers P. W., 2016), (Bochicchio, 2015; Alain Miranda-Quintana et al., 2016; Miranda-Quintana R. A., 2016, Miranda-Quintana, R. A. 2017; Gonzalez et al., 2018). This simple perspective has been successfully applied to the study of redox reactions, the study of solvation processes, and to the reactivity of aromatic species (Miranda-Quintana and Smiatek, 2020; Miranda-Quintana et al., 2022).

Chemistry is not always this simple. Certainly some reactions happen (or at least happen more readily than expected) because the reactant molecules are well-matched. Classic cases include the hard/soft acid/base principle and the Woodward-Hoffmann rules (e.g., ethene and butadiene are chemically similar, yet ethene is susceptible to cycloaddition with butadiene but not with itself). In the case of the hard/soft acid/base principle, treating both reagents as simplified “model reagents” suffices (Ayers, 2005; Ayers et al., 2006; Anderson et al., 2007a; Ayers, 2007b; Cardenas and Ayers, 2013; Miranda-Quintana, 2017a). The Woodward-Hoffmann rules can be treated using the DFT-CR framework, (De Proft et al., 2006; Ayers et al., 2007; De Proft et al., 2008; Sablon et al., 2009; Geerlings et al., 2012) but the fundamental reason that the dual descriptor and the initial hardness response suffice is still somewhat mysterious (Cardenas et al., 2009b). A more complete model for the interaction energy is clearly needed in this case. In the next section we will define just such a model. Section 3 discusses the key results. In follow-up papers we will derive an expression for the initial hardness response and reveal a new link between this quantity and the dual descriptor.

2 An Interaction Energy Model From Density Functional Theory Perturbation Theory

Our analysis is based on Eq. 1, including only the terms that are shown explicitly. Higher-order derivatives with respect to the number of electrons are usually small. (Fuentealba and Parr, 1991; Geerlings and De Proft, 2008; Morell et al., 2013; Heidar-Zadeh et al., 2016b; Munoz and Cardenas, 2017), and although they have attracted some recent interest, we will not take them into account here (Hoffmann et al., 2020; Miranda-Quintana R. A. et al., 2021; Miranda-Quintana R. A. et al., 2021). One may argue that since there are only one-body and two-body terms in the Hamiltonian, the energy and its functional derivatives with respect to v(r) should be nearly quadratic in the number of electrons (Ayers and Parr, 2008). Higher order derivatives with respect to the external potential correspond to hyperpolarizabilities; (Fuentealba and Parr, 1991; Senet, 1996; Cardenas et al., 2009a) these effects are difficult to compute accurately and are hopefully negligible. We assume that all the derivatives exist. Treating them exactly requires working with the system at finite temperature: (Franco-Perez et al., 2015a; Franco-Perez et al., 2015b; Malek and Balawender, 2015; Franco-Perez et al., 2016; Franco-Perez et al., 2017a; b;Franco-Perez et al., 2017c; Franco-Perez et al., 2017d; Franco-Perez et al., 2017e; Miranda-Quintana, 2017c; Miranda-Quintana RA. et al., 2017; Miranda-Quintana R. A. et al., 2017; Polanco-Ramirez et al., 2017; Franco-Pérez et al., 2018; Robles et al., 2018; Gázquez et al., 2019): since there are derivative discontinuities at the zero-temperature limit (Perdew et al., 1982; Yang et al., 2000; Ayers, 2008; Yang et al., 2016). One may argue that, for interacting molecular fragments, the “effective temperature” is nonzero; cf. refs. (Ayers, 2007a; Alain Miranda-Quintana R. and Ayers P. W., 2016).

The perturbation expansion, (1), gives us information about the energy of A in the presence of B and the energy of B in the presence of A . We need to add these two interaction-energy expression and make corrections for 1) the double-counting of interactions and 2) nonadditive energy contributions like exchange and correlation between A and B.

The electron density of the interacting reagents, ρ (r), is divided into electron densities for the subsystems,

It is convenient to also define the nuclear charge density of the reactants, i.e., The total charge density of reactant A is then . There is a similar expression for reactant B. The external potentials of the reactants are given by expressions like

The interaction energy can then be written in terms of the fundamental density functional for the energy, (Hohenberg and Kohn, 1964) Here and denote the electron densities of the separated reactants and denotes the change in nuclear-nuclear repulsion energy when the reactants are brought together. The contribution from rearrangements of the nuclei within the fragments is neglected for simplicity; this assumes that in the initial approach between the reagents, the internal molecular geometries of the individual reagents are preserved. This approximation was suggested, for example, in ref. (Ayers and Parr, 2001). This approximation can be relaxed if a term corresponding to the “preparation” energy of the reagents is added (Morokuma, 1971; Kitaura and Morokuma, 1976; Ziegler and Rauk, 1977, 1979; Wu et al., 2009; Azar and Head-Gordon, 2012).

The electronic energy functional is then decomposed,

Here is the Hohenberg-Kohn functional, (Hohenberg and Kohn, 1964) is the Kohn-Sham kinetic energy functional, (Kohn and Sham, 1965) and is the exchange-correlation energy functional. is the variational energy functional in DFT (Hohenberg and Kohn, 1964; Kohn and Sham, 1965)

Equation 7 can be simplified by breaking each contribution to the energy into its subsystem contributions:

Using the definition of the electronic energy functional and Eq. 4, this simplifies to

The terms in the first two lines are identified as the non-additive kinetic energy and the non-additive exchange-correlation energy,

The next term in Eq. 11 is the “non-additive Coulomb energy” this models the electrostatic interactions between the systems. The final terms in Eq. 11 include the polarization and electron-transfer energies. When the reactants approach each other, the electron densities of the reactants are polarized and the number of electrons in the reactants changes; the terms on the last line correct for these effects. In this approach, dispersion is treated as a facet of polarization.

The electrostatic interactions are readily evaluated in any quantum chemistry program, so the only “difficult” terms are the electron-transfer terms on the last line of Eq. 11. These could be approximated using a quadratic expansion in the electron density, i.e., (Liu and Parr, 1997; Ayers and Parr, 2000)

The second-order derivative is the hardness kernel, (Berkowitz and Parr, 1988)

The chemical potential and the hardness kernel of the isolated reagent could be evaluated exactly from the output of a Kohn-Sham program or using one of the many approximate forms that have been proposed in the literature (Liu et al., 1997; De Proft et al., 1998; Ayers, 2001; Chattaraj and Maiti, 2001; Torrent-Sucarrat et al., 2005; Torrent-Sucarrat et al., 2007; Gómez et al., 2021). Inserting (14) into (11) gives

In order to obtain a simple and intuitive form in Eq. 16, we have defined the electrostatic potential of the reactants in the usual way, e.g., with

The actual interaction energy is then approximated by minimizing the expression in Eq. 16 with respect to all changes in density that preserve the total number of electrons and which do not cause the electron density of either fragment to become negative. i.e., subject to

Equation 19 extends the variational principle for the Fukui function to include effects other than electron transfer (Chattaraj et al., 1995; Ayers and Parr, 2000).

For computational purposes, the preceding approach is probably the most useful. For conceptual purposes, however, an alternative perspective that makes references to the most commonly used reactivity indicators may be preferred. Notice that many of the commonly employed reactivity indicators appear in the expression for the change in electron density. e.g.,

Note that we have defined the change in the electron density of A due to polarization as .

Equation 22 is directly useful only if we know . However, we do not know how the external potential of A changes due to the presence of the other reagent. (There is an exact formulation to compute , but it is computationally demanding (Ayers, 2000; Ayers and Parr, 2001; Kiewisch et al., 2008; Roncero et al., 2008; Fux et al., 2010)) Among the various approximations that have been proposed, the simplest is to identify with the electrostatic potential of molecule B (Ayers and Parr, 2001). (This is certainly acceptable when the reactants are reasonably far apart, so that the electrons in reactant A and reactant B can be considered distinguishable particles.) This simple approximate form may also be rationalized from Eq. 16, where the electrostatic potential of B serves as an “effective external potential” for the electrons in A.

Using and Eq. 22, we can expand the electrostatic terms in Eq. 11, giving,

In addition to the contributions to the interaction energy from electrostatic interactions of the reactants and the electron transfer between them, there are contributions due to mutual polarization of the reactants. The leading order term of this type is

This term is negative: the electrostatic potential of reactant B polarizes the electron density of A in a way that increases the attraction between the reagents.

Some readers may question whether should be multiplied by a factor of 1/2. After all, the change in energy due to polarization is

When the two reactants approach each other, they polarize each other. This polarization increases the attraction of the reactants for each other. However, polarization also mixes excited-state wavefunctions with the ground-state wavefunction; this increases the electronic energy of the reactants(Ayers, 2007a). The total change in energy, cf. Eq. 11, arises because the electronic energy of the reactants, and , increases by while the interaction energy between the reactants decreases by twice that amount. This explains why there is no factor of 1/2 in Eq. 24 and draws our attention to the final contributions to the fundamental expression for the interaction energy, Eq. 11. When the electron density of the reactants changes (due to polarization or electron transfer), the energy of the reactants also changes. One has

The terms in the first line represent the change in system energy due to electron transfer; the terms in the second line shows how the energy of the system increases as the density is polarized away from the ground-state density.

Combining Eqs 11, 22, 24, 27 gives an alternative method for computing the interaction energy between reagents.

A. Determine the non-additive kinetic energy and the non-additive exchange-correlation energy. Sometimes simple Thomas-Fermi-Dirac-like functionals will suffice here. One has

B. Compute the electrostatic potential of the isolated reactants. E.g.,

Use the electrostatic potential to estimate how the density of each reactant is polarized by the electrostatic potential of the other reactant,

Also compute the polarization energy,

C. Determine the extent of electron transfer, , by minimizing the expression for the interaction energy with respect to .

D. The electron density of the reagents is now changed due to electron transfer and mutual polarization. Use the updated electron densities, to compute an updated electrostatic potential . This can then be used to update the density change due to polarization and the polarization energy (using Eqs 31, 32; step B). The extent of electron transfer is now computed by substituting the updated quantities into (Eq. 33) and minimizing the energy with respect to (step C). The electron density of the reactants can then be updated and steps B and C can be iterated until value of the interaction energy has converged to the desired level of precision. (Notice how the supra-indices in Eqs 34, 35 distinguish between the magnitudes that are updated, from those that are kept constant.)

Step A could also be included in the iteration cycle, but if that is done then one needs to use the linear response kernel for the perturbed molecular fragment when evaluating Eq. 31; this may not be practical.

In this form, we are always expanding the Taylor series given in Eq. 33 around a fixed point that corresponds to the isolated reagents. We can improve this procedure by centering the expansion at a point that already contains the initial interaction between the molecules or, even more generally, the influence of the molecular environement (e.g., solvent) (Alain Miranda-Quintana R. and Ayers P. W., 2016, Alain Miranda-Quintana R. and Ayers, P. W. 2016; Alain Miranda-Quintana et al., 2016). For example, we can correct the energies using: where the sub-index M indicates the number of electrons, and the supra-index 0 identifies the isolated-system properties. We can also update the density by considering the effect of the perturbation on the Kohn-Sham orbitals:

In these cases we take the perturbation, , according to Eqs 23, 30. The modified energies and densities can then be used to estimate the reactivity descriptors present in Eq. 33.

An even more realistic solution will be to change the point around which we will expand the Taylor expansion after each iteration step. To do so, we will have to recalculate the perturbed energies and orbitals (e.g., densities) after each iteration, using the electrostatic potentials associated to Eq. 35. But, additionally, we must correct the expressions used to calculate the chemical potentials and Fukui functions to accomodate fractionally-charged systems:

Extending the preceding mathematical framework to include reagents with degenerate ground states is somewhat difficult, but is not impossible (Cardenas et al., 2011a; Bultinck et al., 2013, 2014; Bultinck et al., 2015). Extending this procedure to include open-shell reagents and spin-specific reaction pathways is quite straightforward: (Galvan et al., 1988; Galvan and Vargas, 1992; Ghanty and Ghosh, 1994; Ghosh, 1994; Ghanty and Ghosh, 1996; Vargas and Galvan, 1996; Vargas et al., 1998; Melin et al., 2003; Miranda-Quintana R. A. and Ayers P. W., 2016): none of the formulas change if the matrix-vector notation for spin-resolved reactivity indicator functions is invoked (Garza et al., 2006; Perez et al., 2008).

3 Discussion

3.1 Connections to Other Theoretical Methods

This method of computing the interaction energy from its fragment decomposition is fundamentally related to the method of density-functional embedding (Cortona, 1991; Vaidehi et al., 1992; Wesolowski and Warshel, 1993; Govind et al., 1999; Wesolowski, 2004; Wesolowski and Leszczynski, 2006; Elliott et al., 2010; Fux et al., 2010; Goodpaster et al., 2010; Goodpaster et al., 2011; Laricchia et al., 2011). That method, unlike this one, performs iterative ab initio calculations on the fragments rather than using the perturbation expansion.

This method is also fundamentally related to electronegativity equalization approaches to molecular mechanics (Mortier et al., 1985; Mortier et al., 1986; Yang and Mortier, 1986; Mortier, 1987; Rappe and Goddard, 1991; Bultinck et al., 2002a; Bultinck et al., 2002b; Verstraelen et al., 2009; Verstraelen et al., 2011). In those methods one replaces the non-additive kinetic and exchange-correlation energies with an empirically parameterized “repulsive wall.” One accelerates the method by approximating the Coulomb-like integrals by using atom-centered charges and (sometimes) dipoles, with the charges and dipoles varying in response to external fields (polarization) and electron transfer to/from the molecule; this amounts to using a coarse-grained version of Eq. 16. The result is a relatively fast and accurate method for computing molecular interaction energies. In this latter context it bears mention that both of the fundamental expressions for the interaction energy, Eqs 16, 33, can be adapted to cases where more than two molecules need to be considered (e.g., two reactants and a catalyst, several reactant molecules and the surrounding solvent). Assuming that three-body forces are negligible, the interaction energy in these more complicated situations is merely the sum over all the pairwise interaction energies.

This approach is fundamentally related to the electron-density-based energy decomposition analysis (DFT-EDA) (Wu et al., 2009). In DFT-EDA, two fragments are brought together with frozen density (corresponding to the first line in Eqs 16, 33) and then the exact (nonperturbative) electron-transfer energy and polarization-energies are defined and computed using constrained search computations. Like the approach in this paper, DFT-EDA is intended to elucidate chemistry, not to provide numbers. This approach, however, is mainly intended to predict chemistry in the perturbative regime, while DFT-EDA is generally employed for understanding known chemistry using an exact (nonperturbative) expression for the energy.

Finally, this approach is partly justified by the nearsightedness principle in density-functional theory (Li et al., 1993; Kohn, 1996; Baer and Head-Gordon, 1997; Prodan and Kohn, 2005; Prodan, 2006). It is safe to assume that when two reactants are far from each other, the “long-range” effects are mainly electron transfer and electrostatics, and that large-scale induced electron rearrangements are minimal (Fias et al., 2017). This means, in particular, that the estimate of the polarization density in step B, which implicitly assumes that the response of one subsystem to the other is small, and that the (total system) linear response kernel, , and higher-order responses is negligible (except perhaps for features described by electron transfer and the Fukui function/dual descriptor) when and are located on different subsystems (Cardenas et al., 2009b; Fias et al., 2017). That is, in step (B) we assume that the response kernel is nearsighted enough for the approximation to be valid. A similar nearsightedness argument applies to the hardness kernel in the alternative expression for the interaction energy in Eq. 16.

3.2 Connections to Previous Results in Density Functional Theory for Chemical Reactivity

The two expressions for the interaction energy, Eqs 16, 33, are the key result of this paper. In certain limits, these expressions recover known results from the DFT-CR literature. For example, suppose that the reactivity trends in fragment A do not depend on the detailed properties of fragment B. Then we can neglect the terms in Eq. 16 that depend on fragment B, obtaining,

Minimizing this expression with respect to all with ΔN electrons is equivalent to the variational principle for the Fukui function (Chattaraj et al., 1995; Ayers and Parr, 2000). This establishes the Fukui function as the key nonspecific reactivity indicator for electron-transfer reactivity. Recall that the lowest-curvature normal modes on a potential energy surface correspond to the most facile molecular rearrangements. From Eq. 40 we can approximate the elements of the Hessian matrix as,

This suggests that eigenvectors of the hardness kernel with small eigenvalues indicate favorable reactive modes. Or, to phrase this in the more conventional way: the eigenvectors of the softness kernel with the highest eigenvalues indicate favorable molecular rearrangements (Nalewajski and Koninski, 1987; Cohen et al., 1995; Cardenas et al., 2006). These simple guiding principles based on Eq. 16 obviously fail when semiquantitative results are needed or when detailed “matching” of reactants plays a pivotal role in the chemical reaction.

Equation 33 may seem very intimidating, but its interpretation is relatively simple. There are strong links to the previous two-reagent expressions in the literature, especially from the work of Parr and Pearson (for electron transfer), (Parr and Pearson, 1983) Berkowitz (for Frontier-orbital interactions), (Berkowitz, 1987) Anderson et al. (for the balance between electrostatics and electron-transfer terms), (Anderson et al., 2007a) and Morell et al. (for electrophile/nucleophile matching) (Morell et al., 2005; Ayers et al., 2007; Morell et al., 2007; Morell et al., 2008; Cardenas et al., 2009b; Labet et al., 2009). For example, if one assumes that all regioselective interactions may be neglected, then the entire expression (33) reduces to the conventional Parr-Pearson quadratic expansion for the electron-transfer energy, (Parr and Pearson, 1983; Ayers et al., 2006)

If one furthermore chooses one of the reagents to be a perfect electron donor, the change in energy is equal to the electrophilicity (Parr et al., 1999; Chattaraj et al., 2003; Chattaraj et al., 2006; Liu, 2009b).

If electron transfer is dominant and one wishes to consider regioselectivity, then the leading-order term is the Coulomb attraction between the Fukui functions,

This term plays a decisive role in the theory of Frontier-molecular orbital interactions put forth by Berkowitz (Berkowitz, 1987) and it justifies the use of the Fukui potential as a fundamental reactivity indicator (Cardenas et al., 2011b; Cardenas, 2011; Osorio et al., 2011; Yañez et al., 2021). When only the potential is varied, but not the number of electrons, one obtains potential-based reactivity indicators (Ayers and Parr, 2001; Cedillo et al., 2007; Ayers P. W. et al., 2009; Liu et al., 2009; Muñoz et al., 2020).

3.3 Term-By-Term Interpretation of Eq. 33

The Taylor series expansion in Eq. 1 converges subject to reasonable assumptions (Ayers et al., 2005a). Equation 33 is therefore approximate only in two ways: 1) higher order terms in the Taylor series are neglected and 2) only electrostatic effects are included in the model for the external potential change. Both assumptions are accurate for distant, weakly-interacting, reagents with nondegenerate ground states; they will be useful in all cases where the initial approach of the reagents guides the reactivity preferences (e.g., when the transition state is early). Both assumptions can be relaxed. Many higher-order reactivity indicators are known, (Fuentealba and Parr, 1991; Senet, 1996; Geerlings and De Proft, 2008; Cardenas et al., 2009a; Heidar-Zadeh et al., 2016b) though they are difficult to compute and inconvenient to work with. The exact “effective external potential” that a molecule feels in the presence of a reagent can be computed, and various levels of approximation are known (Ayers and Parr, 2001). Some of these approximations are not that much more complicated than the simple electrostatic approximation. (For example, including an “exchange-correlation charge density” is relatively straightforward (Liu et al., 1999; Ayers and Levy, 2001; Ayers and Parr, 2001; Andrade and Aspuru-Guzik, 2011)) Even with the approximations implicit in Eq. 33, however, the main chemical effects that drive reactions are included.

3.3.1 Steric Repulsion and Electron Pairing

The non-additive kinetic and exchange-correlation energy terms, capture the effects of electon pairing (from combining the electron densities of the fragments) and, when closed shells are pushed together, capture the “Pauli” and “kinetic” portions of the repulsive wall in potential energy curves. The electrostatic repulsion between the fragments also contributes to steric hindrance; this contribution is included in the second line.

We note that this measure of the steric repulsion is different, both physically and mathematically, from the quantification of the steric effect proposed by one of us, (Liu, 2007) and then studied and extended by us and others (Liu, 2007; Nagy, 2007; Liu and Govind, 2008; Torrent-Sucarrat et al., 2009; Ugur et al., 2009; Ess et al., 2010; Liu et al., 2010; Tsirelson et al., 2010; Esquivel et al., 2011; Huang et al., 2011; Fang et al., 2014). Nonetheless, the two measures have some mathematical similarities (in both cases the dominant term is associated with the kinetic energy) and are expected to give similar insights. Note also that Eq. 44 also includes some long-range contributions associated with exchange and correlation and short-ranged contributions from electron pairing that are not true “steric” effects.

3.3.2 Electrostatic Interactions

The electrostatic attractions/repulsions between the isolated fragments are accounted for in the term

Notice that, Eqs 44, 45 together can already describe weak interactions; they are the key elements in the Gordon-Kim theory, for example (Gordon and Kim, 1972).

3.3.3 Polarization and Dispersion

The term accounts for the electrostatic stabilization that occurs because the fragments are polarized. (The very last line in Eq. 33 corrects for the “promotion energy” that is required to deform the electron density to make these favorable electrostatic interactions.) The second-order change in electrostatics due to the mutual polarization of the fragments is captured in the term

If the interaction-energy expression is iteratively converged (step D in Section 2), this term includes the “induced multipole-induced multipole” interactions that are implicated in dispersion.

3.3.4 Electron Transfer; Fukui Function

The Coulomb interaction between the Fukui functions (Parr and Yang, 1984; Yang et al., 1984; Ayers and Levy, 2000; Ayers P. W. et al., 2009; Osorio et al., 2011; Yañez et al., 2021) is the leading-order purely electron-transfer term,

The importance of this term had already been noted by Berkowitz (Berkowitz, 1987). Electron transfer will be favorable when the Coulomb repulsion between the donor Fukui function and the acceptor Fukui function is large; this is associated with strong overlap of the Frontier orbitals. Alternatively, Eq. 48 can be viewed as the attraction between the donor orbital and the acceptor orbital, or between the electron and the hole (Ayers and Levy, 2000). Similar terms are found to play a role in Frontier orbital theories, especially when expansions around the separated atom limit are considered (De Proft et al., 2006; Cardenas et al., 2009b).

Electron transfer also changes the electrostatic interactions between reagents at both the isolated reagent, and polarized reagent levels. Terms similar to Eqs 48, 49 can be important even when electron transfer is small because the asymptotic behavior of the linear response function is governed by the Fukui function (Cardenas et al., 2009b; Fias et al., 2017).

3.3.5 Electron Transfer; Dual Descriptor

The dual descriptor (Morell et al., 2005, 2006; Ayers et al., 2007; Morell et al., 2007; Geerlings et al., 2012) tends to play a role when the reagents are electrostatically flat and the Fukui functions do not establish a strong reactivity preference. The leading-order purely dual-descriptor term is

This expression has appeared previously in the literature, but without derivation (Ayers et al., 2007). Notice that, unlike the interaction between the Fukui functions (Eq. 48), it is favorable to have a large electrostatic attraction between the dual descriptors of the fragments; this corresponds to aligning the nucleophilic portions of one fragment with the electrophilic portions of a different fragment, and vice versa (Morell et al., 2005, 2006; Ayers et al., 2007; Morell et al., 2007; Geerlings et al., 2012).

The change in electron density due to second-order electron transfer also causes a shift in the electrostatic interactions with the isolated (compare Eq. 49) and polarized (compare Eq. 50) reagents. There is also a coupling to the first-order electron transfer (Fukui function) term,

If any of these terms are large, it indicates that the electrostatic potential, polarization, or first-order electron-transfer (Fukui function) terms are likely to be dominant. The dual descriptor is also important when describing the asymptotic behavior of the quadratic and cubic density response functions; it has been argued that the dual descriptor’s chemical relevance may be mostly attributed to those nonlinear polarization terms (Cardenas et al., 2009b).

3.3.6 Energy Change From Electron Transfer; Chemical Potential and Hardness

The remaining terms in Eq. 33 represent the change in energy of the isolated reagents. First of all, the isolated reagents change in energy because the number of electrons in the reagents changes. This leads to the term

The electronic chemical potential is minus one times the electronegativity, (Parr et al., 1978) and measures the intrinsic strength of a Lewis acid or base (Ayers et al., 2006). The chemical hardness is the key to describing the HSAB principle (Parr and Pearson, 1983). The electron-transfer energy expression in Eq. 55 suffices, all by itself, to derive the HSAB principle (Chattaraj et al., 1991; Ayers, 2005; Chattaraj and Ayers, 2005; Ayers and Cardenas, 2013; Cardenas and Ayers, 2013).

3.3.7 Promotion Energy

Polarizing the reagents so that they can form more favorable electrostatic interactions requires mixing excited-state character into the ground-state fragment densities. This raises the electronic energy of the fragments by, though this term is doubly-cancelled out (cf. Eq. 46) by the favorable electrostatic contacts that are created (Ayers, 2007a).

3.3.8 Condensed Reactivity Indicators

It is often convenient to replace the pointwise interactions between reactants in the interaction energy model with interactions between the atoms in the reactants. This gives an atom-condensed version of the reaction energy expression (Yang and Mortier, 1986) which strongly resembles the structure of electronegativity equalization molecular mechanics force fields. Letting denote the atomic weight function associated with the α th atom in reactant A, and let and denote its nuclear charge and position (consistent with the notation in Eq. 3). The expressions for the atom-condensed charges, condensed Fukui function, (Fuentealba et al., 2016) condensed dual descriptor, and condensed polarization density are, respectively

Here we have elected to use the fragment-of-molecular response method, though the response-of-molecular-fragment method could be used instead (Fuentealba et al., 2000; Ayers et al., 2002; Bultinck et al., 2007; Padmanabhan et al., 2007; Zielinski et al., 2012; Miranda-Quintana R. A., 2016). We could also include dipole (and higher-order multipole) components of the descriptors in Eq. 57, but this would complicate the formula for the polarization energy and the interaction energy

This equation is more convenient as a “working formula” than Eq. 33 because it removes the need to perform six-dimensional integrals. We note, however, that once one chooses a suitable energy model, (Parr and Bartolotti, 1982; Parr and Pearson, 1983; Noorizadeh and Shakerzadeh, 2008; Fuentealba and Cardenas, 2013; Heidar-Zadeh et al., 2016a; Miranda-Quintana R. A. and Ayers P. W., 2016) all the integrations in Eq. 33 can be performed analytically by combining the basis-set expansion of one-electron density matrices and the linear response function with the appropriate one- and two-electron integrals.

4 Summary

In the density functional theory approach to chemical reactivity (DFT-CR), it is traditional to focus on the reactivity indicators of a single reagent. This is not always adequate, as interactions between reagents can be decisive when, for example, partner-specific interactions are critically important for determining reactivity preferences. In this work, we have derived mathematical formulas for the interaction energy between reagents in terms of chemical reactivity indicators (Eqs. 16, 33) and shown how these expressions are related to existing computational approaches (density-functional embedding, density-functional energy decomposition analysis) and previous reactivity rules (Section 3.2). We also showed that all of the commonly accepted chemical interactions that influence chemical reactivity are contained in the key Eq. 33 (Section 3.3). We believe that this work provides not only a basis for future work in DFT-CR, but that we can rationalize existing DFT-CR success stories by noting that whenever one term in Eq. 33 is decisive, the conventional one-reagent-one-reactivity-indicator approach is likely to be successful. In further contributions in this series we will explore how the two-reagent picture can be used to provide further analytical arguments in favor of several reactivity principles.