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 Wave Function Level, Part II: Error Sources and Performance Tests Jin Cheng, Kuang Yu, Florian Libisch, Johannes M. Dieterich, and Emily A. Carter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01011 • 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 39

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 Wave Function Level, Part II: Error Sources and Performance Tests Jin Cheng,1 Kuang Yu,2 Florian Libisch,2,3 Johannes M. Dieterich2 and Emily A. Carter4a) 1

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

4

School of Engineering and Applied Science, Princeton University, Princeton, New Jersey

08544 3

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

Austria, EU ABSTRACT Quantum mechanical embedding theories partition a complex system into multiple spatial regions that can use different electronic structure methods within each, in order to optimize tradeoffs between accuracy and cost. The present work incorporates accurate but expensive correlated wavefunction (CW) methods for a subsystem containing the phenomenon or feature of greatest interest, while selfconsistently capturing quantum effects of the surroundings using fast but less accurate density functional theory (DFT) approximations. We recently proposed two embedding methods [for a review, see Acc. Chem. Res. 47, 2768 (2014)]: density functional embedding theory (DFET) and potential functional embedding theory (PFET). DFET provides a fast but non-self-consistent density-based embedding scheme, while PFET offers a more rigorous theoretical framework to perform fully self-consistent, variational CW/DFT calculations [as defined in Part I, CW/DFT means subsystem 1(2) is treated with CW(DFT) methods]. When originally presented, PFET was only tested at the DFT/DFT level of theory as a proof of principle within a planewave (PW) basis. Part I of this two-part series demonstrated that PFET can be made to work well with mixed Gaussian type orbital (GTO)/PW bases, as long as optimized GTO bases and consistent electron-ion potentials are employed throughout. Here in part II we 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

Page 2 of 39

conduct the first PFET calculations at the CW/DFT level, and compare them to DFET and full CW benchmarks. We test the performance of PFET at the CW/DFT level for a variety of types of interactions (hydrogen bonding, metallic, and ionic). By introducing an intermediate CW/DFT embedding scheme denoted DFET/PFET, we show how PFET remedies different types of errors in DFET, serving as a more robust type of embedding theory. _____________________________ a)

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

1 Introduction Accurate electronic structure determination is a challenging yet foundational task in many chemistry and material science applications. A variety of electronic structure theories have been developed that achieve different tradeoffs between accuracy and computational cost. For example, various approximations to density functional theory (DFT) comprise the most common means to study quantum aspects of condensed matter,1 typically with periodic boundary conditions (PBCs) imposed.2 While DFT approximations (DFAs) are fairly efficient, unacceptably large errors emerge in a number of situations (e.g., dispersion interactions3, charge transfer reactions,4 excited states,4 etc.). Correlated wavefunction (CW) methods are more accurate than DFAs but the former’s high computational cost and scaling with system size hinder their applications,5 especially with PBCs and the attendant large planewave (PW) basis sets.6 The general idea of embedding theory is to partition a complex system into multiple subsystems and then treat each one at the physically most suitable level of theory.7 The embedded correlated wavefunction (ECW) method is one way to incorporate a CW-level embedded region (subsystem I) into a periodic DFT-level environment (subsystem II).8 ECW enables the use of an accurate yet computationally expensive CW method to capture the correct exchange-correlation (XC)

2


Page 3 of 39

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


interaction within the embedded region, while still appropriately accounting for the influence of the environment at the DFT level.9 The general formalism to evaluate the total energy of such a hybrid quantum mechanical (QM) system can be expressed as:

Etot = EICW + EIIDFT + Eint

(1)

Here, EICW is the CW energy of the embedded region, E IIDFT is the DFT energy of the environment, and

Eint is the interaction energy between the two subsystems, which can be evaluated using either orbitalfree DFT or an optimized effective potential (OEP) formalism.10 An effective embedding potential Vemb is introduced as a common practice to represent the interaction between the two subsystems, and the electronic structure of the subsystem is solved for in the presence of Vemb . Our group recently proposed two different embedding schemes: density functional embedding theory (DFET)11 and potential functional embedding theory (PFET).10 Both schemes assume a common embedding potential shared by all subsystems, which uniquely determines the density partitioning of the total system. Under this assumption, DFET aims to reproduce the density of the total system, while PFET aims to minimize the total energy. Theoretically, when using a consistent basis set, the two methods are equivalent and should generate an identical embedding potential at the DFT/DFT level (as defined in Part I, “method 1/method 2” respectively denotes the methods used in subsystem 1/subsystem 2). Compared to DFET, which computes the mutual polarization of the two subsystems at the DFT level, PFET provides a much more flexible theoretical framework to enable fully self-consistent hybrid CW/DFT calculations such that the interactions between the subsystems go beyond the DFT level of DFET. Up to now, PFET had only been verified at the PW-DFT/PW-DFT level of theory10 and its performance at the CW/DFT level with mixed basis sets had not been examined.

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

Page 4 of 39

Successful implementation of CW/DFT PFET depends on overcoming two separate challenges: the first is PFET performance with mixed basis sets and the second is PFET performance upon moving beyond DFT. It is common practice to use Gaussian type orbitals (GTO) for CW calculations and PWs in a DFT description of an extended environment. Straightforward utilization of mixed GTO/PW basis sets causes the OEP procedure, a crucial step in PFET, to become numerically unstable. In Part I,12 we addressed this issue via re-optimization of the contraction coefficients of the GTO basis sets using principal component analysis. Test calculations showed that the re-optimized GTO basis sets are compatible with the PW basis sets in the OEP procedure and thus generate reliable results at the GTODFT/PW-DFT PFET level of theory. With a robust basis set representation in hand, we investigate the second challenge in the present paper, namely the accuracy of the hybrid CW/DFT PFET. Our goals here are three-fold: 1) to confirm the flexibility of the PFET formalism by demonstrating PFET calculations at the CW/DFT level; 2) to assess the accuracy of CW/DFT PFET by comparing GTOCW/PW-DFT PFET results with full GTO-CW benchmark calculations; and 3) to identify sources of errors in PFET and DFET and establish their respective ranges of application, by analyzing their fundamental theoretical differences. The present paper is organized as follows: first we briefly review DFET and PFET frameworks, specifically with respect to how they differently evaluate total energies. We will also introduce a mixed DFET/PFET formalism to compare with DFET and PFET. Computational details are provided thereafter. We then investigate several specific applications, including water dimers, and adsorption of small molecules on a metal (Al) as well as on metal oxide (MgO and Al2O3) surfaces. These results will clarify how pure DFT calculations can be systematically corrected by an accurate hybrid description at the CW/DFT level, as discussed in our concluding remarks. 2 Theories 4


Page 5 of 39

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.1 Density Functional Embedding Theory Here we briefly introduce the DFET method, the details of which can be found elsewhere.11 As usual, for notational simplicity we consider two subsystems, but the method is generalizable to an arbitrary number of subsystems. In DFET, we first perform a PW-DFT calculation for the total system to obtain the total DFT density as a reference ( ρtot ). Then for an arbitrary embedding potential Vemb shared by both subsystems, we define the following extended Wu-Yang functional:11, 13

r r r W [Vemb ] = E% I + E% II − ∫ Vemb (r ) ⋅ ρtot (r )dr

(2)

Here, the tilde sign denotes inclusion of the density-potential integral term in the energy, i.e.,

r E% K = EK [ ρ K [Vemb ]] + ∫ ρ K [Vemb ] ⋅Vemb dr with K ( K = I , II ), and ρ K [Vemb ] the ground state density of subsystem K in the presence of the embedding potential Vemb . Maximizing W with respect to Vemb DFET , DFT / DFT yields a unique embedding potential (denoted as Vemb ) which reproduces the total density as the

sum of the two subsystem densities ( ρ I + ρ II = ρtot ). Because both subsystem densities and the total DFET , DFT / DFT system density ( ρ I , ρ II , and ρtot ) are computed using DFT while constructing Vemb , Vemb

effectively captures the mutual polarization between subsystem I and II at the pure DFT level. More specifically, the environment DFT density ρIIDFT is computed in the presence of ρIDFT , instead of ρ ICW , and is unchanged while we perform the subsequent embedded region calculations at the CW level. This is effectively a frozen environment approximation since we fix the environment density at the DFT/DFT level while we perform CW/DFT calculations. The frozen environment approximation is one of the two major approximations employed in DFET.

5



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 39

DFET , DFT / DFT To obtain the DFET energy, Vemb is then applied to the embedded region calculations as DFET , DFT / DFT ] and the embedded DFT an external potential to obtain the embedded CW energy E% ICW [Vemb DFET , DFT / DFT ] . The total energy is evaluated via energy E% IDFT [Vemb

(

DFET , DFT / DFT DFET , DFT / DFT  − E% IDFT Vemb  EtotDFET = EtotDFT + E% ICW Vemb

in an a posteriori way. Here,

( E%

CW I

DFT Etot

)

(3)

is the DFT energy of the entire system. The

)

DFET , DFT / DFT DFET , DFT / DFT Vemb  − E% IDFT Vemb  term in Eq. (3) represents the CW correction for the XC

energy within the embedded region. The DFET method is computationally more efficient compared to the DFET/PFET or PFET methods introduced below. The computational cost scaling of DFET with respect to system size is mainly limited by the DFT calculation of the total system, and a higher prefactor compared to standard DFT is introduced by the iterative process for the maximization of the extended Wu-Yang functional. To theoretically connect the working equation for the DFET energy (Eq. (3)) and the general DFT in an exact energy expression in ECW schemes [Eq. (1)], we decompose the total DFT energy Etot

way as: DFET , DFT / DFT DFET , DFT / DFT  + EIIDFT Vemb  + Eint  ρ IDFT , ρ IIDFT  EtotDFT = EIDFT Vemb

(4)

DFET , DFT / DFT DFET , DFT / DFT DFET , DFT / DFT DFET , DFT / DFT r  = EIDFT Vemb  + ∫ ρ IDFT [Vemb E% IDFT Vemb ] ⋅Vemb dr

(5)

DFET , DFT / DFT DFET , DFT / DFT DFET , DFT / DFT DFET , DFT / DFT r  = EICW Vemb  + ∫ ρ ICW [Vemb E% ICW Vemb ] ⋅ Vemb dr

(6)

We further have:

6


Page 7 of 39

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


Inserting Eqs. (4), (5), and (6) into Eq. (3) leads to:

(

DFET DFET , DFT / DFT DFET , DFT / DFT  + E IIDFT Vemb  + Eint  ρ IDFT , ρ IIDFT  + ∆ Etot = EICW Vemb

)

(7)

Where ∆ is: DFET , DFT / DFT r ∆ = ∫ ( ρ ICW − ρ IDFT ) ⋅Vemb dr

(8)

Eq. (7) indicates that the DFET energy is comprised of the CW embedded region energy ( EICW ), the DFT energy of the environment ( E IIDFT ) and the interaction between the two subsystems (

EintDFET = Eint  ρ IDFT , ρ IIDFT  + ∆ ). A frozen environment approximation is employed and thus the environment density here ( ρIIDFT ) is in principle inaccurate. Moreover, in Eint  ρ IDFT , ρ IIDFT  , the interaction energy is computed using the DFT embedded region density ρIDFT , instead of the CW embedded region density ρICW as required for a self-consistent CW/DFT method. As we will see later (in PFET), using ρICW should require not only ρICW but also a set of corresponding Kohn-Sham orbitals to obtain the exact kinetic energy of the interacting subsystems. The interaction energy error

(E

int

)

 ρ ICW , ρ IIDFT  − Eint  ρ IDFT , ρ IIDFT  can be large if ρICW differs significantly from ρIDFT . However, as

we show below, the ∆ term partly corrects for this error with a linear approximation to the interaction energy. In a converged DFET calculation, the total DFT density is exactly reproduced within the specific DFET functional used, hence Etot at the DFT/DFT level is exactly the ground state DFT energy of the total DFT system ( Etot ), which is a minimum with respect to variations in ρ I and ρ II . Therefore, we have:

7



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

δ EtotDFT δ EI δ Eint = + =0 δρ I δρ I δρ I

Page 8 of 39

(9)

Considering the variational principle for the embedded DFT calculation for subsystem I, we further have: r

DFET , DFT / DFT dr ) δ E δ E% I δ ( EI + ∫ ρ I ⋅Vemb DFET , DFT / DFT I = = + Vemb =0 δρ I δρ I δρ I

(10)

Combining Eq. (9) and (10), we find: DFET , DFT / DFT Vemb =

δ Eint δρ I

(11)

and inserting Eq. (11) into Eq. (8) yields: ∆ = ∫ ( ρ ICW − ρ IDFT ) ⋅

δ Eint ( ρ I , ρ II ) r dr δρ I

(12)

≈ Eint ( ρ ICW , ρ II ) − Eint ( ρ IDFT , ρ II )

Eq. (12) clearly shows the role of ∆ , which employs a linear approximation to account for the method-induced interaction energy error. However, the linearity of the interaction with respect to ρ I may break down, especially when the change in ρ I (i.e., ρ ICW − ρIDFT ) is large. This linear interaction energy approximation is the second major source of error in the DFET scheme.

2.2 Potential Functional Embedding Theory A more rigorous method to conduct the ECW calculation is within PFET, which was developed and tested only at the pure DFT level in previous work.10, 12 As in DFET, PFET assumes a common embedding potential Vemb for all subsystems and computes the densities of the embedded region and the 8


Page 9 of 39

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


environment with Vemb , but here PFET will self-consistently employ CW for the former region and DFT CW for the latter. Then, the exact total energy at the CW/DFT level can be derived using ρI [Vemb ] and

ρIIDFT [Vemb ] : CW / DFT Etot [Vemb ] = EICW [Vemb ] + EIDFT [Vemb ] + Eint  ρ ICW , ρ IIDFT  I

(13)

Note that all three energy terms on the right hand side of Eq. (13) do not include the density times embedding potential integral term; Eq. (13) is equivalent to Eq. (5) in Part I5 after a simple CW rearrangement. It is straightforward to perform embedded CW or DFT calculations to obtain EI [Vemb ] DFT and EII [Vemb ] . To evaluate the interaction energy Eint  ρ ICW , ρ IIDFT  we invert ρ tot = ρ ICW + ρ IIDFT and

ρ ICW using an (unfortunately numerically quite demanding) OEP approach10 to determine the associated exact kinetic energy potentials. These can then be used together with the electrostatic and exchangecorrelation contributions to reconstruct at the DFT level the interaction energy as the difference

EKS [ ρ tot ] − EKS  ρ ICW  − EKS  ρ IIDFT  . The PFET scheme thus aims to directly minimize EtotCW / DFT with PFET ,CW / DFT respect to Vemb . The optimized embedding potential is denoted as Vemb with the final minimized

energy defined as the PFET energy EtotPFET ,CW / DFT . In principle, PFET also allows for optimization of the number of electrons in each subsystem, subject to the constraint that the total number of electrons is fixed. However, because our ultimate aim is to perform CW/DFT PFET calculations, we do not permit fractional subsystem occupation numbers here, in anticipation of the CW requirement of integer numbers of electrons. Note that Nafziger and Wasserman have shown that complete (e.g., PW) basis sets do not require optimization of subsystem electron numbers,14 substantiating our choice. 9



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 39

In contrast to DFET, no linear interaction energy approximation is made in PFET and the PFET ,CW / DFT evaluation of Eint is therefore exact at the CW/DFT level. Furthermore, Vemb is optimized with PFET ,CW / DFT interactions between ρICW and ρIIDFT , instead of ρIDFT and ρIIDFT . Therefore, Vemb polarizes the

environment in the correct way, and the environment density ρIIDFT is obtained self-consistently in presence of the CW density in the embedded region. Such a self-consistent CW-in-DFT embedding was first performed by Dresselhaus et al. using freeze-and-thaw cycles for a density matrix renormalization group CW approach in DFT.15 Note that PFET does not require a freeze-and-thaw approach, as all subsystem densities are optimized in parallel. The PFET scheme likewise removes both approximations/error sources in DFET, and lays out an exact formalism for CW/DFT embedding calculations. However, the exact evaluation of Eint is time-consuming, making PFET a much more expensive method compared to DFET.

2.3 DFET/PFET As outlined above, there are two major approximations in DFET, both of which are removed in PFET. To isolate the errors incurred by these two approximations, we propose a hybrid method combining DFET and PFET to facilitate analysis of the results. In this DFET/PFET scheme, we first DFET , DFT / DFT obtain Vemb using DFET, then use this DFET potential in the PFET energy formula [Eq. (13)]

without

further

optimization.

The

DFET/PFET

energy

can

thus

be

denoted

as

CW / DFT DFET , DFT / DFT DFET , DFT / DFT Vemb  . On one hand, because Vemb EtotDFET / PFET = Etot is used to evaluate ρIIDFT , the

environment is only polarized by a DFT description of the embedded region and the frozen environment CW / DFT approximation is still in force. On the other hand, the exact Eint is utilized in the evaluation of Etot ,

thus eliminating the linear interaction energy approximation used in DFET. By comparing DFET/PFET 10


Page 11 of 39

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


results with DFET and PFET, we therefore can separate the error sources and assess the quality of the linear interaction energy and frozen environment approximations, respectively. The hierarchy of the three ECW schemes (DFET, DFET/PFET, PFET) allows us to quantitatively analyze inaccuracies introduced by the underlying approximations and to systematically evaluate the performance of each. Below we apply these three ECW schemes to representative test cases in order to gain insight into their respective applicability. These four cases consider hydrogen bonding in the water dimer, dative and metallic bonding via CO adsorption and diffusion on an Al12 cluster, and dative and ionic bonding via water adsorption on magnesia (MgO) and alumina (Al2O3) clusters. Except for water dimer, all other cases aim to simulate local interactions between small molecules and extended surfaces, for which a hybrid GTO-CW/PW-DFT description would be needed. Here, however, we use large cluster models instead of extended surfaces as our total system so that we can afford to perform full GTO-CW benchmark calculations on the entire system. Even though for these tests we conduct effectively cluster-in-cluster calculations, we still adopt a mixed GTO/PW basis set setup consistent with the cluster-in-extended surface scenarios expected to be used in actual applications. Computational details for these test cases are introduced in the next section.

3 Computational Details For each test case and at each geometry, a DFET calculation is performed to generate DFET , DFT / DFT DFET , DFT / DFT Vemb . We then use Vemb to compute both the DFET and the DFET/PFET energies. To

DFET , DFT / DFT accelerate convergence, Vemb is also used as the initial guess for the PFET calculation to obtain

PFET ,CW / DFT Vemb and EtotPFET ,CW / DFT . The basis set type, the level of theory, and codes used to perform the

calculations are listed in Table 1.

11



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 39

Table 1. Basis set type, level of theory, and codes used in the DFET and PFET calculations. a) ABINIT 6.0.3/7.0.516 and MOLCAS 7.817 codes were modified to perform different types of calculations in the embedding theories: the prefix “NonSc-” indicates the code was modified to enable DFET potential optimization; the prefix “emb-” indicates the code was modified to perform subsystem calculations with an external embedding potential; and the prefix “invKS-” indicates the code was changed to conduct OEP calculations.

PFET

calculation type

Basis set type

level of theory

Code/version

DFET

PW

DFT/DFT

NonSc-ABINIT 6.0.311

subsystem (embedded region)

GTO

HF, CW

emb-MOLCAS 7.811

subsystem (environment)

PW

DFT

emb-ABINIT 7.0.510

OEP

PW

DFT

invKS-ABINIT 7.0.510

GTO

HF, CW

MOLCAS 7.811

Total system benchmark

All DFT calculations here employ the Perdew-Burke-Ernzerhof (PBE) XC functional.18 All atoms except O utilize the Trail-Needs ECP/pseudopotential pairs19 while we use the OPIUM-derived20 soft ECP/pseudopotential for O, as discussed in part I. All PW-based calculations employ a PW kinetic energy cutoff of 1400 eV to converge the total energy to within 10 meV/atom. All atoms in GTO-based calculations are described by the Dunning-style optimized VTZ basis sets of part I,12 except for GTOCW benchmark calculations of H2O(MgO)4/(MgO)12 and H2O(Al4O6H6O3)/(Al4O6H6O3), in which Dunning-style optimized VDZ basis sets12 are used for the Mg, Al, and O atoms in the environment to reduce computational cost. In all DFET and PFET calculations, the extended Wu-Yang functional value is converged to within 1.0×10-5 hartree.

12


Page 13 of 39

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


For the hydrogen bonding test calculations, the geometry of the water dimer was first optimized at the DFT level within the MOLCAS code using the optimized GTO basis set12 and the B3LYP21 hybrid XC functional. We then place this optimized water dimer in the center of a 12 Å × 16 Å × 12 Å periodic cell. We consider each water molecule as a subsystem in the embedding calculations, with the whole system denoted as H2O/H2O. See also supplementary material for Part I. For the CO on Al12 test calculations, we place an Al12 cluster in the center of a 14.3 Å × 14.3 Å × 17 Å hexagonal periodic cell. We carve the Al12 cluster out of a face-centered cubic (fcc) Al slab to mimic an 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. The six central Al atoms and the CO molecule comprise one subsystem and the six surrounding Al atoms are the other subsystem, denoted as (CO)Al6/Al6. This is the same test case discussed in part I. For the H2O-(MgO)16 test calculations, we place a (MgO)16 cluster in the center of a 16 Å × 16 Å × 20 Å box. We carve the (MgO)16 cluster out of the MgO rock-salt structure with a lattice constant of 4.28 Å corresponding to the DFT-PBE equilibrium lattice constant with the HF-based pseudopotentials. The central (MgO)4 cube of the (MgO)16 cluster and the H2O molecule comprise one subsystem and the surrounding (MgO)12 is the other subsystem, denoted as H2O(MgO)4/(MgO)12. See also supplementary material for Part I for more tests on this system. For the H2O-(Al4O6H6O3)2 test calculations, we place an (Al4O6H6O3)2 cluster in the center of a 19.0 Å × 19.0 Å × 18.1 Å hexagonal supercell. We first carve an (Al2O3)4 cluster out of the Al2O3 corundum phase with experimental lattice constants of a = b = 4.75 Å and c = 12.99 Å.22 We terminate the unsaturated O atoms with H, and the unsaturated Al atoms with OH, and optimize within ABINIT the geometries of the terminating groups with the original (Al2O3)4 part fixed. This results in an 13



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 39

(Al4O6H6O3)2 cluster retaining the corundum phase geometry. We then adsorb a water molecule onto the rigid (Al4O6H6O3)2 cluster and optimize the geometries of the molecular- and dissociative-adsorption states at the DFT-PBE level within ABINIT. The resulting H2O-(Al4O6H6O3)2 cluster is then partitioned into an embedded region comprised of the adsorbate and upper part of the metal oxide cluster, with the lower part of the metal oxide cluster as the environment (vide infra). This partitioning is denoted as H2O(Al4O6H6O3)/(Al4O6H6O3). All of the specific geometries are shown in figures accompanying the results and discussion of the respective systems. Preliminary test calculations for these four cases are carried out at the Hartree-Fock (HF)/DFT level. However, we still use the terminology “ECW scheme” to refer to these calculations. With the insights gained from the HF/DFT embedding calculations, we then proceed to perform complete active space self-consistent-field (CASSCF)23 / DFT embedding calculations.

4 Results and Discussion 4.1 Water Dimer and CO-Al12 One important application of the ECW method is to study chemical phenomena or reactions in aqueous solution.24 Typically, the system is partitioned into two subsystems, an embedded region comprised of the solute molecules and nearby solvent water molecules and the environment region being the rest of the solvent molecules. The two subsystems influence each other via various interactions, including electrostatic, polarization, dispersion, and especially hydrogen bonding. We examine the water dimer as the simplest example of these interactions in order to test the quality of the ECW method in describing them. We investigate the stretching mode of the O-H bond involved in hydrogen bonding for a pair of interacting water molecules (inset of Figure 1), where the rigid water is “the environment” and

14


Page 15 of 39

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 vibrating water is the “embedded” one. The energy profiles of this motion predicted at five different levels of theory are displayed in Figure 1. All the curves are shifted so that they share the same energy value at lΟ1−Η1 = 0.98 Å for ease of comparison. The DFT results for the entire system are given as a baseline, and the pure HF results for the entire dimer are provided as the benchmark for the three HF/DFT embedding methods. The HF and DFT energy curves (red and gray, respectively) feature significantly different equilibrium bond lengths (0.95 Å versus 0.98 Å), demonstrating a qualitative difference between the two methods. This deviation is consistently removed by all three HF/DFT embedding methods (see Figure 1). The equilibrium bond length predicted by all three embedding schemes is 0.95 Å, the same as the full HF result. Furthermore, all three ECW methods generate similar results, indicating that neither DFET/PFET nor PFET show significant improvements on DFET in this case.

Figure 1. Energy changes induced by stretching the O1-H1 bond in an otherwise rigid water dimer (see inset). Green, blue, and black curves are respectively HF/DFT DFET, HF/DFT DFET/PFET, and HF/DFT PFET energies, while the red and gray curves are respectively the HF and DFT results for the entire system. The curves are shifted to have the same total energy at an O1-H1 distance of 0.98 Å. 15



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 39

The fair agreement between the DFET result (blue curve in Figure 1) and the HF benchmark

(

)

DFET , DFT / DFT DFET , DFT / DFT  − E% IDFT Vemb  in Eq. (3). indicates the effectiveness of the correction term E% IHF Vemb

In the presence of the embedding potential, the method-induced difference within the embedded region accounts for the major portion of the difference compared to the full system, proving the validity of the DFET scheme. The improvement upon the DFET results is negligible when an explicit interaction energy term ( Eint  ρ IHF , ρ IIDFT  ) is introduced in the DFET/PFET scheme (green curve in Figure 1). We take the calculations at lΟ1−Η1 = 0.98 Å as representative examples for all geometries along the potential energy curves for a more detailed analysis. Embedded DFT predicts a Mulliken charge of 0.57 e on O1 and 0.32 e on H1, while embedded HF predicts -0.69 e on O1 and 0.38 e on H1, indicating a noticeable difference between ρIDFT and ρ IHF . Explicit evaluation of the interaction energies in DFET (

Eint  ρ IDFT , ρ IIDFT  ) and DFET/PFET ( Eint  ρ IHF , ρ IIDFT  )10 reveals they differ only slightly (by 3.6 meV). Such a small difference proves that, although ρIDFT and ρ IHF differ nontrivially, the linear approximation employed in Eq. (12) of DFET is accurate in this case. The small deviation between DFET and DFET/PFET is fairly constant for different geometries, even if an error associated with the linear approximation exists, and cancels out when considering only relative energy profiles. The effectiveness of the frozen environment approximation of DFET is also proven by the lack of improvement when employing the fully self-consistent PFET, showing that ρIDFT and ρ IHF polarize the environment density in a fairly similar way for this particular system. Thus, the three ECW schemes perform equally well for the water dimer test case at the HF level of theory for the embedded region. The residual deviation from the HF benchmark stems from the DFTlevel description of the environment and the interaction between it and the embedded region, rather than 16


Page 17 of 39

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 lack of environment density relaxation. We also performed a PFET calculation at the HF/HF level, with the interaction energy still determined at the DFT level, which results in an energy curve very similar to the ones produced by the three embedding schemes at the HF/DFT level (not shown). This suggests that the residual mismatch in energy is due to the DFT-level description of the interaction energy. These findings provide insights into the method choice for applications, although the test calculation results for the water dimer do not yield any theoretical advantage of the PFET scheme. The insensitivity of the total energy with respect to different methods in this case is attributed to the closeshell nature of each subsystem. There are no actual chemical bonds formed between the subsystems for this test case. The subsystem interaction therefore is relatively weak, and the curvature of Eint and polarization differences are not observable. We conclude that DFET is a cheap yet accurate embedding method to study weak (hydrogen) bonding interactions. We proceed to test the different ECW schemes on systems with strong chemical bonds between the subsystems. Specifically, the DFET-based ECW scheme has been applied to study a variety of metallic surface phenomena.11,

25

It therefore is important to compare the performance of the higher-

level DFET/PFET and PFET schemes to the DFET scheme for metallic environment partitioning. We consider here the nearly-free-electron-like metal Al and study a finite Al12 cluster (Figure 2(a)), so that benchmark calculations on the full system are affordable.

17



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 39

Figure 2. (a) Geometry of (CO)Al6/Al6. (b) Binding curve of a CO molecule adsorbed on an Al12 cluster. Green, blue, and black curves are respectively the results obtained from HF/DFT DFET, HF/DFT DFET/PFET, and HF/DFT PFET, while the red and gray curves are respectively the benchmark HF and DFT results for the entire system. The energies have been shifted to match at d Al −C = 2.2 Å. Note that DFT is expected to provide a much better description of the metal than HF does, but a worse description of the metal-CO interaction due to the artificially low CO 2π* level predicted by DFT that subsequently over-binds the CO molecule (vide infra). It is therefore a bit unclear as to whether DFT or HF will give a better description of the entire system. As these HF calculations are just meant for testing, with the more physically realistic embedded CW calculations to follow, we still refer to these HF calculations as benchmarks despite their likely inaccuracy. Here, we adsorb a CO molecule on the fcc hollow site (Figure 2(a)) and vary the distance between the top layer of the Al12 cluster and the C atom ( d Al-C ) from 1.0 to 2.2 Å to investigate the resulting energy change, while keeping all other geometrical parameters fixed. The energy curves predicted by the three ECW schemes are displayed in Figure 2(b)), along with the HF and DFT results for the entire system. Note the significant binding of CO to the cluster predicted by DFT, which is probably an artifact, as mentioned above. All HF/DFT 18


Page 19 of 39

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


embedding schemes remove this over-binding and approach the HF benchmark, although none of them reach perfect agreement. Once again, no significant difference among the three HF/DFT embedding approaches is observed. We investigate this result in more detail next. The system with d Al-C = 1.4 Å is taken as a representative point on the curve to analyze the single point energy calculation. The nonlinear interaction energy correction term ( EintDFET / PFET − EintDFET ) is -0.093 eV when comparing DFET and DFET/PFET results, much more than an order of magnitude larger than in the previous hydrogen-bonding case (3.6 meV). The energy evaluation is thus more sensitive to the embedded region’s electron density distribution in a metallic environment compared to the hydrogenPFET , HF / DFT bonding case. The total energy Etot [Vemb ] is also lowered by 0.02 eV with a fully optimized

PFET , HF / DFT

embedding potential Vemb

DFET , DFT / DFT

evaluated with Vemb

in the PFET calculation, compared to the DFET/PFET energy

. The improvement by PFET is not large enough to substantially further

approach the HF benchmark compared to DFET or DFET/PFET, even though the differences between DFET, DFET/PFET, and PFET (on the order of 10 meV) are significantly larger than the water dimer case (on the order of 1 meV). The difference between the HF and the HF/DFT embedding methods are on the order of few tenths (~0.2) eV and cannot be amended by a better description of the density of the embedded region. The residual error therefore is more likely due to the inaccurate DFT description of the environment itself. This error will endure due to the DFT nature of the environment in a hybrid CW/DFT method and thus would require either a larger embedded region or use of PFET with more than one subsystem treated with CW theory. We next perform embedding calculations at the CASSCF/DFT level for the above two cases and consider the energy changes with the same geometry variations. Figure 3 displays the curves evaluated with the embedding schemes at the CASSCF/DFT level, compared to DFT and CASSCF results for the 19



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 39

entire system. We include the σ and σ* orbitals of the two O-H bonds in the embedded H2O, and the 2p lone electrons and orbitals of O1 to form a (6e, 5o) active space for water dimer. We include the two π and two π* orbitals of the CO molecule, along with the three highest occupied molecular orbitals and three lowest unoccupied molecular orbitals from the Al6 cluster to form a (10e, 10o) active space for (CO)Al6/Al6. The CASSCF calculations for the entire system adopt the same active space composition for a fair comparison. We observe that the results of the embedding schemes (the green, blue and black curves) improve upon the pure DFT results (gray curves) in the right direction towards the CASSCF energy curves (in red) in both cases. In the water dimer case, there is only a slight difference on the order of 10 meV between DFET results and the other two, while the DFET/PFET and PFET schemes perform equivalently. In the (CO)Al6/Al6 case, we observe no significant difference among the three embedding schemes. The comparisons are similar to what was observed at the HF/DFT level. The underlying reason also holds here. The bonding type at the interface of the two subsystems affects the relative performance of the embedding schemes, regardless of the specific non-DFT method used to treat the embedded region.

20


Page 21 of 39

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 3. (a) Energy changes induced by displacing the O1-H1 bond away from equilibrium in an otherwise rigid water dimer (inset of Figure 1). The curves are shifted to have the same total energy at an O1-H1 distance of 0.98 Å. (b) Binding energy curve for a CO molecule adsorbed on an Al12 cluster (Figure 2(a)). The energies have been shifted to match at d Al −C = 2.2 Å. Green, blue, and black curves are respectively the results of the CASSCF/DFT DFET, CASSCF/DFT DFET/PFET, and CASSCF/DFT PFET schemes, while the red and gray curves are respectively the CASSCF and DFT results for the entire system. In summary, we observe equivalent performances of the three schemes at the HF/DFT and CASSCF/DFT levels in both test cases analyzed above. This indicates that changes in polarization of the environment induced by density variations in the embedded region (due to introduction of a non-DFT treatment there) are too small in both hydrogen bonding and metallic bonding to require the selfconsistent PFET embedding scheme. This conclusion validates the extensive use thus far of DFET for chemisorption on metal surfaces.

The short screening length for metals and the weak bonding

associated with hydrogen bonds are likely the reason that DFET proves sufficient.

4.2 Adsorption of H2O on Metal Oxide Clusters 21



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 39

Adsorption and reaction of water on metal oxide surfaces occurs in catalysis and geochemistry, among other applications. To this end, we examine water adsorption on two simple ionic materials, magnesia and alumina. We first investigate water adsorption on MgO(100). We use a single water molecule and an (MgO)16 cluster (Figure 4), instead of an extended surface as our total system so as to be able to carry out wavefunction-based benchmark calculations on the entire system. Two orientations of the water molecule are studied, as depicted in Figure 4. In geometry A, the water molecule is perpendicular to the surface with its O atom (Ow) above the center of the embedded (MgO)4 cluster pointing towards vacuum, with its two H atoms pointing toward the O atoms in the top layer of the embedded (MgO)4 cluster (the Oc atoms). The distance between Ow and the MgO surface ( dMgO-Ow ) is varied from 2.0 to 5.0 Å, while keeping all other geometrical parameters fixed (the geometry of the rigid water molecule is the same as that used in the water dimer case). In geometry B, the water molecule is parallel to the MgO surface. The Ow atom is along the diagonal line connecting the two Mg atoms in the MgO lattice, 0.806 Å away from the embedded cluster center (see Figure 4(b2)), and the two H atoms are pointing symmetrically toward the two Oc atoms.

22


Page 23 of 39

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 4. (a1) Geometry A of H2O(MgO)4/(MgO)12; (a2) top and (a3) side views of geometry A. (b1) Geometry B of H2O(MgO)4/(MgO)12; (b2) top and (b3) side views of geometry B. Solid circles surround the O atom in the water molecule (Ow), while dashed circles surround O atoms in the top layer of the (MgO)4 cluster (Oc). DFET , DFT / DFT

We construct Vemb

via DFET calculations and subsequently perform DFET and DFET , DFT / DFT

DFET/PFET calculations. We then use Vemb

as the initial guess to perform self-consistent

HF/DFT PFET and CASSCF/DFT PFET calculations. We take geometry A of H2O(MgO)4/(MgO)12 with dMgO-Ow = 2.6 Å (the approximate equilibrium distance) as an example for single-point calculation analysis. We first examine the difference in HF and DFT embedded region densities ρI − ρI HF

DFT

(Figure

DFET , DFT / DFT

5 (a1) and (a2)), with both densities optimized in the presence of the embedding potential Vemb

.

The positive values around the embedded cluster oxygen atoms and the negative values around the embedded cluster magnesium atoms indicate that HF unsurprisingly predicts more ionic bonding for (MgO)4 compared to DFT. The self-interaction error in DFT leads to over-delocalization of the electrons, resulting in a less ionic electron distribution compared to HF. The nonlinear interaction energy 23



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

correction corresponding to this density difference

(E

DFET / PFET int

Page 24 of 39

DFET − Eint ) is 1.45 eV, much larger than in

the previous types of bonding cases considered.

Figure 5. (a1) Isosurface and (a2) contour plot - crossing the top layer of embedded region atoms - of DFET , DFT / DFT DFET , DFT / DFT  − ρ IDFT Vemb  . (b1) Isosurface and (b2) contour plot - crossing the top layer ρ IHF Vemb PFET , HF / DFT DFET , DFT / DFT  − ρ IIDFT Vemb  . (c1) Isosurface and (c2) contour of environment atoms – of ρ IIDFT Vemb

24


Page 25 of 39

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


PFET , HF / DFT DFET , DFT / DFT  − ρ IHF Vemb  . plot - crossing the top layer of embedded region atoms – of ρ IHF Vemb

Geometry A of Figure 4 (a1)-(a3) is used, with dMgO-Ow = 2.6 Å. Yellow (blue) represents positive (negative) density differences at the isosurfaces. Circles (crosses) represent O (Mg) atoms in the contour plots. All isosurfaces are at ±1.7 × 10-3 a.u. Different embedded region densities ( ρ IHF and ρIDFT ) will produce distinct environment densities ( ρIIDFT ) because of changes in polarization. This effect can be observed by comparing ρIIDFT obtained with the self-consistent PFET and the non-self-consistent DFET schemes. Figure 5 (b1) and (b2) display this difference between the environment densities evaluated within the two schemes using the two PFET , HF / DFT DFET , DFT / DFT embedding potentials ( Vemb and Vemb ).

DFET , DFT / DFT  corresponds to the ρ IIDFT Vemb

PFET , HF / DFT  environment density polarized by the embedded region’s DFT density, while ρIIDFT Vemb

corresponds to the environment density polarized by the embedded region’s HF density. The major difference occurs around the O atoms in the environment [denoted by circles in Figure 5(b2)]. The electrons are attracted to the interface region between the embedded region and the environment. This phenomenon can be understood from electrostatics: the Mg atoms in the embedded region are more positively charged in PFET by HF than in DFET by DFT (+1.5 e at the HF level versus +1.2 e at the DFT level). When we compute the polarization of the environment at the HF/DFT level via PFET, the Mg atoms in the embedded region attract more electron density that accumulates on the environmental O atoms in the vicinity, producing a more ionic environment. Figure 5 (c1) and (c2) display the difference between the embedded region densities evaluated at the HF level with the two embedding potentials ( ρI [Vemb HF

PFET , HF / DFT

, DFT / DFT ] and ρIHF [VemDFET ] ). This b

difference represents the higher-order back polarization of the embedded region caused by the different environment densities evaluated at DFET and PFET levels. We also observe a decrease of 1.57 eV in the 25



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 39

PFET , HF / DFT

PFET energy for the selected geometry when using the optimized Vemb DFET , DFT / DFT

DFET/PFET energy evaluated with Vemb

compared to the

. This large change indicates that the embedding

potential correctly accounts for the mutual polarization between the HF-embedded region and the DFTenvironment and does indeed lower the total energy. In this case, the frozen environment approximation also fails and the fully self-consistent PFET scheme is necessary to capture all the physics. Note that no benchmark exists to directly verify the correctness of the calculated density polarizations in the hybrid HF/DFT embedding scheme. We therefore can examine only relative energy curves and use energy differences as our gauge to compare different embedding methods. PFET waterMgO surface binding energy curves are plotted in Figure 6, along with DFET and DFET/PFET results for comparison. It is surprising that the HF/DFT DFET (green curves) results and full HF benchmarks (red curves) are in reasonable agreement (maximum errors of only 50 meV for geometry A and 3 meV for geometry B). Despite the linear interaction energy approximation and the frozen environment approximation employed in DFET, the a posteriori DFET energy formalism seems to do a decent job in reproducing the HF benchmark energy. However, the DFET/PFET scheme predicts deeper wells (blue curves in Figure 6) for both geometries, incurring large errors. The significant difference between DFET/PFET and DFET represents the breakdown of the linear interaction energy approximation, in line with our previous analysis. The DFET/PFET scheme, which removes only one of the two error sources, leads to worse agreement with the full HF benchmark. By removing the errors from both approximations, the fully-self-consistent PFET results (black curves) agree better with the HF benchmark than DFET, with a smaller residual error (30 meV for geometry A and 0.7 meV for geometry B). This validates the performance of the HF/DFT PFET scheme and proves that the more rigorous energy expression also requires use of the fully-self-consistent HF/DFT embedding potential constructed via the PFET scheme. Comparison of the three levels of embedded wavefunction schemes implies that 26


Page 27 of 39

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 fair performance of the DFET scheme may be due to fortuitous error cancellation here. The two sources of inaccuracy – the linear interaction energy approximation Δ [Eq. (8)] and the frozen environment approximation - cancel each other out except for a fairly constant residual error, which cannot be observed in the shifted energy curve.

Figure 6. Water-MgO surface binding energy curves using the H2O(MgO)4/(MgO)12 cluster model in (a) geometry A and (b) geometry B (Figure 4). Green, blue, and black curves respectively are the results of the HF/DFT DFET, HF/DFT DFET/PFET, and HF/DFT PFET schemes, while red and gray curves are respectively the HF benchmark and DFT results for the entire system. All energies are shifted to match at 5.0 Å. We next proceed to test PFET performance for these types of bonds at an actual CW level for the PFET , HF / DFT

embedded region, namely, at the CASSCF/DFT PFET level. We use Vemb

as the initial guess for

the GTO-CASSCF/PW-DFT PFET calculation. Hydrogen bonds or dative bonds are formed as the water molecule approaches the MgO cluster, depending on the orientation of the water molecule. All electrons and orbitals involved in bond formation comprise the active space for each orientation. For geometry A, 27



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 39

we include the two σ and two σ* orbitals of the O-H bonds in H2O, the 2p lone pair of the Ow atom (Figure 4 (a1)), and the 2pz lone pairs of the two Oc atoms (Figure 4 (a1)) to comprise a (10e, 7o) active space. For geometry B, we include the two σ and two σ* orbitals of the O-H bonds in H2O, the 2s and 2p lone pair orbitals of the Ow (Figure 4 (b1)) and two 3s orbitals of the surface Mg atoms to comprise a (8e, 8o) active space. In the case of geometry B, the shortest distance between the H atom and the Oc (Figure 4 (b1)) is 2.2 Å. At this distance, the hydrogen bond is not yet formed and therefore we do not include the 2p orbitals of the Oc. The energy curves predicted by CASSCF/DFT PFET are displayed in Figure 7, along with benchmark CASSCF results for the full system and DFT results on the full system for comparison. We verified that the active space for the embedded cluster calculation has exactly the same composition as the active space for the benchmark calculation on the full system, so both calculations are depicting the same electronic states and thus the results are comparable. The CASSCF/DFT PFET results agree quite well with the CASSCF benchmark for both geometries.

Figure 7. Water-MgO surface binding energy curves using the H2O(MgO)4/(MgO)12 cluster model in (a) geometry A and (b) geometry B (Figure 4). Black curve is the result of the CASSCF/DFT PFET 28


Page 29 of 39

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


scheme, red is the CASSCF benchmark, and gray is the DFT result for the entire system. All energies are shifted to match at 5.0 Å. PFET , HF / DFT

We noticed that the initial guess Vemb

in the CASSCF/DFT PFET embedding calculations PFET ,CASSCF / DFT

is already a very good approximation of the optimal embedding potential Vemb

. For the

geometry with dMgO-Ow = 2.6 Å, the single point total energy is lowered by only 4.9 meV after the PFET , HF / DFT

optimization, compared to the initial value. This small energy change indicates that Vemb PFET ,CASSCF / DFT

already fairly close to the final converged Vemb

. We also compare the two relative

PFET , HF / DFT

CASSCF/DFT PFET energy curves evaluated with Vemb

and

PFET ,CASSCF / DFT Vemb , respectively PFET , HF / DFT

(Figure 8). The good agreement between the two curves shows again that Vemb PFET ,CASSCF / DFT

approximation to Vemb

is

is a good

. The equivalent performance of the two embedding potentials is

attributed to the similar density distributions predicted by HF and CASSCF in the embedded region. The interaction energy difference term

(E

int

 ρ ICASSCF , ρ IIDFT  − Eint  ρ IHF , ρ IIDFT 

)

is on the meV scale,

indicating that a CASSCF embedded region polarizes the environment similarly to the HF embedded region. This observation indicates that a more efficient HF/DFT PFET calculation should be performed PFET , HF / DFT

to get Vemb

prior to a more expensive CW/DFT PFET calculation. The embedding potential

PFET , HF / DFT Vemb then can be used as a good starting point for further optimization at a higher level of theory.

Such a strategy will have a significant performance advantage, as it reduces the number of expensive CW calculations needed to achieve an optimized PFET embedding potential. It thus could partially alleviate the performance penalty associated by PFET when compared to DFET and make it a more feasible option. 29



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 39

Figure 8. Water-MgO surface binding energy curves using the H2O(MgO)4/(MgO)12 cluster model in geometries A and B (Figure 4). The CASSCF/DFT energies are evaluated with a PFET energy expression using optimized embedding potentials (either HF or CASSCF) at each different geometry. PFET ,CASSCF / DFT PFET , HF / DFT The solid red (black) curve is for Vemb ( Vemb ) in geometry A, while the corresponding

dashed curves are for geometry B. All energies are shifted to match at 5.0 Å. The H2O(MgO)4/(MgO)12 results demonstrate the advantages of PFET. To further validate its performance, we conduct embedding calculations with the three schemes (DFET, DFET/PFET, and PFET) for another ionic environment: water adsorption on the Al-terminated α-Al2O3(0001) surface. We specifically employ H2O-(Al4O6H6O3)2 as a cluster model to investigate the energetics of the molecularand dissociative-adsorption states of the H2O molecule (Figure 9).

30


Page 31 of 39

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 9. Water adsorption on (Al4O6H6O3)2. (a1) top and (a2) side views of the molecular-adsorption state; (b1) top and (b2) side views of the dissociative-adsorption state; (c1) top and (c2) side views of the desorption state. Solid circles surround O atoms in the water molecule (Ow), while dashed circles surround the two O atoms adjacent to the water molecule on the Al4O6H6O3 embedded cluster (Oc). We evaluate the adsorption energy as:

∆Eads = Eads − Edes

(14)

in which Eads is the energy of the adsorbed state and Edes is the energy of the desorbed state. A more negative adsorption energy thus indicates a more energetically favored adsorption state. We consider first the molecular- ( ∆Emol ) and dissociative- ( ∆Ediss ) adsorption energies produced by the three embedding schemes at the HF/DFT level, along with HF and DFT results for comparison ( Table 2). 31



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 39

Table 2. Molecular and dissociative adsorption energies predicted by different levels of theory. molecular-adsorption energy ∆Emol (eV)

dissociative-adsorption energy ∆Ediss (eV)

HF/DFT DFET

-2.00

-2.35

HF/DFT DFET/PFET

-2.16

-2.59

HF/DFT PFET

-2.13

-2.52

CASSCF/DFT DFET CASSCF/DFT DFET/PFET CASSCF/DFT PFET

-1.91

-2.27

-2.07

-2.60

-2.03

-2.51

HF

-2.18

-2.56

CASSCF

-2.08

-2.56

DFT

-1.81

-2.13

HF predicts more negative adsorption energies than DFT for both states. All three embedding schemes correctly predict stronger adsorption (more negative adsorption energies) at the HF/DFT level compared to DFT. HF/DFT DFET underestimates the adsorption energy compared to the HF benchmark by 0.18 eV for ∆Emol and by 0.21 eV for ∆Ediss . Excellent agreement is reached with HF/DFT DFET/PFET due to the exact interaction energy formula in the DFET/PFET scheme. Further optimizing the embedding potential at a hybrid HF/DFT level with the PFET scheme does not lead to a significant energetic difference (-2.13 eV vs. -2.16 eV for ∆Emol and -2.52 eV vs. -2.59 eV for ∆Ediss ). We therefore conclude that in this case the DFET/PFET scheme strikes the best balance of accuracy and computational cost at the HF/DFT level. The lack of improvement when using the full PFET signifies that the response of the embedding potential to changes in the environment density causes a consistent shift in energy for the three geometries, and thus can be neglected here. 32


Page 33 of 39

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


We proceed to carry out PFET calculations at the CASSCF/DFT level to assess if this conclusion PFET , HF / DFT

holds at a CW level of theory, with Vemb

as the initial guess for the embedding potential.

Dissociative adsorption breaks the Ow-H bond of the water molecule and forms an Al-Ow bond and Oc-H bond (Figure 9). We therefore include the two filled σ and two empty σ* orbitals of the Ow-H bonds, the filled 2s and 2p lone pair orbitals of the Ow atom, and the filled 2pz orbitals of the two adjacent Oc atoms (Figure 9), as well as one empty 3s orbital of the surface Al atom to comprise a (12e, 9o) active space. We keep the composition of the active space consistent for the three studied geometries. We again verified that the embedded cluster calculations and the benchmark calculations on the entire system correspond to the same active space for each geometry. The CASSCF/DFT PFET scheme predicts a molecular adsorption energy of -2.03 eV and a dissociative adsorption energy of -2.51 eV, both of which are in excellent agreement with the CASSCF benchmark (-2.08 eV for ∆Emol and -2.56 eV for ∆Ediss ). We also list the CASSCF/DFT energies evaluated with the DFET and DFET/PFET schemes respectively for comparison in Table 2. The results again confirm that both the DFET/PFET and PFET schemes perform very well at the hybrid CASSCF/DFT level, while the DFET scheme is less accurate when compared to the CASSCF benchmark. The good correspondence between PFET in both the HF and CASSCF cases also indicates that benefits from error cancellation between overestimated delocalization in the DFT subsystem and overbinding in the HF part is small. In summary, we observe for both metal oxides examined that all three embedding schemes improve upon pure DFT, with varying levels of agreement with the CW benchmark. In the H2O(MgO)16 case, the relative performance of the three schemes can be summarized as DFET ≈ PFET > DFET/PFET in terms of energy evaluation. In the H2O-(Al4O6H6O3)2 case, the relation becomes 33



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 39

DFET/PFET ≈ PFET > DFET. The inconsistency of the relative performances observed in the two cases indicates that the better performance of either the DFET scheme in the former case or the DFET/PFET scheme in the latter may be attributed to different levels of fortuitous error cancellation. In contrast, the excellent agreement of PFET with the benchmark in both cases demonstrates its robustness and accuracy, both of which stem from PFET’s theoretical rigorousness. We currently limit our test calculations to HF and CASSCF. Since PFET is based on variational minimization of the energy with respect to the embedding potential, we refrain from using nonvariational CW methods to ensure the validity of the PFET approach, thereby eschewing non-variational methods such as Møller-Plesset perturbation theory or coupled-cluster methods. One would first need to investigate the validity of PFET with non-variational methods, obtain a valid evaluation of ρ for perturbative approaches, and examine the relation

δE = ρ , in order to implement non-variational CW δ Vemb

methods with the PFET embedding scheme. Our focus here instead is on basis set compatibility and presenting the first benchmark calculations of PFET using CW techniques. Although the variational multi-reference configuration interaction (MRCI) method in principle could be used straightforwardly, its expense is prohibitive within the self-consistent PFET framework at present. Research in this direction is ongoing.

5 Conclusions We analyzed three ECW schemes, DFET, DFET/PFET, and PFET, to understand their sources of error and how each systematically offers improvement. The DFET scheme employs a DFT/DFT embedding potential to polarize the subsystems, using a linear approximation to evaluate the interaction between the CW embedded region and the DFT environment. The DFET/PFET scheme instead exactly evaluates this interaction energy at the DFT level but does not give a self-consistent description of the 34


Page 35 of 39

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


mutual polarization between the embedded region and the environment. The PFET scheme ameliorates both approximations in DFET, thus representing the most rigorous of the ECW methods. We conducted test calculations on four representative bonding cases [hydrogen bonding in water dimer; dative and metallic bonding in CO-Al12; and dative, hydrogen, and ionic bonding in H2O(MgO)16 and H2O-(Al4O6H6O3)2] to explore the performance of all three ECW schemes. We found that ECW performance depends on whether the type of interaction between subsystems is weak, screened, or strong. For hydrogen-bonded complexes and metals, the nonlinear effect of the interaction energy is weak and the polarization difference between non-DFT and DFT descriptions of the embedded region is small. Thus DFET already is a good approximation, and all three methods perform equally well. For ionic materials, we find case-dependent performances of the DFET and DFET/PFET schemes. In one case, the two errors associated with DFET compensate for each other, while the intermediate DFET/PFET scheme fails because it only fixes one of the two errors (the linear approximation of the interaction energy is fixed but not the not-fully optimized embedding potential at the CW/DFT level). In another case, the intermediate DFET/PFET scheme in contrast showed superior performance to DFET. This indicates that the fair agreement reached by either DFET or DFET/PFET is subject to error cancellation. The most rigorous PFET worked in all cases and led to excellent agreement with benchmarks. In summary, we have shown that CW/DFT PFET is the most rigorous and robust of the ECW schemes, while the accuracy of DFET and DFET/PFET is case-dependent. The intermediate DFET/PFET scheme also serves as a good diagnostic tool to truly understand the physics of the approximations underlying the DFET method.

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. 35



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 39

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 Dunning-style optimized VDZ basis sets used for the benchmark calculations. This information is available free of charge via the Internet at http://pubs.acs.org/.

References 1. Parr, R. G.; Yang, W., Density-Functional Theory of Atoms and Molecules. Oxford University Press: New York, 1989. 2. (a) Kresse, G.; Furthmüller, J., Efficient Iterative Schemes for ab initio Total-Energy Calculations using a Plane-Wave Basis Set. Phys Rev B 1996, 54, 11169-11186; (b) Kresse, G.; Furthmüller, J., Efficiency of ab-initio Total Energy Calculations for Metals and Semiconductors using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15-50; (c) 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; (d) 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; (e) 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: First-principles approach to material and nanosystem properties. Comput Phys Commun 2009, 180, 2582-2615; (f) 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.

36


Page 37 of 39

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


3. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J Chem Phys 2010, 132 (15), 154104. 4. Cohen, A. J.; Mori-Sánchez, P.; Yang, W., Insights into Current Limitations of Density Functional Theory. Science 2008, 321 (5890), 792-794. 5. (a) Szabo, A.; Ostlund, N. S., Modern Quantum Chemistry. Dover Publications Inc.: Mineola, 1996; (b) Helgaker, T.; Jorgensen, P.; Olsen, J., Molecular Electronic-Structure Theory. Wiley 2000. 6. (a) Marsman, M.; Grüneis, A.; Paier, J.; Kresse, G., Second-order Møller-Plesset perturbation theory applied to extended systems. I. Within the projector-augmented-wave formalism using a plane wave basis set. J Chem Phys 2009, 130, 184103; (b) Marsman, M.; Grüneis, A.; Paier, J.; Kresse, G., SecondOrder Møller–Plesset Perturbation Theory Applied to Extended Systems. I. Within the ProjectorAugmented-Wave Formalism Using a Plane Wave Basis Set. J Chem Phys 2009, 130 (18), 184103; (c) Booth, G. H.; Gruneis, A.; Kresse, G.; Alavi, A., Towards an Exact Description of Electronic Wavefunctions in Real Solids. Nature 2013, 493 (7432), 365-370. 7. (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 self-consistent 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; (c) 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; (d) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesolowski, T. A., The merits of the frozen-density embedding scheme to model solvatochromic shifts. J Chem Phys 2005, 122, 094115; (e) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F., Density functional theory embedding for correlated wavefunctions: improved methods for open-shell systems and transition metal complexes. J Chem Phys 2012, 137, 224113. 8. Jacob, C. R.; Neugebauer, J., Subsystem density-functional theory. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (4), 325-362. 9. Libisch, F.; Huang, C.; Carter, E. A., Embedded Correlated Wavefunction Schemes: Theory and Applications. Acc. Chem. Res. 2014, 47 (9), 2768-2775. 10. Huang, C.; Carter, E. A., Potential-functional embedding theory for molecules and materials. J. Chem. Phys. 2011, 135 (19), 194104. 11. Huang, C.; Pavone, M.; Carter, E. A., Quantum mechanical embedding theory based on a unique embedding potential. J. Chem. Phys. 2011, 134 (15), 154110. 12. Cheng, J.; Libisch, F.; Yu, K.; Chen, M.; Dieterich, J. M.; Carter, E. A., Potential Functional Embedding Theory at the Correlated Wavefunction Level, Part I: Mixed Basis Set Embedding. submitted 2016. 13. Yang, W. T.; Wu, Q., Direct method for optimized effective potentials in density-functional theory. Phys. Rev. Lett. 2002, 89 (14), 143002. 14. Nafziger, J.; Wasserman, A., Density-Based Partitioning Methods for Ground-State Molecular Calculations. J. Phys. Chem. A 2014, 118 (36), 7623-7639. 15. Dresselhaus, T.; Neugebauer, J.; Knecht, S.; Keller, S.; Ma, Y.; Reiher, M., Self-consistent embedding of density-matrix renormalization group wavefunctions in a density functional environment. J Chem Phys 2015, 142 (4), 044111. 16. (a) Gonze, X.; Amadon, B.; Anglade, P. M.; Beuken, J. M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Cote, 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.; 37



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 38 of 39

Verstraete, M. J.; Zerah, G.; Zwanziger, J. W., ABINIT: First-principles approach to material and nanosystem properties. Comput. Phys. Commun. 2009, 180 (12), 2582-2615; (b) Gonze, X.; Rignanese, G. M.; Verstraete, M.; Beuken, J. M.; 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.; Hamann, D. R.; Allan, D. C., A brief introduction to the ABINIT software package. Z. Kristallogr. 2005, 220 (5-6), 558-562. 17. Aquilante, F.; De Vico, L.; Ferre, N.; Ghigo, G.; Malmqvist, P. A.; Neogrady, P.; Pedersen, T. B.; Pitonak, M.; Reiher, M.; Roos, B. O.; Serrano-Andres, L.; Urban, M.; Veryazov, V.; Lindh, R., Software News and Update MOLCAS 7: The Next Generation. J. Comput. Chem. 2010, 31 (1), 224-247. 18. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77 (18), 3865-3868. 19. (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. 20. http://opium.sourceforge.net/sci.html. (accessed January 4, 2017). 21. (a) Becke, A. D., Density‐functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98 (7), 5648-5652; (b) 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 (2), 785-789; (c) 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 (8), 1200-1211; (d) 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 (45), 11623-11627. 22. Ishizawa, N.; Miyata, T.; Minato, I.; Marumo, F.; Iwai, S., A structural investigation of α-Al2O3 at 2170 K. Act. Cryst. B 1980, 36 (2), 228-230. 23. (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. 24. (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 (9), 094115; (b) Gomes, A. S. P.; Jacob, C. R.; Visscher, L., Calculation of local excitations in large systems by embedding wave-function theory in density-functional theory. Phys. Chem. Chem. Phys. 2008, 10 (35), 5353-5362; (c) Daday, C.; König, C.; Valsson, O.; Neugebauer, J.; Filippi, C., State-Specific Embedding Potentials for Excitation-Energy Calculations. J. Chem. Theory Comput. 2013, 9 (5), 2355-2367. 25. (a) 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; (b) 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: PlasmonInduced Dissociation of H2 on Au. Nano Lett. 2013, 13 (1), 240-247; (c) Libisch, F.; Cheng, J.; Carter, E. A., Electron-Transfer-Induced Dissociation of H-2 on Gold Nanoparticles: Excited-State Potential Energy Surfaces via Embedded Correlated Wavefunction Theory. Z. Phys. Chem. 2013, 227 (9-11), 1455-1466; (d) Libisch, F.; Krauter, C. M.; Carter, E. A., Corrigendum to: Plasmon-Driven Dissociation of H2 on Gold Nanoclusters. Z. Phys. Chem. 2016, 230 (1), 131-132; (e) Cheng, J.; Libisch, F.; Carter, E. A., Dissociative Adsorption of O2 on Al(111): The Role of Orientational Degrees of Freedom. J. Phys. Chem. Lett. 2015, 6 (9), 1661-1665; (f) Zhou, L.; Zhang, C.; McClain, M. J.; Manjavacas, A.; 38


Page 39 of 39

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


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), 14781484.

TOC graphic

39


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

Recommend Documents