Observing the 3D chemical bond and its energy distribution in a projected space

TL;DR

This paper introduces a novel method for visualizing and analyzing chemical bonds and their energy distributions in a two-dimensional projected space, enabling detailed 3D bond analysis across different types of materials.

Contribution

The authors present a new approach to analyze chemical bonds in a condensed charge density, allowing for objective identification of bond regions and energy distribution in a projected space.

Findings

Bond silhouettes can be reverse-projected to reveal 3D bonding regions.

The method applies to both metallic and covalent bonds, showing formation and structure.

It reproduces known organic chemistry results from charge density data.

Abstract

Our curiosity-driven desire to "see" chemical bonds dates back at least one-hundred years, perhaps to antiquity. Sweeping improvements in the accuracy of measured and predicted electron charge densities, alongside our largely bondcentric understanding of molecules and materials, heighten this desire with means and significance. Here we present a method for analyzing chemical bonds and their energy distributions in a two-dimensional projected space called the condensed charge density. Bond "silhouettes" in the condensed charge density can be reverse-projected to reveal precise three-dimensional bonding regions we call bond bundles. We show that delocalized metallic bonds and organic covalent bonds alike can be objectively analyzed, the formation of bonds observed, and that the crystallographic structure of simple metals can be rationalized in terms of bond bundle structure. Our method also reproduces the expected results of organic chemistry, enabling the recontextualization of existing bond models from a charge density perspective.

Equations

$$T_G(\Omega )=\frac{{\hslash }^2}{4m}N{\int }_{\Omega }d\mathbf{r}\int d{\tau }^{\prime }\text{leavevmode}\nabla {\Psi }^{\ast }\cdot \nabla \Psi .$$

$$T_L(\Omega )=-\frac{{\hslash }^2}{4m}N{\int }_{\Omega }d\mathbf{r}\int d{\tau }^{\prime }\text{leavevmode}[\Psi {\nabla }^2{\Psi }^{\ast }+{\Psi }^{\ast }{\nabla }^2\Psi ].$$

$$T_L(\Omega )-T_G(\Omega )=-\frac{{\hslash }^2}{4m}N\int dS(\Omega ,\mathbf{r})\nabla \rho (\mathbf{r})\cdot \mathbf{n}(\mathbf{r}),$$

$$E(\Omega )=\frac{1}{2}V(\Omega )=-T(\Omega ).$$

$$\mathcal{F}[f_i]\equiv {\mathcal{F}}_i[\theta ,\varphi ,f_i(\theta ,\varphi ,s)]=\underset{{\mathrm{G}}_i(\theta ,\varphi )}{\int }{f}_i(s)dA(s)ds.$$

$$\bm{\tau}(p)=\frac{\partial{\rm{G}}_{i}(\theta,\phi,s)}{\partial s}\bigg{\rvert}_{p}\\$$

$$l_{1}(s_{i+1})=l_{1}(s_{i})\big{(}1+\kappa_{1}(s_{i})\Delta s\big{)},$$

$$l_1(s_n)=$$

$$=$$

$$l_{1}(s)=\displaystyle{d\theta\,s_{0}\exp\bigg{(}\int\limits_{s_{0}}^{s}\kappa_{1}(s)ds}\bigg{)}.$$

$$l_{2}(s)=\displaystyle{d\phi\,s_{0}\exp\bigg{(}\int\limits_{s_{0}}^{s}\kappa_{2}(s)ds}\bigg{)},$$

$$dA(s)=$$

$$=$$

$${\mathcal{V}}_i(\theta ,\varphi )=\underset{{\mathrm{G}}_i(\theta ,\varphi )}{\int }dA(s)ds$$

Full text

\dates

This manuscript was compiled on

\acresetall\acuse

ncp,bcp,rcp,ccp

1. 1 Introduction

2. 2 Charge density partitioning

3. 2.1 Condensed charge density space

4. 3 Bonds in molecules

5. 3.1 The topology of the condensed charge density

6. 3.2 Bond formation

7. 3.3 Spurious bond points and paths

8. 3.4 Missing bond points and paths

9. 4 Condensed charge density and local kinetic energy

10. 5 Metallic bonds

11. 6 Summary

1 Introduction

\dropcap

The great statistician George Box reputedly remarked, “All models are wrong but some are useful” [1]. Box’s observation has been referenced as relevant to scientific models in general and is particularly germane to the chemical sciences [2, 3, 4, 5, 6]. After all, much of chemistry relies on empirical models that have survived by proving useful to those creating new molecules and materials. Foremost among these are representations of chemical bonding that now undergird all of chemistry.

Bonding models are useful when providing a framework from which to estimate energy differences due to subtle changes in the arrangement or composition of an atomic system. Such useful models have a venerable history, arguably beginning with Gilbert Lewis’ century old insights regarding electron sharing [7], which later formed the kernel of the valence bond theory of Slater and Pauling [8]. The evolving perspectives of chemical bonding are almost too numerous to mention and constitute a significant portion of the chemical literature. However, the contributions to this corpus by such luminaries as Mulliken, Hammond, Coulson, Fukui, Hoffmann, Ruedenberg, Hückel, Goddard, Pople, Parr, Peyerimhoff, Karplus, Levitt, and Warshel [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22] demand recognition.

Despite the advantages current bonding models confer, they have proven difficult to apply broadly, that is, to all the stuff held together by “bonds.” For example, to metals and alloys where suitable chemical bonding models might prove just as useful as those that have been employed in the design and synthesis of organic polymers [23].

There are forces at work that may lead to a change in this situation. Advances in computational density functional methods [24] coupled to ever more accurate measurements of the electronic charge density [25] are providing the impetus for the development of charge density based bonding models [26]. Because the electronic charge density is an observable—existing independently from the methods used for its calculation or measurement—such models should prove useful across all classes of molecules and solids. An observable chemical bond in particular could similarly be used to extend existing bonding models to new fields by connecting them with their roots in the charge density.

Foremost among the efforts to frame chemical principles around the charge density is Bader’s \acqtaim [27, 28]. With its clearly constructed formalism through which to identify an atom’s boundary, \acqtaim brings clarity and consistency to a number of previously ill-defined chemical concepts such as the energies, sizes and electron counts of the atoms comprising a molecule or solid. However, the topological representation of chemical bonding ensuing from this theory is plagued by an ongoing debate—both questioning [29, 30, 6, 31, 32] and supporting [33, 34, 35, 36, 37, 38, 39] its rigor and utility. It should be noted, however, that the topological approach to bond analysis due to \acqtaim neglects the central attribute of the theory—the partitioning of charge density into regions with well-characterized energies. We show here that this omission is not intrinsic to \acqtaim-based bond analysis approaches.

We proceed by reviewing the consequences and rationale behind charge density partitioning. Building on this rationale, we define a maximally partitioned charge density space—the condensed charge density—in which the kinetic energy resulting from the Laplacian family of kinetic energy operators is everywhere well-defined [40]. While faithfully recovering the essential elements of \acqtaim, the topology of this space reveals a charge density volume with the properties of a chemical bond, which is designated a bond bundle. We then examine several instances in which QTAIM’s topological bond has been asserted to be faulty and demonstrate that these assertions are made moot by the bond bundle construct. We additionally argue that the condensed charge density space is ideally suited to describe chemical phenomena as its structure derives from a preferred moving coordinate frame giving primacy to charge density isosurfaces. Supported by these findings and arguments, we apply our approach to some simple metals and demonstrate that their crystallographic structure can be rationalized as a consequence of bond bundle structure.

2 Charge density partitioning

The significance of the electronic kinetic energy as a mediator of chemical bonding was recognized as early as the 1930s, for example by Hellman [41] and Slater [42]. The central role of the kinetic energy took on further chemical import in Ruedenberg’s classic 1962 paper, The Physical Nature of the Chemical Bond [14], which prompted subsequent efforts to capture changes to the local kinetic energy as necessary for useful theories of bonding. These efforts were confounded, however, by the the multiple representations for local kinetic energy [43, 40]. Among the forms commonly used is one referred to here as the gradient representation, $T_{G}$ [44, 45], in which the total kinetic energy over a region $\Omega$ appears as,

[TABLE]

An alternative form for the kinetic energy of the same region may be expressed in terms of the Schrödinger (Laplacian) kinetic energy [46] as,

[TABLE]

In general $T_{L}(\Omega)$ and $T_{G}(\Omega)$ differ. It is straightforward to show [43] that

[TABLE]

where $\nabla\rho(\bm{r})$ is the gradient of the charge density, $S(\Omega,\bm{r})$ is the surface bounding the region $\Omega$ and $\bm{n}(\bm{r})$ is the unit vector normal to this surface at $\bm{r}$. The integral gives the net flux of the charge density gradient through $S(\Omega,\bm{r})$. Where this flux is zero, the two kinetic energy representations give the same value—as is required if $\Omega$ spans an entire molecule or the unit cell of a periodic solid. Thus, the kinetic energy of a molecule or a crystal’s unit cell is unambiguous. However, smaller regions contained in a molecule or unit cell may also possess well-defined energies if the net flux of the gradient of the charge density is zero over their boundaries [47], which necessarily includes regions bounded by surfaces over which the flux of $\nabla\rho(\bm{r})$ is everywhere zero. Such regions are said to be bounded by zero flux surfaces and it was believed that each possess a well-defined kinetic energy [27].111 Anderson et al. [40] have demonstrated that the kinetic energy ambiguity is much broader than embodied in Eqs. 1 and 2. We will address the broader implications of kinetic energy ambiguity in Section 4 of this paper.

Bader noted a unique zero flux surface surrounding the nucleus of every atom of a molecule, making it possible to associate a kinetic, potential and hence total energy ($T(\Omega)$, $V(\Omega$), $E(\Omega)$ respectively) with an atom in a molecule. Consistent with chemical models that assume energies to be additive, the molecular energy is given as the sum of these atomic energies. In addition to their well-characterized energies, Bader showed that these atomic regions (actually all zero flux surface bounded regions) satisfy the virial theorem [27]. Accordingly, for molecules or solids in structural equilibrium or at a transition state [42, 48],

[TABLE]

The atoms in a molecule (or solid)—called Bader atoms or atomic basins—owe their energetic significance to the zero flux surfaces that bound and define them. Obviously, The \aclqtaim derives its name from this fundamental association.

Explicit atomic boundaries impose a connectivity on the atoms of a molecule. Specifically, two atomic basins that touch over some finite region of their mutual boundary may be connected by an observable ridge of maximum charge density that extends from one atomic nucleus to the other. Bader called this ridge a “bond path.” A necessary condition for the existence of such a path is a charge density saddle point of index $-1$ on their shared boundary, which Bader named a \acbcp.

For the molecules Bader originally investigated, their sets of bond paths recovered Lewis’ molecular graphs. However, subsequent and ongoing investigations (see for example References 30 and 33) have recovered bond paths and \acpbcp between atoms that the Lewis model does not predict to be “bound.” In some instances bond paths have been discovered between atoms that are argued to be repulsive [37]. And in still other instances bond paths are absent between atoms that chemical reasoning predicts to be bound [32, 49, 50, 51, 52]. As mentioned, such revelations have prompted questions regarding the extent to which a bond path should be taken as an indicator of a bond and the indicated interaction as “bonding" (see particularly References 29, 30 and 31).

Pendás et al. used a wave function-based interacting quantum atoms approach to investigate the selection mechanism for bond paths. They showed that bond paths indicate the channels along which quantum exchange energy is most lowered [39]. A bond path may still correspond to stable, unstable, or metastable interactions depending on the role played by exchange, and their work offers a clue as to why bond paths may be found in some instances but not others.

The larger debate, however, assumes that connections between the atoms of a molecule must be stabilizing with respect to total energy. For the vast majority of organic molecules, which conform to the Lewis model, stabilizing effects proceed through electron sharing and hence connectivity. However, we see no a priori reason that this must be the case and believe that structures giving rise to stability must be embedded in a space in which energy is everywhere well defined. This is the space of all charge density volumes bound by zero flux surfaces called the \acP.

2.1 Condensed charge density space

As we summarized in an earlier article [53], the space \acP is constructed through a mapping of \aclpgp \acresetgp of \acrho to points in \acP. Every \acgp in \acrho originates from a local minimum called a \acccp and terminates at a maximum, almost always coincident with a nucleus and hence denoted a \acncp. We assume these \acpgp to be parameterized by arc length, $s$.

Sufficiently close to their terminus \acpgp are radial, making it convenient to imagine every \acncp as the center of a sphere Si of radius $dr$ (in practice $dr\approx 0.2$Å). Passing through every point on the surface of these spheres is a \acgp. The points on such a sphere may be specified by a polar and an azimuthal angle, so each of the molecule’s \acpgp may be specified by a pair of coordinates and the index of the \acncp at its terminus, i.e. \acgp${}_{i}(\theta,\phi)$. In this reference system, the charge density of the atomic basin $i$ is a function of $\theta,\phi$ and $s$, where $\theta$ and $\phi$ pick out a unique gradient path in the atomic basin $i$, and $s$ is the euclidean distance along this path.

Imagine covering each Si with a set of nonintersecting differential elements of area $dA=dr^{2}\,d\theta\,d\phi$ (Fig. 1). The \acpgp passing through the points interior to each of these area elements give rise to a family of infinitesimal volume elements called \acusedgb\aclpdgb, \acadgb [54, 55], each of which is bounded by a zero flux surface. Significantly (see Section 4), the cross sectional areas of these \acadgb (Fig. 2) change throughout their length and thus $dA$ is a function of $s$, $\theta$, and $\phi$. The union of all \acdgbi is equivalent to the union of all \acpgp terminating at \acncp $i$ and hence recovers Bader’s atomic basin. And significantly, these are the smallest structures bounded by zero flux surfaces and accordingly possessing well-defined energies and energy related properties.

Explicitly, for any scalar field, $f_{i}$, there exists a corresponding condensed scalar field, ${\cal F}_{i}$, that is a function of $\theta$ and $\phi$ and a functional of $f_{i}$, such that

[TABLE]

In particular, the charge density yields the condensed charge density (${\cal F}[\rho]={\cal P}$), the gradient or Laplacian kinetic energy densities yield the condensed kinetic energy density (${\cal F}[T_{G}]={\cal F}[T_{L}]={\cal T}$), and so forth. For the special case where $f_{i}(\theta,\phi,s)=1$, the \aclV is produced (${\cal F}[1]={\cal V}$).222In open systems a step function—defined to be one inside some charge density isosurface (commonly taken to be $\rho=0.001\leavevmode\nobreak\ au$) and zero outside this region—is used to obtain \acV.

As an illustrative example, Fig. 3 depicts \acP for each of the $i$ atoms of vinyl alcohol, where every point on the sphere surrounding an atom is colored to represent the magnitude of the integrated electron density in the \acdgb originating at that point. Borrowing terminology from differential geometry, each \acPi is referred to as an atomic chart and the set of all atomic charts comprising a molecule is termed its molecular atlas.

As an alternative to representing atomic charts on spheres, they may be projected onto a flat space as in Fig. 4 where the alcohol carbon (C1) atomic chart is represented with a stereographic projection centered at the C = C bond path.

By construction, every point in \acP maps to a \acgp in \acrho. Hence every trajectory through \acP maps to a zero flux surface in \acrho, and any closed loop in \acP maps to a volume in \acrho bounded by a zero flux surface and thus characterized by a well-defined energy. Such volumes are called gradient bundles [54, 55, 56]. All previously noted zero flux surface-bounded volumes, e.g. atomic basins, are proper subsets of the space of gradient bundles.

3 Bonds in molecules

3.1 The topology of the condensed charge density

Maxima in \acP typically map to bond paths in \acrho, as is the case for vinyl alcohol and evident in Fig. 3 where each of the molecule’s bond paths (shown as black lines) intersects an atomic chart at a maximum in \acP.

Just as all \acpgp terminating at the same maximum in \acrho delineate an atomic basin as a unique volume, \acfPgp terminating at the same maximum define a unique two-dimensional basin in \acP and hence a unique zero flux surface-bounded volume in \acrho.

Continuing with the vinyl alcohol example, Fig. 5(a) shows stereographic projections of the C1 atomic chart centered on each of its three maxima. Also shown are the three sets of \acspPgp terminating at each maximum and delineating corresponding \acP-basins. The image of all \acP-basins of the molecular atlas partitions \acrho into space-filling regions bounded by zero flux surfaces, each containing a bond path. The energy of each region is well defined and additive to give the molecular energy. These are the characteristics associated with a chemical bond.

Accordingly, we offer the following definitions: i) A bond wedge is the image in \acrho of \aclPgp terminating at a common maximum; and ii) a bond bundle is the union of bond wedges sharing a common intersection with an interatomic surface.333In the vast majority of cases the general definition provided here recovers the same regions as those resulting from an earlier bond bundle definition [57, 58]. In the case of lone electron pairs, the terms lone pair wedge and lone pair bundle can be used. Some of the bond bundles and lone pair wedges of vinyl alcohol are shown in Fig. 6.

Bond wedges and bond bundles are attractive in their one-to-one mapping with organic bonds as understood through valence bond theory, but they are merely the demarcation of charge density made possible through a gradient bundle decomposition of the system. The condensed charge and energy distribution within a bond wedge contains vastly more information, and we anticipate that these distributions will be significant in discovering the charge density structure underlying subtle differences between systems otherwise thought to operate via the same chemical mechanisms.

3.2 Bond formation

Though maxima in \acP typically map to bond paths, this is not always the case—a feature that allows one to “see” bond formation during chemical reactions and to assess the relative energies associated with such processes. Consider for example the Lewis acid base reaction between borane and ammonia,

\ce

BH3 + $:$NH3 -> H3B-NH3

Fig. 7(a) depicts the bond bundles and bond wedges of \ceNH3 and \ceBH3 alongside their respective \acpP.

Ammonia’s nitrogen atom is distinguished by four maxima in \acP and hence four bond wedges. Three of the these share bond wedges with those from hydrogen atoms, yielding three N – H bond bundles. The bond wedge centered on what we identify as the nitrogen atom’s south pole does not share an interatomic surface with another atom, as one would expect of a lone electron pair. Borane, on the other hand, is characterized by three shared bond wedges with hydrogen atoms yielding three B – H bond bundles. These bond wedges intersect along gradient paths that map to minima in \acP located at the boron atomic poles.

The reaction between borane and ammonia molecules was simulated by aligning the south pole of the nitrogen atom with the north pole of the boron atom at an initial B – N distance of 3.0Å. The reaction was allowed to proceed to the equilibrium B – N distance of 1.65Å, forming the B – N bond bundle depicted in Fig. 7(b).

The evolution of the bond bundle and its associated condensed properties along the reaction profile are depicted in Fig. 8. Inspection of this figure reveals that at 3.0Å separation, while there is a \acbcp and a bond path connecting the boron and nitrogen \acpncp, the two atoms do not contribute to a common bond bundle because the B – H bond wedges account for the entirety of \acP on the boron. The lone pair bond wedge on the nitrogen atom persists but there is no corresponding bond wedge on boron, which, at its pole, instead hosts a \acP minimum—as well as minima for both \acV, and \acT.

Bond bundle formation begins at a B – N distance of approximately 2.70Å with the near simultaneous development of maxima in \acP, \acV, and \acT at the boron pole. The maximum in \acP necessitates a bond wedge on the boron atom and an accompanying B – N bond bundle that continues to evolve over the course of the reaction to encompass a greater portion of the interatomic region between the boron and nitrogen atoms.

Across all calculations we find a correlation between \acP, \acV, and \acT, hinting at an energy-charge density structure relationship we will discuss in more detail in Section 4. However, for now recall that when the forces acting on the molecules are small, i.e. when far apart and at the equilibrium separation, the virial theorem asserts that $E(\Omega)=-T(\Omega)$. Taking $\Omega$ to be the B – N bond bundle, in principle the energy of the region corresponding to B – N bond is given by the integral of \acT over the appropriate bond wedges on the boron and nitrogen atoms. Without evaluating the integral, however, it is clear from Fig. 8 that the contribution to the total energy from the boron bond wedge is stabilizing, as one would expect. Moreover, the simultaneous formation of basins in \acP and \acT emphasizes that a bond wedge is also an energetic stabilization basin, in agreement with the common conception of bonding as the concentration of charge that results in stabilization of the system.

3.3 Spurious bond points and paths

Another example in which the bond bundle picture provides more information than the conventional \acqtaim representation is provided through an analysis of “hydrogen-hydrogen” bonding, as distinct from hydrogen or dihydrogen bonding [33]. Two cases are considered here: i) the planar biphenyl molecule (Fig. 9) characterized by \acpbcp between ortho hydrogen atoms that are not found in its twisted, lower energy conformation, and ii) dibenzene (Fig. 10), also displaying unexpected H … H bond paths.

The first instance, planar biphenyl, has been thoroughly investigated. Some have argued that the H … H interaction is destabilizing and hence the occurrence of a bond path between the two hydrogen atoms is a failure of \acqtaim [6, 30]. While others assert that the H … H bond path in fact lowers the energy of the metastable configuration [33, 36, 59, 39].

When analyzed from the bond bundle perspective, these arguments become moot. The H … H bond paths of biphenyl (Fig. 9) do not map to maxima in \acP. These bond paths lie within the ortho C – H bond bundles. Accordingly, there is no zero flux surface partitioning that exclusively constitutes the H … H bond path region. Hence the energy of this region is ill-defined.

Second, the case of dibenzene in Fig. 10, representative of the general complex of two molecules, R – H … H – R, between which H … H bond paths are found within some threshold intermolecular separation [6]. Here we do find a H … H bond bundle, allowing for a quantitative assessment of the H … H interaction.

The electron and kinetic energy counts of the H … H bond bundle are approximately $4\times 10^{-3}$ electrons and 2 kcal/mol (less than $0.1\leavevmode\nobreak\ eV$), which amounts respectively to a 0.20% and 0.26% share of each hydrogen’s 0.98 electrons and 380 kcal/mol of kinetic energy. The energies of these bond bundles are on the same order as Van der Waals interactions.

For such weak interactions, it is important to note the relationship between the analysis resolution, i.e. the number of differential gradient bundles used, and the minimum solid angle that a \acP-basin must occupy in order to be resolved. For example, all the calculations for this work used (approximately) twenty-thousand \acpdgb, which means that a \acP-basin whose solid angle is less than $1/20{,}000^{\rm{th}}$ of a steradian will not be recovered. It is possible that every bond path will map to a maximum in \acP if a sufficiently high resolution is used. Drawing from the H … H bond bundle of dibenzene (Fig. 10), we expect the energy and electron counts of similar tiny bond bundles to be negligibly small.

If all bond paths do map to maxima in \acP, one could speculate that \acpbcp may function as nucleation sites about which bond bundles grow, not dissimilar to the nucleation sites in crystal growth. Though this would change nothing about the instantaneous energetic significance of \acpbcp and bond paths associated with very small bond bundles, it would cast the \acbcp itself as a spacial indicator of regions likely to become energetically significant. The nature and extent of that significance could then be anticipated by the behavior of \acP in the regions that correspond to bond paths.

3.4 Missing bond points and paths

Equally disconcerting to those seeking to rationalize molecular properties based on a \acqtaim analysis are instances where \acpbcp are absent between atoms that, based on chemical reasoning, should be bound [32, 49, 50, 51, 52].

One such molecule is an iron trimethylenemethane complex represented in Fig. 11. Based on simple bond counting schemes, interatomic distances, source function studies, delocalization indices, and spectroscopic evidence, the Fe atom and the methylene carbon atoms (C2 in Fig. 11) should be bound [32]. Yet there is no \acbcp between these atoms, and hence no bond path linking their \acpncp. As shown in Fig. 11 and consistent with the preponderance of the evidence and with chemical intuition, from a bond bundle perspective there is indeed a bond between these atoms.

4 Condensed charge density and local kinetic energy

Though the local kinetic energy expressed in Eqs. 1 and 2 are those most discussed in the literature [43, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69], Anderson et al. [40] have shown that these are but two representations drawn from an apparently infinite number of definitions for the local kinetic energy and correctly assert that “Regardless of how one partitions the system, it seems that the kinetic energy of an atom in a molecule is not uniquely defined in quantum mechanics: for any choice of subsystem partitioning, one can always find two $\dots$ functions that give different values for the regional kinetic energy.” While the authors go on to comment on the conceptual utility afforded through the use of Eqs. 1 and 2, or linear combinations of the two—constituting the Laplacian class of kinetic energy operators—they conclude, “[We] cannot think of any physical or intuitive justification for excluding local kinetic energies from outside the Laplacian family.”

It is imperative to note that Anderson et al. were commenting on the quantum mechanical ambiguity associated with defining local kinetic energy and not reflecting directly on the value of approaches using one or another kinetic energy form. Still, the implicit, and in some cases explicit [70] expectation that QTAIM must be quantum mechanically rigorous has, in our opinion, compromised its utility.

Useful models often enlarge an existing conceptual framework, which for models of bonding means enhancing the chemist’s view of charge density expansion and polarization. These concepts are fundamental to many modern theories of bonding. For example, ab initio basis sets are specifically designed to capture the response of the charge density—and its underlying orbital basis—to expansion and polarization [71]. Basis sets are even named in a manner that allows the user to quickly evaluate their ability to recover these charge density responses, e.g. double zeta with polarization (DZP). The degree to which \acrho expands and polarizes upon bond formation is often used to categorize bonds as ionic, polar, covalent or polar covalent, and is frequently illustrated through the depiction of changes to charge density isosurfaces. It is the central role of charge density isosurfaces in mediating the structure of \acP that provides a rationale for representing the kinetic energy within the Laplacian family of energy operators.

Clearly, \acP derives it structure most directly from the gradient field. Consider an arbitrary gradient path \acgp and a point $p$ on this path. Through this point passes a charge density isosurface with its normal in the direction of the unit tangent vector to \acgp at $p$, which we represent as $\bm{\tau}(p)$, i.e.

[TABLE]

There are an infinite number of planes containing $\bm{\tau}(p)$. Each of these normal planes intersects the isosurface along a plane curve, and in general each is characterized by a different curvature (not to be confused with the curvature of \acrho). The principal curvatures at $p$, denoted $\kappa_{1}(p)$ and $\kappa_{2}(p)$, are the maximum and minimum values of these curvatures. The directions of these curves in the tangent plane at $p$ are referred to as principal directions, designated here as $\bm{e}_{1}(p)$ and $\bm{e}_{2}(p)$, with $\bm{e}_{1}(p)\cdot\bm{e}_{2}(p)=0$. An isosurface’s lines of curvature are surface geodesics that are by definition always tangent to a principal direction. There will be two orthogonal lines of curvature through every point on the surface that together act as a natural isosurface coordinate system.

Thus $\bm{\tau}$, $\bm{e}_{1}$ and $\bm{e}_{2}$, form an orthogonal moving frame ideally suited to describing charge density structure from a chemical perspective. The fields $\bm{e}_{1}$ and $\bm{e}_{2}$ recover the shape of isosurfaces and $\bm{\tau}$ provides information as to the nesting of these surfaces.

Using the $(\bm{\tau},\bm{e}_{1},\bm{e}_{2})$ moving frame, it is possible to derive the form of the differential area and volume elements introduced in Eq. 5. We begin by picking an arbitrary point $a_{0}$ on the arbitrary gradient path \acgp0 depicted in Fig. 12(a). (Here, consider the point to be located at roughly one third of the distance to the nearest atom.) This point lies on an isosurface ${\cal{S}}_{a}$. Around this point construct a planar square patch by placing one of its vertices at $a_{0}$ and the remaining vertices at $a_{1}=a_{0}+\epsilon\,\bm{e}_{1}(a_{0})$, $a_{2}=a_{0}+\epsilon\,\bm{e}_{2}(a_{0})$ and a_{3}=a_{0}+\epsilon\,\big{(}\bm{e}_{1}(a_{0})+\bm{e}_{2}(a_{0})\big{)}. This patch is normal to $\bm{\tau}(a_{0})$ and hence is a tangent plane to ${\cal{S}}_{a}$ at $a_{0}$. For sufficiently small $\epsilon$ the edges of this patch are principal curves. Passing through each vertex $a_{i}$ is a gradient path \acgpi.

Now pick another point on \acgp0, say $b_{0}$, which lies on the isosurface designated ${\cal{S}}_{b}$. A vector originating at $b_{0}$ and lying in the direction $\bm{e}_{1}(b_{0})$ will intersect \acgp1 at a point designated $b_{1}$ and in a similar fashion the vectors originating at $b_{0}$ and lying in the directions $\bm{e}_{2}(b_{0})$ and $\bm{e}_{1}(b_{0})+\bm{e}_{2}(b_{0})$ will intersect \acgp2 and \acgp3 respectively at points designated $b_{2}$ and $b_{3}$. The planar patch with vertices $b_{i}$ is necessarily a tangent plane to the isosurface ${\cal{S}}_{b}$ at $b_{0}$ with edges along principal directions.

At every point $s_{0}$ along \acgp there exists a coordinate patch (with vertices $s_{i}$) that is normal to $\bm{\tau}(s_{0})$ and spanned by the vectors $\bm{e}_{1}(s_{0})$ and $\bm{e}_{2}(s_{0})$ such that \acgpi is the union of all $s_{i}$.444This construction generates a set of triply orthogonal surfaces. As with any such set, they necessarily intersect along lines of curvature.

Through every point in any of these patches there is a gradient path that maps one-to-one and onto each of these coordinate patches. The differential gradient bundle may be taken as the union of all such coordinate patches.

The charge contained in a differential gradient bundle can be found by integrating over the set of volumes formed by coordinate patches separated by $ds$, yielding a volume element equal to $ds\times dl_{1}\times dl_{2}$ where $dl_{1}$ and $dl_{2}$ are the lengths of the edges of the coordinate patch.

Obviously, $dl_{1}$ and $dl_{2}$ are functions of $s$, varying along the extent of the differential gradient bundle. The form of this functionality may be discerned by considering a set of coordinate patches separated by a small distance $\Delta s$ (Fig. 12(b)). It is straightforward to show that for sufficiently small $\Delta s$,

[TABLE]

where $\kappa_{1}(s_{i})$ is the isosurface curvature in the $\bm{e}_{1}$ direction at $s_{i}$. It follows that,

[TABLE]

Noting that $l_{1}(s_{0})=d\theta\,s_{0}$ (see Fig. 1) we have,

[TABLE]

Similarly,

[TABLE]

and the area element, $dA(s)=l_{1}(s)\times l_{2}(s)$ is,

[TABLE]

where $H(s)$ is the mean curvature, \big{(}\kappa_{1}(s)+\kappa_{2}(s)\big{)}/2, of the isosurface at $s$.

Combining Eqs. 5 and 4, it is apparent that all gradient bundle condensed properties depend on the isosurface mean curvature along \acsgp. In particular, the gradient bundle condensed volume at a point $(\theta,\phi)$ of atom $i$ is,

[TABLE]

which is a characteristic of the charge density’s extrinsic structure and may be thought of as measuring the total mean curvature of the isosurfaces along a gradient path. The change in the magnitude of this quantity provides a direct measure of charge density expansion along a particular \acgp, while its derivative properties with respect to $\theta$ and $\phi$ measure polarization. As an illustration consider the evolution of \acV through the Lewis acid-base reaction discussed in Section 3.2.

Fig. 13 depicts the change in the charge density contours of the developing B – N bond bundle through the chemical reaction (recall Fig. 8). Before the bond bundle forms, the bond path intersects isosurfaces about the boron atom at points of both positive and negative mean curvature. In fact, the majority of boron isosurfaces are concave along the bond path. The fraction of concave isosurfaces decreases along with the B – N internuclear distance. At the equilibrium separation all boron isosurfaces are convex and a noticeable smooth corner has developed along the full length of the B – N bond path.

The behavior represented in Fig. 13 results from an increase in $p_{z}$ admixture from the boron atom. Initially, the $p_{z}$ atomic orbital is essentially empty, becoming partially occupied through overlap with the lone pair orbital of the nitrogen atom. Simply, the changing mean curvature along the bond path results from charge density expansion and polarization, and the total mean curvature as given by \acV and \acP serves as a measure of this response.

5 Metallic bonds

Real space models of metallic bonding date back to the earliest days of quantum mechanics. Pauling proposed that interatomic forces and metallic structure could be rationalized from a resonating-valence-bond perspective [72, 73]. Altmann et al. [74] employed directed valence bond approaches in an attempt to explain the preferred crystal structure of the non-magnetic transition metals. At nearly the same time, Engel and subsequently Brewer [75], based largely on correlations, suggested that the spherically averaged number of valence $s$-$p$–electrons was the determining factor favoring one metal structure over another, an assertion deemed deficient by none other than Hume-Rothery [76]. Much more recently, Hoffmann [23] unfolded band diagrams as a means of explaining the crystal structure preference of transition metals.

Setting aside the highly focused work of Hoffmann, none of the other investigations produced insights as useful as the chemical models from which the various researchers took inspiration. However, it is notable that in all cases the goal was to further chemical understanding of crystallographic structure. Historically and practically, the ability of chemical models to account for molecular and solid state structure has served as an acid test for their continued exploration. It is for this reason we compare the bond bundles of the BCC metal Nb with the FCC metal Cu as a way to demonstrate the potential of gradient bundle analysis. A more complete comparison of bond bundles across the transition metals will be the subject of an upcoming paper.

Fig. 14 depicts \acP and \acT for both Cu and Nb in FCC and BCC crystal structures. Though the stable structure of Nb is BCC and that of Cu is FCC, both structures of each metal were modeled as a way of comparing stable and unstable gradient bundle condensed properties.

As with the previous molecular calculations, there is a strong correlation between \acP and \acT. (A weaker correlation with \acV is not shown.) Remarkably, \acP appears to reveal the underling $d$-orbital character of \acrho, with Cu and Nb maxima capturing the nodal character of atomic orbitals that transform as the irreducible representations T2g and Eg respectively. Though the magnitude of condensed property features for both \acP and \acT change dramatically with crystal structure, their topologies do not. A fact that is readily explained by the observation regarding the orbital characteristics of gradient bundles. On the one hand, the topology of \acP is principally controlled by the distribution of electrons between Eg and T2g atomic orbitals—a symmetry property (OH for both BCC and FCC). On the other hand, the magnitude of \acT and hence \acP is controlled by the nature of the $\sigma$, $\pi$ and $\delta$ overlap between these atomic orbitals—a crystallographic property.

With the virial theorem satisfied at every point of the condensed charge density, local maxima in kinetic energy correspond to regions of low total energy and hence are stabilizing relative to neighboring regions. Accordingly, Fig. 14(b) reveals that the Eg-like gradient bundles stabilize the BCC structure of Nb, while T2g-like gradient bundles stabilize FCC Cu.

One “lobe” of the high kinetic energy T2g-like gradient bundles most responsible for stabilizing the FCC structure of Cu is represented in Fig. 15 alongside the gradient bundles emanating from the same region of FCC Nb. As expected, these gradient bundles are similar in shape. However, for Nb, with five valence electrons, they emerge from local minima in \acP and hence are “under occupied” compared to neighboring gradient bundles. For Cu, with eleven valence electrons, they emerge from local maxima in \acP, indicating that they are preferentially occupied by the valence electrons added when traversing the transition metal series from group 5 to group 11. Of course, such preferential occupation is driven by a lower total energy—higher kinetic energy—with respect to all other gradient bundle occupations for the same number of electrons. All in all, condensed charge density accumulates in the gradient bundles that will maximally lower the system total energy. The correlation between \acP and \acT merely confirms the well-defined energy of gradient bundles.

In an analogous fashion, the high kinetic energy, Eg-like, gradient bundles most responsible for stabilizing BCC Nb are represented in Fig. 16 alongside those emanating from the same region of BCC Cu. Again these gradient bundles are similarly shaped though obviously they are preferentially filled by the first valence electrons of the transition metal series.

We speculate that the preferential filling of gradient bundles is due to their relative content of bonding and antibonding character. To illustrate, Fig. 17 shows early and late filling gradient bundles superimposed on representative bonding and antibonding molecular orbital contour diagrams.555The orbital contours are from large cluster calculations simulating the FCC and BCC structures. It has been demonstrated that the central atom of such large clusters possesses an energy within a fraction of an eV of that from a bulk calculation and that the resulting charge densities are indistinguishable from those resulting from bulk calculations [77].

The variation in the contour curvature is controlled by orbital nodal character. The contour lines near nodes have a small mean curvature and experience their maximum mean curvature along gradient paths that are at greatest distance from nodes.

As a consequence, in the case of bonding orbitals it is the orbital gradient paths (distinct from \acrho gradient paths) more or less aligned with orbital antinodes that intersect the interatomic surface along orbital contours of maximum curvature (Fig. 17(a)). For antibonding orbitals the situation is quite different. The orbital contours near the interatomic surface possess low mean curvature, as do those near the orbital’s intrinsic angular nodes. Such constraints force the formation of a curvature corner along the orbital gradient path roughly bisecting the interatomic and angular nodes (Fig. 17(b)).

The gradient field of \acrho must reflect the character of its orbital basis. Gradient bundles containing predominately bonding character will approach the interatomic surface along nearly normal directions, while gradient bundles containing antibonding character will approach the surface with a tangential component.

As illustrated in Fig. 17, the gradient bundles stabilizing BCC Nb intersect regions that are predominately bonding and actually increase their volumes in these regions by intersecting orbital contours along paths near high mean curvature maxima. The gradient bundles stabilizing FCC Cu maximize their volume in antibonding regions. Leading to the seemingly paradoxical argument that FCC Cu is stabilized by antibonding interactions.

The paradox is easily resolved by noting that the virial theorem requires there be a corresponding decrease in the potential energy within the stabilizing gradient bundle. This decrease results from a radial contraction of the charge density in the part of the gradient bundle closer to the nucleus, which screens that non-radial—curving—part of the bundle and allows it to expand into the high kinetic energy antibonding region. Basically the gradient bundle of Fig. 17(b) is harvesting excess kinetic energy from the antibonding region to stabilize the FCC structure.

The previous conjecture is supported by considering the shape of the FCC stabilizing gradient bundle when Cu is forced BCC, as shown in Fig. 18. Just as for the FCC structure, the gradient bundle projects into the crystallographic “octahedral” hole. However, this interstice is less tightly packed when BCC, and in fact hosts a \acbcp and bond path between second neighbors. Obviously the larger expanse of the BCC versus FCC octahedral holes diminishes the intensity of antibonding interactions in the former, and accordingly the amount of kinetic energy that can be harvested from this hole. The effect on the gradient bundle is evident, it does not curve and expand into the octahedral hole as it does in the FCC structure.

We conclude that the stability of the FCC transition metals results from maximizing antibonding interactions, which is why the structure occurs late in the series where antibonding orbitals of the series are being filled. The BCC structure is stabilized by gradient bundles that contain predominantly bonding character and therefore occurs early in the series.666 The reader may wonder why the HCP structure occurs early (when predominately bonding orbitals are filled) and also just after the midpoint of the transition metal series when antibonding orbitals are predominately filled. Though these two groups of transition metals have the same crystal structure, their charge density topologies are different and hence they possess different gradient fields and gradient bundle structures. [78]

The consequence of each structure is that its respective stabilizing gradient bundles contain maximal kinetic energy, thus maximizing the redistribution of charge and energy, yielding the minimum system total energy.

6 Summary

We have introduced a projected space that appears to be ideally suited to the investigation and discovery of charge density-property relationships in all classes of molecules and materials. This space, denoted \acP, constitutes all volumes bounded by zero flux surfaces in the gradient of the charge density. As a consequence, when confined to the conventional Laplacian family of kinetic energy operators, each point of \acP is endowed with a well-defined energy. As a subset of the chemically meaningful structures of \acP, all of those intrinsic to the \aclqtaim are recovered, e.g. Bader atoms. In addition, we also recover a class of nonstandard \acqtaim structures that maximize charge density concentration and minimize energy, properties commonly associated with chemical bonds. These structures are called bond bundles. We show that bond bundles are recoverable from \acrho for any class of atomic system, provide more chemical information, and are not subject to the ambiguities and misgivings that have been associated with the \acqtaim bond path. We argue that the robust properties of bond bundles are rooted in the structure of \acP, which is endowed with a natural reference frame that is determined by the system’s charge density isosurfaces. We demonstrate this fact by successfully analyzing transition metal structure within the bond bundle construct.

This work should serve to further attempts at approaching traditional metallurgical problems from a chemical perspective. And also provides a new and possibly fruitful approach to analyze charge density in terms of its natural reference frame and the variations in this frame associated with chemical reactions and mechanical perturbations.

Computational methods

All simulations were performed with the \acadf [79, 80, 81] ab initio software using the Perdew-Burke-Ernzerhof (PBE) functional [82] and a triple-zeta with polarization (TZP) all-electron basis set. Calculation of \acP was performed within the Tecplot 360 visualization package [83] using the \aclgba tool of the in house Bondalyzer package by the \acmtg at \accsm.

Acknowledgments

Support of this work under Office of Naval Research Grant No. N00014-10-1-0838 is gratefully acknowledged.