通过价层电子密度分析展现分子电子结构

卢天，陈沁雪

北京科音自然科学研究中心，北京 100022

1 In troduction

Electron density lies at the heart of density functional theory(DFT). The first Hohenberg-Kohn theorem shows that electron density distribution of a system at its ground state essentially carries all information about it1. Despite of this, simple analysis on molecular electron density is incapable of uncovering any information of chemical interest. The reason is that, for almost all kinds of chemical systems, electron density always decays exponentially from nucleus of heavy atoms1. It is clear that from the monotonic pattern of the electron density,it is incapable of recovering well known chem ical concepts,such as covalent bond, hybrid of atomic orbitals, lone pair and multi-center electronic conjugation. However, if some special treatments are imposed on the electron density, fruitful and useful information may be gained. For example, the Fukui function2, which is defined as derivative of electron density w ith respect to the number of electrons, is shown to be a nice indicator of regioselectivity of organic reactions3,4. Moreover,if electron density is properly partitioned into atomic contributions, then atomic charge can be derived, which has great significance in theoretical and practical aspects5.

The nature and formation region of chemical bonds, as well as the locations where lone pair electrons occur, are problems which theoretical chem ists pay great attentions to. Sofar numerous methods have been purposed aim ing for shedding light on these points. The so-called Laplacian of electron density (∇2ρ) analysis is such a method, and it is one of the key ingredients of the famous atoms in molecules (AIM) theory6,7.∇2ρ is a three-dimension real space function w ith follow ing explicit expression8

where ∇2is Laplacian operator, ρ corresponds to electron density, r is position vector, and x, y, z are its Cartesian components. Positive and negative signs of ∇2ρ at a point indicate that electron density is locally depleted and concentrated,respectively. It is well known that if covalent interaction occurs between two atoms, then their valence electrons must be concentrated toward each other. This phenomenon is known as valence bonding charge concentration8,9. If in a valence region,electron density is locally concentrated as indicated by negative sign of ∇2ρ, but no chem ical bond is formed, then this region should be regarded as containing lone pair electrons.

The electron localization function (ELF)10is another real space function frequently used to study chem ical bonds, while it is not defined based on electron density like ∇2ρ but based on orbital wavefunctions. Its expression for closed-shell system is11

in which φ and η denote orbital wavefunction and orbital occupation number, respectively. ELF essentially is an indicator of measuring electron localization degree in local region, its value range is [0,1]. A large ELF value at a point suggests electrons are greatly localized around it, implying presence of covalent bond, lone pair or atomic inner shell. ELF has been substantially used for a w ide variety of systems, such as organic molecules, atom ic crystals, coordination compounds,clusters; and for different problems, such as studying aromaticity12, portraying favorable electronic delocalization path13,14and investigating variation of electronic structure during chem ical reaction15.

Aside from ∇2ρ and ELF, there are many other real space functions possessing sim ilar capacity of revealing molecular electronic structure, such as localized orbital locator (LOL)16,phase-space-defined Fisher information density (PS-FID)17,density overlap regions indicator (DORI)18and region of slow electrons (RoSE)19. Since they are currently not as popular as∇2ρ and ELF, they w ill not be mentioned in detail.

It is worth noting that deformation density (ρdef)20,21is a very useful concept and frequently employed in discussing interatomic interactions in crystallology studies22–25. ρdef is defined as difference between actual electron density and promolecular density (ρpro), the latter is constructed by simple superposition of spherically averaged density of all atoms constituting the system at their free-states. Clearly, ρdefis a direct reflection of electron re-distribution in the course of molecular formation and thus carries wealth of chemically interesting information.

The exponential decay behavior of total electron density irrespective of actual electronic structure makes it useless for discussing chem ical bonds. The essential reason leading to the monotonic decay behavior, as w ill be shown later, is that the spherically symmetric core electron density dominates the distribution character of total electron density in core regions.As a result, the informative details of valence electron density distribution are completely shielded. It is easy to imagine that if the density corresponding to core electrons is discarded, in other words, if one only focuses on distribution of valence electrons, then the electron density itself may also behave as a useful indicator of formation of chemical bond or lone pair in analogous w ith ∇2ρ, ELF and ρdef. Although this idea is very simple and prom ising, as far as we know no studies have systematically explored its viability. Actually, we found valence electron density rarely received deserved attentions in existing literatures.

In this work, we w ill examine whether or not valence electron density has adequate ability to exhibit actual characteristics of electronic structure of various kinds of systems. In the meantime, sim ilarities and differences between analysis result of valence electron density, ∇2ρ, ELF and ρdefw ill be compared, and some cautions of using the valence density analysis w ill be mentioned. We hope that this work w ill make readers better recognize the great value of valence electron density.

2 Com pu tational details and m ethodo logies

In this work, all quantum chem istry calculations were finished by Gaussian 16 A.03 program26. The popular B3LYP exchange-correlation functional27in conjunction w ith def2-TZVP basis set28were used throughout all tasks,including geometry optim izations, frequency analyses and generation of electronic wavefunctions. A ll geometries used in the studies were fully optim ized and confirmed to be m inimum since no imaginary frequencies were found. The Multiw fn program29version 3.4 developed by one of us was used to conduct all kinds of wavefunction analyses, which is open-source and freely available on Internet30.

There is no unique way to separate density of core electrons and valence electrons. In our work, in order to investigate valence electron density, we first carried out regular calculation to yield all-electrons wavefunction, and then set occupation number of the molecular orbitals composed of core atomic orbitals to zero; after that, the total electron density analyzed by Multiw fn just corresponds to valence electron density. Since energy level of core atom ic orbitals are much lower than the valence ones and thus unable to mix w ith them, assuming that there are Ncorecore electrons, valence electron density can be easily obtained by setting occupation number of the lowest Ncore/2 closed-shell orbitals or Ncoreopen-shell orbitals to zero.Unless otherw ise specified, only electrons at outermost shell are conserved.

It is worth mentioning that it is also possible to assign pseudopotential to all atoms except for hydrogen, so that only valence electrons are explicitly presented in the resulting wavefunction, which makes analysis of valence electron density more straightforward. However, this strategy was not employed by us, because quality of electron distribution at valence region may be deteriorated due to the approximate representation of core electrons.

The density of atoms in their free-states involved in the computation of deformation density were calculated at the same level as molecular calculations, namely B3LYP/def2-TZVP. The detailed procedure of realizing the spherical average of atomic density has been detailedly documented in the Multiw fn manual31and thus w ill not be described here. Since at the B3LYP/def2-TZVP level, the restrictive open-shell Kohn-Sham wavefunction of boron atom at its ground state is unable to get converged due to unsolvable numerical reason, while it is needed by Multiw fn during plotting deformation map for B+13 cluster, therefore, for this map another popular exchange-correlation functional M 06-2X32was adapted instead of B3LYP. According to our tests, the difference between deformation map at B3LYP and M 06-2X levels is generally negligible.

3 Resu lts and discussion

Fig.1 Valence, core and total electron density along the axis of HCN molecule.color online.

3.1 Hyd rogen cyanide

Before discussing more involved systems, we first use the simple molecule HCN to briefly show why valence electron density may be usable for exhibiting molecular electronic structure. In Fig.1, the black curve corresponds to total electron density profile, it decays exponentially from nuclear positions as mentioned earlier. It is clear that chemical bond and lone pair are undetectable from this curve, while they should exist in present system according to common chem ical know ledge. The blue curve in Fig.1 portrays core electron density, it can be seen that it dominates total electron density in the region near to the nuclei, thus it is naturally expected that if core electron density is shielded, then the remainder density which carries rich information about electronic structure w ill be exposed. Indeed,the valence density, as shown by red curve, clearly indicates that electrons accumulate between carbon and nitrogen atoms,rendering formation of the C―N covalent bond. In addition, at the right side of the nitrogen there is an evident peak of valence density, implying occurrence of lone pair electrons (however,decisive conclusion about this point should be resorted to exam ination of plane map, as shown below). Noticeably, red curve has respective peaks at nuclear positions of carbon and nitrogen, they must arise from the small amount of 2s valence electrons that penetrated into the core regions. The electrons lying at atom ic orbitals w ith higher angular quantum number obviously do not contribute to the peaks because these atom ic orbitals have nodal character at atom ic center.

It is important to notice that, if one attempts to discuss bonding problem by means of valence electron density analysis, the geometry must be fully optimized. For example,experimentally the equilibrium distance of Ne2dimer is 0.3094 nm33, and it is known that the two atoms are not covalently bonded. As expected, the valence electron map corresponding to the equilibrium separation (Fig.S1 in Supporting Information) does not show detectable electron density distribution between the two atoms. However, if the two neon atoms are artificially placed closely w ith interval of 0.10 nm,from the corresponding curve map of valence electron density(Fig.S2), it is seem ingly that electrons aggregate towards intermediate region between the two atoms and thus contribute to covalent bonding of Ne―Ne, which is apparently untrue.

Next, we draw plane map of valence electron density in the cutting plane of HCN molecule, see Fig.2. As a comparison, the maps of ELF, ∇2ρ and ρdefare also given together. From the color-filled map, the feature of electronic structure can be captured more easily than the curve map. The chem ical bonds of C―H and C―N are vividly seen in the valence electron density map in Fig.2, it is seen that the valence electron density at corresponding region is remarkably higher than surrounding regions. A lso, it is visible that nitrogen forms a lone pair, as there is a marked red zone at the terminal of HCN axis. The overall profile of ∇2ρ map is highly analogous to the valence density map. By comparing the two maps, it is found that negative regions of the ∇2ρ map, namely electron concentration zones, typically occur where valence density distribution is relatively high. The negative zones of ∇2ρ map clearly reveal the fact that covalent bonds and lone pair exist in the present system. Moreover, both the two maps show that C―N is a highly polar bond, as the greatest electron concentration and the highest valence electron distribution along the C―N bond are at the nitrogen side.

Most above-mentioned findings from valence density and∇2ρ analyses can also be extracted from the ELF and ρdefmaps in Fig.2. Their features are roughly identical, though many details are different due to their very different physical definitions. Noticeably, ELF and ρdef do not tend to exhibit the high polarity of the C―N bond like valence density and ∇2ρ,but more emphasize to unveil its multiplicity nature, as it can be seen that in the ELF and ρdefmaps, the high localization region and the main positive distribution of deformation density corresponding to the C―N bond are quite broad in the direction perpendicular to the molecular axis.

3.2 Hyd rocarbons

In this section, we employ ethane and cyclopropane as instances to exam ine features of ELF, ∇2ρ, ρdef, and especially,valence electron density. First we look at ethane. In order to represent distribution character of valence electron density and∇2ρ, their isosurface maps are plotted, see Fig.3. The two plots again are very similar w ith each other, implying that the very simple valence electron density in many cases is adequate to replace the role of ∇2ρ, which is more involved in the definition and more difficult to be calculated. In the figure, the C―H bonds can be vividly recognized, and it can be concluded that C―C bond exists in present system, since not only the valence density map shows that valence electrons of a carbon are strongly polarized towards another carbon ow ing to the bonding interaction, but also the ∇2ρ map clearly indicates that corresponding valence electron concentration has formed. In addition, the common chem ical know ledge that carbon is in sp3hybridization state in ethane can be well recovered from Fig.3,since the figure shows that valence electrons of carbon atom are mainly located at the four vertices evenly distributed around the nucleus.

Fig.2 Color-filled map of valence electron density and ELF, as well as contour line map of ∇2ρ and ρdef of HCN molecule.In the color-filled map, the region w ith value higher than the upper lim it of color scale is depicted as white. The X and Y axes correspond to coordinate (since it does not affect discussion, for brevity the scale is hidden, sim ilarly hereinafter). In the contour line map, solid red line and blue dash line represent positive and negative regions, respectively. color online.

Fig.3 Isosurface map of valence electron density and∇2ρ of ethane.

Then we consider cyclopropane, which contains a threemembered ring accompanying large ring strain effect. A reliable and versatile electronic structure analysis method should be able tofaithfully reveal the influence on C―C bonding due to the remarkable strain effect. In order to explore this point, plane maps of valence electron density, ELF, ρdefand∇2ρ are plotted for cyclopropane in the plane defined by the three carbon atoms, see Fig.4. Observations from all the four real space function maps are basically identical, that is the severe ring strain makes the C―C bonds bulge out w ith respect to ring center. From Fig.4a one can clearly see that the centroid of the valence electrons used toform the C―C bonds deviates appreciably from the straight line linking neighbouring carbons. It is interesting to check whether total electron density is also as useful as valence electron density in revealing bonding character of cyclopropane, therefore we plotted total electron density map, see Fig.S3. This figure is apparently not as intuitive as valence electron density map, rendering the importance of discarding core electron density.

3.3 Atom ic c lusters

Fig.4 Color-filled map of valence electron density (a) and ELF (b), as well as contour line map of ρdef (c) and ∇2ρ (d) of cyclop ropane in the p lane defined by the three carbon atom s.In the contour line map, solid red line and blue dash line correspond to positive and negative, respectively.For clarity, the nuclei of carbons are linked by purple straight line. color online.

Fig.5 Color-filled map of (a) valence electron density, (b) ELF, (c) ∇2ρ and (d) ρdef of cluster in its p lane.

In this section, we test if valence electron density analysis is viable for deciphering electronic structure of atomic clusters,whose bonding character is usually much more involved than organic species and difficult to be predicted according to classic chem ical concepts, mostly due to the various possible types of multi-center delocalization.

Color-filled maps of valence electron density, ELF, ∇2ρ and ρdefof B+13cluster in its plane are provided as Fig.5. Electronic structure characters unveiled from the four subgraphs are qualitatively the same, namely strong covalent bonds formed between each pair of neighboring boron atoms lying at boundary of the cluster. Moreover, we can draw conclusion that there is an evident three-center bond at the center of the cluster,the evidence is quite sufficient: Fig.5a displays that valence electron is heavily distributed in the three-center bonding region; In Fig.5b, ELF shows that electron localization among the three central boron atoms is prom inent; In Fig.5c, ∇2ρ indicates that electron density of the three boron atoms significantly concentrates towards the ring center, and finally,the ρdefin Fig.5d manifests notable increase in electron density in the three-center bonding region. Another electronic structure character can be found from Fig.5 is that boron atoms lying at inner shell of the cluster have conspicuous bonding interaction w ith some boron atoms located at outer shell. For example, the three-center bond of B1―B7―B8 and the four-center bond of B2―B3―B4―B5 can be clearly identified from valence electron density map. This finding can be further consolidated by referring to the ELF, ∇2ρ and ρdef maps.

It is noteworthy that the valence density graph shown in Fig.5 do not convey full information about electronic structure of thecluster, because the graph is plotted on the cluster plane, only σ type of bonding interactions can be revealed,while π interaction is not reflected in the graph since π orbitals have nodal character along the plane and thus have vanished contribution to the graph. However, if one ignores all σ molecular orbitals and plots valence density map slightly above or below the cluster plane, then the π interactions can also be studied. Such a map is provided as Fig.S4. From the graph one can see that the multi-center bonds among the boron atoms at inner shell and those at outer shell of the cluster can also be found in π space; however, the plot suggests that no π bond is formed between the three central boron atoms. It is worth mentioning that if π and σ molecular orbitals are ignored when calculating ELF, then ELF-σ and ELF-π are obtained,respectively, which have been extensively used to study σ and π aromaticities34. The ELF-π map ofcluster slightly deviated to its plane is presented in Fig.S5, its distribution character is basically consistent w ith Fig.S4, implying the reliability and usefulness of employing valence electron density to analyze σ and π types of interatom ic interactions separately.

Another instance iscluster. Fig.6 presents isosurface maps of its valence electron density, ELF, ∇2ρ and ρdef. The isovalues are properly chosen so that distribution characteristics of the real space functions can be maximally exhibited.A lthough the isosurface shapes of various real space functions are different, all of them reveal the same truth that there exists a pair of four-center bonds respectively located at the center of Li1―Li2―Li3―Li5 and the center of Li1―Li2―Li3―Li4.Relatively speaking, due to the complex shape, the isosurface map of ELF is not as intuitive as that of other real space functions. On the other hand, the isosurfaces of valence electron density offer a very concise and clear description on the multi-center interactions.

Fig.6 Isosurface maps of valence electron density, ELF, ∇2ρ and ρdef of cluster.Isovalues are indicated in the parentheses. In the ρdef map, the blue isosurfaces correspond to negative regions. color online.

Fig.7 Color-filled map of valence electron density of CH3Br (a) and that of NH3 (b). Optim ized geometry of CH3Br ··· NH3 dimer is shown in (c).

3.4 CH3Br ···NH3 dim er

In this section, we illustrate possibility of using valence electron density to discuss non-covalent interaction,CH3Br···NH3dimer is taken as an example. This is a typical dimer bounded by halogen bond35,36. Since it is demonstrated that M 06-2X performs much better than B3LYP for this kind of interaction37, we optimized the dimer geometry at M 06-2X/def2-TZVP level, the final structure is presented in Fig.7c. How to interpret the underlying driven force of formation of this dimer? In fact if we look into the valence electron density of the monomers, the answer w ill become clear. Valence electron density of CH3Br and NH3 are shown as Fig.7(a, b),respectively. Note that when plotting Fig.7a, only the outermost shell of electrons of the brom ine atom, namely 4s and 4p electrons, are kept. If the subvalence 3s, 3p or 3d electrons are not discarded, the unique characteristic of valence electron distribution around the bromine atom w ill be obscured. From Fig.7a we can easily identify the lone pairs of the bromine atom, their distributions are highly anisotropic. As can be seen in the figure, almost no lone pair electrons occur at term inal of the C―Br σ-bond, this character is known as σ-hole36. The vacancy of valence electron at this place leads to significant deshielding of the positive charge carried by brom ine nucleus,and hence, the σ-hole zone effectively behaves as a local Lew is acid. In Fig.7b, the appearance of the lone pair of nitrogen atom in NH3 is clearly visible from the valence density distribution,the corresponding region can act as a local Lew is base due to the rich negative charge carried by the lone pair electrons.Since it is known that Lew is acid always tends to approach Lew is base as much as possible to maxim ize electrostatic attractive interaction, the valence electron density analysis nicely elucidates the underlying reason of form ing the CH3Br···NH3 dimer, and explains why in the dimer configuration the terminals of CH3Br and NH3are directly facing each other.

Fig.8 Energy variation along intrinsic reaction path (IRC) of Diels-Alder reaction between 1,3-butadiene and ethene.Isosurfaces of valence electron density of reactant (a), transition state (b), a featured IRC point (c) and product (d) are plotted as insets. Isovalues of all graphs are set to 0.25. Indices of carbon atoms are labelled in (a).

3.5 Study variation of elec tronic struc tu re du ring chem ical reac tion

Tracing variation of electronic structure along intrinsic reaction path (IRC)38can provide deeper insight into a reaction. An illustrative example is given as Fig.8. In the figure,energy variation in the process of a typical Diels-A lder reaction between 1,3-butadiene and ethene is shown, along w ith valence electron density maps corresponding tofour representative IRC points. From the graphs, below notable features can be found:

(1) In reactant complex, namely the initial point of the IRC,the character of the three C―C bonds in the butadiene are unequal, the C1―C2 and C3―C4 bonds are evidently stronger than the C2―C3 bond. This can be easily understood from different thickness of the cylindrical isosurfaces of valence density in the bonding regions. The isosurface profile of C1―C2 or C3―C4 is close to the C―C bond in ethene, which is a prototype of double bond. This observation is in line w ith bond order analysis, which suggests the boundary C―C bonds should be approximately regarded as double bond, while the central C―C bond is much weaker8.

(2) At transition state (TS), all the four C―C bonds in the Diels-A lder system have similar strength, as can be seen that the characters of corresponding valence density isosurfaces are almost identical. Actually, this conclusion can also be justified quantitatively by means of bond order calculation; the Laplacian bond orders of C1―C2, C2―C3 and C5―C6 are 1.61, 1.50 and 1.61, respectively. It is worth mentioning that the Laplacian bond order we proposed previously and employed herein is defined based on integrating negative part of ∇2ρ in fuzzy overlap space between the two interacting atoms. It is demonstrated that this bond order definition performs better than any other kinds of bond order in measuring bonding strength of organic system8.

(3) Inset (c) in Fig.8 corresponds to midpoint between TS and product. At this moment, the uniform ity of the C―C bonds at TS has been severely broken. From the corresponding valence density isosurface maps, one can readily recognize that C5―C6 is the current weakest C―C bond, while the C2―C3 has already become the strongest one, though it has the lowest strength at initial geometry.

(4) The final point of the IRC in Fig.8 is the product,cyclohexene. Formation of the new C1―C5 and C4―C6 bonds has finished. From inset (d) of Fig.8 it can be seen that C1―C5, C4―C6 and C5―C6 show equalized character, all of them can be attributed to typical single bond, this is in line w ith the Lew is structure shown in upper left side of Fig.8. From the isosurface map it is found that C1―C2 or C3―C4 is slightly stronger than the just mentioned C―C single bonds, this is indeed true and can be further confirmed by bond order analysis; the Laplacian bond order of C1―C2 is higher than C1―C5 by 0.17.

Overall, this example shows that smooth variation of bonding characteristic in the course of organic reaction can be successfully revealed by tracing change of valence density map along IRC, further demonstrating the powerfulness of the valence density analysis.

3.6 Basin analysis of valence electron density

Basin analysis is a popular technique of analyzing quantitative distribution character of real space functions. It was firstly introduced by Bader and extensively used in his AIM theory6. Then pioneered by Silvi and Savin, the basin analysis method has also been prevalently utilized in combination w ith ELF39,40. In the basin analysis, the whole molecular space is partitioned into individual spaces separated by zero-flux surfaces of the real space function under study.Each individual space is known as basin and contains one and only one maximum, which is called as attractor. Attractors and basins are one-to-one correspondence. Quantities of chem ical interest may be gained by integrating specific real space functions in the basins. In principle, basin analysis can be applied to all kinds of real space functions. In this section, we tentatively employ the basin analysis technique to valence electron density and exam ine which information could be extracted.

The results of basin analysis on valence electron density of ethane and acetonitrile are shown in Fig.9(a, b), respectively.The green spheres denote position of attractors, while the values stand for integral of valence electron density in corresponding basins. If the readers feel confusing on understanding the meaning of the attractors, they are suggested to refer to Fig.S6, in which the attractors are shown along w ith properly selected isosurfaces. From the data in Fig.9(a),follow ing findings can be obtained: Each carbon utilizes one of its valence electrons toform covalent bond w ith one of adjacent hydrogens, and approximately uses one valence electron to build the C―C single bond. In addition, four attractors are distributed at tetrahedron vertices around each carbon atom, indicating that carbon is in very standard sp3hybridization state. Aside from the basins at valence region,there are also basins found at carbon core regions, each of them contains 0.06 electrons, which could be interpreted as the amount of 2s electrons that penetrated into the core.

Next, we turn our attention to Fig.9b and explain the basins from top to bottom of the figure. The uppermost basin corresponds to lone pair of nitrogen, thus one may conclude that the total amount of electron that can be attributed to lone pair is 2.73. There is only one basin between the carbon and nitrogen, which contains 4.82 electrons and may be regarded as actual number of electrons contributed to the C―N multiple bonds. Like the ethane case, there are also two basins between the two carbons, however, the number of electrons in the two basins are no longer equal. In the present system, because of the larger electronegativity of cyano group w ith respect to methyl group41, the basin close to the former has evidently greater number of electrons (1.63) than the one close to the latter (0.66), reflecting the high polarity of the C―C bond. The electron w ithdraw ing nature of the cyano group also indirectly affects the C―H bonds in the methyl group. By comparing the methyl group in Fig.9a and in 9b, it can be found that the C―H bonds in acetonitrile are slightly more polar than those in ethane, since in the case of acetonitrile the basin at the carbon side has detectably higher number of electrons (1.06) than that at the hydrogen side (0.92), while this phenomenon does not appear in ethane, where both the two kinds of basins have 1.00 electrons. Besides, it is found that character of core basins is quite stable, the number of electrons contained in this kind of basin kept unchanged even external chem ical environment has varied a lot.

Fig.9 Attractors of valence electron density of ethane (a) and acetonitrile (b). The integral of valence electron density in thecorresponding basins are labelled. (c) Attractors of ELF of acetonitrile along w ith integral of total electron density in the ELF basins.color online.

For comparison purpose, basin analysis of ELF of acetonitrile was also carried out, the result is presented in Fig.9c. Notice that the real space function integrated in the ELF basins is total electron density rather than valence electron density. There are some degrees of sim ilarities and discrepancies between Figs.9b, 9c. Like the situation of basin analysis on valence density, basin analysis on ELF also yields basin corresponding to the lone pair of nitrogen and that corresponding to the C―N bond. According to the ELF convention, they can be referred to as V(N) and V(N,C),respectively, where V denotes valence. The number of electrons found in the ELF V(N) and V(N,C) basins are somewhat different to the counterpart obtained from valence density basin analysis. It is not possible to judge which analysis result is more reasonable, because of the physically unobservable nature of the quantity in question. However, it is seem ingly that the result of valence electron analysis is slightly more chem ically meaningful, because according to classical know ledge about chemical bond, the lone pair of cyano group should have exactly two electrons, and there should be three pairs and thus totally six electrons involved in the C≡N triple bonds, the basin analysis result of valence density is closer to this conventional belief than ELF analysis. The ELF basin analysis describes each C―H and C―C bond as a single basin, which is controversial to the two-basin representation of the valence density analysis and does not provide any information about bond polarity. A very common phenomenon observed in ELF analysis is that the result of integrating total electron density in core basin is slightly higher than the actual number of core electrons. For example, our ELF analysis on acetonitrile shows that carbon core basin C(C) has 2.10 electrons, which is ostensibly unphysical because it is larger than the expected 2.0.However, if we recall the findings in our valence density analysis, the reason of the overestimation of core electrons become a bit clear, it can be partially explained as consequence of penetration effect of valence electrons into core region.

3.7 Som e no tices on valence elec tron density analysis

Although above we have illustrated many successful use of valence electron density in shedding light on electronic structure of various kinds of systems, some notices should be kept in m ind when using this analysis method.

In our point of view, the major drawback of graphical analysis of valence density is that it is sometimes difficult or impossible tofind a proper color scale for color-filled map or isovalue for isosurface map to simultaneously exhibit electronic structure characters at various regions, the fluoroethane in Fig.10 is a typical example. In Fig.10a, the color scale is set to 0.00–0.33, the electronic structure characters of hydrogen and carbon atoms are clearly visible, however, in this setting no information about fluorine atom can be gained visually. When the color scale is changed to 0.00–2.00, as depicted in Fig.10b,the lone pair of fluorine atom is discrim inable, but information in other regions are completely lost. As a result, tofully investigate fluoroethane, one has to examine two plots of valence electron density, it is more or less cumbersome. The underlying reason leading to this contradiction is that fluorine atom possesses evidently more electrons than other atoms in this system. A possible strategy to overcome this drawback is to transform the valence electron density in a proper mathematical manner, so that its magnitude in different areas can always fall into a small range and thus can be graphically studied using a single color scale or isovalue.

Fig.10 Color-filled map of valence electron density of fluoroethane in its F-C-C-H p lane.Different color scales are used for (a) and (b). color online.

Another point should be noticed is that valence electron density analysis is uncapable of revealing character of certain type of interaction; for example, the Re―Re bond in [Re2Cl8]2−anion. It is w idely accepted that this is a quadruple bond w ith configuration of (σ2π4δ2)42. The σ bond results from overlap of dz2-pzhybrid orbitals of the two rhenium atoms, the two π bonds stem from overlap of their dxzand dyzorbitals, while the δ bond is formed by overlap of their dxyorbitals. The optim ized geometry of [Re2Cl8]2−is shown in Fig.S7, the valence density map and ELF map are given as Figs.S8 and S9, respectively.Notice that for rhenium, the def2-TZVP we used is a pseudopotential basis set w ith 4s, 4p, 5d, 6s electrons explicitly expressed, when plotting the valence electron density map we discarded molecular orbitals that composed of the subvalence 4s and 4p atom ic orbitals. From Fig.S8 it can be seen that the multiplicity character of the quadruple bond is undetectable, in other words, from the valence density map it is difficult to discriminate the quadruple bond against common σ bond.However, if we resort to ELF analysis, from Fig.S9 one can easily identify the unique feature of the quadruple bond,namely there are four relatively high localization regions linking the two rhenium atoms. Hence, for this and sim ilar cases, ELF is able to discover more information about bonding than valence electron density. In order to assess whether or not∇2ρ and ρdefare also capable of exhibiting special feature for the quadruple bond, the corresponding plane maps were plotted as Figs.S10 and S11, respectively. We find ∇2ρ is completely useless in current situation, it is even failed to show the existence of the Re―Re bonding, since there is no negative region between the two atoms. A lthough the conspicuous positive region of ρdefbetween the two rhenium atoms definitively proves the strong Re―Re covalent interaction,unfortunately, no additional information about the multiplicity nature of the quadruple bond is explicitly revealed.

4 Conc lusions

Electron density is a real space function that can be easily obtained theoretically, or experimentally via high-resolution crystallographic diffraction. It has long been thought that electron density is of little use on studying molecular electronic structure due to its monotonic character. However, this is no longer true if the density corresponding to core electron is removed. Despite its simplicity, valence electron density shows significant application value of revealing electronic structure of chemical systems. However, this point has not been substantially explored and recognized by theoretical chemists.In this contribution, we present a detailed investigation on valence electron density using a few chem ical systems of various kinds as examples. It is proven that graphical analysis of valence electron density is able to successfully unveil information about covalent bonds, lone pairs, multi-center interactions and so on. Non-covalent interaction and variation of electronic structure in the course of chemical reaction can also be nicely analyzed by means of valence electron density. If advanced analysis technique, namely basin analysis, is adapted,even more content of chem ical interest could be extracted from distribution of the valence electron density.

The prevalently used real space functions ELF and ∇2ρ, as well as the less popular but very useful function ρdefare also included in our discussion for comparison purpose. It is found that information conveyed by these functions are largely identical, despite different functions have very different physical meanings. For most situations, valence electron density can replace the role of other functions, w ith additional advantage of reduced computational complexity. However, as noticed and illustrated in the text, one should also recognize that the powerful valence electron density analysis is not completely free of drawbacks. In the circumstances when valence electron density analysis does not work well, other real space functions such as ELF should be employed together to supply complementary information.

In summary, the usefulness and reliability of valence electron density analysis for revealing molecular electronic structure is strongly demonstrated, we hope this work w ill make valence electron density receive more attention from chemists as it deserves.

Suppo rting In fo rm ation:available free of charge via the internet at http://www.whxb.pku.edu.cn.

References

(1) Koch, W.; Holthausen, M. C. A Chemist's Guide to Density Functional Theory; Wiley-VCH Verlag GmbH: Weinheim, Germany,2001; pp. 24–28.

(2) Parr, R. G.; Yang, W. J. Am. Chem. Soc. 1984, 106, 4049.doi: 10.1021/ja00326a036

(3) Fu, R.; Lu, T.; Chen, F. W. Acta Phys. -Chim. Sin. 2014, 30, 628.[付蓉, 卢天, 陈飞武. 物理化学学报, 2014, 30, 628.]doi: 10.3866/PKU.WHXB201401211

(4) Cao, J. S.; Ren, Q.; Chen, F. W.; Lu, T. Sci. China Chem. 2015, 45,1281. [曹静思, 任庆, 陈飞武, 卢天. 中国科学: 化学, 2015, 45,1281.] doi: 10.1007/s11426-015-5494-7

(5) Lu, T.; Chen, F. W. Acta Phys. -Chim. Sin. 2012, 28, 1. [卢天, 陈飞武. 物理化学学报, 2012, 28, 1.]doi: 10.3866/PKU.WHXB2012281

(6) Bader, F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: New York, USA, 1994.

(7) Matta, C. F.; Boyd, R. J. The Quantum Theory of Atoms in Molecules-from Solid State to DNA and Drug Design; W ILEY-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2007.

(8) Lu, T.; Chen, F. J. Phys. Chem. A 2013, 117, 3100.doi: 10.1021/jp4010345

(9) Gillespie, R. J.; Popelier, P. L. A. Chemical Bonding and Molecular Geometry-From Lewis to Electron Densities; Oxford University Press: New York, NY, USA, 2001; pp. 163–180.

(10) Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990, 92, 5397.doi: 10.1063/1.458517

(11) Lu, T.; Chen, F. W. Acta Phys. -Chim. Sin. 2011, 27, 2786. [卢天,陈飞武. 物理化学学报, 2011, 27, 2786.]doi: 10.3866/PKU.WHXB20112786

(12) Poater, J.; Duran, M.; Solà, M.; Silvi, B. Chem. Rev. 2005, 105,3911. doi: 10.1021/cr030085x

(13) Manzetti, S.; Lu, T. RSC Adv. 2013, 3, 25881.doi: 10.1039/c3ra41572d

(14) Manzetti, S.; Lu, T.; Behzadi, H.; Estrafili, M. D.; Thi, H. L. T.;Vach, H. RSC Adv. 2015, 5, 78192. doi: 10.1039/C5RA17148B

(15) Emam ian, S.; Lu, T.; Moeinpour, F. RSC Adv. 2015, 5, 62248.doi: 10.1039/C5RA08614K

(16) Schm ider, H. L.; Becke, A. D. J. Mol. Struct. THEOCHEM 2000,527, 51. doi: 10.1016/s0166-1280(00)00477-2

(17) Astakhov, A. A.; Tsirelson, V. G. Chem. Phys. 2014, 435, 49.doi: 10.1016/j.chemphys.2014.03.006

(18) De Silva, P.; Corm inboeuf, C. J. Chem. Theory Comput. 2014, 10,3745. doi: 10.1021/ct500490b

(19) Jacobsen, H. Chem. Phys. Lett. 2013, 582, 144.doi: 10.1016/j.cplett.2013.07.027

(20) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129.doi: 10.1007/bf00549096

(21) Lu, T.; Chen, F. J. Theor. Comp. Chem. 2012, 11, 163.doi: 10.1142/S0219633612500113

(22) Dunitz, J. D.; Schweizer, W. B.; Seiler, P. Helv. Chim. Acta 1983, 66,123. doi: 10.1002/hlca.19830660113

(23) Eisenstein, M. Acta Crystallogr. Sect. B 1979, 35, 2614.doi: 10.1107/S0567740879010001

(24) Cai, Y.; Luo, S.; Wang, Z.; Xiong, J.; Gu, H. J. Materiomics 2017, 3,130. doi: 10.1016/j.jmat.2016.12.005

(25) Cameron, T. S.; Borecka, B.; Kw iatkowski, W. J. Am. Chem. Soc.1994, 116, 1211. doi: 10.1021/ja00083a006

(26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb,M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.;Nakatsuji, H.; et al. Gaussian 16; Gaussian Inc.: Wallingford, CT,2016.

(27) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J.J. Phys. Chem. 1994, 98, 11623. doi: 10.1021/j100096a001

(28) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297.doi: 10.1039/b508541a

(29) Lu, T.; Chen, F. J. Comput. Chem. 2012, 33, 580.doi: 10.1002/jcc.22885

(30) Website of Multiw fn program. http://sobereva.com/multiw fn(accessed Aug 15, 2017).

(31) Multiw fn Manual, version 3.4, Section 3.7.4, available at http://sobereva.com/multiw fn (accessed Aug 15, 2017).

(32) Zhao, Y.; Truhlar, D. Theor. Chem. Acc. 2008, 120, 215.doi: 10.1007/s00214-007-0310-x

(33) Wüest, A.; Merkt, F. J. Chem. Phys. 2003, 118, 8807.doi: 10.1063/1.1566944

(34) Santos, J. C.; Andres, J.; Aizman, A.; Fuentealba, P. J. Chem. Theory Comput. 2005, 1, 83. doi: 10.1021/ct0499276

(35) Politzer, P.; Lane, P.; Concha, M. C.; Ma, Y.; Murray, J. S. J. Mol.Model. 2007, 13, 305. doi: 10.1007/s00894-006-0154-7

(36) Clark, T. WIREs: Comp. Mol. Sci. 2013, 3, 13.doi: 10.1002/wcms.1113

(37) Kozuch, S.; Martin, J. M. L. J. Chem. Theory Comput. 2013, 9,1918. doi: 10.1021/ct301064t

(38) Jensen, F. Introduction to Computational Chemistry; John Wiley &Sons: West Sussex, UK, 2007; pp. 487–492.

(39) Silvi, B.; Savin, A. Nature 1994, 371, 683. doi: 10.1038/371683a0

(40) Fuentealba, P.; Chamorro, E.; Santos, J. C. Understanding and Using the Electron Localization Function. In Theoretical Aspects of Chemical Reactivity; Toro-Labbé, A., Ed.; Elsevier B.V.:Amsterdam, The Netherlands, 2007; p. 57.

(41) Bratsch, S. G. J. Chem. Educ. 1985, 62, 101.doi: 10.1021/ed062p101

(42) Cotton, F. A.; Murillo, C. A.; Walton, R. A. Multiple Bonds between Metal Atoms; Springer Science and Business Media, Inc.: New York,USA, 2005.