Potential Functional Embedding Theory at the ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Potential Functional Embedding Theory at the Correlated Wavefunction Level, Part I: Mixed Basis Set Embedding Jin Cheng, Florian Libisch, Kuang Yu, Mohan Chen, Johannes M. Dieterich, and Emily A. Carter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01010 • Publication Date (Web): 26 Jan 2017 Downloaded from http://pubs.acs.org on January 30, 2017

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.

Journal of Chemical Theory and Computation 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 41

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

Journal of Chemical Theory and Computation

Potential Functional Embedding Theory at the Correlated Wavefunction Level, Part I: Mixed Basis Set Embedding

Jin Cheng,1 Florian Libisch,2,4 Kuang Yu,2 Mohan Chen,2 Johannes M. Dieterich,2 and Emily A. Carter3,a) 1

Department of Chemistry, 2Department of Mechanical and Aerospace Engineering, and 3School of

Engineering and Applied Science, Princeton University, Princeton, New Jersey 08544 4

Institute for Theoretical Physics, Vienna University of Technology, 1040 Vienna, Austria, EU

ABSTRACT Embedding theories offer an elegant solution to overcome intrinsic algorithmic scaling and accuracy limitations of simulation methods. These theories also promise to achieve the accuracy of high-level electronic structure techniques at near the computational cost of much less accurate levels of theory by exploiting positive traits of multiple methods. Of crucial importance to fulfilling this promise is the ability to combine diverse theories in an embedding simulation. However, these methods may utilize different basis set and electron-ion potential representations. In this first part of a two-part account of implementing potential functional embedding theory (PFET) at a correlated wavefunction level, we discuss remedies to basis set and electron-ion potential discrepancies and assess the performance of the PFET scheme with mixed basis sets. _____________________________ a)

Author to whom correspondence should be addressed. Electronic mail: [email protected].

1 ACS Paragon Plus Environment


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

1 Introduction Embedding theories promise simulations with the accuracy of a high-level electronic structure method at significantly reduced computational effort.1 These theories divide a system into two or more subsystems, each potentially studied at different levels of theory, and approximate the interactions between the subsystems. Both the division into subsystems and the choice of theory for the subsystems are based on their physical characteristics, e.g., for chemisorption, the adsorbate and atoms comprising the adsorption site might qualify as one subsystem and the rest of the surface and subsurface as another. One of the more interesting method combinations is a correlated wavefunction (CW) for the subsystem of physical interest, which we will refer to henceforth as the embedded region, and density functional theory (DFT) for the other subsystem, the environment.1e, 2 This combination allows for wavefunction-quality descriptions of localized physical and chemical phenomena, while retaining the ability to include the influence of extended condensed matter. A defining aspect of embedding theories is the representation of the interaction between the subsystems.1e Frozen density embedding (FDE) approaches have been employed successfully;2e,

2h

however, as the name indicates, they lack response of the environment to

density changes in the embedded region. A more accurate representation requires inclusion of these effects. Density functional embedding theory (DFET)2i and its antecedents2b-d, 2f, 2g have been used with considerable success for embedded CW studies of, e.g., ground and excited state molecule-metal surface interactions.1c, 2d, 2g, 3 The partition-DFT proposed by Wasserman and coworkers shares a great similarity with DFET applied at the DFT level, and likewise does not rely on the frozen density approximation.4 However, embedded CW within DFET still neglects the higher-order mutual polarization between subsystems induced by the CW description of the


Page 2 of 41

Page 3 of 41

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


embedded region. The self-consistent potential functional embedding theory (PFET),5 however, provides a much more flexible formalism that theoretically accounts for such mutual polarization. When first introduced, PFET was verified within a planewave (PW) pseudopotential DFT-only framework.5 However, practical extension of PFET towards general applicability and broader method choices, including CW methods, requires significant method and program developments, as well as added capabilities for every step of the PFET calculation. In this two-part contribution, we extend the implementation of PFET from a purely PWpseudopotential DFT-based scheme5 to a more general one that permits use of mixed basis sets and multiple electronic structure methods. Verification of the PFET implementation at the DFTonly level with mixed basis set and electron-ion potential representations is an essential prerequisite prior to implementing full-blown CW/DFT PFET. (Here “ion” refers to a nucleus plus core electrons and “CW/DFT” refers to treatment of one subsystem with CW and one with DFT; hence “DFT/DFT” refers to a DFT-only treatment, i.e., DFT treatment within all subsystems.) In part I, we describe a consistent implementation of electron-ion potentials and mixed basis sets primarily at the DFT/DFT level. Tests of accuracy and consistency of these components lay the foundation for part II, which focuses on the performance of hybrid CW/DFT methods within PFET. After briefly reviewing the PFET and DFET formalisms (Section 2.1), we introduce and analyze numerical issues caused by basis set incompatibility (Subsection 2.2.1) and introduce the concept of spillage as a way to assess basis set performance in Subsection 2.2.2. We then use principal component analysis (PCA) to optimize Gaussian-type orbital (GTO) basis sets appropriate for the PW-based optimized effective potential (OEP) calculations used in PFET (Subsection 2.2.3). We close Section 2 by a brief discussion of the numerical details of our



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 41

formalism (Section 2.3). In Section 3, we present our results on representative test cases (Section 3.1) to verify that the optimized GTO basis sets (Section 3.2) are satisfactory for PW-based OEP calculations (Section 3.3) and therefore suitable for use in PFET simulations. PFET calculations at the DFT/DFT level with mixed basis sets for CO adsorption on an aluminum cluster are then used to benchmark PFET performance under different chemical bonding scenarios (metallic and donor-acceptor bonds, Section 3.4).1 We also compare PFET with DFET, and extract insights anticipating application of PFET at the CW/DFT level. 2 Theory and implementation 2.1 Embedding theories We only discuss here relevant aspects of DFET2i and PFET;5 an in-depth discussion is given in Ref. 1e. A detailed analysis of approximations made and associated errors is presented in part II. Both DFET and PFET are orbital-free embedding theories (OFETs) that rely solely on an electron density representation. These theories partition the density of the total system ρtot into a sum of subsystem densities ρi : N

ρtot = ∑ρi

(1)

i =1

v

where interactions between the subsystems are described by an embedding potential Vemb ( r ) . Explicit calculation of couplings between wavefunctions are not required.6 This approach is

1

More test cases of different bonding types are discussed in the supporting information, Sections 3.1 through 3.3. 4 ACS Paragon Plus Environment

Page 5 of 41

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


compatible with a wide range of theoretical treatments of the individual subsystems because it does not require any quantities besides the density, such as orbitals or density matrices. We consider only two regions I and II for simplicity. We further identify subsystem I as the physically more challenging (embedded) region while subsystem II comprises the environment. The Schrödinger equation for the embedded region of all OFETs then becomes:

(H I + Vemb [ ρI , ρtot ])ΨI = EI ΨI

(2)

where H I is the electronic Hamiltonian, EI the energy, and Ψ I the wavefunction of the embedded region, and the embedding potential Vemb [ ρI , ρtot ] is designed to represent the environment’s influence on the embedded region. We typically further partition this interaction potential into its kinetic Ts , exchange-correlation (XC) Vxc , Hartree VJ , and electron-ion VIE components:1b  δT [ρ ] δT [ρ ]  Vemb [ ρ I , ρtot ] =  s tot − s I  + Vxc [ ρtot ] − Vxc [ ρ I ] + VJ [ ρtot ] − VJ [ ρ I ] + VIE , tot − VIE , I δρ I   δρtot

(3)

Both DFET and PFET assume a common embedding potential to represent the influence of subsystem I on II and vice-versa. For DFET we also demand that the sum of DFT subsystem densities optimized in the presence of the embedding potential reproduce the total DFT density. The embedding potential satisfying both of these requirements is unique7 and can be solved for by maximizing an extended Wu-Yang functional known from OEP theory:2i, 8 r r r W [Vemb ] = E% I [ ρ I , Vemb ] + E% II [ ρ II , Vemb ] − ∫ ρtot ( r ) ⋅ Vemb ( r ) dr


(4)


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 41

where the subsystem energy E% K [ ρ K ,Vemb ], ( K = I , II ) includes the interaction with the r r r embedding potential ( E% K [ ρ K , Vemb ] = E K [ ρ K ] + ∫ ρ K ( r ) ⋅ Vemb ( r ) dr ). The total reference density

is obtained from a DFT calculation of the entire system, typically in a PW basis representation. Embedded subsystem densities and energies are obtained with PW-DFT during the OEP procedure. The optimized embedding potential thus obtained can then be used within higherlevel electronic structure (e.g., CW) calculations on subsystem I using Eq. (2). In DFET, the embedded region’s interaction with the environment is only evaluated at the DFT level and the environment density is fixed in subsequent embedded CW calculations. PFET addresses this shortcoming by allowing the mutual polarization between subsystems to account for the more refined CW description of the embedded region (or, more generally, of any subsystem). Instead of deriving the embedding potential by requiring reproduction of the total DFT density, PFET variationally optimizes the total energy as a functional of the embedding potential:

Etot [Vemb ] = E% I [Vemb ] + E% II [Vemb ] + E%int  ρI [Vemb ] , ρII [Vemb ]

(5)

This functional is directly minimized with the subsystem energies able to be calculated at different levels of theory.5 Hence, polarization of the environment by density changes in the embedded region can occur (and vice-versa, until self-consistency is reached). The interaction between subsystems in PFET ( E% int  ρ I [Vemb ] , ρ II [Vemb ] ) is separated as usual into exchangecorrelation, Coulomb, cross-subsystem electron-ion, kinetic energy terms, and a densityembedding potential integral to cancel out the corresponding terms in the subsystem energies. The interaction kinetic energy is computed using the Wu-Yang OEP method.8a PFET is


Page 7 of 41

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


equivalent to DFET in the infinite basis set limit with DFT treatment of all subsystems. At mixed levels of theory, the polarization feedback built into PFET provides a self-consistent, more refined description. In DFET, the density and OEP optimizations in Eq. (4) are performed using only PW basis sets. By contrast, the Wu-Yang OEP8a in PFET employs a PW basis set to invert the subsystem densities, each of which could be obtained with GTO basis sets as discussed below. Here, inverting a density means to find, for a given density ρ , the effective Kohn-Sham potential that reproduces ρ using, e.g., the Wu-Yang OEP algorithm. Differing basis set representations introduce numerical challenges for this step, which we discuss, analyze, and ultimately remedy in the following sections. 2.2 GTO basis set optimization Implementation of CW/DFT PFET requires a careful look into the underlying basis representations. CW calculations usually employ GTO basis sets9 whereas PW-based DFT calculations are typically used to treat extended systems. Accordingly, we will use GTO basis functions for high-level treatments of the region of interest while relying on PWs to describe the environment and for OEP calculations of the exact kinetic energy Ts [ ρ ] and kinetic potential VTs

=,

δ Ts [ ρ ] , to ensure numerical accuracy.5, δρ

10

These mixed bases introduce possible

incompatibilities. For example, a GTO-CW/PW-DFT PFET iteration involves a PW-based OEP calculation, in which the Kohn-Sham (KS) equations are inverted using the GTO-based embedded region CW density to obtain Ts [ ρGTO −CW ] and

δTs [ ρGTO−CW ] . Previous studies suggest δρGTO−CW

that GTO-based OEP calculations are prone to numerical instabilities in the final potential.2h, 11 7 ACS Paragon Plus Environment


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 41

To avoid these problems, we implement all OEP calculations in a PW-basis throughout our embedding calculations. Differences also exist between the electron-ion effective core potentials (ECPs) used with GTO bases and electron-ion pseudopotentials used with PW bases. Conventional ECPs and pseudopotentials, although similar in physical concept, are mostly not interchangeable. We must construct GTO bases and use electron-ion potentials consistent across subsystems as a prerequisite to accurate GTO-CW/PW-DFT embedding. 2.2.1 Inconsistencies between GTO and PW basis sets GTO basis sets are localized, atom-centered bases widely used in most quantum chemistry programs.9 PW bases inherently include periodicity and are mainly used in solid-state DFT software.12 Here, we first introduce these two basis set types and then discuss basis set incompatibility in CW/DFT PFET. v GTO basis functions are written as Φ ( r ) = Rl ( r ) Ylm (θ , φ ) , with the radial part Rl ( r )

represented by a contraction of primitive Gaussian functions P

Rl ( r ) = r l ∑c p A ( l , α p ) e

−α p r 2

(6)

p =1

in which cp is a contraction coefficient, A ( l , α p ) a normalization factor, and αp an exponent. Many GTO basis sets suitable for quantum chemistry have been constructed.13 The geometric shape of a GTO basis function resembles an atomic orbital, except close to the nucleus where the cusp is missing; these basis sets thus describe electron distributions in atoms well, even with relatively few functions. GTO basis sets also straightforwardly represent molecular orbitals in a multi-atom system with the “linear combination of atomic orbitals” ansatz.14


Page 9 of 41

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


In contrast to the locality of atom-centered GTO bases, the completely delocalized PW bases were introduced because their inherently periodic nature makes them a natural choice for computations of the electronic structure of crystals. PW basis functions are written as

1 iGvi1i2i3 ⋅rv e Ω

v

ϕi i i ( r ) = 123

(7)

where i j indexes from 1 to the number of grid points in j-th direction, Nj , and Ω is the volume

v

v

v

of the unit cell defined by the three lattice vectors a1 , a2 , and a3 . By introducing the reciprocal lattice vectors v v v a2 × a3 b1 = 2π Ω v v v a ×a b2 = 2π 3 1 Ω v v v a ×a b3 = 2π 1 2 Ω

(8)

v N v  N v  N v  Gi1i2i3 =  i1 − 1  b1 +  i2 − 2  b2 +  i3 − 3  b3 2  2  2    

(9)

v the wave-vector Gi1i2i3 may be defined as

and the PW expansion of the wavefunction becomes N1 N 2 N3 v v v v 1 N1 N2 N3 v iGi1i2i3 ⋅r % un = ∑∑∑ u%n Gi1i2i3 ϕi1i2i3 ( r ) = u G e ∑∑∑ n i1i2i3 Ω i1 =1 i2 =1 i3 =1 i1 =1 i2 =1 i3 =1

(

v

with u% n ( G i i i

1 2 3

)

)

(

)

(10)

the expansion coefficient for the PW basis function ϕi1i2i3 . The number of PW

basis functions ( N1 × N 2 × N 3 ) is controlled by the PW kinetic energy cutoff Ecut with the v constraint 1 G 2 ≤ E cut in atomic units. The advantage of PW basis sets is that they can be 2

systematically improved and made effectively complete simply by increasing Ecut . By contrast,



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 41

the extension of GTOs to the complete basis set limit is quite challenging: although several proposals have been made, the general use of GTOs remains restricted to existing systematic basis sets.15 Fortunately, a much smaller number of GTOs can typically be used to obtain accuracy similar to a PW description. Both types of basis sets serve their own respective purposes well, although they are not designed to be compatible with each other. One symptom of this incompatibility is that the density derived from a PW-DFT calculation typically does not match the density of a GTO-DFT calculation with the same XC functional. Further, inversion of the KS equations in a PW-OEP scheme using a GTO-DFT density produces an OEP potential Voep that unphysically deviates from the KS effective potential VKS . In fact, Voep derived from a ground-state KS-DFT density should exactly equal VKS in the complete basis set limit. The mismatch in the density between the two basis sets manifests in a spurious, unphysical contribution to Voep . Inversion of the KS equations for a GTO-CW density therefore would yield a potential Voep of questionable value, despite containing physically appropriate deviations from VKS accounting for the difference in CW versus DFT densities. L. Füsti-Molnar and P. Pulay recognized the incompatibility between these two bases and proposed a mixed-basis augmented PW method to address the problem.16 However, employing such a mixed basis set in the embedding method at hand leads to various technical difficulties and thus becomes impractical. We must therefore devise a way to optimize GTO bases so as to alleviate the basis-set-induced difference in the density and hence eliminate spurious contributions to Voep .


Page 11 of 41

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


2.2.2 Spillage Spillage is a metric to quantify how well a localized basis set represents a set of PWbased wavefunctions.17 We use this metric to optimize GTO basis sets for our PFET calculations. Below, φi denotes atomic basis functions and ψ i

the target, orthonormal, PW-based

wavefunctions. The spillage of a given basis set { φi } is then defined as:17a

s=

1 n ∑ ψ i 1 − pˆ ψ i n i =1

(11)

Here, n is the number of wavefunctions and pˆ is the projection operator spanned by the atomic basis functions: −1 φν pˆ = ∑ φµ Sµν

(12)

µν

where S is the overlap matrix between the basis functions, evaluated as Sij = φi φ j . We re-write the above equations in matrix form for simplicity. The matrix W , representing the overlap between an atomic basis function and a target wavefunction, is evaluated numerically as Wij = φi ψ j . With this matrix, we re-arrange Eq. (11) as:

1 n ∑ ψ i 1 − pˆ ψ i n i =1 1 n 1 n = ∑ ψ i I ψ i − ∑ ψ i pˆ ψ i n i =1 n i =1 1 n −1 = 1 − ∑∑ ψ i φµ S µν φν ψ i n i =1 µν

s=

= 1− Here the { ψ i

(13)

1 n 1 W ′S −1W ) = 1 − Tr (W ′S −1W ) ( ∑ ii n i =1 n

} are orthonormal and

ܹ ᇱ is the transpose of ܹ . We similarly derive the spillage

of the uncontracted Gaussian functions

{ χ } as: i



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

(

1 s0 = 1 − Tr W% ′S% −1W% n

)

Page 12 of 41

(14)

Here, S% and W% are respectively the overlap matrices between primitive basis functions (

S%ij = χi χ j ), and between primitive basis functions and target wavefunctions ( Wij = χi ψ j ). Eqs. (13) and (14) are related via the contraction coefficients matrix C (see Eq. (6)), so % ' . We evaluate s via Eq. (14) by numerically evaluating S% and W% , W = CW% and S = C SC 0

(

)

computing W% ′S% −1W% , and evaluating s0 from the trace of the latter matrix. The subspace spanned by a contracted basis for a fixed set of exponents is always smaller than that spanned by the uncontracted basis. A larger subspace will give rise to lower spillage. The spillage of the uncontracted basis s0 therefore always bounds that of any contracted version and is the best possible spillage for a given set of exponents. We thus can optimize the GTO basis by optimizing first the exponents to minimize s0 (vide infra) and then the matrix C to minimize the difference between s and s0 . We do not directly minimize the spillage value s for the second step; we instead employ a PCA method to search for the optimal coefficient matrix C that best represents the uncontracted basis set, as described next.

2.2.3 PCA-based optimization PCA is a well-established method for dimensionality reduction,18 employed here to derive an expression for the contraction coefficient matrix C for a given set of exponents. We briefly summarize the PCA scheme and how it applies to our problem. We first outline how to obtain the coefficients of the primitive basis functions to fit the target wavefunctions. We denote these coefficients by X% , a g 0 by n matrix with g 0 the number of primitive Gaussian functions and n the number of wavefunctions. It can be proven that the 12 ACS Paragon Plus Environment

Page 13 of 41

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


least squares estimate of X% can be computed via the ordinary least squares formula

(

X% = S% ' S%

)

−1

S% 'W% . We then solve for the contraction coefficient matrix C using X% . C is a g 0

by g matrix with g the number of contracted Gaussian functions. We generate a set of target wavefunctions for a given element (vide infra), and the number of contracted Gaussian functions is usually much smaller than the number of target wavefunctions, i.e., g < n . We thus convert the optimization problem to a search for a g 0 by g matrix C that best represents the g 0 by n matrix X% , with the constraint that g < n . This is a typical problem solvable by a standard

% % ' and performing a singular value PCA,18 by evaluating the covariance matrix A = XX decomposition of A : V − 1 DV = A , with D = diag ({di } ) , and ݀௜ ∈ ℝା the singular values of A. The first g 0 rows of V comprise the matrix C if the di are sorted in descending order. For each element we consider, we prepare a set of PW-DFT one-electron wavefunctions for homonuclear diatomics of that element over a range of bond lengths. We then evaluate s0 for a set of primitive Gaussians with given exponents at each optimization step. Data from all the bond lengths are included in the evaluation of s0 simultaneously, with equal weights. We next conduct a brute force scan to search for the optimal exponents that minimize s0 . In practice, we always found that the s0 value for the exponents of the existing basis set was sufficiently low (vide infra) to render further optimization of exponents unnecessary. We thus simply obtain the coefficient matrix C via PCA, using the exponents from the original primitive Gaussian basis set.



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 41

2.3 Numerical details We employ the ABINIT 7.0.5 code12a, 19 with the Perdew-Burke-Ernzerhof (PBE)20 XC functional to perform PW-basis, spin-polarized DFT calculations. PW-based OEP calculations are carried out using a modified ABINIT 7.0.5 code (“invKS-ABINIT”)5, 12a, 19 with the same XC functional.5 We use an OPIUM-generated21 soft Hartree-Fock (HF)-based pseudopotential for O to achieve reasonable PW kinetic energy cutoffs and Trail-Needs (TN)-HF-based pseudopotentials for all other elements.22 The ECP form22 of each of these pseudopotentials (see supporting information) is employed in GTO-DFT-PBE calculations, within the MOLCAS 7.8 code.9a, 23 Diatomic molecules are used to break the spherical symmetry of atoms so as to provide further basis set flexibility when optimizing the GTO basis set. We generate target wavefunctions for ground-state Al, H, C, O, Mg, and Na dimers at varying bond lengths (Table 1) via spin-polarized PW-DFT-PBE calculations for ground state triplet Al2, singlet H2, singlet C2, triplet O2, singlet Mg2, and singlet Na2, each placed in a 22 Å × 16 Å × 16 Å supercell, using

Ecut s of 1000 eV for all elements except O, which requires 1400 eV. The bond length values are chosen around the equilibrium distance for each dimer. We include all occupied orbitals (and the first unoccupied p-type orbitals for the metals that we consider) as target wavefunctions to achieve a more effective basis set optimization (vide infra). Table 1. Dimer bond lengths for generating target wavefunctions for Al, H, C, O, Mg, and Na basis set optimization. Dimers Al2 H2 C2 O2 Mg2

Bond lengths for generating target wavefunctions (Å) 2.0, 2.5, 3.0, 3.75, 4.5 0.6, 0.75, 0.9, 1.2, 1.5 1.0, 1.25, 1.5, 2.0, 3.0 1.0, 1.208, 1.5, 2.0, 3.0 2.125, 2.375, 2.875, 3.375 14 ACS Paragon Plus Environment

Page 15 of 41

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


Na2

2.0, 2.5, 3.0, 3.75, 4.5

Test cases to evaluate performance of the optimized basis sets included singlet states of a linear chain of four Al atoms, H2O, CO, and NaH, each placed in a 22 Å × 10 Å × 10 Å supercell. Tests for Mg are given in the supporting information. Al-Al bond lengths were taken to be 3.0 Å for the linear Al4 chain. DFT-PBE optimized geometries were used for H2O and CO, while NaH’s experimental bond length of 1.88 Å was used.24 An Ecut of 1000 eV was used for Al4 and NaH, and 1400 eV for H2O and CO. The corresponding GTO-DFT calculations used the same geometries as the PW-DFT ones. We compare densities for the Al4 chain, H2O, and CO, while test OEP calculations are done for the Al4 chain and NaH. For the latter molecule, we employed HF25 and complete-active-space self-consistent field (CASSCF)26 densities for testing beyond DFT. We perform DFET calculations in a PW basis using a modified ABINIT code (“NonScABINIT”).2i, 12a, 19 For PFET calculations, the PW-based embedded DFT calculations are carried out with a different modified ABINIT code5,

12a, 19

(“emb-ABINIT”), while the GTO-based

embedded calculations are performed with a modified MOLCAS code23 (“emb-MOLCAS”). The two modified software packages read in the embedding potential on a uniform real-space grid and output the electron density on the same grid. PW-based OEP calculations in PFET are performed with “invKS-ABINIT” as mentioned above,5, 12a, 19 which reads in a reference density on a uniform real-space grid and outputs the optimized embedding potential, as well as the kinetic potential/kinetic energy on the same grid. Benchmark DFT results for the entire system are all from GTO-DFT calculations within the MOLCAS code,23 except for the H2O adsorption on an MgO cluster example discussed in the supporting information, where also ABINIT PWDFT results are presented. DFET, PFET, and benchmark DFT calculations are done with the 15 ACS Paragon Plus Environment


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 41

PBE XC functional.20 We use an Ecut of 1400 eV for all PW-based tests, to converge total energies to within 10 meV/atom. We use Dunning-style optimized VTZ basis sets (vide infra) for all atoms in the GTO-DFT calculations encompassing the entire system. The main goal of part I is to assess the performance of PFET with mixed basis sets; to do so we must first verify the fidelity of the pseudopotentials/ECPs and optimized GTO basis sets used in this mixed representation. As is commonly done (without rigorous justification), we perform DFT test calculations using HF-based ECPs, in anticipation of CW calculations that properly build on HF theory. Each atom is treated with the same ECP/pseudopotential appropriate for the respective GTO/PW basis set to keep physical consistency throughout. In this way, setups of PW-based KS-DFT, GTO-based KS-DFT, GTO-DFT/GTO-DFT PFET, and GTO-DFT/PW-DFT PFET or DFET are consistent with each other, thereby simplifying our analysis of PFET. We do not recommend using HF-based pseudopotentials for DFT calculations, but we employ them here simply for method verification. The comparisons between different levels of theory are meaningful as long as the PFET DFT/DFT calculations are set up consistently with the corresponding benchmark KS-DFT calculations. Note that the PFET formalism is general in nature, independent of the nature of the subsystems or the basis sets used for the subsystems. Although our test cases are all performed with mixed basis sets, the environment does not have to be an extended periodic one. Using Gaussian basis sets for the environment is also an option. We label henceforth each embedding calculation “((basis set type in subsystem 1)-(level of theory in subsystem 1)) /((basis set type in subsystem 2)-(level of theory in subsystem 2)) embedding scheme”. For example, GTO-DFT/PW-DFT PFET represents a PFET calculation in which subsystem 1 is described with a GTO basis set at the DFT level and subsystem 2 treated with a PW basis set at the DFT level.


Page 17 of 41

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


In the Al12 and (CO)Al12 test calculations, we place the Al12 cluster in the center of a 14.3 Å × 14.3 Å × 17 Å hexagonal supercell. We carve the cluster out of a face-centered cubic (fcc) Al crystal to mimic the Al (111) surface, with a lattice constant of 4.05 Å (the equilibrium lattice constant from periodic DFT-PBE calculations with the HF-based pseudopotentials). We fix the CO bond length at 1.134 Å, the GTO-DFT-PBE equilibrium value. We partition the bare Al12 cluster into a central Al6 cluster and an environment with six Al atoms, denoted as Al6/Al6. Note that the choice of Al6/Al6 as the test case mitigates complexity introduced by the basis. For this system, we can perform both pure-GTO and pure-PW calculations on the entire system. These are the benchmarks for the subsequent analysis of mixed basis sets. Only then can we understand what differences are caused by the basis set choices versus other numerical/algorithmic choices. Such a comparison would be difficult if we use an extended environment, as we do not have the means to perform a consistent, full-scale GTO calculation with periodic boundary conditions. When we include the CO adsorbate, the central Al6 cluster and the CO molecule comprise one subsystem and the other six Al atoms the other. The total system geometry is thus denoted as (CO)Al6/Al6. We suppress unphysical oscillations in the final embedding potential with a penalty

(

)

(

)

v function of λ ∫ Vemb∇2Vemb dr ( λ = 10−5 hartree-1·bohr-7) in the DFET calculations and a penalty v function of λ ∫ Vemb∇2Vemb dr ( λ = 10−4 hartree-1·bohr-7) for the PFET calculations. The WuYang functional value change is required to be smaller than 1.0×10-5 hartree for convergence in both schemes. We use electron densities and total energies as two metrics to evaluate the quality of the embedding potential. For each Vemb , we evaluate the electron density as:



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

ρ [Vemb ] = ρ I [Vemb ] + ρ II [Vemb ]

Page 18 of 41

(15)

and use the PFET energy expression: DFT Etot [Vemb ] = EIDFT [Vemb ] + EIIDFT [Vemb ] + Eint [ ρ I , ρ II ]

(16)

to obtain the total energy. Note that for embedding with a GTO-based subsystem (e.g., GTODFT/GTO-DFT or GTO-DFT/PW-DFT), the GTO basis set is explicitly used to represent the subsystem density in Eq. (16). A PW-based OEP calculation with the GTO-derived density provides the non-additive kinetic energy in Eint [ ρ I , ρ II ] (Eq. (16)). We denote total energies evaluated this way as PFET energies with a designated embedding potential (vide infra).

3 Results and discussion 3.1 ECP construction and tests As discussed in the supporting information (Sections 1.1 and 1.2), the TN ECPs/pseudopotentials used here are soft enough for PW-based calculations for all elements considered except oxygen. These ECPs/pseudopotentials have been used extensively in quantum Monte Carlo CW calculations, hence they should be valid for use here. The DFT-PBE total energy difference for a 3P O atom in a 22 Å × 10 Å × 10 Å supercell when using kinetic energy cutoffs Ecut of 1600 versus 1800 eV is 0.36 eV with the TN22 pseudopotential, whereas this difference is only 0.005 eV for Ecut of 1400 versus 1600 eV using a soft HF-based pseudopotential derived from the OPIUM21 code. We therefore employ the soft HF-based pseudopotential for O henceforth. We use the OPIUM21 code to map the soft HF-based pseudopotential onto the ECP form, as discussed in the supporting information. We take the {n k } powers of r (see Eqs. SI.8 and SI.9) from the TN-ECP for O27 and add one more Gaussian function with nk = 5 to channels l 18 ACS Paragon Plus Environment

Page 19 of 41

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


=1 and l =2. This initial guess for {n k } had a smaller starting residual before fitting than the original {n k } taken from Ref. 27. We set the exponents {ξ k } in a range of 5 to 40. We optimize the parameters such that the residuals defined in Eqs. (SI.8) and (SI.9) are below 10-2 a.u. for each channel and the pseudo-wavefunction-based error function (Eq. (SI.11)) is below 10-5 a.u. We then examine how well the fitted ECP form matches the target numerical pseudopotential (Figure 1(a)) and likewise for the pseudo wavefunctions corresponding to the two potentials (Figure 1 (b)). We find almost perfect agreement for all target quantities. Likewise, comparison with benchmark all-electron calculations and experimental data show excellent agreement, further corroborating that our ECPs are well suited for future CW calculations (see supporting information, Section 2).

Figure 1. Oxygen atom (a) target pseudopotentials and fitted ECPs and (b) pseudo wavefunctions for different angular momentum channels. For each color, the light solid line represents the target, while the dark dashed line represents the fitting result. The targets and fitting results are almost on top of each other.



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 41

3.2 Results of GTO basis set optimization Al is used to illustrate the GTO basis set optimization procedure; we then apply the optimization protocol to H, C, O, Mg, and Na, and evaluate performance of the optimized GTO bases in a number of applications. The initial GTO guesses are Burkatzki et al.’s QMC-VXZ bases, X being the cardinal number of the basis set.28 Their general form follows the Dunning style.13g For example, the QMC-VTZ basis set for Al is 11s11p2d1f contracted to 3s3p2d1f. Basis functions for s and p channels are similar in structure, with a contracted coefficient matrix as schematically represented in Figure 2(a). Gray blocks denote non-zero coefficients while white blocks denote zeroes. Each column in the matrix represents the coefficients for a particular contraction. In their original construction,28 exponents of primitive Gaussians were optimized to minimize the coupled cluster with single, double, and perturbative triples excitations total energy for the atom, while contraction coefficients were set by the HF reference state’s orbital coefficients. We only optimized parameters for s and p basis functions because their shapes dominate the electron density profiles for main group elements considered here.


Page 21 of 41

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


Figure 2. Schematic diagram of an 11 × 3 contraction coefficient matrix for the s and p basis functions in a QMC-VTZ basis set: (a) the original Dunning style QMC-VTZ basis set and (b) the fully optimized (f-opt) QMC-VTZ basis set. Gray blocks denote non-zero values in the matrix; white blocks denote zeroes. We first examine the extent to which the spillage s0 depends on the exponents. The original Al basis has exponents defined by

{a

0, p

{a

0, s

= 0.0455180, q s = 2.2037} for the s and

= 0.0309670, q p = 2.0856} for the p channel; see Table 2 and Table 3 for definitions. We

do a brute force scan of those four parameters to roughly locate the optimal exponents that minimize ‫ݏ‬଴ . Target wavefunctions used to compute the spillage came from the five Al2 configurations mentioned earlier. A combination of small a and q values means a narrow span of exponents, leading to a high ‫ݏ‬଴ value. At the other extreme, large a and q values lack sufficient contributions from diffuse functions to properly describe wavefunctions in the valence region, also leading to a high

‫ݏ‬଴ value. Within a certain range (denoted in bold in Table 2 and Table 3), the values for a and q 21 ACS Paragon Plus Environment


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 41

are balanced, with ‫ݏ‬଴ less sensitive to the parameters. The original basis set values

{a

0, s

= 0.0455180, q s = 2.2037} and {a0, p = 0.0309670, q p = 2.0856} , lie just in this region and

almost hit the minimal spillage obtainable with this basis set size. Thus the exponents as provided in the original QMC-VTZ basis set are suitable as is. We next proceed to optimize the coefficients, employing the original exponents. Table 2. Spillage s0 for Al2 at the DFT-PBE level with different combinations of {a0,s , qs } (sexponents of the Gaussians are α k = a0,s qsk −1 , k = 1, 2,L , 9 ). Values in bold are optimal.

a0,s = 0.0015

a0,s = 0.005

a0,s = 0.015

a0,s = 0.05

a0,s = 0.15

qs = 1.1

0.12

0.053

0.011

0.0064

0.0030

qs = 1.5

0.020

0.0097

0.0031

0.0029

0.0047

qs = 2.0

0.0040

0.0029

0.0029

0.0031

0.0078

qs = 3.0

0.0030

0.0031

0.0032

0.0033

0.011

{

}

Table 3. Spillage s0 for Al2 at the DFT-PBE level with different combinations of a0,p , q p (pexponents of the Gaussians are α k = a0, p q kp −1 , k = 1, 2,L , 9 ). Values in bold are optimal.

a0,p = 0.0015

a0,p = 0.005

a0,p = 0.015

a0,p = 0.05

a0,p = 0.15

qp = 1.1

0.32

0.18

0.017

0.0031

0.0062

qp = 1.5

0.035

0.0033

0.0030

0.0030

0.015

qp = 2.0

0.0036

0.0029

0.0029

0.0031

0.032

qp = 3.0

0.0030

0.0030

0.0031

0.0031

0.045

One straightforward way to optimize an “11 contracted 3” basis set is to first perform the PCA and then select the first three columns (principal components) of the matrix V as the contraction coefficients. The optimized coefficient matrix is comprised of the three most dominant vectors spanning the target wavefunctions and is an 11 × 3 matrix with all non-zero values (Figure 2(b)). We denote this fully optimized basis set as the f-opt QMC-VTZ basis. 22 ACS Paragon Plus Environment

Page 23 of 41

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


Another approach is to optimize the coefficients just of the first contraction with PCA, using a 9×3 matrix, and then leave the other two exponents the same as in the original QMC-VTZ basis set. This optimized basis set retains the Dunning style as in the original QMC-VTZ basis (Figure 2(a)). We denote this Dunning-style optimized basis set as the D-opt QMC-VTZ basis. We next obtain PW-DFT wavefunctions for spin-polarized triplet Al dimers. There are six occupied orbitals for each dimer, four s-like and two p-like, at five bond lengths. Having only ten p-like wavefunctions is insufficient for an effective optimization of the parameters for the p basis function. We therefore also include the unoccupied p wavefunctions as target wavefunctions. We follow the above two PCA schemes to optimize the basis set coefficients and evaluate the spillage for those new bases as well as the original and uncontracted QMC-VTZ bases for comparison. These spillage values are given for each dimer configuration in Table 4.

Table 4. Spillage of Al GTO basis sets for different coefficient optimization schemes, for Al dimers with five varying bond lengths. Al-Al bond length (Å) Optimization scheme original QMC-VTZ uncontracted QMC-VTZ f-opt QMC-VTZ D-opt QMC-VTZ

2.00

2.50

3.00

3.75

4.50

7.1×10-3 5.6×10-3 6.5×10-3 6.8×10-3

4.3×10-3 3.5×10-3 3.8×10-3 4.1×10-3

4.1×10-3 3.3×10-3 3.7×10-3 3.9×10-3

1.8×10-3 1.3×10-3 1.7×10-3 1.6×10-3

9.7×10-4 6.9×10-4 1.1×10-3 7.9×10-4

The PCA-optimized basis sets generally decrease the spillage when compared to the original basis set. The only exception is the f-opt QMC-VTZ result for the most stretched Al-Al bond. The optimization scheme is designed to improve basis set performance in terms of decreasing overall spillage across dimer configurations; a decrease in spillage for every configuration is not guaranteed. The spillage for the largest bond length (4.50 Å) is already



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 41

smaller than for all other configurations and thus is not dominant in the optimization procedure. The f-opt QMC-VTZ basis yields smaller spillage at shorter Al-Al bond lengths whereas the Dopt QMC-VTZ basis performs better at longer ones. Averaging over the five configurations, the f-opt QMC-VTZ leads to a lower average spillage of 3.36 × 10-3 than the D-opt QMC-VTZ basis set’s value, 3.44 × 10-3, though the difference is very small. The data in Table 4 suggest that the variance in spillage produced by different GTO bases is not very significant. We therefore need another metric in addition to spillage to effectively assess basis set performance. We take an Al4 chain as an example and directly compare the densities obtained with different GTO basis sets. We take the PW density as a reference to plot density deviations along the Al-Al bond (Figure 3). Optimization of the GTO basis set successfully alleviates much of the deviation from the PW density. The maximum absolute density difference drops from 5×10-3 to 5×10-4 a.u. for either of the optimized basis sets. The density distribution greatly improves in both the core and bonding regions. The improvements corresponding to the two optimization schemes are almost equivalent. The first column extracted from the matrix V in the PCA procedure of the Dunning-style optimization covers a large portion of the subspace spanned by the first nine primitive Gaussian functions. We form the Dopt QMC-VTZ basis set by adding two more primitive Gaussian functions (the most diffuse of the original basis) to this extracted column. The resulting D-opt QMC-VTZ basis set therefore covers the majority of the space spanned by the overall 11 primitive Gaussian functions. This explains why the D-opt QMC-VTZ and f-opt QMC-VTZ basis sets exhibit similar improvements in the spillage. The oscillation observed in the density difference profile is inevitable, however, because of the intrinsic difference between the GTO and PW basis sets. However, both PCAbased optimization schemes work well to reduce the density discrepancy with respect to the PW


Page 25 of 41

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


density. To retain the flexibility of the GTO basis set for further possible optimization and to guarantee CW quality, we employ the Dunning-style optimization procedure so that the exponents of the two diffuse primitive Gaussians can be optimized independently following the “correlation-consistent” approach proposed by Dunning and coworkers.13g

The agreement

between our optimized basis sets and available Dunning-style basis sets for CW calculations is satisfactory (see supporting information, Section 2).

Figure 3. Density deviation of an Al4 chain along the bond axis, evaluated as the difference between the optimized GTO and PW basis densities. The aluminum nuclei are represented by gray spheres located at z = 6.5, 9.5, 12.5, and 14.5 Å. We optimize basis sets for O, H, C, Mg, and Na following the same procedure as described above for Al, but this time using only the Dunning-style optimization. In each case, the resulting average spillage from the optimized bases (Table 5) was lower than for the original bases, and are properly bounded from below by the fully uncontracted bases. Table 5. Average spillage for O, H, C, Mg, and Na dimers with different GTO basis sets.

original QMC-VTZ

O 2.2×10-3

H 1.3×10-3


C 1.3×10-3

Mg 5.8×10-2

Na 3.5×10-3


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

uncontracted QMC-VTZ D-opt QMC-VTZ

9.6×10-4 1.3×10-3

1.2×10-5 5.4×10-4

8.1×10-4 1.1×10-3

Page 26 of 41

5.0×10-3 1.1×10-2

1.4×10-3 1.8×10-3

We next examine the densities of H2O and CO molecules to test transferability of the optimized basis sets. We evaluate the electron density on a 3D uniform grid and compare the density at each grid point by plotting the GTO-based density versus the PW-based density (Figure 4). If the density produced by the two basis sets matches perfectly, such a plot should result in a straight line with a slope of 1. Improvement after optimization is evident in both cases: the data points cluster around the optimal line when using the optimized basis sets.

Figure 4. Comparison between GTO- and PW-derived densities for (a) H2O and (b) CO. The horizontal axis has the PW-derived density value at each grid point, while the vertical axis gives the GTO-derived density at the same grid point. The green dashed line represents the ideal case of perfect agreement (x=y).

In summary, test calculations show that a simple PCA-based optimization is effective at decreasing the spillage and at improving agreement between GTO- and PW-derived densities. 26 ACS Paragon Plus Environment

Page 27 of 41

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


Performance tests of the optimized basis sets when used in CW calculations of atomic and dimer properties can be found in the supporting information, Section 2.

3.3 OEP calculations with optimized GTO basis sets Beyond comparisons of the GTO-derived density, we must test whether such densities behave well in PW-based OEP calculations without introducing severe numerical noise. This metric is more practical and important since Voep [ ρ ] is the quantity directly involved in optimization of the embedding potential. At the complete basis set limit, an OEP calculation with a given DFT density in theory should yield a Voep [ ρ ] identical to the KS effective potential ( VKS [ ρ ] ). VKS [ ρ ] can be evaluated as a functional of the reference density by summing up the Coulomb, XC, and ionic potentials. We use this metric to examine whether the optimized GTO basis sets are sufficiently compatible to ensure a well-behaved OEP calculation. We use the DFT densities of an Al4 chain evaluated with the original and the optimized GTO basis sets as reference densities for PW-based OEP calculations. We compare VKS [ ρ ] and Voep [ ρ ] along the Al-Al bond axis (Figure 5). Use of the original basis set density as input to the

OEP algorithm shows a significant discrepancy between the two potential profiles whereas the two potentials evaluated with the optimized basis set agree very well. Minimization of the difference between the two potentials is essential, as the self-consistency cycle in PFET will

(

)

amplify any errors induced by inconsistencies between the basis sets. Voep [ ρ ] − VKS [ ρ ] should be zero when all subsystems are treated with DFT and indeed Voep [ ρ ] does not deviate from VKS [ ρ ] at the DFT level when using the optimized GTO basis set. We therefore can expect

(V [ ρ ] − V [ ρ ]) oep

KS

with the OEP derived from a CW density will nicely capture only the



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 41

physical difference between the CW and DFT descriptions of the subsystem, excluding any unphysical artifacts induced by basis set incompatibility, as desired.

Figure 5. Profiles of KS effective potentials VKS [ ρ ] and OEPs Voep [ ρ ] for linear Al4, employing (a) QMC-VTZ and (b) Dunning-style optimized QMC-VTZ basis sets. (a) Black and (b) red solid lines are VKS [ ρ ] ; (a) green and (b) blue dashed lines are Voep [ ρ ] . The aluminum nuclei are represented by gray spheres located at z = 6.5, 9.5, 12.5, and 14.5 Å. A more stringent test of basis set compatibility is to perform the OEP calculation with a non-DFT density. HF and CASSCF densities were computed for NaH to consider an alternative (ionic) bonding scenario. The complete valence active space of the CASSCF calculation consists of two electrons in five orbitals (Na 3s3p and H 1s). Figure 6(a) displays the difference between the non-DFT and DFT reference densities, as well as the analogous density differences produced non − DFT DFT by the OEP. The perfect match between ( ρ ref − ρ ref ) and ( ρ ononep − DFT − ρ oeDFp T ) proves that the

PW-based OEP algorithm is sufficiently robust to employ non-DFT GTO-based densities in its

(

)

non − DFT  − Voep  ρreDFT  KS equation inversion scheme. Figure 6(b) also illustrates that Voep  ρref f  is


Page 29 of 41

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


negatively correlated with the density difference

(ρ

non − DFT ref

DFT − ρ ref

)

despite its oscillatory

behavior. Such a negative correlation between the two quantities is physical, as a more positive potential prevents electron accumulation, while a potential well attracts electrons, yielding a higher electron density.

CASSCF DFT − PBE Figure 6. (a) Density differences along the Na-H bond. Black solid line is ( ρ ref − ρ ref ) CASSCF DFT − PBE and green dashed line is ( ρ oep − ρ oep ) ; red solid line is

line is

(ρ

HF oep

(ρ

HF ref

DFT − PBE − ρ ref ) and blue dashed

− PBE − ρ oeDFT ) . (b) OEP differences along the Na-H bond. Black curve is p

CASSCF DFT − PBE HF DFT − PBE (Voep − Voep ) and red curve is (Voep − Voep ).

3.4 Full PFET with a mixed basis at the DFT/DFT level We carry out PFET calculations on an Al6/Al6 cluster to investigate PFET’s performance for metallic bonding and then adsorb CO on this cluster to test PFET for donor-acceptor bonding.



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 30 of 41

We assess performance of GTO-DFT/PW-DFT PFET with regard to reproducing the total electron density and evaluating the total energy. We first examine the topology of the converged PFET embedding potential (Figure 7). The embedding potential has wells in the interface region, promoting metallic bonding, and is positive around the Al nuclei, as expected for a valence electron distribution (see Figure 7(c) and (d) for red spots coinciding with Al atom locations).

Figure 7. (a) Geometry of Al6/Al6 test case. (b) Isosurfaces at ±0.041 a.u. of embedding potential Vemb { PFET , ( GTO − DFT ) / ( PW − DFT ) } (yellow represents positive values, blue negative).

Contour plots of the embedding potential crossing the first (c) and second (d) atomic layers. Density deviations in the first atomic layer produced by different embedding potentials, relative to the ground-state PW-DFT density of the entire system, are displayed in Figure 8. All embedding potentials are applied within a GTO-DFT/PW-DFT embedding framework to evaluate the density. Recall that the DFET and PFET PW-based embedding potentials are formally equivalent at the DFT/DFT level, despite being numerically different due to their 30 ACS Paragon Plus Environment

Page 31 of 41

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


distinct evaluation methods (see supporting information, Section 3.1 for a detailed discussion). As the algorithm for its optimization is simpler and leads to less numerical noise than the OEP algorithm in PFET, we take the DFET PW-based embedding potential as the embedding potential benchmark. We thus compare results using a converged PW basis in each subsystem and the embedding potential Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} to all other schemes. We first use Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} to derive the density at the GTODFT/PW-DFT level (Figure 8(a)). The resulting small density deviation shows that the DFET embedding potential constructed with PW-based calculations is transferable to the mixed basis case, with GTOs in the embedded region and PWs elsewhere. By contrast, the sum of the isolated subsystem densities (with GTOs in the embedded region, PWs in the environment), i.e., when Vemb = 0 is applied, exhibits large density deviations in the interface region (Figure 8(b)), exactly where an embedding potential is needed to restore correct behavior. Finally, the density deviation again almost vanishes everywhere with a fully optimized embedding potential Vemb { PFET , ( GTO − DFT ) / ( PW − DFT )} [Figure 8(c)]. The similarity between Figure 8 (a)

and (c) demonstrates that construction of embedding potentials with PFET at the GTO-DFT/PWDFT level is nearly as accurate as that of Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} in terms of reproducing the density, and that use of a mixed basis set compared to all PWs does not introduce numerical noise in the density.



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 32 of 41

Figure 8. Contour plots of Al6/Al6 GTO-DFT/PW-DFT density deviations across the first atomic layer derived from (a) Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} , (b) Vemb = 0 , and (c) Vemb {PFET , ( GTO − DFT ) / ( PW − DFT )} . Circles represent the Al nuclei. Deviations are

relative to a PW-DFT computation on the entire system. We next consider how accurately PFET evaluates energetics. For this purpose, we adsorb a CO molecule perpendicularly onto the Al12 cluster with the C atom closest to the metal. We examine two types of relative energy curves: binding of vertically oriented CO molecule onto the fcc hollow site and then moving this CO molecule parallel to the surface along the 112  direction with dAl-C fixed at 1.400 Å and dC-O fixed at 1.134 Å (Figure 9). We use the embedding potential Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} constructed for the Al6/Al6 without the CO 0 adsorbate (denoted as Vemb ) as the PFET initial guess for all adsorption geometries, to accelerate

convergence.


Page 33 of 41

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


0 Figure 10 displays the PFET energy curves produced by the initial guess Vemb , the

optimized

embedding

potential

Vemb {PFET , ( GTO − DFT ) / ( PW − DFT )}

,

and

Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} calculated for each adsorption geometry, as well as the energy curve for the GTO-DFT benchmark on the entire system. All embedding potentials are applied within a GTO-DFT/PW-DFT embedding framework to derive the PFET energies.

Figure 9. Geometry of CO adsorbed on the Al12 cluster. The two crosses represent the fcc and bridge sites on the Al(111) surface.



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 34 of 41

Figure 10. PFET energy curves for (a) CO binding perpendicularly to the fcc hollow on the Al12 cluster and (b) parallel motion of perpendicularly adsorbed CO along the 112  direction on the Al12 cluster (Figure 9). Black is the GTO-DFT benchmark for the entire system; blue is the PFET 0 energy with fixed Vemb ( Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} constructed for Al6/Al6); red

is the PFET energy with Vemb {PFET , ( GTO − DFT ) / ( PW − DFT )} optimized at each geometry; and green is the PFET energy with Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} optimized at each geometry. The PFET energy curves produced by Vemb { DFET , ( PW − DFT ) / ( PW − DFT )} are in excellent agreement with the GTO-DFT benchmark (green versus black curves in Figure 10). The good agreement proves that the subsystem GTO basis is sufficient, with no need to employ a system-wide one (see supporting information, Section 3.2 for a detailed definition). Further, the GTO basis behaves well within the PW-based OEP procedure and does not induce significant numerical error in the energy evaluation. Thus, the PW-derived embedding potential is


Page 35 of 41

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


transferable to the embedding framework with mixed basis sets in terms of accurate energy evaluation. 0 PFET using a fixed Vemb captures the general shape of the (CO)Al12 binding energy curve

(blue curve in Figure 10(a)) but overbinds by up to 0.2 eV. Upon further optimization of the embedding potential with PFET, the energy curve (red curve in Figure 10(a)) agrees with the GTO-DFT benchmark (black curve in Figure 10(a)). For CO motion parallel to the cluster 0 surface, the validity of a fixed Vemb depends on the position of the CO molecule. When the

molecule is above the fcc site at the center of the cluster, it only induces a local density distortion within the Al6 cluster. In this case, the embedding potential constructed for Al6/Al6 is still sufficiently accurate to describe the interaction between (CO)Al6 and the environmental Al atoms. However, errors accumulate when the molecule approaches and moves beyond the bridge site (see blue versus black curves in Figure 10 (b)). The origin is a substantial density distortion at the interface between the central Al6 cluster and its Al environment induced by the CO molecule when it is close to the boundary of the Al6 cluster at the bridge site. The fixed embedding potential constructed for the isolated Al6/Al6 system then no longer appropriately describes the interaction therein. Only further DFET or PFET optimization of the embedding potential in the presence of the CO molecule, followed by evaluation of the PFET energy, successfully produces good agreement with the GTO-DFT benchmark (see red/green versus black curves in Figure 10 (b)). These results demonstrate the good performance of the PFET embedding potential constructed with a mixed basis set in terms of energy evaluation. Moreover, the results here also suggest that use of the more numerically robust DFET scheme to optimize the embedding potential followed by evaluation of the total energy within the PFET framework could provide an optimal balance of computational efficiency and accuracy. Although in CW-in35 ACS Paragon Plus Environment


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 36 of 41

DFT embedding, such an approach forfeits the advantage of PFET to self-consistently treat the polarization of the environment induced by changes in the embedded region, accuracy is still improved over DFET alone, as well as over a pure DFT calculation. Part II will consider this point in more detail.

4 Conclusions Because of the contrasts in GTO-based and PW-based electronic structure calculations, we constructed optimized GTO basis sets and used consistent electron-ion potentials to achieve compatibility between the two types of calculations. We used PCA to optimize GTO basis sets in order to remove discrepancies between GTO- and PW-derived densities. Densities derived from the optimized GTO basis sets better match the PW-derived electron densities than densities derived using the original GTO basis sets. Inverting the KS equations using a non-DFT GTO density within a PW-based code further corroborates that the optimized GTO basis sets yield the numerically stable OEP calculations needed for PFET. Test calculations at the CW level indicate that the quality of CW calculations is well-preserved with the optimized basis sets. These systematic tests allow us to conclude that our optimized GTO basis sets are valid and robust, and suitable for the OEP computations required in CW/DFT PFET calculations with mixed GTO/PW basis sets. Further tests proved that PFET now works with a purely PW basis set as well as with mixed GTO/PW basis sets, at the pure DFT level. PFET constructs a physically appropriate embedding potential through variational optimization. Because of the different algorithms in DFET and PFET, the uniqueness of the embedding potential cannot be numerically guaranteed, even when using only a PW basis set. This breakdown of formal uniqueness is a direct result of the nature of the PFET algorithm. First, the functional value ( Etot [Vemb ] ) is evaluated with a 36 ACS Paragon Plus Environment

Page 37 of 41

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


numerical OEP procedure; second, the evaluation of the functional derivative (

δ Etot [Vemb ] ) δ Vemb

employs a first-order finite-difference approximation. These two features decrease the numerical stability of PFET even when only a PW basis set is used. The intrinsic differences between PW and GTO basis sets introduce further numerical deviations in the embedding potential. However, test calculations demonstrate that basic properties such as densities and energies are not sensitive to these numerical discrepancies. These conclusions set a solid foundation to implement GTOCW/PW-DFT PFET for more complicated and interesting phenomena, as done in the companion article presented next.

Acknowledgements We are grateful to the National Science Foundation (Award No. 1265700) and the Department of Energy, Office of Science, Basic Energy Sciences (Award No. DE-SC0002120) for support of this work. The authors are pleased to acknowledge computational support from the TIGRESS high performance computer center at Princeton University, which is jointly administered by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology. We also thank Nari Baughman for her help in editing the paper.

Supporting Information The supporting information provides the details of pseudopotential and ECP construction, the CW results for atomic and diatomic properties with optimized GTO basis sets, and other PFET test cases at the DFT/DFT level. This information is available free of charge via the Internet at http://pubs.acs.org/.



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

References 1. (a) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesolowski, T. A., The merits of the frozendensity embedding scheme to model solvatochromic shifts. J Chem Phys 2005, 122, 094115; (b) Huang, P.; Carter, E. A., Advances in correlated electronic structure methods for solids, surfaces, and nanostructures. Annu. Rev. Phys. Chem. 2008, 59, 261-90; (c) Libisch, F.; Huang, C.; Carter, E. A., Embedded correlated wavefunction schemes: theory and applications. Acc. Chem. Res. 2014, 47, 276875; (d) Jacob, C. R.; Neugebauer, J., Subsystem density-functional theory. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 325-362; (e) Yu, K.; Krauter, C. M.; Dieterich, J. M.; Carter, E. A., Density and potential functional embedding: theory and practice. In Fragmentation: Toward Accurate Calculations on Complex Molecular Systems, Gordon, M., Ed. in press 2016. 2. (a) Wesolowski, T. A.; Warshel, A., Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050-8053; (b) Govind, N.; Wang, Y. A.; da Silva, A. J. R.; Carter, E. A., Accurate ab initio energetics of extended systems via explicit correlation embedded in a density functional environment. Chem. Phys. Lett. 1998, 295, 129-134; (c) Govind, N.; Wang, Y. A.; Carter, E. A., Electronic-structure calculations by first-principles density-based embedding of explicitly correlated systems. J Chem Phys 1999, 110, 7677; (d) Klüner, T.; Govind, N.; Wang, Y. A.; Carter, E. A., Periodic density functional embedding theory for complete active space self-consistent field and configuration interaction calculations: Ground and excited states. J Chem Phys 2002, 116, 42; (e) Jacob, C. R.; Wesolowski, T. A.; Visscher, L., Orbital-free embedding applied to the calculation of induced dipole moments in CO2...X (X = He, Ne, Ar, Kr, Xe, Hg) van der Waals complexes. J Chem Phys 2005, 123, 174104; (f) Huang, P.; Carter, E. A., Self-consistent embedding theory for locally correlated configuration interaction wave functions in condensed matter. J Chem Phys 2006, 125, 084102; (g) Sharifzadeh, S.; Huang, P.; Carter, E. A., All-electron embedded correlated wavefunction theory for condensed matter electronic structure. Chem. Phys. Lett. 2009, 470, 347-352; (h) Fux, S.; Jacob, C. 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; (i) Huang, C.; Pavone, M.; Carter, E. A., Quantum mechanical embedding theory based on a unique embedding potential. J Chem Phys 2011, 134, 154110; (j) Knizia, G.; Chan, G. K.-L., Density matrix embedding: a simple alternative to dynamical mean-field theory. Phys Rev Lett 2012, 109, 186404. 3. (a) Klüner, T.; Govind, N.; Wang, Y.; Carter, E., Prediction of electronic excited states of adsorbates on metal surfaces from first principles. Phys Rev Lett 2001, 86, 5954-5957; (b) Sharifzadeh, S.; Huang, P.; Carter, E., Embedded configuration interaction description of CO on Cu(111): resolution of the site preference conundrum. J. Phys. Chem. C 2008, 112, 4649-4657; (c) Libisch, F.; Huang, C.; Liao, P.; Pavone, M.; Carter, E. A., Origin of the energy barrier to chemical reactions of O2 on Al(111): evidence for charge transfer, not spin selection. Phys Rev Lett 2012, 109, 198303; (d) Mukherjee, S.; Libisch, F.; Large, N.; Neumann, O.; Brown, L. V.; Cheng, J.; Lassiter, J. B.; Carter, E. A.; Nordlander, P.; Halas, N. J., Hot electrons do the impossible: plasmon-induced dissociation of H2 on Au. Nano Lett. 2013, 13, 240-7; (e) Libisch, F.; Cheng, J.; Carter, E. A., Electron-transfer-induced dissociation of H2 on gold nanoparticles: excited-state potential energy surfaces via embedded correlated wavefunction theory. Z. Phys. Chem. 2013, 227, 1455; (f) Libisch, F.; Krauter, C. M.; Carter, E. A., Corrigendum to: Plasmon-driven dissociation of H2 on gold nanoclusters. Z. Phys. Chem. 2016, 230, 131-132; (g) Zhou, L.; Zhang, C.; McClain, M. J.; Manjavacas, A.; Krauter, C. M.; Tian, S.; Berg, F.; Everitt, H. O.; Carter, E. A.; Nordlander, P.; Halas, N. J., Aluminum nanocrystals as a plasmonic photocatalyst for hydrogen dissociation. Nano Lett. 2016, 16 (2), 1478-1484. 4. Elliott, P.; Burke, K.; Cohen, M. H.; Wasserman, A., Partition density-functional theory. Phys. Rev. A 2010, 82, 024501.


Page 38 of 41

Page 39 of 41

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


5. Huang, C.; Carter, E. A., Potential-functional embedding theory for molecules and materials. J Chem Phys 2011, 135, 194104. 6. (a) Cortona, P., Self-consistently determined properties of solids without band-structure calculations. Phys Rev B 1991, 44 (16), 8454-8458; (b) Cortona, P., Direct determination of selfconsistent total energies and charge densities of solids: A study of the cohesive properties of the alkali halides. Phys Rev B 1992, 46 (4), 2008-2014. 7. (a) Cohen, M. H.; Wasserman, A., On Hardness and Electronegativity Equalization in Chemical Reactivity Theory. J. Stat. Phys. 2006, 125 (5), 1121-1139; (b) Nafziger, J.; Wasserman, A., Density-Based Partitioning Methods for Ground-State Molecular Calculations. J. Phys. Chem. A 2014, 118 (36), 76237639. 8. (a) Wu, Q.; Yang, W., A direct optimization method for calculating density functionals and exchange–correlation potentials from electron densities. J Chem Phys 2003, 118, 2498; (b) Bulat, F. A.; Heaton-Burgess, T.; Cohen, A. J.; Yang, W., Optimized effective potentials from electron densities in finite basis sets. J Chem Phys 2007, 127, 174101. 9. (a) Karlström, G.; Lindh, R.; Malmqvist, P.-Å.; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P.-O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L., MOLCAS: a program package for computational chemistry. Comput. Mater. Sci. 2003, 28, 222-239; (b) Gordon, M. S.; Schmidt, M. W., Advances in electronic structure theory: GAMESS a decade later. In Theory and Applications of Computational Chemistry: the first forty years, Dykstra, C. E.; Frenking, G.; Kim, K. S.; Scuseria, G. E., Eds. Elsevier: Elsevier, 2005; pp 1167-1189. 10. (a) Städele, M.; Majewski, J. A.; Vogl, P.; Görling, A., Exact Kohn-Sham exchange potential in semiconductors. Phys Rev Lett 1997, 79 (11), 2089-2092; (b) Görling, A., Orbital- and state-dependent functionals in density-functional theory. J Chem Phys 2005, 123 (6), 062203. 11. Jacob, C. R., Unambiguous optimization of effective potentials in finite basis sets. J Chem Phys 2011, 135, 244102. 12. (a) Gonze, X.; Beuken, J.-M.; Caracas, R.; Detraux, F.; Fuchs, M.; Rignanese, G.-M.; Sindic, L.; Verstraete, M.; Zerah, G.; Jollet, F.; Torrent, M.; Roy, A.; Mikami, M.; Ghosez, P.; Raty, J.-Y.; Allan, D. C., First-principles computation of material properties: the ABINIT software project. Comput. Mat. Sci. 2002, 25 (3), 478 - 492; (b) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. 13. (a) Hehre, W. J.; Ditchfield, R.; Pople, J. A., Self-consistent molecular orbital methods. XII. further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules. J Chem Phys 1972, 56, 2257; (b) Dill, J. D.; Pople, J. A., Self-consistent molecular orbital methods. XV. Extended Gaussian-type basis sets for lithium, beryllium, and boron. J Chem Phys 1975, 62, 2921; (c) Francl, M. M.; Petro, W. J.; Hehre, W. J.; Binkely, J. S.; Gordon, M. S.; Defrees, D. J.; Pople, J. A., Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. J Chem Phys 1982, 77, 3654; (d) Dunning, T. H., Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J Chem Phys 1989, 90, 1007; (e) Widmark, P.-O.; Malmqvist, P.-Å.; Roos, B. O., Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions. Theor. Chim. Acta 1990, 77, 291-306; (f) Widmark, P.-O.; Persson, B. J.; Roos, B. O., Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions. Theor. Chim. Acta 1991, 79, 419-432; (g) Kendall, R. A.; Dunning, T. H.; Harrison, R. J., Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J Chem Phys 1992, 39 ACS Paragon Plus Environment


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

96, 6796; (h) Schäfer, A.; Horn, H.; Ahlrichs, R., Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J Chem Phys 1992, 97, 2571; (i) Woon, D. E.; Dunning, T. H., Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J Chem Phys 1993, 98, 1358; (j) Woon, D. E.; Dunning, T. H., Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties. J Chem Phys 1994, 100, 2975; (k) Schäfer, A.; Huber, C.; Ahlrichs, R., Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J Chem Phys 1994, 100, 5829; (l) Rassolov, V. A.; Pople, J. A.; Ratner, M. A.; Windus, T. L., 6-31G∗ basis set for atoms K through Zn. J Chem Phys 1998, 109, 1223; (m) Weigend, F.; Ahlrichs, R., Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297-305. 14. Lennard-Jones, J. E., The electronic structure of some diatomic molecules. Trans. Faraday Soc. 1929, 25, 668. 15. (a) Feller, D. F.; Ruedenberg, K., Systematic approach to extended even-tempered orbital bases for atomic and molecular calculations. Theor. Chim. Acta 1979, 52 (3), 231-251; (b) Schmidt, M. W.; Ruedenberg, K., Effective convergence to complete orbital bases and to the atomic Hartree–Fock limit through systematic sequences of Gaussian primitives. J Chem Phys 1979, 71 (10), 3951-3962; (c) Kutzelnigg, W.; Klopper, W., Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General theory. J Chem Phys 1991, 94 (3), 1985-2001. 16. Füsti-Molnar, L.; Pulay, P., Accurate molecular integrals and energies using combined plane wave and Gaussian basis sets in molecular electronic structure theory. J Chem Phys 2002, 116 (18), 7795-7805. 17. (a) Sánchez-Portal, D.; Artacho, E.; Soler, J. M., Projection of plane-wave calculations into atomic orbitals. Solid State Commun. 1995, 95, 685-690; (b) Sánchez-Portal, D.; Artacho, E.; Soler, J. M., Analysis of atomic orbital basis sets from the projection of plane-wave results. J. Phys.: Condens. Matter 1996, 8, 3859-3880; (c) Kenny, S. D.; Horsfield, A. P.; Fujitani, H., Transferable atomic-type orbital basis sets for solids. Phys Rev B 2000, 62 (8), 4899-4905; (d) Chen, M.; Guo, G.-C.; He, L., Systematically improvable optimized atomic basis sets for ab initio calculations. J. Phys.: Condens. Matter 2010, 22, 445501. 18. Wold, S.; Esbensen, K.; Geladi, P., Principal component analysis. Chemom. Intell. Lab. Syst. 1987, 2, 37-52. 19. (a) Gonze, X.; Rignanese, G.-M.; Verstraete, M.; Beuken, J.-M.; Pouillon, Y.; Pouillon, Y.; Caracas, R.; Jollet, F.; Torrent, M.; Zerah, G.; Mikami, M.; Ghosez, P.; Veithen, M.; Raty, J.-Y.; Olevano, V.; Bruneval, F.; Reining, L.; Godby, R.; Onida, G.; AllanX, D.; Hamann, R.; Allan, D. C., A brief introduction to the ABINIT software package. Z. Kristallogr. 2005, 220, 558-562; (b) Gonze, X.; Amadon, B.; Anglade, P.M.; Beuken, J.-M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Côté, M.; Deutsch, T.; Genovese, L.; Ghosez, P.; Giantomassi, M.; Goedecker, S.; Hamann, D. R.; Hermet, P.; Jollet, F.; Jomard, G.; Leroux, S.; Mancini, M.; Mazevet, S.; Oliveira, M. J. T.; Onida, G.; Pouillon, Y.; Rangel, T.; Rignanese, G.-M.; Sangalli, D.; Shaltaf, R.; Torrent, M.; Verstraete, M. J.; Zerah, G.; Zwanziger, J. W., ABINIT: Firstprinciples approach to material and nanosystem properties. Comput Phys Commun 2009, 180, 25822615. 20. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximation made simple. Phys Rev Lett 1996, 77, 3865-3868. http://opium.sourceforge.net/sci.html. (accessed January 4, 2017). 21. 22. (a) Trail, J. R.; Needs, R. J., Norm-conserving Hartree-Fock pseudopotentials and their asymptotic behavior. J Chem Phys 2005, 122, 14112; (b) Trail, J. R.; Needs, R. J., Smooth relativistic Hartree-Fock pseudopotentials for H to Ba and Lu to Hg. J Chem Phys 2005, 122, 174109. 23. Aquilante, F.; De Vico, L.; Ferre, N.; Ghigo, G.; Malmqvist, P.-A.; Neogrady, P.; Pedersen, T. B.; Pitovnak, M.; Reiher, M.; Roos, B. O.; Serrano-Andres, L.; Urban, M.; Veryazov, V.; Lindh, R., MOLCAS 7: the next generation. J. Comput. Chem. 2010, 31 (1), 224-247. 40 ACS Paragon Plus Environment

Page 40 of 41

Page 41 of 41

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


24. Huber, K. P.; Herzberg, G., Molecular Spectra and Molecular Structure, IV. Constants of Diatomic Molecules. 1979, 440. 25. Roothaan, C. C. J., New developments in molecular orbital theory. Rev. Mod. Phys. 1951, 23, 6989. 26. (a) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M., A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach. Chem. Phys. 1980, 48, 157-173; (b) Roos, B. O., The complete active space SCF method in a fock-matrix-based super-CI formulation. Int. J. Quantum Chem. 2009, 18, 175-189. https://vallico.net/casinoqmc/pplib/. (accessed January 4, 2017). 27. 28. Burkatzki, M.; Filippi, C.; Dolg, M., Energy-consistent pseudopotentials for quantum Monte Carlo calculations. J Chem Phys 2007, 126, 234105.

TOC graphic


Potential Functional Embedding Theory at the ... - ACS Publications

Recommend Documents