Electron Delocalization Range in Atoms and on Molecular Surfaces

Jun 10, 2016 - The electron delocalization range function EDR(r⃗; d) (J. Chem. Phys. 2014, 141, 144104) quantifies the extent to which an electron a...
1 downloads 8 Views 5MB Size
Article pubs.acs.org/JCTC

Electron Delocalization Range in Atoms and on Molecular Surfaces Benjamin G. Janesko,*,† Kenneth B. Wiberg,‡ Giovanni Scalmani,¶ and Michael J. Frisch¶ †

Texas Christian University, Fort Worth, Texas 76129, United States Yale University, New Haven, Connecticut 06520, United States ¶ Gaussian, Inc., 340 Quinnipiac St., Bldg. 40, Wallingford, Connecticut 06492, United States ‡

S Supporting Information *

ABSTRACT: The electron delocalization range function EDR(r;⃗ d) (J. Chem. Phys. 2014, 141, 144104) quantifies the extent to which an electron at point r ⃗ in a calculated wave function delocalizes over distance d. This work illustrates how atomic averages of the EDR, and plots of the EDR on molecule surfaces, provide a chemically intuitive picture of the sizes of occupied orbital lobes in different regions. We show how the surface and atomic delocalization distinguish aminophosphine ligand’s hard N and soft P lone pairs, distinguish the site preference for Ag+ cation binding to conjugated oligomers, and provide information that is different from and complementary to conjugation lengths. Applications to strong correlation and the prototropic tautomerization of phosphinylidenes R1R2HPO illustrates how the surface and atomic delocalization can work with other tools to provide a nuanced picture of reactivity.

1. INTRODUCTION Electron delocalization is widely used to rationalize ideas of bonding, conjugation, and aromaticity.1 Computed molecular orbitals are often used to interpret delocalization. However, orbital-based interpretations have some limitations.2 Orbital localization can change the shape of individual orbitals without changing computed observables.3−5 Single molecular orbital configurations (single Slater determinant wave functions) provide an incomplete picture of strongly correlated systems, such as dissociating covalent bonds.6 One may instead focus on the one-particle density matrix γ(r,⃗ r′⃗ ) of N-electron wave function Ψ, γ( r ⃗ , r ′⃗ ) ≡ N

interactions give negative contributions. Off-diagonal elements of the Hilbert-space density matrix give rise to quantities such as Hückel,11,12 Wiberg,13 and Mayer14 bond orders. More broadly, the off-diagonal two-particle reduced density matrix is connected to superconductivity,15 while the off-diagonal elements of transition density matrices quantify excited-state coherence.16−18 It is difficult to interpret a function of six variables such as γ(r,⃗ r′⃗ ). Many investigators have worked to extract chemically relevant information out of the density matrix. As discussed above, the Hilbert-space off-diagonal density matrix quantifies bond orders.11−13 The real-space density matrix γ(x, x′) plotted for points x along a symmetry axis or after spherical averaging provides insight into atoms and diatomics.19,20 Integrating the one-particle density matrix, density matrix square root, or Fermi hole over atomic21 or other22 domains provide orbital- and basis-set-independent bond orders.10,23−29 (These contributions to bond orders are arguably more important than contributions from the irreducible two-particle density matrix.30) Domain-averaged Fermi holes31,32 and plots of γ(r,⃗ rR⃗ ) or the Fermi hole about reference point rR⃗ 10,33 give finer details of electron−electron interactions. Deviations from idempotency 1 − ρ−1(r)⃗ ∫ d3 r′⃗ |γ(r,⃗ r′⃗ )|2 measure the spatial distribution of effectively unpaired electrons.34,35 Recent extensions use the B05 exchange hole model36−38 to identify effectively unpaired electrons within idempotent density matrices.39 Partitioning the one-particle density matrix into

∫ d3 r2⃗ , ..., d3 rN⃗ Ψ( r ⃗ , r2⃗ , ..., rN⃗ )Ψ*( r ′⃗ , r2⃗

, ..., rN⃗ )

(1)

(spin dependence is suppressed for conciseness). The density matrix is invariant to unitary transforms of integer-occupied orbitals and robust to noninteger orbital occupancy. Diagonal elements ρ(r)⃗ = γ(r,⃗ r)⃗ give the probability density for finding an electron at r.⃗ Off-diagonal elements, like molecular orbitals themselves, are not quantum-mechanical observables. Offdiagonal elements contribute to the electronic kinetic energy, whose expectation value is an explicit functional of the oneparticle density matrix7 but not of the electron density.8,9 Offdiagonal elements contribute to the pair probability of finding like-spin electrons at points r ⃗ and r′⃗ via Fermi hole terms ∝ |γ(r,⃗ r′⃗ )|2.10 Off-diagonal elements are connected to chemical bonding: covalent bonds between atoms A and B typically give positive contributions to γ(r ⃗ ∈ A, r′⃗ ∈ B), while antibonding © 2016 American Chemical Society

Received: April 5, 2016 Published: June 10, 2016 3185

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation

The first new interpretive tool is the surface delocalization, defined as D(r)⃗ at points r ⃗ on a molecular surface. This complements the common practice of plotting electrostatic potentials on molecular surfaces.65 Among the many definitions of “molecule surface”, we use the (spin)density isosurface ρ(r)⃗ = 0.001 electrons/bohr3.66 D(r)⃗ is sensitive to the choice of isosurface, with ref 62. proving that a uniform electron gas with constant spin density ρ has D ∝ ρ−1/3. However, we will show that variations in D(r)⃗ on a particular isosurface provide chemically useful information. (Figure SI-5 in the Supporting Information illustrates that trends in the benzene D(r)⃗ values are robust when plotted on different density isosurfaces.) The second new interpretive tool is the atom delocalization DA, defined as the average delocalization length of electrons assigned to atom A:

pairs of atomic domains γAB(r,⃗ r′⃗ ) provides additional insights into covalent bonding and electron sharing.40−42 The kinetic energy density8,43 provides the electron localization function (ELF) and several related descriptors.44−46 (Note that the ELF can be expressed in terms of the full two-particle reduced density matrix,47,48 and has been interpreted in terms of the covariance of the electron pair distribution46,49,50 and the nodal structure of occupied orbitals.51) Partitioning the kinetic energy into “classical” and “exchange” contributions gives an alternative perspective on bonding.52,53 The localized orbital locator (LOL) and parity function provide complementary perspectives on the off-diagonal structure about rR⃗ .54−56 Some of us recently proposed an electron delocalization range function57 EDR(r;⃗ d), which complements these other tools by focusing on the distance between points r ⃗ and r′⃗ in γ(r,⃗ r′⃗ ). The EDR projects the density matrix onto a reference density matrix spanning distance d, to quantify the degree to which electrons at point r ⃗ delocalize over distance d: EDR( r ⃗ ; d) =

∫ d r ′⃗ g( r ⃗ , r ′⃗ ; d)γ( r ′⃗ , r ⃗)

(2)

g ( r ⃗ , r ′⃗ ; d) ≡

⎛ | r ⃗ − r ′⃗ |2 ⎞ ⎛ 2 ⎞3/4 1 ⎜ ⎟ exp⎜ − ⎟ 2 ⎝ πd ⎠ ⎝ d2 ⎠ ρ( r ⃗)

(3)

DA =

g(r,⃗ r′⃗ ; d) is constructed from a Gaussian fit to the one-particle density matrix of a uniform electron gas (UEG) with density ρ ∝ d−3.58 This reference density matrix ensures |EDR(r;⃗ d)|2 ≤ 1. Equation 2 is connected to the Fermi hole normalization34 and to the B05 model’s fit of the exchange hole to a localized model,36−39 and is reminiscent of quantities appearing in the electron pair localization function.59 The EDR is based on earlier work projecting the density matrix onto model density matrices found in (semi)local model exchange holes.60 We have previously explored four methods to visualize and quantify the EDR: system averages ⟨EDR(d)⟩ = ∫ d3 rρ⃗ (r)⃗ EDR(r;⃗ d), plots of EDR(x; d) for points x along symmetry axes,57 contour or isosurface plots of EDR(r;⃗ dspecial) at special length scales dspecial,61,62 and topological analyses.63 We find that the delocalization length dmax maximizing the change in ⟨EDR(d)⟩ upon ionization characterizes weakly bound electrons in solvent clusters and electrides,61,62 extending existing orbital-based tools64 and allowing facile studies of strong (nondynamical) correlation among multiple weakly bound electrons. Local maxima in the EDR distinguish, e.g., core vs valence regions and normal vs stretched bonds.63 The present contribution introduces two new interpretive tools built from the EDR. Both are based on the “characteristic delocalization length” of electrons at point r,⃗ d

(5)

Among the many methods21 for determining weights wA(r)⃗ assigning point r ⃗ to atom A, we use the Hirshfeld partitioning as implemented in Gaussian 09.67 This work presents numerical illustrations of surface and atomic delocalization. These orbital-independent tools prove to characterize the size of orbital lobes, providing a convenient extension of ideas from molecular orbital theory. Applications to amphiphilic ligands, conjugated systems, strong correlation, and a chemical reaction highlight these new tools’ utility for building chemically intuitive interpretations.

3

D( r ⃗) = arg max EDR( r ⃗ ; d)

∫ d3 r ρ⃗ ( r ⃗)D( r ⃗)wA( r ⃗)

2. COMPUTATIONAL METHODS Calculations use the EDR as implemented in the development version of the Gaussian suite of programs.57,58,68 Like our previous studies,57,61,62 this work uses one-particle density matrices calculated with Hartree−Fock theory, Kohn−Sham density functional theory (DFT),69 generalized Kohn−Sham DFT,70 coupled cluster singles and doubles (CCSD) theory, and complete active space self-consistent field (CASSCF) approaches. DFT calculations use the local density approximation71 (LDA), B3LYP,72−75 or ωB97X-D76 approximate exchange-correlation functionals. CCSD density matrices are evaluated using the Z-vector method.77,78 Calculations use standard atom-centered Gaussian basis sets.79−81 Calculations treat the nuclei as classical point charges with zero kinetic energy. Calculations on closed-shell systems show the density matrix of only the alpha spins. Calculations on one-electron systems show the alpha-spin density matrix. Atomic averages of the valence delocalization length or kinetic energy are evaluated by first removing the core orbitals from the one-particle density matrix, then using this modified density matrix to evaluate the electron density, kinetic energy density, and EDR. Figure SI-1 in the Supporting Information illustrates this approach for the Be atom. This work follows the perspective of refs 70, 82, and 83. The one-particle density matrix of a Kohn−Sham reference system is not treated as an approximation to the interacting system’s density matrix, but as a rigorously defined quantity obeying its own set of exact conditions. As discussed above, D(r)⃗ is generally plotted on the 0.001 electrons/bohr3 isosurface of the majority-spin electron density. Individual molecular orbitals ψ(r)⃗ are plotted on the corresponding |ψ ( r ⃗)| = 0.001 bohr−3/2 isosurface, and colored by the orbital sign. The D(r)⃗ color scale is such that delocalization increases in the direction red → orange → yellow

(4)

Here, “argmaxd” denotes the value of the argument variable d that maximizes the function of interest. The characteristic delocalization length D(r)⃗ is thus defined as the distance d, which maximizes EDR(r;⃗ d) at point r.⃗ Operationally, at each point r,⃗ we evaluate EDR(r;⃗ di) on a grid of ∼30 different di values, determine the value of di that gives the maximum EDR, then perform a three-point numerical fit about that value. (Section SI-II in the Supporting Information highlights how EDR(r;⃗ d), in special circumstances, may have multiple local maxima in d. For simplicity, the present work considers only the global maximum at each r.⃗ ) 3186

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation

Figure 1. H atom (left), ground-state H+2 (middle), and first excited-state H+2 (right). The top row represents molecular orbital |ψ ( r ⃗)| = 0.001 bohr−3/2. The bottom row represents density isosurfaces ρ(r⃗) = 0.001 electrons/bohr−3 painted with delocalization lengths D(r⃗) from 1.8 bohr (red) to 2.5 bohr (blue).

→ green → blue. Many pictures of molecular density isosurfaces show only the back half of the isosurface. This convention allows us to display the molecular geometry and density isosurface in a single figure. We follow the conventional practice of representing molecular geometries with spheres at the atomic positions, connected by single, double, or broken lines. While these lines generally correspond to chemically reasonable bond orders, they are not evaluated from quantummechanical calculations, and are included merely as a guide to the eye.

orbital lobes. The H atom has the largest and most diffuse orbital, with the lowest electronic kinetic energy, the largest atomic delocalization DH, and the largest (dark blue) surface delocalization. Ground-state H+2 is modestly more localized, with a larger electronic kinetic energy consistent with the virial theorem,84 a smaller DH, and a smaller (light blue) surface delocalization. Excited-state H+2 has an internal node giving two relatively small orbital lobes. This increases the electronic kinetic energy and reduces the atom and surface delocalization lengths. Figure SI-2 in the Supporting Information shows the positive-semidefinite electronic kinetic energy density,8,43

3. RESULTS 3.1. Model Systems. We begin by illustrating the surface and atom delocalization of three one-electron systems, the ground-state H atom, and the ground and first excited states of H+2 at a bond length of 1.057 Å. Calculations use the aug-ccpVQZ basis set. Figure 1 shows the molecular orbitals and surface delocalization. (For illustration purposes, the excitedstate orbital is taken as the UHF ground-state lowest unoccupied molecular orbital (LUMO).) Table 1 presents

τ ( r ⃗) =

on these density isosurfaces, plotted as the localized orbital locator (LOL),55,56 LOL =

τ0 =

H+2

kinetic energy, KE (Hartree) atomic delocalization length, DH (bohr)

H atom

ground state, GS

excited state, ES

0.50 2.6

0.60 2.3

0.80 1.7

τ0/τ 1 + (τ0/τ )

where τ0 is the kinetic energy density of a uniform electron gas with electron spin density ρ(r)⃗ :

Table 1. Computed Kinetic Energies and Atomic Delocalization Lengths for the Systems Depicted in Figure 1

property

1 [∇r ·∇r γ( r ⃗ , r ′⃗ )]r ⃗ = r ⃗ ′ 2 ⃗ ⃗′

3 (6π 2)2/3 ρ5/3 ( r ⃗) 10

(LOL variations on a density isosurface where τ0 is constant highlight variations in τ.) While the LOL is relatively large on H atom and small on ground-state H+2 , it is rather large (blue) on the weakly curved periphery of the excited-state orbital. This “short-sighted” behavior of τ has been discussed previously.56 We suggest that the surface delocalization introduced in this work measures the scale of individual orbital lobes, a scale that has a tendency to be intermediate between the orbital curvature measured by the kinetic energy density and the systemaveraged off-diagonal delocalization measured by the parity function.56 Table 2 provides a more detailed analysis of how surface and atomic delocalization change with orbital lobe size. The table shows the surface and atomic delocalization of spherically symmetric, filled s, p, d, f, and higher shells. Calculations all use

electronic kinetic energies and atom-averaged delocalization lengths DH. For reference, ref 20 plots the one-particle density matrix of ground-state H2 along the bond axis, ref 61 presents a direct comparison of the density matrix and EDR of onedimensional systems, and ref 63 shows local maxima of the EDR in H2. The results show that the surface and atom delocalization lengths are qualitatively consistent with the size of individual 3187

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation Table 2. Results for Filled Atomic Shellsa shell

KE/e (Hartree)

Ex/e (Hartree)

DA (bohr)

rsurf (bohr)

102 LOL(rsurf)

D(rsurf) (bohr)

s p d f g h i

1.5 2.5 3.5 4.5 5.5 6.5 7.5

−0.564 −0.517 −0.458 −0.416 −0.385 −0.361 −0.341

1.27 0.86 0.70 0.60 0.54 0.50 0.46

1.77 2.14 2.39 2.59 2.76 2.92 3.06

0.73 0.61 0.55 0.51 0.47 0.44 0.42

2.10 1.50 1.24 1.08 0.97 0.88 0.82

characteristic delocalization lengths of hydrated electrons (H2O)−N match the HOMO radius of gyration.61,64 These trends in valence delocalization lengths also match KE−1/2, consistent with Table 2. The HOMO radius for helium is relatively small, compared to DA and KE−1/2, while the HOMO radii for d-block atoms are relatively large, compared to DA and KE −1/2 . We speculate this difference arises from the delocalization length’s strong dependence on angular nodes seen in Table 2. The valence shells of heavy elements have a tendency to have relatively long atomic delocalization lengths. For example, the ratio of DA to KE−1/2 in Figure 2 is 1 for Na and K atoms. This appears to be because the atomic delocalization is relatively sensitive to lowdensity regions of delocalized electrons, like those in the outer lobe of potassium’s valence 3s orbital. Section SI-II analyzes this effect in a model system. Table SI-1 illustrates the basis set dependence of the surface and atomic delocalization of the Li atom. Both are reasonably well-converged in double-ζ basis sets, broadly consistent with Hirshfeld populations86 and ELF basin populations (however, see ref 87). 3.3. Small Molecules. We next illustrate the surface delocalization of a representative small molecule, the aminophosphine ligand Me2N−CH2−CH2−PMe2. Aminophosphines combine a chemically hard88 nitrogen lone pair with a softer phosphorus lone pair. They are widely used in catalytic asymmetric hydrogenation,89 and their amphiphilic hard−soft character is important to their coordination chemistry.90 While tools such as computed molecular electrostatic potentials capture the local acid−base character of such molecules, measuring this ”local hardness” remains a challenge.91,92 Figure 3 shows the aminophosphine surface delocalization evaluated from the CCSD/6-311++G(2d,2p) density matrix at

a

Kinetic energy per electron (KE/e), exchange energy per electron ((Ex/e)), atom-averaged delocalization DA, radius rsurf of the 0.001 electrons/bohr−3 density isosurface in these spherically symmetric systems, and LOL(rsurf) and surface delocalization D(rsurf).

the same radial functions, a single Gaussian function with an exponent of 1 bohr−2, such that the only differences between the results are the number of internal angular nodes in the orbitals. The kinetic energy per electron smoothly increases with the number of angular nodes. Atom and surface delocalization lengths decrease with increasing angular nodes, scaling approximately as (KE/e)−1/2. The exchange energy per electron varies more slowly with the number of angular nodes, because of its dependence on the positive-semidefinite absolute value squared of the density matrix. The LOL on the molecule surface also varies relatively slowly with the number of angular nodes, which is consistent with Figure SI-2. 3.2. Isolated Atoms. We next illustrate periodic trends in atom-averaged delocalization lengths. Figure 2 shows results

Figure 2. Valence atomic delocalization lengths DA and kinetic energies KE−1/2 of isolated atoms H−Kr, and HOMO radii from ref 85. Delocalization lengths and HOMO radii given in bohr, kinetic energies given in units of 3(Hartree−1/2).

from Hartree−Fock calculations in the large all-electron UGBS basis set. Calculations sum majority- and minority-spin contributions, and include only contributions from the valence orbitals. Atomic cores are [He] for Li−Ne, [Ne] for Na−Ar, [Ar] for K−Zn, and [Ar]3d10 for Ga−Kr. Figure SI-3 in the Supporting Information shows all-electron results. The figure also includes the inverse square root of the atomic kinetic energy (KE−1/2) and the atomic HOMO radius (from ref 85). The results in Figure 2 follow the expected periodic trends. Valence delocalization decrease along rows of the periodic table, and increase down columns. p-block elements are generally more localized than s-block elements, and d-block elements are quite localized, consistent with Table 2. These trends in valence delocalization lengths match the atomic HOMO radii,85 just as our previous study showed that the

Figure 3. Graphic depictions of aminophosphine: (left) computed geometry and (right) surface delocalization D(r⃗) from 2.8 bohr (red) to 3.5 bohr (blue).

the B3LYP/6-311++G(2d,2p) equilibrium geometry. Electrons on the molecular surface are relatively localized (red) in the nitrogen lone pair region, relatively delocalized (blue) in the phosphorus lone pair region, and modestly delocalized (yellow) in the alkyl regions. The surface delocalization distinguishes the delocalized, soft Lewis base site from the localized, hard Lewis base site. Figure SI-4 in the Supporting Information reports the B3LYP/6-31+G(d,p) surface delocalization of other representative small molecules. The results largely mirror Figure 3: 3188

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation

smaller than that of Li+. Despite this caveat, we suggest that, based on Figure 3, the surface delocalization can help distinguish the relative hardness of different sites. 3.4. Orbital Lobe Size Versus Conjugation Length. Before continuing, we present a result that requires careful explanation. The word ”delocalization” can sometimes mean different things in different contexts. Such distinctions are common in the literature on excited states, where transition density matrices are characterized in terms of diagonal and offdiagonal delocalization.16−18 Here, we emphasize (cf. Figure 1) that the ”delocalization” measured by the EDR corresponds to the size of individual orbital lobes, not to conjugation lengths. Figure 4 illustrates this with the HOMO and surface delocalization of ethylene and linear π-conjugated 1,3,5,7,9pentadecene. (These calculations use B3LYP/6-31+G(d,p) density matrices. Figure 5 shows that correlated wave functions

CC bonds, and the lone pairs of O and N atoms, are more localized than normal C−C or C−H bonds. An important caveat is that surface delocalization does not capture all aspects of chemical hardness. Chemical hardness is formally defined as the second derivative of energy, with respect to electron number, and can be approximated from the HOMO−LUMO gap.88,93−95 Hardness is also related to system size.85,96 The surface delocalization is dependent on the occupied orbitals, but not on the HOMO−LUMO gap or the virtual orbitals; therefore, it cannot capture all aspects of chemical hardness. Table 3 reports the hardness parameter93 Table 3. Spherical Atomic Lewis Acids’ Valence Shell Angular Momentum, Hardness Parameter (ηA), Radius (rsurf) of the 0.001 electrons/bohr−3 Density Isosurface, and Surface Delocalization Lengths (D(rsurf))

a

acid

valence

ηAa (eV)

rsurf (bohr)

D(rsurf) (bohr)

Li+ Na+ K+ Rb+ Cs+ Cu+ Ag+ Au+

s p p p p d d d

35.1 21.1

1.68 2.31 3.07 3.37 3.76 3.00 3.22 3.45

2.04 1.83 2.37 2.59 2.88 1.61 1.75 1.85

11.7 6.3 6.9 5.7

Data taken from ref 93.

and surface delocalization of some spherically symmetric Lewisacidic cations. Surface delocalization is computed from Hartree−Fock calculations in the def2-TZVP basis set. For atoms with a given type of valence shell, surface delocalization is approximately inversely correlated with hardness: Na+ is harder and more localized than Rb+, and Ag+ is harder and more localized than Au+. However, trends across different valence shells are not consistent: Li+ is harder and more delocalized than Na+, while Rb+ is harder and more delocalized than the Group 11 cations. This is readily rationalized by Table 2: the Na+ valence 2p orbital lobes are simply smaller than the Li+ valence 1s orbital lobe, thus the Na+ surface delocalization is

Figure 5. Graphic depiction of benzene surface delocalization. CCSD/ 6-311++G(2d,2p) density matrices, D(r)⃗ from 2.8 bohr (red) to 3.1 bohr (blue).

give similar results.) Each of the 10 lobes of the pentadecene HOMO are somewhat smaller than the lobes of the ethylene HOMO in this visualization. The surface delocalization

Figure 4. Graphic depictions of ethylene (left) and 1,3,5,7,9-pentadecene (right): (top) HOMO, viewed at tilted orientation, and (bottom) surface delocalization D(r⃗) from D = 2.8 bohr (red) to D = 3.2 bohr (blue). 3189

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation captures this, with pentadecene’s central CC bonds being more localized (more red) than the ethylene π bond. This does not imply that pentadecene is “not conjugated”, merely that the individual lobes of pentadecene’s occupied orbitals had a tendency to be rather small. Figure 5 shows the surface delocalization of benzene evaluated with CCSD/6-311++G(2d,2p) density matrices and B3LYP/6-31+G(d,p) geometries. Consistent with Figure 4, the relatively small orbital lobes in the π electron region are highlighted with relatively small (red) delocalization lengths. Again, this does not imply that benzene is ”not conjugated” or ”not aromatic”, merely that the individual lobes of benzene’s conjugated π orbitals are rather small. Figures SI-5 and SI-6 in the Supporting Information show that qualitatively similar results are found on other density isosurfaces, and in calculations with other basis sets or Hartree−Fock or DFT wave functions. Figure SI-7 in the Supporting Information shows the EDR at other length scales EDR(x;d) for points x along benzene’s symmetry axes, confirming the presence of a single maximum in d at these points. Note that regions between adjacent C−H bonds in Figures 4 and 5 are dark blue, indicating relatively large delocalization lengths. This is arguably an artifact of the test function g(r,⃗ r′⃗ ; d) in eq 3, arising when points r ⃗ and r′⃗ are on adjacent nonbonded atoms. Figure 1 of ref 61 discusses this effect in detail. We next argue that the orbital lobe size measured by the EDR is chemically relevant, and is both distinct f rom and complementary to measures of conjugation length. (This is analogous to the distinct and complementary diagonal and offdiagonal delocalization of transition density matrices mentioned above.16−18) Of course, many properties of conjugated oligomers can be explained in terms of conjugation length; including polarizabilities, conductivities, and UV/visible absorbance spectra.97 However, some other properties are relatively insensitive to conjugation length. For example, the computed Li+ affinities of butadiene, hexatriene, and octatetraene are −26.4, −29.2, and −31.5 kcal/mol; showing relatively weak dependence on conjugation length.98 (Organic chemists will no doubt agree that the terminal H2CCH− groups of these molecules are, in many ways, chemically similar, which is a transferability or nearsightedness often discussed in the literature.99,100) We find that the relatively small surface delocalization of pentadencene’s central CC bonds, highlighted in Figure 4, is consistent with their computed hard−soft base behavior. Pentadecene’s binding affinities for the soft Ag+ cation are relatively strong (−2.10 eV) at the terminal CC bond, and relatively weak (−1.92 eV) at the central CC bond. In contrast, the hard Li+ cation binds to pentadecene’s C−C bonds (as in ref 98) with little preference for central vs terminal bonds (−1.50 eV vs −1.53 eV). (These binding affinities are computed with the dispersion-corrected ωB97X-D exchange-correlation functional, the LANL2DZ basis set and effective core potential, and no counterpoise corrections. Figure SI-8 in the Supporting Information shows the computed adsorption geometries and total energies.) Table 4 shows ωB97X-D/LANL2DZ computed binding energies ΔE for Ag+ and Li+ binding to the terminal H2CCH− group of all-trans polyacetylene oligomers. Individual binding energies increase with conjugation length, consistent with the charge-transfer arguments of ref 98. However, the relative affinity for soft Ag+ is almost constant, which is consistent with the correspondingly small variations in the terminal group’s atomic delocalization.

Table 4. Results for All-trans Polyacetylene Oligomers H−(HCCH−)NH, Including the HOMO−LUMO Gap, Total (Valence + Core) Atomic Delocalization DC and DH of the Terminal CH2 Group, Binding Energies ΔELi+ and ΔEAg+ for Li+ and Ag+ Binding to the Chain End, and the Ratio ΔEAg+/ΔELi+ Giving the Relative Affinity for Ag+ Total (Valence + Core) Atomic Delocalization (bohr)

Binding Energy (eV)

N

gap (eV)

DC

DH

ΔEAg+

ΔELi+

ΔEAg+/ΔELi+ ratio

2 3 4 5 6

9.5 8.2 7.4 6.9 6.5

1.340 1.340 1.340 1.340 1.340

1.962 1.963 1.963 1.963 1.963

−1.66 −1.85 −1.99 −2.10 −2.18

−1.20 −1.34 −1.45 −1.53 −1.59

1.38 1.39 1.38 1.38 1.38

We conclude this section with two comments. First, we reiterate that we use “orbital lobe” language to aid understanding, not because the EDR is, in some way, “orbitaldependent”. Like other density-matrix-based descriptors, the EDR is independent of unitary transforms among integeroccupied orbitals. Localizing pentadecene’s occupied orbitals will not change the surface delocalization in Figure 4, just as it will not change the electron density, density Laplacian, ELF, or LOL. Moreover, Figure 5, as well as Figure SI-6 in the Supporting Information, confirm that different levels of theory give comparable surface delocalization. Second, we note that the distinction between conjugation length and orbital lobe size is analogous to our previous study of the EDR in the uniform electron gas.62 D(r)⃗ in the UEG, evaluated from the exact Kohn−Sham density matrix, is a constant proportional to ρ−1/3, giving delocalization lengths of ∼1−10 bohr for typical valence densities. We argued previously that this is chemically reasonable, in that metallic bonding does not imply large bond orders between distant metal atoms. Here, we argue that Figure 4 is also chemically reasonable, in that pentadecene’s large conjugation length does not imply large covalent bond orders between, e.g., the α and ω C atoms. 3.5. Strong Correlation. A strength of density-matrixbased tools such as the EDR is their ability to treat mean-field and correlated calculations on an equal theoretical footing. We next illustrate how our tools quantify the effects of ”strong” nondynamical correlation. We follow recent work by Feixas and co-workers,48 who used carbon monoxide CO stretched to a bond length of 2 Å to highlight correlation effects on the ELF. Our topological analysis of the EDR also considered this system.63 Figure 6 shows the surface delocalization evaluated for equilibrium and stretched CO. Table 5 shows the corresponding atomic Hirshfeld charges QC = −QO, atomic delocalization lengths DC and DO, and atom-averaged kinetic energies per electron KEC and KEO. Results are presented for Hartree−Fock theory, DFT, CCSD, and multireference CASSCF(10,10) calculations in the aug-cc-pVTZ basis set. Figure SI-9 in the Supporting Information confirms that the surface delocalization at the equilibrium bond length is less sensitive to the method used. The dipole moment of equilibrium CO is computed to be +0.24 D for Hartree− Fock theory and −0.11 D for CASSCF, which is consistent with previous work.101 At the equilibrium bond length, the surface delocalization distinguishes the more localized electrons on the oxygen atom (orange). Hartree−Fock theory overpolarizes the molecule, 3190

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation

Figure 6. Graphic depiction of surface delocalization of carbon monoxide: (left) equilibrium geometry, CASSCF(10,10); (middle) stretched bond, CASSCF(10,10); and (right) stretched bond, Hartree−Fock. D(r⃗) is plotted from 2.5 bohr (red) to 3.5 bohr (blue).

Table 5. Atomic Hirshfeld Charges (QC),a Delocalization Lengths (DC and DO), and Kinetic Energies Per Electron (KEC and KEO) for Equilibrium and Stretched Carbon Monoxide Delocalization Length (bohr) method

a

QC (e)

HF LDA B3LYP CCSD CASSCF

0.162 0.082 0.105 0.106 0.105

HF LDA B3LYP CCSD CASSCF

0.359 0.215 0.266 0.181 0.124

DC

DO

Equilibrium Bond 1.371 1.062 1.392 1.066 1.383 1.064 1.368 1.058 1.368 1.059 Stretched Bond 1.462 1.142 1.518 1.156 1.508 1.157 1.491 1.133 1.502 1.129

3.6. Chemical Reactions. We conclude by illustrating the surface delocalization of a chemical reaction. We revisit our recent theoretical and experimental study103 of phosphinylidenes R1R2HPO. These compounds display prototropic tautomerism between a stable P(V) form R1R2P(O)H and a reactive P(III) form R1R2P−OH.104 The uncatalyzed reaction’s high barrier105 is reduced by dimerization106 or in complexes with bicyclic guanidines.107 Figure 7 shows the computed reactant, transition state, and product of uncatalyzed P(V) → P(III) tautomerization of dimethyl phosphine oxide. These B3LYP/6-311++G(3df,3pd) calculations follow ref 103. The figure compares the molecular electrostatic potential and surface delocalization. The reactant’s electrostatic potential highlights the negatively charged oxygen atom (red), and the surface delocalization highlights the localized (red) electrons on this electronegative atom. The product’s electrostatic potential highlights the positively charged O−H hydrogen (blue) and modestly negatively charged phosphorus lone pair (yellow). In contrast, the product’s surface delocalization clearly highlights the delocalized phosphorus lone pair (blue), while showing that both O and H atoms on the OH group are quite localized (red), despite their different electrostatic potentials. Tautomerization thus converts the reactant’s modestly delocalized (green) P−H bond into a highly delocalized (blue) phosphorus lone pair and a localized (red) O−H bond. The transition state’s electrostatic potential shows that the forming O−H and breaking P−H bonds only exhibit modest charges (yellow and light blue). In contrast, the transition state’s surface delocalization clearly highlights the delocalized, stretched P−H and O−H bonds (blue). Table 6 shows the atomic delocalization and Hirshfeld charges of selected atoms in Figure 7. The results quantify the trends in Figure 7. The reaction localizes the transferred hydrogen, with DH being relatively small in the product. The reaction delocalizes the phosphorus, consistent with the product’s delocalized lone pair. In contrast, the DO value of oxygen is largely unchanged during the reaction, despite the substantial changes in charge QO. As expected, the reaction does not significantly change the charge or delocalization of the spectator CH3 groups. Overall, we conclude that the surface and atomic delocalization complement surface electrostatic potentials and atomic charges to provide a nuanced picture of this reaction.

Kinetic Energies Per Electron (Hartree) KEC

KEO

6.53 6.37 6.46 6.49 6.48

9.13 9.16 9.21 9.23 9.23

6.64 6.40 6.51 6.46 6.39

8.91 9.00 9.00 9.13 9.18

QC = − QO.

giving QC values larger than CASSCF or CCSD. The LDA overcorrects the HF atomic charges, while B3LYP charges approach the correlated methods. Hartree−Fock and DFT atomic delocalization lengths are larger than the accurate calculations. This is consistent with an admixture of highly oscillatory “virtual” orbitals into the accurate one-particle density matrices, and the resulting positive kinetic energy of correlation. Atomic kinetic energies per electron generally follow this trend, with the correlated calculations having relatively large kinetic energies per atom. At the stretched bond length, the surface delocalization shows how multireference calculations ”smooth out” the Hartree−Fock delocalization lengths. Atomic delocalization lengths and kinetic energies are larger and smaller, respectively, than the equilibrium values, which is consistent with the increased delocalization expected in stretched bonds. Hartree− Fock theory again overpolarizes the bond, and even the LDA is insufficient to correct the overpolarization. CCSD results are quite different from CASSCF, consistent with the T 1 diagnostic102 (>0.25), indicating significant multireference character. Just as for equilibrium bond lengths, the oxygen atom delocalization DO is significantly overestimated by the mean-field calculations. Mean-field theory’s overpolarization of carbon leads to tighter binding of the remaining electrons, an error cancellation improving agreement with accurate DC. We conclude that surface and atomic delocalization lengths provide a convenient way to quantify correlation effects in stretched bonds.

4. CONCLUSION This work illustrates how the electron delocalization range can visualize and quantify delocalization in computed wave functions. The surface delocalization provides an “orbital lobe size”, capturing aspects of delocalization and chemical hardness, 3191

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation

Figure 7. Graphical depiction of the reactant, transition state, and product of uncatalyzed P(V) → P(III) tautomerization of dimethyl phosphine oxide: (top) molecular geometry; (middle) density isosurface ρ(r⃗) = 0.001 electrons/bohr3 colored with molecular electrostatic potential from −0.08 a.u. (red) to +0.09 a.u. (blue); and (bottom) same density isosurface colored with delocalization D(r⃗) from 2.6 bohr (red) to 3.6 bohr (blue).

new route to building chemical insight from electronic structure calculations.

Table 6. Atomic Charges (QA) and Valence Delocalization Lengths (DA) for Selected Atoms (See Figure 7) propertya Atomic Charge (e) QP QO QC QH Valence Delocalization Length (bohr) DP DO DC DH a

reactant

transition state

product

0.455 −0.449 −0.166 −0.031

0.304 −0.388 −0.170 0.076

0.208 −0.294 −0.183 0.169

1.837 1.348 1.734 2.074

1.920 1.351 1.739 2.048

1.994 1.327 1.742 1.705



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00343. Basis set effects on Li atom delocalization, valence vs allelectron delocalization of Be atom, H and H+2 density isosurfaces colored with LOL, all-electron atomic delocalization lengths and kinetic energies, surface delocalization of small molecules, benzene surface delocalization with other methods and density isosurfaces, benzene EDR along symmetry axes, geometries and energies of cation-pentadecene complexes, CO equilibrium surface delocalization, and analysis of model systems with diffuse electron densities (PDF)

“H” denotes the transferred hydrogen.



capable of distinguishing hard and soft sites on a single molecule (Figure 3) and visualizing effects of electron correlation (Figure 6). The atomic delocalization quantifies these trends. Analysis of a chemical reaction shows how the surface and atomic delocalization can be combined with other tools to provide a richer picture of chemical reactivity. We suggest that the surface and atomic delocalization provide a

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 3192

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation



(39) Proynov, E.; Liu, F.; Kong, J. Phys. Rev. A: At., Mol., Opt. Phys. 2013, 88, 032510-1−032510-8. (40) Alcoba, D. R.; Lain, L.; Torre, A.; Bochicchio, R. C. J. Chem. Phys. 2005, 123, 144113-1−144113-7. (41) Vanfleteren, D.; Van Neck, D.; Bultinck, P.; Ayers, P. W.; Waroquier, M. J. Chem. Phys. 2010, 132, 164111-1−164111-10. (42) Vanfleteren, D.; Van Neck, D.; Bultinck, P.; Ayers, P. W.; Waroquier, M. J. Chem. Phys. 2012, 136, 014107-1−014107-13. (43) Bader, R. W. F.; Preston, H. J. T. Int. J. Quantum Chem. 1969, 3, 327−347. (44) Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990, 92, 5397− 5403. (45) Savin, A. J. Mol. Struct.: THEOCHEM 2005, 727, 127−131. (46) Kohout, M. Int. J. Quantum Chem. 2004, 97, 651−658. (47) Matito, E.; Duran, M.; Solá, M. J. Chem. Educ. 2006, 83, 1243− 1248. (48) Feixas, F.; Matito, E.; Duran, M.; Solá, M.; Silvi, B. J. Chem. Theory Comput. 2010, 6, 2736−2742. (49) Ayers, P. W. Proc.Indian Acad. Sci., Chem. Sci. 2005, 117, 441−454. (50) Savin, A. J. Phys. Chem. Solids 2004, 65, 2025−2029. (51) Burdett, J. K.; McCormick, T. A. J. Phys. Chem. A 1998, 102, 6366−6372. (52) Wilson, C. W., Jr.; Goddard, W. A., III. Theor. Chim. Acta 1972, 26, 195−210. (53) Goddard, W. A., III; Wilson, C. W., Jr. Theor. Chim. Acta 1972, 26, 211−230. (54) Schmider, H. L. J. Chem. Phys. 1996, 105, 11134−11142. (55) Schmider, H. L.; Becke, A. D. J. Mol. Struct.: THEOCHEM 2000, 527, 51−61. (56) Schmider, H. L.; Becke, A. D. J. Chem. Phys. 2002, 116, 3184− 3193. (57) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Phys. 2014, 141, 144104-1−144104-13. (58) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Phys. 2014, 141, 034103-1−034103-14. (59) Scemama, A.; Caffarel, M.; Chaudret, R.; Piquemal, J.-P. J. Chem. Theory Comput. 2011, 7, 618−624. (60) Janesko, B. G.; Scuseria, G. E. J. Chem. Phys. 2007, 127, 164117. (61) Janesko, B. G.; Scalmani, G.; Frisch, M. J. Phys. Chem. Chem. Phys. 2015, 17, 18305−18317. (62) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Theory Comput. 2016, 12, 79−91. (63) Janesko, B. G. Topological analysis of the electron delocalization range, J. Comput. Chem., 2016, DOI: 10.1002/jcc.24421. (64) Williams, C. F.; Herbert, J. M. J. Phys. Chem. A 2008, 112, 6171−6178. (65) Politzer, P. A.; Murray, J. Theor. Chem. Acc. 2002, 108, 134− 142. (66) Bader, R. W. F.; Preston, H. J. T. Theor. Chim. Acta 1970, 17, 384−395. (67) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129−138. (68) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Bloino, J.; Janesko, B. G.; Izmaylov, A. F.; Marenich, A.; Lipparini, F.; Zheng, G.; Sonnenberg, J. L.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Parandekar, P. V.; Mayhall, N. J.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.;

ACKNOWLEDGMENTS B.G.J. was supported by the National Science Foundation (No. DMR-1505343).



REFERENCES

(1) Nordholm, S. J. Chem. Educ. 1988, 65, 581−584. (2) Autschbach, J. J. Chem. Educ. 2012, 89, 1032−1040. (3) Boys, S. F. In Quantum Theory of Atoms, Molecules, and the Solid State; Löwdin, P. O., Ed.; Academic Press: New York, 1966; pp 253− 262. (4) Edmiston, C.; Ruedenberg, K. In Quantum Theory of Atoms, Molecules, and the Solid State; Löwdin, P. O., Ed.; Academic Press: New York, 1966; pp 263−280. (5) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735−746. (6) Cremer, D. Mol. Phys. 2001, 99, 1899−1940. (7) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (8) Abramov, Y. A. Acta Crystallogr., Sect. A: Found. Crystallogr. 1997, A53, 264−272. (9) Karasiev, V. V.; Jones, R. S.; Trickey, S. B.; Harris, F. E. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 245120-1−245120-17. (10) Bader, R. F. W.; Streitwieser, A.; Neuhaus, A.; Laidig, K. E.; Speers, P. J. Am. Chem. Soc. 1996, 118, 4959−4965. (11) Coulson, C. A. Proc. R. Soc. London, Ser. A 1939, 169, 413−428. (12) Mulliken, R. S. J. Chem. Phys. 1955, 23, 2343−2346. (13) Wiberg, K. B. Tetrahedron 1968, 24, 1083. (14) Mayer, I. Chem. Phys. Lett. 1983, 97, 270−274. (15) Yang, C. N. Rev. Mod. Phys. 1962, 34, 694−704. (16) Tretiak, S.; Mukamel, S. Chem. Rev. 2002, 102, 3171−3212. (17) Head-Gordon, M.; Graña, A. M.; Maurice, D.; White, C. A. J. Phys. Chem. 1995, 99, 14261−14270. (18) Plasser, F.; Wormit, M.; Dreuw, A. J. Chem. Phys. 2014, 141, 024106-1−024106-13. (19) Gadre, S. R.; Kulkarni, S. A.; Pathak, R. K. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 40, 4224−4231. (20) Schmider, H.; Edgecombe, K. E.; Smith, V. H., Jr.; Weyrich, W. J. Chem. Phys. 1992, 96, 8411−8419. (21) Bader, R. F. W. Atoms In MoleculesA Quantum Theory; Oxford University Press: Oxford, U.K., 1990; pp 1−20. (22) Silvi, B. Phys. Chem. Chem. Phys. 2004, 6, 256−260. (23) Bader, R. F. W.; Stephens, M. E. J. Am. Chem. Soc. 1975, 97, 7391−7399. (24) Cioslowski, J.; Mixon, S. T. J. Am. Chem. Soc. 1991, 113, 4142− 4145. (25) Fulton, R. L.; Mixon, S. T. J. Phys. Chem. 1993, 97, 7530−7534. (26) Ponec, R.; Cooper, D. L. J. Mol. Struct.: THEOCHEM 2005, 727, 133−138. (27) Bader, R. W. F.; Heard, G. L. J. Chem. Phys. 1999, 111, 8789− 8798. (28) Ponec, R.; Roithová, J. Theor. Chem. Acc. 2001, 105, 383−392. (29) Ponec, R.; Cooper, D. L.; Savin, A. Chem.Eur. J. 2008, 14, 3338−3345. (30) Bochicchio, R. C.; Lain, L.; Torre, A. Chem. Phys. Lett. 2003, 374, 567−571. (31) Fulton, R. L.; Mixon, S. T. J. Phys. Chem. 1995, 99, 9768−9783. (32) Bultinck, P.; Cooper, D. L.; Ponec, R. J. Phys. Chem. A 2010, 114, 8754−8763. (33) Kent, P. R. C.; Hood, R. Q.; Towler, M. D.; Needs, R. J.; Rajagopal, G. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 15293−15302. (34) Takatsuka, K.; Fueno, T.; Yamaguchi, K. Theor. Chim. Acta (Berl.) 1978, 48, 175−183. (35) Proynov, E. I. J. Mol. Struct.: THEOCHEM 2006, 762, 159−163. (36) Becke, A. D.; Roussel, M. R. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 39, 3761−3767. (37) Becke, A. D. J. Chem. Phys. 2003, 119, 2972−2977. (38) Becke, A. D. J. Chem. Phys. 2005, 122, 064101-1−064101-6. 3193

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194

Article

Journal of Chemical Theory and Computation Fox, D. J. Gaussian Development Version, Revision I.02+; Gaussian, Inc.: Wallingford, CT, 2014. (69) Kohn, W.; Sham, L. Phys. Rev. 1965, 140, A1133−A1138. (70) Seidl, A.; Görling, A.; Vogl, P.; Majewski, J. A.; Levy, M. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 3764−3774. (71) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (72) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (73) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (74) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (75) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (76) Chai, J.-D.; Head-Gordon, M. J. Chem. Phys. 2008, 128, 0841061−084106-15. (77) Handy, N. C.; Schaefer, H. F., III. J. Chem. Phys. 1984, 81, 5031−5033. (78) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. J.; Frisch, M. J. J. Phys. Chem. 1992, 96, 671−679. (79) Krishnan, R.; Binkley, J.; Seeger, R.; Pople, J. J. Chem. Phys. 1980, 72, 650−654. (80) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007−1023. (81) Jorge, F. E.; de Castro, E. V. R.; da Silva, A. B. F. Chem. Phys. 1997, 216, 317−321. (82) Perdew, J. P.; Levy, M. Phys. Rev. Lett. 1983, 51, 1884−1887. (83) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. J. Chem. Phys. 1997, 107, 5007−5015. (84) Ruedenberg, K. Rev. Mod. Phys. 1962, 34, 326−376. (85) Ghanty, T. K.; Ghosh, S. K. J. Phys. Chem. 1996, 100, 17429− 17433. (86) Davidson, E. R.; Chakravorty, S. Theor. Chem. Acc. 1992, 83, 319−330. (87) Berski, S.; Latajka, Z.; Gordon, A. J. J. Comput. Chem. 2010, 31, 2555−2567. (88) Pearson, R. G. J. Am. Chem. Soc. 1963, 85, 3533−3539. (89) Xie, J.-H.; Liu, X.-Y.; Yang, X.-H.; Xie, J.-B.; Wang, L.-X.; Zhou, Q.-L. Angew. Chem., Int. Ed. 2012, 51, 201−203. (90) Williams, D. B. G.; Traut, T.; Kriel, F. H.; van Zyl, W. E. Inorg. Chem. Commun. 2007, 10, 538−542. (91) Torrent-Sucarrat, M.; De Proft, F.; Geerlings, P.; Ayers, P. W. Chem.Eur. J. 2008, 14, 8652−8660. (92) Torrent-Sucarrat, M.; De Proft, F.; Ayers, P. W.; Geerlings, P. Phys. Chem. Chem. Phys. 2010, 12, 1072−1080. (93) Parr, R. G.; Pearson, R. G. J. Am. Chem. Soc. 1983, 105, 7512− 7516. (94) Pearson, R. G. J. Chem. Educ. 1968, 45, 581−587. (95) Pearson, R. G. J. Chem. Educ. 1968, 45, 643−648. (96) Ghanty, T. K.; Ghosh, S. K. J. Phys. Chem. 1993, 97, 4951−4953. (97) Brédas, J.-L.; Beljonne, D.; Coropceanu, V.; Cornil, J. Chem. Rev. 2004, 104, 4971−5003. (98) Vijay, D.; Sastry, G. N. Phys. Chem. Chem. Phys. 2008, 10, 582− 590. (99) Prodan, E.; Kohn, W. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 11635−11638. (100) Bader, R. W. F. J. Phys. Chem. A 2008, 112, 13717−13728. (101) Grimaldi, F.; Lecourt, A.; Moser, C. Int. J. Quantum Chem. 1967, 1, 153−161. (102) Lee, T. J.; Taylor, P. R. Int. J. Quantum Chem. 1989, 36, 199− 207. (103) Janesko, B. G.; Fisher, H. C.; Bridle, M. A.; Montchamp, J.-L. J. Org. Chem. 2015, 80, 10025−10032. (104) Montchamp, J.-L. Acc. Chem. Res. 2014, 47, 77−87. (105) Wesolowski, S. S.; Brinkmann, N. R.; Valeev, E. F.; Schaefer, H. F., III; Repasky, M. P.; Jorgensen, W. L. J. Chem. Phys. 2002, 116, 112−122. (106) Mamaev, V. M.; Prisyajnuk, A. V.; Laikov, D. N.; Logutenko, L. S.; Babin, Y. V. Mendeleev Commun. 1999, 9, 240−241.

(107) Cho, B.; Tan, C.-H.; Wong, M. W. Org. Biomol. Chem. 2011, 9, 4550−4557.

3194

DOI: 10.1021/acs.jctc.6b00343 J. Chem. Theory Comput. 2016, 12, 3185−3194