One-and two-center energy components in the atoms in molecules theory

The energy decomposition scheme proposed in a recent paper has been realized by performing numerical integrations. The sample calculations carried out for some simple molecules show excellent agreement with the chemical picture of molecules, indicating that such an energy decomposition analysis can be useful from the point of view of connecting quantum mechanics with the genuine chemical concepts. © 2001 American Institute of Physics. @DOI: 10.1063/1.1381407 #


I. INTRODUCTION
Most recently, one of us demonstrated 1 that the Hartree-Fock ͑HF͒ energy of a molecule can be presented exactly as a sum of atomic and diatomic contributions in the framework of the Atoms in Molecules ͑AIM͒ theory, 2 introduced by Bader, or in any other scheme in which the threedimensional physical space is decomposed into disjunct atomic domains ͑basins͒.This result follows from the simple facts that ͑a͒ the integrals over the whole space are equal to a sum of integrals over the individual domains, ͑b͒ the nuclei and the atomic basins can usually be put into one-to-one correspondence with each other.
The HF energy contains both one-and two-electron terms, described by single and two-fold integrals over the 3D space, respectively.The kinetic energy integrals and the electrostatic interaction of the electronic charge within a given basin with its own nucleus contribute to the atomic ͑one-center͒ energy component, while the electrostatic interaction of the electronic charge in one domain with a nucleus in another one represents a diatomic effect.Each two-electron ͑Coulomb or exchange͒ integral over the molecular orbitals ͑MO's͒ will contain both monoatomic and diatomic components, depending on whether the two integration variables are actually in the same or different basins.The total molecular energy contains also the intermolecular repulsions, which are obviously of diatomic nature.
The authors of Ref. 1 have also demonstrated the close conceptual correspondence of such a decomposition in the AIM theory with the recent chemical energy decomposition analysis ͑CECA͒ performed in the linear combination of atomic molecules ͑LCAO͒ framework, 3,4 which seems to be a promising simple tool for the a posteriori chemical analysis of the results obtained in the ab initio calculations.The equations of these two schemes are connected with a mathematical mapping which, however, does not mean that the numerical values should coincide; moreover, the energy decomposition in the AIM case is, in principle, exact while in CECA it is only a ͑good͒ approximation.We consider the possibility of presenting the molecular energy as a sum of atomic and diatomic contributions to be a result of upmost conceptual importance from the point of view of connecting quantum mechanics with the genuine chemical concepts.In fact, such an energy decomposition can help us to obtain a deeper insight into the physical content which is behind the chemical structural formulas.Of course, the abstract mathematical results of Ref. 1 will appear truly useful from the practical chemical point of view only if the decomposition leads to numerical results which correlate well with the chemical picture of the molecules.The aim of the present work is to demonstrate that this is really the case.As the integration of molecular orbitals over the atomic basins cannot be performed-at least at time being-analytically, we have had to recur to numerical integrations.
It is worthwhile to point out that this decomposition is not restricted to the HF level of theory, as the expectation value of the energy can always be written in the form of sums in which the one-and two-electron integrals over the molecular orbitals are multiplied by the respective first or second order density matrix elements.
It should also be mentioned that the AIM theory itself contains an energy decomposition scheme which is based on the atomic virial theorem. 2 This scheme presents the molecular energy as a sum of atomic energies, ⌺ A E(A), integrated for each atomic basin.That type of decomposition is extremely important from a physical point of view, but it is lacking a direct connection with the chemical notion of atoms interacting with each other.This aspect of the original AIM energy decomposition motivated Sierraalta and Frenking to assign part of the atomic energies to the diatomic interactions, by using some overlap integrals over the atomic basins as proportionality factors. 5Their theory, however, always gives negative diatomic energy components, thus it is unable to distinguish between bonding and non-bonding-or more generally, between attractive and repulsiveinteractions between the individual pairs of atoms.

II. THEORY
As discussed in Ref. 1, the self-consistent field ͑SCF͒ energy can be strictly decomposed into one-center and two-center components: where N/2 is the number of doubly filled MO's and the ͓12͉12͔ convention is used for the two-electron integrals.In the one-electron integrals, the subscript ''A'' denotes that the integration is performed over the Ath atomic basin For the two-electron integrals, the subscript ''A,B'' indicates that the integration for the electrons 1 and 2 are carried out over the atomic domains ⍀ A and ⍀ B , respectively:

III. COMPUTATIONAL ASPECTS
The integrals defined over the atomic basins must be computed numerically.The identification of the individual AIM basins and the numerical calculation of the one-electron ͑kinetic energy and nuclear attraction͒ integrals over them have been performed by using a slightly modified version of the PROAIMV 6 program.The same program was utilized to generate some arrays of data, which were saved on disk and have been used to calculate the two-electron integrals as follows: The PROAIMV program performs the numerical integration over the atomic domains by generating a grid of points with the radius vectors r ជ and assigning them respec- tive weight factors w , so that the integral of some scalar function f (r ជ ) over the basin ⍀ A is approximated as Here and further on, the notation A indicates the set of points belonging to the basin of the atom A.
We have used such numerical integrations to calculate the kinetic energy and nuclear attraction integrals for each nucleus C and each occupied MO i , by defining f (r ជ ) ϭϪ 1 2 i (r ជ )ٌ 2 i (r ជ ) and f (r 2 for each atom C, respectively.Only real orbitals have actually been considered in the calculations.
The numerical calculation of the two-electron integrals can be performed in a straightforward manner by a repeated use of the integration rule given in Eq. ͑7͒.As a preparatory step to such calculations, we have generated by the properly modified program PROAIMV the set of data defined as for each atomic basin A, and each pair of occupied MO's i and j ͑including the case iϭ j) in every point r ជ , A.
Using these quantities, the repeated use of Eq. ͑7͒ for calculating the integrals of exchange type leads simply to A similar formula holds for the individual Coulomb integrals, too.However, one can compute the sum of the Coulomb type contributions at once, because For this case, the functions G i j A with iϭ j have been utilized.By introducing the auxiliary quantities the sum of the integrals of Coulomb-type over the pairs of atomic basins can be written as All the above formulas also hold in the case when basins ⍀ A and ⍀ B coincide.The numerical calculation of the two-electron integrals over the pairs of atomic basins according to Eqs. ͑8͒-͑13͒ has been performed by a small FORTRAN program written by us.It also handles the symmetry of the molecule in order to avoid repetitive calculation of identical quantities.This program has been used on different Linux machines in a single node regime and on a SGI Origin 2000 in a parallel way, with a very good scaling with the number of nodes.Moreover, since we are only using the set of data generated by the PROAIMV program, the timing of the two-electron energy contribution depends only on the size of the grid, but is independent of the basis set used.
As given by the PROAIMV program, the total number of points is generated from three input parameters, the number of planes , planes and the number of radial points.A calculation using some standard integration parameters ͑64,48,96͒ involves more than 140.000 points per atomic basin for the H 2 molecule.That means that to compute the two-electron integrals ca.2ϫ140.000 2 Ϸ4•10 10 points would appear in the numerical integration, which makes the two-electron integrals extremely costly from a computational point of view.Therefore, we are forced to use the smallest grid of points for each atomic basin which still gives an acceptable accuracy.
It has been found that, in general, about 40.000 points are needed to get a good accuracy.Table I gathers the results of the integration for the model N 2 molecule with respect to the number of points used for each atomic basin.It can be seen that the one-electron part converges faster to the exact value than the two-electron one.Our integration shows, in general, similar accuracy to the sum of the AIM atomic energies.However, when a relatively small grid is used, a better accuracy is obtained.This may be due to the fact that AIM atomic energies are computed assuming the fulfilment of the atomic virial theorem within the atomic domains, so that a more accurate numerical integration is desirable.The energy partitioning we propose is only based on the assumption that the whole space is partitioned in domains, and each domain must contain one nucleus.

IV. ILLUSTRATIVE CALCULATIONS
We have performed test calculations for several diatomic, i.e., H 2 , N 2 , BH, and HF, and more complex molecular systems such as C 2 H 2 , C 2 H 4 , C 2 H 6 , and B 2 H 6 .We have used the 6 -31G(d,p) basis set whenever it was possible.It is to be mentioned that the AIM analysis sometimes yields basins with so-called non-nuclear attractors ͑NNA͒.These cases can be included in the above frame by assigning to the NNA a dummy nucleus with zero nuclear charge.However, it is not chemically very appealing that then some atomic and diatomic ͑bonding͒ energy will also be assigned to the dummy atoms.In such cases it may be worthwhile either to introduce another decomposition of the space or to check whether the appearance of the non-nuclear attractor is not an artifact of the basis set applied.For instance, for the acetylene, the 6 -31G(d, p) basis set and many others like 6 -311G(d,p) or even the 6 -311G(2d f ,2p) and 6 -311ϩϩG(3d f ,2pd) exhibited the undesirable nonnuclear attractors in the middle of the C-C bond.Hence, we have been forced to use the 6 -31G for this specific case.Nevertheless, the results of the integration for the 6 -31G(d,p) basis set using a dummy nucleus will be discussed in more detail later.Also, the molecule of diborane presented integration problems with the 6 -31G(d,p) basis set, due to a bad location of the bond critical points in the molecular plane.For this reason we have used the 6 -31G(d) basis set instead.The number of planes , planes , and the number of radial points of the atomic grid has been set to 32, 24, and 32, respectively for all the calculations described in Tables II-V.
Tables II, III, and IV collect the results obtained for several molecules.The accuracy of the integration is very good and comparable to the AIM integration.The dominating energy components are the large negative atomic ͑one-center͒ ones, while the diatomic ͑two-center͒ interaction energy components are smaller in absolute value-except for the H-atom energies in polar molecules such as HF or BH.The diatomic energy components are negative for chemically bonded atoms and can be of either sign for the nonbonded interactions.It is to be stressed that these energy components are static parameters corresponding to the given geometry and wave function, so the diatomic energy components are not directly related to the dissociation energies.The onecenter components are somewhat higher than the free atomic energies, reflecting the promotion of atoms during the bond formation.͑This energy is then regained in the form of bond energy.͒Accordingly, the static diatomic energy components are more negative than the respective dissociation energies, as the latter give the net energy effect with respect to the sum of free atomic energies.
The changes in the two-center energy components in the hydrocarbon series are in good agreement with the chemical intuition.The C-C energy monotonically increases ͑in absolute value͒ from ethane to ethylene and acethylene, and the opposite trend is observed, in turn, for the C-H energy.
The accuracy of the energy expansion which can be practically achieved by performing the numerical integrations in the present-conceptually exact-method is approximately the same as that which one gets in the recent approximate LCAO energy decomposition scheme 3,4 ͑CECA͒.The results of the two schemes agree qualitatively but not quantitatively: the AIM bonding energy components are usually less negative, and the one-center components are not as negative as observed for the CECA method.This behavior seems to be closer to the intuitive chemical picture than that of CECA, as CECA gives numbers which may be considered somewhat exaggerated.͑At the same time, CECA is a computationally very cheap method which can be routinely applied to large systems, too.͒ The only surprising result corresponds to the large repulsive interaction between the two boron atoms in the diborane molecule.A previous result obtained with the CECA method 3 produced attractive interactions between the boron atoms ͑Ϫ0.227 a.u.͒ and also for the B-H interactions ͑Ϫ0.279 and Ϫ0.517 a.u.for the B-H bridge and B-H terminal, respec-  tively͒.In the present method, the B-H energies are much more negative ͑Ϫ0.746 and Ϫ0.829, respectively͒ and compensate for the repulsive B-B and H-H energy components.Finally, the results for the acethylene molecule exhibiting a non-nuclear attractor in the middle of the C-C bond are collected in Table V.The energy of the dummy atom X corresponding to the NNA is positive.This is due to the fact that the kinetic energy is not compensated by any electronnuclear interaction, because there is no nucleus to be assigned to the NNA.The C-C interaction is also repulsive whereas the interaction between the carbon atom and the dummy atom at the NNA ͑C-X͒, is strongly attractive due to the proximity of both basins.The one-center carbon energy components are less negative than the values obtained with the 6 -31G basis set as shown in Table III.An effective C-C interaction could be computed by summing up the twocenter C-C and C-X and the one-center X contributions.The computed (C-C) eff energy contribution is still more negative than the corresponding C-C value in Table III by more than 0.5 a.u.This difference can be assigned to the decrease of the C energy component mentioned above.

V. SUMMARY
The partitioning of the SCF energy in terms of one-and two-center interactions in the AIM framework has been carried out.The one-electron contributions have been calculated by using a slightly modified version of the PROAIMV program.The two-electron integration over the atomic basins has been computed with a small program written by us.These large scale numerical integrations have proved to be extremely costly from a computational point of view and large supercomputation resources have been necessary.The results obtained in the test calculations are in excellent agree-ment with the chemical notion of molecules consisting of interacting atoms.Also, a generally good qualitative agreement has been observed with the recently introduced approximate LCAO energy decomposition scheme ͑CECA͒.Nevertheless, important differences have also been found for some molecular systems such as diborane.The presence in some cases of non-nuclear attractors destroys the chemical picture of the molecule and is, therefore, undesirable.

TABLE II .
One and two-center energy ͑a.u.͒ components for H 2 , N 2 , BH, and HF molecules using the 6 -31G(d,p) basis set.

TABLE I .
Integration parameters and energies ͑a.u.͒ computed for the N 2 molecule using the 6 -31G(d,p) basis set.Total energy stands for the sum of one and two-center energy components and AIM energy corresponds to the sum of the AIM atomic energies ͓⌺ A E(A)͔.Exact one-electron energy: Ϫ194.963 3; exact two-electron energy: 86.019 3; exact total energy: Ϫ108.9439.

TABLE III .
One and two-center energy ͑a.u.͒ components for ethane, ethene and acetylene.

TABLE IV .
One and two-center energy ͑a.u.͒ components for diborane using the 6 -31G(d) basis set.Subscripts b and t in H atoms hold for bridge and terminal, respectively.

TABLE V .
Energy ͑a.u.͒ components for acetylene with non-nuclear attractor ͑X͒ in the middle of the C-C bond using the 6 -31G(d,p) basis set.The (C-C) eff value is computed as E(C-C)ϩ2E(C-X)ϩE(X).