A one-electron perspective on dative and ionic bonding

Abstract

The distinction between electron-sharing and dative bonding remains a matter of continued debate more than a century after Lewis's foundational contributions to chemical bonding. Here we show that part of the apparent dichotomy may arise from an overemphasis on electron pairs as the fundamental units of bonding. Using a combination of Generalized Valence Bond (GVB) theory, real-space descriptors (Laplacian of the electron density, delocalization indices, and IQA energy decomposition), electron distribution functions, and variational quantum Monte Carlo Born maxima, we find that many bonds traditionally considered either dative or ionic can be consistently interpreted within a peculiar single-electron bonding framework. In these interactions, one electron of the nominal bonding pair remains largely localized on one fragment, whereas its spin-paired companion delocalizes between fragments, providing a dominant contribution to bonding. We illustrate this motif in prototypical ionic, polar covalent, and donor–acceptor systems. While this perspective does not fully resolve the conceptual distinctions between dative and ionic bonding, it offers a coherent interpretation across different bonding regimes and suggests that the proposed kind of single-electron bonding may be rather common and chemically relevant.

---

Surprisingly, the chemical bond, widely accepted as the basic building block of the edifice of chemistry, remains a subject of academic debate even 110 years after Lewis's seminal work.¹ Besides the electron-pair bond, Lewis² also introduced an electronic definition of acids and bases as acceptors or donors of pairs of electrons, respectively, and Sidgwick³ proposed an arrow notation for the donor–acceptor, known today as dative, bonds. Although the arrow was found not convenient by Pauling,⁴ Haaland⁵ noticed in 1989 that, regardless of how ethane and ammonia borane are formed, the molecules differ in the nature of the fragments obtained when the central bonds are broken. He proposed that, for dative bonds, minimum-energy rupture proceeds heterolytically, yielding either species with or without net spin and net charge. Also, many dative bonds tend to be considerably weaker and to display longer bond lengths than those typical of single bonds, a phenomenon that has not yet been conclusively explained. After Haaland's paper, the donor–acceptor perspective has been revitalized by many researchers to provide a wealth of new insights into main-group chemistry.⁶ Not unexpectedly, some have led to heated debates^7,8 that have slowly decayed over the years. For many, dative bonds are not different from electron-sharing interactions, for how can one tell where the electrons of a pair come from in a bound complex? For others, dative bonds display unique, distinct features.

We believe that the root of many of these disagreements runs deep, going back to the chemical imaginary of the electron pair. As has been emphasized many times (see, for example, a recent discussion in ref. 9), bonding is not directly related to the existence of a spin-coupled pair of electrons. In fact, the pervasiveness of the pair stems from Pauli's exclusion of same-spin electrons: the Fermi hole around an electron is leveraged by an opposite-spin counterpart, which fills it. Actually, the Fermi hole has been used to construct one of the most successful visual descriptors of pairing, the electron localization function (ELF) of Becke and Edgecombe.¹⁰ Since the pair is not the essence of the chemical glue, one naturally turns to single-electron bonding. Indeed, the H⁺₂ molecular ion¹¹ has served as a testbed for developing chemical bonding models for decades. Interestingly, the binding energy of H₂, 109.5 kcal·mol⁻¹, is less than twice that of H⁺₂, 64.4 kcal·mol⁻¹. This well-known fact, seldom discussed in basic chemistry courses, prompts us to reexamine the role of one-electron bonding. An even more striking case is that of Li₂ (D_e = 24.5 kcal mol⁻¹),¹² and Li⁺₂ (D_e = 30.0 kcal·mol⁻¹).¹³

In conventional computational chemistry, decoupling the role of the two electrons of a Lewis pair is not easy. For example, restricted Hartree–Fock (HF) and density functional theory (DFT) calculations assign the same one-electron function to a couple of opposite-spin electrons, which are thus effectively transparent to each other from the statistical point of view. Similarly, explicitly correlated wavefunction methods are usually built from the HF reference, resulting in two well-separated manifolds of highly and lowly occupied natural orbitals that are difficult to associate with specific electron pairs. As we will see, valence-bond theory (VB) can be used to remedy this situation, vide infra.

The central proposition advanced in this work is that the two electrons forming an otherwise conventional spin-coupled Lewis pair need not contribute equivalently to bonding. As will be shown, in certain ionic and dative compounds, only one member of the electron pair actively engages in bonding, giving rise to what we label as a single-electron bond. These situations must be distinguished from the canonical one-electron bonds found in species such as H⁺₂ or Ng⁺₂, in which only a single electron is effectively localized in the bonding region, resulting in spin-doublet ground states. By contrast, the delocalized electron discussed herein remains spin-paired with its partner electron at all times.

Evidence that this kind of one-electron bonds can be hidden in rather standard bonding patterns is easy to find. A 2007 Faraday Discussions study,¹⁴ which examined in detail the real-space evolution of the charge transfer process leading from Li + H to the LiH diatomic, showed that the distribution of electron counts in the two atoms (as described in real space by quantum chemical topology¹⁵ (QCT)) is overwhelmingly dominated by only two possibilities at all interatomic distances: either three electrons are found in Li and one in H, n_Li = 3, n_H = 1 (the neutral structure), or two and two, n_Li = 2, n_H = 2 (the ionic structure). At large distances the probability of the first distribution, p(3, 1) ≈ 1.00, while at equilibrium p(2, 2) ≈ 1 – p(3, 1) ≈ 0.90. Although not explicitly stated at that moment, this charge-transfer model implies that only one electron of the Li·+ H·pair participates in the chemical bonding of the species, whose binding energy is about 58.0 kcal mol⁻¹.¹⁶ The situation is similar in other ionic molecules. More recently,¹⁷ an equivalent pattern was found in the nominally dative bond present in LiBe⁺, where unlike in LiH, not only two electron distributions dominate, but also an instability toward an unrestricted solution was discovered. In this solution, one of the two electrons of the bonding pair remains as an unperturbed 2s-like Be function, while the other is described by a heavily hybridized orbital that delocalizes over the Li moiety.

Since no bonding exists without electron delocalization,¹⁸ the possibility that one electron of a pair remains localized over an atom or fragment while the other engages in more extensive delocalization raises the question of whether the weakness of many dative and ionic bonds is due to this phenomenon. In that case, both dative and several ionic bonds would be better described as single- or one-electron bonds, such that one electron of the pair is in an inactive (I), passive, non-participating state while the other is an active (A), participating in delocalization, state. Dative and ionic bonds would differ in their dissociation pathway. Upon rupture, the active electron in dative links would remain in the same fragment as the spectator electron. In ionic interactions, the delocalizing electron would end up in the other fragment, possibly after an avoided crossing. The dissociation of dihydrogen would imply an AA-to-II transition, wherein the two electrons participate in bonding, while the dissociation of LiH would be an AI-to-II transition, in which only one electron dominates the bonding landscape.

Here, we present compelling evidence that this is indeed the case. To this end, we use a battery of real-space descriptors: the Laplacian of the electron density,¹⁹ electron distribution functions (EDFs),²⁰ real-space bond orders,²¹ and the interacting quantum atoms (IQA) energy decompositions²² obtained from real-space atomic/fragment partitioning.¹⁵ At the same time, a lesser-known procedure that locates the most likely positions of all the electrons in a system by finding the maxima of the square of the wavefunction,^23,24 will also reveal additional independent information. EDFs are invariant under orbital transformations and are very stable with respect to the basis set. They also do not change substantially with the specific real-space partitioning method used, at least for the ionic and dative bonds under study. Born maxima and Laplacian plots, on the other hand, are fully method-independent. They can be obtained whenever a wavefunction or electron density is available, respectively. Thus, real-space evidence strongly supports our claims, and the existence of spectator/delocalized electrons is directly inferred, in the absence of any orbital representation, from the presence of only two non-negligible probability values in the two-electron EDFs. Our protocol provides a convenient framework for discovering or recognizing that some well-known pair-electron bonds are better understood as spectator-plus-active pair bonds than active-plus-active pair bonds.

To satisfy those interested in orbitals while maintaining simplicity, we also employ the Generalized Valence Bond (GVB) method of Goddard and coworkers.²⁵ In this ansatz, selected electron pairs are described by self-consistent, overlapping, singly occupied orbitals. GVB is possibly the simplest approach that allows one to monitor the function of each electron in a bonding pair individually. Since each function of a GVB pair is expanded in the full molecular basis set, the transition from localized to delocalized regimes can be easily observed. This is difficult with conventional valence-bond procedures, where one-electron functions are spanned within a subset of the full basis set. Note that GVB is able to correctly dissociate bonds. Its use provides visual clues to the otherwise neat results obtained from orbital-invariant real-space descriptors. The close agreement between the real-space and GVB outcomes demonstrates that our conclusions possess a high degree of generality and are not contingent upon a specific orbital representation.

Our work does not alter the traditional polar-covalent or donor–acceptor bonding mechanisms of the electron-pair bonds examined. However, it introduces a clear conceptual shift in the standard paradigm, where both electrons in a bonding pair are described by the same spatial function by construction. Thus, the spatial polarity of each electron is identical, and both electrons contribute equally to bonding. While this mean-field description can provide good energetics, it can be fundamentally flawed in other respects. Our analysis reveals different spatial polarities for each electron in a pair. This results in significant energy reduction and a framework that is conceptually distinct from and inaccessible to MO theory. This is not merely a reinterpretation of spatial polarization; it unveils the individual behavior of electrons hidden in pair-centric textbook representations.

The GVB ansatz, the real-space artillery, and the appearance of single-electron bonds

In this section, we describe the methodology using three basic examples while keeping the theoretical content to a minimum. Very succinctly, the GVB method involves replacing a limited set of the ϕ_i(1)α(1)ϕ_i(2)β(2) doubly occupied orbitals in a HF determinant, each describing an electron pair, with a pair of non-orthogonal valence bond functions (NOGVBs). These two NOGVBs are coupled to a singlet as (ϕ_ia(1)ϕ_ib(2) + ϕ_ib(1)ϕ_ia(2)) (α(1)β(2) – β(1)α(2)). Since the two non-orthogonal orbitals can be orthogonalized, this ansatz is equivalent to the linear combination (λϕ_iab(1)ϕ_iab(2) + µχ_iab(1)χ_iab(2)) (α(1)β(2) – β(1)α(2)), where ϕ_iab, χ_iab are two orthonormal functions. Typically, |λ|≫|µ| and the pair resembles a standard bonding (ϕ_ab)-antibonding (χ_ab) set. Any GVB expansion is equivalent to a complete active space (CAS) calculation with 2 electrons and 2 active orbitals per GVB pair. Since our goal is to discover new bonding motifs rather than to produce very accurate results, we just report in this work GVB calculations with one active pair, from which we build the one- and two-particle GVB densities, which are further analyzed with real-space techniques. The Born maxima will be computed from variational quantum Monte Carlo calculations with optimized Jastrow factors that start from the cc-pVTZ HF wavefunction, as described, for example, in ref. 24. This allows to separate correctly the positions of the electrons forming a pair.

A few notes of caution are in order. Firstly, GVB wavefunctions are strict eigenfunctions of 〈S²〉. The asymmetry of the two distinct NOGVB functions of a GVB pair is a consequence of the greater variational flexibility of the GVB ansatz compared to the single determinant approximation. A similar phenomenon occurs when the restriction of α and β electrons being described by the same spatial function is removed in unrestricted schemes. In some closed-shell molecules, it is possible to reach a lower energy by breaking spin symmetry. This is the case of H₂ at internuclear distances larger than the Coulson-Fischer point.²⁶ Secondly, Hund's rule does not apply to the two NOGVBs. Typically, the greater variational flexibility of the singlet space of the pair compared to that of the triplet situates the former below the latter in energy. Finally, the NOGVB asymmetry discussed here is fundamentally different from the naïve dissimilarity of the atomic orbital coefficients of polar-covalent bonds. In textbook MO theory, a doubly occupied bonding orbital between fragments A and B is written as ϕ ≈ c_aχ_a + c_bχ_b, where the square of these coefficients measures spatial polarity. This is the same for each electron in the pair. Although our systems engage in polar-covalent interactions, NOGVB asymmetry refers to different behavior of each electron in the pair, and the two functions ϕ_ia, ϕ_ib can each display their own particular spatial polarity.

Fig. 1 sketches the results of GVB calculations on the LiH, H₂, and HF molecules at their equilibrium geometries. Although these results have been known to the chemistry community for decades, their interpretation has not permeated enough. Beginning with dihydrogen, the Hartree-Fock description provides a doubly occupied symmetric 1σ_g orbital that contains two statistically independent electrons, each with a 50% probability of being found within each H atom. This yields²⁷ the binomial distribution p(1, 1) = 0.5, p(2, 0) = p(0, 2) = 0.25. The two H valence shells fuse together, and the Laplacian is strongly negative in the internuclear region. The maximum value of |Ψ|² occurs when one electron is located on each nucleus. When the two electrons are correlated at the GVB or CAS(2,2) level (the full configuration interaction wavefunction displays essentially the same behavior), the antibonding orbital of the GVB pair (equivalent to the low-occupation CAS natural orbital) has a population n* = 0.024. In the GVB description, the two NOGVBs, shown in the Figure, can be understood as broken-symmetry σ functions: one is more localized around the left nucleus and the other around the right, with an overlap S = 0.80. This description results in an energy decrease of −11.7 kcal mol⁻¹ with respect to the Hartree–Fock solution. As expected, electron correlation separates the electrons of the Lewis pair. At equilibrium, the Laplacian shows two minima at the nuclear positions, precisely where the electrons are located at the Born maximum. Sure enough, when the H atoms are taken to the dissociation limit, the two NOGVB localize further and further until they end up as 1s H functions.

Fig. 1 Top panel: orbital densities (arbitrary units, along the internuclear axis, z, in a.u.) of the two non-orthogonal functions forming the GVB pair in LiH (left), H₂ (middle), and HF (right) molecules at their cc-pVTZ//GVB equilibrium geometries (do not confuse the HF molecule with the Hartree–Fock approximation). Nuclei are marked with vertical lines, with the H atom in LiH and HF on the right. Medium panel: ± 0.1 a.u. Isosurfaces of the two GVB non-orthogonal orbitals. Bottom panel: ∇²ρ isocontours (positive in blue, negative in red), with local minima highlighted with a yellow sphere. Below the representation, the position of the α (white spheres) and β (pale-yellow spheres) of all valence electrons at the Born maximum are also indicated. Doubly occupied cores are represented by magenta (Li) or green (F) spheres. Finally, the probabilities of the dominant electron counts are also shown.

The chemical interpretation of the bond formation process in H₂ is that of two 1s electrons on each H atom that progressively delocalize over the other as the internuclear distance R decreases. Very importantly, at the HF level, delocalization is complete at any R, and each electron is described by a σ_g spinorbital. At the correlated level, however, delocalization is incomplete due to the presence of the other electron (correlation) and each particle behaves as an effective, only partially delocalized, electron described by an NOGVB. Bonding results from the independent delocalization of each of these correlated electrons, and the pair-electron bond is formed from two correlated one-electron bonds. These ideas are perfectly compatible with real-space analyses. The correlated EDF yields p(1, 1) = 0.586, p(2, 0) = p(0, 2) = 0.207, and a bond order DI = 0.828, which is smaller than the DI = 1.000 obtained at the HF level. This distribution can be reproduced from two effective electrons: the first with a probability p¹_a of being found in the left H_a atom and a probability p¹_b = 1 − p¹_a of being located in the right H_b atom. The second electron has probabilities p²_a and p²_b, of being found in H_a and H_b, respectively. Solving for the direct product that yields the full EDF (see Section S3 of the SI for details) yields p¹_a = p²_b = 0.708. Each of these effectively independent electrons contributes 2p(1 – p) to the bond order, where p is either p¹_a or p²_b. The bond order provided by each of these electrons obviously peaks at 0.5 when each p = 0.5, as occurs in a HF description. As we evolve towards dissociation and localize each of the electrons, p¹_a and p²_b approach 1, and the DI approaches 0. An IQA decomposition of the GVB/CAS wavefunction leads to an interatomic exchange-correlation energy V^ab_xc = −138.9 kcal mol⁻¹. This is the covalent component of the interatomic interaction energy. It vanishes at dissociation. When a bond forms from a pair of atoms, each of which provides an electron, the deviation of p¹_a,p²_b from 1 tells us about the independent delocalization of each constituent of the pair.

Now, we turn to the LiH and HF molecules. In LiH, a textbook example, the Laplacian shows how the valence shell of Li has disappeared. This leads to a strongly charged quasi-spherical Li atom (with QTAIM charge Q = 0.889) and a highly polarized H moiety. A recent discussion about whether this should or should not be called an ionic system can be found in ref. 28. The Born maximum shows two electrons on both the Li and H nuclei, and, as previously mentioned, the probability distribution is dominated by just two electron counts. Looking at the NOGVBs (mid-panel of Fig. 1) is illuminating. One of the orbitals is nearly a pristine spherical 1s H function. Unlike its H₂ counterpart, the other is also centered at the H and has a polarization pattern that resembles that of the K Laplacian shell of the H basin shown in the bottom-left panel of the Figure. The orthogonal antibonding GVB orbital has n* = 0.043, and the energy lowering of the pair is −11.04 kcal·mol⁻¹. When the molecule dissociates, the first NOGVB remains almost intact, while the second becomes the 2s Li orbital. This can be seen immediately as a pair formed by a localized spectator H electron and a 2s Li function that first delocalizes and is finally transferred to the H atom in the so-called harpoon mechanism.²⁹ This is basically a one-electron transfer and thus a one-electron bond. The GVB bond order of the system is 0.229, and the decomposition of the EDF shown in the Figure leads to p¹_a = 0.114 and p²_b = 0.991, where a is associated with the Li atom. The initial Li electron of the pair, with p¹_b = 1 − 0.114 = 0.886, has about a 90% probability of being found on the H atom, while the initial H electron, with p²_b = 0.991, does not delocalize at all. This results in 82% of the bond order coming from the delocalizing electron and a V^ab_xc value of −25.8 kcal·mol⁻¹.

A reverse, less extreme case is hydrogen fluoride, which exhibits a shared-shell Laplacian with three minima: two at the nuclei and one on the internuclear axis, which is closer to the F than to the H atom. The Born maximum displays a clear pair of opposite-spin electrons, one of which is at the H nucleus. The two NOGVBs resemble 2p_z F functions considerably, though one of them has a much more spherical lobe centered around the H atom. The rationale is similar. Upon dissociation, the upper NOGVB collapses onto a 1s H function, and the lower onto a pure 2p_z F orbital. Note again that it is one of the two NOGVBs that changes character upon bond formation; the other remains rather localized around its original center. This behavior is completely general in heteropolar bonds. The EDF decomposition yields p¹_a = 0.946, p²_b = 0.212, where a now refers to the F atom. The bond order is 0.438, and V^ab_xc = −91.6 kcal·mol⁻¹.

In these three examples, dissociation leads to full separation of the two electrons of the pair. However, this is not always the case.

Dative bonds

We have examined a set of prototypical systems based on the archetypal NH₃BH₃ molecule. Table 1 contains their binding energies, as well as the preparation energy required to deform the donor and acceptor fragments from their in vacuo to their in-the-molecule geometries. All of the systems satisfy Haaland's criterion; upon dissociation, the two electrons of the bonding pair remain on the donor fragment. Since it is well known that not all density functionals provide good descriptions of dative bonds,³⁰ we report M06-2X results together with GVB single-point calculations at the DFT equilibrium geometries. We also followed DFT dissociation paths by constraining donor(D)-acceptor(A) distances (e.g. the d(B–N) distance) while optimizing all other degrees of freedom. We will use ammonia borane and the CH₃NH⁺₃ cation to exemplify our results, then comment on the differences among the other systems.

Table 1 Preparation, E_p, and dissociation, E_b, energies of the dative systems under study. All data obtained at M06-2X/cc-pVTZ level, in the gas phase, and in kcal·mol⁻¹

System	E^Dp	E^Ap	Eb
BH3NH3	0.1	13.3	33.2
BH3NMe3	1.2	14.1	41.0
BMe3NH3	0.2	16.0	18.5
BF3NH3	0.2	23.1	24.9
CH3NH^+3	0.2	25.6	113.0
BH3PH3	1.7	10.3	24.0
BH3Cl^−	0.0	16.5	39.2
AlH3NH3	0.1	4.0	31.9
SiH3NH^+3	0.2	8.2	81.1

---

The DFT optimized geometry of BH₃NH₃ shows an interatomic N–B distance of 1.650 Å, in good agreement with the value obtained from microwave spectroscopy (1.658 Å).^5,31 Its computed dissociation energy, 33.2 kcal·mol⁻¹ also lies within the range reported by Haaland. For comparison, GVB provides a somewhat longer equilibrium distance, 1.692 Å. As previously mentioned, we consider a single GVB pair for the NH₃ – BH₃ σ bonding orbital at the DFT geometry. The occupation of the antibonding orthogonal GVB orbital is 0.013e, providing a pair stabilization (measured with respect to the HF reference) of −9.0 kcal·mol⁻¹. Fig. 2 shows the two NOGVBs, which display an overlap S = 0.849. A quick look at the plots reveals a familiar pattern now: one of the orbitals, labeled ϕ₁ in the following, is largely localized on the ammonia (D, donor) moiety, resembling its canonical lone-pair orbital. The other orbital, labeled as ϕ₂, delocalizes considerably more over the borane (A, acceptor) fragment, showing two maxima along the B–N axis. However, the maximum closer to the B atom is rather shallow. It is easy to justify the identification of function ϕ₁, which is directed along a bond axis and resembles a conventional σ-like function, with a lone pair of the ammonia fragment. Firstly, by not maintaining orthogonality constraints with ϕ₂, it decays smoothly towards the BH₃ fragment, showing no orthogonalization tails. Secondly, it can be quantitatively compared with the lone pair orbital of the free ammonia molecule, or with its equivalent in stretched ammonia borane, provided that the B–N distance is sufficiently large. Fig. S18 illustrates this latter situation: the overlaps, moments, and other quantitative orbital descriptors are remarkably similar in both cases.

Fig. 2 The BH₃NH₃ molecule at its DFT equilibrium geometry. First row:density of the NOGVB ϕ₁ (a), and ϕ₂ (b), along the B–N axis, denoted as z. Note that the donor atom (N) is placed at the bond coordinate z = 0.0, while the acceptor atom (B) is placed at the DFT equilibrium distance of 1.650 Å (nuclear positions are marked with vertical black lines). In the diagram for ϕ₂ in (b), the N–B internuclear region is zoomed in the inset to better assess the existence of two maxima. Second row: NOGVB |ϕ| = 0.1 a.u. Isosurfaces of ϕ₁ (c), and ϕ₂ (d). Third row: (e) snapshot of the ∇²ρ isocontours (with an electron density threshold of 0.001 a.u.) on a plane containing the B–N axis (positive/negative contours in blue/red), with the minima of ∇²ρ indicated as yellow spheres, and bond critical points as green spheres; and (f) positions of all α (yellow) and β (white) electrons at the Born maximum (hydrogens have been omitted, N nucleus in blue, B nucleus in pink).

Interestingly, the Laplacian of the DFT full density shows two minima located at 0.404 and 0.918 Å away from the N nucleus, much like in the ordinary CH₃–CH₃ bond. These minima lie within the N basin before the bond critical point (the small green sphere in the Figure). The Born maximum shows a lone pair around the N atom that points toward the boron atom in what can be understood as a direct mapping of the H₃N → BH₃ representation. The two electrons of this lone pair are located at similar, but not equal distances from the N atom, at positions very close to those of the main maxima of ϕ₁ and ϕ₂ along the B–N internuclear axis. The NOGVB orbitals thus provide a mapping of the Born maximum. A more covalent-like Born maximum where the first electron is located close to the internuclear maximum of ϕ₁ in Fig. 2a while the second is close to the secondary maximum of ϕ₂ in the inset of 2b is also possible.

QCT provides a DFT charge transfer of 0.075e from the donor to the acceptor, with a B–N bond order of 0.32 and V^BN_xc = −54.9 kcal·mol⁻¹. As previously reported,³² this system exhibits a considerable electron delocalization among the negatively charged H's of BH₃ and the N atom. This can be understood as π-backdonation or as an incipient case of collective bonding.³³ Regardless, the total donor–acceptor exchange-correlation energy is large, V^DA_xc = −110.9 kcal·mol⁻¹, and a significant electrostatic interaction V^DA_cl = −75.9 kcal·mol⁻¹ also exists, with a total DI(D, A) = 0.83. Only the BH₃ fragment leads to a relevant preparation energy due to the pyramidalization cost, about 13.3 kcal·mol⁻¹. Discounting this from E_b, V^BN_xc turns out to be, as usual, a reasonable measure of the corrected bond dissociation energy. When the GVB pair is separated over the QCT basins, the two-electron EDF provides p(2, 0) = 0.750, p(1, 1) = 0.240, p(0, 2) = 0.010 using the donor–acceptor (a ≡ D, b ≡ A) order. These values can be transformed into p¹_a = 0.958, p²_b = 0.217, which shows that one of the pair's electrons is primarily localized in the donor fragment while the other engages in more extensive delocalization. Note that these p values are reasonably similar to those found in the HF molecule.

Dissociation, modeled through a relaxed scan during which all the degrees of freedom are optimized except for the d(B–N) distance, proceeds towards full planarization of the BH₃ fragment. During this process, the ϕ₁ NOGVB undergoes minimal changes, while the ϕ₂ NOGVB localizes over the NH₃ fragment and the GVB pair collapses towards the configuration observed in free ammonia. A full account can be found in Section S8.1 of the SI. It is important to note that the maximum value of µ² occurs at a B–N distance of approximately 2.0 Å. Fig. 3 displays vividly how only ϕ₂ changes substantially throughout the process. As the N–B distance increases, ϕ²₂ passes from displaying one maximum, well separated from that of ϕ²₁, to two via a standard splitting catastrophe at d(N–B) ≈ 1.60 Å. The second maximum, represented by the red curve in the Figure, lies much closer to the N atom and the single maximum of ϕ²₁. At d(N–B) ≈ 1.65 Å, a second catastrophe destroys the first maximum of ϕ²₂, leaving only the second one, which quickly evolves together with the static maximum of ϕ²₁ marked by the black curve, towards the lone pair of free ammonia.

Fig. 3 Distance between the donor atom (N) and the maximum critical points (CP) of the NOGVBs densities as the N–B bond distance is stretched from 1.20 to 5.00 Å in BH₃NH₃. The black line shows the evolution of the unique maximum of ϕ²₁, also displayed as a small black sphere in the embedded molecular plot that depicts the BH₃NH₃ equilibrium situation. The red and blue lines show the positions of the maxima of ϕ²₂, also depicted as spheres of the same color in the embedded molecule. Grey dashed lines in the zoomed region delimit the range of N–B distances in which three maxima coexist. All distances are given in Å.

The CH₃NH₃⁺ system dissociates into a planar methyl cation and ammonia. However, as shown in Table 1, its binding energy is more than three times larger than that of ammonia borane, even though the pyramidalization cost of the CH⁺₃ acceptor is 12.3 kcal·mol⁻¹ larger than that of BH₃ in BH₃NH₃. The DFT equilibrium geometry shows an N–C distance of 1.505 Å, at which the occupation of the orthogonal antibonding GVB orbital is n* = 0.019e. This value is slightly more relevant to the σ-bonding between the CH₃–NH₃ fragments than it was in the BH₃–NH₃ case. This translates into a pair stabilization of −11.6 kcal·mol⁻¹ and a NOGVB orbital overlap of 0.823.

A look at Fig. 4 and 5 reveals a pattern similar to that in ammonia borane, though there are also clear differences. First, ∇²ρ at the C–N bond critical point (BCP) is now clearly negative (−0.461 a.u.) and displays two minima along the C–N axis. The minimum closest to the C atom lies very near to the BCP, slightly within the N basin. The first NOGVB, ϕ₁, resembles the ammonia lone pair, and its density displays just one maximum along the N–C line, 0.285 Å away from the N nucleus. ϕ₂ is a more symmetric σ function with a density displaying only one maximum, 0.593 Å from the C atom. In this sense, the cation's behavior is closer to H₂, yielding similar GVB stabilization and NOGVB orbital overlap.

Fig. 4 The CH₃NH⁺₃ molecule at the DFT equilibrium geometry. First row: density of the NOGVB ϕ₁ (a) and ϕ₂ (b) along the N–C axis, denoted as z. Note that the donor atom (N) is placed at the bond coordinate z = 0.0, while the acceptor atom (C) is placed at the DFT equilibrium distance of 1.505 Å (nuclei are marked with vertical black lines). Second row: NOGVB |ϕ| = 0.1 a.u. Isosurfaces of ϕ₁ (c) and ϕ₂ (d). Third row: (e) isocontours of ∇²ρ (with an electron density threshold of 0.001 a.u.) on a plane containing the N–C axis (positive/negative contours in blue/red, respectively) with the minima of ∇²ρ depicted as yellow spheres, and the bond critical points as green spheres; and (f) positions of all the valence α (yellow), and β (white) electrons at the Born maximum, with N and C cores depicted as grey and blue spheres, respectively (hydrogen atoms have been omitted).

Fig. 5 Distance between the donor atom (N) and the maxima critical points (CP) of the NOGVB densities as the N–C bond distance is stretched from 1.20 Å to 5.00 Å in the CH₃NH⁺₃ system. The black line shows how the unique maximum of ϕ²₁ evolve. This maximum is also represented as a black sphere in the inset image. The red and blue lines correspond to the positions of the maxima of ϕ²₂, also displayed as spheres of the same color in the embedded molecular plot, which corresponds to an interatomic N–C bond distance of 2.00 Å.

This is corroborated by delocalization and IQA energetic data. The very large donor–acceptor bond order, 1.14, is dominated by the C–N DI of 0.89. Similarly, the V^DA_xc = −183.8 kcal·mol⁻¹ value is also controlled by its N–C component, which is −161.6 kcal·mol⁻¹. Unlike ammonia borane, V^DA_cl = −3.6 kcal·mol⁻¹ is very small because the positive charge is almost equally shared among the fragments due to the large QCT charge transfer of 0.420 e from the donor to the acceptor. All this points toward a very strong covalency.

Even in this case, the signatures of one-electron bonding are easily identifiable. First, the bonding pair at the Born maximum exhibits a distinct lone-pair distribution. Second, when we separate the GVB pair across the QCT basins, as before, the two-electron EDF yields p(2, 0) = 0.491, p(1, 1) = 0.459, p(0, 2) = 0.050 with the same donor–acceptor (a, b) order. These values transform into p¹_a = 0.888 and p²_b = 0.447, reflecting a pair formed by one rather localized electron at the D fragment and one almost fully delocalized electron. We stress that even though p¹_a remains significant, there is appreaciable delocalization of ϕ₁ over the methyl fragment. Furthermore, the evolution of the NOGVBs upon dissociation reveals that the evolution of the ϕ₁ function is minimal. The ϕ₂ function, however, changes from an N–C σ orbital polarized towards the methyl group at equilibrium to an ammonia lone pair at dissociation.

Discussion

Building on the insights revealed in the previous section, we performed a series of geometry optimization and dissociation scans on the systems listed in Table 1. See Section S8 of the SI for complete details. Our previous findings remain unchanged, and the following paragraphs discuss them succinctly. First, Table 1 reveals that the geometrical preparation energy of the donor fragments (NH₃, Cl⁻, PH₃) is very small, while that of the acceptor fragments (BH₃, CH⁺₃, BF₃, AlH₃, SiH⁺₃) is considerably larger due to the pyramidalization cost, which peaks at about 26 kcal·mol⁻¹ for the methyl cation. Note that the binding energy of the fragments is not correlated with the preparation energy. In all cases, a ϕ₁ orbital that remains practically unchanged and considerably localized on D upon bond formation is found. Its ϕ₂ companion, on the contrary, delocalizes as the donor–acceptor (DA) distance decreases. Its square shows a consistent trend: one maximum close to D in the bonding region at large DA distances, which evolves into two maxima as the distance decreases, and which can even collapse into just one maximum located closer to A in strongly bound systems such as CH₃NH⁺₃.

A deeper understanding of the energetics of these systems is provided through the IQA partitioning, summarized in Table 2 and 3. While these data (see Section S4 in the SI) contain a wealth of information, we focus on the broad picture they provide.

Table 2 Deformation energies of donor and acceptor fragments (E^D_def and E^A_def, respectively) of the systems under study. All energies were obtained at M06-2X/cc-pVTZ level in the gas phase and are given in kcal·mol⁻¹

System	E^Ddef	E^Adef
BH3NH3	101.4	52.3
BH3NMe3	112.1	67.6
BMe3NH3	100.6	70.2
BF3NH3	96.3	69.5
CH3NH^+3	153.3	−79.0
BH3PH3	46.7	37.4
BH3Cl^−	97.0	28.5
AlH3NH3	99.4	49.3
SiH3NH^+3	61.2	49.3

---

Table 3 IQA Interaction energies between donor (D) and acceptor (A) fragments. The total interaction energy, E_int, as well as its electrostatic, V_cl and exchange-correlation contributions, V_xc, are shown in (kcal·mol⁻¹). Also shown are the D–A delocalization indices (δ^DA, in a.u.), and fragment QTAIM charges (Q_D, Q_A in a.u.). All calculations were performed at DFT equilibrium geometries

System	D	A	E^DAint	V^DAcl	V^DAxc	δ^DA	QD	QA
BH3NH3	NH3	BH3	−186.75	−75.88	−110.87	0.83	0.075	−0.075
	N	B	−550.54	−495.68	−54.86	0.32	−1.142	1.891
BH3NMe3	NMe3	BH3	−220.80	−86.39	−134.41	1.05	0.061	−0.062
	N	B	−578.33	−521.21	−57.12	0.32	−1.137	1.890
BMe3NH3	NH3	BMe3	−189.44	−73.76	−115.68	0.89	0.077	−0.077
	N	B	−541.24	−493.04	−48.21	0.28	−1.128	1.955
BF3NH3	NH3	BF3	−190.52	−90.92	−99.60	0.72	0.109	−0.109
	N	B	−669.35	−635.51	−33.84	0.19	−1.164	2.481
CH3NH^+3	NH3	CH3	−187.33	−3.58	−183.75	1.14	0.420	0.580
	N	C	−281.19	−119.55	−161.64	0.89	−1.013	0.289
BH3PH3	PH3	BH3	−172.73	−29.90	−142.83	1.17	0.113	−0.114
	P	B	247.85	317.87	−70.02	0.43	1.732	1.705
BH3Cl^−	Cl	BH3	−149.63	−45.47	−104.16	0.87	−0.804	−0.196
	Cl	B	−324.88	−267.30	−57.58	0.40	−0.804	1.842
AlH3NH3	NH3	AlH3	−115.90	−60.07	−55.83	0.50	0.021	−0.020
	N	Al	−518.88	−493.08	−25.80	0.18	−1.202	2.305
SiH3NH^+3	NH3	SiH3	−206.38	−120.58	−85.81	0.69	0.100	0.900
	N	Si	−791.32	−740.85	−50.47	0.32	−1.305	2.847

---

First, let us concentrate on the first three systems. They display rather different binding energies, which have traditionally been explained via the electron-donating inductive effect of the methyl groups. However, we observe that E^A_p does not follow the E_b trend. Examining the IQA deformation and interaction energies reveals that the inter-fragment interaction energies and mutual charge transfer and bond orders are similar in BH₃NH₃ and BMe₃NH₃. It is the large deformation of the BMe₃ group that can be tracked (Tables S2 and S4 in the SI) to a weakening of the B–Me links, which justifies the small binding energy of the latter system. Similarly, although the B–N covalency is slightly larger in BH₃NMe₃ than in BH₃NH₃, it is the 1,3 and 1,4 collective-like interactions that provide a considerably larger inter-fragment bond order and increased interaction energy in the former (see Table S6 in the SI).

We also noticed a clear correlation between V^DA_cl and the donor's deformation energy, which appears to be electrostatically driven. Similarly, ϕ²₂ displays two maxima in BH₃NH₃ and BH₃NMe₃, at 0.361, 0.807 and 0.415, 0.917 Å from the N nucleus, respectively, while it displays only one maximum, at 0.362 Å away from N in BMe₃NH₃. These results further confirm that the covalency of the B–N link is firmly tied to the extent of delocalization of ϕ₂.

In CH₃NH⁺₃, the even distribution of the cationic charge among the fragments yields a very large N–C bond order of about 0.90 and a smaller role of collective interactions. E^A_def is negative due to the extensive neutralization of the methyl cation in the dimer. This also justifies the large E^D_def value resulting from the concomitant loss of electron population. All this results in electrostatics playing a very small role in the final dimer, which is mainly covalently bonded. Recall that now ϕ²₂ has only one maximum close to the C nucleus, and that ϕ₁ is still more localized around N.

The energetic features of the other systems can be equally rationalized by a combination of deformations, collective interactions and dative bonds, which are mainly due to just one delocalization channel.

Conclusions

Our analysis indicates that a significant fraction of bonds conventionally labelled as dative or ionic can be fruitfully interpreted within a peculiar type of single-electron bonding picture, rather than solely as classical two-electron polar-covalent interactions. Across the systems examined, from LiH and HF to NX₃BY₃ complexes, the GVB description frequently reveals a “spectator” electron described by a strongly localized orbital spin-paired with a companion orbital whose delocalization provides the main bonding contribution. Real-space descriptors such as EDF decompositions, Laplacians of the electron density, delocalization indices, and IQA energy components provide mutually consistent evidence for this asymmetric electronic behaviour. This leads to a chemically intuitive picture in which dative and many ionic interactions can be viewed as arising from a dominant single-electron delocalization channel, in contrast with the more balanced two-electron delocalization typical of covalent shared-shell bonds.

This perspective also helps rationalize trends in donor–acceptor strength. Fragment deformation costs, collective interactions, and the spatial evolution of the predominantly delocalizing orbital play a central role in determining the final binding energy, while the more localized electron contributes less directly to bonding. Notably, even strongly bound species such as CH₃NH⁺₃ retain identifiable signatures of this asymmetric bonding pattern, despite their substantial covalent character.

Summarizing, adoption of the single-electron bonding point of view provides physical grounding for several long-standing observations. (i) Explanatory insight: our model rationalizes the Haaland criterion regarding the weakness and longer bond lengths of dative bonds. By demonstrating that only one electron of the pair effectively participates in delocalization, the half-pair nature of the bond explains its inherent reduced strength relative to standard single bonds, naturally mapped onto delocalization indices and exchange-correlation (covalent) energies. No satisfactory explanation for this phenomenon has been reported to our knowledge. (ii) Predictive insight: the framework provides a simple mathematical mapping of dissociation. In ammonia borane, the model predicts a topological splitting catastrophe in the GVB density as the bond stretches. MO theory lacks the mathematical vocabulary to describe this event, where the active orbital undergoes a sudden dramatic character change to collapse back onto the donor while the spectator remains pristine. (iii) Dissociation dynamics: the model dissects whether a bond breaks homolytically or heterolytically based on the real-space dynamics of the active versus spectator electron, allowing us to map the transition from localized to delocalized regimes with individual-electron resolution.

Taken together, these results provide a unified interpretative framework that connects diverse bonding motifs and highlights the relevance of electron-by-electron analysis. While not settling fully the broader conceptual debate on the nature of dative bonding, this approach offers a complementary viewpoint that may help clarify underlying electronic mechanisms and guide further discussion.

Computational details

All systems were optimized at the M06-2X³⁴ DFT level of theory using Gaussian16.³⁵ The M06-2X functional has been shown to be well suited for systems bearing dative bonds.^30,36 Generalized valence bond (GVB) calculations were performed with GAMESS.³⁷ In both DFT and GVB calculations, the cc-pVTZ basis set³⁸ was used (see the SI, Sections S1 and S2 for further details). Real space analyses were done with AIMAll³⁹ and Promolden.⁴⁰ EDFs were computed with the EDF program²⁰ at the GVB level. Born maxima for selected DFT geometries were calculated with AMOLQC⁴¹ through variational Quantum Monte Carlo (VMC) simulations in which a reference HF/cc-pVTZ wavefunction was optimized via the generalized power expansion of isotropic Jastrow factors (see the SI, Table S10).

Author contributions

A. M. P.: Funding acquisition, investigation, methodology, project administration, resources, supervision, validation, writing, review & editing. D. B-E.: data curation, formal analysis, investigation, methodology, writing – original draft. E. F.: data curation, formal analysis, investigation, methodology, software, supervision, validation, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting the investigation reported in this manuscript, including molecular geometries and input and relevant output files, have been deposited in the Zenodo repository (https://doi.org/10.5281/zenodo.18380848).

Supplementary information (SI): section S1, methodology; S2, basis set evaluation; S3, decomposition of the 2e electron distribution functions; S4, energy decomposition; S5, Jastrow factors; S6, electron distribution functions; S7, GVB features; S8, topology of the single-pair GVB and full DFT wavefunctions. See DOI: https://doi.org/10.1039/d6sc01276k.

Acknowledgements

The authors acknowledge the Spanish MICINN (Grant No. PID2024-155569NB-I00, https://doi.org/10.13039/501100011033) and the European Union “ERDF A way of making Europe” for financial support. D. B.-E. acknowledges the Spanish FICyT for a predoctoral grant (PA-23-BP22-168).

References

1. G. N. Lewis, J. Am. Chem. Soc., 1916, 38, 762–785.
2. G. N. Lewis, Valence and the Structure of Atoms and Molecules, Chemical Catalog Company, Incorporated, 1923.
3. N. V. Sidgwick, The electronic theory of valency, Oxford University Press, 1929.
4. L. Pauling, The nature of the chemical bond and the structure of molecules and crystals: an introduction to modern structural chemistry, Cornell university press, 1960, vol. 18.
5. A. Haaland, Angew Chem. Int. Ed. Engl., 1989, 28, 992–1007.
6. L. Zhao, M. Hermann, N. Holzmann and G. Frenking, Coord. Chem. Rev., 2017, 344, 163–204.
7. D. Himmel, I. Krossing and A. Schnepf, Angew Chem. Int. Ed. Engl., 2014, 53, 370–374.
8. G. Frenking, Angew Chem. Int. Ed. Engl., 2014, 53, 6040–6046.
9. L. Zhao, W. E. Schwarz and G. Frenking, Nat. Rev. Chem., 2019, 3, 35–47.
10. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397–5403.
11. K. Ruedenberg, Rev. Mod. Phys., 1962, 34, 326–376.
12. C. H. Wu, J. Chem. Phys., 1976, 65, 2040.
13. H. Bouzouita, C. Ghanmi and H. Berriche, J. Mol. Struct. THEOCHEM, 2006, 777, 75–80.
14. A. Martín Pendás, E. Francisco and M. A. Blanco, Faraday Discuss., 2007, 135, 423–438.
15. R. F. W. Bader, Atoms in Molecules, Oxford University Press, Oxford, 1990.
16. W. C. Stwalley and W. T. Zemke, J. Phys. Chem. Ref. Data, 1993, 22, 87–112.
17. A. Martín Pendás and E. Francisco, Nat. Commun., 2022, 13, 3327.
18. A. Martín Pendás and E. Francisco, Phys. Chem. Chem. Phys., 2022, 24, 639–652.
19. R. F. W. Bader, P. J. MacDougall and C. D. H. Lau, J. Am. Chem. Soc., 1984, 106, 1594–1605.
20. E. Francisco, A. Martín Pendás and M. A. Blanco, Comput. Phys. Commun., 2008, 178, 621–634.
21. C. Outeiral, M. A. Vincent, Á. Martín Pendás and P. L. A. Popelier, Chem. Sci., 2018, 9, 5517–5529.
22. M. A. Blanco, A. Martín Pendás and E. Francisco, J. Chem. Theory Comput., 2005, 1, 1096–1109.
23. A. Lüchow, J. Comput. Chem., 2014, 35, 854–864.
24. M. Menéndez-Herrero, E. Francisco and Á. Martín Pendás, J. Chem. Theory Comput., 2025, 21, 2448–2461.
25. F. W. Bobrowicz and I. W. A. Goddard, in The Self-Consistent Field Equations for Generalized Valence Bond and Open-Shell Hartree-Fock Wave Functions, ed. I. H. F. Schaefer, Plenum Press, New York, 1977, pp. 79–126.
26. C. A. Coulson and I. Fischer, Philos. Mag., 1949, 40, 386–393.
27. A. Martín Pendás, E. Francisco and M. A. Blanco, Phys. Chem. Chem. Phys., 2007, 9, 1087–1092.
28. A. Martín Pendás, D. J. L. Rodrigues and E. Francisco, Phys. Chem. Chem. Phys., 2025, 27, 1789–1793.
29. M. Rodríguez-Mayorga, E. Ramos-Cordoba, P. Salvador, M. Solà and E. Matito, Mol. Phys., 2015, 114, 1345–1355.
30. B. G. Janesko, J. Chem. Theory Comput., 2010, 6, 1825–1833.
31. L. R. Thorne, R. D. Suenram and F. J. Lovas, J. Chem. Phys., 1983, 78, 167–171.
32. A. Martín Pendás, E. Francisco and J. Casals-Sainz, Intermolecular Interactions in Crystals: Fundamentals of Crystal Engineering, ed. J. A. Novoa, Royal Society of Chemistry, London, 1st edn, 2018, ch. 6, pp. 178–221.
33. S. Sowlati-Hashjin, V. Šadek, S. Sadjadi, M. Karttunen, Á. Martín Pendás and C. Foroutan-Nejad, Nat. Commun., 2022, 13, 2069.
34. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
35. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian ∼16 Revision C.01, Gaussian Inc. 2016, Wallingford CT.
36. F. Huang, J. Jiang, M. Wen and Z.-X. Wang, J. Theor. Comput. Chem., 2014, 13, 1350074.
37. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis and J. A. Montgomery Jr, J. Comput. Chem., 1993, 14, 1347–1363.
38. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
39. T. A. Keith, AIMAll (Version 19.02.13), TK Gristmill Software, 2019, Overland Park KS, USA.
40. A. Martín Pendás and E. Francisco, PROMOLDEN. A QTAIM/IQA code, Available from the authors upon request at ampendas@uniovi.es and at the https://github.com/QTCOVI%20github%20repository.
41. A. Lüchow, S. Manten, C. Diedrich, A. Bande, T. C. Scott, A. Schwarz, R. Berner, R. Petz, A. Sturm, M. Hermsen, K. H. Mood, C. Schulte, L. Reuter, M. A. Heuer, J. Ludovicy, V. Terzi and M. V. Heinz, Amolqc (v8.0.0), 2022.

---