Decomposition of Molecular Integrals into Atomic Contributions via Becke Partitioning Scheme: a Caveat

: Decomposition of molecular integrals into physically meaningful atomic contributions by means of the Becke integration scheme requires some care with respect to the choice of suitable atomic size adjustments. Using a simple illustrative example, it is shown that the adjustment of cell boundaries, as originally proposed by considering Bragg-Slater atomic radii, does not pro vide reliable results. Alternatively, the positions of the bond critical points of the electron density can be adopted to define heteronuclear cutoff profiles which allow for a more reasonable atomic partition of the molecular electron density.


INTRODUCTION
HE theoretical resolution of molecular properties into atomic contributions is a long standing research field which dates back to the early days of quantum chemistry. Among all various methods presented so far, the most popular are surely those based on the integration of density functions connected to the properties of interest. [1] However, other approaches based on the atomic resolution of quantum mechanical operators might be recalled, as, for example, the breaking down of the molecular electric polarizability, [2] magnetizability, [3] and electric hyperpolarizability, [4] which have been obtained using force and torque gauges.
In the context of density functions, the quantum theory of atoms in molecules (QTAIM) has played a fundamental role, [5] see for example the atomic partitioning of polarizabilities [6] and magnetizabilities [7] presented by the Bader group. More recently, Hirshfeldbased partitioning scheme for the intrinsic polarizability density [8] and origin-independent decomposition of the static polarizability [9] have been presented.
The breaking down of molecular integrals into atomic terms provided by the Becke numerical integration scheme [10] has come recently into play, see for example the numerical study on the partitioning of the molecular polarizability by Mei et al. [11] and the resolution into atomic contributions of the nuclear magnetic shielding constants proposed by the Sundholm group. [12] In this work we wish to examine some aspects of the Becke scheme, when applied to the determination of atomic contributions, which deserve a bit of attention not paid so far.
As it is very well known, Becke scheme [10] deals with the evaluation of three-dimensional integrals of the type and assumes that a relative weight function wA(r) can be assigned to each atom A in the molecule such that and such that each wA(r) is equal to unity in the vicinity of its own nucleus and vanishes in a continuous and wellbehaved manner near any other nucleus. The molecular function F(r) is decomposed into single-center components since, as it can be easily proven, Therefore (1) reduces to a sum of single-center integrations IA over each atom in the molecule: Evaluation of an atomic integral can be best performed adopting spherical polar coordinates with origin on the nucleus of atom A, even if other numerical techniques can be adopted.
A key step within the Becke scheme is the evaluation of the nuclear weights wA(r), which are obtained as normalized cell functions. [10] These can be determined assuming a homonuclear system or adopting some kind of atomic size adjustments. In the original formulation of the method, Bragg-Slater (BS) atomic radii [13] were used in the adjustment of cell boundaries and only recently a new alternative [14] has been presented which adopts the bond critical point (BCP) positions of the electron density. [5] However, it should be noted that the result of summation (5), i.e., the final value of the three-dimensional integral, is independent of the way the nuclear weights are determined. The reason to use atomic size adjustments is that they provide a faster convergence to the required accuracy.
Having said that, however, it can be observed that the integration of molecular density functions, carried out with or without atomic size adjustments, provides very different atomic partitions, as given in Eq. (6). For example, atomic contributions to the nuclear magnetic shielding constants in HF and H2O molecules, calculated via the Becke scheme for the homonuclear partition and atomic size adjustments (ASA) based on Bragg-Slater (BS) radii and bond critical points (BCP) of the electron density, are reported in Table 1 (for calculation details see later).
As can be observed, total values are invariant with respect to ASA, whilst atomic contributions undergo considerable variations. Remarkably, the contribution of hydrogen to its own shielding pass from being prevalent for the homonuclear partition and BS to be minor for BCP in both molecules.
This opens a question about the physical significance, if any, of the atomic breakdown obtained with the Becke partitioning scheme. To answer this question we require a term for comparison, i.e., some theoretical or experimental quantity to compare with, which is hard to devise for the magnetic shielding example given above, as well as for many other molecular properties. As we will see in the following, a much simpler example will permit to draw some conclusions. However, before introducing the subject, let us first argue that, as it concerns the democratic subdivision of the molecular space given by the homonuclear partition, there is no reason to believe that it can provide meaningful atomic quantities, with the exception of highly symmetric cases. Instead, when atomic size adjustments are taken into account, some improvement in this regard is expected to occur. Using Bragg-Slater radii as guide in the adjustment of cell boundaries, neither the specific action that one atom has on the electron density nor the different substitutions that the same atom can have in different molecules seem properly taken into account. In this respect, the actual BCP positions of the electron density in each molecule under investigation can provide an optimal choice to guide the cell boundary adjustments. This can be easily implemented substituting the ratio of the Bragg-Slater radii χ = R i/Rj (see Eq. (A4) and omitting the restriction imposed by Eq. (A3) of Ref. [10]) with the ratio ri/rj of the distances from nuclei i and j to the BCP located between them. [14] METHOD In order to illustrate all these aspects, we have devised a very simple example in which the molecular function F(r) is taken as the ground state electron density function ρ(r) which integrates to the electron number N in the molecule (see for example Ref. [15], chapter 5) Defining a net charge on each atom as where ZA is the atomic number, it is straightforward to consider each qA as the charge transferred among the atoms in the molecule, since for a neutral molecule.
The α component of the total electric dipole moment for the point charge distribution corresponding to the qA' s can be calculated as where A R is the vector position of atom A (as can be easily verified, the dipole moment for the charge distribution does not depend on the origin of the coordinate system).
Comparing the calculated dipole moment for the charge distribution obtained via Eq. (12) with the expectation value of the dipole moment operator for the ground state electron density, will allow us to ascertain the physical significance of the different atomic partitions obtained by means of the Becke scheme. This has been done for a number of simple heteronuclear, polar molecules such as: HF, LiH, LiF, H2O, NH3. Owing to their modest size, very accurate computations have been carried out using the BHandHLYP functional, [16] adopting basis sets containing function types of high angular momentum taken from BSE. [17] In particular, on H and Li atoms we have adopted the aug-cc-pV5Z, which is a (6s5p4d3f2g) basis set on H and a (7s6p5d4f3g2h) basis set on Li. For N, O, and F atoms we have adopted the aug-cc-pV7Z, which corresponds to a (9s8p7d6f5g4h3i2j) basis set. Molecular geometries were optimized at BHandHLYP/aug-cc-pVTZ level using the Gaussian 16 program package. [18] Numerical integrations using the Becke scheme have been included within the SYSMOIC system of programs. [19] Nuclear magnetic shielding constants in Table 1 have been determined using the CSGT method.

RESULTS AND DISCUSSION
As it can be observed, the net atomic charges resulting from the partition of the electron density function via the Becke scheme, reported in columns 3 and 4 of Table 2, are in qualitative agreement with the electronegativity of the elements only for BCP based atomic size adjustments. For the Bragg-Slater based atomic size adjustments there is the notable exception of LiF, since the large χ = RLi/RF ratio provides a wrong direction of the charge transfer. The homonuclear partition provides counter-intuitive results for HF, H2O, and NH3 as well. At any rate, notable discrepancies can be observed comparing BCP and BS net charges, a fact which obviously depends on the different ratios of BS radii on the one hand, and distances between two joined nuclei, on the other hand. In this regard, we remark that in our calculations we followed exactly the Becke recipe, [10] which uses a modified radius for hydrogen atom of 0.35 Å. Adopting the original Bragg-Slater [13] radius of 0.25 Å, a choice made also in Ref. [11], the BS results move a little toward the BCP ones.  The comparison of the dipole moments obtained for the point charge distributions and the expectation values of the dipole moment operator permits an even clearer analysis. For the homonuclear partition the comparison is always rather poor and wrong in sign in three cases over five. The BS based electron partition provides dipole moments correct in sign (with the notable exception of LiF, as already remarked) but in large disagreement with the expectation value. On the contrary, the BCP based partition provides dipole moments which are in fairly nice agreement with μ   , with the only exception of water for which the comparison -3.2 vs. -1.9 Debye is not as good (even though net charges are in very good agreement with the QTAIM results, see Table 3.1 of Ref. [5]). It should also be considered that a dipole calculation in terms of point charges, can be a reasonable approximation only if the dipolar polarization of the atomic charge can be neglected, as exemplified by Bader for diatomic molecules. [5]

CONCLUSION
In conclusion, there is no doubt that the resolution of molecular properties into atomic contributions via the Becke scheme is sensitive to the atomic size adjustment, and this sensitivity can change the results even qualitatively. In the case of the electric dipole charges derived from cell boundaries either adjusted using Bragg-Slater radii or not at all adjusted, do not provide reliable results the former and lack of any solid physical ground the latter. The alternative way provided by the bond critical point positions to define heteronuclear cutoff profiles offers a promising solution for estimating atomic contributions to molecular properties. However, some care must be always paid by checking if the electron partition is reasonable on a physical ground, at least.