Energy decompositions according to physical space partitioning schemes.

This work describes simple decompositions of the energy of molecular systems according to schemes that partition the three-dimensional space. The components of those decompositions depend on one and two atomic domains thus providing a meaningful chemical information about the nature of different bondings among the atoms which compose the system. Our algorithms can be applied at any level of theory (correlated or uncorrelated wave functions). The results reported here, obtained at the Hartree-Fock level in selected molecules, show a good agreement with the chemical picture of molecules and require a low computational cost in comparison with other previously reported decompositions.


I. INTRODUCTION
The techniques of population analysis have proven to be powerful devices for the insight of molecular structures and chemical bondings. The use of these treatments has allowed one to carry out evaluations of classical chemical quantities as atomic charges, bond indices, valences, free valence indices, etc., in a satisfactory way, which turns out to be extremely useful for chemists. As is well known, the studies of population analysis are based on the partitioning of a certain molecular property, usually described by means of reduced density matrices, into contributions associated with each atom or group of atoms in the molecule. Hence, these procedures constitute an appropriate connection between the quantum mechanics and the intuitive chemical concepts. The most popular techniques in population analysis are those of Mulliken type, 1 in which the partitioning is carried out in the Hilbert space spanned by the basis function set, and those of topological type, 2 in which the ordinary physical space is decomposed into atomic domains.
The study of bond orders provides an important chemical information and consequently the numerical determination of these quantities from different population analysis methods has aroused the interest of many authors. [3][4][5][6][7][8][9][10][11] Similar to the evaluation of bond orders, the study of the decomposition of the total molecular energy into one-and two-center contributions has also been tackled. This kind of partitioning presents yet a greater interest since an appropriate partitioning scheme may enable one to identify the intramolecular bondings and to perform calculations of the bonding strengths, provided that the two-center components are related to the interactions between two atoms. Recently, several works related with energy partitionings have been reported as within the Mulliken scheme 12,13 as within methods that partition the three-dimensional space. [14][15][16][17] Although the last procedures have been regarded as more realistic to describe chemical features, the approaches proposed so far have been extremely CPU demanding. 15,16 The main aim of this paper is to describe simpler partitioning schemes, which exactly decompose the molecular energy into one-and twodomain terms. As the former partitionings, 16 our treatments allow us to draw results in the chemical scale but require a lower computational expense. Although the framework of our formulas can be developed at any level of theory ͑correlated and uncorrelated wave functions͒, we have limited our calculations, performed in several molecules, to the Hartree-Fock ͓self-consistent field ͑SCF͔͒ case to establish an appropriate comparison with the results derived from previously reported methods.
The paper is organized as follows. Section II describes the derivation of an algorithm which exactly decomposes the energy of a molecule, according to a partitioning of the real space. Two different ways have been selected to develop this algorithm: in the first procedure the real space is decomposed into disjunct atomic domains according to Bader's atoms in molecules ͑AIM͒ theory; in the second procedure a fuzzy atom approach 16,[18][19][20] is used, in which the space is decomposed into overlapping atomic domains. Section III reports the computational details and the results found in selected molecules as well as the discussion of these results. Finally, we have dedicated Sec. IV to point out the remarks and conclusions of this work.

II. THE ENERGY PARTITIONING
We will refer to an N-electron molecule with clamped nuclei; the nonrelativistic electronic energy corresponding to a determined state of this system is in which i , j , k , l ,... are orbital functions of an orthonormal basis set-e.g., SCF molecular orbitals, natural orbitals, etc.,-A , B ,... are the nuclei of the molecule, R AB the distance between those nulcei and Z A , Z B ,... are the corresponding nuclear charges. T j i mean the matrix elements of the kinetic energy operator ͑ − 1 2 ٌ 2 ͒ . A V j i are the matrix elements corresponding to the ͑−Z A / r A ͒ operator and B jl ik = ͗ik ͉ jl͘ are the standard two-electron integrals in the ͑12͉12͒ convention. 1 D j i and 2 D jl ik denote the matrix elements of the spin-free first-order and second-order reduced density matrices, respectively.
In order to decompose the energy expressed in Eq. ͑1͒, let us rewrite this equation as The Kronecker deltas in Eq. ͑2͒ can be formulated in different manner according to the selected partitioning of the real space. These formulations lead to treatments which are developed in the next two sections.

A. AIM treatment
As is well known, Bader's AIM theory divides the whole physical space into disjunct atomic domains ⍀ A , which are defined by surfaces having zero flux in the gradient vector field of the electron density. 2 In this treatment each domain ⍀ A is generally associated with a determined nucleus A.
Following the procedure used in population analysis studies to determine bond orders, 10,21 let us express the Kronecker deltas as in which ͗i ͉ j͘ are the standard overlap integrals ͑where the integration is performed over the total space͒ and ͗i ͉ j͘ ⍀ A are the overlap integrals over Bader atomic domains ⍀ A ͑where the integration is limited to this kind of domains͒. The substitution of the Kronecker deltas in Eq. ͑2͒ according to formula ͑3͒ leads to a topological partitioning of the electronic energy E into terms of monoatomic and diatomic character, as a function of nuclei A , B ,... and their corresponding atomic domains ⍀ A , ⍀ B ,..., so that in which The partitioning of the electronic energy expressed in Eq. ͑4͒ is exact and valid at any level of theory. Its use only needs the standard one-electron and two-electron integrals, that is, the matrix elements T j i , A V j i , B jl ik , the overlap integrals ͗i ͉ j͘ ⍀ A , as well as the elements of the second-order reduced density matrix 2 D jl ik . According to previous interpretations, 17 the overlap integrals over the atomic region ⍀ A in Eq. ͑5͒ act as weight factors over the one-and two-electron terms thus providing the one-center energy corresponding to that region. Similarly, two-center contributions in Eq. ͑6͒ are obtained through the corresponding terms weighted by overlap integrals over one or two different atomic regions ⍀ A and ⍀ B . The total molecular energy also contains the internuclear repulsions which are, obviously, of diatomic nature and consequently the corresponding term is also included in Eq. ͑6͒.

B. Fuzzy atom treatment
The second procedure to partition the energy E in Eq. ͑2͒ is based on the so-called fuzzy atom approach. 16,[18][19][20] According to this method, a non-negative continuous weight function w A ͑r͒ is introduced for each atom A. These weight functions measure the degree in which a given point of space r is considered to belong to atom A, fulfilling the conditions Consequently, in this approach there are not any sharp boundaries between the atomic domains but a continuous transition from one to another. Using these tools, the Kronecker deltas can be expressed as which requires to carry out numerical integrations for evaluating the ͗i͉w A ͑r͉͒j͘ integrals. The introduction of these deltas in Eq. ͑2͒ leads to in which and where the overlap integrals act in similar way to the partitioning reported in preceding section.
As is well known, the AIM treatment could be regarded as a limit case of the fuzzy atom formalism in which the weight functions w A ͑r͒ are zero or one. However, we have developed both formulations in an independent way, in order to be able to compare numerical results arising from both versions. The closed shell SCF formulas of these treatments can straightforwardly be derived expressing the second-order reduced density matrix elements in Eqs. ͑4͒ and ͑10͒ by the well-known relationship According to the above equations, in our proposals the one-and two-electron integrals have to be evaluated in the usual way ͑over all space͒ and only the overlap integrals ͗i ͉ j͘ ⍀ A or ͗i͉w A ͉j͘ have to be calculated over the Bader domains or through the weight functions, respectively. Other partitionings related to the AIM decomposition of the physical space previously reported 15,17 use the one-electron integrals T j i and A V j i and two-electron ones B jl ik calculated over the atomic domains ⍀ A instead of the whole space. Although the computational effort of two-electron integrals over arbitrary regions of the space may be reduced using some reported methods, 22,23 it still requires to perform large-scale numerical integrations. Similarly, a recently reported treatment in terms of the fuzzy atom scheme 16 is based on the use of weight functions w A to calculate all type of integrals ͑the one-and two-electron ones͒ which obviously may be only carried out in a numerical way, increasing also the computational cost.
The formulation of Kronecker deltas and their decompositions in terms of the matrix elements of the atomic overlap integrals play a fundamental role within the framework of our proposals. As mentioned in Secs. II A and II B, these matrix elements act as proportionality factors in the assignment of energy to the different nuclei and nucleus pairs, thus providing natural partitionings of the total energy in terms of energy components of different physical nature, that is, oneor two-center character terms. Similar treatments have successfully been used both in population analysis studies 10,21 and AIM energy decompositions. 17 Another property which deserves to be commented is the invariance of our expressions under a unitary transformation. In fact, each term within our partitionings is a trace, in mathematical sense, of product of matrices that transform according to the general rule. Hence, each of these terms turn out to be invariant under this kind of transformations. This property has been checked computationally proving the theoretical predictions in our both topological and fuzzy atom procedures.
In the following section we report results derived from the two above described treatments as well as a discussion on their agreement with the chemical features of the studied molecular systems. In that section we also point out that each of these simple treatments provides similar or even better results than the corresponding previously reported energy partitioning developed within the same decomposition of the physical space ͑AIM type 15 or fuzzy atom type 16 ͒ requiring significant lower computational cost.

III. COMPUTATIONAL DETAILS, RESULTS, AND DISCUSSION
As mentioned in the Introduction, our energy decompositions can be applied at any level of theory although the numerical determinations performed in this paper have been limited to the SCF case in order to be able to compare our results with the previous reported ones. Hence, our calculations have been carried out in the SCF molecular orbital basis sets. The computational implementation of our partitionings requires the calculation of the one-electron integrals T j i and A V j i and the two-electron integrals B jl ik . These integrals, as well as the SCF molecular orbitals, have been computed using a modified version of GAMESS program, 24 while the overlap integrals calculated over the Bader regions ͗i ͉ j͘ ⍀ A have been obtained from a modified GAUSSIAN94 code. 25 The values of weight functions w A and the numerical integration of the expressions ͗i͉w A ͉j͘ have been determined with the code cited in Ref. 26 which follows a Becke integration scheme 27 based on the weight functions originally proposed by this author. These weight functions, which satisfy Eqs. ͑7͒ and ͑8͒, depend on the empirical Slater-Bragg atomic radii of the atoms composing the molecule under study. According to the Ref. 16 we have increased the radius of hydrogen to the value 0.35 Å but for fluorine atom we have used its covalent radius instead of the average of covalent and ionic radii suggested by those authors. The reported results have been obtained with the basis sets 6-31G and 6-31G͑d , p͒. For all systems, the geometries were optimized for the corresponding basis set within SCF wave functions. Table I gathers the results obtained within the SCF approximation, with the basis set 6-31G, for simple systems ͑H 2 ͒, hydrocarbons with different bond orders ͑CH 4 , C 2 H 6 , C 6 H 6 , C 2 H 4 , and C 2 H 2 ͒ and some second-row hydrides ͑NH 3 , H 2 O, and HF͒. The results reported in the third and fourth columns correspond to one-and two-center energy components derived from the partitionings described by Eqs. ͑4͒ ͑AIM scheme͒ and ͑10͒ ͑fuzzy atom scheme͒, respectively. All these two-center components refer to classical bondings and are negative; other two-center values are low and have not been included. Although it is not obvious to establish a simple correspondence between the dissociation energies and the results derived from the energy partitionings proposed in this work or in other reported treatments, [15][16][17] we have included in the fifth column the experimental values 28 in order to take into account a reference. Similarly, as we are dealing with procedures that partition the three-dimensional space the concept of "atom" is not well defined and consequently the comparisons between monoatomic energies and energies of isolated atoms are not direct. For the sake of comparison, the last column reports results calculated according to the alternative energy partitioning within the fuzzy atom approach reported in Ref. 16, in which the oneand two-electron integrals are evaluated by numerical integrations through the use of the weight functions w A . These results have been calculated with the code picked up in Ref. 29, where the kinetic energy has been dealt providing only one-center energy components. Another treatment assigning diatomic components from the kinetic energy pointed out in Ref. 16 has not been considered in this work.
The values described in columns three and four of Table I show that the energies corresponding to the A -H bondings ͑A =H, C, N, O, F͒ are closer to the experimental ones in our fuzzy atom treatment ͑E b values͒, except in the simpler systems ͑H 2 , CH 4 , and C 2 H 6 ͒ in which our AIM procedure ͑E a values͒ presents a better approximation. It means that probably the fuzzy atom method is more appropriate to describe the polarity of the A -H bondings. In relation with the energies of the C-C bondings, both approaches fulfill the tendency of the series C 2 H 6 , C 6 H 6 , C 2 H 4 , and C 2 H 2 , according to the bond orders ͑multiplicity͒ but the values obtained within the AIM treatment are closer to the reference ones in column five. A comparison between the values obtained from our fuzzy atom treatment ͑E b energies͒ and those obtained from the fuzzy atom approach in previous report, using the "simple" energy decomposition scheme, ͑E d energies͒ shows that our values are considerably closer to the experimental ones. Obviously, the computational cost of our fuzzy atom treatment is much lower, since fewer number of numerical integrations are required. Moreover, the error of the overall integration ͑the difference between the sum of all one-and two-center computed contributions and the exact SCF energy͒ is negligible in our treatment while is of the order of 1 -5 Kcal/ mol in the results E d . Table II collects results obtained with identical treatments to those mentioned in Table I but using the wider basis set 6-31G͑d , p͒. No significant changes have been found in the behavior of our two methods in relation with the A -H and C-C bondings from what has been commented above. The comparison between values of energies E b and E d ͑see Table I of Ref. 16͒ arising from both fuzzy atom procedures also deserves similar comments to those reported in Table I. However, the description of acetylene molecule in the AIM approach presents a non-nuclear attractor ͑identified as a dummy atom with zero nuclear charge within this approach and denoted as X in Table II͒, which is in agreement with previously reported results. 15 These authors have used a treatment within the AIM framework which requires the evaluation of integrals which are extremely costly from a computational point of view. However, the two-center energy values derived from Eq. ͑6͒ are closer to the experimental ones than those described in Ref. 15 and the results related with the non-nuclear attractor are similar in both treatments. A survey to Tables I and II shows that the influence of the basis set in the treatments proposed in this work do not change the qualitative description of chemical features in the studied compounds and the quantitative differences between the energies of the fragments derived from both basis sets are similar to those obtained in previous studies.

IV. CONCLUDING REMARKS
In conclusion, in this work we have proposed a scheme to perform decompositions of the electronic molecular energy at any level of theory ͑correlated or uncorrelated one͒ according to partitionings of the three-dimensional space. This scheme can be adapted to any decomposition of the physical space through an appropriate formulation of the overlap integrals. Two methods have been developed providing suitable mathematical algorithms in which the usual oneelectron and two-electron integrals are calculated over the whole space. The first one belongs to the AIM approach and carries out the calculation of overlap integrals according to that technique. In the second method the overlap integrals are calculated through a fuzzy atom procedure. The results arising from both methods provide an adequate description of the chemical features of the studied molecules and are comparable with those obtained in previous reported studies that require a much higher computational effort.