Quantifying Electron Delocalization in Electrides - ACS Publications

Article pubs.acs.org/JCTC

Quantifying Electron Delocalization in Electrides Benjamin G. Janesko,*,† Giovanni Scalmani,‡ and Michael J. Frisch‡ †

Department of Chemistry, Texas Christian University, Fort Worth, Texas 76129, United States Gaussian, Inc., 340 Quinnipiac Street Building 40, Wallingford, Connecticut 06492, United States

‡

S Supporting Information *

ABSTRACT: Electrides are ionic solids whose anions are electrons confined to crystal voids. We show that our electron delocalization range function EDR(r;⃗ d), which quantifies the extent to which an electron at point r ⃗ in a calculated wave function delocalizes over distance d, provides useful insights into electrides. The EDR quantifies the characteristic delocalization length of electride electrons and provides a chemically intuitive real-space picture of the electrons’ distribution. It also gives a potential diagnostic for whether a given formula unit will form a solid electride at ambient pressure, quantifies the effects of electron−electron correlation on confined electrons’ interactions, and highlights analogies between covalent bonding and the interaction of interstitial quasi-atoms in high-pressure electrides. These results motivate adding the EDR to the toolbox of theoretical methods applied to electrides.

1. INTRODUCTION Electron delocalization is a central idea in chemistry, rationalizing a huge body of experimental work. The interplay of electron delocalization and confinement in anionic water clusters,1−3 solvated electrons,4,5 F-center defects,6 Rydberg states,7 and other systems provide unique and technologically relevant properties. Electrides, ionic salts whose anions are electrons trapped in crystal voids, are another important example of this interplay.8 Electrides’ trapped electrons can produce strongly variable conductivities 9 and magnetic stabilities8 and have potential applications including catalysis10−12 and thermionic electron emitters.13 Recent experimental demonstrations of room-temperature-stable electrides14,15 and recent computational16−24 and experimental25−28 evidence for electride behavior at high pressure make fundamental investigations of electrides particularly timely. 1.1. Previous Theoretical Studies of Electrides. Quantum-mechanical calculations are essential for providing insight into electrides, as the anions are generally too diffuse to be located in X-ray crystal structures.8 There is extensive literature on such calculations. Here we briefly highlight some studies relevant to the present work. Singh and co-workers simulated Cs+(15-crown-5)2e− using periodic local density approximation calculations, characterizing the electron distribution in interstitial regions.29 Dye and co-workers used van der Waals isosurfaces to visualize electrides’ void spaces.30 Several groups have studied high-pressure electrides using periodic density functional calculations.16−22 Other groups have proposed candidate electrides based on simulations of single formula units of alkali metal adducts.31−35 Rousseau and Ashcroft used an electron-and-hard-sphere model to build insight into high-pressure electrides.36 Miao and Hoffmann considered a related model, in which interstitial quasiatoms (ISQ) are modeled as anionic vacancies in solid helium under pressure.23,24 Dale and co-workers performed a systematic evaluation of known ambient-pressure electrides’ band © XXXX American Chemical Society

structures and characterized the trapped electrons’ structure using valence band densities, van der Waals isosurfaces, the noncovalent interaction index (NCI, ref 37), and the quantum theory of atoms in molecules (QTAIM, ref 38).39 Kimoon Lee and co-workers used band structure calculations and the electron localization function (ELF, refs 40 and 41) to assign Ca2N as a two-dimensional electride.42,43 Postils and coworkers recently defined ”molecular” electrides in terms of the topology of the electron density and ELF.44 Kumar and Gadre performed a topological analysis of the electrostatic potential in single formula units of experimentally verified solid electrides and other alkali adducts.45 1.2. The Electron Delocalization Range Function. We have proposed a new theoretical tool that explicitly addresses the question ”How far do electrons delocalize?”.46 Our electron delocalization range function EDR(r;⃗ d) quantifies the degree to which an electron at point r ⃗ in a computed wave function delocalizes over distance d (referred to as u in previous work46,47). The EDR is based on the nonlocal one-particle density matrix γ(r,⃗ r′)⃗ of a computed N-electron wave function Ψ(r1⃗ , r2⃗ ...rN⃗ )48 → ⎯ γ( r ⃗ , r′) ≡ N

→ ⎯

∫ d3r2⃗ ... d3rN⃗ Ψ( r ⃗ , r2⃗ ... rN⃗ )Ψ*(r′, r2⃗ ... rN⃗ ) (1)

γ(r,⃗ r′)⃗ characterizes delocalization between points r ⃗ and r′,⃗ as discussed in Section 1.3 below. Diagonal elements limr′→r⃗ ⃗ γ(r,⃗ r′)⃗ = ρ(r)⃗ give the probability density for finding an electron at r.⃗ The EDR quantifies the degree to which an electron at r ⃗ delocalizes over distance d. It does so by contracting γ(r,⃗ r′)⃗ with a test function that decays in distance |r− ⃗ r′|⃗ with characteristic decay length d: Received: August 12, 2015

A

DOI: 10.1021/acs.jctc.5b00993 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation EDR( r ⃗ ; d) =

→ ⎯

→ ⎯ → ⎯

∫ d3 r′gd( r ⃗ , r′)γ(r′, r ⃗)

→ ⎯ ⎞ ⎛ ⎛ 2 ⎞3/4 → ⎯ | r ⃗ − r′|2 ⎟ gd ( r ⃗ , r′) ≡ ρ−1/2 ( r ⃗)⎜ 2 ⎟ exp⎜⎜ − ⎝ πd ⎠ d 2 ⎟⎠ ⎝

covalent bond between atoms A and B typically gives an off⃗ B) large and diagonal density matrix element γ(r∈ ⃗ A, r′∈ positive, while an antibonding interaction between A and B ⃗ B) negative.57 Off-diagonal typically gives γ(r∈ ⃗ A, r′ ∈ delocalization is also connected to the electronic kinetic energy expectation value 2⟨T⟩ = ∫ d3r ⃗ limr′→r⃗ ⃗ ∇r⃗·∇r′γ ⃗ (r,⃗ r′)⃗ . Many existing theoretical tools quantify the extent of offdiagonal delocalization between pairs of atoms. These include Hückel bond order matrices,58,59 Mulliken overlap populations,60−62 Wiberg bond orders,63 and the delocalization index64−67 and related quantities.68,69 Other tools provide a real-space picture of delocalization through 2τ(r )⃗ ≡ limr′→r⃗ ⃗ ∇r⃗·∇r′γ ⃗ (r,⃗ r′)⃗ . These include the electron localization function evaluated from single-determinant wave functions40,41 and related functions.70−72 The EDR complements these two classes of tools by quantifying the length scale of off-diagonal delocalization. Appendix A shows how the EDR complements the delocalization index by predicting different characteristic delocalization lengths dmax for bonding and antibonding states evaluated from single-determinant wave functions. Appendix B shows that the EDR predicts that bulk metals’ dmax scales as the Fermi length k−1 F .

(2)

(3)

The prefactors in eq 3 ensure |EDR(r;⃗ d)| ≤ 1. Our choice of a Gaussian test function enables analytic integration over r′ ⃗ in eq 2, when the orbitals are expanded in standard atom-centered Gaussian basis sets.49,50 Global descriptors of delocalization may be obtained from the system average 2

⟨EDR(d)⟩ =

∫ d3r ρ⃗ ( r ⃗)EDR( r ⃗; d)

(4)

We use such averages to define the characteristic delocalization dmax of a finite N-electron system’s most weakly bound electron. dmax is the length scale d that maximizes the difference between ⟨EDR(d)⟩ of the N-electron and (N-1)-electron systems: dmax = arg max Δ⟨EDR(d)⟩

(5)

Δ⟨EDR(d)⟩ = ⟨EDR(d)⟩(N ) − ⟨EDR(d)⟩(N − 1)

(6)

d

The corresponding difference between total energies defines the vertical detachment energy VDE. We have shown46,47 that the EDR provides a novel perspective on electron delocalization and complements other theoretical tools.51 The EDR gives a chemically intuitive picture of delocalization in molecules: EDR(r;⃗ d) at small distances d ∼ 0.2 bohr peaks at points r ⃗ in the cores of first-row atoms, while EDR(r;⃗ d) at larger distances d ∼ 1 bohr peaks at points r ⃗ near covalent bonds. We recently applied the EDR to electrons solvated in finite clusters, e.g., (H2O)−N.47 Solvated electrons’ characteristic delocalization lengths dmax, calculated from Hartree−Fock density matrices, have a one-to-one correspondence with the size of the Hartree−Fock highest occupied molecular orbital (HOMO).52 Real-space pictures of EDR(r;⃗ dmax) show that it highlights the same region of space as the major lobe of the HOMO and spin density. The EDR also quantifies the relations between delocalization and confinement: for example, dmax of an electron confined in the octahedral (H2O)−6 Kevan53 cavity increases asymptotically linearly with cavity radius. The density-matrix-based EDR is independent of rotations among occupied orbitals and can be evaluated for any calculation that produces a density matrix. Our study of ammoniated electrons’ insulator-to-metal transition illustrated that the EDR quantifies relationships between delocalization and strong (multireference) correlation.47 These results motivate the present study, which uses the EDR to visualize and quantify the interplay of delocalization, confinement, and correlation in electrides. 1.3. Diagonal vs Off-Diagonal Delocalization. Before continuing, we specify the ”delocalization” considered here. We follow ref 54 in distinguishing the diagonal delocalization of charge or spin over multiple centers, from the electron coherence (off-diagonal delocalization) distinguishing bonding from antibonding interactions. Diagonal delocalization is a property of the diagonal density matrix γ(r,⃗ r)⃗ = ρ(r)⃗ . Examples include the diagonally delocalized charge of sodium acetate’s two crystallographically equivalent Oδ− atoms55 and the diagonally delocalized electron spin of C6H6− giving six equivalent electron−proton couplings.56 Off-diagonal delocal⃗ r)⃗ . A ization depends on the off-diagonal density matrix γ(r,⃗ r′≠

2. METHODOLOGY We consider four sets of model systems. The first set of model systems is single formula units of four experimentally realized electrides Cs + (15C5) 2 e − (15C5 = 15-crown-5), 7 3 K + (Cryptand-2.2.2)e − , 74 Li + (Cryptand-2.1.1)e − , 75 and Na+(tripip-aza-2.2.2)e−15 and a full unit cell of Li+(Cryptand2.1.1)e−. Geometries are taken from experimental crystal structures. Two CH2 groups in Na+(tripip-aza-2.2.2)e− have C−H bond lengths modified to 1.07 Å; the geometry used is reported in the Supporting Information. Calculations at the B3LYP/6-31++G(d,p) optimized geometry of a single formula unit reduce dmax from 17.3 to 16.4 bohr, a modest shift that maintains the ordering in Table 1. Table 1. B3LYP/6-31++G(d,p) VDE (eV), and dmax (bohr)a system

VDE

dmax

K+(Cryptand-2.2.2)e− Na+(tripip-aza-2.2.2)e− Cs+(15C5)2e− Li+(Cryptand-2.1.1)e−

1.92 1.98 1.73 2.12

18.7 17.2 16.1 15.8

·− Na·+ 2 TCNQ ·− Li·+ TCNQ 2 Na+(Calix)e− Li+(Calix)e− Li(B10H14)

6.16 6.27 3.57 3.88 6.74

5.4 4.1 4.2 3.9 2.4

5.33 4.37 5.62 5.94 11.28 6.96 9.68

7.7 6.3 6.0 5.8 2.3 2.3 1.9

Li2 Li-NCH Li Li-HCN acetylene C10H12 C10H22 a

Results for single formula units of four experimentally realized electrides, five alkali adducts theoretically proposed as electrides, and seven ”sanity tests”.

B


Article

Journal of Chemical Theory and Computation

Figure 1. B3LYP spin density for K+(Cryptand-2.2.2)e−, Na+(tripip-aza-2.2.2)e−, Cs+(15C5)2e−, and Li+(Cryptand-2.1.1)e−. Results plotted at contour 0.000024 bohr−3 for K+(Cryptand-2.2.2)e−, 0.001 bohr−3 for Cs+(15C5)2e−, and 0.004 bohr−3 otherwise.

corresponding to the pressure range from P = 0.001 GPa (He−He separation 2.873 Å) and P = 500 GPa (He−He separation 1.390 Å) reported in ref 23. At He−He separations below 1.3 Å, triplet calculations on our two-ISQ system have one electron escape the finite cluster to reside on the cluster surface. Other computational details are as follows. All calculations use the development version of the GAUSSIAN suite of programs.77 Calculations on open-shell systems are performed spin-unrestricted. The EDR is evaluated as previously described.46,78 We quantify dmax as in eqs 5 and 6. Calculations on systems with D delocalized electrons use ⟨EDR(d)⟩(N−D) in place of ⟨EDR(d)⟩(N−1) in eq 6. We use a Li−C4H4N2··· Na2+ 2 reference (D = 2) for Li−C4H4N2···Na2, a + 6 charged reference (D = 6) for six formula units of Li+(Cryptand2.1.1)e−, a He133 reference (D = 2) for two interacting ISQ in the He2− 133 ”compression chamber”, and a reference with no electrons (D = 2) for H2. Plots of the EDR evaluate it from the majority-spin density matrix. ⟨EDR(d)⟩ are evaluated summing over spins. We fit dmax using a three-point fit to ⟨EDR(di) ⟩ evaluated for multiple length scales {di}. Most calculations use 35 values di = (1.2)i/2 bohr, i = 0,1,2....

The second set of model systems is single formula units of alkali atom adducts treated in other theoretical works.31−35 We consider the pyrroles Li+(Calix)e− and Na+(Calix)e− (Calix = ·− Calix[4]pyrrole),32 the triplet radical ion pair salts Li·+ 2 TCNQ 34 ·+ ·− and Na2 TCNQ (TCNQ = tetracyanoquinodimethane), and Li(B10H14).33 Geometries are taken from ref 35. We also consider the alkali pyrazine Li−C4H4N2···Na2 evaluated at MP2/6-311+G(d,p) geometries.76 The third set of model systems provides ”sanity tests” of our predictions. We consider two lithium cyanide adducts,31 Li atom and Li2,44 the small molecule acetylene, and the large molecules decane and πconjugated 1,3,5,7,9-pentadecene C10H12. Cyanide adducts LiHCN and Li-NCH are evaluated at geometries from ref 31. Others use B3LYP/6-31G(d) geometries. The fourth set of model systems is one or two interstitial quasiatoms (ISQ)23 in pressurized helium. This is similar to the ”helium compression chamber” of refs 23 and 24; however, periodic boundary conditions are not used, and the atoms are frozen to crystal positions. Single ISQ are simulated at the center of a finite 135-atom supercell of fcc He. Two adjacent ISQ are simulated by adding a second ISQ. We scan over He− He separations between 2.8 and 1.3 Å, approximately C


Article


Table 2. Dynamical Correlation Effects on Delocalizationa

Most calculations use the modest 6-31++G(d,p) Pople basis set49,50 used in ref 45. Reference 35 showed that a somewhat smaller basis provides reasonable (hyper)polarizabilities for single formula units of alkali atom adducts. Calculations on Cs atom use the LANL2DZ basis set and effective core potential. Calculations on a full unit cell of Li+(Cryptand-2.1.1)e− use the 3-21++G basis for computational convenience. Calculations on the He compression chamber use the 6-31G(d,p) basis for the He atoms and place the He atom aug-cc-pVQZ79 basis functions at the center of each ISQ. H2 calculations use the aug-cc-pVQZ basis set. Calculations on Li−C4H4N2···Na2 use the 6-31+G(d,p) basis set. Calculations in Appendix A use the minimal STO-3G80 and large aug-cc-pVQZ basis sets and a single Gaussian function with exponent 2 bohr−2. Calculations of dynamical correlation effects compare uncorrelated Hartree−Fock theory (HF) to second-order many-body perturbation theory (MP2) or density functional theory with the B3LYP exchange-correlation functional.81,82 MP2 density matrices are obtained using the Z-vector method83,84 with Gaussian keyword ”Density = Current”. Appendix C discusses some technical issues with these calculations. Calculations of nondynamical correlation effects in singlet Li−C4H4N2···Na2 compare restricted Hartree−Fock (RHF), symmetry-broken spin-unrestricted Hartree−Fock (UHF), and complete active-space self-consistent field calculations with two active electrons in two orbitals (CAS(2,2)). The (N-2)-electron reference system Li−C4H4N2··· Na2+ 2 is treated with RHF. CAS calculations use the converged RHF molecular orbitals as an initial guess. Calculations of nondynamical correlation in two singlet-coupled ISQ compare RHF and UHF. Calculations on the singly- and doubly excited singlet states of minimal-basis H2 in Appendix A use CAS(2,2), i.e., full CI in the minimal basis. We evaluate the bonding energy of a pair of adjacent ISQ using a strain energy correction following ref 24. We define the bonding energy Ebond as the difference between the total ISQISQ interaction energy Etotal and the strain energy Estrain. Denoting the defect-free cluster as He135, the cluster with a singlet ISQ as He134−, the cluster with two adjacent ISQ as He2− 133, and the corresponding neutral clusters (without bound electrons) as He134 and He133, we have

ΔVDE system +

Li (Cryptand-2.1.1)e ·− Na·+ 2 TCNQ ·− ·+ Li2 TCNQ Na+(Calix)e− Li+(Calix)e− Li(B10H14)

B3LYP

MP2

B3LYP

0.20 0.36 0.38 0.48 0.50 0.71

0.48 0.65 0.46 0.67 0.69 0.90

−0.27 −0.23 −0.21 −0.31 −0.23 0.09

1.20 −0.15 −0.63 −0.61 −0.49 0.06

ΔVDE (eV) and Δdmax (bohr) for MP2 and B3LYP calculations, relative to Hartree-Fock theory.

a

and VDE 6.55 eV (Table 2) are close to the MP2/6-31+G(d) values 3.59 au and 6.12 eV reported in ref 33. The UB3LYP/631++G(d,p) average static electronic polarizabilities of Li+(Calix)e−, Na+(Calix)e−, Li(B10H14), Li2·+TCNQ·−, and ·− Na·+ are 407.08, 450.73, 150.07, 340.91, and 376.26 2 TCNQ au, largely consistent with the hybrid basis results in Table 1 of ·− ref 35. UB3LYP calculations on Li·+ with the mixed 2 TCNQ basis set of ref 35 (6-31+G(d) on Li, 6-31+G on other atoms), at geometries optimized with that basis set, yield a polarizability 371.94 au matching Table 1 of that reference. The Li− C4H4N2···Na2 Na−N distances are 2.529 and 2.587 Å, and the Na−Na distance is 3.204 Å, within 0.001 Å of the values in ref 76. The highest occupied orbital energy for a single ISQ at He− He separation 1.4 Å is +8.0 eV, notably lower than the ∼+28 eV shown for P = 500 GPa (He−He separation 1.39 Å) in Figure 3 of ref 23. The strain energy of two ISQ at this separation is 0.06 eV, lower than the 0.3 eV at P = 500 GPA in Figure 7 of ref 24. We attribute these to our different choice of basis set, periodic boundary conditions, background charge, and exchange-correlation functional.

3. RESULTS 3.1. Introduction and Comparison with Other Theoretical Tools. Figures 1 and 2 show the spin density and majority-spin EDR(r;⃗ dmax) from B3LYP/6-31++G(d,p) calculations on single formula units of the experimentally realized solid electrides Cs+(15C5)2e−, K+(Cryptand-2.2.2)e−, Li+(Cryptand-2.1.1)e−, and Na+(tripip-aza-2.2.2)e−. Supporting Information Figure SI-1 shows that the HOMO is essentially indistinguishable from the spin density. To improve readability, particularly diffuse electrons are plotted with transparent isosurfaces. The most important result is that EDR(r;⃗ dmax), the real-space distribution of electrons delocalized over characteristic distance dmax, highlights the same region of space as the major lobe of the spin density. This confirms our previous results.47 The EDR predicts that the most delocalized electron is delocalized over the outer surfaces of K+(Cryptand2.2.2)e−, Na+(tripip-aza-2.2.2)e−, and Cs+(15C5)2e− and bound at the top of Li+(Cryptand-2.1.1)e−. 3.2. Delocalization in Electrides, Alkali Adducts, and Normal Molecules. Table 1 shows the B3LYP VDE and characteristic delocalization length dmax for single formula units of the four experimentally realized electrides in Section 3.1. Results are also presented for five alkali metal adducts previously considered for electride behavior31−35 and seven ”sanity tests” discussed above. Supporting Information Figure SI-2 illustrates several of these systems’ Δ⟨EDR(d)⟩. The computed dmax follow chemically reasonable trends. Weakly bound electrons with small VDE are quite delocalized and have large dmax. Tightly bound electrons with large VDE are more

2− − Etotal = (E(He133 ) − E(He135)) − 2(E(He134 ) − E(He135))

(7)

Estrain = (E(He133) − E(He135)) − 2(E(He134) − E(He135)) (8)

Ebond = Etotal − Estrain

−

Δdmax

MP2

(9)

Before discussing our results, we confirm that our calculations are consistent with previous work. The computed HOMO of Li+(Calix)e− (Supporting Information Figure SI-3) is qualitatively similar to Figure 1 of ref 32. The computed B3LYP/6-31++G(d,p) IP 3.88 eV is qualitatively consistent with the 4.16 eV of ref 32. B3LYP calculations in that reference’s basis set (6-311++G(3df,3pd) on Li, 6-311++G elsewhere) give a IP 3.92 eV. The MP2 VDE 6.19 eV ·− ·+ ·− (Li·+ 2 TCNQ ) and 5.87 eV (Na2 TCNQ ) (Table 2) are quite close to the 6.19 and 5.80 eV restricted open-shell values reported in ref 34. Calculations with additional diffuse basis functions on hydrogen52,85,86 give negligible changes to the MP2 dmax, changing the Li+(Calix)e− dmax from 2.218 to 2.226 Å and changing the Na+(Calix)e− dmax from 2.418 to 2.429 Å. The Li(B10H14) MP2/6-31++G(d,p) dipole moment 3.11 au D


Article


Figure 2. B3LYP EDR(r⃗; dmax) for the systems in Figure 1. Results plotted at contour 0.7 for Li+(Cryptand-2.1.1)e−, 0.6 otherwise.

localized and have smaller dmax. Replacing Li with Na in the calixarene and TCNQ electrides increases dmax consistent with Na atom’s increased atomic radius and larger valence delocalization. One important result in Table 1 is that the chemically stable ”normal” molecules acetylene, decane, and pentadecene have the largest VDE and the smallest dmax. This result is consistent with the virial theorem, which predicts that the total and electronic kinetic energy expectation values obey ⟨E⟩ = −⟨T⟩ for molecules at optimized geometries. Chemically stable, lowenergy systems will thus tend to have relatively large groundstate electronic kinetic energies, with corresponding small delocalization lengths dmax. (For example, formation of the stable covalent bond in H2 increases the electronic kinetic energy relative to separated H atoms, a result consistent with the coherent off-diagonal delocalization of the bonding interaction.)57 Conjugated pentadecene, like the metals in Appendix B, has a small dmax consistent with its chemical stability. Another important result in Table 1 is that single formula units of the four experimentally realized electrides have substantially smaller VDE, and substantially larger dmax, than any of the other tested systems. Their dmax > 8 Å are in the

range we reported for electrons on the surface of water clusters.47 This difference between the experimentally realized electrides and the other alkali adducts in Table 1 is consistent with the fact that, to our knowledge, crystals of the latter have not been demonstrated to contain interstitial electrons at ·− ambient pressure. For example, Na·+ 2 TCNQ is a colorless salt 87 without low-energy excited states. This difference is also consistent with the recent study by Kumar and Gadre, who computed the molecular electrostatic potential on single formula units of the experimentally realized electrides Li+(Cryptand-2.1.1)e−, K+(Cryptand-2.2.2)e−, and Na+(tripipaza-2.2.2)e − and the alkali adducts Li + (Calix)e − and Na+(Calix)e−. These authors found that the distance between the alkali atom and the molecular electrostatic potential minimum is large (>3 Å) for the experimentally realized electrides and ”relatively small” (0.14.) However, the excess electron is also bound quite tightly to the borane, giving a large VDE and a small coherence length (off-diagonal delocalization length) dmax. We conclude that, like antibonding orbitals (Appendix A) and bulk metals (Appendix B), the tightly bound

excess electron of Li(B10H14) combines diagonal delocalization and off-diagonal localization. 3.3. Multiple Formula Units. Calculations on bulk electrides find that the electrons are confined to voids between adjacent formula units.30,39 Figure 4 explores this effect by showing the EDR in the six formula units of the Li+(Cryptand2.1.1)e− unit cell. While such cluster model calculations are not a perfect representation of crystalline electrides, they can provide a reasonable picture of trends inside the central void spaces.45 The left panel of Figure 4 shows the three highestoccupied orbitals containing the six most weakly bound electrons. The right panel shows EDR(r;⃗ dmax = 16.4 bohr). As in Figure 2, the EDR highlights the same regions of space as the frontier orbitals, showing the interstitial electrons. The computed dmax = 16.4 bohr is somewhat larger than the dmax 14.2 bohr evaluated for a single Li+(Cryptand-2.1.1)e− formula F


Article


Figure 5. Isosurfaces of the Li−C4H4N2···Na2 EDR(r;⃗ dmax) = 0.75. RHF (left, dmax = 8.15 bohr), UHF ↑- and ↓-spin (center, dmax = 6.74 bohr), and CAS(2,2) (right, dmax = 7.24 bohr) wave functions.

Figure 6. dmax (left) and Ebond (right) as a function of He−He separation for the ”helium compression chamber” model of one and two ISQ in highpressure electrides. Horizontal line denotes the limit of one isolated He atom. dmax from singlet and triplet UHF calculations on two ISQ are nearly identical to dmax of a single ISQ.

tional studies of electrides,32−35,45 treated the simple case of single formula units with only one weakly bound electron. However, many electride properties depend on interactions of electride electrons. Nondynamical correlation effects can be essential in these interactions.21 A recent review highlighted ”A new world of electrons moving off their atoms and new modes of correlation” in high-pressure electrides.18 Our earlier study46 of stretched singlet H2, a textbook example of nondynamical correlation,57,88−93 illustrated the EDR’s utility for analyzing such systems. Here we present an initial exploration of nondynamical correlation in electrides. We focus on the singlet-coupled electron pair in a single formula unit of the proposed76 electride Li−C4H4N2···Na2. Our calculations suggest that nondynamical correlation plays some role in this system. The brokensymmetry unrestricted Hartree−Fock (UHF) singlet and the symmetric CAS(2,2) singlet states are a modest 11 kcal/mol more stable than the RHF singlet state. Test CAS(4,4) calculations lower the total energy by 0. Supporting Information Figure SI-4 illustrates that a single value of dmax does not completely characterize this system’s B3LYP ΔEDR. 3.5. Delocalization and Nondynamical Correlation. Most of the above calculations, like many previous computaG


Article


Figure 7. dmax (left) and E (right) as a function of H−H separation for H2. Horizontal line denotes the limit of two isolated H atoms.

Figure 8. EDR(r;⃗ dmax) = 0.8 isosurfaces for the He compression chamber with zero (top left), one (top right), and two singlet-coupled (RHF and UHF, bottom left and right) ISQ. He atoms and vacancy sites are shown as spheres. Results at He−He separation 2.0 Å.

CAS(2,2) state is spin-unpolarized and lacks a well-defined ”HOMO”. CAS(2,2) calculations give two natural orbitals with significant noninteger occupancies 1.70 and 0.30. This result motivates further application of the EDR to electron−electron coupling in electrides.

CAS(2,2) calculations localize the electrons and give overall smaller values for EDR(r;⃗ dmax) due to the fractional occupancy effects discussed in ref 47. We emphasize that the spin density analysis in Figure 1 and the HOMO analysis in Supporting Information Figure SI-3 are unavailable for this system, as the H


Article

Journal of Chemical Theory and Computation 3.6. Interstitial Quasi-Atoms in High-Pressure Electrides. The above results show how the EDR illustrates the interplay of electrides’ structure, bonding, correlation, and delocalization. We conclude by exploring high-pressure electrides, which have seen enormous recent computational16−24 and experimental25−28 interest. Figure 6 illustrates the interplay between pressure and delocalization for a ”helium compression chamber” model of interstitial quasi-atoms (ISQ).23,24 The figure shows how the characteristic delocalization length dmax of a single ISQ and two adjacent singlet- and triplet-coupled ISQ vary as the surrounding lattice is compressed. Test calculations on a neutral defect-free He cluster are included for comparison. Symmetry-broken UHF calculations on the singlet-coupled ISQ pair provide a preliminary exploration of nondynamical correlation. The figure also includes the ISQ pairs’ computed bonding energies (eqs 7−9). Figure 6 shows that the EDR provides a rich representation of the interplay of delocalization, confinement, chemical bonding, and electron correlation in high-pressure electrides. As expected, the defect-free He cluster has dmax equal to isolated He atom at large lattice parameter. Compressing the defect-free cluster localizes the electrons due to increasing Pauli repulsion, decreasing dmax. The single ISQ is relatively delocalized, giving dmax larger than the defect-free cluster at all tested lattice parameters. Expanding the lattice increases the single ISQ’s dmax to very large values,95 while compressing the lattice reduces the cavity size and decreases dmax. A notable result in Figure 6 is that the EDR highlights analogies between the chemical-bonding-type interactions of two singlet-coupled ISQs24 and the classic problem of H2 dissociation.57,88−93 We illustrate this by comparing Figure 6 to the ground-state singlet and triplet H2 dmax and binding energy curves in Figure 7. Both systems’ RHF wave functions are excessively delocalized at intermediate separation. For example, the RHF singlet ISQ pair’s dmax is 4.07 bohr at lattice parameter 2 Å, significantly larger than the single ISQ’s dmax 3.541 bohr; while the RHF singlet H2 molecule’s dmax is >4 bohr at H−H separation 3.0 Å, significantly larger than a single H atom’s dmax 2.538 bohr. For both systems, symmetry-broken singlet UHF calculations at intermediate separation lower the energy and dramatically reduce dmax. The triplet states are even more localized, consistent with their additional Pauli repulsion. For example, the UHF singlet and triplet ISQ pair’s dmax are 3.564 and 3.545 bohr at lattice parameter 2 Å, comparable to the single ISQ’s dmax 3.541 bohr; while the UHF singlet and triplet H2 molecule’s dmax are 2.543 and 2.528 bohr at H−H separation 3.0 Å, respectively slightly larger and slightly smaller than the H atom dmax 2.538 bohr. Correlated calculations on singlet H2 yield a similar localization, with dmax 2.615 bohr at H−H separation 3.0 Å. (This difference between UHF and full configuration interaction dmax again highlights the UHF overlocalization mentioned above.94) Another notable result in Figure 7 is that dmax of equilibrium H2 is smaller than that of isolated H atom. This localization upon bonding is consistent with the virial theorem as discussed in ref 57. The large RHF dmax in Figure 6 suggest that our model of two adjacent ISQ is beyond the crossover between RHF singlet and UHF triplet states. The computed energies support this. Like H2 beyond this crossover, Figure 6 shows that the RHF interaction energy of the singlet ISQ pair is more positive (less favorable) than the corresponding UHF singlet or triplet energies. Unlike ref 24, we find that the ISQ pair has positive (thermodynamically unstable) bonding energies at all tested

lattice parameters. This is consistent with differences in the model systems. The uniform background potential used in ref 24 ensures that the 108-atom supercell containing two ”extra” electrons is electrically neutral. In contrast, our model of two adjacent ISQ places two ”extra” electrons into neutral He133, such that the energy difference in eq 7 is dominated by the additional electron−electron repulsion in He2− 133. Supporting Information Figure SI-5 confirms that a single electron bound by two adjacent vacancies (the quasi-atom equivalent of H2+) is thermodynamically stable relative to a single electron in one vacancy. We conclude by reconsidering the real-space distribution of electrons in high-pressure electrides. Figure 8 illustrates EDR(r;⃗ dmax) for the helium compression chamber at lattice parameter 2 Å. Results for the defect-free He crystal chamber (top left, dmax = 1.86 bohr) show that this length scale is characteristic of the valence electrons, such that EDR(r;⃗ dmax) highlights every He atom in the lattice. In contrast, results for the single ISQ (top right, dmax = 3.54 bohr) show that this length scale is characteristic of the defect, such that EDR(r;⃗ dmax) highlights only the ISQ electron. RHF calculations on the two singletcoupled ISQ (bottom left, dmax = 4.07 bohr) have EDR(r;⃗ dmax) highlight the bonding region between the two ISQ, while UHF calculations (bottom right, dmax = 3.56 bohr) have the ↑- and ↓ -spin EDR(r;⃗ dmax) separately highlight the two ISQ centers, illustrating the symmetry breaking and ”relocalization” of electrons. Supporting Information Figure SI-6 shows the corresponding HOMO of the one- and two-ISQ systems, confirming that the results are qualitatively reasonable.

4. CONCLUSIONS The results presented here illustrate that the electron delocalization range function effectively quantifies the delocalization of confined electrons in electrides. The EDR provides both a global measure dmax of the most weakly bound electrons’ average delocalization and a real-space picture EDR(r;⃗ dmax) of these electrons. Our results confirm ref 47: dmax follow chemically reasonable trends, and EDR(r;⃗ dmax) maps out the same region of space as the highest occupied orbitals and spin density. The large dmax computed for experimentally realized electrides suggests that dmax may provide a diagnostic for whether a given formula unit will form a solid electride at ambient pressure. Calculations at various levels of theory illustrate that the EDR can visualize and quantify the interplay of delocalization and correlation critical to electride behavior. Finally, the EDR illustrates analogies23,24 between covalent bonding of atoms and bonding-type interactions between interstital quasi-atoms in high-pressure electrides. Overall, these results motivate adding the EDR to the ”toolbox” of theoretical methods applied to electrides.

■

APPENDIX

A. Diagonal and Off-Diagonal Delocalization in H2+ and H2

The distinction between diagonal delocalization of spin or charge and coherent (off-diagonal) delocalization of the oneparticle density matrix may be illustrated by comparing the H2+ bonding σg and antibonding σu states. We first consider a minimal basis of 1s Gaussian-type functions centered on atoms A and B: χi ( r ⃗) = I

⎛ 2α ⎞3/4 ⎜ ⎟ exp( −α| r ⃗ − R⃗ i|2 ) ⎝π ⎠

, i = {A , B}

(10)


Article


Figure 9. Electron delocalization in the bonding and antibonding states of H+2 , plotted along the bond axis. (Top left) electron density ρ(r)⃗ . (Top ⃗ with point r⃗ fixed to the left H atom (vertical line). (Bottom left) EDR(r⃗;d) at two right) density matrix γ(r⃗,r′)⃗ and EDR test function gd(r⃗−r′), representative length scales d. (Bottom right) System-averaged delocalization ⟨EDR(d)⟩.

σg( r ⃗) = (χA ( r ⃗) + χB ( r ⃗))(2 + 2S)−1

(11)

σu( r ⃗) = (χA ( r ⃗) − χB ( r ⃗))(2 − 2S)−1

(12)

(rA⃗ ;dbond), but has a negative overlap with atom B for the antibonding orbital, leading to a small EDR(rA⃗ ,dbond). The bottom left panel of Figure 9 shows the bonding and antibonding states’ EDR(r;⃗ d) for points r ⃗ along the bond axis, at two representative length scales d. The bonding and antibonding EDR(r;⃗ d) are very similar at small d, where the test function generally probes one atom at a time. However, for d on the order of a covalent bond length, EDR(r;⃗ d) is substantially larger for the off-diagonally delocalized bonding state. The bottom right panel of Figure 9 shows that the system average of eq 4 quantifies the off-diagonal delocalization of the bonding state. Both states’ dmax are marked with vertical lines. The bonding state peaks at dmax larger than the antibonding state. The results in Figure 9 are not an artifact of the basis set or bond length. To illustrate, we consider nearly exact HartreeFock/aug-cc-pVQZ calculations on H+2 at the ground-state optimized bond length 1.057 Å. These calculations give dmax 2.229 bohr for the bonding σg state and 1.602 bohr for the antibonding σu state. We conclude that the EDR quantifies the off-diagonal delocalization that differentiates bonding interactions from antibonding interactions. In chemists’ terms, dmax differentiates the single large lobe of a bonding orbital from the two smaller lobes of an antibonding orbital.

We choose α = 2 bohr−2, RAB = 1.5 bohr. Minimal basis calculations use a Mathematica worksheet provided as Supporting Information. Reference 47 presents a similar analysis of noninteracting Fermions in a 1-D box. The top left panel of Figure 9 shows the bonding and antibonding states’ ρ(r)⃗ for points r ⃗ along the bond axis. Both states exhibit diagonal delocalization of the electron across the two atoms. Each atom has atomic charge +1/2. The top right panel of Figure 9 shows the one-particle density matrix γ(rA⃗ ,r′)⃗ about a point rA⃗ in the atomic basin65 of atom A, for points r′ ⃗ along the bond axis. rA⃗ is marked with a vertical line. While the bonding state’s density matrix is positive for points r′ ⃗ on atom B, the antibonding state’s density matrix is negative on atom B. The bonding orbital exhibits off-diagonal delocalization yielding a positive H−H bond order, while the antibonding orbital exhibits off-diagonal localization yielding a negative H−H bond order. The top right panel of Figure 9 also shows the EDR test function gd(rA⃗ ,r′)⃗ of eq 3 at a distance dbond = 1.3 bohr characteristic of covalent bond lengths, used to evaluate EDR(rA⃗ ,d). The test function has a positive overlap with atom B for the bonding orbital, leading to a large EDRJ


Article

Journal of Chemical Theory and Computation The delocalization index δ(A,B) between atoms A and B is obtained by integrating the exchange-correlation contributions to the diagonal two-particle reduced density matrix ⃗ 64−67 2Γ(r,⃗ r′)⃗ − ρ(r)⃗ ρ(r′)⃗ over atomic domains r∈ ⃗ A, r′∈B. For single-determinant calculations, this yields integration over 2 ⃗ |γ(r∈ . The EDR complements the delocalization ⃗ A,r′∈B)| index in two ways. First, δ(A,B) quantifies delocalization between particular atoms A and B, while EDR(r;⃗ d) quantifies delocalization over a particular distance d. Second, the EDR’s use of γ in place of γ2 distinguishes bonding from antibonding interactions in single-determinant wavefunctions. To illustrate this, Table 3 shows our predictions for the minimal-basis H2

Substituting γUEG into eq 2 gives that the EDR is constant at every point in space ⎛ 9 ⎞1/4 ⎛ erf(d )̃ 2 ⎞ EDR UEG(d) = ⎜ − 2 exp( −d ̃ )⎟ ⎜ π 2⎟ d̃ ⎝ 2πd ̃ ⎠ ⎝ ⎠ d̃ =

(14)

This function peaks at numerical value d̃ = 1.2357..., or d = −1/3 2.4712 k−1 . F = 0.634ρ It is not feasible to directly evaluate dmax from the difference between N and (N-1)-electron uniform electron gases, as these systems are infinite. However, we can consider the delocalization of the most weakly bound frontier electrons whose Fermi vector is between kF and kF(1−δ):

Table 3. H2 HF/STO-3G dmax (bohr), Kinetic Energy Expectation Value ⟨T⟩ (Hartree), and Delocalization Index δ(A,B) (unitless) between Atoms A and Ba system

dmax

⟨T⟩

δ(A,B)

singlet ground-state H2 singly-excited singlet state H2 singly-excited triplet state H2 doubly-excited singlet state H2

2.24 1.72 1.72 1.18

1.20 2.14 2.14 3.05

1 1 1−4(Sab(A))2 1

→ ⎯ γfrontier( r ⃗ , r′) = (2π )−3

(15)

(Note that γfrontier is normalized to ρδ(δ − 3δ+3), not to ρ.) Substituting γfrontier into eq 2 gives a result similar to eq 14 2

⎛ 9 ⎞1/4 ̃ EDR frontier(d) = ⎜ ⎟ [ π (erf(d )̃ + erf((δ − 1)δ )) ⎝ 2πd 2̃ ⎠ ̃ 2(exp(−d 2̃ ) + (δ − 1)exp(− (δ − 1)2 d 2̃ ))] /d −

(16)

As δ → 0, EDRfrontier(d) peaks at a somewhat smaller numerical value d̃ = 0.866.... This is consistent with the fact that the frontier electrons with k → kF have the highest kinetic energy and the smallest off-diagonal delocalization. We conclude, based on both the full density matrix and the density matrix of the most weakly bound electrons, that metals’ coherence lengths (off-diagonal delocalization lengths) dmax scale as the Fermi length k−1 F . This length scale is typically on the order of a covalent bond length. For example, alkali metals have Fermi vectors on the order of 108 cm−1, giving Fermi 97 lengths k−1 F on the order of one Angstrom. In chemists’ terms, metallic bonding does not imply that distant atoms in a metal have large bond orders.

B. Bulk Metals

C. Evaluating the EDR from MP2 Density Matrices

Bulk metals combine diagonal delocalization of charge or spin with off-diagonal localization, yielding characteristic delocalization lengths on the order of the Fermi length k−1 F . This Appendix proves this for the uniform electron gas with constant ↑-spin and ↓-spin electron density ρ(r)⃗ = k3F/(6π2) and noninteracting (exact Kohn−Sham) spin density matrix kF

→ ⎯

⎛ sin(k s) − sin((1 − δ)k s) F F = 3ρ⎜ 3 (kFs) ⎝ cos(kFs) − (1 − δ) cos((1 − δ)kFs) ⎞ ⎟ − (kFs)2 ⎠

states treated using single-determinant wavefunctions in ref 66. Our calculations use H−H bond length 0.74 Å. Unlike the delocalization index, dmax distinguishes the ground, singly, and doubly occupied singlet states of single-determinant wavefunctions. Excited states have higher electronic kinetic energies and smaller characteristic delocalization lengths. We note for completeness that the sharing index proposed by Fulton predicts that ground and singly excited states of H2 both have bond index 1, while the bond index of Cioslowski and Mixton is 1 for the ground state and 1/2 for the first excited state.68,69 We also note that both the delocalization index and the EDR capture the RHF approximation’s over-delocalization of the singlet ground state of stretched H2. Figure 7 shows that dmax of stretched RHF H2 is unrealistically large at long bond lengths (e.g., dmax > 4 bohr at H−H separation 3.5 Å), while dmax of the full configuration interaction wavefunction is comparable to that of H atom. Similarly, ref 96 shows that the delocalization index of stretched singlet RHF H2 is identically 1 at all bond lengths, while the delocalization index from full configuration interaction calculations smoothly converges to zero as the bond is stretched.

∫

kF

∫k (1−δ) d3k ⃗ exp(−ik ⃗· r ⃗)exp(ik ⃗·r′) F

a Delocalization indices from ref 66, other values are our calculations at H−H bond length 0.74 Å. Sab(A) is the overlap of bonding and antibonding orbitals integrated over the domain of atom A.

→ ⎯ γUEG( r ⃗ , r′) = (2π )−3

kFd 2

It is worth noting an issue that can arise in evaluating the EDR from approximate MP2 response density matrices. Response density matrices from nonvariational methods such as MP2 are not guaranteed to have eigenvalues between zero and one.98 These negative natural orbital occupation numbers can produce unphysical negative electron densities. For example, Figure 2 of ref 98 illustrates that the MP2/6-31G(d) response density matrix of stretched N2 yields a region of negative electron density between the atoms. The EDR in these unphysical regions is imaginary due to the ρ−1/2(r)⃗ prefactor in eq 3. While our calculations only evaluate the EDR at points where ρ(r)⃗ > (a small positive threshold), we find that the EDR can become unphysical at points close to the unphysical negative density regions. In particular, we find that the MP2/6-31++G(d,p)

→ ⎯ d3k ⃗ exp(−ik ⃗· r ⃗)exp(ik ⃗·r′)

⎛ sin(k s) cos(kFs) ⎞ F ⎟ = 3ρ⎜ − 3 (kFs)2 ⎠ ⎝ (kFs) → ⎯ s ≡ | r ⃗ − r′| (13) K


Article

Journal of Chemical Theory and Computation density matrix of a single formula unit of K+(Cryptand-2.2.2)e− has minimum eigenvalue −0.092, minimum ρ(r)⃗ = −1.96 × 10−6 electrons/bohr3, and an imaginary EDR in regions near the solvated electron. We also note that our computational setup cannot evaluate the MP2/6-31++G(d,p) response density matrix of single formula units of Na+(tripip-aza-2.2.2)e− and Cs+(15C5)2e−.

■

(21) Gatti, M.; Tokatly, I. V.; Rubio, A. Phys. Rev. Lett. 2010, 104, 216404. (22) Zhu, Q.; Oganov, A. R.; Lyakhov, A. O. Phys. Chem. Chem. Phys. 2013, 15, 7696−7700. (23) Miao, M.-S.; Hoffmann, R. Acc. Chem. Res. 2014, 47, 1311− 1317. (24) Miao, M.-S.; Hoffmann, R. J. Am. Chem. Soc. 2015, 137, 3631− 3637. (25) Matsuoka, T.; Shimizu, K. Nature 2009, 458, 186−189. (26) Ma, Y.; Eremets, M.; Oganov, A. R.; Xie, Y.; Trojan, I.; Medvedev, S.; Lyakhov, A. O.; Valle, M.; Prakapenka, V. Nature 2009, 458, 182−186. (27) Lazicki, A.; Goncharov, A. F.; Struzhkin, V. V.; Cohen, R. E.; Liu, Z.; Gregoryanz, E.; Guillaume, C.; Mao, H.-K.; Hemley, R. J. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 6525−6528. (28) Marqués, M.; Santoro, M.; Guillaume, C. L.; Gorelli, F. A.; Contreras-García, J.; Howie, R. T.; Goncharov, A. F.; Gregoryanz, E. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 184106. (29) Singh, D. J.; Krakauer, H.; Haas, C.; Pickett, W. E. Nature 1993, 365, 39−42. (30) Dye, J. L.; Wagner, M. J.; Overny, G.; Huang, R. H.; Nagy, T. F.; Tománek, D. J. Am. Chem. Soc. 1996, 118, 7329−7336. (31) Chen, W.; Li, Z.-R.; Wu, D.; Li, R.-Y.; Sun, C.-C. J. Phys. Chem. B 2005, 109, 601−608. (32) Chen, W.; Li, Z.-R.; Wu, D.; Li, Y.; Sun, C.-C.; Gu, F. L. J. Am. Chem. Soc. 2005, 127, 10977−10981. (33) Muhammad, S.; Xu, H.; Liao, Y.; Kan, Y.; Su, Z. J. Am. Chem. Soc. 2009, 131, 11833−11840. (34) Li, Z.-J.; Wang, F.-F.; Li, Z.-R.; Xu, H.-L.; Huang, X.-R.; Wu, D.; Chen, W.; Yu, G.-T.; Gu, F. L.; Aoki, Y. Phys. Chem. Chem. Phys. 2009, 11, 402−408. (35) Garcia-Borrás, M.; Solà, M.; Luis, J. M.; Kirtman, B. J. Chem. Theory Comput. 2012, 8, 2688−2697. (36) Rousseau, B.; Ashcroft, N. W. Phys. Rev. Lett. 2008, 101, 046407. (37) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, P.; Cohen, A. J.; Yang, W. J. Am. Chem. Soc. 2010, 132, 6498. (38) Bader, R. F. W. Atoms In Molecules - A Quantum Theory; Oxford University Press: Oxford, 1990; pp 1−20. (39) Dale, S. G.; Otero-de-la Roza, A.; Johnson, E. R. Phys. Chem. Chem. Phys. 2014, 16, 14584−14593. (40) Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990, 92, 5397. (41) Savin, A. J. Mol. Struct.: THEOCHEM 2005, 727, 127. (42) Baker, C. F.; Barker, M. G.; Blake, A. J. Acta Crystallogr., Sect. E: Struct. Rep. Online 2001, E57, i6−i7. (43) Lee, K.; Kim, S. W.; Toda, Y.; Matsuishi, S.; Hosono, H. Nature 2013, 494, 336−341. (44) Postils, V.; Garcia-Borràs, M.; Solà, M.; Luis, J. M.; Matito, E. Chem. Commun. 2015, 51, 4865. (45) Kumar, A.; Gadre, S. R. Phys. Chem. Chem. Phys. 2015, 17, 15030−15035. (46) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Phys. 2014, 141, 144104. (47) Janesko, B. G.; Scalmani, G.; Frisch, M. J. Phys. Chem. Chem. Phys. 2015, 17, 18305−18317. (48) All formulas use Hartree atomic units. Unless noted otherwise, electron spin and time dependence are suppressed throughout for conciseness. All calculations treat approximate stationary state solutions of Schrödinger’s equation. (49) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257−2261. (50) Hariharan, P.; Pople, J. Theor. Chem. Acc. 1973, 28, 213−222. (51) Johnson, E. R.; Otero-de-la Roza, A.; Dale, S. G. J. Chem. Phys. 2013, 139, 184116. (52) Williams, C. F.; Herbert, J. M. J. Phys. Chem. A 2008, 112, 6171−6178. (53) Kevan, L. Acc. Chem. Res. 1981, 14, 138. (54) Tretiak, S.; Mukamel, S. Chem. Rev. 2002, 102, 3171. (55) Zachariasen, W. H. J. Am. Chem. Soc. 1940, 62, 1011−1013. (56) Bolton, J. R.; Carrington, A. Mol. Phys. 1961, 4, 497−503.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00993. Figures SI-1−SI-6 (PDF) All calculated geometries, total energies, and tables of ⟨EDR(d)⟩ at representative d values (TXT) Mathematica worksheet for the H2+ calculations in Appendix A and the UEG calculations in Appendix B (TXT)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS B.G.J. acknowledges support from the National Science Fundation (DMR-1505343). REFERENCES

(1) Jacobson, L. D.; Herbert, J. M. Int. Rev. Phys. Chem. 2011, 30, 1− 48. (2) Turi, L.; Rossky, P. Chem. Rev. 2012, 112, 5641. (3) Casey, J. R.; Kahros, A.; Schwartz, B. J. J. Phys. Chem. B 2013, 117, 14173−14182. (4) Schwartz, B. J.; Rossky, P. J. J. Chem. Phys. 1994, 101, 6902− 6916. (5) Zurek, E.; Edwards, P. P.; Hoffmann, R. Angew. Chem., Int. Ed. 2009, 48, 8191−8232. (6) Kölmel, C.; Ewig, C. S. J. Phys. Chem. B 2001, 105, 8538−8543. (7) Gudmundsdottir, H.; Zhang, Y.; Weber, P. M.; Jonsson, H. J. Chem. Phys. 2014, 141, 234308. (8) Dye, J. M. Acc. Chem. Res. 2009, 42, 1564. (9) Moeggenborg, K. J.; Papaioannou, J.; Dye, J. L. Chem. Mater. 1991, 3, 514−520. (10) Kitano, M.; Inoue, Y.; Yamazaki, Y.; Hayashi, F.; Kanbara, S.; Matsuishi, S.; Yokoyama, T.; Kim, S.-W.; Hara, M.; Hosono, H. Nat. Chem. 2012, 4, 934−940. (11) Kuganathan, N.; Hosono, H.; Shluger, A. L.; Sushko, P. V. J. Am. Chem. Soc. 2014, 136, 2216−2219. (12) Kitano, M.; Kanbara, S.; Inoue, Y.; Kuganathan, N.; Sushko, P. V.; Yokoyama, T.; Hara, M.; Hosono, H. Nat. Commun. 2015, 6, 6731. (13) Huang, R. H.; Dye, J. L. Chem. Phys. Lett. 1990, 166, 133−136. (14) Matsuishi, S.; Toda, Y.; Miyakawa, M.; Hayashi, K.; Kamiya, T.; Hirano, M.; Tanaka, I.; Hosono, H. Science 2003, 301, 626. (15) Redko, M. Y.; Jackson, J. E.; Huang, R. H.; Dye, J. L. J. Am. Chem. Soc. 2005, 127, 12416−124212. (16) Neaton, J. B.; Ashcroft, N. W. Nature 1999, 400, 141−144. (17) Christensen, N. E.; Novikov, D. L. Solid State Commun. 2001, 119, 477−490. (18) Grochala, W.; Hoffmann, R.; Feng, J.; Ashcroft, N. W. Angew. Chem., Int. Ed. 2007, 46, 3620−3642. (19) Pickard, C. J.; Needs, R. J. Phys. Rev. Lett. 2009, 102, 146401. (20) Yao, Y.; Tse, J. S.; Klug, D. D. Phys. Rev. Lett. 2009, 102, 115503. L


Article

Journal of Chemical Theory and Computation (57) Ruedenberg, K. Rev. Mod. Phys. 1962, 34, 326−376. (58) Hückel, E. Eur. Phys. J. A 1931, 70, 204. (59) Coulson, C. A. Proc. R. Soc. London, Ser. A 1939, 169, 413−428. (60) McWeeny, R. J. Chem. Phys. 1951, 19, 1614−1615. (61) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (62) Mulliken, R. S. J. Chem. Phys. 1955, 23, 2343. (63) Wiberg, K. A. Tetrahedron 1968, 24, 1083. (64) Bader, R. F. W.; Streitwieser, A.; Neuhaus, A.; Laidig, K. E.; Speers, P. J. Am. Chem. Soc. 1996, 118, 4959−4965. (65) Bader, R. W. F.; Heard, G. L. J. Chem. Phys. 1999, 111, 8789. (66) Fradera, X.; Solà, M. J. Comput. Chem. 2002, 23, 1347−1356. (67) Poater, J.; Duran, M.; Solá, M.; Silvi, B. Chem. Rev. 2005, 105, 3911. (68) Cioslowski, J.; Mixon, S. T. J. Am. Chem. Soc. 1991, 113, 4142. (69) Fulton, R. L.; Mixton, S. T. J. Phys. Chem. 1993, 97, 7530−7534. (70) Schmider, H. L. J. Chem. Phys. 1996, 105, 11134. (71) Schmider, H. L.; Becke, A. D. J. Mol. Struct.: THEOCHEM 2000, 527, 51. (72) Kohout, M. Int. J. Quantum Chem. 2004, 97, 651. (73) Dawes, S. B.; Eglin, J. L.; Moeggenborg, K. J.; Kim, J.; Dye, J. L. J. Am. Chem. Soc. 1991, 113, 1605−1609. (74) Ward, D. L.; Huang, R. H.; Dye, J. L. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1988, C44, 1374−1376. (75) Huang, R. H.; Wagner, M. J.; Gilbert, D. J.; Reidy-Cedergren, K. A.; Ward, D. L.; Farber, M. K.; Dye, J. L. J. Am. Chem. Soc. 1997, 119, 3765−3772. (76) Ma, F.; Li, Z.-R.; Xu, H.-L.; Li, Z.-J.; Li, Z.-S.; Aoki, Y.; Gu, F. L. J. Phys. Chem. A 2008, 112, 11462−11467. (77) 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.; Fox, D. J. Gaussian Development Version, Revision I.02+; Gaussian, Inc.: Wallingford, CT, 2014. (78) Janesko, B. G.; Scalmani, G.; Frisch, M. J. J. Chem. Phys. 2014, 141, 034103. (79) Peterson, K.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. J. Chem. Phys. 2003, 119, 11113. (80) Hehre, W. J.; Stewart, R. F.; Pople, J. A. J. Chem. Phys. 1969, 51, 2657. (81) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (82) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (83) Handy, N. C.; Schaefer, H. F., III J. Chem. Phys. 1984, 81, 5031. (84) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. J.; Frisch, M. J. J. Phys. Chem. 1992, 96, 671−679. (85) Herbert, J. M.; Head-Gordon, M. J. Phys. Chem. A 2005, 109, 5217−5229. (86) Herbert, J. M.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2006, 8, 68−78. (87) Khatkale, M. S.; Devlin, J. P. J. Chem. Phys. 1979, 70, 1851. (88) Gunnarsson, O.; Lundqvist, B. I. Phys. Rev. B 1976, 13, 4274. (89) Cremer, D. Mol. Phys. 2001, 99, 1899. (90) Becke, A. D. J. Chem. Phys. 2003, 119, 2972. (91) Ponec, R.; Cooper, D. L. J. Phys. Chem. A 2007, 111, 11294. (92) Perdew, J. P.; Ruzsinszky, A.; Constantin, L. A.; Sun, J.; Csonka, G. I. J. Chem. Theory Comput. 2009, 5, 902.

(93) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. J. Chem. Phys. 2008, 129, 121104. (94) Hollett, J. W.; McKemmish, L. K.; Gill, P. M. W. J. Chem. Phys. 2011, 134, 224103. (95) Our use of the finite aug-cc-pVQZ basis set for ISQs ensures that a single ISQ’s dmax goes to 11.4 bohr in the limit of large lattice parameter. Calculations with a complete basis set would have dmax → ∞. (96) Matito, E.; Duran, M.; Solá, M. J. Chem. Educ. 2006, 83, 1243− 1248. (97) Ashcroft, N. R.; Mermin, N. D. Solid State Physics; Harcourt College Publishers: Fort Worth, 1976; p 38. (98) Gordon, M. S.; Schmidt, M. W.; Chaban, G. M.; Glaesemann, K. R.; Stevens, W. J.; Gonzalez, C. J. Chem. Phys. 1999, 110, 4199−4207.

M


Quantifying Electron Delocalization in Electrides - ACS Publications

Recommend Documents