Accurate Dissociation of Chemical Bonds Using DFT-in-DFT

Chemistry Department, University of North Dakota, Grand Forks, North Dakota 58202, United States. J. Phys. Chem. A , 2017, 121 (1), pp 256–264. DOI:...
3 downloads 12 Views 516KB Size
Subscriber access provided by University of Otago Library

Article

Accurate Dissociation of Chemical Bonds Using DFT-inDFT Embedding Theory with External Orbital Orthogonality Patrick K. Tamukong, Yuriy G. Khait, and Mark R. Hoffmann J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b09909 • Publication Date (Web): 09 Dec 2016 Downloaded from http://pubs.acs.org on December 10, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Accurate Dissociation of Chemical Bonds Using DFT-in-DFT Embedding Theory with External Orbital Orthogonality Patrick K. Tamukong, Yuriy G. Khait, and Mark R. Hoffmann* Chemistry Department, University of North Dakota, Grand Forks, North Dakota 58202

Abstract Our recent density functional theory (DFT)-in-DFT embedding protocol, which enforces inter-subsystem (or external orbital) orthogonality, is used for the first time to investigate covalent bond dissociation and is shown to do so accurately. Full potential energy curves for the dissociation of an H-O bond in H2O and the C-C bond in H 3 C − CH 3 have been constructed using the new embedding method, as have the challenging ionic bonds in LiH and LiF, and were found to match the reference Kohn-Sham (KS)-DFT curves to at least one part in 106. The added constraint of external orbital orthogonality allows for the formulation of an embedding protocol that neither relies on approximate kinetic energy functionals for the evaluation of the so-called nonadditive kinetic potential, nor introduces compensatory potentials, nor requires a total system calculation at any stage. The present work extends the demonstrated applicability of the external orthogonality variant of embedding theory by more than a factor of two to the interaction strength range of strong single bonds. In particular, it is demonstrated that homolytic cleavage of both covalent and ionic bonds into radicals can be accomplished.

*

Author for correspondence. E-Mail: [email protected], Phone: (+1) 701-777-2742, FAX: (+1) 701-7772331

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 34

I. Introduction DFT-in-DFT embedding theories seek to describe relatively large systems by partitioning such systems into fragments that are both described using density functional theory (DFT) (cf. e.g. refs. 1–7 and recent reviews 8 and 9). Conventional DFT scales cubically with system size, making it incapable of simulating some of the larger systems investigated in laboratory experiments, especially those of biological or materials relevance. Moreover, so-called environmental effects (e.g., solvent effects, inter-chromophore effects) are included in many computer simulations only implicitly due to size limitations. Examples of implicit solvent inclusion include the conductor-like screening model (COSMO),10 and polarizable continuum model (PCM) methods, such as the integral equation formalism PCM (IEF-PCM)11 and the continuous surface charge formalism.12 Although such models describe well scenarios in which solute-solvent interactions are mainly electrostatic in nature, they are pressed for use in many situations in which they are not ideal, such as those in which explicit solvent inclusion would be appropriate but would lead to prohibitively high computation costs. Effective fragment potential (EFP) methods13-15 are a significant improvement wherein a system is partitioned into fragments and the fragment of interest described using an ab initio quantum chemistry method, while approximating the rest of the fragments by potentials and assuming non-bonded interactions between surrounding fragments and the fragment of interest. Such potentials are, however, most appropriate when quantum mechanical bonding interactions between fragments are not significant. Since many chemical reactions cannot be well described by neglecting environmental effects (whether inter- or intra-molecular), the development of methods for describing large realistic systems has been a significant subject of theoretical chemistry research. In this regard,

2 ACS Paragon Plus Environment

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

new methods16,17 have been developed for describing relatively large systems, based on the idea of “from fragments to molecule”. Embedding theory similarly uses a fragment-based or divideand-conquer approach to problem solving. Since many properties of interest depend only on local effects, such an approach is not only computationally efficient, but captures the essential physics in a natural way. Specifically, a system can be partitioned into a small fragment of interest (the embedded subsystem) while the effects of the rest of the system (the environment) on the embedded subsystem are considered at the same or another level of approximation. Furthermore, embedding theory allows for multiple disjoint fragments. In DFT-in-DFT, a system is partitioned into an embedded subsystem (hereafter referred to as subsystem A) and an environment subsystem (hereafter referred to as subsystem B), although B could also be, e.g., another molecule or part of a molecule. The fragments are

r

r

assigned integer numbers of electrons, NA and NB, and densities ρ A ( r ) and ρ B (r ) , such that

r r r ρ tot (r ) = ρ A (r ) + ρ B (r ) ,

(r )

(1)

r

where ρ tot r is the total density. Thus partitioned, the environment density, ρ B (r ) , could be

r

determined once and subsequently frozen during optimization of the embedded density, ρ A ( r ) , (referred to as frozen density embedding (FDE)9) or both densities can be optimized selfconsistently in a so-called freeze-and-thaw cyclic procedure. In either case, the coupled equations are generally termed Kohn-Sham equations with constrained electron density (KSCED). If the densities are disjoint, the nonadditive kinetic potential (i.e. vT) is rigorously zero;9,18 since approximation to vT is a significant source of error, embedding methods in which vT=0 is enforced have been pursued.7,19–22 Partitioning the density of a system into a sum of fragment densities is conveniently done, within an orbital approximation, if orbitals of those

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 34

fragments are strictly orthogonal to each other.7,19–22 Not taking into account external orbital orthogonality leads to increased complexity in formulations of DFT-in-DFT. In very recent work, Unsleber et al.23 demonstrate formally that density summability (i.e., eq 1) can be achieved without enforcing external orbital orthogonality. Through analytical and numerical study, they show that the external orbital pair orthogonality is the relevant criterion. But, they also observe that for modest basis sets, including those of most DFT calculations, external orthogonality is essentially enforced by their derived exact conditions on density. Their study shows that the statement of Hoffmann and coworkers19,21,22 on the necessity of external orthogonality needs to be considered as a sufficiency, especially in the case of very large basis sets. Most DFT-in-DFT research efforts (cf. refs 8 and 9) were (and still are) directed at obtaining improved approximations to the nonadditive kinetic potential,24 recasting vT in an exact form,25 or seeking a unique embedding potential6,26,27 common to interacting subsystems. Most of these developments result in methods that described weakly interacting subsystems well, but prove to be less successful for more strongly coupled subsystems. For example, Beyhan et al.28 used several kinetic energy functionals within the frozen density variant of DFT-in-DFT (i.e., frozen density embedding, FDE) to study the weak Ng-AuF bond (where Ng = Ar, Kr, Xe) and concluded that DFT-in-DFT embedding was incapable of describing metal-ligand interactions with covalent character. Besides the possibility of artificial charge leaks,29 the

r

quality or accuracy of FDE depends on the choice of the frozen density, ρ B (r ) .30 Jacob et al. concluded that DFT-in-DFT was inadequate for strongly coupled subsystems, pending the development of more accurate approximations to vT that are applicable in cases of strong interactions.31 Fux et al.,32 as well as Kiewisch et al.,33 also focused on the connection between

4 ACS Paragon Plus Environment

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the quality of FDE and exact exchange-correlation, E xc [ρ], and kinetic energy, Ts [ρ] , functionals. These authors compared electron densities of selected systems obtained from FDE embedding with reference KS-DFT densities and found non-negligible deviations, particularly at the interface between the subsystems. In particular, the study of Fux et al. found FDE to fail for cases involving strong inter-subsystem interactions and connected such inaccuracies with inaccuracies in vT.32 Alternative approaches to DFT-in-DFT embedding methods that seek a unique embedding potential common to both subsystems have yielded promising results on the tested systems.6,26,27 However, they involve initial KS-DFT calculations on the total system; such approaches are intrinsically limited to systems that are amenable to whole system DFT calculations. And since they result in more computational cost in comparison with conventional DFT, the strategy is best suited for wavefunction-in-DFT approaches, in which the preliminary whole system DFT calculation precedes a computationally expensive subsequent step. Similarly, the work of Miller and coworkers that seeks to evaluate vT in a more accurate form involves an initial total system DFT calculation.25 More recent work by Miller and coworkers7,20 introduced a variant of DFT-in-DFT (and wavefunction-in-DFT) which enforces the orthogonality of orbitals between subsystems by use of a level shift projection operator that shifts the energies of the KS orbitals of the complementary subsystem to higher values; of course, an initial total system KS-DFT calculation is performed. Perhaps one of the strongest interactions to have been well described previously using a variant of DFT-in-DFT is the heterolytic cleavage of the highly polarized C-C bond of

H3C − CF3 into H 3C + and CF3- fragments.34 In that study, the authors considered H 3C + and CF3-

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 34

fragments as separate subsystems with 8 and 34 electrons, respectively. Their method was able to overcome large deviations (at least 730 kcal/mol) obtained using conventional DFT-in-DFT. In contrast to other efforts, the method of DFT-in-DFT that this laboratory introduced is agnostic of the total system, as is traditional DFT-in-DFT embedding, but includes orbital orthogonality between subsystems as an added constraint, using a Lagrangian formulation.19,21,22 And since the method does not require a total system calculation at any stage, it is potentially attractive for calculations of large systems such as biological molecules. Enforcing external orbital orthogonality eliminates the nonadditive kinetic potential (i.e. vT(r)=0 everywhere),9,18 and therefore results in a method that does not rely on approximate kinetic energy functionals. This method was recently shown to accurately describe subsystems that were more strongly interacting than conventional DFT-in-DFT could reproduce.21,22 Interaction energies agreeing nearly exactly with reference KS-DFT values were obtained. Moreover, electron densities were also well reproduced by the new embedding method, with only negligible deviations at the interface between subsystems, irrespective of interaction strength. Nonetheless, the question of whether the embedding-with-external-orbital-orthogonality method could address the yet more complicated electron densities encountered when dissociating actual chemical bonds into radical moieties had not been asked. Here, we report full potential energy curves (PECs) for the dissociation of the covalent H-O bond in H 2 O and the C-C bond in

H 3 C − CH 3 as well as the difficult to characterize ionic bonds in LiH and LiF, and compare to those of reference KS-DFT. The rest of this article is organized as follows. In Section II, the external orthogonality variant of DFT-in-DFT is reviewed and other computational details are described. In Section III, results are presented and discussed; conclusions are given in Section IV. 6 ACS Paragon Plus Environment

Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

II. Theory a. DFT-in-DFT with External Orbital Orthogonality Details of the embedding-with-external-orbital-orthogonality DFT-in-DFT method can be found elsewhere,19,21,22 so only the key equations relevant to the present study are reviewed. To enforce orbital orthogonality between subsystems within DFT-in-DFT embedding theory, the following Lagrangian is constructed

[

]

[

] ∑ ∑Θ

Ω ρA , ρB = E v ρA + ρB −

NI

I cc′

φ cI φ cI ′

I = A, B c, c′ = 1



NA NB

∑ ∑α

a =1b =1

ba

φ aA φ Bb −

(2)

NA NB

∑ ∑β

ab

φ Bb φ aA

a =1b =1

where Ev[ρA+ρB] is the energy functional for the total system expressed as a functional of the densities of subsystems; the second term on the right hand side expresses the constraint of orbital orthonormality ( ϕ iI ϕ jI = δij where I = A or B while δij is Kronecker’s delta) within a subsystem, which is convenient although not essential; while the last two terms express the constraint of inter-subsystem or external orbital orthogonality. The coefficients ϴcc’, αba and βab are Lagrange multipliers associated with these constraints. Note that for a canonical set of orbitals, ϕ A

or ϕ B , one has Θ aAa ′ = δ aa ′ ε aA and Θ Bbb′ = δ bb′ ε Bb , respectively. The energy

functional, Ev[ρA+ρB], is partitioned into

E v [ρ A + ρB ] = E v A [ρA ] + E v B [ρB ] + Eint [ρA ,ρB ] ,

(3)

where E v A [ρA ] and E v B [ρB ] are the functionals for subsystems A and B, respectively (with the same definitions as in usual DFT), while E int [ρ A , ρ B ] describes inter-subsystem interactions, and is defined as 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

E int [ρ A , ρ B ] = F[ρ A + ρ B ] − F[ρ A ] − F[ρ B ] + ∫ (v Aρ B + v Bρ A )dr ,

(4) r

where F [ρ ] is the kinetic energy and electron-electron interaction functional, and v A (r ) and r v B (r ) represent the potentials due to nuclei assigned to the individual subsystems.

Minimization of the Lagrangian in eq 2 with respect to the orbitals of the subsystems (and hence, allowable densities), subject to the charge constraints

r r

∫ ρ (r )d r = N I

I

(I = A or B),

leads to the set of coupled Euler-Lagrange equations

[h [h

KS

KS

]

NA

NB

a′

b =1

+ vTA ϕaA = ∑ ϕaA′ ΘaA′a + ∑ ϕbB α ba

+v

B T



B b

NB

=∑ϕ

B b′

Θ

b′

B b′b

NA

(a ∈ [1, N A ])

+ ∑ ϕ aA α ab

(b ∈ [1, N B ])

(5a)

(5b)

a

where hKS is the Kohn-Sham Hamiltonian and v TI (I = A, B) is the nonadditive kinetic energy potential (N.B. v TI = 0 in the case of external orthogonality). Assuming canonical orbital sets

ϕ A and ϕ B leads to

[h

KS

]

NB

+ v TA ϕ aA = ϕ aA ε aA + ∑ ϕ bB α ba

(a ∈ [1, N A ])

(6a)

(b ∈ [1, N B ])

(6b)

b =1

and

[h

KS

]

NA

+ v TB ϕ bB = ϕ bB ε Bb + ∑ ϕ aA β ab a

With the constraint ϕ aA ϕ bB = ϕ bB ϕ aA = 0 , one has

α ab = α ba = ϕ bB h KS + v TA ϕ aA

and β ab = β ba = ϕaA h KS + vTB ϕbB

(7)

Hence, eq 6 may be rewritten as 8 ACS Paragon Plus Environment

Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(I − P )[h B

(I − P )[h A

KS

KS

]

+ v TA ϕ aA = ϕ aA ε aA

(8a)

]

+ v TB ϕ bB = ϕ bB ε Bb

(8b)

where NA

NB

a =1

b =1

P A = ϕ A ϕ A = ∑ ϕ aA ϕ aA and P B = ϕ B ϕ B = ∑ ϕ bB ϕ bB

(9)

are projection operators with respect to subsystems A and B, respectively. Noting that

(I − P )ϕ B

A

(

)

= ϕ A and I − P A ϕ bB = ϕ bB , eq 8 may alternatively be written in hermitian form

as

(I − P )[h

KS

+ v TA I − P B ϕ aA = ϕ aA ε aA ,

(I − P )[h

KS

+ v TB I − P A ϕ bB = ϕ bB ε Bb

B

A

](

)

(10a)

](

)

(10b)

where I represents the identity operator within the total one-electron space. Thus, the embeddingwith-external-orbital-orthogonality equations require the KS orbitals of each subsystem to be

(

)[

](

)

(

)[

](

eigenvectors of the reduced Hamiltonians I − P B h KS + v TA I − P B and I − P A h KS + v TB I − P A

)

(of course, with v TI = 0 ). Additional details on the computational implementation of this method are available in Refs. 21 and 22.

b. Computational Details Calculations of the potential energy curves (PECs) of the dissociation of LiH, LiF, an HO bond in H 2 O , and the C-C bond in H 3 C − CH 3 into radical fragments were performed using locally developed computer programs based on the above described embedding method. These calculations are designated by KSCED(e), KSCED(s), KSCED(e, Ext. Orth.) and KSCED(s, Ext. Orth.), where the letters e and s within parentheses denote that either an extended monomer 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 34

basis21,22 or a supermolecular basis was used to expand the orbitals of the subsystems, respectively; Ext. Orth. denotes that inter-subsystem orbital orthogonality was enforced and its absence that it was not. Note that the first two cases (i.e., KSCED(e), KSCED(s)) represent conventional DFT-in-DFT calculations, i.e. without orthogonality; a Lembarki-Chermette (LC94) kinetic functional35 was used to calculate vT for these comparative calculations. In calculations of the dissociation of LiH, H • and Li • composed the two subsystems; for LiF, Li • and F • ; for H 2 O , H • and HO • were the two subsystems, while for H 3 C − CH 3 , each H 3 C • radical was treated as a subsystem. The unpaired electron on one of the embedded fragments was constrained to have beta spin projection while that on the other was constrained to have alpha spin projection (i.e., Ms = 1/2). The program specifications were set as follows: integration threshold = 5.749 × 10-11; self-consistent field energy convergence criterion = 10-6; and convergence criterion for orbital rotations = 10-8. A Mura-Knowles log3 grid type36 with 96 radial and 302 angular grid points was used. The VWN5,37 PW91,38,39 and B3LYP40-43 functionals with the cc-pVTZ44 basis set were used in the calculations. Specifically, the H2O and LiH calculations utilized the VWN5 and B3LYP functionals, H 3 C − CH 3 used PW91, while LiF used B3LYP. Use of different functionals was intended to investigate the sensitivity of the embedding-with-external-orbital-orthogonality

DFT-in-DFT

method

to

different

DFT

functionals. In calculations of H 2 O and H 3 C − CH 3 , the experimental geometries of the molecules were used as the reference geometries, and the O-H and C-C bonds being dissociated were subsequently varied while all other bonds and angles were kept fixed. Vibrational frequencies were computed using the harmonic oscillator approximation and three quadrature points, i.e.

10 ACS Paragon Plus Environment

Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1 k 1 ω= = 2πc µ 2πc

∂2E

∂R 2

µ

∂2E 1 E ( R + ∆R ) − 2 E ( R ) + E ( R − ∆R )  = 2  2 ∂R ( ∆R ) 

(11a)

(11b)

where ω is the harmonic frequency; k is the force constant and µ the reduced mass; c is the speed of light; E is the total energy; R is the bond length; and a displacement of ∆R=0.01 Å was used.

III. Results and Discussion a. LiF The PECs of the dissociation of LiF, obtained with KSCED(s, Ext. Orth.) and KS-DFT using B3LYP and the cc-pVTZ basis set, are shown in Figure 1. As can be seen, the two curves match each other exactly (i.e., to within 10-5 Hartree) at all geometries and lead essentially to the same spectroscopic constants as shown in Table 1. The Li-F equilibrium bond length of 1.57 Å and binding energy of 5.43 eV obtained by both methods compare favorably with the experimental values of 1.564 Å45 and 5.898±0.35 eV,46 respectively. The obtained frequency of 917.9 cm-1 is also close to the experimental value of 910.57 cm-1.47 As shown in Table 1, the KSCED(s) method (i.e., without orbital orthogonality and using LC94 KE functional) leads to results that are not even in qualitative agreement with reference KS-DFT and experimental data for LiF. The KSCED(s) bond length of 3.30 Å is 1.736 Å longer than the experimental value (i.e., 1.564 Å), whereas the bond energy of 1.18 eV is 4.718 eV less than the experimental value of 5.898 eV. Thus, DFT-in-DFT, without additional modifications, proves utterly incapable of describing bond cleavage in LiF. This is consistent with the

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

observation of Goodpaster et al.34 (and our earlier work21) that even for the weak hydrogen bond in the water dimer (N.B. the experimental value of the hydrogen bond in water is only 0.24 eV,48 quite weak compared to > 5 eV for LiF), deviations of 20% in bond energy were obtained using KSCED(s) (and a Lembarki-Chermette (LC94) kinetic functional35). When using the ThomasFermi kinetic energy functional, yet larger deviations of 40% were observed in the hydrogen bond energy of the water dimer.34 Of course, their exact embedding method7, which requires a total system KS-DFT calculation at the initial stage, reproduced KS-DFT results.

b. LiH The PECs of the dissociation of LiH, obtained with KSCED(s, Ext. Orth.) and KS-DFT using the cc-pVTZ basis set together with the B3LYP or VWN5 functional, are shown in Figure 2. The data describing the curves are in Table 1, as are corresponding experimental results. As can be seen in Figure 2, the KSCED(s, Ext. Orth.) and KS-DFT curves are indistinguishable to within the width of the line for both functionals (B3LYP and VWN5) at all geometries, indicating the insensitivity of the embedding-with-external-orbital-orthogonality DFT-in-DFT method to DFT functional type. It should be noted that DFT-in-DFT, without orbital orthogonality and using LC94 KE functional, has previously been demonstrated (e.g., ref 21) to show dependency on type of DFT functional. Whereas the quality of KS-DFT results generally depends on functional type49 (i.e., higher rungs on the DFT ladder generally yield more accurate results than lower rungs50), DFT-in-DFT (without orbital orthogonality and using LC94 KE functional), lacks predictability in accuracy of results in relation to functional type. For example, for a variety of weakly interacting complexes for which LDA functionals proved adequate in DFT-in-DFT, GGA functionals used with DFT-in-DFT led to worse results whereas KS-DFT

12 ACS Paragon Plus Environment

Page 13 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

with those GGA functionals improved the results in comparison with reference data.51 The accuracy of embedding-with-external-orbital-orthogonality DFT-in-DFT is insensitive to choice of DFT functional tested thus far and should be expected to essentially reproduce reference KSDFT results with any functional. As can be seen in Table 1, the B3LYP results are in better agreement with experimental data than those obtained with VWN5. The bond length of 1.59 Å obtained with B3LYP in KSCED(s, Ext. Orth.) and KS-DFT agrees quite well with the experimental value of 1.596 Å.52 Likewise the B3LYP binding energy of 2.53 eV and harmonic frequencies of 1416.7 cm-1 and 1417.0 cm-1, obtained with KSCED(s, Ext. Orth.) and KS-DFT, respectively, are in good agreement with the experimental values of 2.52 eV53 and 1405.50 cm-1.54 On the other hand, the bond length of 1.61 Å, binding energy of 2.63 eV, and harmonic frequency of 1361.0 cm-1 obtained with both KSCED(s, Ext. Orth.) and KS-DFT calculations with VWN5 deviate from experimental values by 0.014 Å (or 0.88%), 0.11 eV (or 4.37%), and 44.5 cm-1 (or 3.2%). Perhaps because of the much lower electron density of H vis-à-vis F, the KSCED(s) method with B3LYP performs much better for LiH than for LiF. In particular, the KSCED(s)/B3LYP harmonic frequency of 1482.1 cm-1 is larger than the experimental value by only 76.6 cm-1, and the bond length and energy are 0.166 Å and 0.87 eV less than reference experimental values, respectively. This difference between LiH and LiF in terms of KSCED(s)/B3LYP performance may also be understood in terms of their much different interaction strengths (i.e., the experimental binding energy of LiF is more than twice that of LiH). In general, LDA functionals tend to over bind in KS-DFT calculations.55 This is the case in the present calculations when comparing the VWN5 binding energy of LiH to that obtained with

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

B3LYP and to experiment. The key observation is that KSCED with external orthogonality reproduces KS-DFT results whether KS-DFT is accurate or is not for a given problem.

c. H2O The PECs of the dissociation of H2O, obtained with KSCED(s, Ext. Orth.) and KS-DFT using the cc-pVTZ basis set together with the VWN5 functional, are shown in Figure 3A. Also shown in Figure 3A is the KSCED(e, Ext. Orth.) PEC. The spectroscopic data obtained from the curves are given in Table 2. As can be seen in Figure 3A, the KSCED(s, Ext. Orth.) and KS-DFT PECs are indistinguishable at all internuclear distances, whereas the KSCED(e, Ext. Orth.) curve is close, but distinguishable at the resolution of the plot, from the other two. In calculations with the extended monomer basis, the one-particle basis for the expansion of KS orbitals of the hydrogen radical ( H • ) subsystem was extended by the inclusion of basis functions centered at the oxygen atom of the hydroxyl radical ( HO • ) fragment. The effect of additional basis functions, perhaps not centered on atoms, has not yet been investigated but might alleviate the discrepancy, albeit small between bases. As shown in Table 2, the KSCED(s, Ext. Orth.) bond length (0.95 Å) and dissociation energy (5.27 eV) obtained with B3LYP agree exactly with reference KS-DFT results to the reported precision. Moreover, these values agree well with the experimental values of 0.9587 Å56 and 5.101 eV,57 respectively. Also, the KSCED(s, Ext. Orth.) bond length (0.81 Å) and dissociation energy (6.13 eV) obtained with VWN5 exactly match reference KS-DFT values, to the given level of precision. Again, as noted previously, the VWN5 functional tends to over bind and leads to a shorter bond length, as is seen in the present case. The KSCED(e, Ext. Orth.)

14 ACS Paragon Plus Environment

Page 15 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

values obtained with VWN5 are also relatively close to the reference KS-DFT ones, differing by 0.06 Å in bond length and 0.08 eV in bond energy. In sharp contrast to the embedding-with-external-orbital-orthogonality DFT-in-DFT results, the conventional KSCED(s) (i.e. DFT-in-DFT, without orbital orthogonality and using LC94 KE functional) results are very poor. To even be displayable, KSCED(s) relative energies in Figure 3B were multiplied by a factor of 10 for comparison with the reference KS-DFT PEC. The obtained spectroscopic data is similarly compromised: the KSCED(s) bond length from VWN5 is 0.2213 Å larger than the experimental value, whereas the KSCED(s) predicted H-O bond energy is negligible (i.e., 0.47 eV compared to 5.101 eV from experiment). This is the first demonstration of the capabilities of the external orthogonality DFT-in-DFT formalism for a covalent bond, even when conventional DFT-in-DFT is unable to describe chemical bonding.

d. H3C-CH3 The PECs of the dissociation of H 3 C − CH 3 from KSCED(s, Ext. Orth.) and KSCED(e, Ext. Orth.) calculations, obtained using PW91 and cc-pVTZ, are shown in Figure 4 and compared with the reference KS-DFT curve obtained with the same functional and basis set. As can be seen in the figure, the KSCED(s, Ext. Orth.) curve matches the reference KS-DFT curve to the width of the line at all internuclear distances, whereas KSCED(e, Ext. Orth.) underbinds. In KSCED(e, Ext. Orth.) calculations, the monomer basis of each methyl radical was extended by the inclusion of basis functions centered at the C-atom of the complementary methyl radical. That is, basis functions of all H-atoms of the complementary subsystem were excluded in expanding the KS orbitals of a given H 3 C • fragment. As noted for H 2 O , it is anticipated that using additional basis functions within the external orthogonality embedding formalism would 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 34

systematically improve the quality of results obtained with the extended monomer basis set in comparison with KSCED(s, Ext. Orth.) and reference KS-DFT, but that study is outside the scope of the present work. Shown in Table 2 are the spectroscopic data corresponding to the H 3 C − CH 3 curves in Figure 4. As can be seen, the same values of equilibrium bond length (1.53 Å) and dissociation energy (4.93 eV) are obtained from the KSCED(s, Ext. Orth.) and KS-DFT methods. The computed bond length is in good agreement with the experimental value of 1.54 Å,58 obtained from electron diffraction studies. On the other hand, both the obtained KSCED(s, Ext. Orth.) and KS-DFT bond energies are 1.02 eV in excess of the experimental value of 3.91 eV.59 The obtained KSCED(e, Ext. Orth.) values of bond length and bond energy compare reasonably well with the reference KS-DFT; the bond length is 0.03 Å too long whereas the binding energy is 0.50 eV too low. Conversely, conventional DFT-in-DFT (i.e., KSCED(s)) predicts H 3 C − CH 3 to be unbound. Since the orbitals, and subsequent densities, from a KSCED calculation will not, in general, agree exactly with those from a KS-DFT calculation of the whole system, the differences in densities were obtained and examined.

Figures 5 and 6 show the density

difference (i.e. ∆ρ = ρKS-DFT – ρKSCED) for the CH3-CH3 molecule at the equilibrium C-C bond length (i.e., 1.53 Å) as a function of distance along the C-C axis and a perpendicular axis, using the supermolecular basis and the extended basis respectively. Several observations can be made; the deviations of KSCED densities from KS-DFT densities: 1) are small; 2) are smooth, with minimal variations in gradient; and 3) are larger for the extended monomer basis than for the supermolecular basis. Moreover, the largest differences are concentrated near the center of the

16 ACS Paragon Plus Environment

Page 17 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

C-C bond, which corroborates our hypothesis that basis set support near the interface between subsystems is an area for improvement.

IV. Conclusions In this article, the results of the first application of the embedding-with-external-orbitalorthogonality DFT-in-DFT method19,21,22 to the cleavage of strong bonds, including first-row covalent bonds, to form radical moieties were investigated. By studies on the ionic bond species LiH and LiF, as well as on the polarized O-H bond in H2O and the unpolarized C-C bond in CH3-CH3, it was found that the embedding-with-external-orbital-orthogonality DFT-in-DFT method is well suited not only for the description of interacting molecules and ions21,22 but also of actual chemical bonds. To our knowledge, for the first time, homolytic bond cleavage into radicals has been accurately described (accuracy assessed with respect to KS-DFT) using a variant of DFT-in-DFT embedding theory that does not require any full system calculations. The PECs of the dissociated bonds obtained from the external orthogonality embedding protocol using the supermolecular basis, i.e., KSCED(s, Ext. Orth.), were found to be indistinguishable (to within the width of a line in the graphs) at all internuclear distances when compared with reference KS-DFT curves. In fact, total energies from KSCED(s, Ext. Orth.) agreed with those from KS-DFT at least one part in 106 at all geometries of the studied bonds. Consequently, the computed equilibrium bond lengths, harmonic frequencies, and dissociation energies from KSCED(s, Ext. Orth.) were essentially indistinguishable from the reference KSDFT values to the reported precision. Of course, although results herein were compared with experimental data, no variant of DFT-in-DFT should be expected to be more accurate than KSDFT itself excepting fortuitous cases. The new embedding method is thus accurate with respect

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

to DFT, and seems to be mainly by inadequacies of DFT exchange-correlation functionals. It was also observed that when deviations in subsystem densities (i.e., relative to whole molecule KS-DFT) occur, as a result of confining one electron (e.g., alpha electron) to one subsystem and the other to another subsystem, they are small as are their gradients. Parenthetically, this suggests, if the observations for H3C-CH3 prove to be more general that the external orthogonality protocol is an effective means of localizing orbitals and densities. The results of this work corroborate in part our previous finding that the extended monomer basis approach to expanding the KS orbitals of a subsystem is a suitable alternative to the more computationally demanding supermolecular basis approach.21 However, it should be noted that there is sometimes a fall-off of accuracy with use of extended monomer basis, as previously construed, with respect to the supermolecular basis. Although it must be appreciated that even when the results are not as accurate as those with the supermolecular basis, they are vastly superior to DFT-in-DFT, without orbital orthogonality and using LC94 KE functional. We anticipate that the accuracy of results obtained with the extended monomer approach within the new embedding method should be expected to be systematically improvable with additional basis functions in boundary regions that are distant from atom-centered basis functions; investigation of this is underway. Because of the small molecules used in this study, variations in the size of the buffer region (i.e., the “extended” of extended monomer basis) were not able to be explored, since extending the extended monomer basis to include those of the next-nearestneighbor atoms was equivalent to using the supermolecular basis; systematic investigation of buffer size is warranted and will be investigated with second-generation computer codes. Although computational efficiency was not a consideration in our pilot computer code, it is perhaps valuable to remark that enforcing external orbital orthogonality did not increase

18 ACS Paragon Plus Environment

Page 19 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

computer time; in fact, in our largest study, of the dissociation of H3C-CH3, the external orthogonality calculations took slightly less time than those without, as a result of quicker convergence of the density. Lastly, we recall that enforcing external orthogonality in wave function-in-DFT embedding theory follows the same logic as in DFT-in-DFT.19,21 Consequently, it should be anticipated that enforcing the external orthogonality constraint within wave function-in-DFT will produce chemically accurate results in situations in which one of the subsystems cannot be adequately described using DFT (e.g., some excited states19). We anticipate future research efforts in this direction.

Acknowledgements The authors gratefully acknowledge the National Science Foundation (Grant No. EPS-0814442) for funding this project.

References 1. Wesolowski, T. A.; Warshel, A. Frozen Density Functional Approach for Ab Initio Calculations of Solvated Molecules. J. Phys. Chem. 1993, 97, 8050–8053. 2. Wesolowski, T. A.; Weber, J. Kohn-Sham Equations with Constrained Electron Density: An Iterative Evaluation of the Ground-State Electron Density of Interacting Molecules. Chem. Phys. Lett. 1996, 248, 71–76. 3. Jacob, C. R.; Neugebauer, J.; Visscher, L. A. Flexible Implementation of Frozen‐Density Embedding for Use in Multilevel Simulations. J. Comput. Chem. 2008, 29, 1011–1018. 4. Wesolowski, T. A. Embedding a Multi-determinantal Wavefunction in Orbital-Free Environment. Phys. Rev. A 2008, 77, 012504/1–9. 5. Khait, Y. G.; Hoffmann, M. R. Embedding Theory for Excited States. J. Chem. Phys. 2010, 133, 044107/1–6. 6. Huang, C.; Pavone, M.; Carter, E. A. Quantum Mechanical Embedding Theory Based on a Unique Embedding Potential. J. Chem. Phys. 2011, 134, 154110/1–11. 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

7. Manby, F. R.; Stella, M.; Goodpaster, J. D.; T. F. Miller, T. F. A Simple, Exact DensityFunctional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564–2568. 8. Krishtal, A.; Sinha, D.; Genova, A.; Pavanello, M. Subsystem Density-Functional Theory as an Effective Tool for Modeling Ground and Excited States, Their Dynamics and Many-Body Interactions. J. Phys.: Condens. Matter 2015, 27, 183202/1–27. 9. Wesolowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891–5928. 10. Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and its Gradient. J. Chem. Soc., Perkin Trans. 1993, 2, 799–805. 11. Tomasi, J.; Mennucci, B.; Cances, E. The IEF Version of the PCM Solvation Method: An Overview of a New Method Addressed to Study Molecular Solutes at the QM Ab Initio Level. Theochem 1999, 464, 211–226. 12. Scalmani, G.; Frisch, M. J. Continuous Surface Charge Polarizable Continuum Models of Solvation. I. General formalism. J. Chem. Phys. 2010, 132, 114110/1–15. 13. Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Chemistry. J. Chem. Phys. 1996, 105, 1968–1986. 14. Kairys, V.; Jensen, J. H. QM/MM Boundaries Across Covalent Bonds:  A Frozen Localized Molecular Orbital-Based Approach for the Effective Fragment Potential Method. J. Phys. Chem. A 2000, 104, 6656–6665. 15. Adamovic, I.; Freitag, M. A.; Gordon, M. S. Density Functional Theory Based Effective Fragment Potential. J. Chem. Phys. 2003, 118, 6725–6732. 16. Wu, F.; Liu, W.; Zhang, Y.; Li, Z. Linear-Scaling Time-Dependent Density Functional Theory Based on the Idea of “From Fragments to Molecule”. J. Chem. Theory Comput. 2011, 7, 3643–3660. 17. Song, L.; Han, J.; Lin, Y.; Xie, W.; Gao, J. Explicit Polarization (X-Pol) Potential Using Ab Initio Molecular Orbital Theory and Density Functional Theory. J. Phys. Chem. A 2009, 113, 11656–11664. 18. Wesolowski, T. A. One Electron Equations for Embedded Electron Density: Challenge for Theory and Practical Payoffs in Multi-level Modelling of Complex Polyatomic Systems. In Computational Chemistry: Reviews of Current Trends, Vol. 10; Leszczynski, J., Ed.; World Scientific: Singapore, 2006; pp 1–82. 19. Khait, Y. G.; Hoffmann, M. R. On the Orthogonality of Orbitals in Subsystem Kohn-Sham Density Functional Theory. Annu. Rep. Comput. Chem. 2012, 8, 53–70. 20. Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F., III. Density Functional Theory Embedding for Correlated Wavefunctions: Improved Methods for Open-Shell Systems and 20 ACS Paragon Plus Environment

Page 21 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Transition Metal Complexes. J. Chem. Phys. 2012, 137, 224113/1–10. 21. Tamukong, P. K.; Khait, Y. G.; Hoffmann, M. R. Density Differences in Embedding Theory with External Orbital Orthogonality. J. Phys. Chem. A 2014, 118, 9182–9200. 22. Tamukong, P. K. Extension and Applications of the GVVPT2 Method to the Study of Transition Metals. Ph.D. Dissertation, University of North Dakota, Grand Forks, ND, 2014. 23. Unsleber, J. P.; Neugebauer, J.; Jacob, C. R. No Need for External Orthogonality in Subsystem Density-Functional Theory. Phys. Chem. Chem. Phys. 2016, 18, 21001–21009. 24. Fux, S.; Jacob, Ch. R.; Neugebauer, J.; Visscher, L.; Reiher, M. Accurate Frozen-Density Embedding Potentials as a First Step towards a Subsystem Description of Covalent Bonds. J. Chem. Phys. 2010, 132, 164101/1–18. 25. Goodpaster, J. D.; Ananth, N.; Manby, F. R.; Miller, T. F. Exact Nonadditive Kinetic Potentials for Embedded Density Functional Theory. J. Chem. Phys. 2010, 133, 084103/1– 10. 26. Huang, C.; Carter, E. A. Potential-Functional Embedding Theory for Molecules and Materials. J. Chem. Phys. 2011, 135, 194104/1–17. 27. Elliott, P.; Burke, K.; Cohen, M. H.; Wasserman, A. Partition Density Functional Theory. Phys. Rev. A 2010, 82, 024501/1–4. 28. Beyhan, S. M.; Götz, A. W.; Jacob, Ch. R.; Visscher, L. The Weak Covalent Bond in NgAuF (Ng= Ar, Kr, Xe): A Challenge for Subsystem Density Functional Theory. J. Chem. Phys. 2010, 132, 044114/1–9. 29. Dułak, M.; Wesołowski, T. A. On the Electron Leak Problem in Orbital-Free Embedding Calculations. J. Chem. Phys. 2006, 124, 164101/1–9. 30. Humbert-Droz, M.; Zhou, X.; Shedge, S. V.; Wesolowski, T. A. How to Choose the Frozen Density in Frozen-Density Embedding Theory-Based Numerical Simulations of Local Excitations? Theor. Chem. Acc. 2014, 133, 1405/1–20. 31. Jacob, C. R.; Beyhan, S. M.; Visscher, L. Exact Functional Derivative of the Nonadditive Kinetic-Energy Bifunctional in the Long-Distance Limit. J. Chem. Phys. 2007, 126, 234116/1–13. 32. Fux, S.; Kiewisch, K.; Jacob, Ch. R.; Neugebauer, J.; Reiher, M. Analysis of Electron Density Distributions from Subsystem Density Functional Theory Applied to Coordination Bonds. Chem. Phys. Lett. 2008, 461, 353–359. 33. Kiewisch, K.; Eickerling, G.; Reiher, M.; Neugebauer, J. Topological Analysis of Electron Densities from Kohn-Sham and Subsystem Density Functional Theory. J. Chem. Phys. 2008, 128, 044114/1–15. 34. Goodpaster, J. D.; Barnes, T. A.; Miller, T. F., III. Embedded Density Functional Theory for Covalently Bonded and Strongly Interacting Subsystems. J. Chem. Phys. 2011, 134, 164108/1–9. 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 34

35. Lembarki, A.; Chermette, H. Obtaining a Gradient-Corrected Kinetic-Energy Functional from the Perdew-Wang Exchange Functional. Phys. Rev. A 1994, 50, 5328–5331. 36. Mura, M. E.; Knowles, P. J. Improved Radial Grids for Quadrature in Molecular DensityFunctional Calculations. J. Chem. Phys. 1996, 104, 9848–9858. 37. Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200–1211. 38. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671–6687. 39. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Erratum: Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1993, 48, 4978–4978. 40. Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange J. Chem. Phys. 1993, 98, 5648–5652. 41. Becke, A. D. A New Mixing of Hartree-Fock and Local Density-Functional Theories. J. Chem. Phys. 1993, 98, 1372–1377. 42. Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789. 43. Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. 44. Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. 45. Wharton, L.; Klemperer, W.; Gold, L. P.; Strauch, R.; Gallagher, J.; Derr, V. Microwave Spectrum, Spectroscopic Constants, and Electric Dipole Moment of Li6F19. J. Chem. Phys. 1963, 38, 1203–1210. 46. Bulewicz, E.; Phillips, L.; Sugden, T. Determination of Dissociation Constants and Heats of Formation of Simple Molecules by Flame Photometry. Part 8.—Stabilities of the Gaseous Diatomic Halides of Certain Metals. Trans. Faraday Soc. 1961, 57, 921–931. 47. Hedderich, H.; Frum, C.; Engleman, R., Jr.; Bernath, P. The Infrared Emission Spectra of LiF and HF. Can. J. Chem. 1991, 69, 1659–1671. 48. Suresh, S.; Naik, V. Hydrogen Bond Thermodynamic Properties of Water from Dielectric Constant Data. J. Chem. Phys. 2000, 113, 9727–9732. 49. Huang, P.; Carter, E. A. Advances in Correlated Electronic Structure Methods for Solids, Surfaces, and Nanostructures. Annu. Rev. Phys. Chem. 2008, 59, 261–290. 22 ACS Paragon Plus Environment

Page 23 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

50. Perdew, J. P.; Schmidt, K. Jacob's Ladder of Density Functional Approximations for the Exchange-Correlation Energy. AIP Conf. Proc. 2001, 577, 1–20. 51. Dulak, M.; Kaminski, J. W.; Wesolowski, T. A. Equilibrium Geometries of Non-Covalently Bound Intermolecular Complexes Derived from Subsystem Formulation of Density Functional Theory. J. Chem. Theory Comput. 2007, 3, 735–745. 52. Chan, Y.; Harding, D.; Stwalley, W.; Vidal, C. Inverted Perturbation Approach (IPA) Potentials and Adiabatic Corrections of the X  1Σ+ State of the Lithium Hydrides near the Dissociation Limits. J. Chem. Phys. 1986, 85, 2436–2444. 53. Vidal, C.; Stwalley, W. Potential Energy Curves and Adiabatic Corrections of Weakly Bound States: Application to the LiH B 1Π State. J. Chem. Phys. 1984, 80, 2697–2703. 54. Dulick, M.; Zhang, K.-Q.; Guo, B.; Bernath, P. Far- and Mid-Infrared Emission Spectroscopy of LiH and LiD. J. Mol. Spectrosc. 1998, 188, 14–26. 55. Kohn, W.; Becke, A. D.; Parr, R. G. Density Functional Theory of Electronic Structure. J. Phys. Chem. 1996, 100, 12974–12980. 56. Cook, R. L.; De Lucia, F. C.; Helminger, P. Microwave Spectra of Molecules of Astrophysical Interest V. Water Vapor. J. Mol. Spectrosc. 1974, 53, 62–70. 57. Maksyutenko, P.; Rizzo, T. R.; Boyarkin, O. V. A Direct Measurement of the Dissociation Energy of Water. J. Chem. Phys. 2006, 125, 181101/1–3. 58. Pauling, L.; Brockway, L. O. Carbon—Carbon Bond Distances. The Electron Diffraction Investigation of Ethane, Propane, Isobutane, Neopentane, Cyclopropane, Cyclopentane, Cyclohexane, Allene, Ethylene, Isobutene, Tetramethylethylene, Mesitylene, and Hexamethylbenzene. Revised Values of Covalent Radii. J. Am. Chem. Soc. 1937, 59, 1223– 1236. 59. Blanksby, S. J.; Ellison, G. B. Bond Dissociation Energies of Organic Molecules. Acc. Chem. Res. 2003, 36, 255–263.

23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 34

Table 1. Equilibrium bond lengths, dissociation energies, and harmonic frequencies of the LiF and LiH molecules, obtained with the external orthogonality embedding method compared with reference KS-DFT, DFT-in-DFT (using LC94 KE functional) and with experimental values. KS-DFT

KSCED(s, Ext. Orth.)

KSCED(s)

Expt.

LiF (B3LYP Functional) Re (Å)

1.57

1.57

3.30

1.564a

De (eV)

5.43

5.43

1.18

5.898±0.35b

ωe (cm-1)

917.9

917.9

144.2

910.57c

LiH (B3LYP Functional) Re (Å)

1.59

1.59

1.43

1.596d

De (eV)

2.53

2.53

1.65

2.52e

ωe (cm-1)

1416.7

1417.0

1482.1

1405.50f

LiH (VWN5 Functional) Re (Å)

1.61

1.61

1.49

De (eV)

2.63

2.63

1.87

ωe (cm-1)

1361.0

1361.0

1422.9

a

Ref. 45; bRef. 46; cRef. 47; dRef. 52; eRef. 53; fRef. 54.

24 ACS Paragon Plus Environment

Page 25 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 2. Equilibrium bond lengths and dissociation energies of the O-H bond in H2O and the CC bond in ethane, obtained with the external orthogonality embedding method compared with reference KS-DFT, DFT-in-DFT (using LC94 KE functional) and with experimental values. KS-DFT

KSCED(s, Ext. Orth.)

KSCED(e, Ext. Orth.)

KSCED(s)

Expt.

H 2 O (B3LYP Functional) Re (Å)

0.95

0.95

0.96

1.30

0.9587a

De (eV)

5.27

5.27

5.24

0.22

5.101b

H 2 O (VWN5 Functional) Re (Å)

0.81

0.81

0.75

1.18

De (eV)

6.13

6.13

6.05

0.47

H 3 C − CH 3 (PW91 Functional) Re (Å)

1.53

1.53

1.56

unbound

1.54±0.02c

De (eV)

4.93

4.93

4.43

-

3.91d

a

Ref. 56; bRef. 57; cRef. 58; dRef. 59 (Note that this is the dissociation enthalpy at 298 K).

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 34

Figure Captions Figure 1. Relative energy dissociation curves of LiF obtained with the embedding-with-externalorbital-orthogonality DFT-in-DFT (green) and KS-DFT calculations (black). The green curve completely obscures the KS-DFT curve.

Figure 2. Relative energy dissociation curves of LiH obtained using B3LYP or VWN5 with the embedding-with-external-orbital-orthogonality DFT-in-DFT (olive: B3LYP; blue: VWN5) and KS-DFT calculations (black: B3LYP; cyan: VWN5). The olive and blue curves completely obscure the KS-DFT curves.

Figure 3. (a) Relative energy curves of the dissociation of H2O into OH and H fragments, obtained

with

the

embedding-with-external-orbital-orthogonality

DFT-in-DFT

(green:

supermolecular basis; magenta: extended monomer basis) and KS-DFT calculations (black). The green curve completely obscures the KS-DFT curve. (b) Comparison of minimum of conventional DFT-in-DFT (red: KSCED(s)) curve with KS-DFT (black). Note that the KSCED(s) relative energies have been multiplied by a factor of 10 to allow comparison of the two curves.

Figure 4. Relative energy curves of the symmetric dissociation of H3C–CH3 obtained with the embedding-with-external-orbital-orthogonality DFT-in-DFT (green: supermolecular basis; magenta: extended monomer basis) and KS-DFT calculations (black). The green curve completely obscures the KS-DFT curve.

26 ACS Paragon Plus Environment

Page 27 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Density differences between reference KS-DFT and embedding-with-external-orbitalorthogonality DFT-in-DFT densities obtained using the supermolecular basis set (i.e., ∆ρ = ρ DFT − ρ KSCED (super ) ) for the CH3-CH3 molecule at the equilibrium C-C bond length (i.e.,

1.53 Å). Both sets of calculations used the PW91 functional and the cc-pVTZ basis set.

Figure 6. Density differences between reference KS-DFT and embedding-with-external-orbitalorthogonality DFT-in-DFT densities obtained using the extended monomer basis set (i.e., ∆ρ = ρ DFT − ρ KSCED (ext ) ) for the CH3-CH3 molecule at the equilibrium C-C bond length (i.e., 1.53

Å). Both sets of calculations used the PW91 functional and the cc-pVTZ basis set.

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 34

TOC Graphic

28 ACS Paragon Plus Environment

Page 29 of 34

5

4

R e la tiv e E n e rg y (e V )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

The Journal of Physical Chemistry

3

2

1

0 1

2

3

4

5

B o n d L e n g th (Å

6

7 )

ACS Paragon Plus Environment

8

9

1 0

The Journal of Physical Chemistry

3 .0

2 .5

R e la tiv e E n e r g y ( e V )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 30 of 34

2 .0

1 .5

1 .0

0 .5

0 .0 1

2

3

4 L i- H

5

B o n d L e n g th (Å )

ACS Paragon Plus Environment

6

7

8

9

Page 31 of 34

1 2 1 0

(a ) 8

R e la tiv e E n e rg y (e V )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

The Journal of Physical Chemistry

6 4 2 0 1 0

(b ) 8 6 4 2 0 0 .0

1 .0

2 .0

3 .0

4 .0

H -O B o n d L e n g th (Å ) ACS Paragon Plus Environment

5 .0

6 .0

The Journal of Physical Chemistry

5 .5 5 .0 4 .5 4 .0

R e la tiv e E n e rg y (e V )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 32 of 34

3 .5 3 .0 2 .5 2 .0 1 .5 1 .0 0 .5 0 .0

- 0 .5 1

2

3

4

5

C -C B o n d L e n g th (Å ) ACS Paragon Plus Environment

6

7

8

Page 33 of 34

4

0 .0 0 4 4 0 0 - 0 .0 0 2 6 7 5

2

- 0 .0 0 9 7 5 0

C - C A x is ( a .u .)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

The Journal of Physical Chemistry

- 0 .0 1 6 8 3

0

- 0 .0 2 3 9 0 - 0 .0 3 0 9 8

-2

- 0 .0 3 8 0 5 - 0 .0 4 5 1 3 - 0 .0 5 2 2 0

-4 -4

-2

0

2

⊥ ACS A x Paragon i s ( a . Plus u . ) Environment

4

The Journal of Physical Chemistry

4

0 .0 8 8 0 0 0 .0 7 6 9 8

2

0 .0 6 5 9 5

C - C a x is ( a .u .)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

Page 34 of 34

0 .0 5 4 9 3

0

0 .0 4 3 9 0 0 .0 3 2 8 8

-2

0 .0 2 1 8 5 0 .0 1 0 8 3 - 2 .0 0 0 E - 0 4

-4 -4

-2

0

2

⊥ACS a x Paragon i s ( a . Plus u . ) Environment

4