Electron-density topology in molecular systems : Paired and unpaired densities

This work studies the partitioning of the electron density into two contributions which are interpreted as the paired and the effectively unpaired electron densities. The topological features of each density field as well as of the total density are described localizing the corresponding critical points in simple selected molecules (local formalism). The results show that unpaired electron-density concentrations occur out of the topological bonding regions whereas the paired electron densities present accumulations inside those regions. A comparison of these results with those arising from population analysis techniques (nonlocal or integrated formalisms) is reported.


I. INTRODUCTION
The theory of atoms in molecules ͑AIM͒, introduced long time ago, 1 describes the topology of electron density in molecular systems providing important information which is essentially twofold.On one side, the electron density and its associated Laplacian field are characterized by the localization and classification of their critical points ͑cp's͒, that is, maxima, minima, and saddle ones.This aspect will be denominated the local formulation of the theory.On the other side, the topological definition of an atom or atomic region as the link of a density attractor and a zero flux surface basin enables one to perform studies of population analysis.As is well known, these studies allow one to evaluate classical quantities of paramount interest in quantum chemistry such as bonding populations ͑bond multiplicities͒, atomic populations, valences, free valences, etc.This second aspect will be called the nonlocal or integrated formulation of the theory.The localization of the cp's ͑of type bond, ring, and cage͒ and the value of the Laplacian of the charge density at these points, ٌ 2 ͑r͒, turn out to be essential in describing the electronic structure of the molecules. 1,2The cp's of the charge distribution ͑points of the distribution where ٌ =0͒ are classified according to the characteristics of the Hessian matrix of the function . 1,2The cp's that exhibit three non-null eigenvalues, two of them negative and one positive, are said to be bond critical points ͑bcp's͒ ͑3, −1͒ and they are of main interest in this work.The first two eigenvalues correspond to the perpendicular curvature and the latter one provides a curvature along the internuclear axis.Thus, a covalent bond is characterized by the electron cloud possessing two large negative curvatures perpendicular to the bond and a relatively small positive curvature along the bond at the position of the ͑3, −1͒ critical point. 1,2Another derived quantity is the Laplacian function, that is, the sum of the curvatures in the electron density along any orthogonal coordinate axis at the point r.Its sign indicates whether the charge density is locally depleted ͑positive͒ or locally concentrated ͑negative͒.This is a basic information to classify the behavior of the density around a local point. 1,2Moreover, other critical points under consideration in this work are the nuclear attractors ͑3, −3͒, that is, cp's which exhibit three non-null negative eigenvalues describing the electron-density concentration around the nuclear position.These points are used to define the topological bonding region as the line joining two nuclei involved in a bond.
The pth-order reduced density matrix ͑p-RDM͒ corresponding to an N-electron system is derived as an average of the N-electron-density matrix ͑nonreduced density matrix͒. 3,4he procedure to obtain a lower-order RDM from a higher one is called contraction mapping.Thus, the electron density, expressed as the diagonal part of the first-order reduced density matrix ͑1-RDM͒, can be derived from the contraction mapping from the second-order reduced density matrix ͑2-RDM͒ ͑pair density͒ to the first-order density matrix ͑1-RDM͒ ͑particle density͒.6][7][8][9] This result provides the appropriate method for studying the topology of each component in an independent way and consequently to draw out all the features of the electron density.The main goal of this work is to perform the partitioning of the electron density within the pth-RDM framework into two terms so that the properties of both of them can topologically be analyzed.For this purpose, a matrix formulation of RDM is briefly introduced and then the coordinate representation is used to describe the electron density and their components.Within this approach the relationships between correlation or many-body effects of the 1-RDM, related to cumulant matrices of 2-RDM, 10,11 are directly taken into account.This treatment allows us to incorporate correlation effects in a natural way, as well as to analyze the influence of many-body effects over this quantity and over its topology.
The organization of this article is as follows.Section II reports the theoretical framework of the partitioning of the electron density and the associated topological quantities.Section III is devoted to the description of the computational details of the calculations performed in selected systems as well as the discussion of the results.Finally, some remarks and conclusions are presented in Sec.IV.

A. The decomposition of the electron density
The spin-free N-electron density matrix of an N-electron system is defined by the kernel integral of the state function as where the integral symbol acts over the spin variables .The reduced density matrices of pth-order or p-electron matrices ͑p-RDM's͒ are obtained by contraction mapping from the N-electron density matrix, p D͑r 1 …r p ͉r 1 Ј…rЈ p ͒ = ͵ dr ͑p+1͒ …dr N N D͑r 1 …r p r ͑p+1͒ …r N ͉r 1 Ј…r p Јr ͑p+1͒ …r N ͒ ͑ 2͒ or in matrix notation, which is also called partial trace ͓tr ͑p+1͒…N ͑ N D͒ = p D, where the "tr" symbol means the mathematical trace operation͔.
The main interest in this work is to search for topological analysis of the electron density ͑r͒, which is the diagonal element of the spin-free 1-RDM, 1 D͑r ͉ rЈ͒, ͑r͒ = 1 D͑r͉r͒.͑4͒ This equation can be written down from the contracting mapping of the 2-RDM to 1-RDM, The introduction of the 2-RDM structure 10-12 where ⌳ jl ik are the matrix elements of the cumulant matrix of the 2-RDM ͓note that tr͑ 2 D͒ = ͑ N 2 ͒ and tr͑ 1 D͒ = N͔.
In this equation each contribution has a definite physical meaning, The 1 D ͑exch͒ matrix comes from the second term in 2-RDM decomposition ͑6͒, the exchange term, describing the paired electron populations, while the 1 D j ͑cumul͒ i term comes from the cumulant density and reads as where the u matrix is the effectively unpaired density matrix describing the unpaired electron cloud. 8Thus, the 1-RDM has two contributions, the paired and unpaired density matrices, respectively.In the coordinate representation Eq. ͑8͒ is Therefore Eq. ͑4͒ which defines the electron density reads as ͑r͒ = ͑p͒ ͑r͒ + ͑u͒ ͑r͒, ͑12͒ where ͑p͒ ͑r͒ and ͑u͒ ͑r͒ mean the paired and unpaired density contributions to the total density ͑r͒.These contributions are defined by

B. The topological study
The above-described electron-density decomposition allows one to study the topology of the two different fields in order to characterize bonding concepts.Interesting information is obtained from the localization of ͑3, −3͒ and ͑3, −1͒ cp's for each contribution of the electron density.cp's for the total density are found throughout the gradient of the field by ٌ͑r͉͒ r c = 0, ͑15͒ or according to Eq. ͑12͒, where r c = ͕r i c ; i =1,…M͖ indicates the critical point set of the total density ͑r͒.Then the straightforward relation holds, which means the opposite direction in which each field component increases/decreases the value of the density, i.e., when ͑p͒ ͑r͒ increases its value, the remaining part, the unpaired contribution, properly does it in the opposite direction.Then no simultaneous increase/decrease is possible in the cp's of the density components.As has been mentioned in the Introduction, the aim of this work is to search critical points of each field, to determine their shifts from those of the total electron density, and to interpret their meaning.6][7][8] In the case of a double occupation numbers such as the closed-shell Hartree Fock or the density-functional theory ͑DFT͒, in which ͑u͒ ͑r͒ is intrinsically zero, are fulfilled, and therefore no compensation is possible as in Eq. ͑17͒.

III. COMPUTATIONAL DETAILS, RESULTS, AND DISCUSSION
The AIMPAC ͑Ref.13͒ modules were appropriately modified in order to deal with paired and unpaired densities.The search for cp's, their characterization from the ͑p͒ and ͑u͒ densities, and the relieves and contour plots of the ͑p͒ and ͑u͒ functions were applied to selected simple molecular systems.The molecular wave functions were obtained from the GAUSSIAN98 package 14 using the basis sets 6-31G ** .For all systems, the geometries were optimized for these basis sets within the configuration-interaction ͑CI͒ wave functions with single and double excitations ͑SDCI͒.The numerical examples reported in this section attempt to use the present decomposition of the total electron density ͑r͒ into paired ͑p͒ ͑r͒ and unpaired ͑u͒ ͑r͒ components in order to show that the bonding phenomena have a pairing nature from the point of view of the local formulation of the topological description.For such a goal the results are reported in two tables, some selected maps and relieves of the densities and a schematic graphic showing the nuclear and critical point positions for each density.
Table I reports the topological properties of the total density as nuclear and bonding cp's, the distance between the nuclear attractor cp and the actual nuclei and the values of Distance between the nucleus and the corresponding nuclear attractor.

144116-3
each component of the cpys.Table II shows the distances between cp's of the total density and of those of each component.This table must jointly be read with Fig. 1 ͓scheme of the relative shifts of the nuclear ͑3, −3͒ cps͔.Table I describes the results as follows.The range of the paired component of the density varies from one ͑H 2 molecule͒ to four orders of magnitude ͑NaCl molecule͒ greater than the unpaired counterpart.Thus, the paired density is close to the total electron density.The greater nonpairing effect, i.e., presence of appreciable unpaired density, may be observed over H atoms in all compounds.This effect decreases as a function of the environment from H 2 , H 2 O, HF, CH 4 , and C 2 H 4 ; it can be interpreted as follows: the presence of a heavy atom linked to H induces an "enlargement" of basins of the former ones compelling the H electron cloud to be closer to the valence electrons of heavy atoms, leading them to be paired.The values of and ͑p͒ at the nuclear cp's increase as a function of the atomic number, indicating the concentration of ͑p͒ over them.No such trend may be found for ͑u͒ ; their values at cp of are much smaller than those of ͑p͒ .Another important feature is the shift of the location of the total electron-density cp from the nucleus positions: all of them are practically zero except that corresponding to H atoms which exhibits a small but appreciable displacement.This is a measure of the polarization of the electron cloud around the actual nucleus which is more pronounced in H atoms, i.e., the nonpairing effects are enlarged.The direction of these displacements are shown schematically in Fig. 1.
The results of the topological analysis also reveal the positions of the cp's of the density components.Table II shows the relative position of each critical point of ͑p͒ and ͑u͒ from those of , jointly with Fig. 1, in which the direction of the ͑3, −3͒ displacements is shown.cp's are in all systems placed between those of ͑p͒ and ͑u͒ showing that ͑p͒ cp's are shifted towards the topological bonding region ͑space between nuclear attractors͒ while those of ͑u͒ do inversely, as indicated by Eqs.͑16͒ and ͑17͒.Heavy atoms exhibit a practically null displacement for the paired density cp while H atoms do not.However, all them are very small and of the same order as those of bcp.The same trends are followed by ͑u͒ cp despite all the displacements being greater by some orders of magnitude than that of ͑p͒ .The greater displacements are noted for ͑u͒ at ͑3, −1͒ cp.From the values in this table and Fig. 1 it may be noted that ͑p͒ cp's always are displaced to the topological bonding region while ͑u͒ ͑3, −3͒ follows the inverse displacement.Thus, the unpaired cloud moves towards the atomic or nonbonding region while the paired one to the bonding region.For all H atoms the nuclear position is "outer" than that of cp, i.e., the cp's are inside the internuclear region.For all of the other atoms except the C and Na ones in C 2 H 4 and NaCl molecules, the nuclear position is located between the ͑u͒ and cp's.This means that the cp's of ͑p͒ and are placed in the internuclear region, while that of ͑u͒ does not.All these results indicate the accumulation of paired density over the bonding region with the corresponding decrement of the unpaired one.
It is worthy to note that the ͑3, −1͒ cp displacements in diatomic heteronuclear systems do not vanish and are of importance for ͑u͒ .All their ͑3, −1͒ cp's are shifted towards the nuclear position of the H atom while the behavior of ͑p͒ is the opposite, i.e., towards the bonding region.This result means that unpaired density is accumulated very close to the H nucleus.In diatomic homonuclear systems such displacements vanish identically because of the symmetry of the systems ͑see the tables in which these numbers are not re-ported͒.The contour and relief maps for both components of the density are shown in Figs.2͑a͒-6͑d͒ for N 2 , F 2 , HF, NaCl, and C 2 H 4 systems.It may be noted for the systems that the ͑p͒ field is too much stronger than the ͑u͒ field ͑see

144116-5
Electron-density topology J. Chem.Phys.123, 144116 ͑2005͒ Table I and contour maps but not the relief maps because they have been drawn in different scales in order to be properly appreciated͒.The paired density for heavy atoms preserves the spherical form for their electron distribution while the unpaired one does not, showing the polarization effect also in the valence region as may be seen for F, N, and C atoms in F 2 , HF, N 2 , and C 2 H 4 systems.In a previous work 15 the electronic correlation effects have been considered by means of the topological analysis of the Laplacian field of the electronic charge density in C-n-butonium cations which yield the lowering of the charge concentration at the valence shell and its contraction near the nucleus.Thus, the present work allows to interpret the earlier results as a natural consequence of the nonpairing effects as shown from the importance of the unpaired density over the nuclear positions and the valence region.Following the contour maps over the lines of greater density closing two bonded atoms it may be observed that the ͑u͒ field exhibits an important increment of the curvature in the interatomic region following the order HFϽ C 2 H 4 Ͻ N 2 Ͻ F 2 ͑see the corresponding relief plots͒.This feature indicates that ͑u͒ when integrated in space ͑nonlocal or integrated population analysis formalism͒ has an important contribution to the cumulant populations the results of which are appreciable, as can be seen from Ref. 12, and follow the same trend.
In order to complete our discussion we choose two further well-known examples of molecular systems lying out of equilibrium: N 2 molecule at a different conformation around equilibrium geometry and ethylene molecule with a CH 2 group rotated a / 2 angle about its principal axis.These systems allow us to show how to interpret the information contained in both contributions of the density.The N 2 mol-ecule as in the equilibrium examples has been calculated in the SDCI approximation while in the case of the C 2 H 4 molecule and its rotated conformation a simple correlated generalized valence bond-perfect pair 16 ͑GVB-PP͒ ͑ − * ͒ has been used.It is well-known that the SDCI approximation does not properly work in this case 17 and therefore it would need to go beyond this size expansion and consider triples and fourth excitations with the corresponding computational cost.Besides the GVB-PP approach includes adequate correlation effects to describe this situation. 17The geometrical parameters have been chosen from the SDCI calculation for the planar conformation while the rotated structure has been taken from the literature. 17igure 7 show ͑p͒ and ͑u͒ maps for the nitrogen molecule as a function of R − R eq relative bonding length where R and R eq stand for the internuclear distance and the equilibrium bond length, respectively.The minus/plus sign below each of the subfigures indicates shortened or enlarged length from equilibrium conformation, respectively.Figures 8͑a͒,  8͑b͒, and 9͑a͒, 9͑b͒ show both the contributions of the density at the nuclear and bond critical points as a function of the relative distance, as indicated above.The graphics in Figs. 8 and 9 permit to note the depletion of both densities in the region in which the nucleus becomes far apart ͑dissociation region͒ while an appreciable concentration in the region where they are compelled to be closer than in equilibrium is noted.The value of ͑p͒ at the nuclear cp position ͓see Fig. 8͑a͔͒ is decreased from the compressed conformations towards an intermediated situation in which it shows a "frustrated" maxima as an attempt to reorder the electron cloud before reaching the equilibrium conformation.Then it reaches a constant value for long distances where the nucleus are far apart.This result may be interpreted as a predissocia- tion state.For ͑u͒ ͓see Fig. 8͑b͔͒ this path begins with an increment towards the "reordering" stage while finally tending to a small value for long distances.The densities at the bcp in Figs.9͑a͒ and 9͑b͒ show a uniform decrease of ͑p͒ vanishing at large distances, while the unpaired density exhibits maxima in the "compressed" region and a sharp decrement before reaching the equilibrium conformation to go on decreasing from it until it becomes null at large distances.This behavior is expected according to the physical meaning we have discussed when introducing these contributions.
The results for the C 2 H 4 system are shown in Figs. 10 and 11.We choose the perpendicular plane to the planar conformation containing the carbon atoms because the results can be observed in both conformations, i.e., planar and rotated.Figures 10͑a͒ and 10͑b͒ show the paired density relief and the unpaired one, respectively, on the carbon atoms in the planar conformation.The unpaired part is viewed as "loges" at both sides of the C nucleus.In Figs.10͑c͒ and 10͑d͒ the associated maps are presented.The bond can be clearly appreciated from the ͑u͒ map, and thus it may be noted that the main contribution to the unpaired density comes from this bond.The rotated structure, Figs.11͑a͒-11͑d͒, show that the CC bond is weaker than that of the planar one.This result can be interpreted by comparison of Figs. 10, 11͑c͒, and 11͑d͒, which show the contour maps for each component of the total density .The region between the carbon atoms ͓Figs.10 and 11͑c͔͒ shows a similar ͑p͒ contribution ͑density of contour lines in each graph͒ in the planar and rotated conformations while on the contrary the ͑u͒ component is lowered in the rotated conformation, as may be noted from Figs. 10 and 11͑d͒.This fact is due to bond which is present in the planar conformation ͓see Fig. 10͑d͔͒ and is missing in the rotated one, leading to a lower density in the bonding region ͓see Fig. 11͑d͔͒.Therefore at the region between carbon atoms the only remaining contribution comes from the paired component which stands for the character of the bond, leading the unpaired density to be concentrated near the C nucleus.This effect under the rotation of the CH 2 group shows the strong diradical character of this conformation.These results are in complete agreement with previous population analysis studies on ethylene molecule. 17

IV. FINAL REMARKS AND CONCLUSIONS
The numerical results found provide a clear conclusion: bonding phenomena are only supported by the paired component of the density and therefore once more the Lewis model is successfully revisited.The effectively unpaired electron density, arising from the corresponding density matrix, collects two effects, that of the spin density and that of the nonuniform occupation of molecular or natural orbitals, i.e., the departure of the ␣and ␤-electron clouds from the idempotency property.Only the second effect is present throughout the numerical examples in this work since all the systems studied are in their closed-shell singlet ground states.Excited states, radicals, and ionic systems also possess the other effect which will be considered elsewhere.The fact that the unpaired density is not involved in bonding phenomena does not mean that it vanishes in the internuclear or bonding region.As mentioned in Sec.III, some systems ͑see maps and relieves͒ possess in this zone a non-negligible unpaired density contribution.This result allows to interpret by means of this local formalism the onset of high cumulant populations in nonlocal or integrated population analysis formalism. 9However, it may be noted that these populations are very dependent on the basis set as well as on the size of the CI expansion used in the calculations, while the populations arising from the paired density turn out to be quite the opposite. 18Also we must mention that these quantities are useful to detect formation and bond breaking in the course of molecular regrouping in nonequilibrium conformations.Finally, we wish to emphasize that the topology of ͑u͒ is a complex problem.Only the features necessary to achieve the goals of this paper have been considered here.Other aspects will be studied in a forthcoming work.

FIG. 3 .
FIG. 3. F 2 molecule: ͑a͒ contour map of the paired density, ͑b͒ relief map of the paired density, ͑c͒ contour map of the unpaired density, and ͑d͒ relief map of the unpaired density.

FIG. 2 .
FIG. 2. N 2 molecule: ͑a͒ contour map of the paired density, ͑b͒ relief map of the paired density, ͑c͒ contour map of the unpaired density, and ͑d͒ relief map of the unpaired density.

FIG. 6 .
FIG. 6. C 2 H 4 molecule: ͑a͒ contour map of the paired density, ͑b͒ relief map of the paired density, ͑c͒ contour map of the unpaired density, and ͑d͒ relief map of the unpaired density.FIG. 5. NaCl molecule: ͑a͒ contour map of the paired density, ͑b͒ relief map of the paired density, ͑c͒ contour map of the unpaired density, and ͑d͒ relief map of the unpaired density.

FIG. 8 .
FIG.8.͑a͒ ͑p͒ at the nuclear cp for N 2 at different molecular conformations around equilibrium geometry; ͑b͒ ͑u͒ at the nuclear cp for N 2 at different molecular conformations around equilibrium geometry.

FIG. 11 .
FIG. 11.C 2 H 4 relieves and maps of ͑p͒ ͓͑a͒ and ͑b͔͒ and of ͑u͒ ͓͑c͒ and ͑d͔͒ for the 90°rotated CH 2 group conformation viewed from the perpendicular plane containing the CC axis.

TABLE I .
Topological local properties of the total density ͑͒: components ͑ ͑p͒ , ͑u͒ ͒, critical points at the CISD/ 6-31G ** level of calculation.All quantities are in atomic units.

TABLE II .
Topological local properties of the paired, ͑p͒ , and unpaired, ͑u͒ , densities at the CISD/ 6-31G ** level of calculation.All quantities are in atomic units.
a Distance between paired and total density cp's.b Distance between unpaired and total density cp's.c No ͑u͒ ͑3, −1͒ cp's are found in the topological bonding region.FIG. 1. Relative nuclear coordinates ͑capital letters A and B͒; total electron density ͑circle͒, paired electron density ͑concentric circles͒, and unpaired electron density ͑star͒ ͑3, −3͒ cp scheme locations.