A new, more flexible definition of fuzzy Voronoi cells is proposed as a computationally efficient alternative to Bader's Quantum Theory of Atoms in Molecules (QTAIM) partitioning of the physical space for large-scale routine calculations. The new fuzzy scheme provides atomic charges, delocalization indices, and molecular energy components very close to those obtained using QTAIM. The method is flexible enough to either ignore the presence of spurious non-nuclear attractors or to readily incorporate them by introducing additional fuzzy Voronoi cells.

The identification of an atom within the molecule is a crucial concept that permits the decomposition of physical quantities into atomic (and often also diatomic) contributions. Bader's Quantum Theory of Atoms in Molecules (QTAIM)^{1} represents nowadays one of the most widely used definitions of atom in the molecule. In QTAIM, the topology of the electron density is cleverly used to derive the atomic basins and the sharp boundaries that define them. The resulting atoms in the molecule possess some particular properties, namely, the virial theorem is fulfilled within each basin of the system. Such topological atoms are extensively used in the literature for different applications, from a number of different electron population/distribution^{2} analysis to molecular energy components^{3,4} among others.

There are two disadvantages of the method. First, the analysis of the electron density may lead to so called non-nuclear attractors (NNA), i.e., basins that have no nucleus associated to it. Whereas such objects do exist in condensed phase^{5,6} or molecules such as Li_{2},^{7} a particular combination of level of theory and basis set may give rise to spurious NNAs in dipole-bound water cluster anions^{8} or simple molecules such as C_{2}H_{2}^{3} or Si_{2}H_{2}.^{9} The existence of such spurious NNAs makes any analysis in chemical terms rather difficult and ambiguous.^{9,10} Second, the detailed exploration of the topology of the electron density and specially the cumbersome integration over the complex-shaped atomic domains may be a rather time consuming procedure, as compared with alternative definitions of atom in the molecule. The systematic use of the method in quantum chemical applications of systems of growing complexity such as fullerenes may be compromised.

An alternative to QTAIM within 3D-space analysis are atomic definitions that introduce overlapping or “fuzzy” boundaries. Such “fuzzy” atoms are best represented by introducing a non-negative weight function

*A*and each point of the 3D space, which satisfies the requirement

everywhere. It is assumed that the atomic weight function is large “inside” of atom *A*and small “outside.” Of course, in the special case of QTAIM,

*A*and

Over the last years we have shown that many quantities such as bond orders, overlap populations,^{11} energy components,^{12} effective atomic orbitals, or local spins could be generalized and computed in the framework of the “fuzzy atoms.” For that mere purpose we have often made use of the simplest Becke's atoms; i.e., the fuzzy atomic Voronoi cells introduced by Becke^{13} in the context of his celebrated multicenter numerical integration scheme.

Becke's atoms were originally devised for effective numerical integration of three-dimensional functions of marked atomic character. Becke's atomic weights,

*B*≠

*A*; that is, the region (point) of the space occupied by the nucleus is fully associated to its corresponding atom.

The scheme can be formulated as follows. Let us consider two nuclei, *A* on the left-hand side and *B* on the right-hand side, at a distance *R*_{AB}. For any point of the space,

where *r*_{A} and *r*_{B} represent the distance of that point to nucleus *A* and *B*. The surface μ_{AB} = 0 corresponds to the perpendicular plane that bisects the internuclear axis (i.e., the face of the Voronoi cell), whereas μ_{AB} = −1 and μ_{AB} = 1 conform with the line extending from nucleus *A* to infinity and from nucleus *B* to infinity, respectively. Thus, the points of the space for which μ_{AB} < 0 are within the Voronoi cell of atom *A*. The following step function

can be used to define the Voronoi cell of atom *A*.

In order to generalize the scheme for a system formed by a set of atoms, one can define the atomic Voronoi cell of atom *A* as

The step function (3) can be replaced by a continuous, monotonically decreasing function in the range (−1, 1), fulfilling the requirements *s*(−1) = 1 and *s*(+1) = 0 in order to define the fuzzy Voronoi cells. Becke suggested the simple polynomial function

where

*f*

_{k}(μ) =

*f*[

*f*

_{k−1}(μ)]. The integer

*k*is known as the

*stiffness*parameter and controls the shape of the cutoff profile. The larger the value of

*k*the steeper is the cutoff profile.

The fuzzy Voronoi cells thus defined do not account for the different atomic sizes in heteronuclear systems, as the faces of the Voronoi cells exactly bisect the internuclear axis between the two neighboring atoms. A shifted cutoff profile can be obtained with the cutoff function (5) using as an argument the following transformed coordinate:

where −1/2 ⩽ *a*_{AB} ⩽ 1/2 in order to ensure that −1 ⩽ ν_{AB} ⩽ 1. The position of the cell boundary is now controlled by the value of the parameter *a*_{AB}, which in the new coordinates is given by the condition ν_{AB} = 0. At the point where the shifted cell boundary intersects the interatomic plane, one has

Becke suggested to relate the ratio *r*_{A}/*r*_{B} with the relative size of the atoms *A* and *B* making use of some set of reference *fixed* atomic radii

In his original paper, Becke used the set of empirical atomic radii of Bragg and Slater^{14} and the value *k* = 3, on the basis of a better performance of his numerical integration scheme.

Even though Becke's atoms are only mathematical constructs, we have shown that, with some special provisions, they can be readily used as atoms in molecules models, from which chemically sound values of atomic populations, bond orders, or one and two-center energy components can be obtained. However, the use of a set of *fixed* atomic radii is a clear limitation of this model, basically because same atoms are treated on equal footing in different chemical environments, i.e., partial ionic character of atoms may not be properly taken into account.

One way to sort out this problem is to use an internal criterion to determine the position of the cell boundaries between all neighboring atoms. Thus, instead of a set of fixed radii, one needs to establish a strategy to define the set of χ_{AB} values, where *A* and *B* represent any pair of atoms sharing a cell boundary. In Ref. 11, some of us proposed as a criterion to use the position of the extremum (typically a minimum) of the density along the internuclear axis connecting every two neighboring atoms to locate their cell boundary. That scheme was later referred to as Becke-ρ in Ref. 15. A similar scheme making explicit use of the bond critical points of the density was also proposed by Francisco *et al.*^{16} Not surprisingly, the Becke-ρ scheme typically yielded electron populations, bond orders, and energy components close to those obtained by applying the disjoint QTAIM atoms.^{15–17}

There is, however, a limitation of the original formulation of Becke-ρ, that also applies to the original Becke scheme itself. We just need to recall that the transformation (6) does not permit an arbitrary shift of the cell boundary because if

In order to be able to fully accommodate the criterion behind the original definition of the Becke-ρ scheme, we propose here the following alternative transformation:

One can easily verify that

_{AB}= −1 and μ

_{AB}= 1, and also if χ

_{AB}= 1, i.e., in the homonuclear case.

The alternative new transformation is monotonic for any value of χ_{AB}, so it can be applied to shift the cell boundary to any position along the internuclear axis. It is worth to note that the transformations (6) and (9) lead to different cutoff profiles, even within the limits of applicability of the former. The differences are minimal yet the new definition exhibits a somewhat less sharp profile, as shown in Figure S1 of the supplementary material.^{18}

We have implemented this new scheme, henceforth, topological fuzzy Voronoi cells (TFVC), in our APOST-3D program.^{19} For TFVC, we have adopted a stiffness parameter of *k* = 4. We have observed a better overall agreement with QTAIM atomic populations than for *k* = 3 without a significant loss of numerical accuracy in the overall integration. The agreement with the QTAIM bond orders or delocalization indices (DIs) also improves with *k* = 4, in particular for non-bonded atoms. The observation that in benzene the DI between carbon atoms in para position is larger than in meta position^{20} is better reproduced with this TFVC model.

In our implementation, a rather modest grid of 40 radial and 146 angular points per atom typically yields the total number of electrons with an error of ca. 10^{−3}. For the numerical integration of other quantities such as one-electron energy components, we use a larger grid of 70 × 434 points per atom. In a most recent implementation of QTAIM integrations,^{21} the author uses up to 0.87M grid points per atom for a typical run. The reduction of the atomic grid size in fuzzy approaches is dramatic due to the lack of strict atomic boundaries. These methods can integrate all centers at once and scale linearly with the total number grid points, essentially for any level of theory. Most recently, delocalization indices for a set of fullerenes were determined using Becke-ρ.^{22} We have included in the supplementary material^{18} a comparative analysis of the accuracy of TFVC and QTAIM for a challenging endohedral complex.

Another advantage of the TFVC scheme is that it is flexible enough to either ignore the presence of spurious NNAs, as described above, or to readily incorporate them if necessary by defining extra fuzzy Voronoi cells centered on the position of the attractors. In the supplementary material,^{18} we provide illustrative examples of the TFVC method with NNAs.

Finally, it is worth to note that, strictly speaking, there are no boundaries in the fuzzy Voronoi cells. Thus far we have referred to cell boundary between atoms *A* and *B* to the intersecting plane for which μ_{AB} = 0. Hence, each pair of Voronoi cells (atoms), whether they are chemically bonded or not, has a formal boundary, whose location must be established. In the Becke-ρ method,^{11} if two atoms were too far to be considered as *neighbors* (the distance was larger than twice the sum of their fixed atomic radii), the ratio given in Eq. (8) was used to locate the fuzzy cell boundary. This leads to some inconsistencies when dealing with intermolecular interactions of bond dissociation curves, where atoms sharing cell boundaries may be too far away to be neighbors. In TFVC, we have adopted another criterion that also avoids any use of arbitrary fixed atomic radii. Thus, a given pair of atoms are not considered as neighbors if their midpoint is closer to a third nucleus or attractor of the molecular system. In that case, the formal cell boundary is placed at the midpoint, that is, a χ_{AB} = 1 value is taken.

We have computed with TFVC the atomic populations and delocalization indices of the molecular set described in Ref. 2, which features a number of hydrides as well as prototypical organic molecules. The hydrides have been chosen because the limits of applicability of the original Becke-ρ scheme are easily reached when one of the involved atoms is a hydrogen. The wavefunctions and densities have been obtained at the fully optimized Hartree-Fock and B3LYP levels of theory with the cc-pVTZ basis set. The agreement between the TFVC and QTAIM atomic charges at B3LYP level for the molecular set is very good, as indicated by the values of R^{2} and the slopes and intercepts of Fig. 1. Homonuclear diatomics and equivalent atoms have been excluded in the comparison. The results at the Hartree-Fock level are very similar and are not discussed. The mean unsigned error between TFVC and QTAIM charges is below 0.12, and the larger discrepancies are observed for hydrides involving 3rd period atoms with electronegativity similar to that of the H atom. The agreement between the TFVC and QTAIM delocalization indices is also strikingly good (see Fig. 2). The mean unsigned error is 0.06 and the largest one is just 0.25. The corresponding plots using *k* = 3 are included in the supplementary material.^{18}

An obvious physical limitation of the TFVC method is that the cell boundaries are described by planes, whereas in QTAIM they are typically curved. However, neither the overlapping character of the boundaries nor their curvature seems to be critical for TFVC integrations to be close to QTAIM. In our experience, it is when the electron density on the interatomic region is very flat (and large) that the QTAIM and TFVC values are more likely to diverge (see the supplementary material^{18} for a deeper analysis).

Finally, for illustration purposes we have also compared the behavior of the TFVC scheme in the decomposition of the Hartree-Fock molecular energy.^{3,16} Table I gathers the different contributions (namely, kinetic, electron-nuclear repulsion, Coulomb, and exchange) to the one and two-center energy components for hydrogen fluoride at the RHF/cc-pVTZ level of theory. The QTAIM results have been obtained with the program AIMALL.^{23} The overall accuracy of the numerical integration, expressed as the difference between the total energy and the sum of the one- and two-center terms separately, is similar for all methods. The two-electron numerical integrations in TFVC and Becke-ρ schemes have been carried out using solely a pair of rotated grids of 40 × 146 points per atom, as described in Ref. 12. On the other hand, the value of the Laplacian of the density integrated over the atomic domain is a typical measure of the fulfillment of the zero-flux condition within QTAIM. The value obtained with TFVC (ca. 10^{−2}) is one order of magnitude smaller than that obtained with Becke-ρ scheme. This value is somewhat too large compared with QTAIM, however, it does not seem to influence too much the values of the different one and two-center contributions to the energy. The TFVC energy contributions are much closer to QTAIM that to the Becke-ρ ones, especially the purely electrostatic ones. It is worth to note that for the sake of better comparison, the Becke-ρ values on Table I have been obtained using *k* = 4. Thus, the main difference between TFVC and Becke-ρ in this case is the proper location of the fuzzy cell boundary.

In conclusion, the TFVC method represents fast and simple atoms in molecules scheme that can be routinely used to extract chemical information from large-scale *ab initio* calculations. The method presented here can be regarded as a general purpose computationally more efficient alternative to Bader's QTAIM.

Financial help has been furnished by the Spanish MICINN Project No. CTQ2011-23441/BQU, the FEDER UNGI08-4E-003, and the Generalitat de Catalunya SGR528 grants. Support from the Xarxa de Referència en Química Teòrica i Computacional is also acknowledged. E.R.-C. acknowledges support from the Spanish FPU program (Grant No. AP2008-01231).