Covalent Bonding in the Hydrogen Molecule - ACS Publications

Article pubs.acs.org/JPCA

Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

Covalent Bonding in the Hydrogen Molecule George B. Bacskay*,† and Sture Nordholm*,‡ †

School of Chemistry, The University of Sydney, Sydney, NSW 2006, Australia Department of Chemistry and Molecular Biology, The University of Gothenburg, SE-412 96 Göteborg, Sweden

‡

S Supporting Information *

ABSTRACT: This work addresses the continuing disagreement between two schools of thought concerning the mechanism of covalent bonding. According to Hellmann, Ruedenberg, and Kutzelnigg, covalent bonding is a quantum mechanical phenomenon whereby lowering of the kinetic energy associated with electron sharing, i.e., delocalization, is the key stabilization mechanism. The opposing view of Slater, Feynman, and Bader has maintained that the source of stabilization is electrostatic potential energy lowering due to electron density redistribution to binding regions between nuclei. Following our study of H2+ we present an analogous detailed study of H2 where bonding involves an electron pair with repulsion and correlation playing a significant role in its properties. We use a range of different computational approaches to study and reveal the relevant contributions to bonding as seen in the electron density and corresponding kinetic and potential energy distributions. The energetics associated with the more complex electronic structure of H2, when examined in detail, clearly agrees with the analysis of Ruedenberg; i.e., covalent bonding is due to a decrease in the interatomic kinetic energy resulting from electronic delocalization. Our results support the view that covalent bonding is a quantum dynamical phenomenon requiring a properly quantized kinetic energy to be used in its description.

1. INTRODUCTION The concept of chemical bonds is undoubtedly central to our understanding of chemistry and as such it is an issue about which most chemists have strong and long held views and theories. The physical description of covalent bonding, i.e. the mechanism and energetic consequences of valence electron sharing among the atoms of a molecule, has been an especially contentious issue since the first quantum mechanical calculations on the H2+ and H2 molecules in 1927, by Burrau,1 Condon,2 and Heitler and London3 respectively. The Heitler− London treatment of H2 has become standard textbook material and the source of the archetypal valence-bond (VB) description of a molecule.4−6 It provides a quantum mechanical interpretation of Lewis’ theory,7 whereby a shared pair of electrons from different atoms corresponds to a single covalent bond. The above calculations,1−3 in addition to validating the applicability of wave mechanics to molecules (by yielding accurate predictions of the bond lengths and dissociation energies of H2+ and H2) provided valuable insight and understanding of the quantum mechanical description of shared electrons. © XXXX American Chemical Society

Ninety years on, the lack of consensus concerning the physical interpretation and description of the above quantum mechanical results is quite remarkable, especially in light of the many thousands of sophisticated quantum chemical calculations of ever increasing accuracy, which have since been carried out. The widely adopted view that is prevalent in most general chemistry textbooks is that covalent bond formation in a molecule such as H2 is accompanied by an increased electron density in the interatomic region where it experiences a stronger attractive field due to two nuclei, rather than one in the original atomic environment. As noted by Slater8 in 1933, such an essentially electrostatic mechanism is consistent with the predictions of the virial theorem, whereby bonding is associated with a drop in the potential energy of the molecule, despite a concomitant rise in the kinetic energy. This view was reiterated by Feynman9 and also by Coulson whose famous book Valence10 from 1952 on has had strong influence on the Received: September 8, 2017 Revised: November 11, 2017

A

DOI: 10.1021/acs.jpca.7b08963 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

energies that are below the atomic value and repulsive potential energies in accord with Hellmann’s14 view of a covalent bond. However, orbital optimization (via ζ) results in qualitative changes in the kinetic and potential energies, although the effect on the total energy is quite small. While these energy changes due to orbital optimization result in the virial theorem being satisfied, do they also bring about a qualitative change in the covalent bonding mechanism? The electrostatic description of covalent bonding was challenged by Ruedenberg18 in 1962. In his landmark paper, The Physical Nature of the Chemical Bond, Ruedenberg concluded that covalent bonding is a quantum effect, since the critical component of bonding is electron sharing, i.e. delocalization, which is associated with constructive interference and results in a net decrease in the kinetic energy of the molecule. However, covalent bonding is accompanied by orbital contraction that is essentially an intra-atomic effect,18,19 as it rescales the molecular wave function (in particular its AO building blocks), decreasing the potential energy while the kinetic energy is increased. In addition to reducing the total energy, the rescaling process ensures that the virial theorem is satisfied.20−22 In summary, Ruedenberg18 emphasizes that “any explanation of chemical binding based on an electrostatic, or any other nonkinetic concept, misses the very reason why quantum mechanics can explain chemical binding, whereas classical mechanics cannot.” Subsequent papers by Ruedenberg and co-workers,23−32 have further elaborated and extended the original ideas, culminating in the physical explanation for orbital contractions31,32 in that “they further enhance the delocalization and, hence, the associated kinetic energy lowering. This kinetic energy lowering through enhanced interatomic delocalization drives the intraatomic contractions, which per se are antibonding, even though they generate the negative potential binding energy component at the equilibrium geometry.”32 Our own work over a period of 25 years,33−41 while in broad agreement with Ruedenberg’s views, has departed significantly from earlier mechanistic interpretations of covalent bonding as it has pointed to a deeper quantum dynamical origin of covalent bonding. Noting that Thomas−Fermi (TF) theory,42,43 the original and simplest type of density functional theory (DFT),44,45 is unable to describe covalent bonding,46,47 we analyzed the reasons for this failure in an effort to better understand the basic physics of bonding.33,36,39,40 We have found that, at least in simple systems, the major source of the problem is the simplified semiclassical form of the TF kinetic energy functional, resulting in a theory that is unable to account for nonergodicity and slow interatomic electron transfer and thus for any hindered internal electron dynamics in atoms and molecules. As it is the relaxation of these dynamical constraints that facilitates electron transfer in molecules, i.e., electron sharing, it is clear that covalent bonding is a dynamical process, which is implicit in the phenomenon of delocalization. This underlying dynamical origin is, we believe, the reason why analysis in terms of potential or kinetic energy can appear paradoxical. It is also an underlying motivation for us to verify that, allowing both potential and kinetic components to be resolved, a detailed energy analysis confirms the leading role of kinetic energy in the bonding mechanism. Kutzelnigg48−50 and Goddard51 have also embraced and extended Ruedenberg’s original ideas while Bader11,52,53 has strongly criticized Ruedenberg’s description of chemical bonding, reiterating that chemical bonding is “a result of the

chemical community in matters of bonding. As the virial theorem is satisfied by all accurate quantum chemical calculations, the above view has gained considerable general support, most notably, and quite recently, by Bader,11 in a personal tribute to Slater, emphasizing that Slater’s original views on bonding are totally correct. Thus, bonding came to be understood to be an essentially electrostatic phenomenon, inasmuch as the movement of electron charge into the interatomic “binding” regions results in increased attraction to the nuclei. The underlying role of quantum mechanics is seen as the prediction of the appropriate wave functions and hence charge distributions. Recent decades saw the development of “chemical topology” and its application to a wide variety of chemical systems, mostly to the description of electron density, largely stemming from and inspired by Bader’s quantum theory of atoms in molecules (QTAIM).12 Yet, the status of chemical concepts within a topological approach is open to question.13 In particular, as noted by Ruedenberg,13 “topology per se cannot yield physical concepts,” e.g. the description of covalent bonding at a fundamental physical level. In direct contrast with the electrostatic model of bonding, in the same year (1933) that Slater first proposed it, a different interpretation was offered by Hellmann.14 He suggested that covalent bonding should be understood as a quantum mechanical effect, brought about by the lowering of the ground state kinetic energy associated with the delocalization of valence electron motion between the atoms of a covalently bound molecule. This explanation differs radically from the prevailing electrostatic view, and as it appeared to violate the virial theorem, it failed to gain acceptance. The source of argument and disagreement about the origin of covalent bonding in H2, as well as in other molecules, is neatly summarized in Figure 1, showing the internuclear

Figure 1. Total VB energies (in black) and their kinetic (in red) and potential (blue) components computed with ζ = 1 (dashed lines) and with ζ optimized (solid lines) relative to the energies of two H atoms.

distance dependence of the binding energy and its kinetic and potential components, computed using VB wave functions (constructed from 1s atomic orbitals) with and without orbital optimization. The former calculations (with), originally formulated by Wang15 and Rosen,16 satisfy the virial theorem, while the latter (without orbital optimization), i.e. the original Heitler-London calculations, do not. The energy curves in Figure 1 are well-known, of course, and have been reproduced in the literature on numerous occasions, most notably in Slater’s Quantum Theory of Matter.17 The Heitler−London calculations (with fixed orbital exponent ζ of 1) yield kinetic B


Article


one-electron system is very fortunate, it is all the more surprising that experts on complex many-electron systems cannot agree on what H2+ is telling us about covalent bonding. Apparently there is still a need for a clearer and more complete analysis of the bonding phenomenon in this simplest of molecules. Even with all facts laid bare, several interpretations may, of course, seem possible and defensible. The choice must then be made on the basis of the clarity and consistency with which they explain covalent bonding. In light of the disagreement that still exists between these two points of view we carried out detailed analyses of the quantum chemical predictions for H2+ based on the spatial characteristics of the electron density and associated energy changes that characterize the bonding process.41 The analysis has now been extended to H2. As our work concentrates on the spatial analysis of electron and energy densities, we are reexamining the bonding concepts of Ruedenberg using methods that are closer to those of Bader, presenting the former ideas in new light. Our work thus represents a bridge between the two approaches and we hope that it will help to resolve issues that are still contentious. In the case of H2+ we concluded41 that, in agreement with Ruedenberg et al.,18,19,23−32,61 electron sharing, i.e. delocalization, is the dominant mechanism of binding and therefore an essential component of any reasonable model of covalent bonding. These conclusions were underscored by the observation that simple models of bonding that include electron delocalization and little else, such as Hü ckel theory,64−67 have long been successfully used by chemists to interpret and understand a range of chemical phenomena. The mechanism of orbital contraction, an intra-atomic effect, yet a consequence of delocalization,31,32 is responsible for a significant increase (∼20%) in the stability of H2+ and a reduction in its bond length, since it results in a more compact molecule. Orbital contraction is effectively responsible for ensuring that the molecular wave function and energies satisfy the virial theorem, although that is not a critical aspect of covalent bonding. We also found that in order to understand the binding process and the overall effect of several contributing mechanisms it is necessary to investigate each in detail. Inspection of the total electron and energy density changes alone means that the information on individual effects of sharing, orbital contraction and polarization are ignored which can easily lead to the wrong conclusions. This difference in degree of decomposition and detail is a source of disagreement between our comprehensive analysis and that behind the electrostatic models, as the latter concentrate on the spatial characteristics of the total density and density differences, without attempting to account for the effects of the sharing, orbital contraction and polarization. Consequently, in our view, the electrostatic explanation of bonding (whereby bonding in a diatomic molecule is simply due to increased electron density in the internuclear “binding region” where it experiences an increased attractive field due to two nuclei) is an oversimplification of the quantum mechanical results where the variation in the kinetic energy is of fundamental importance. Unfortunately, to this day, most undergraduate chemistry texts opt for the electrostatic explanation when discussing the physics of covalent bonding. One may well question the need for a paper such as this, especially as the main proponents of the electrostatic view are no longer with us. Admittedly, Ruedenberg’s description of covalent bonding is widely accepted by the quantum chemical

lowering of the potential energy in the bonding region caused by the accumulation of density that attracts the nuclei”,53 in complete agreement with Slater’s8 original explanation of bonding. While in broad agreement with Ruedenberg’s ideas, Frenking and co-workers,54,55 building on the earlier work of Morokuma and Kitaura56,57 as well as of Ziegler and Rauk,58 have also proposed an energy decomposition analysis (EDA) that highlights the importance of electrostatic and Pauli repulsion effects in a range of covalently bound molecules. The results have an obvious bearing on the question of the origin of the bonding mechanism in a range of molecules but do not directly contribute to the central question pursued here, i.e., whether the bonding is predominantly an electrostatic or kinetic phenomenon. Quite recently, Levine and Head-Gordon,59 using an EDA method,60 have found that orbital contraction, the chief “bone of contention” in the debate between proponents of electrostatic or kinetic mechanisms of bonding, indeed has a major effect on bonds to hydrogen, but a very minor role in single bonds between heavy atoms. Noting that in a range of diatomic molecules, in particular B2, C2, N2, O2, and F2, Ruedenberg et al.31,32,61 found that the intra-atomic components of the kinetic and potential energies are positive and negative respectively, it is possible that these trends are not entirely due to orbital contraction but to polarization and rehybridization.59 In any case, the important result that has emerged from the calculations of Ruedenberg et al.31,32,61 is that the attractive component of the interatomic energy is kinetic, exactly as in H2+ and H2. Thus, in all these diatomics, the covalent bonds are associated with delocalization, i.e., electron sharing. The density based EDA scheme of Wu et al.,62 implemented within the DFT framework, is another potentially useful method of analysis, although thus far it has been applied only to intermolecular interactions. Clearly, to this day there are two very different and strongly held views of the physical origin of covalent bonding. Much of the difference, in our opinion, stems from the different methods of analyzing quantum mechanical results. Ruedenberg’s approach18,19,23−32,61 is essentially quantum chemical, using Hilbert space methods to obtain quantitative estimates of the intra- and interatomic energy components of covalent bonding. The dominant bonding component is interatomic sharing, i.e., electron delocalization, which is invariably accompanied by a decrease in kinetic energy. Bader et al.11,52,53 and other proponents of the electrostatic model focus on the electron density and its topological features and in particular the density differences between a molecule and its constituent atoms. Bonding is said to be “the result of the accumulation of electron density in the binding region as required by Feynman’s electrostatic theorem.”11 The concept of “binding” and “antibinding” regions, as proposed by Berlin,63 is central to Bader’s model, whereby the movement of density among certain spatial domains and the associated changes in the electrostatic energy of the molecule are the descriptors of a covalent bond. The manifestations of covalent bonding at issue here are ubiquitous. Nevertheless, it is generally agreed that, whatever the physical mechanism is, it should already be expressed in the one-electron molecule H2+. Thus, in principle, we do not need to resolve the complex character of correlated electron dynamics in order to identify the basic mechanism of covalent bonding. While this projection of the discussion onto a simple C


Article


co-workers, hence a direct comparison, as will be done here, is instructive and useful.

community, but less so by the more general chemical community or indeed by the authors of undergraduate texts. It is especially the future generation of chemists that our paper is primarily aimed at, using concepts and language that is widely accessible. Moreover, our dynamical view and description of covalent bonding is relatively new, although we aim to publish more about it. In our H2+ paper,41 we emphasized that the rigorous quantum mechanical representation of kinetic energy, largely absent in electrostatic theories, is essential to the description of covalent bonding. The failure of Thomas−Fermi theory,42,43,46,47 the original form of electronic DFT,44,45 to predict covalent bonding graphically underscores this. Indeed, the correction that produced the successful modern density functional theory was to replace the semiclassical kinetic energy with the quantum mechanical kinetic energy operator of wave function theories.68,44,45 Moreover, we concluded41 that no real understanding of the mechanism of covalent bonding is provided by the Hellmann−Feynman9,14 and virial theorems,8 or by the concept of “binding” and “antibinding” regions,63 all of which appear to have crucial roles in the electrostatic theory of bonding. An extension of the above studies of covalent bonding based on wave function quantum mechanics is to start with the nonbonding Thomas-Fermi form42,43 of DFT. The use of the von Weizsäcker69 (vW) gradient correction of the TF kinetic energy would then introduce bonding and, in its pure form, the vW kinetic energy would give us exact or nearly exact results for H2+ and H2but for these molecules alone. For larger molecules with Fermi (exchange) correlations among the electrons the kinetic energy becomes a mixture of TF and vW components and there is no exact kinetic energy functional available. This approach would cast additional light on the dynamical origin of the covalent bond and its representation in DFT but it would not change the basic premise of our analysis or, we believe, our conclusions about the key role of kinetic energy in an energy decomposition of the bonding mechanism. We plan to expand upon these basic facts in future publications. In this work we extend our studies on the H2+ system to H2. The former is “the amazingly simple prototype” of covalent bonding. Perhaps it is even so simple that there is doubt that it can resolve an argument of interpretation such as the one we have been engaged in. H2 is, of course, a real molecule in all respects. It has a pair of electrons which repel each other, giving rise to electron correlation effects, while its wave function must also account for the antisymmetry condition of electron statistics (which leads to the Pauli principle and aufbau structure of atoms and molecules). Thus, the original Lewis picture of covalent bonding,7 regarded as the “basis for the subsequent construction and generalization of VB theory”,6 applies to H2.4,5 While the “independent electron picture” of the electronic motion and corresponding electron structure must apply to the bond in H2+ it can only be an approximation of limited validity for H2. In fact, we shall instead use the VB wave function (which is exact in the limit of infinite atomic separation) as a starting point for our studies of H2. As in our previous work,41 we explore the spatial dependence of the kinetic and potential energies as well as of the density itself, in particular the spatial dependence of these quantities in regions that have been identified as associated with interference that is at the heart of delocalization and hence covalent bonding. The approach used differs from that developed by Ruedenberg and

2. THEORETICAL AND COMPUTATIONAL DETAILS The majority of H2 wave functions in this work are constructed using a minimal 1s(ζ) Slater atomic orbital (AO) basis set {ϕa, ϕb} with exponent ζ that is either optimized (yielding quasi atomic orbitals, QUAO) or set to their atomic values of 1. The spatial H2 wave function in terms of these QUAOs can be written in a general form as combinations of atomic (or covalent) and ionic terms as31 Ψ=

⎛ N ⎞1/2 ⎜ ⎟ {cos γ[ϕa(1)ϕb(2) + ϕb(1)ϕa(2)] ⎝2⎠ + sin γ[ϕa(1)ϕa(2) + ϕb(1)ϕb(2)]}

(1)

where N = (1 + Sab 2 + 2Sab sin 2γ )−1 , with Sab = ⟨ϕa|ϕb⟩

(2)

The electron coordinates, specifying the occupancy of the 1s QUAO on nucleus a or b, are simply written as 1 and 2. With γ = 0, the resulting wave function is known as valence-bond (VB) type, while γ = π/4, i.e., equal contributions of atomic and ionic terms, corresponds to the molecular orbital (MO) wave function, as it can be shown to be equivalent to the doubly occupied bonding MO (ϕa + ϕb). Allowing maximum flexibility, i.e. choosing γ variationally (at any internuclear distance) yields the configuration interaction (CI) type wave function. These three types of wave functions neatly represent three different types of electron dynamics in the H2 molecule: independent (MO), strongly correlated (VB) and optimally correlated (CI) electron dynamics (noting that the term correlation, as used here, is essentially static correlation, i.e. the longer range mechanism that ensures the correct asymptotic character of the wave function as the molecule dissociates.) Comparison of results obtained with these three different methods will therefore allow us to explore the role of electron correlation dynamics in covalent bonding. Allowing the occupied orbital basis greater flexibility, polarization in particular, via a complete active space SCF70−73 (CASSCF) calculation, will yield more accurate wave functions of the above types (1), with the QUAOs {ϕa, ϕb} replaced by {ψa, ψb}. The latter are expanded in terms of a chosen basis, such as Dunning’s aug-cc-pVQZ set.74,75 A robust and reliable approach for the generation of QUAOs, based on the Singular Value Decomposition (SVD) method, have been developed by Lu et al.76,77 and by West et al.78 It has been successfully utilized in the work of Schmidt et al.31,32 and also in the more recent work of West et al.61 In this work the QUAOs are the “modified AOs” (MAO) developed by Erhardt and Ahlrichs.79 Briefly, the desired MAO ψa is the eigenfunction with the highest eigenvalue da (typically ∼0.6) of the density operator D̂ A, which is the appropriate intra-atomic part of the total molecular density matrix D̂ . Thus, DÂ ψa = daψa

(3)

where DÂ =

⎧ Dfg χ and χ on center A f g ⎪

∑ |χf ⟩DfgA ⟨χg |, DfgA = ⎨

⎩0 ⎪

f ,g

otherwise (4)

D


Article

The Journal of Physical Chemistry A and {|χf⟩} is the basis set. In this work the molecular densities are obtained at the CASSCF/aug-cc-pVQZ level of theory. Dynamic electron correlation can be accounted for by extending the CASSCF wave function to include single and double excitations from it, i.e., via a multireference CI (MRCI) calculation.80−82 We note here that the ground state (two-electron) wave function of H2 is separable into a product Ψ(r)Θ(σ) of a symmetric spatial wave function Ψ(r) with the antisymmetric singlet spin function Θ(σ) (where r and σ represent space and spin coordinates). As the molecular Hamiltonian is nonrelativistic in these studies the (normalized) spin component integrates to one and, consequently, the energy changes that accompany the bonding process result entirely from variations in the (normalized) spatial wave function. The interference density, ρI, associated with bond formation at a given geometry R is defined as the electron density difference between the H2 molecule and its constituent (quasi) atoms at the same geometry,31 i.e.

Figure 2. Contour map of the interference density ρI in the xz plane of H2 based on the VB density at R = 1.4 a0 and ζ = 1.1695 (where z is the internuclear axis). The positive/negative contours in red/blue are ±0.002, ±0.004, ±0.008, ±0.016, and −0.032. The dashed lines represent zero ρI contours. The nuclei (a, b) are at (0 ± 0.7a0).

ρI = Ψ2 − (ϕa 2 + ϕb 2) = p[2ϕϕ − Sab(ϕa 2 + ϕb 2)] a b

ment, i.e., their values in the positive and negative ρI regions. An alternative measure of the charge accumulation in the bond is to compute ρI at the bond midpoint.31 We have done this along with the computation of QI for a range of internuclear distances. In our previous work41 on H2+, the spatial distributions of the electron density and the expectation values of the kinetic and potential energies and hence total energies were computed as functions of the internuclear coordinate z. In this work, specifically by way of comparison, the same type of calculations have been carried out for H2 as well. The basic idea is that a given property P is expressed in terms of contributions, as summarized by the following equations in the case of oneelectron properties:

(5)

where the bond order p is defined as p = N (Sab + sin 2γ )

(6)

(Note that other workers have used the term deformation density for ρI in eq 5.83,84) Similarly, the interference contribution to the kinetic energy, TI, (that is equivalent to the interatomic component, Tinter) is defined as TI = Tinter = ⟨Ψ|T̂ |Ψ⟩ − [⟨ϕa|T̂ |ϕa⟩ + ⟨ϕb|T̂ |ϕb⟩]

(7)

where T̂ represents the appropriate kinetic energy operators. In the case of the potential energy the interatomic component, Vinter, (in atomic units) is simply Vinter = ⟨Ψ|Vâ + Vb̂ + r12

−1

⟨P⟩ =

∫ Ψ*P̂Ψ dx dy dz = ∑ Pn

−1

− [⟨ϕa|Vâ |ϕa⟩ + ⟨ϕb|Vb̂ |ϕb⟩]

where (8)

where V̂ a and V̂ b are the Coulomb operators for the two nuclei, r12−1 is the electron repulsion operator and R is the internuclear distance. The square bracketed terms in eqs 7 and 8 are the intra-atomic contributions to the total kinetic and potential energies, Tintra and Vintra. A detailed analysis31,32 of Vinter shows that it consist of a quasi classical term, Vqc, defined as Vqc = ⟨ϕa|Vb̂ |ϕa⟩ + ⟨ϕb|Vâ |ϕb⟩ + [ϕϕ |ϕ ϕ ] + R−1 a a b b

(10)

n

+ R |Ψ⟩

Pn =

∫z

zn + 1 n

{∫

∞

∞

∫ Ψ*P̂Ψ dx dy} dz −∞ −∞

(11)

The Pn integrals within a cylindrical volume element (with x and y chosen sufficiently large to guarantee convergence) are computed numerically, using a 48-point Gaussian quadrature.85 In the case of two-electron integrals

(9)

Pn =

as well as a Coulombic sharing, Vsc, and two distinct types of interference terms, VI and VII. The two-electron integral [ϕaϕa| ϕbϕb] represents the Coulomb interaction of the electrons on centers a and b respectively. By way of illustration a contour map of a calculated interference density is shown in Figure 2 that clearly indicates the movement of charge relative to the quasi atoms as well as the sheet of zero ρI that separates the regions of density buildup and loss. The total charge QI that is actually moved into the interatomic (bonding) region is computed by numerical integration, using a 48-point Gaussian quadrature.85 (In the case of H2+, Schmidt et al.31 derived and used an analytic expression.) Using the numerical information on the location of the zero ρI sheet, we have also computed the kinetic and nuclear attraction energies associated with the charge move-

∫z

zn + 1 n

⎧ dz1⎨ ⎩

∞

∞

∫−∞ dx1 ∫−∞ dy1 ∫r dr2 2

ρ1(r1)ρ2 (r2) ⎫ ⎬ r12 ⎭ (12)

where ρ1 and ρ2 are the densities of electron 1 and 2 respectively, with coordinates r1 (x1, y1, z1) and r2 (x2, y2, z2). Thus, e.g. the electron repulsion integral [ϕaϕb|ϕaϕb] is obtained using ρ1(r1) = ϕa(r1)ϕb(r1) and ρ2(r2) = ϕa(r2)ϕb(r2). The numerical integration in eq 12 has been carried out using the VEGAS (Monte Carlo) algorithm.85 This approach of studying densities and energies, as described in this paragraph, while in practice simpler because of the rectilinear integration, is complementary to the method described above, with reference to Figure 2, where the x,y limits are z-dependent, i.e., where the effect of interference on the properties is studied by contributions to and from the regions of density increase and loss, respectively. E


Article


the wave functions do not account for dynamic correlation. Therefore, the question is whether accounting for these effects is necessary from the point of view of bonding in H2 and its basic mechanism. In other words, is it sufficient to rely on the simple minimal basis VB and CI theories to reach a definitive description of covalent bonding in H2? The diagrams in Figures S1, S2, and S3 of the Supporting Information provide a comparison of the energetics of the minimal basis calculations with those obtained by CASSCF and MRCI methods, in addition to the restricted (RHF) and unrestricted (UHF) Hartree−Fock calculations, using the augcc-pVQZ basis set. For the sake of comparison, results obtained by DFT44,45,68 using the UB3LYP functional92−94 are also presented. Clearly, dynamical correlation does not change the qualitative behavior of the energies (as judged by the close agreement of the CASSCF and MRCI curves), hence it is not expected to have a significant effect on the basic physical description of covalent bonding. This is in agreement with the findings of Schmidt et al.31,32 for H2 and other diatomic molecules, concluding that correlation does not change the physical interactions that lead to bonding, even though its effects on dissociation energies are significant. The close agreement between the minimal basis VB and CI energies and those from the large basis CASSCF and MRCI calculations similarly suggest, in agreement again with Schmidt et al.,31,32 that the essential elements of covalent bonding are captured by the minimal QUAO calculations, i.e., polarization functions are not essential to describe the fundamental aspects of bonding. It is well-known that RHF theory, including the simple MO method, is unable to describe the dissociation of H2 into atoms, but as the current plots suggest, for distances less than ∼2.5a0, i.e., in the region that is critical for the study of bonding, these approaches too are adequate. Although in this work DFT44,45,68,92−94 is not utilized beyond the plots in Figures S1−S3, the current results suggest that DFT too would describe bonding the same way as VB. 3.2. Scaling and Wave Function Contraction. Orbital contraction and expansion, that lead to contraction and expansion of the wave function of the molecule, in the case of minimal basis MO, VB and CI calculations are described by the optimization of the QUAO orbital exponent, ζ. A plot of the computed exponents at a range of internuclear separations are shown in Figure S4 of the Supporting Information. In the case of the VB wave function, orbital contraction occurs as R decreases from ∼2.8a0, while CI predicts its onset at the larger distance of ∼3.6 a0. In all cases, the maximum ζ is 1.6875 as R → 0, corresponding to the screened 1s2 configuration of the helium atom. As the VB and CI wave functions correctly describe the dissociation process, their QUAO exponents approach the limiting value of 1.0 as R → ∞ after dropping below 1.0, signifying orbital expansion. In the case of the MO calculation, the large orbital expansion that occurs is nonphysical, as it is associated with the incorrect dissociation. The orbital optimization process, carried out at each internuclear distance, ensures that the Virial Theorem is satisfied at all distances, i.e.

The choice of local kinetic energy functional, however, is nonunique.86−89 In this work (in atomic units), we use the 1 positive definite form 2 |∇Ψ|2 that yields positive values everywhere. Such a choice is consistent also with the recommendation of Ayers et al.,88 as “the most acceptable definition of the local kinetic energy” in the context of local temperature in DFT. In our previous work41 on H2+ we mostly 1 used the alternate Laplacian form − 2 ∇2 Ψ that in some regions does result in negative values of the local kinetic energy. The total kinetic energy, i.e., integral of the local kinetic energy, is of course independent of the choice of functional. While we recognize the loss of stringency associated with the nonunique nature of the local kinetic energy functional, making a sensible choice is all one can do if we want to understand and quantify the topological nature of the kinetic energy and hence total energy. The minimal basis calculations were carried out using our own programs, while the CASSCF, MRCI, and density functional calculations using the aug-cc-PVQZ (AVQZ) basis were done using the Gaussian0990 and MOLPRO2015 codes.91

3. RESULTS AND DISCUSSION 3.1. H2 Wave Functions. The simplest (spatial) wave function for H2, derived from that of H2+, is the MO wave function, ΨMO = ψg(1)ψg(2)

(13)

where the (normalized) bonding MO ψg is doubly occupied. Following Burrau’s calculations on H2+, the first MO calculation for H2 was reported by Condon2 in the same year, in 1927. While a MO wave function adequately describes the electronic structure of H2 in the region of its equilibrium geometry, it is well-known that its asymptotic behavior with increasing internuclear distance is incorrect, as it predicts dissociation to H atoms and ions. By contrast, the Heitler−London wave function (also published in 1927), being formulated in terms of 1s hydrogen AOs, is exact in the limit of infinite separation. Written in terms of MOs the Heitler−London VB wave function is ΨVB = [2(1 + Sab 2)]−1/2 [(1 + Sab)ψg (1)ψg (2) − (1 − Sab)ψu(1)ψu(2)]

(14)

where ψu is the antibonding MO, i.e. ψg , u = [2(1 ± Sab)]−1/2 (ϕa ± ϕb)

(15)

where {ϕa, ϕb} is the minimal AO (or QUAO) basis. The second configuration in eq 14, with both electrons occupying the antibonding MO, cancels the ionic terms that are implicit in the MO wave function. A more general and flexible wave function is of the CI type ΨCI = c1ψg(1)ψg(2) + c 2ψu(1)ψu(2)

(16)

where c1 and c2 are variational coefficients. These three wave functions account for covalent as well as ionic terms, fixed or variable contributions, as present in the MO and CI wave functions. This is especially clear when the wave functions are expanded in terms of the QUAOs, as in eq 1. The limitations of the above wave functions are 2-fold. First, the orbital basis contains no polarization functions and second,

1 1 T = − V − R(dE /dR ) 2 2

(17)

As shown by Hirschfelder and Kincaid,20 Coulson and Bell21 and by Löwdin,22 given a N-electron wave function Ψ(r1,r2,...; R), scaling (i.e., stretching or contracting) of the electron F


Article


densities and (intra-atomic) energies is a natural extension of the scaling process. Consequently, the description and explanation of the driving force for bonding is independent of the choice of orbital exponents or indeed of the choice of wave function, as concluded by Ruedenberg and co-workers on many occasions.18,19,23−32,61 Figure 3 presents a comparison of the interatomic contributions to the kinetic, potential and total energies of

position vectors relative to the appropriate nuclei, r1, r2,..., and of the internuclear distance R by the scale factor ζ yields Ψζ = ζ 3N /2 Ψ(ζ r1, ζ r2, ...; ρ)

(18)

where

ρ = ζR

(19)

(Note that in this section alone ρ represents a scaled internuclear distance, rather than density.) Consequently, given a wave function computed with ζ = 1 at a given internuclear distance ρ and its kinetic and potential energy expectation values T(1, ρ) and V(1, ρ), the scaled energies T(ζ, R) and V(ζ, R) are readily obtained. The scale factor ζ is given by the equation22 ζ=−

V (1, ρ) + ρVρ(1, ρ) 2T (1, ρ) + ρTρ(1, ρ)

(20)

where Tρ and Vρ are the derivatives of T and V with respect to ρ. The effect of scaling on the kinetic and potential energies at the internuclear distance R are then simply given as

T (ζ , R ) = ζ 2T (1, ρ)

(21)

Figure 3. Interatomic kinetic (Tinter), potential (Vinter) and total energies (Einter) of H2 computed with ζ optimized (full lines) and with ζ = 1 (dashed lines) computed at the VB level. The curve in bold black is the total binding energy (Einter + Eintra) with ζ optimized.

and V (ζ , R ) = ζV (1, ρ)

(22)

the total energy being just their sum, i.e. E ( ζ , R ) = T (ζ , R ) + V ( ζ , R )

H2. While the interatomic potential energies are quite insensitive to the choice of ζ, there is a large effect on the interatomic kinetic and hence total energies. As the total binding energy includes the intra-atomic contributions, the effect of orbital optimization on the former is not nearly as large as on its kinetic and potential components. As a numerical example of the scaling process for H2 consider the VB results, where the equilibrium internuclear distance obtained with ζ = 1 is 1.643 a0 and the kinetic, potential and total energies are 0.8377, −1.9537 and −1.1160 Eh respectively. Clearly, the virial theorem is not satisfied. Substituting these numbers (and the appropriate derivatives) into eqs 20−23 yields ζ = 1.1678, so that R(a0) = 1.643/1.1678 = 1.407, the resulting kinetic, potential and total energies being 1.1425, −2.2816, and −1.1391 Eh respectively. Although these energies satisfy the Virial Theorem as given in eq 17 (noting that at this point dE/dR = −0.0045), the calculation does not exactly predict the equilibrium geometry. In fact, the equilibrium bond length and energies are predicted by the calculation at ρ = 1.6489 that yields ζ = 1.1661, Re = 1.414 and kinetic, potential and total energies of 1.1391, −2.2782, and −1.1391 Eh respectively, exactly those obtained by direct orbital optimization at R = 1.414. While the differences between the two sets of results are practically negligible, we believe them to be instructive as they provide some insight into the virial theorem and how it works. A simpler estimate22,41 of the effects of orbital contraction can be obtained by neglecting the derivatives in eq 20. Starting with the results obtained with ζ = 1 at ρ = 1.643 we get ζ = 1.9537/(2 × 0.8377) = 1.1661. Hence R = 1.409 a0, while the kinetic, potential and total energies are 1.1391, −2.2782, and −1.1391 Eh, respectively, seemingly satisfying the virial theorem at the predicted equilibrium, i.e. 1 T = − V = −E (30) 2

(23)

Thus, by use of eqs 20−23 the wave functions, energies and orbital exponents (that minimize the total energy and satisfy the virial theorem at all geometries) can be obtained from the knowledge of the analogous unscaled (ζ = 1) quantities. In other words, orbital optimization in the case of the current VB and CI wave functions is formally equivalent to scaling the wave functions and energies computed with a fixed exponent. Further, if we partition the unscaled energies into intra- and interatomic contributions, i.e. write T (1, ρ) = Tintra(1, ρ) + Tinter(1, ρ)

(24)

and V (1, ρ) = Vintra(1, ρ) + Vinter(1, ρ)

(25)

the scaled interatomic energies (from eqs 19 and 20) are simply Tinter(ζ , R ) = ζ 2Tinter(1, ρ)

(26)

and Vinter(ζ , R ) = ζVinter(1, ρ)

(27)

As the intra-atomic energies scale the same way, it is clear that a change of orbital exponent by scaling will result in intraand interatomic components in exactly the original ratio, i.e. Tinter(ζ , R )/Tintra(ζ , R ) = Tinter(1, ρ)/Tintra(1, ρ)

(28)

and Vinter(ζ , R )/Vintra(ζ , R ) = Vinter(1, ρ)/Vintra(1, ρ)

(29)

According to this analysis, the definition of intra- and interatomic energies at a specific internuclear distance R and optimized exponent ζopt is exactly the same as at the scaled distance ρ(= ζoptR) and ζ = 1. Thus, in our view, an analysis of covalent bonding among atoms with contracted (or expanded) G


Article


Figure 4. Spatial distributions of the total integrated density (ρ) and kinetic (T), potential (V), electron repulsion (Vee) and total energies (E) in H2 at R = 1.4 a0 calculated with the VB method (1s, ζ = 1) (panel I) and corresponding difference maps (panel II), i.e. {H2 − [H(a) + H(b)]} where H(a) and H(b) are hydrogen atoms. Total energies (Eh) in panel I: T = 0.8452, V = −1.9506, Vee = 0.5277 E = −1.1055. Total energy changes (Eh) in panel II: ΔT = −0.1548, ΔV = 0.0494, ΔE = −0.1055.

The difference maps in panel II of Figure 4, which display the interatomic contributions to the density and energies, are much more informative from the bonding point of view. The transfer of density into the internuclear region is very obvious as well as substantial, leading to a uniform drop in kinetic energy. By contrast, while there is an analogous decrease in the potential energy in the interatomic region, the increase in the outer regions (where the interference component of the density is negative), more than cancels the decrease, resultig in a net increase in potential energy. The overall result, as discussed above, is that binding, i.e. a drop in total energy of the molecule relative to the atoms, is due to the lower kinetic energy of the molecule. Electron repulsion is obviously an important component of the total potential energy, resulting in substantial cancellation of the attractive electron−nucleus Coulomb energy. (A map of the two-electron integrals that contribute to the electron repulsion in a VB calculation at R(a0) = 1.4 are shown in Figure S8, along with the z-dependence of the total repulsion.) An alternative view to the explicit treatment of interelectron repulsion is of course to regard the electrons as screening the nuclei. The analogous maps to those in Figure 4, at R = 1.4 a0, obtained with the optimized exponent of 1.1695 are shown in Figure S9. At a qualitative level the results in Figure S9 are essentially the same as those in Figure 4. To display the effects of orbital optimization we computed spatial maps of the change in density and energy distributions in H2 resulting from a change in the orbital exponent from 1.0 to its optimized value ζopt. As we want to compare the effects of orbital contraction on the molecule with the analogous contraction on the (noninteracting) atoms, it is useful to provide a mathematical summary of the process. Given a property P, the effect of contraction in the molecule and free atoms,ΔPH2 and ΔP2H, are defined as

While this simple calculation accurately reproduces the scale factor and energies, its prediction of the equilibrium separation is slightly in error. Contour maps of the densities obtained with ζ = 1 and ζ = 1.1661 at their appropriate equilibrium geometries and the resulting difference map are shown in Figures S5 and S6, illustrating the contraction in density that is largely described by the scaling process discussed above. One dimensional analogues of the contour maps are shown in Figure S7. Clearly, there is a marked density increase between z ≅ − 1.5 and 1.5 a0, mostly around the nuclei, as expected, since the degree of contraction depends on the wave function itself and hence the density (see eq 18). By the same token, the spatial contributions to the energy (as considered in detail in the next section) will be scaled (according to eqs 21 and 22) and a change in ζ will not change the relative balance of intra- and interatomic energies. 3.3. Spatial Analysis of the Density and Energy Changes on Covalent Bonding in H2. The bond length dependence of the binding energy and its kinetic and potential components, as shown in Figure 1, summarizes the source of argument and disagreement about the origin of covalent bonding in H2 as well as other molecules, including H2+. In our previous work41 on H2+ we generated and studied density maps and density difference maps, including those of the kinetic and potential energy integrals so as to test the essential features of the electrostatic and the contrasting Ruedenberg theories. We used an approach that focuses on the spatial dependence of the density and energy integrals, as discussed in section 2, whereby a given property P is expressed in terms of contributions from each element, as summarized by eq 10. The results of such an analysis for H2 at the experimental internuclear separation of 1.4 a0, using a VB wave function constructed from H AOs (ζ = 1) are shown in Figure 4. In panel I, the contributions to the density as well as the kinetic, potential and total energies are shown, while in panel II the results of the same analysis are shown for the difference H2 − [H(a) + H(b)], i.e., the contributions, as defined by eqs 5−8. The sharp minima in the potential and hence total energies clearly identify the nuclear positions, where the density and kinetic energy have maxima, although in the case of the former such maxima are barely discernible due to the short bond length. The electron repulsion is very smooth by comparison, with a flat maximum.

ΔPH2 = PH2(ζopt) − PH2(ζ = 1)

(31)

and ΔP2H = P2H(ζopt) − P2H(ζ = 1)

(32)

The difference between ΔPH2 and ΔP2H can then be written as ΔΔP = ΔPH2 − ΔP2H = Pinter(ζopt) − Pinter(ζ = 1) H

(33)


Article


Figure 5. Panel I: Molecular contraction at R = 1.4 a0, i.e. difference maps of integrated densities (ρ) and kinetic (T), total potential (V), nucleuselectron interaction (Vne), electron repulsion (Vee) and total energies (E) between H2 with ζopt = 1.1695 and ζ = 1.0. Total energy changes (Eh): ΔT = 0.3007, ΔV = −0.3342, ΔVee = 0.0562, ΔVne = −0.3904, ΔE = −0.0336. Panel II: Difference between molecular and atomic contractions at R = 1.4 a0 (as defined by eq 33). Total energy differences (Eh): ΔΔT = −0.0670, ΔΔV = 0.0048, ΔΔE = −0.0623.

since the interatomic components of P with ζopt and ζ = 1 are defined as Pinter(ζopt) = PH 2(ζopt) − P2H(ζopt)

due to the interference component of the kinetic energy. These observations are in complete accord with the findings of Schmidt et al.,31,32 whereby orbital contraction enhances delocalization, i.e. the transfer of electron density to the interatomic region that results in a lowering of the kinetic energy and hence a lower total energy. The above results clearly demonstrate that orbital contraction affects the density and energy contribution from the immediate neighborhood of the nuclei and are contrary to the notion that the increased stabilization due to orbital contraction is brought about by the interaction between the increased charge density in the bonding region and the nuclei. The model, first suggested by Ruedenberg,18,19,23−32,61 that we should consider electron sharing between atoms with contracted densities, i.e., quasi atoms, as the source of covalent bonding, is strongly supported by the current studies. The results of the above type of calculations at internuclear distances of 1.2, 1.643, 2.0, and 3.0 a0 are summarized in Figures S10−S17, demonstrating the universal nature of orbital contraction and its effect on the density and energies of H2. Note, however, that at 3.0 a0 an orbital expansion occurs that changes the appearance of the difference maps in Figures S16 and S17 when compared with those at the shorter bond lengths where orbital contraction occurs. A question that arises concerns the role of electron correlation, especially as the results above were obtained at the VB level of theory. As noted already, the VB and CI methods account for static correlation and are therefore exact in the dissociation limit. In the equilibrium region, however, they resolve only a relatively small degree of dynamic (left−right) correlation. Thus, the uncorrelated MO method yields potential energy curves that agree well with those obtained at the VB and CI levels in the equilibrium region but not at stretched geometries (see Figures S1−S3). The inability of (restricted) Hartree−Fock methods (i.e., the more general versions of the MO method applied to H2 in this work) to correctly describe dissociation is of course widely appreciated. In order to demonstrate the similarities and differences in the methods we have used, we compare the spatial dependence of the interatomic components of the electron density and of the kinetic, potential and total energies at the equilibrium (R = 1.4 a0) and a stretched geometry (3.0 a0) in Figure S18. A numerical summary of the computed energies is given in Table S1. While the MO, VB and CI methods yield very similar results at R = 1.4 a0, at 3.0 a0 the MO density and potential

(34)

and Pinter(ζ = 1) = PH2(ζ = 1) − P2H(ζ = 1)

(35)

The molecular contraction results, as defined by eq 31, computed at the VB level at the internuclear distance of 1.4 a0, where ζopt = 1.1695, are displayed in Panel I of Figure 5. The results clearly indicate that the effects of orbital contraction on the density of the molecule and the corresponding kinetic, potential and total energy changes are greatest near the nuclei. As a conservative estimate of the atomic contributions to the energy changes of the H2 molecule on contraction we consider the nonatomic or “binding” region to be between z = −0.4 and 0.4. By excluding this region from the computation of energies we obtain the following atomic contributions to the kinetic, potential and total energies (in Eh, showing percentages of the total in parentheses): ΔT = 0.226 (75%), ΔV = −0.261 (78%), ΔVne = −0.279 (72%), ΔE = −0.035 (104%). Furthermore, the orbital contraction also results in an increase in the density in the interatomic region, computed to be 0.026e in the [−0.4, 0.4] region, as well as a small positive contribution of ∼0.001 Eh to the total energy. Thus, 75 and 78% of the kinetic and potential energy changes respectively that occur on contraction are estimated to be atomic in nature. The net result is a small drop in the total energy, almost an order of magnitude less than the decrease in the potential energy. By comparison, the effect of orbital contraction for two H atoms are ζ2 − 1 and −2(ζ + 1) for the kinetic and potential components of the binding energy respectively, i.e., 0.368 and −0.339 Eh, and hence an increase in the total energy by 0.029 Eh. The effect of orbital contraction on the interatomic components of the properties,Pinter(ζopt) and Pinter(ζ = 1), as defined in eqs 33-35 are shown in panel II of Figure 5. The plots convincingly show that orbital contraction results in a larger degree of density transfer to the binding region along with a corresponding uniform decrease in the kinetic energy. The changes in density are mirrored by changes in the potential energy (that is dominated by the attraction of the electrons to the nuclei), resulting in a net antibonding effect. The decrease in total energy, brought about by orbital contraction, is clearly I


Article


essentially the same results, inasmuch as the interatomic (i.e., interference) component of the kinetic energy is strongly bonding, while the interatomic potential energy is mildly antibonding, as noted by Schmidt et al.31,32 Moreover, the close agreement between the minimal basis CI, CASSCF and MRCI results indicates that accounting for polarization, by virtue of the extensive nature of the aug-cc-pVQZ basis, has a minor effect on the interatomic energies, as does dynamic electron correlation. The interatomic energies, as discussed above, are very similar to what had been calculated with a fixed value of 1 for the exponent ζ. The net result of orbital optimization is a lower total energy E, i.e., increased stability and shorter bond length: 1.414 a0 (obtained by VB). The corresponding equilibrium binding energy is −0.1391 Eh and its kinetic and potential components are 0.1391 and −0.2762 Eh respectively, clearly satisfying the virial theorem, eq 30. Figure 6 displays the various components of the potential energy of binding in H2. The effects of orbital expansion (at R >

energy are already quite different from what has been obtained using VB or CI, since the MO, i.e. Hartree−Fock type approach, does not account for static correlation. The ability of the electrons to move between the two atomic centers, even if tightly correlated, is clearly still the driving force for bonding. As a visual demonstration of the electron correlation, or its absence, we refer the reader to Figure S19. Thus, irrespective of the presence of correlation in the theory applied, the kinetic nature of the bonding mechanism is unaltered for a molecule with a pair of correlated electrons. The meaning of “delocalization” is refined by correlation but in a clear and plausible way so that the main features are readily seen to remain. 3.4. Ruedenberg Analysis of Interference Contributions to Covalent Bonding in H2. The quantitative bonding analysis, as proposed by Schmidt, Ivanic and Ruedenberg,31,32 starts with the notion of quasi atoms in the molecule which are bonded, as a result of electron sharing between them. In the simplest applications, such as explored here for H2, the quasi atoms are hydrogens whose wave functions are QUAOs. In the case of a complex multiconfiguration (correlated) wave function a transformation is required so that the molecular wave function is re-expressed in terms of subunits that exhibit atomic character. While the essential elements of this analysis have been published by Schmidt et al.,31,32 we present extra detail that we believe are useful in the context of covalent bonding, in particular the energetic consequences of interference, as illustrated in Figure 2. Using a minimal 1s QUAO basis (with exponent ζ) the kinetic and potential contributions to the binding energies of H2 at any internuclear distance R are written as sums of intraand interatomic contributions,31,32 as defined in eqs 7 and 8. Note, however, that the interatomic contributions to the potential energy, as defined by Schmidt et al.,31,32 are Vinter = Vqc + Vsc + VI + VII

Figure 6. Intra- and interatomic contributions to the potential energy component (Vbind) of the binding energy of H2 calculated using a VB wave function with with ζ optimized.

(36)

where Vqc is the quasi-classical Coulomb interaction, Vsc is the Coulombic sharing contribution31 or interpenetration associated with electron sharing (denoted Vpen),32 VI is the potential interference energy (with one- and two-electron contributions) and VII is the interelectron interference energy. The total binding energy is thus

2.85 a0) and contraction (at R < 2.85 a0) on the intra-atomic component, as discussed above, are very obvious given the small energy scale of this plot. All the interatomic components, with the exception of the quasi classical potential energy, are repulsive, although the VII interference contribution is very small. The VI component is the dominant interatomic term. It is almost entirely due to the one-electron interference component, VI(1), of the potential energy due to the charge movement into the bond, i.e.

E bind = (Tintra + Vintra) + (Tinter + Vinter) − 2(TH + VH) (37)

where TH and VH are the kinetic and potential energies of a hydrogen atom, 0.5 and −1.0 Eh respectively. The interatomic energies for H2, computed with a VB wave function, are shown in Figure 3, while a comparison of intraand interatomic energies is shown in Figure S20. In the region of R < 2.85 a0, ζ monotonically increases with decreasing R, as shown in Figure S4. In this region therefore Tintra and Vintra are monotonically increasing and decreasing functions respectively with decreasing R. On the other hand, for R > 2.85 a0 ζ decreases with increasing R to a minimum of 0.9901 at R = 3.5 a0 and then increases to its asymptotic value of 1.0 (reaching 0.9999 at R = 7.2 a0). The corresponding variations in Tintra and Vintra are hardly discernible though in the plots in Figure S20. A comparison of the interatomic energies computed by the minimal basis VB, CI and MO methods is shown in Figure S21, while in Figure S22 the CI results are compared with CASSCF and MRCI energies calculated with the aug-cc-pVQZ basis. These results demonstrate that qualitatively all methods yield

VI(1) = ⟨Ψ|Vâ + Vb̂ |Ψ⟩ + R−1 − Va(ζ ) − Vb(ζ ) − Vqc(1) (38)

∫ dr (Vâ + Vb̂ )ρI

=

(39)

The interference density (ρ I) associated with bond formation, as defined in eq 1, describes the transfer of charge into the bond, i.e. the internuclear region of the molecule, as shown in the contour map of Figure 2. The total charge transferred and the location of the zero ρI sheet are readily calculated by numerical integration as described above. Similarly, the effect of charge transfer on the kinetic and potential energies can be computed by integration over the regions of positive and negative ρI. The results of such calculations are shown in Figure 7. In the interatomic region where the density change is positive (red contours in Figure 1) J


Article


Figure 7. Charge transfer ρ+ into the bond of H2, and kinetic and nuclear attraction potential energy changes (T+, V+) in the region defined by ρ+ (full lines). Similarly, T− and V− are the analogous changes in the regions of negative density change ρ− (dashed lines). Panel I: VB wave function with ζ optimized. Panel II: Comparison of VB results with ζ = 1 (dotted lines) and ζ optimized (full and dashed lines).

analytic calculations on the H2+ molecule, that for R < ∼4 a0, orbital interference is enhanced (as indicated by an increase of the charge density in the bond) by orbital contraction. That, in turn, implies enhancement of the kinetic interference effect, i.e., a lower kinetic energy relative to atoms in contracted states. For R > ∼4 a0, orbital expansion, viz. a swelling, occurs, as indeed can be verified computationaly. Thus, orbital contraction and expansion are consequences of electron delocalization. Simple measures of the total electron density, as well as its interatomic (interference) component in the interatomic region, are the density ρm and its interference component ρIm at the bond midpoint. The R-dependence of these quantities and the effect of orbital optimization is shown in Figure 9. The

both kinetic (T+) and one-electron potential energies (V+) decrease. On the other hand, in the regions outside the interatomic region, where the density change is negative (blue contours in Figure 1), the potential energy (V−) increases, while the kinetic energy (T − ) decreases. The partial cancellation of V+ and V− results in a small positive VI, whereas TI (= T+ + T−) is negative and in magnitude roughly twice that of VI, as shown in Figure 7. The dominant contributions to the interatomic potential energy are the interference term VI and the quasi classical term Vqc, defined as Vqc = ⟨ϕa|Vb̂ |ϕa⟩ + ⟨ϕb|Vâ |ϕb⟩ + [ϕϕ |ϕ ϕ ] + R−1 a a b b

(40)

where the two-electron term represents the Coulomb interaction of the charge distributions |ϕa|2 and |ϕb|.2 The effects of orbital optimization on charge transfer and on the interference contributions to the kinetic and potential energies are obvious when the results in panel II of Figure 7 are compared. The distance dependence of the interatomic kinetic and potential terms and the dominant contributions to the latter are displayed in Figure 8. As remarked already, only the quasi-

Figure 9. Total density ρm (in red) and interference density ρIm (in black) at bond midpoint obtained using VB wave function with ζ optimized (full lines) and ζ = 1 (dotted lines).

plot indicates that as the internuclear distance decreases the density at the bond midpoint monotonically increases and once R < ∼2.5 a0, orbital contraction accelerates the process. Moreover, as noted by Schmidt et al.31,32 for H2+, the interference density is enhanced by orbital contraction, reaching a maximum value in the case of H2 at R ≅ 0.7 a0. Orbital expansion occurs for R > 2.85 a0 as shown in Figure S23 and both the charge at the bond midpoint as well as the density change due to interference are directly affected by the degree of orbital expansion, i.e. the magnitude of ζopt − 1. 3.5. Bonding in H2 Viewed in a Quantum Chemical Perspective. A simple, yet physically correct explanation of covalent bonding is to this day an issue that remains contentious. While we can routinely predict chemically accurate binding energies using the techniques of modern quantum

Figure 8. Interatomic contributions to the kinetic and potential energies of H2 computed with a VB wave function with ζ optimized.

classical component of the potential energy is attractive (in addition to the dominant intra-atomic part). Thus, at all distances the total interatomic potential energy is repulsive in contrast with the strongly attractive interatomic kinetic energy. The relationship between delocalization and orbital contraction has been put on firm theoretical footing by Schmidt, Ivanic, and Ruedenberg.31,32 These authors have shown, via K


Article


of contracted (or expanded) quasi atoms in the analysis of covalent bonding is, in our view, a totally reasonable process, especially in light of the fact that the ratios of the inter- and intra-atomic components of the kinetic and potential energies are independent of the scaling factor. Orbital contraction, as a distinct energetic component of the covalent bonding process, is not recognized in the electrostatic model. Polarization, i.e. the energetic consequence of including polarization functions in the basis set that yields polarized AOs in the construction of the wave functions, results in relatively minor changes in the interatomic energies. Thus, the essential elements of interference and hence bonding are readily recovered from the simple minimal basis calculations, such as the VB method. Similarly, dynamical correlation has only a minor effect on the magnitudes of the interatomic components of the energies. It is widely appreciated that even low level quantum chemical calculations are capable of accounting for chemical bonding, inasmuch as reasonably accurate molecular geometries and force constants are readily predictable by ab initio or semiempirical wave function as well as (Kohn−Sham) density functional methods44,45,68,92−95 using small, even minimal, such as STO-3G, basis sets. However, Thomas−Fermi (TF) theory,42,43 the original and simplest type of DFT,44,45,68 is unable to describe covalent bonding.46,47 Even though it contains all the electrostatic interaction terms between the electrons and nuclei, its application to molecules produces repulsive potential energy surfaces. This repulsion arises despite the complete flexibility of the density, which is free to reproduce the sort of variations that are generally associated with bond formation. This fundamental failure, encapsulated in the nonbinding theorem,44 is inherent in the mathematical structure of the TF kinetic energy functional, obtained by a statistically uniform sampling of phase space between bounding energy surfaces determined by the correspondence principle. It is hardly surprising therefore that the kinetic energy in modern DFT is evaluated by a wave function method, i.e., using MOs, as proposed by Kohn and Sham.68 The critical conclusion that emerges from the above observations is that the correct quantum mechanical treatment of kinetic energy is crucial for the description of covalent bonding. We note here also that the virial theorem holds in TF theory,96 yet it fails to predict covalent bonding. While the proper quantum mechanical treatment of kinetic energy is of crucial importance in the bonding mechanism, it does not imply that the energetic stabilization of a covalently bonded molecule is associated with a decrease in the total kinetic energy. As illustrated above, delocalization is an interatomic process related to wave function interference and causes a drop in the interatomic, i.e., molecular part of the kinetic energy, while the intra-atomic orbital contraction process results in an increase in the intra-atomic part of the kinetic energy as well as a decrease in the potential energy that is dominated by the electron−nucleus interaction. The end result for the molecule, as predicted by accurate quantum chemical calculations, is actually a net increase of kinetic energy when we compare the molecule at its equilibrium geometry with its constituent atoms, in agreement with the requirement of the virial theorem. However, covalent bonding is predicted by significantly less accurate calculations as well, where the kinetic-potential balance can be quite different from the requirements of the virial theorem, e.g., by neglecting orbital contraction. Thus, there are two processes happening as the

chemistry, there is still disagreement on the nature of bonding, even in the simplest of molecules, such as H2+ and H2. The simple Lewis theory of covalent bonds as related to shared electrons is however universally accepted and as such it represents the natural starting point for the interpretation of the quantum mechanical results. The resolution of the controversy between kinetic and electrostatic interpretations of covalent bonding is obviously an important pedagogical issue. It is also an important aspect of the progress of quantum chemistry since it relates directly to the early failure and recent success of the DFT methods and the search for ever better functionals for the prediction of ground state energies from orbital wave functions and the corresponding electron densities. The basic theoretical and computational approaches used in this work on the H2 molecule are essentially the same as those for H2+41 with the appropriate major extensions to accommodate electron repulsion and correlation, both static and dynamic. In particular, we have carried out an extensive spatial analysis of the changes in the integrated electron density and energy densities that occur on covalent bonding, in an effort to make a connection between the energetics of H2, as used in the work of Ruedenberg and co-workers,18,19,23−32,61 and the spatial density and energy partitioning which is implicit in the electrostatic models of bonding.8−11,52,53 As the latter rely on the electron density alone, the spatial dependence of the kinetic energy density has been largely ignored. Yet, our detailed investigation shows that it is the decrease in the interatomic, i.e., molecular, component of the kinetic energy which is responsible for the stabilization that occurs on delocalization, i.e., the basic mechanism of covalent bonding. Using the methodology recently developed by Schmidt et al.,31,32 we have shown that the interference of QUAOs not only results in an electronic density transfer into the bonding region of the molecule from the outside regions, but also in a decrease in the interatomic, i.e., interference, component of the kinetic energy in both regions. The movement of charge, however, results in a cancellation in the corresponding potential energy changes, so that the net result is actually antibonding. The above effects become more pronounced with increasing orbital contraction, as it enhances the degree of charge transfer and interatomic kinetic energy drop, while the effect on the interatomic potential energy remains quite small. This is a powerful demonstration that bonding is essentially due to electron sharing, i.e., delocalization, that results from the interference of the QUAOs, leading to a drop in the interatomic component of the kinetic energy as well as charge transfer to the bond. Orbital contraction, the concept that has proved most contentious to understand and interpret or indeed recognize, provides an extra degree of stabilization as well as resulting in a shorter bond length, i.e., a more compact molecule. As the process of orbital contraction or expansion is essentially the scaling of distances and wave function, it results in large changes in the kinetic and potential energies, but only a modest decrease in the total energy. It does not alter, however, the fundamental features of bonding, i.e., the interference components of the density and energies. Orbital contraction and its role in covalent bonding represent one of the main differences between the electrostatic and Ruedenberg’s models.18,19,23−32,61 In the latter, orbital contraction is regarded as a process distinct from electron sharing, although it is induced by sharing.31,32 Moreover, it is understood to be a primarily intra-atomic effect. The concept L


Article


(when carried out in a minimal basis), produces a potential energy surface for H2 where the dissociation products include H+ and H− along with H atoms, as shown in Figure S1. Such a calculation leads to a gross overestimation of the bond energy. This well-known “dissociation error” rules out the “independent electron approximation” for H2 except in the equilibrium region (and smaller bond lengths). This is the reason why in this work the strongly correlated VB wave function has been taken as the simplest consistent representation of covalent bonding in H2. The structure of the VB wave function is readily seen to imply that the motion of the two electrons between the atomic centers is fully correlated, i.e. the electron transfers are simultaneous and in opposite direction. The motion of the electrons is therefore regarded as opposite of “independent”. This is certainly true in the region where the nature of correlation is essentially static, i.e. at bond lengths larger than ∼2.5 a0, as discussed in section 5. At smaller distances the nature of correlation is mostly dynamic, i.e. short-range. Its role is the resolution of the Coulomb hole around each electron, i.e. to ensure that the wave function satisfies the electron cusp condition. While the VB and minimal basis CI calculations allow for a degree of left−right correlation, that is only a minor part of the total. (At R = 1.4 a0 the VB, CI and MRCI/aug-ccpVQZ correlation energies are −0.0109, −0.0196, and −0.0404 Eh, respectively.) In any case, as found in this work, as well as by Schmidt et al.31,32 the basic aspects of covalent bonding are not affected by dynamic correlation. Moreover, in light of the accuracy and reliability of Hartree−Fock theory in the equilibrium region, we can safely conclude that it captures the essential elements of covalent bonding in that region. The two-electron probabilities computed at the bond lengths of 1.4 and 3.0 a 0 (shown in Figure S19) provide a visual demonstration of the change in the degree of left−right electron correlation at the equilibrium geometry to that at the elongated bond length. We note however that Handy and Cohen97,98 have argued that the term “static” correlation should be replaced by the more precise “left−right” correlation. Moreover, the component of left−right correlation resolved by the VB wave function (as used in this work) is an exchange term in the DFT terminology, accurately resolved by the OPTX functional. Thus, irrespective of the presence of correlation in the theory applied, the dynamical bonding mechanism in H2 remains fundamentally the same as in H2+, where the motion of the electron is “independent” by definition. As our work clearly demonstrates, the kinetic nature of the bonding mechanism is unaltered for a molecule with a pair of correlated electrons. The meaning of “delocalization” is refined by correlation but in a clear and plausible way so that the main features are readily seen to remain. In summary, we emphasize that covalent bonding is critically dependent on the proper quantum mechanical treatment of kinetic energy as it is responsible for the correct description of the variation in dynamical constraints during the formation of a molecule. Electron correlation, such as between the pair of electrons in H2 in the VB or CI methods of computation, is a type of dynamical constraint always present to some degree in molecules with two or more electrons. It redefines the meaning of electron delocalization and limits its extent but leaves the fundamental driving force the same: molecules form due to the stabilization associated with facilitated interatomic electron motion. Both the quantum chemical energy decomposition and

molecule forms: electron delocalization and density contraction around the nuclei. The common factor in all calculations that successfully predict covalent bonding is the phenomenon of electron delocalization that will invariably lead to stabilization, while the correct kinetic−potential energy balance depends on both processes being accurately described that requires a more complex calculation. 3.6. The Dynamical View of Bonding in H2. From a fundamental quantum mechanical point of view, covalent bond formation can be interpreted as the result of the relaxation of dynamical constraints, which make the nonbonded atoms or molecular fragments dynamically decoupled or nonergodic.39,40 In H2+ for example, the important constraint is the spatial confinement of the electron to the vicinity of one of the protons when the internuclear distance is large. A potential barrier centered at the bond midpoint prevents the electron from transferring between the two Coulombic potential wells except by extremely slow tunneling. As the bond length decreases toward the equilibrium value the localization constraint is relaxed and the electron is able to more rapidly delocalize over the two protons. The covalent bonding mechanism in H2+ as well as other molecules is consistent with a more ergodic, or spatially delocalized, form of electron dynamics in the molecules than in the fragments that form them. Nonergodicity in quantum mechanics is manifested by the presence of degeneracies (or near-degeneracies) in the energy levels and in the localized character of the corresponding energy eigenfunctions of the Hamiltonian. The spherical symmetry of atoms in the one-electron picture results in the conservation of angular momentum, which in turn is related to the shell structure of atoms. The constraint of angular momentum conservation is notably absent in the Thomas− Fermi theory of atoms, which therefore lacks shell structure and atomic reactivity. As argued in our previous work,39−41 the absence of dynamical constraints in TF theory (due to its failure to describe kinetic energy correctly) means that it is unable to account for nonergodicity and hindered internal electron dynamics in atoms as well as in molecules. It is the relaxation of these constraints, resulting in rapid delocalized electronic motion, which allows for covalent bond formation. Thus, covalent bonding is fundamentally related to facilitated internal electron dynamics, i.e. to interatomic electron transfer in molecules. In the quantum dynamical picture of covalent bonding there is a significant difference between the mechanisms seen in H2+ and in H2. In H2, and most multielectron molecules, the motion of the electrons is “correlated”, i.e., not “independent” as in the one-electron H2+. Nevertheless, as demonstrated by a large number of quantum chemical calculations, Hartree−Fock theory is able to predict a range of molecular properties quite accurately, especially equilibrium geometries. The level of accuracy is significantly higher when modern forms of density functional theories are used. Both Hartree−Fock and Kohn− Sham DFT use the self-consistent field (SCF) approximation where the motion of the electrons is treated as being independent for the purpose of evaluating kinetic energy. Thus, both schemes rely on the covalent bonding mechanism being contained within or, at least, fully representable within, the “independent electron” approximation. However, the reality in a multielectron molecule, even in the simplest of them, i.e. H2, as explored in detail in this work, appears very different. A Hartree−Fock calculation, referred to as MO in this work M


Article


We conclude with the observation that 55 years ago Ruedenberg15 accurately captured the physical nature of the chemical bond, and his theory has more than withstood the test of time.

the quantum dynamical analysis thus reach the same conclusion: covalent bonding is the manifestation of electron delocalization.

■

4. CONCLUSIONS We studied and analyzed the covalent bond in H2 by three distinct quantum chemical approaches, concentrating in particular on the nature and importance of orbital contraction that has been the most contentious issue in the description of the physical nature of covalent bonds. Through a detailed study of the spatial characteristics of the kinetic and potential energy contributions to the binding energy of H2 we conclude that, as in H2+, the effect of orbital contraction on the molecular wave function is such that the resulting electron density and energy changes are largely restricted to the atomic regions. The above conclusions are in agreement with the results of a simple scaling analysis that shows that orbital optimization is equivalent to scaling the molecular wave function and energies so that the virial theorem is obeyed. As the ratios of the interand intra-atomic components of the kinetic and potential energies are independent of the scaling factor, we concluded that the concept of contracted (or expanded) quasi atoms in the analysis of covalent bonding is a totally reasonable process. Delocalization is the consequence of the interference between the QUAOs and the changes in density and kinetic energy that occur in the process of bonding are due to interference, i.e. due to the quantum mechanical nature of electrons. By contrast, interference has a small antibonding effect on the potential energy. Irrespective of the nature of the H2 wave functions or of the QUAOs that are used in their construction, the picture that consistently emerges is that the transfer of density into the internuclear region, while destabilizing with respect to the potential energy, results in a large decrease in the kinetic energy. Thus, accounting for orbital contraction and satisfying the virial theorem are important issues but not the basic elements of covalent bonding. The essential findings of this work, in that covalent bonding is a consequence of quantum mechanical interference and its effect on the kinetic energy of the molecule, contradict the electrostatic model that explains bonding in a diatomic molecule in terms of increased attraction of the electron density in the bonding region to two nuclei. The same conclusions were reached in a number of studies on the H2+ molecule,18,19,23−41 indicating that the basic physics of covalent bonding are perfectly discernible in the simplest of molecules. Thus, covalent bonding is essentially a one-electron phenomenon. Further, we have argued that covalent bonding in H2 is fundamentally related to electron dynamics, specifically to a relaxation of nonergodic constraints that are manifested in quantum mechanical interference and orbital delocalization that lead to a reduction of the “kinetic energy pressure”. Orbital contraction and electron correlation play important subsidiary roles, particularly at short and long bond lengths, respectively. However, it is our view that the issues raised by the virial theorem, especially in the equilibrium region, and the dissociation error implicit in the Hartree−Fock descriptions of molecules such as H2, should not be allowed to hide or confuse the underlying one-electron mechanism and kinetic character of the covalent bond.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.7b08963. Figures S1−S24 and Table S1 (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*(G.B.B.) E-mail: [email protected]. *(S.N.) E-mail: [email protected]. ORCID

George B. Bacskay: 0000-0001-5117-9440 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI), via INTERSECT Australia. NCI is supported by the Australian Government.

■

REFERENCES

(1) Burrau, Ø. Berechnung des Energiewertes des Wasserstoffmolekel-Ions (H2+) im Normalzustand. Det Kgl. Danske Videnskab Selskab. 1927, 7 (14), 1−18. (2) Condon, E. U. Wave Mechanics and the Normal State of the Hydrogen Molecule. Proc. Natl. Acad. Sci. U. S. A. 1927, 13, 466−470. (3) Heitler, W.; London, F. Wechselwirkung neutraler Atome und homöopolare Bindung nach der Quantenmechanik. Z. Phys. 1927, 44, 455−472. (4) Pauling, L. The Shared-Electron Chemical Bond. Proc. Natl. Acad. Sci. U. S. A. 1928, 14, 359−362. (5) Pauling, L. The Nature of the Chemical Bond, 3rd ed.; Cornell University Press: Ithaca, NY, 1960. (6) Shaik, S.; Hiberty, P. C. A Chemist’s Guide to Valence Bond Theory; John Wiley & Sons, Inc.: Hoboken, NJ, 2007. (7) Lewis, G. N. The Atom and the Molecule. J. Am. Chem. Soc. 1916, 38, 762−785. (8) Slater, J. C. The Virial and Molecular Structure. J. Chem. Phys. 1933, 1, 687−691. (9) Feynman, R. P. Forces in Molecules. Phys. Rev. 1939, 56, 340− 343. (10) Coulson, C. A. Valence, 2nd ed.; Oxford University Press: London, 1961. (11) Bader, R. F. W. Worlds Apart in Chemistry: A Personal Tribute to J. C. Slater. J. Phys. Chem. A 2011, 115, 12667−12676. (12) Bader, R. W. F. Atoms in Molecules: A Quantum Theory; Oxford University Press: New York, USA, 1994. (13) Ayers, P. W.; Boyd, R. J.; Bultinck, P.; Caffarel, M.; CarboDorca, R.; Causa, M.; Cioslowski, J.; Contreras-Garcia, J.; Cooper, D. L.; Coppens, P.; Gatti, C.; Grabowsky, S.; Lazzeretti, P.; Macchi, P.; Pendas, A. M.; Popelier, P. L. A.; Ruedenberg, K.; Rzepa, H.; Savin, A.; Sax, A.; Schwarz, W. H. E.; Shahbazian, S.; Silvi, B.; Sola, M.; Tsirelson, V. Six questions on topology in theoretical chemistry. Comput. Theor. Chem. 2015, 1053, 2−16. (14) Hellmann, H. Zur Rolle der kinetischen Elektronenenergie für die Zwischenatomaren Kräfte. Z. Phys. 1933, 85, 180−190. (15) Wang, S. C. The Problem of The Normal Hydrogen Molecule in The New Quantum Mechanics. Phys. Rev. 1928, 31, 579−586.

N


Article

The Journal of Physical Chemistry A (16) Rosen, N. The Normal State of The Hydrogen Molecule. Phys. Rev. 1931, 38, 2099−2114. (17) Slater, J. C. Quantum Theory of Matter, 2nd ed.; McGraw-Hill, Inc.: New York, 1968. (18) Ruedenberg, K. The Physical Nature of the Chemical Bond. Rev. Mod. Phys. 1962, 34, 326−376. (19) Gordon, M. S.; Jensen, J. H. Perspective on “The physical nature of the Chemical Bond. Theor. Chem. Acc. 2000, 103, 248−251. (20) Hirschfelder, J. O.; Kincaid, J. F. Application of the Virial Theorem to Approximate Molecular and Metallic Eigenfunctions. Phys. Rev. 1937, 52, 658−661. (21) Coulson, C. A.; Bell, R. Kinetic Energy, Potential Energy and Force in Molecule Formation. Trans. Faraday Soc. 1945, 41, 141−149. (22) Löwdin, P.-O. Scaling Problem, Virial Theorem, and Connected Relations in Quantum Mechanics. J. Mol. Spectrosc. 1959, 3, 46−66. (23) Feinberg, M. J.; Ruedenberg, K.; Mehler, E. The Origin of Binding and Antibinding in the Hydrogen Molecule Ion. Adv. Quantum Chem. 1970, 5, 27−98. (24) Feinberg, M. J.; Ruedenberg, K. Paradoxical Role of the KineticEnergy Operator in the Formation of the Covalent Bond. J. Chem. Phys. 1971, 54, 1495−1511. (25) Feinberg, M. J.; Ruedenberg, K. The Heteropolar One-Electron Bond. J. Chem. Phys. 1971, 55, 5804−5818. (26) Ruedenberg, K. The Nature of the Chemical Bond, an Energetic View. In Localization and Delocalization in Quantum Chemistry; Daudel, R., Ed.; Reidel Publishing Co.: Dordrecht, The Netherlands, 1975; Vol. 1, pp 222−245. (27) Ruedenberg, K.; Schmidt, M. W. Why Does Electron Sharing Lead to Covalent Bonding? A Variational Analysis. J. Comput. Chem. 2007, 28, 391−410. (28) Bitter, Th.; Ruedenberg, K.; Schwarz, W. H. E. Toward a Physical Understanding of Electron-Sharing Two-Center Bonds. I. General Aspects. J. Comput. Chem. 2007, 28, 411−422. (29) Ruedenberg, K.; Schmidt, M. W. Physical Understanding through Variational Reasoning: Electron Sharing and Covalent Bonding. J. Phys. Chem. A 2009, 113, 1954−1968. (30) Bitter, T.; Wang, S. G.; Ruedenberg, K.; Schwarz, W. H. E. Toward a Physical Understanding of Electron-Sharing Two-Center Bonds. II. Pseudo-Potential Based Analysis of Diatomic Molecules. Theor. Chem. Acc. 2010, 127, 237−257. (31) Schmidt, M. W.; Ivanic, J.; Ruedenberg, K. The Physical Origin of Covalent Bonding. In The Chemical Bond: Fundamental Aspects of Chemical Bonding; Frenking, G.; Shaik, S., Eds.; Wiley-VCH Verlag: Weinheim, Germany, 2014. (32) Schmidt, M. W.; Ivanic, J.; Ruedenberg, K. Covalent Bonds Are Created by the Drive of Electron Waves to Lower their Kinetic Energy through Expansion. J. Chem. Phys. 2014, 140, 204104. (33) Nordholm, S. Analysis of Covalent Bonding by Nonergodic Thomas-Fermi Theory. J. Chem. Phys. 1987, 86, 363−369. (34) Nordholm, S. Delocalization − The Key Concept of Covalent Bonding. J. Chem. Educ. 1988, 65, 581−584. (35) Bacskay, G. B.; Reimers, J. R.; Nordholm, S. The Mechanism of Covalent Bonding. J. Chem. Educ. 1997, 74, 1494−1502. (36) Eek, W.; Nordholm, S. Simple Analysis of Atomic Reactivity: Thomas-Fermi Theory with Nonergodicity and Gradient Correction. Theor. Chem. Acc. 2006, 115, 266−273. (37) Nordholm, S.; Bäck, A.; Bacskay, G. B. The Mechanism of Covalent Bonding: Analysis within the Hückel Model of Electronic Structure. J. Chem. Educ. 2007, 84, 1201−120310.1021/ed084p1201 and Supplementary Information: http://www.jce.divched.org. ezproxy2.library.usyd.edu.au/Journal/Issues/2007/Jul/abs1201. html#supplement. (38) Bacskay, G. B.; Eek, W.; Nordholm, S. Is Covalent Bonding a One-Electron Phenomenon? Analysis of a Simple Potential Model of Molecular Structure. Chem. Educator 2010, 15, 42−54. (39) Nordholm, S.; Eek, W. Ergodicity and Rapid Electron Delocalization - The Dynamical Mechanism of Atomic Reactivity and Covalent Bonding. Int. J. Quantum Chem. 2011, 111, 2072−2088.

(40) Nordholm, S.; Bacskay, G. B. The Role of Quantum Dynamics in Covalent Bonding − A Comparison of Thomas-Fermi and Hückel Models. In Advances in Quantum Theory; Cotaescu, I. I., Ed.; InTech: Rijeka, Croatia, 2012; pp 107−152. http://www.intechopen.com/ books/advances-in-quantum-theory/. (41) Bacskay, G. B.; Nordholm, S. Covalent Bonding: The Fundamental Role of the Kinetic Energy. J. Phys. Chem. A 2013, 117, 7946−7958. (42) Thomas, L. H. The Calculation of Atomic Fields. Math. Proc. Cambridge Philos. Soc. 1927, 23, 542−548. (43) Fermi, E. Un Metodo Statistico per la Determinazione di alcune Prioprietà dell’Atomo. Rend. Accad. Lincei 1927, 6, 602−607. (44) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1994. (45) Jones, R. O. Density Functional Theory: Its Origins, Rise to Prominence and Future. Rev. Mod. Phys. 2015, 87, 897−923. (46) Teller, E. On the Stability of Molecules in the Thomas-Fermi Theory. Rev. Mod. Phys. 1962, 34, 627−631. (47) Balázs, N. L. Formation of stable molecules within the statistical theory of atoms. Phys. Rev. 1967, 156, 42−47. (48) Kutzelnigg, W. The Physical Mechanism of the Chemical Bond. Angew. Chem., Int. Ed. Engl. 1973, 12, 546−562. (49) Kutzelnigg, W.; Schwarz, W. H. E. Formation of the Chemical Bond and Orbital Contraction. Phys. Rev. A: At., Mol., Opt. Phys. 1982, 26, 2361−2367. (50) Kutzelnigg, W. The Physical Origin of the Chemical Bond. In The Concept of the Chemical Bond; Maksic, Z. B., Ed.; Springer: Berlin, 1990; Vol. 2, pp 1−43. (51) Wilson, C. W.; Goddard, W. A. The Role of Kinetic Energy in Binding. I. The Nonclassical or Exchange Kinetic Energy. Theor. Chim. Acta 1972, 26, 195−210. (52) Bader, R. F. W. On the Non-Existence of Parallel Universes in Chemistry. Found. Chem. 2011, 13, 11−37. (53) Bader, R. F. W.; Hernandez-Trujillo, J.; Cortes-Guzman, F. Chemical Bonding: From Lewis to Atoms in Molecules. J. Comput. Chem. 2007, 28, 4−14. (54) Esterhuysen, C.; Frenking, G. The Nature of the Chemical Bond Revisited. An Energy Partitioning Analysis of Diatomic Molecules E2 (EN−Bi, F−I), CO and BF. Theor. Chem. Acc. 2004, 111, 381−389. (55) Frenking, G.; Krapp, A. Essay. Unicorns in the World of Chemical Bonding Models. J. Comput. Chem. 2007, 28, 15−24. (56) Morokuma, K. Molecular Orbital Studies of Hydrogen Bonds. III. CO···H-O Hydrogen Bond in H2CO···H2O and H2CO···2H2O. J. Chem. Phys. 1971, 55, 1236−1244. (57) Kitaura, K.; Morokuma, K. A New Energy Decomposition Scheme for Molecular Interactions within the Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 10, 325−340. (58) Ziegler, T.; Rauk, A. On the Calculation of Bonding Energies by the Hartree-Fock Slater method. I. The Transition State Method. Theoret. Chim. Acta 1977, 46, 1−10. (59) Levine, D. S.; Head-Gordon, M. Quantifying the Role of Orbital Contraction in Chemical Bonding. J. Phys. Chem. Lett. 2017, 8, 1967− 1972. (60) Levine, D. S.; Horn, P. R.; Mao, Y.; Head-Gordon, M. Variational Energy Decomposition Analysis of Chemical Bonding. 1. Spin-Pure Analysis of Single Bonds. J. Chem. Theory Comput. 2016, 12, 4812−4820. (61) West, A. C.; Schmidt, M. W.; Gordon, M. S.; Ruedenberg, K. Intrinsic Resolution of Molecular Electronic Wave Functions and Energies in Terms of Quasi-atoms and Their Interactions. J. Phys. Chem. A 2017, 121, 1086−1105. (62) Wu, Q.; Ayers, P. W.; Zhang, Y. K. Density-based energy decomposition analysis for intermolecular interactions with variationally determined intermediate state energies. J. Chem. Phys. 2009, 131, 164112. (63) Berlin, T. Binding Regions in Diatomic Molecules. J. Chem. Phys. 1951, 19, 208−213. (64) Hückel, E. Zur Quantentheorie der Doppelbindung. Z. Phys. 1930, 60, 423−456. O


Article

The Journal of Physical Chemistry A (65) Hückel, E. Quantentheoretische Beiträge zum Benzolproblem. Z. Phys. 1931, 70, 204−286. (66) Hückel, E. Quantentheoretische Beiträge zum Benzolproblem. Z. Phys. 1931, 72, 310−337. (67) Hückel, E. Quantentheoretische Beiträge zum Problem der Aromatischen und Ungesättigten Verbindungen. III. Z. Phys. 1932, 76, 628−648. (68) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (69) Weizsäcker, v. C. F. Zur Theorie der Kernmassen. Z. Phys. 1935, 96, 431−458. (70) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. S. A Complete Active Space SCF Method (CASSCF) Using a Density-Matrix Formulated Super-CI Approach. Chem. Phys. 1980, 48, 157−173. (71) Roos, B. O. The Complete Active Space Self-Consistent Field Method and its Applications in Electronic Structure Calculations. Adv. Chem. Phys. 1987, 69, 399−445. (72) Werner, H.-J.; Knowles, P. J. A second order multiconfiguration SCF procedure with optimimum convergence. J. Chem. Phys. 1985, 82, 5053−5063. (73) Knowles, P. J.; Werner, H.-J. An efficient second-order MCSCF method for long configuration expansions. Chem. Phys. Lett. 1985, 115, 259−267. (74) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (75) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. Electron Affinities of the First-Row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (76) Lu, W. C.; Wang, C. Z.; Schmidt, M. W.; Bytautas, L.; Ho, K. M.; Ruedenberg, K. Molecule intrinsic minimal basis sets. I. Exact resolution of ab initio optimized molecular orbitals in terms of deformed atomic minimal-basis orbitals. J. Chem. Phys. 2004, 120, 2629−2637. (77) Lu, W. C.; Wang, C. Z.; Schmidt, M. W.; Bytautas, L.; Ho, K. M.; Ruedenberg, K. Molecule intrinsic minimal basis sets. II. Bonding analyses for Si4H6 and Si2 to Si10. J. Chem. Phys. 2004, 120, 2638− 2651. (78) West, A. C.; Schmidt, M. W.; Gordon, M. S.; Ruedenberg, K. A comprehensive analysis of molecule-intrinsic quasi-atomic, bonding, and correlating orbitals. I. Hartree-Fock wave functions. J. Chem. Phys. 2013, 139, 234107. (79) Ehrhardt, C.; Ahlrichs, R. Population Analysis Based on Occupation Numbers II. Relationship between Shared Electron Numbers and Bond Energies and Characterization of Hypervalent Contributions. Theor. Chim. Acta 1985, 68, 231−245. (80) Werner, H.-J.; Knowles, P. J. An Efficient Internally Contracted Multiconfiguration-Reference Configuration Interaction Method. J. Chem. Phys. 1988, 89, 5803−5814. (81) Knowles, P. J.; Werner, H.-J. An Efficient Method for the Evaluation of Coupling Coefficients in Configuration Interaction Calculations. Chem. Phys. Lett. 1988, 145, 514−522. (82) Shamasundar, K. R.; Knizia, G.; Werner, H.-J. A New Internally Contracted Multi-Reference Configuration Interaction Method. J. Chem. Phys. 2011, 135, 054101. (83) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge Densities. Theoret. Chim. Acta 1977, 44, 129−138. (84) Davidson, E. R.; Chakravorty, S. A test of the Hirshfeld definition of atomic charges and moments. Theoret. Chim. Acta 1992, 83, 319−330. (85) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77. The Art of Scientific Computing, 2nd ed.; Cambridge University Press: New York, 1992. (86) Cohen, L. Local Kinetic Energy in Quantum Mechanics. J. Chem. Phys. 1984, 80, 4277−4279. (87) Cohen, L. Representable Local Kinetic Energy. J. Chem. Phys. 1979, 70, 788−789.

(88) Ayers, P. W.; Parr, R. G.; Nagy, A. Local Kinetic Energy and Local Temperature in the Density-Functional Theory of Electronic Structure. Int. J. Quantum Chem. 2002, 90, 309−326. (89) Ayers, P. W.; Parr, R. G. Sufficient Condition for Monotonic Electron Density Decay in Many-Electron Systems. Int. J. Quantum Chem. 2003, 95, 877−881. (90) Frisch, M. J.; et al. Gaussian 09, revision A.1.; Gaussian Inc.: Wallingford, CT, 2009. (91) MOLPRO, version 2015.1; a package of ab initio programs; Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. and et al., 2015. See http://www.molpro.net. (92) Becke, A. D. Density Functional Thermochemistry. 3. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (93) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (94) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200−1211. (95) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (96) March, N. H. The Virial Theorem in the Thomas-Fermi Theory. Philos. Mag. Series 7 1952, 43, 1042−1046. (97) Handy, N. C.; Cohen, A. J. Left-right correlation energy. Mol. Phys. 2001, 99, 403−412. (98) Handy, N. C. Exchange and dynamic correlation. Mol. Phys. 2009, 107, 721−726.

P


Covalent Bonding in the Hydrogen Molecule - ACS Publications

Recommend Documents