Electronic Density Approaches to the Energetics of Noncovalent Interactions

We present an overview of procedures that have been developed to compute several energetic quantities associated with noncovalent interactions. These formulations involve numerical integration over appropriate electronic densities. Our focus is upon the electrostatic interaction between two unperturbed molecules, the effect of the polarization of each charge distribution by the other, and the total energy of interaction. The expression for the latter is based upon the Hellmann-Feynman theorem. Applications to a number of systems are discussed; among them are dimers of uracil and interacting pairs of molecules in the crystal lattice of the energetic compound RDX.


Introduction
Noncovalent interactions are ubiquitous: between enzymes and substrates, in hydrogen bonding, physical adsorption, solvation, condensation processes, etc. Calculating the energy associated with such an interaction is, in principle, straightforward; if a system M is formed from N components M i , then the stabilization energy of M is, where E M and are the respective equilibrium, ground-state energies.Eq. ( 1) is exact.However it suffers from the fact that ∆E i M E stab is typically several orders of magnitude smaller than E M and the ; thus, barring fortuitous cancellation, any errors in these quantities will be greatly magnified in This problem can of course be minimized by computing E M and the at high levels of accuracy, but this is likely to be prohibitively expensive in terms of processing resources for many systems of practical interest.
i M E Another difficulty with eq. ( 1) is the so-called basis set superposition error (BSSE).This refers to the size-imbalance between the basis sets used for M and the M i , which are smaller.The result is an artificial stabilization of M [1][2][3].The effect diminishes with the use of larger basis sets, but this solution again involves increased computational cost.Some time ago, Boys and Bernardi suggested addressing BSSE by introducing "ghost" orbitals in treating the M i [4]; this is known as the counterpoise procedure.It is now widely used, although there has been considerable controversy in the past concerning its effectiveness [1,3].
The use of eq. ( 1) in connection with noncovalent interactions is sometimes denoted the ab initio or supermolecular approach.For more extensive discussions, see Chalasinski and Szczesniak [5] and Rappe and Bernstein [6].
The problem of achieving sufficient accuracy with eq. ( 1) has made perturbation theory an attractive alternative [7][8][9].This directly produces the interaction energy E int between the components M i , without the need to take differences between large numbers.However the M i are normally assumed to retain their ground-state geometries; no account is taken of any changes in these that may accompany the interaction [5,7].In contrast, ∆E stab , eq. ( 1), is obtained using the optimized structures of M and the M i .Thus ∆E stab and E int differ by the energy involved in any geometry changes that occur.This may be quite small, however, as shall be shown later.
In the perturbation theory formulation of E int , it is expressed as the sum of a series of terms.This is frequently considered to be an advantage, since these can be assigned physical interpretations.For example, E int is often viewed as composed of four elements: 1) the electrostatic interaction between the unperturbed M i ; 2) their mutual polarization of each other's charge distributions; 3) dispersion effects, involving intermolecular electronic correlation; and 4) exchange/repulsion, reflecting the overlapping of electronic distributions [8,10].E int is thus written as, Higher levels of theory yield more elaborate representations of E int [5,7,9].The various contributions to E int , as in eq. ( 2), can all be formulated in terms of the perturbation operator and the wave functions of the M i [7][8][9].However simplified expressions are frequently used [8,[10][11][12], for instance in molecular dynamics simulations [8,11].In the latter, a point-charge approximation is commonly employed for the electrostatic interaction, E es , and a Lennard-Jones or Buckingham-type potential for E disp + E ex-rep .E pol is not taken into account, although it could be done, for example by periodically changing the magnitudes of the point charges in the course of the simulation [13].
In the remainder of this paper, we shall focus upon E es , E pol and E int .For convenience in notation, we shall treat the case of two components A and B forming a noncovalently-bound complex AB.Extension to more components, on a pair-by-pair basis, is straightforward.

Methodology
The energy E es of the electrostatic interaction between the unperturbed components A and B is given rigorously in terms of their electronic densities ρ A (r) and ρ B (r) by, In eq. ( 3), Z M,A and Z N,B are the charges on nuclei M and N of components A and B; R M and R N are their locations.
Since the charge distributions of A and B do not remain unperturbed as the interaction proceeds, but rather have a polarizing effect upon each other, it is necessary to include the associated energy, E pol , in E int .Several ways of determining E pol have been proposed [7,8,12,14].We have recently introduced another approach [15], which involves expressing the electronic densities of A and B after interaction as ρ A * (r) = ρ A (r) + ∆ρ A (r) and ρ B * (r) = ρ B (r) + ∆ρ B (r) ; ∆ρ A (r) and ∆ρ B (r) are the changes due to polarization.Replacing ρ A (r) and ρ B (r) in eq. ( 3) by ρ A * (r) and ρ B * (r) produces E es + E pol ; subtracting eq. ( 3) from both sides leaves, We represent ∆ρ A (r) and ∆ρ B (r) by partitioning them into overlapping and nonoverlapping portions [15], which are then treated separately.They are obtained from the electronic density of the complex, ρ AB (r) .
Eqs. ( 3) and ( 4) give E es and E pol in terms of integrals over electronic densities.We evaluate these numerically [15], using a technique modeled after that of Gavezzotti [16].The electronic charge distributions in the integrands are divided, by means of three-dimensional grids [17], into large numbers of uniform volume units, "e-voxels."Each of these has a negative charge with magnitude equal to its volume times the electonic density at its center.We discard those e-voxels that are outside of the assigned boundary surface of the charge distribution, which is defined to be a specific value of its electronic density, ρ min .To facilitate the computations, blocks of n x n x n e-voxels are "condensed" into "super-e-voxels," having charges equal to the sums of those of their constituents.The integrals in eqs.( 3) and ( 4) are then evaluated by calculating the electrostatic interactions of the super e-voxels of each molecule with those of the other, or with its nuclei.More detailed discussions of this integration procedure can be found in Gavezzotti [16] and in Ma and Politzer [15].

Applications
We used the methodology that has been outlined to determine E es for several molecular dimers [15]: (H 2 O) 2 , (CH 3 OH) 2 , (CH 2 Cl 2 ) 2 , (CH 3 CN) 2 , (CH 3 COCH 3 ) 2 , (CH 3 SOCH 3 ) 2 .This was done primarily at several Hartree-Fock levels.One of our objectives was to ascertain the number of e-voxels, the condensation number n, and the ρ min that would be most effective.We found that E es converges for approximately 10 6 e-voxels and, for the smaller systems, when ρ min ≤ 10 −5 au (electrons/bohr 3 ).It was also our experience that at least 2000 super e-voxels are needed; thus, for 1.0 x 10 6 e-voxels, n should be no larger than 7. Our E es agree well with those obtained by the Morokuma-Kitaura scheme for partitioning interaction energies [18,19].E pol was computed for the water dimer [15], for which we could compare the results with those from the Morokuma-Kitaura and also the reduced variational space self-consistent-field methods [20].Hartree-Fock electronic densities were used, corresponding to ten different basis set combinations.The three sets of E pol were in good accord when ρ min was taken to be 0.01 au.It is not surprising that the optimum ρ min is not the same for E pol as for E es ; the extent of overlap between the components, which is determined by ρ min , influences E pol more directly than E es .E pol was observed to have a relatively low sensitivity to basis set, less than that of E es [15].Among the ten (H 2 O) 2 calculations, the largest difference in E pol was 0.41 kcal/mole, between HF/6-31G(d,p)//6-31G(d,p) and HF/cc-pVDZ//6-31G(d,p).
We have also extended these studies to some larger systems, the first of which was the dimer of uracil (1).Since this involves bigger molecules than any of the other dimers for which we computed E es , we tested whether ρ min ≤ 10 −5 au is still sufficient for E es to converge [15].Two dimer structures were investigated, which differ in the relative orientations of the uracil molecules: face-toface and face-to-back [21].We found that ρ min ≤ 10 −6 au is now required; this can be seen in Table I.The fact that the magnitude of E es for the face-to-face dimer is more than double that for the faceto-back was attributed to the former having more intermolecular N−H---O hydrogen bonds [15].
Finally, we investigated the intermolecular interactions in the crystal lattice of RDX (2, hexahydro-1,3,5-trinitro-s-triazine) [22], which is of considerable interest as an important energetic compound.RDX is frequently the subject of molecular dynamics simulations [11,23], which generally obtain E es by a point-charge approximation.One of our objectives was to examine how this compares with our E es from eq. ( 3).We considered (a) the interaction within an interlocked pair of molecules, and (b) that between two molecules in neighboring interlocked pairs, at the Hartree-Fock and B3PW91 levels, with 6-311+G** basis sets.The electrostatic interaction energies E es were calculated by two pointcharge models, involving (a) Mulliken and (b) CHelpG atomic charges [17]; the latter are derived from electrostatic potentials.We also determined E es from the electronic densities by means of eq. ( 3), using 1.4 x 10 6 e-voxels and ρ min = 1.0 x 10 −6 au.The Mulliken charges produced very poor E es , positive for both pairs of interacting molecules.The CHelpG were negative but significantly smaller in magnitude than the E es from eq. ( 3) (e.g.−8 vs. −3 kcal/mole for the interlocked pair).We conclude that these point-charge models do not satisfactorily reproduce E es .E pol was also computed for the two RDX pairs, with eq. ( 4); it was in the neighborhood of −1 kcal/mole in each case.

Methodology
We shall not approach E int in terms of eq. ( 2), but rather from the standpoint of the Hellmann-Feynman theorem [24][25][26].In its general form, this states that, for a system described by, In eq. ( 6), λ is any parameter appearing in the Hamiltonian ˆ H . Thus, the theorem can be expressed in various ways, depending upon the choice of λ.For example, letting λ = Z, the nuclear charge of an N-electron atom with electronic density ρ(r), produces [27,28], E Z N (r) r dr V 0 (7) in which V 0 represents the electrostatic potential at the nucleus due to the electrons.This leads to exact relationships between E and V 0 [27,29,30], including, Analogous equations can be derived for molecules [29][30][31].
For our present purpose, we take λ to be R M , the position of nucleus M in the system of interest.

Since
gives the force M exerted upon M by the electrons and other nuclei, then eq. ( 6) becomes, where ρ(r) is the electronic density of the system.Eq. (11) shows that the force upon any nucleus is given by classical electrostatics.
The binding energy of a nucleus can in principle be determined by using eq.( 11) to calculate the work done in bringing it from infinity to its equilibrium position in the force field of the remainder of the system [24,[32][33][34].Extending this approach, we have recently formulated the stabilization energy of a noncovalently-bound complex AB as the work done upon the nuclei and electrons of component A as it is brought from infinity to its position in AB in the force fields of the nuclei and electrons of B [35].A complicating factor, which was pointed out by Bader in a different context [34], is that the electronic densities of A and B change somewhat as they approach.We neglect this, and use the geometries and electronic densities of A and B as they are in AB.In view of this, the quantity that we obtain is the interaction energy of A and B as they are in the complex, and shall be designated .
It is given by [35], In deriving eq. ( 12), it is assumed that the components A and B retain their identities in AB.

Conceptually, differs from E
* int E int , eq. ( 2), in that the latter is obtained using the ground-state geometries of isolated A and B. Since these often remain essentially the same during the formation of AB, it can be anticipated that , E * int E int and ∆E stab will frequently be quite similar, if evaluated at comparable levels of accuracy.Indeed, the total energies required to transform the component molecules in (H 2 O) 2 and (HF) 2 from their ground states to their forms in the dimers were found to be 0.09 kcal/mole [36] and 0.03 kcal/mole [37], respectively.For the face-to-face dimer of uracil [21], a larger molecule, this energy is 0.79 kcal/mole at the MP2/6-31+G* level [35].It can be significantly greater, however, for ion-molecule interactions, e.g.F − (H 2 O) [36].
Eq. ( 12) shows, in the spirit of Feynman [25], that the total interaction energy is classically electrostatic.Eq. ( 12) is in fact formally identical to eq. ( 3), differing only in that it involves the electronic densities of A and B as they are in the complex rather than in the free states.Thus the quantities E pol , E disp and E ex-rep in eq. ( 2) are simply compensating for E es not being in terms of the appropriate electronic densities.If it were, then E es alone would suffice in eq.(2).
A practical concern with regard to eq. ( 12) is partitioning the total electronic density ρ AB (r) into ρ A (r) and ρ B (r).To do this, we first establish a boundary surface ρ min for the complex; this may differ from those used to obtain E es and E pol .At each point r within this surface, we determine the ratios of its distance from each nucleus divided by the van der Waals radius of that atom.The point r and the corresponding ρ AB (r) are then assigned to the atom (and therefore the component A or B) for which this ratio has the lowest value.We perform the integrations in eq. ( 12) by means of the numerical technique described in an earlier section of this paper.

Applications
The number of e-voxels that we use in computing depends upon the size of the molecules that are involved, but it continues to be of the order of 10 * int E 6 .Thus, it was 1.0 x 10 6 for (H 2 O) 2 [35], but 3 x 10 6 for the pair interactions in the crystal lattice of RDX (2) [22].With regard to ρ min , we found that converges for ρ * int E min ≤ 10 −4 au [35].
Our initial calculations of were for four molecular dimers for which reasonable computational/experimental estimates of the stabilization energies ∆E * int E stab are available in the literature, to which our results could be compared.The systems studied included (H 2 O) 2 , (HF) 2 , (H 3 COH) 2 and (HCOOH) 2 [35]; the calculations were carried out with the Hartree-Fock, MP2, B3LYP and B3PW91 procedures, and three different correlation-consistent basis sets.
For (H 2 O) 2 and (HF)  [38].Several factors may contribute to this (besides the approximations in our procedure), one being a degree of uncertainty in the literature values.It is also likely that our computed dimer structures differ somewhat from the experimental ones upon which these ∆E stab are based.Finally, the calculated electronic densities are of course not exact.The Hartree-Fock were invariably more negative than the others, but the spread was less than 1 kcal/mole, except for (HCOOH) * int E 2 , for which it was about 3 kcal/mole.For a given computational method, the three basis sets usually gave quite similar results, particularly the two larger ones, cc-pVTZ and cc-pVQZ.
We also determined for the two pairs of molecules in the RDX crystal lattice for which, earlier in this paper, we discussed E * int

E
es and E pol .Our predicted for the interlocked pair was * int E about −8 kcal/mole, and −2 to −3 kcal/mole for the interaction between interlocked pairs [22].These values are very similar to what we obtained for the corresponding E es , approximately −8 and −3 kcal/mole.A similar situation was found by Bukowski et al [39] in a symmetry-adapted perturbation theory analysis of dimers of dimethylnitramine, (H 3 C) 2 N−NO 2 , a molecule with the same structural elements as RDX.They found E es and E int to differ by ≤ 1 kcal/mole for each of the three most stable dimer configurations.

Discussion and summary
The fact that E es is a good approximation to for two pairs of RDX molecules, i.e., the electrostatic interations between the separate components nearly match the total interaction energies, suggests that (barring fortuitous cancellation) * int E ρ A (r) and ρ B (r) , eq. ( 3), are similar to ρ A (r) and ρ B (r), eq. ( 12).This appears to be the case for the dimethylnitramine dimers as well, in view of the findings of Bukowski et al [39] (vide supra).On the other hand, for the two uracil dimers mentioned earlier, our computed E es differed in each instance by about 3 kcal/mole from the best estimate of ∆E stab [15], being more negative for the face-to-face dimer and more positive for the face-to-back.
The point-charge model was not successful in reproducing E es for the two pairs of RDX molecules, for the charge definitions investigated.We tested the possibility that the model might be more effective if the charges were obtained for each pair after interaction, perhaps yielding reasonable approximations to , but there was, in general, no improvement.

E
We believe that our results overall support the formulations in terms of electronic densities that have been given for E es , E pol and , and the validity of the numerical integration technique that is used to evaluate them.However there is certainly a need for continuing efforts to optimize the assignments of the parameters − i.e., number of e-voxels, ρ * int E min and condensation number n − taking into account the energy quantity being sought and the sizes and shapes of the molecules.The effects of various computational methods (e.g., Hartree-Fock, MP2, density functional) and basis sets should also be further explored.
2, there was overall very good agreement between and ∆E