Projection-Based Correlated Wave Function in Density Functional

Subscriber access provided by Kaohsiung Medical University

Article

Projection-Based Correlated Wave Function in Density Functional Theory Embedding for Periodic Systems Dhabih Vishaal Chulhai, and Jason D. Goodpaster J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01154 • Publication Date (Web): 01 Mar 2018 Downloaded from http://pubs.acs.org on March 3, 2018

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 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 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 52 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

Projection-Based Correlated Wave Function in Density Functional Theory Embedding for Periodic Systems Dhabih V. Chulhai and Jason D. Goodpaster∗ Department of Chemistry, University of Minnesota. 207 Pleasant St. SE, Minneapolis MN, 55455, USA. E-mail: [email protected]

Abstract We present a level shift projection operator-based embedding method for systems with periodic boundary conditions—where the “active” subsystem can be described using either density functional theory (DFT) or correlated wave function (WF) methods, and the “environment” is described using DFT. Our method allows for k-point sampling, is shown to be exactly equal to the canonical DFT solution of the full system under the limit that we use the full system basis to describe each subsystem, and can treat the active subsystem either with periodic boundary conditions—in what we term “periodic-in-periodic” embedding—or as a molecular cluster—in “cluster-in-periodic” embedding. We explore each of these methods, and show that cluster WF-in-periodic DFT embedding can accurately calculate the absorption energy of CO on to a Si(100)2×1 surface.

1

ACS Paragon Plus Environment

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

Introduction The choice of electronic structure method in computational chemistry is often a balancing act between accuracy and computational cost. For complex chemical systems, embedding methods promise the accuracy of a higher level (and more costly) method, at near the computational cost of a lower level (and cheaper) method. This is achieved by describing only a small region or subsystem of interest at the higher level method, while the remainder— a majority of the total system—is described at the lower level method. The challenge with embedding methods lies in accurately describing the interaction between subsystems, and Kohn-Sham density functional theory (KS-DFT) 1,2 —in the form of subsystem DFT 3–5 — provides a formally exact framework for embedding one quantum mechanical system in the environment of another. Within this framework, embedding may be performed using either approximations to the nonadditive kinetic potential, 6–8 optimized effective potentials, 9–15 embedded mean field theory, 16–20 or level shift projection operators (termed projection-based embedding). 21–36 Each of these has their own unique benefits and drawbacks. We are interested in using highly accurate embedded correlated wave function (WF) methods to describe chemistry in extended systems—that is, systems described by periodic boundary conditions. There is a long history of using embedded cluster models, where the region of interest is treated as a molecular cluster embedded in a periodic or extended environment, to accurately describe chemistry in periodic systems. 37–40 In such embedding methods, the environment may be described classically using electrostatic potentials 41–44 or mechanical and polarizable force fields, 45–53 or quantum mechanically through Green’s function methods 54–59 or methods based on the aforementioned subsystem DFT framework. In particular, Carter and coworkers have a celebrated history of embedding WF clusters in periodic DFT environments. 13,14,60–74 Their methods rely on using either approximate nonadditive kinetic potentials 60–63,65,66,75 or optimized effective potentials, 13,14,69,71,73,74 the latter of which requires an additional computational step to obtain the embedding potentials. More recently, accurate embedding in DFT environment subsystems 25–28,31–35 have been 2


Page 2 of 52

Page 3 of 52 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


achieved through the use of level shift projection operators. 42,76–81 These operators work by shifting the occupied orbital energies of the environment subsystems—from the perspective of the active subsystem—to higher energies, ensuring that the occupied orbitals of all subsystems are mutually orthogonal; this removes the need for approximate nonadditive kinetic potentials or any additional computational steps to obtain embedding potentials. Projectionbased embedding has been extensively used for molecular systems, 25–28,31–34,82–84 and we are also starting to see applications of projection-based embedding in periodic systems. 35 In that work, Libisch and coworkers first perform a total system KS-DFT calculation and project a subset of the total orbitals on to local basis functions to obtain a localized subsystem, on which they can subsequently perform a DFT calculation with a more expensive exchangecorrelation (XC) functional. However, their method does not yet include k-point sampling nor a WF treatment of the active subsystem. In this paper, we present projection-based DFT-in-DFT and WF-in-DFT embedding for periodic systems; our moniker “method A-in-method B” represents the methods used to describe the active subsystem (method A) and the environment subsystem (method B). Our DFT-in-DFT method describes the active DFT subsystem using either a periodic boundary conditions model or a non-periodic molecular (cluster) model, while our WF-in-DFT method describes the active WF subsystem using a cluster model. Our formulations allow for k-point sampling, and our method has been implemented into a Gaussian orbitals basis periodic code that then allows for easy subsequent description of the cluster WF region using standard molecular quantum chemistry software. This paper is organized as follows: In the theory section, we present projection-based embedding for periodic systems, including how one may approximate the active subsystem as a molecular cluster and still obtain a periodic embedding potential. In the results section, we first show that our formulation of projectionbased DFT-in-DFT periodic embedding is exact when both subsystems are described using the same basis and XC functionals, and how one may obtain accurate results using a more expensive XC functional to describe the subsystem of interest. We also explore how one may

3



Page 4 of 52

approximate the active subsystem as a molecular cluster, using cluster B3LYP-in-periodic LDA or PBE as a benchmark test case. Finally, we present our cluster WF-in-periodic DFT embedding method by calculating the absorption energy of CO on to a Si(100)-2×1 surface.

Theory Projection-Based Embedding in Molecular Systems We first present a summary of projection-based DFT-in-DFT and WF-in-DFT embedding for molecular systems. 25–28,31–34 For simplicity, we will divide the total system into two subsystems labeled A and B, and the total system density γ tot is obtained from the sum of subsystem densities as

γ tot = γ A + γ B

(1)

where γ A and γ B are the densities of subsystems A and B, respectively. The Fock matrix of subsystem A embedded in the “environment” of subsystem B can be written as

FA-in-B [γ A , γ B ] = hA-in-B [γ A , γ B ] + g[γ A ]

(2)

where g contains the two electron terms—that is, Coulomb and exchange for Hartree-Fock or Coulomb and exchange-correlation (XC) for DFT—and the embedded core Hamiltonian is hA-in-B [γ A , γ B ] = (3) A

B

A

B

B

h + g[γ + γ ] − g[γ ] + P [γ ] where h is the total one electron Hamiltonian that contains the kinetic and nuclear potential operators (from both subsystems), and PB is the level shift projection operator that enforces intersubsystem orbital orthogonality. We can analogously construct an embedded 4


Page 5 of 52 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


core Hamiltonian for subsystem B embedded in A, hB-in-A , with a corresponding level shift projection operator PA from subsystem A. Our projection-based embedding solve these embedded Fock matrices FA-in-B and FB-in-A self-consistently, in what is termed DFT-in-DFT embedding. For WF-in-DFT embedding, or DFT-in-DFT embedding using a higher accuracy XC functional for the active subsystem, we simply modify the core Hamiltonian of the active subsystem A by the embedded core Hamiltonian in equation 3. An example of a level shift projection operator PB is the Huzinaga operator 26,33,78,85

PB = −

i 1 h AB B BA F γ S + SAB γ B FBA 2

(4)

where FAB and SAB are elements of the total Fock matrix and overlap matrix described over the basis functions of subsystems A and B; see Refs. 26 and 33 for the implementation of this operator in projection-based embedding for molecular systems. For projection-based embedding in periodic systems, we will introduce a new level shift projection operator that is based on the Huzinaga operator.

Projection-Based Embedding in Periodic Systems For periodic systems, we take advantage of the fact that the Kohn-Sham crystalline orbitals— the periodic analogue of Kohn-Sham molecular orbitals—at each k-point are mutually orthogonal to the crystalline orbitals at all other k-points. Or,

hφ˜i (k)|φ˜j (k0 )i = δij δkk0

(5)

where φ˜i and φ˜j are Kohn-Sham crystalline orbitals i and j at k-points k and k0 , respectively, and δ is the Kronecker delta. In general, we will use the tilde (˜) to represent properties in reciprocal space. In this work, the crystalline orbitals are expanded using a set of Bloch functions {χ(k)} ˜

5



|φ˜i (k)i =

Page 6 of 52

C˜µi (k)|χ˜µ (k)i

X

(6)

µ

where C˜ are the expansion coefficients. These Bloch functions are, in turn, obtained from a set of real space, atom centered Gaussian-type orbitals {χ}

χ˜µ (k; r) =

X

eikT χµ (r − T)

(7)

T

where T is a lattice translation vector. As a result of the inherent orthogonality of crystalline orbitals at different k-points, we only need to ensure that the occupied crystalline orbitals of subsystem A at a particular kpoint (φÃ (k)) are orthogonal to the occupied crystalline orbitals of subsystem B at the same k-point (φ˜B (k)). The required projection operator is therefore independent of the properties of the subsystems at all other k-points. We introduce a new level shift projection operator that is based on the Huzinaga operator, which we termed the “Fermi-shifted Huzinaga operator”. This operator, in the Bloch functions basis matrix representation, is ˜ B (k) = P i 1 h ˜ AB ˜ AB (k) γ˜ B (k)S ˜ BA (k) F (k) − F S − 2

(8)

i

h

˜ AB (k)˜ ˜ BA (k) − F S ˜ BA (k) + S γ B (k) F

where γ˜ B (k) is the density matrix of subsystem B at k-point k, and F is the maximum ˜ AB/BA are elements of the total Fock matrix described Fermi energy among all subsystems. F using basis of of subsystems A and B; for example

AB ˆ tot |χ˜B (k)i F˜µν (k) = hχÃ µ (k)|F ν

(9)

where χÃ and χ˜B are the Bloch functions used to describe subsystems A and B, respec-

6


Page 7 of 52 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


˜ AB/BA is the overlap matrix between tively, and Fˆ tot is the total Fock operator. Similarly, S subsystems A and B, where each element is defined as

AB S˜µν (k) = hχÃ ˜B µ (k)|χ ν (k)i

(10)

This new Fermi-shifted projection operator is required for systems with positive occupied crystalline orbital energies. A more detailed discussion of this operator is provided in the Results section. The Fock matrix of subsystem A embedded in subsystem B becomes ˜ A-in-B [˜ F γ A , γ˜ B ; k] ˜ A-in-B [˜ ˜ [˜ =h γ A , γ˜ B ; k] + g γ A ; k] (11) ˜ ˜ B (k) ˜ [˜ = h(k) +g γ A + γ˜ B ; k] + P ˜ AA (k) + P ˜ B (k) =F ˜ A-in-B is defined within the equation (and is analogous to Equation 3), h ˜ is the core where h ˜ is the two electron potential matrix of the system, Hamiltonian of the periodic system, g ˜ AA are elements of the total Fock matrix spanning the Bloch functions of subsystem and F A. Note that we do not perform a separate calculation for the two electron potentials of subsystem A as we already have, and only need, the elements of the total Fock matrix—of ˜ AA is a sub-block and includes the potential due to both subsystems A and B—from which F calculating the Huzinaga operator. This means that at every step in the self-consistent cycle, the total Fock matrix at the lower level of theory is generated, and our method scales at best as the periodic KS-DFT on the full system at this lower level of theory. Using this total Fock matrix, we can also construct the Fock matrix of subsystem B embedded in A, ˜ B-in-A (k), and perform a self-consistent DFT-in-DFT embedding similar to that described F for molecular systems. Using this projection-based embedding framework, we will describe two types of embed-

7



Total system:

Periodic-in-periodic embedding: Subsystem A

Subsystem B

Cluster-in-periodic embedding: Subsystem A Subsystem B

Figure 1: Pictorial representation of periodic-in-periodic and cluster-in-periodic embedding of a one dimensional system. Both subsystems are described by periodic boundary conditions in periodic-in-periodic embedding, while the active subsystem is described as a non-periodic molecular cluster in cluster-in-periodic embedding.

8


Page 8 of 52

Start

Guess 𝜸 ෥ 𝐴 (𝒌) & 𝜸 ෥𝐵 𝒌 𝜸 ෥𝑡𝑜𝑡 𝒌 = 𝜸 ෥𝐴 𝒌 + 𝜸 ෥𝐵 𝒌 ෩ 𝑡𝑜𝑡 𝒌 Construct 𝑭 ෩ 𝐴-𝑖𝑛-𝐵 (𝒌) from 𝑭 ෩ 𝑡𝑜𝑡 𝒌 Get 𝑭 ෡ 𝐵 (𝒌) using 𝜸 ෥𝐵 𝒌 to construct 𝑷

SCF 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


Diagonalize and update 𝜸 ෥ 𝐴 (𝒌) ෩ 𝐵-𝑖𝑛-𝐴 (𝒌) from 𝑭 ෩ 𝑡𝑜𝑡 𝒌 Get 𝑭 𝐴 ෡ 𝐴 (𝒌) using 𝜸 ෥ 𝒌 to construct 𝑷

Diagonalize and update 𝜸 ෥𝐵 (𝒌) 𝜸 ෥ 𝐴 (𝒌) & 𝜸 ෥𝐵 𝒌 converged?

No

Yes 𝐴-𝑖𝑛-𝐵 ෩ Get 𝑭ℎ𝑖𝑔ℎ-𝑖𝑛-𝑙𝑜𝑤 (𝒌)

SCF II

Page 9 of 52

𝐴 Diagonalize and update 𝜸 ෥ℎ𝑖𝑔ℎ (𝒌)

𝐴 𝜸 ෥ℎ𝑖𝑔ℎ (𝒌) converged?

No

Yes End

Scheme 1: The periodic-in-periodic embedding workflow. We first converge subsystems A and B at the lower level of theory, then use these converged densities to construct an embedding potential for subsystem A at the higher level of theory.

9



Page 10 of 52

ding for periodic systems: 1) A periodic-in-periodic embedding where both subsystems are described using periodic boundary conditions; and 2) a cluster-in-periodic embedding where the active subsystem is described as a non-periodic molecular cluster. Figure 1 is a visual representation of the differences between these two embedding methods. Our periodic-inperiodic embedding algorithm is analogous to the method presented in ref. 33, and is outlined in Scheme 1. We separate our periodic-in-periodic embedding into a low level-in-low level DFT embedding (“SCF I” in Scheme 1), followed by a high level-in-low level DFT embedding (“SCF II”). For each step in the low level-in-low level embedding: 1) we generate the total Fock matrix at the lower level of theory using the sum of subsystem densities; 2) we divide this total Fock matrix into the embedded Fock matrices for each subsystem, applying the appropriate level shift projection operator; and 3) we then diagonalize each of these subsystem Fock matrices to obtain new subsystem density matrices. These steps are repeated until we achieve self-consistency for both subsystem densities—which, when using the full system basis to describe each subsystem, is identical to the canonical KS-DFT solution on the full system. From the converged low level-in-low level DFT embedding, one may perform a high levelin-low level embedding, where one subsystem—the subsystem of interest—may be described using more accurate but costly XC functional (“SCF II” in Scheme 1). The Fock matrix for such a high level-in-low level periodic DFT embedding is constructed as A ˜ A-in-B F γhigh , γ˜ A , γ˜ B ; k] high-in-low [˜

˜ ˜ low [˜ = h(k) +g γ A + γ˜ B ; k]

(12)

A ˜ B (k) ˜ low [˜ ˜ high [˜ −g γ A ; k] + g γhigh ; k] + P

˜ high and g ˜ low are the two electron potential matrices calculated at the high and low where g A levels of theory, respectively, and γ˜high is the density matrix of subsystem A at the high

level of theory. Note that for the high level-in-low level embedding, we do not update the density of subsystem B, which remains fixed at the self-consistent low level-in-low level DFT

10


Page 11 of 52 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


solution. This way, the high level calculation only has to be performed once, yet we obtain highly accurate results as was shown in ref. 33. The density of subsystem A at the high level A is solved self-consistently by using a modified core Hamiltonian as of theory γ˜high

A ˜ A-in-B [˜ h ˜ B ; k] high-in-low γ , γ

˜ ˜ low [˜ = h(k) +g γ A + γ˜ B ; k]

(13)

˜ B (k) ˜ low [˜ −g γ A ; k] + P This leads to an (uncorrected) high-in-low DFT-in-DFT energy, defined as uncorrected A Ehigh-in-low [˜ γhigh , γ˜ A , γ˜ B ]

h

i

˜ · (˜ = Tr h γ A + γ˜ B ) + Glow [˜ γ A + γ˜ B ] + ENN h

(14)

i

˜ A-in-B − Tr h ˜ A − Glow [˜ γ A] high-in-low · γ h

i

A A ˜ A-in-B + Tr h ˜high + Ghigh [˜ γhigh ] high-in-low · γ A where γ˜high is the density of subsystem A at the higher level of theory, ENN is the nuclear-

nuclear repulsion energy, and Glow and Ghigh are the two electron energies at the low and high levels of theory, respectively:

xc Glow(high) [γ] = J[γ] + Elow(high) [γ]

(15)

xc where J is the Coulomb energy, and Elow(high) is the XC energy from the low (high) method.

We also take advantage of the absolute localization scheme, 33 which restricts the movement of the subsystem’s electrons to only the Bloch functions constructed from orbitals that are centered on the nuclei of that subsystem. This type of localization restricts the electron density to the basis functions of each subsystem, and electrons cannot transfer between the subsystems; however, electrons may still be delocalized over multiple unit cells, as each subsystem is still periodic and described using Bloch functions. This is distinctly different from the real space unit cell localization that we will later perform for cluster-in-periodic embed-

11



Page 12 of 52

ding, which further prevents electrons from being delocalized over multiple unit cells. This reduced basis description of the active subsystem significantly reduces the computational cost of any subsequent high level method calculation on that subsystem. For embedding using absolute localization, the low level-in-low level DFT embedding is no longer exactly equal to the canonical KS-DFT solution, and the high level-in-low level energy may be corrected as corrected A Ehigh-in-low [˜ γhigh , γ˜ A , γ˜ B , γ˜ KS-DFT ] uncorrected A = Ehigh-in-low [˜ γhigh , γ˜ A , γ˜ B ]

(16)

+ E˜low [˜ γ KS-DFT ] − E˜low [˜ γ A + γ˜ B ] where γ˜ KS-DFT is the density of the canonical KS-DFT solution, and E˜low is the KS-DFT total energy functional of the periodic system at the lower level of theory. In the case where the embedding is exact—that is, γ˜ A + γ˜ B = γ˜ KS-DFT —this additional correction is numerically zero. In the absolute localization scheme, however, this correction is not negligible as electrons are forbidden to move between subsystems; it therefore becomes necessary to solve the canonical KS-DFT on the full system at the lower level of theory in order to obtain γ˜ KS-DFT .

Embedding a Molecular Cluster in a Periodic Environment We also present a molecular cluster-in-periodic environment embedding method (see Figure 1). In cluster-in-periodic embedding, the subsystem of interest (subsystem A) is described as a molecular (zero dimensional) system in real space, embedded in the potential from the periodic environment (subsystem B) described in reciprocal space. We achieve this by first obtaining the total Fock matrix of a periodic subsystem A in reciprocal space, Fourier transforming this into real space for the origin unit cell, then subtracting out the two electron potentials of the (also Fourier transformed) density of subsystem A—thereby obtaining an embedding potential for a molecular subsystem A in real space. This Fourier 12


Page 13 of 52 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


transformation can be thought of as localizing the electrons to the real space unit cells, and each molecular orbital is now expanded using atom centered Gaussian-type orbitals instead of Bloch functions; this “unit cell localization” is different from, and should not be confused with, the absolute localization method mentioned earlier. We then construct an embedding potential for one of these localized real space unit cells (the origin unit cell), which can then be added to any DFT or WF method with a molecular description of subsystem A. The cluster-in-periodic embedding workflow is shown in Scheme 2, where we start with the converged low level-in-low level periodic-in-periodic solutions for the densities of subsystems A and B. To transform and obtain an embedding potential of subsystem A in real space, we first construct the Fock matrix of subsystem A embedded in B at the lower level of theory ˜ A-in-B ) as in Equation 11. This is then Fourier transformed into real space for the origin (F unit cell from a set of discrete k-points using

FA-in-B [˜ γ A , γ˜ B ] = wk

Nk X

˜ A-in-B [˜ F γ A , γ˜ B ; k]

(17)

k

where Nk is the number of k-points and wk is the weight of each k-point; our method is implemented into PySCF, which uses a Monkhorst-like sampling of the k-space grid where all k-points are equally weighted. 86 We perform a similar transformation to obtain the real space density matrix of subsystem A γ A from its reciprocal space analogue γ˜ A (k). This transformed Fock matrix may then be diagonalized to obtain a new γ A . In order to construct a new Fock matrix in k-space, the new γ A is Fourier transformed back into reciprocal space; since we only have a description of γ A in the origin cell, we will use the transformation γ˜ A (k) = γ A for all k. This process, shown as “SCF III” in Scheme 2, is repeated until γ A is self-consistent. We should point out that one may also update the density of the periodic subsystem B γ˜ B (k), and therefore ensure that both γ A and γ˜ B (k) are self-consistent. However we find that additional update of γ˜ B (k) in the cluster-in-periodic algorithm is unstable and leads to convergence issues, so we keep γ˜ B (k) fixed at the periodic-in-periodic self-consistent 13



Start Converged 𝜸 ෥ 𝐴 (𝒌) & 𝜸 ෥𝐵 𝒌 𝜸 ෥𝑡𝑜𝑡 𝒌 = 𝜸 ෥𝐴 𝒌 + 𝜸 ෥𝐵 𝒌 ෩ 𝑡𝑜𝑡 𝒌 Construct 𝑭

SCF III

෩ 𝐴-𝑖𝑛-𝐵 (𝒌) from 𝑭 ෩ 𝑡𝑜𝑡 𝒌 Get 𝑭 𝐵 ෡ 𝐵 (𝒌) using 𝜸 ෥ 𝒌 to construct 𝑷 ෩ 𝐴-𝑖𝑛-𝐵 𝒌 → 𝑭𝐴-𝑖𝑛-𝐵 Transform 𝑭 Diagonalize and update 𝜸𝐴 Transform 𝜸𝐴 → 𝜸 ෥ 𝐴 (𝒌)

𝜸𝐴 converged?

No

Yes Get

𝐴-𝑖𝑛-𝐵 𝑭ℎ𝑖𝑔ℎ -𝑖𝑛-𝑙𝑜𝑤

SCF IV

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

𝐴 Diagonalize and update 𝜸ℎ𝑖𝑔ℎ

𝐴 𝜸ℎ𝑖𝑔ℎ converged?

No

Yes End

Scheme 2: The finite cluster-in-periodic embedding workflow. We begin by using the converged periodic densities of subsystems A and B at the lower level of theory (see Scheme 1). We first converge subsystem A in real space at the lower level of theory. After convergence, we use the same embedding potential for subsystem A at the higher level of theory.

14


Page 14 of 52

Page 15 of 52 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


solution. Furthermore, when the cluster approximation can be accurately made—that is, ˜ A-in-B (k) → FA-in-B and γ˜ A (k) → γ A transformations without loss when we can perform the F of information—then there is no need to update γ˜ B (k), and “SCF III” is quickly converged; in such cases, the cluster-in-periodic solution is exactly equal to the periodic-in-periodic solution for low level-in-low level DFT embedding. For a cluster high level-in-periodic low level embedding, including WF-in-DFT embedding, we construct the cluster embedded core Hamiltonian for subsystem A at the low level of theory as

A-in-B A hlow [γ , γ˜ A , γ˜ B ] = FA-in-B [˜ γ A , γ˜ B ] − glow [γ A ]

(18)

where glow is the two electron potentials at the lower level of theory for the finite molecular system. The obtained (uncorrected) cluster high level-in-periodic low level energy becomes uncorrected A Ehigh-in-low [γhigh , γ A , γ˜ A , γ˜ B ]

h

i

A-in-B A ˜ low [˜ = Tr hlow · γhigh +G γ A , γ˜ B ]

(19)

A − Glow [γ A ] + Ghigh [γhigh ] + ENN

where Ghigh and Glow are the two electron energies of the finite molecular subsystem at the high and low levels of theory, respectively. In the case of absolute localization, this energy is corrected as corrected A Ehigh-in-low [γhigh , γ A , γ˜ A , γ˜ B , γ˜ KS-DFT ] uncorrected A = Ehigh-in-low [γhigh , γ A , γ˜ A , γ˜ B ]

(20)

+ E˜low [˜ γ KS-DFT ] − E˜low [γ A , γ˜ B ] Note that the last term in this equation evaluates the DFT-in-DFT energy of the periodic system using the finite cluster γ A that has been Fourier transformed back into reciprocal space.

15



Page 16 of 52

Extending to Correlated Wave Function Methods A-in-B in Eq. 18 describes the potential on a finite cluster The embedded core Hamiltonian hlow

due to a periodic DFT environment. We can therefore easily apply this embedding potential A-in-B ˆ A-in-B ; note that hlow already to any WF method to obtain an embedded WF Hamiltonian H

contains the electron kinetic potential and nuclear potential—that is, the core Hamiltonian— of the electrons and nuclei in subsystem A. The (uncorrected) cluster WF-in-periodic DFT energy is thus uncorrected [ΨA , γ A , γ˜ A , γ˜ B ] EWF-in-DFT

ˆ A-in-B |ΨA i = hΨA |H

(21)

˜ low [˜ +G γ A , γ˜ B ] − Glow [γ A ] where ΨA is the WF of the embedded subsystem A. In absolute localization, this energy is corrected as corrected EWF-in-DFT [ΨA , γ A , γ˜ A , γ˜ B , γ˜ KS-DFT ] uncorrected = EWF-in-DFT [ΨA , γ A , γ˜ A , γ˜ B ]

(22)

˜ low [˜ ˜ low [˜ +G γ KS-DFT ] − G γ A , γ˜ B ] where we use the canonical KS-DFT two electron energy of the full periodic system, rather than the two electron energies from the sum of periodic subsystems A and B. Since the projection operator is added to the Hamiltonian of subsystem A, all resulting HF orbitals (occupied or virtual) strive to be orthogonal to the occupied orbitals of subsystem B—with an energy penalty added to the eigenvalues of the orbitals that are not orthogonal. In the case where subsystem A is described using the full system basis, there will be exactly n—where n is the number of occupied orbitals in subsystem B—”fictitious” virtual orbitals that are not orthogonal; all other HF orbitals remain perfectly orthogonal to all the occupied orbitals of subsystem B. For the parameter-dependent operator of Manby, Miller and co-

16


Page 17 of 52 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


workers, 25 the large energy penalty will prevent excitations into these fictitious orbitals in a subsequent WF calculation. Excitations into these fictitious virtual orbitals, which are actually the occupied orbitals of subsystem B, are allowed with the Huzinaga operator, though there is an associated (smaller) energy penalty. In our previous work, we found that this did not lead to larger errors in WF-in-DFT calculations. 33 Additionally, in the case where subsystem A is described using the monomolecular basis—which is used in this paper—it is often mathematically impossible to produce HF orbitals that are perfectly orthogonal to all the occupied orbitals of subsystem B. As such, all orbitals have some (infinitesimal) overlap and energy penalty. The subsequent WF calculation will similarly result in natural orbitals that have some (infinitesimal) overlap with the occupied orbitals of subsystem B; though, once again, the projection operator strives to prevent exciting into orbitals that are occupied (by subsystem B) via the energy penalty. As observed in Ref. 33, the Huzinaga operator—and in particular the Fermi-shifted Huzinaga presented in this work—is attractive for subsystem calculations in this monomolecular basis since the energy penalty for nonorthogonality (which is unavoidable in this monomolecular basis) is on the same order as the eigenvalue of the orbitals themselves.

Computational Details Our embedding methods were implemented into a modified version of the python-based simulations of chemistry framework (PySCF) 94 version 1.4 beta. 95 PySCF uses a dual real space and reciprocal space cubic grid representation of the crystalline orbitals in order to calculate the Coulomb potentials, 86,96 and we used grids with a spacing of 0.33 ˚ A for all systems with only first row elements, and a spacing of 0.25 ˚ A for all systems that included a second row element; these values ensured that our all-electron energies were converged to within 10−6 Eh with respect to grid spacing. We used the Fermi-shifted Huzinaga operator (Equation 8) for all systems with positive Fermi energies (He9 , Diamond, Silicon, and NaCl), and the

17



Page 18 of 52

Table 1: Summary of the computational details for each of the systems examined in this paper. System

Dimensions

k-point grid

Basis

He9 Polyethylene Neoprene Neoprene-2 Diamond Silicon NaCl CO/Si(100)

2 1 1 1 3 3 3 2

10×10×1 6×1×1 4×1×1 2×1×1 2×2×2 2×2×2 2×2×2 1×1×1

6-31G 87 6-311G* 88,89 6-311G* 6-311G* 6-311G* 6-311G(d) 90 8-511G (Na); 91 86-311G (Cl) 92 6-311G(d) (Si); cc-pVTZ (H, C, O) 93

traditional Huzinaga operator for all other systems, unless stated otherwise. Other system specific computational details are given in Table 1, including number of dimensions with periodic boundary conditions, k-point grids, and the basis sets used. Calculations were carried out using DFT with either the LDA, 97,98 PBE 99 or B3LYP 100,101 XC functionals, or using coupled cluster with singles and doubles (CCSD), 102 or with singles, doubles, and perturbative triples (CCSD(T)) 103 WF methods, where stated. Periodic benchmark calculations for the hybrid B3LYP functional were done in the Crystal computational chemistry software, 104 while geometry optimizations for the CO/Si(100) systems were performed in the Vienna abinitio simulation package (VASP). 105 All geometries, including unit cell lattice parameters and subsystem divisions, are included in the Supporting Information.

Results and Discussion Exact Periodic-in-Periodic DFT Embedding Projection-based embedding uses level shift projection operators to change the eigenvalues of the occupied orbitals of an “inactive” subsystem, from the perspective of the “active” subsystem, such that they are above the Fermi energy and would therefore be unoccupied. In their seminal paper, 25 Manby, Miller, and coworkers used a parameter-dependent level 18


Page 19 of 52 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


shift projection operator. For periodic systems, this µ-parameter operator is

˜ B (k) = µS ˜ AB (k)˜ ˜ BA (k) P γ B (k)S

(23)

where µ is some large (positive) number. This operator always shifts the occupied orbitals’ eigenvalues by a constant µ, and therefore µ should be sufficiently large so as to shift the occupied eigenvalues above the Fermi energy; for very large µ (read limit µ → ∞), this method is independent of the value of the Fermi energy. However, the Fock matrix with this level shift projection operator does not commute with the density operator of the “inactive” subsystem—except when each subsystem is described using the full system basis or when µ = ∞—and therefore it is not always possible to obtain a set of self-consistent orbitals for both subsystems. An alternative, which is both independent of user-defined parameters and commutes with the density operator of the inactive subsystem, is the Huzinaga operator. 78,85 For periodic systems, this operator may be defined as ˜ B (k) = P 1 ˜ AB ˜ BA (k) F (k)˜ γ B (k)S − 2

(24)

˜ AB (k)˜ ˜ BA (k) + S γ B (k)F

This operator works by changing the sign of the projected occupied crystalline orbitals’ eigenvalues. If the subsystems have zero or negative Fermi energies, the Huzinaga operator changes the sign of the eigenvalues of the projected occupied orbitals from negative to positive, and moves these orbitals above the Fermi energy, causing them to be unoccupied. However, it is possible for a subsystem to have a positive Fermi energy, with positive orbital eigenvalues. For these subsystems, the Huzinaga projection operator will fail (results not shown), as it shifts some projected orbitals to lower energies by changing their eigenvalues from positive to negative.

19



Our presented Fermi-shifted Huzinaga operator in Equation 8 seeks to resolve this positive Fermi energy issue, while still having the desirable features of the traditional Huzinaga operator. This operator works by flipping the value of the projected occupied orbitals about the Fermi energy F , instead of about “zero energy” as in the traditional Huzinaga operator, thereby always shifting the projected orbitals above the Fermi energy and into the space of unoccupied orbitals. F is taken as the maximum Fermi energy because, in general, the Fermi energy is subsystem dependent, and the operator only commutes with the density operator if ˜ A and P ˜ B . In all of our the same value of F is used in both subsystem projection operators P calculations, we define the Fermi energy as the average energy between the highest occupied and lowest unoccupied crystalline orbitals. This modification of the Huzinaga operator can be thought of as a “fix” for systems with unphysical positive occupied eigenvalues; in the limit that the Fermi energy is zero, the Fermi-shifted Huzinaga operator is identical to the traditional Huzinaga operator. We further recommend the use of this operator for projectionbased embedding in molecular systems, in place of the traditional Huzinaga operator, where it may also be possible to obtain positive eigenvalues for occupied orbitals; a similar energyshifted Huzinaga operator will ensure that these positive eigenvalue occupied orbitals are always shifted to higher energies. Projection-based periodic-in-periodic DFT-in-DFT embedding using this Fermi-shifted Huzinaga projection operator is exact; that is, DFT-in-DFT, with both subsystems and their interaction described using the same XC functional, is identical to the canonical KSDFT on the full system. To show this, we examine the energies and densities obtained from a series of test cases. These test cases include two one dimensional systems—polyethylene and neoprene—and two three dimensional systems—diamond and silicon. The unit cells of these test systems are shown in Figure 2, and their coordinates are included in the Supporting Information. In the images shown in Figure 2, the larger atoms describe one subsystem (subsystem A), while the smaller atoms describe another subsystem (subsystem B). The division of these four systems into subsystems are thus: a) CH2 /CH2 ; b) CHCCl/C2 H4 ; c)

20


Page 20 of 52

Page 21 of 52 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


a) Polyethylene

b) Neoprene

c) Diamond

d) Silicon

Figure 2: The unit cells of four test systems: a) Polyethylene; b) Neoprene; c) Diamond; and d) Silicon. Systems a and b are one dimensional, and c and d are three dimensional. The larger atoms represent subsystem A, and the smaller atoms are subsystem B.

21



Page 22 of 52

C2 /C6 ; and d) Si/Si7 . This particular division of each of the test systems allows us to have subsystem unit cells that are both charge neutral and closed-shell. To show that projectionbased embedding is exact for periodic systems, each subsystem will be described using the basis of the full system. Table 2: Absolute energy difference (|∆E|) in nano-Hartrees (nEh ), and the number of displaced electrons (h∆ρi) of the four test systems with periodic-in-periodic LDA-in-LDA, PBE-in-PBE and B3LYP-in-B3LYP embedding. System Polyethylene Neoprene Diamond Silicon

LDA-in-LDA |∆E| (nEh ) h∆ρi ×10−6 0.001 0.014 0.000 0.009

1 5 1 11

PBE-in-PBE |∆E| (nEh ) h∆ρi ×10−6 0.003 0.009 0.000 0.005

B3LYP-in-B3LYP |∆E| (nEh ) h∆ρi ×10−6

1 4 1 11

0.002 0.007 0.000 0.002

1 5 1 13

Evidence of the exactness of the projection-based embedding method is shown in Table 2, where both the energy and density errors from the subsystem method are presented, and “method 1”-in-“method 2” describe the methods used for subsystems A (method 1) and B (method 2), respectively. In this table, we include the data for periodic-in-periodic LDA-inLDA, PBE-in-PBE, and B3LYP-in-B3LYP—with errors defined with respect to canonical LDA, PBE, and B3LYP, respectively—to show that projection-based embedding is similarly accurate for local, semi-local, and hybrid functionals. We observe errors on the order of 10−11 Eh or better for all four systems examined. To compare the densities of the embedding method to that of the canonical KS-DFT in Table 2, we use the number of displaced electrons h∆ρi, which is defined as 106–109 1 Z A ρ (r) + ρB (r) − ρKS-DFT (r) dr h∆ρi = 2 Cell

(25)

where ρKS-DFT , ρA and ρB are the densities of the canonical KS-DFT system, and of subsystems A and B, respectively. Integration was performed over the unit cell using a cubic grid with a grid spacing of 0.1 ˚ A. We observed h∆ρi values on the order of 10−5 or less for all systems in Table 2; that is, each of these embedding calculations produced the exact 22


Page 23 of 52 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


canonical KS-DFT densities.

Figure 3: Total canonical KS-DFT and embedding energies (top) and embedding energy error (bottom) of bulk NaCl as a function of the lattice constant using both the µ-parameter and Fermi-shifted Huzinaga operators. The inset shows the total unit cell; this system is divided into two equal (Na2 Cl2 ) subsystems, shown as the larger or smaller atoms.

Finally, we show that our projection-based embedding method is exact by examining the lattice constant of bulk (three dimensional) NaCl; these results are shown in Figure 3. This system has been examined before by Carter and coworkers, 14 where the authors found that their optimized effective potential (OEP) embedding method accurately, to within 10−3 Eh , reproduced the canonical KS-DFT results; the small errors were attributed to the use of a penalty function in the OEP process. In our examination, we divided the NaCl unit cell into two subsystems, each containing two Na and two Cl atoms. These two subsystems are shown using the larger or smaller atom models in the inset in Figure 3. Once again, this division of the total system allowed us to maintain charge neutral, closed-shell unit cells for each subsystem. The resulting embedded LDA-in-LDA energies—with both the µ-parameter operator and the Fermi-shifted Huzinaga operator—exactly reproduced the canonical KSDFT LDA energies, with energy errors (shown in the lower panel of Figure 3) that are ∼10−7 23



Eh for the µ-parameter operator (using a µ value of 106 Eh ), and ∼ 10−12 Eh (numerically zero) for the Fermi-shifted Huzinaga operator. Before we end this section, we should note that our current implementation of projectionbased embedding can only exactly project out fully occupied crystalline orbitals. For systems with many states near the Fermi energy, such as in conductors, convergence with integer occupations is difficult, and one often resorts to smearing the occupation of these states. While the various level shift projection operators in their current implementation may still be used for such fractionally occupied orbitals, we may not necessarily obtain the exact canonical KS-DFT results for two reasons: (1) fractionally occupied orbitals eigenvalues are not projected/shifted to the same extent as fully occupied orbitals; and (2) subsystem calculations often result in crystalline band energies that are different from the canonical KS-DFT, even though their total energies are identical, and therefore any smearing of electron distributions will result in different band occupancies in the canonical and subsystem calculations. We are therefore only able to reproduce exact KS-DFT results for semiconductors and insulators using our current implementation. Further study is needed to ascertain the applicability of these level shift projection operators for embedding in metallic systems.

Absolute Localization in Periodic-in-Periodic Embedding So far, we have only described projection-based embedding where both subsystems are described using the same level of theory, and where there is no computational cost saving compared to canonical KS-DFT as the subsystems are still described in the basis of the full system. In ref. 33, we significantly reduced computational cost and maximize error cancellation by restricting the electrons to the basis functions centered on the atoms of their respective subsystems—the subsystems are described in the monomolecular basis. The total energy of such an embedding calculation is then corrected with the canonical KS-DFT energy of the full system as in Equation 16. We referred to these two strategies—the use of the monomolecular basis, and the correction of the total energy with that of the KS-DFT 24


Page 24 of 52

Page 25 of 52 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


energy—as the “absolute localization” approach. This approach significantly reduces the cost of WF-in-DFT calculations, since the WF region is only described by a subset of the full system basis. Similar cost savings may be obtained for hybrid DFT-in-local (or semi-local) DFT, as the exchange operator for the absolutely localized embedded hybrid region spans only a subset of the total basis. To take advantage of this in periodic-in-periodic embedding, we will similarly apply this absolutely localized technique by restricting each subsystem’s electrons to the Bloch functions constructed from atomic orbitals centered on the nuclei of that subsystem only. This type of localization restricts the location of the electrons to within their respective subsystems only, but the electrons are still delocalized over all periodic images of the unit cell as they are described by Bloch functions. Later in this paper, we will also localize the electrons to a single unit cell in our “finite cluster-in-periodic environment” embedding method. Another attractive reason for using subsystems described in the monomolecular basis is that we obtain a unique embedding solution. Consider the total density γ tot that minimizes the total energy and is obtained from a sum of subsystem densities A and B, or γ tot = γ A + γ B . In general, when each subsystem density matrix is described in the full system basis, the solutions of γ A and γ B are not unique, as there may exist many γ A’ = γ A + δ and γ B’ = γ B − δ—where each of γ A’ and γ B’ integrates to the correct number of electrons in that subsystem—that yield the same γ tot . In such cases, the final solution of γ A and γ B will be determined by the initial guesses for these densities, which would have to be carefully chosen. However, in the monomolecular basis, γ tot may be written as 

 A

γ  γ tot =   0

0  γ

B

 

(26)

This results in a γ tot that is block diagonal, and there exists only one γ A and γ B that yields the same γ tot . Of course, it is impossible to obtain the exact canonical KS-DFT density using this block diagonal form; this is why the embedding energies are corrected with the 25



canonical KS-DFT energies, as shown in equation 16, in the absolute localization method.

Figure 4: Embedding energy error, compared to canonical KS-DFT, for a) polyethylene and b) neoprene in the monomolecular basis with the traditional Huzinaga and Fermi-shifted Huzinaga operators. In Figure 4, we show the error introduced in the embedding method by restricting the electrons of each subsystem to the monomolecular basis. In the monomolecular basis, projectionbased embedding is no longer exact—that is, we are not guaranteed to obtain two sets of perfectly orthogonal occupied orbitals—and the results depend on the choice of projection operator. The Huzinaga operator generally results in an order of magnitude smaller errors than the results from the µ-parameter operator (results shown in Figure S1 in the Supporting Information). As discussed before, 33 the µ-parameter operator was never intended to be used in the monomolecular basis, and the error in this basis was attributed to the high penalty cost for non-orthogonality. Figure 4 shows that the Fermi-shifted Huzinaga operator, described in equation 8, results in marginally smaller errors than those obtained from the traditional Huzinaga operator, with an average error of 0.34 Eh with the Fermi-shifted Huzinaga operator compared to 0.37 Eh with the traditional Huzinaga operator. In general, we do not expect to see a significant difference between the traditional Huzinaga operator 26


Page 26 of 52

Page 27 of 52 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


and the Fermi-shifted Huzinaga operator for systems with negative Fermi energies. The energy errors presented in this figure are what is being corrected in equation 16. (II)

(IV)

(III)

(V)

Figure 5: Dissociation energy curves of one of the C-Cl bonds in neoprene-2. Curves are shown for canonical KS-DFT with PBE and LDA, and for PBE-in-LDA embedding for four different divisions of the total systems labeled (II), (III), (IV), and (V), with the larger atoms in the inset describing the “active” subsystem. The bond dissociation energies are given in the legend. While the absolute energy errors shown in Figure 4 are large, we observe much smaller errors when looking at energy differences using the absolute localization approach. To show this, we used a unit cell of neoprene that has two C-Cl bonds, labeled neoprene-2, and examined the bond dissociation energy curve of one of the C-Cl bonds. These results are presented in Figure 5, where we show the embedded PBE-in-LDA binding energies with absolute localization for four subsystem partitions, as well as the canonical PBE and LDA binding energies. The four subsystem partitions include either 2, 3, 4, or 5—labeled (II), (III), (IV), or (V), respectively—carbon atoms, and all other (non-carbon) atoms bonded 27



to them, in the “active” subsystem; this active subsystem is shown as the larger atoms in the inset in the figure. We also show the bond energies for each of these methods in the figure legend. Using canonical PBE as the benchmark, we observe that LDA overestimates the bond energy by 10.1 kcal/mol, but this error is drastically reduced to 2.6 kcal/mol for PBE-in-LDA subsystem division (II)—which includes only four atoms (two carbons, one chlorine, and one hydrogen) in the active region. Intuitively, the errors are reduced as we include more atoms in the active region, and we observe less than 1 kcal/mol errors with subsystem partitions (IV) and (V). With our current implementation, the periodic low level-in-low level self-consistent cycles (“SCF I” in Scheme 1) are difficult to converge when the subsystems are described in the monomolecular basis. Typically, convergence is straight forward in the supermolecular basis, and we often employ the direct inversion of the iterative subspace (DIIS) 110 to speed this up. However, with the monomolecular basis we have to resort to mixing of the current Fock matrix with that of the previous iteration, often with factors as high as 0.9, in order to obtain converged results. Convergence with this mixing, and with DIIS turned off, is very slow, with convergence taking as much as 500 iterations for the systems examined in this paper in order to obtain energies and densities that are converged to 10−8 Eh and 10−4 atomic units, respectively. However, we are converging to a well-defined energy minimum; therefore, in princple fast convergence should be obtainable. We expect that more advanced convergence techniques, such as the ones specifically designed for subsystem DFT, 111 will lead to significantly faster convergence and we are currently exploring these options. We do not see it as a major issue preventing the widespread application of this method, particularly for finite cluster-in-periodic DFT embedding and WF-in-periodic DFT embedding described in the subsequent sections.

28


Page 28 of 52

Page 29 of 52 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


Finite Cluster DFT-in-Periodic DFT Embedding We will now explore our molecular cluster-in-periodic environment embedding method, where the active subsystem (subsystem A) is described by a finite or zero-dimensional molecular cluster, and it is embedded in an environment (subsystem B) with periodic boundary conditions. Such a method is particularly useful if we wish to describe the active subsystem using, for example, a hybrid XC functional where the expensive exchange operator in periodic boundary conditions may be computationally intractable. This finite cluster approximation also makes it possible to describe the active subsystem using any standard WF method, as we will explore in the next section. Before we show the results for our cluster-in-periodic DFT embedding method, we will first discuss its limitations. The periodic-to-cluster transformation in equation 17 is effectively a weighted average of the Fock matrix over all k-points. This transformation of the active subsystem is therefore only “exact”, or we do not lose information about the subsystem, when the Fock matrix is identical at all k-points—that is, the transformation is only accurate when the subsystem does not significantly interact with its periodic images and its band energies are more discrete or molecule-like. This suggests that the application of the cluster-in-periodic embedding method is limited to active regions that are molecule-like, such as studying the interactions of molecules with bulk material. A similar idea has been previously explored by Genova and Pavanello for subsystem DFT, 109 where they explored reducing the k-point sampling in molecule-like subsystems. Our cluster-in-periodic embedding takes this one step further by removing all periodic boundary conditions for the active subsystem, which is then described in real space rather than in reciprocal space. While the cluster approximation can only be made for subsystems with molecule-like band energies, we find that the periodic-in-periodic embedding method naturally results in subsystems with molecule-like band energies—even for systems with rich band structures; that is to say, our periodic-in-periodic embedding forces localization. Consider the two dimensional He9 model system, as shown in the inset in Figure 6, where the atoms are 29



Figure 6: Band energies of the occupied bands of the He9 system from canonical KS-DFT (black solids), subsystem A (blue dashes) and subsystem B (orange dash-dots) from periodicin-periodic embedding, and subsystem A (green dots) from cluster-in-periodic embedding. Inset shows the atoms described as subsystem A (orange) or subsystem B (blue). evenly spaced by a distance of 1.5 ˚ A. We divide the total system such that subsystem A contains the central He atom, while subsystem B describes the remaining eight atoms. Figure 6 shows the band energies of the occupied crystalline orbitals from canonical LDA KS-DFT (solid black lines), which shows a rich band structure with no particular band being molecule-like. However, once we perform a periodic-in-periodic embedding—which returns the exact canonical KS-DFT total energy—we observe that the band energy for subsystem A (dash-dotted orange lines) is molecule-like with a band width of only 0.3 eV, while subsystem B (dashed blue lines) retains much of the band structure from the canonical KS-DFT results. This makes a cluster approximation of subsystem A possible, with little loss of accuracy; the orbital energy from an embedded cluster subsystem A (dotted green lines) is a good approximation to the band energy of the periodic subsystem A, and the cluster-in-periodic method returns an error in the total energy of only 2×10−3 Eh compared to canonical KS-DFT. Periodic-in-periodic embedding is thus able to return molecule-like 30


Page 30 of 52

Page 31 of 52 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


bands for subsystem A as it forces the crystalline orbitals of that subsystem to be localized, even though both subsystems in Figure 6 are still described in the basis of the full system. Another way to interpret this forced localization is that the presence of subsystem B screens subsystem A from significantly interacting with its periodic images, resulting in narrow band widths for subsystem A. Ideally, we will also be using the monomolecular basis description of each subsystem to save on computational cost in the absolute localization approach, which also helps to localize the orbitals. The cluster approximation therefore further benefits from absolute localization, and we observe an energy difference of 2×10−5 Eh for cluster-inperiodic embedding compared to periodic-in-periodic embedding of the He9 system using the monomolecular basis. We stress that the current implementation of periodic-in-periodic and cluster-in-periodic embedding are ground state embedding methods, and therefore should not be used for obtaining band structures; band structures are only used here to show that we can accurately approximate subsystem A as a molecular cluster. To further demonstrate the accuracy of the absolute localization approach in finite clusterin-periodic embedding, we once again examine the C-Cl bond dissociation in neoprene-2 using the hybrid B3LYP XC functional as the benchmark. These results are shown in Figure 7, where we use subsystem partitions (II) and (III) from our earlier discussion (in Figure 5), these partitions are shown as insets in the figure; we also give the bond dissociation energies in the figure legend. The methods labeled as B3LYP, LDA, and PBE are the canonical KSDFT with periodic boundary conditions, while B3LYP-in-LDA and B3LYP-in-PBE are the cluster-in-periodic methods, where subsystem A is a B3LYP molecular cluster described only in the basis of the atoms of the active region, and subsystem B is a carried out using periodic boundary conditions. Therefore, these embedding methods are considerably cheaper than the benchmark canonical B3LYP calculation with periodic boundary conditions. The B3LYP dissociation energy for this C-Cl bond is 70.7 kcal/mol; LDA over binds by 7.1 kcal/mol and PBE under binds by 3.0 kcal/mol. In each of the B3LYP-in-LDA (III), B3LYP-in-PBE (II), and B3LYP-in-PBE (III) methods, we observe bond energies that were within 1 kcal/mol of

31



(II)

Page 32 of 52

(III)

Figure 7: Top: Cluster-in-periodic bond dissociation energy curves for the C-Cl bond in neoprene-2. The subsystem partitions are shown in the insets. Bottom: Bond dissociation energy as a function of the subsystem partition. “None” refers to the canonical KS-DFT dissociation energy.

32


Page 33 of 52 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 B3LYP benchmark. These results show that we can obtain the correct high level energies, once a large enough active region is selected, regardless of what low level XC method we use; although PBE, in this case, is a better low level approximation as we were able to obtain accurate energies with both subsystem partitions (II) and (III), compared to only obtaining accurate energies with subsystem partition (III) with LDA. We should point out that subsystem partitions (IV) and (V) do not yield converged cluster-in-periodic results. This is because the periodic-in-periodic embedding did not sufficiently localize the crystalline orbitals for the active subsystem A; in other words, the sizes of subsystem A for these two partitions are too large, relative to the size of the unit cell, that they are not sufficiently screened from interacting with their periodic images by subsystem B. A back-of-the-envelope method to determine when a subsystem might be too large to result in molecule-like band energies is to examine the maximum of the overlap of the real space AO basis functions of subsystem A in one particular unit cell and that of its immediate neighboring cell. For the (II) and (III) subsystem partitions, the maximum overlaps in the AO functions of subsystem A are on the order of 10−7 and 10−5 , respectively, but this increases to 10−3 and 10−2 for the (IV) and (V) partitions, respectively. For such larger active subsystems, it will be necessary to extend the size of the unit cell in order to perform a cluster-in-periodic embedding. Further study is needed to elucidate exactly when we can or cannot obtain molecule-like bands for the active subsystem from periodic-in-periodic embedding.

Correlated Wave Function-in-Periodic DFT Embedding Having set up the necessary framework in the preceding sections, we now demonstrate our WF-in-periodic DFT embedding method. In WF-in-periodic DFT embedding, we use the same embedding potential as in the finite cluster-in-periodic DFT embedding method to modify the core Hamiltonian of the WF region. This method is therefore similar to the WF-in-DFT embedding with absolute localization in molecular systems, 33 with the only 33



exception being that the embedding potential is due to an environment that is described by periodic boundary conditions. In practice, we first perform a Hartree-Fock (HF)-in-periodic DFT calculation using the transformed periodic embedding potential—that is, an embedded Hartree-Fock calculation for “SCF IV” in Scheme 2—followed by a CCSD(T) calculation using the HF orbitals and using the same core Hamiltonian as outlined in Equation 18. Note that in the current implementation, we do not polarize the periodic DFT environment to the (new) density of the WF subsystem.

(II)

(III)

Figure 8: Top: WF-in-periodic DFT bond dissociation energy curves for the C-Cl bond in neoprene-2. The subsystem partitions are shown in the insets. Bottom: Bond dissociation energy as a function of the subsystem partition. “None” refers to the canonical KS-DFT dissociation energy. We begin by showing that the WF-in-periodic DFT embedding energies converge with the 34


Page 34 of 52

Page 35 of 52 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


size of the active subsystem; these results are shown in Figure 8. We once again explore the bond dissociation of neoprene-2 using subsystem partitions (II) and (III) that were explored in Figure 5; these subsystem partitions are shown as insets in the figure. The WF subsystems were described using CCSD, and the DFT subsystems were described using either the LDA or PBE XC functionals. We do not have a reference periodic CCSD benchmark for this system, but we do observe a convergence of the bond dissociation energy with respect to subsystem partition size (Figure 8 Bottom); CCSD-in-LDA and CCSD-in-PBE yield consistent results and are only within 1.3 kcal/mol for subsystem partition (III), compared to a difference more than 10 kcal/mol between the canonical LDA and PBE results, again showing that the DFT approximation does not significantly influence the WF-in-DFT energies. This convergence is similar to that observed for the B3LYP-in-LDA and B3LYP-in-PBE cluster-in-periodic embedding methods in Figure 7.

Figure 9: Unit cells used in examining the WF-in-periodic DFT binding of CO on to a Si(100)-2×1 surface. The larger atoms represent the WF subsystem. In applying this method to a more realistic system, we look at the binding energy of CO on to the Si(100)-2×1 surface. This system is representative of a typical molecular

35



system interacting with a periodic environment, and is therefore easily studied using our cluster-in-periodic and WF-in-periodic DFT embedding methods. The unit cells for this system are shown in Figure 9, where we include the CO molecule and one Si atom in the WF subsystem, presented as the larger atom models, and the remainder is denoted as the DFT subsystem. The H, C, and O atoms are all described using the cc-pVTZ basis, 93 while the Si atoms are described using the 6-311G(d) basis. 90 The system is periodic in two dimensions with lattice constants (7.680624, 7.680624 ˚ A), and the geometries are given in the Supporting Information. We divide the total system such that the CO molecule and one Si atom make up the WF subsystem; therefore, the WF calculation is only performed over these three atoms and its cost is negligible compared to the preceding DFT-in-DFT and canonical KS-DFT on the full periodic system. Like in our absolute localization WF-in-DFT method for molecular systems, 33 we find that the computational costs for WF-in-periodic DFT embedding are dominated by the DFT-in-DFT and KS-DFT on the full system, and therefore we are paying DFT-like costs for WF-like accuracy. Table 3: Absorption energy of CO on Si(100)-2×1 (in kcal/mol) from canonical KS-DFT and CCSD(T)-in-periodic DFT embedding. Method

Energy (Error)

Experiment 112 LDA PBE CCSD(T)-in-LDA CCSD(T)-in-PBE

11.6 32.1 20.7 13.0 10.3

(20.5) (9.1) (1.4) (1.3)

The CO absorption energies are given in Table 3, where the experimental binding energy of 11.6 kcal/mol was taken from Ref. 112. The canonical LDA and PBE methods are inconsistent, and overestimate the binding energy by 20.5 and 9.1 kcal/mol, respectively. However, both CCSD(T)-in-LDA and CCSD(T)-in-PBE are within 1.5 kcal/mol of the experimental value, and represent significant corrections of 19.1 and 10.4 kcal/mol, respectively, to the canonical KS-DFT values. These WF-in-DFT results show, as was observed in Figures 7 and 36


Page 36 of 52

Page 37 of 52 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


8, that our energies are independent of the choice and accuracy of the XC functional used. The consistent energies obtained by CCSD(T)-in-LDA and CCSD(T)-in-PBE—along with obtaining experimentally accurate energies—suggests that we have selected a sufficiently large active subsystem, even though the active region is made up of only three atoms. These results therefore demonstrate that we can obtain accurate energies for reactions involving bulk systems at DFT-like computational costs using projection-based WF-in-periodic DFT embedding.

Summary and Conclusion We have described a periodic-in-periodic DFT embedding method, a finite cluster-in-periodic DFT embedding method, and a WF-in-periodic DFT embedding method; each method is successively built upon the preceding method. Periodic-in-periodic embedding uses level shift projection operators, where we also proposed an alternative to the Huzinaga operator that solves the positive Fermi energy problem, and exactly reproduces canonical KS-DFT energies in the limit that all subsystems are described in the basis of the full system. To save on computational costs, one can limit the electrons of a subsystem to the basis functions of that subsystem only and correct the total embedding energy with the canonical KS-DFT energy—the so called absolute localization approach—that benefits from error cancellation and accurately reproduces high level energies from a high level-in-low level calculation. We can then take the result from a low level periodic-in-periodic embedding, transform the embedded Fock matrix—that also contains the level shift projection operator—into real space in order to obtain an embedding potential for a zero-dimensional molecular cluster. This embedding potential can then be applied to a more computationally expensive DFT method, as in finite cluster-in-periodic DFT embedding, or to a WF method, as in WF-inperiodic DFT embedding. Our embedding methodology is efficient and highly accurate, therefore we expect signifi-

37



cant adoption of this method for the study of semiconductors, insulators, and a wide variety chemical reactions involving bulk systems. Due to the use of electron occupation smearing in conductors, our method will not reproduce KS-DFT calculations on these systems. However, further study is needed to determine if our method is accurate for energy differences in these systems, a subject of our on-going research. Additionally, we are exploring the use of more advanced convergence methods for further computational efficiency. In conclusion, we showed that projection-based embedding is exact for periodic systems using polyethylene, neoprene, diamond, silicon, and NaCl as model test cases. We further demonstrated—by examining the C-Cl bond dissociation of a model neoprene-2 system—that we can obtain accurate high level-in-low level DFT energies with both periodic-in-periodic and cluster-in-periodic embedding and using the absolute localization approach. Finally, we presented accurate—within 1.5 kcal/mol of the experimental value—WF-in-DFT binding energies of CO on to the Si(100)-2×1 surface using CCSD(T)-in-LDA and CCSD(T)-in-PBE embedding.

Acknowledgement This research was carried out within the Nanoporous Materials Genome Center, which is supported by the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Award DE-FG02-12ER16362. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota, and the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, for providing resources that contributed to the results reported within this paper.

38


Page 38 of 52

Page 39 of 52 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


Supporting Information Available The code used in this work is available at https://github.com/Goodpaster/pbpe.git. The embedding energy errors in the monomolecular basis for all three level shift projection operators, and geometry data for all systems are included in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864– B871. (2) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (3) Cortona, P. Self-consistently determined properties of solids without band-structure calculations. Phys. Rev. B Condens. Matter 1991, 44, 8454–8458. (4) Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050–8053. (5) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. WIREs Comput Mol Sci 2014, 4, 325–362. (6) 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. (7) Götz, A. W.; Beyhan, S. M.; Visscher, L. Performance of Kinetic Energy Functionals for Interaction Energies in a Subsystem Formulation of Density Functional Theory. J. Chem. Theory Comput. 2009, 5, 3161–3174.

39



(8) Xia, J.; Huang, C.; Shin, I.; Carter, E. A. Can orbital-free density functional theory simulate molecules? J. Chem. Phys. 2012, 136, 084102. (9) 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–2509. (10) Jacob, C. R. Unambiguous optimization of effective potentials in finite basis sets. J. Chem. Phys. 2011, 135, 244102. (11) Becke, A. D.; Johnson, E. R. A simple effective potential for exchange. J. Chem. Phys. 2006, 124, 221101. (12) Staroverov, V. N.; Scuseria, G. E.; Davidson, E. R. Optimized effective potentials yielding Hartree-Fock energies and densities. J. Chem. Phys. 2006, 124, 141103. (13) Huang, C.; Pavone, M.; Carter, E. A. Quantum mechanical embedding theory based on a unique embedding potential. J. Chem. Phys. 2011, 134, 154110. (14) Huang, C.; Carter, E. A. Potential-functional embedding theory for molecules and materials. J. Chem. Phys. 2011, 135, 194104. (15) Yu, K.; Carter, E. A. Extending density functional embedding theory for covalently bonded systems. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, E10861–E10870. (16) Ding, F.; Tsuchiya, T.; Manby, F. R.; Miller, T. F., 3rd Linear-Response TimeDependent Embedded Mean-Field Theory. J. Chem. Theory Comput. 2017, 13, 4216– 4227. (17) Ding, F.; Manby, F. R.; Miller, T. F., 3rd Embedded Mean-Field Theory with BlockOrthogonalized Partitioning. J. Chem. Theory Comput. 2017, 13, 1605–1615.

40


Page 40 of 52

Page 41 of 52 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


(18) Miyamoto, K.; Miller, T. F., 3rd; Manby, F. R. Fock-Matrix Corrections in Density Functional Theory and Use in Embedded Mean-Field Theory. J. Chem. Theory Comput. 2016, 12, 5811–5822. (19) Fornace, M. E.; Lee, J.; Miyamoto, K.; Manby, F. R.; Miller, T. F., 3rd Embedded mean-field theory. J. Chem. Theory Comput. 2015, 11, 568–580. (20) Fornace, M. E.; Lee, J.; Miyamoto, K.; Manby, F. R.; Miller, T. F., 3rd Correction to Embedded Mean-Field Theory. J. Chem. Theory Comput. 2015, 11, 3968. (21) Gutowski, M.; Piela, L. Interpretation of the Hartree-Fock interaction energy between closed-shell systems. Mol. Phys. 1988, 64, 337–355. ˙ (22) Rajchel, L.; Zuchowski, P. S.; Szcze´sniak, M. M.; Chalasi´ nski, G. Derivation of the supermolecular interaction energy from the monomer densities in the density functional theory. Chem. Phys. Lett. 2010, 486, 160–165. (23) Rajchel, L.; Zuchowski, P. S.; Szcze´sniak, M. M.; Chalasi´ nski, G. Density functional theory approach to noncovalent interactions via monomer polarization and Pauli blockade. Phys. Rev. Lett. 2010, 104, 163001. (24) Gritsenko, O. V. Recent Progress in Orbital-free Density Functional Theory; Recent Advances in Computational Chemistry; WORLD SCIENTIFIC, 2012; Vol. 6; pp 355– 365. (25) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F., 3rd A Simple, Exact Density-Functional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564–2568. (26) Hégely, B.; Nagy, P. R.; Ferenczy, G. G.; Kállay, M. Exact density functional and wave function embedding schemes based on orbital localization. J. Chem. Phys. 2016, 145, 064107. 41



(27) Barnes, T. A.; Goodpaster, J. D.; Manby, F. R.; Miller, T. F. Accurate basis set truncation for wavefunction embedding. J. Chem. Phys. 2013, 139, 024103. (28) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller, T. F., 3rd Accurate and systematically improvable density functional theory embedding for correlated wavefunctions. J. Chem. Phys. 2014, 140, 18A507. (29) Bennie, S. J.; Stella, M.; Miller, T. F., 3rd; Manby, F. R. Accelerating wavefunction in density-functional-theory embedding by truncating the active basis set. J. Chem. Phys. 2015, 143, 024105. (30) Pennifold, R. C. R.; Bennie, S. J.; Miller, T. F., 3rd; Manby, F. R. Correcting densitydriven errors in projection-based embedding. J. Chem. Phys. 2017, 146, 084113. (31) Tamukong, P. K.; Khait, Y. G.; Hoffmann, M. R. Density differences in embedding theory with external orbital orthogonality. J. Phys. Chem. A 2014, 118, 9182–9200. (32) Chulhai, D. V.; Jensen, L. Frozen Density Embedding with External Orthogonality in Delocalized Covalent Systems. J. Chem. Theory Comput. 2015, 11, 3080–3088. (33) Chulhai, D. V.; Goodpaster, J. D. Improved Accuracy and Efficiency in Quantum Embedding through Absolute Localization. J. Chem. Theory Comput. 2017, 13, 1503– 1508. (34) Culpitt, T.; Brorsen, K. R.; Hammes-Schiffer, S. Communication: Density functional theory embedding with the orthogonality constrained basis set expansion procedure. J. Chem. Phys. 2017, 146, 211101. (35) Libisch, F.; Marsman, M.; Burgdörfer, J.; Kresse, G. Embedding for bulk systems using localized atomic orbitals. J. Chem. Phys. 2017, 147, 034110. (36) Huo, P.; Uyeda, C.; Goodpaster, J. D.; Peters, J. C.; Miller, T. F. Breaking the

42


Page 42 of 52

Page 43 of 52 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


Correlation between Energy Costs and Kinetic Barriers in Hydrogen Evolution via a Cobalt Pyridine-Diimine-Dioxime Catalyst. ACS Catal. 2016, 6, 6114–6123. (37) Sauer, J. Molecular models in ab initio studies of solids and surfaces: from ionic crystals and semiconductors to catalysts. Chem. Rev. 1989, 89, 199–255. (38) Whitten, J. L.; Yang, H. Theory of chemisorption and reactions on metal surfaces. Surf. Sci. Rep. 1996, 24, 55–124. (39) Jug, K.; Bredow, T. Models for the treatment of crystalline solids and surfaces. J. Comput. Chem. 2004, 25, 1551–1567. (40) Huang, P.; Carter, E. A. Advances in correlated electronic structure methods for solids, surfaces, and nanostructures. Annu. Rev. Phys. Chem. 2008, 59, 261–290. (41) Winter, N. W.; Pitzer, R. M.; Temple, D. K. Theoretical study of a Cu+ ion impurity in a NaF host. J. Chem. Phys. 1987, 86, 3549–3556. (42) Barandiarán, Z.; Seijo, L. The abinitio model potential representation of the crystalline environment. Theoretical study of the local distortion on NaCl: Cu+. J. Chem. Phys. 1988, 89, 5739–5746. (43) Sousa, C.; Casanovas, J.; Rubio, J.; Illas, F. Madelung fields from optimized point charges for ab initio cluster model calculations on ionic systems. J. Comput. Chem. 1993, 14, 680–684. (44) Dungsrikaew, V.; Limtrakul, J.; Hermansson, K.; Probst, M. Comparison of methods for point-charge representation of electrostatic fields. Int. J. Quantum Chem. 2004, 96, 17–22. (45) Vail, J. M. Theory of electronic defects: Applications to MgO and alkali halides. J. Phys. Chem. Solids 1990, 51, 589–607.

43



(46) Lin, H.; Truhlar, D. G. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185. (47) Dick, B. G.; Overhauser, A. W. Theory of the Dielectric Constants of Alkali Halide Crystals. Phys. Rev. 1958, 112, 90–103. (48) Pascual, J. L.; Seijo, L. Ab initio model potential embedded cluster calculations including lattice relaxation and polarization: Local distortions on Mn2+-doped CaF2. J. Chem. Phys. 1995, 102, 5368–5376. (49) Nasluzov, V. A.; Rivanenkov, V. V.; Gordienko, A. B.; Neyman, K. M.; Birkenheuer, U.; Rösch, N. Cluster embedding in an elastic polarizable environment: density functional study of Pd atoms adsorbed at oxygen vacancies of MgO (001). J. Chem. Phys. 2001, 115, 8157–8171. (50) Sauer, J.; Sierka, M. Combining quantum mechanics and interatomic potential functions in ab initio studies of extended systems. J. Comput. Chem. 2000, 21, 1470–1493. (51) Shoemaker, J. R.;

Gordon, M. S. SIMOMM: An Integrated Molecular Or-

bital/Molecular Mechanics Optimization Scheme for Surfaces. J. Phys. Chem. A 1999, 103, 3245–3251. (52) Nasluzov, V. A.; Ivanova, E. A.; Shor, A. M.; Vayssilov, G. N.; Birkenheuer, U.; Rösch, N. Elastic Polarizable Environment Cluster Embedding Approach for Covalent Oxides and Zeolites Based on a Density Functional Method. J. Phys. Chem. B 2003, 107, 2228–2241. (53) Liu, Y.; Lu, G.; Chen, Z.; Kioussis, N. An improved QM/MM approach for metals. Modell. Simul. Mater. Sci. Eng. 2007, 15, 275. (54) Pisani, C. Approach to the embedding problem in chemisorption in a self-consistentfield-molecular-orbital formalism. Phys. Rev. B Condens. Matter 1978, 17, 3143–3153. 44


Page 44 of 52

Page 45 of 52 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


(55) Pisani, C.; Dovesi, R.; Carosso, P. Moderately-large-embedded-cluster approach to the study of local defects in solids. Vacancy and substitutional impurities in graphite. Phys. Rev. B Condens. Matter 1979, 20, 5345–5357. (56) Kruger, S.; Rosch, N. The moderately-large-embedded-cluster method for metal surfaces; a density-functional study of atomic adsorption. J. Phys. Condens. Matter 1994, 6, 8149. (57) Inglesfield, J. E. Embedding at surfaces. Comput. Phys. Commun. 2001, 137, 89–107. (58) Bormet, J.; Neugebauer, J.; Scheffler, M. Chemical trends and bonding mechanisms for isloated adsorbates on Al(111). Phys. Rev. B Condens. Matter 1994, 49, 17242– 17252. (59) Whitten, J. L.; Yang, H. Theoretical studies of surface reactions on metals. Int. J. Quantum Chem. 1995, 56, 41–47. (60) 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. (61) Govind, N.; Wang, Y. A.; Carter, E. A. Electronic-structure calculations by firstprinciples density-based embedding of explicitly correlated systems. J. Chem. Phys. 1999, 110, 7677–7688. (62) Kl¨ uner, T.; Govind, N.; Wang, Y. A.; Carter, E. A. Prediction of electronic excited states of adsorbates on metal surfaces from first principles. Phys. Rev. Lett. 2001, 86, 5954–5957. (63) Kl¨ uner, 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. 2001, 116, 42–54. 45



(64) 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. (65) 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. (66) 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. (67) Kanan, D. K.; Sharifzadeh, S.; Carter, E. A. Quantum mechanical modeling of electronic excitations in metal oxides: Magnesia as a prototype. Chem. Phys. Lett. 2012, 519–520, 18–24. (68) 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. Zeitschrift f¨ ur Physikalische Chemie 2013, 3, 1455– 1466. (69) Huang, C.; Libisch, F.; Peng, Q.; Carter, E. A. Time-dependent potential-functional embedding theory. J. Chem. Phys. 2014, 140, 124113. (70) Libisch, F.; Huang, C.; Carter, E. A. Embedded correlated wavefunction schemes: theory and applications. Acc. Chem. Res. 2014, 47, 2768–2775. (71) Yu, K.; Libisch, F.; Carter, E. A. Implementation of density functional embedding theory within the projector-augmented-wave method and applications to semiconductor defect states. J. Chem. Phys. 2015, 143, 102806.

46


Page 46 of 52

Page 47 of 52 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


(72) Yu, K.; Krauter, C. M.; Dieterich, J. M.; Carter, E. A. In Fragmentation: Toward Accurate Calculations on Complex Molecular Systems; Gordon, M. S., Ed.; John Wiley & Sons, Ltd, 2017; pp 81–117. (73) Cheng, J.; Yu, K.; Libisch, F.; Dieterich, J. M.; Carter, E. A. Potential Functional Embedding Theory at the Correlated Wave Function Level. 2. Error Sources and Performance Tests. J. Chem. Theory Comput. 2017, 13, 1081–1093. (74) Cheng, J.; Libisch, F.; Yu, K.; Chen, M.; Dieterich, J. M.; Carter, E. A. Potential Functional Embedding Theory at the Correlated Wave Function Level. 1. Mixed Basis Set Embedding. J. Chem. Theory Comput. 2017, 13, 1067–1080. (75) Huang, P.; Carter, E. A. Local electronic structure around a single Kondo impurity. Nano Lett. 2006, 6, 1146–1150. (76) Lykos, P. G.; Parr, R. G. On the Pi-Electron Approximation and Its Possible Refinement. J. Chem. Phys. 1956, 24, 1166–1173. (77) Phillips, J. C.; Kleinman, L. New Method for Calculating Wave Functions in Crystals and Molecules. Phys. Rev. 1959, 116, 287–294. (78) Huzinaga, S.; Cantu, A. A. Theory of Separability of Many-Electron Systems. J. Chem. Phys. 1971, 55, 5543–5549. (79) Stoll, H.; Paulus, B.; Fulde, P. On the accuracy of correlation-energy expansions in terms of local increments. J. Chem. Phys. 2005, 123, 144108. (80) Henderson, T. M. Embedding wave function theory in density functional theory. J. Chem. Phys. 2006, 125, 014105. (81) Mata, R. A.; Werner, H.-J.; Sch¨ utz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106.

47



(82) Chulhai, D. V.; Jensen, L. External Orthogonality in Subsystem Time-Dependent Density Functional Theory. Phys. Chem. Chem. Phys. 2016, 18, 21032–21039. (83) Tamukong, P. K.; Khait, Y. G.; Hoffmann, M. R. Accurate Dissociation of Chemical Bonds Using DFT-in-DFT Embedding Theory with External Orbital Orthogonality. J. Phys. Chem. A 2017, 121, 256–264. (84) Shimazaki, T.; Kitaura, K.; Fedorov, D. G.; Nakajima, T. Group molecular orbital approach to solve the Huzinaga subsystem self-consistent-field equations. J. Chem. Phys. 2017, 146, 084109. (85) Francisco, E.; Mart´ın Pendás, A.; Adams, W. H. Generalized Huzinaga building-block equations for nonorthogonal electronic groups: Relation to the Adams–Gilbert theory. J. Chem. Phys. 1992, 97, 6504–6508. (86) McClain, J.; Sun, Q.; Chan, G. K.-L.; Berkelbach, T. C. Gaussian-Based CoupledCluster Theory for the Ground-State and Band Structure of Solids. J. Chem. Theory Comput. 2017, 13, 1209–1218. (87) 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–2261. (88) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72, 650–654. (89) McLean, A. D.; Chandler, G. S. Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, Z= 11–18. J. Chem. Phys. 1980, 72, 5639–5648. (90) Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. Energy band gaps and lattice

48


Page 48 of 52

Page 49 of 52 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


parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional. J. Chem. Phys. 2005, 123, 174101. (91) Dovesi, R.; Roetti, C.; Freyria-Fava, C.; Prencipe, M.; Saunders, V. R. On the elastic properties of lithium, sodium and potassium oxide. An ab initio study. Chem. Phys. 1991, 156, 11–19. (92) Apra, E.; Causa, M.; Prencipe, M.; Dovesi, R.; Saunders, V. R. On the structural properties of NaCl: an ab initio study of the B1-B2 phase transition. J. Phys. Condens. Matter 1993, 5, 2969. (93) 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–1023. (94) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J.; Sharma, S.; Wouters, S.; Chan, G. K.-L. The Python-based Simulations of Chemistry Framework (PySCF v1.2). arXiv, 2017. (95) Sun, Q. PySCF v1.4 (beta). https://github.com/sunqm/pyscf, 2017; Accessed: 2017-9-22. (96) Sun, Q.; Berkelbach, T. C.; McClain, J. D.; Chan, G. K.-L. Gaussian and plane-wave mixed density fitting for periodic systems. J. Chem. Phys. 2017, 147, 164119. (97) Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Cambridge Philos. Soc. 1930, 26, 376–385. (98) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can. J. Phys. 1980, 58, 1200–1211. (99) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. 49



(100) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (101) 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. JPC J. Planar Chromatogr. - Mod. TLC 1994, 98, 11623–11627. (102) Hirata, S.; Podeszwa, R.; Tobita, M.; Bartlett, R. J. Coupled-cluster singles and doubles for extended systems. J. Chem. Phys. 2004, 120, 2581–2592. (103) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479–483. (104) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; Others, CRYSTAL14: A program for the ab initio investigation of crystalline solids. Int. J. Quantum Chem. 2014, 114, 1287–1317. (105) Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B Condens. Matter 1996, 54, 11169– 11186. (106) Jacob, C. R.; Visscher, L. A subsystem density-functional theory approach for the quantum chemical treatment of proteins. J. Chem. Phys. 2008, 128, 155102. (107) Kiewisch, K.; Eickerling, G.; Reiher, M.; Neugebauer, J. Topological analysis of electron densities from Kohn-Sham and subsystem density functional theory. J. Chem. Phys. 2008, 128, 044114. (108) Genova, A.; Ceresoli, D.; Pavanello, M. Periodic subsystem density-functional theory. J. Chem. Phys. 2014, 141, 174101. 50


Page 50 of 52

Page 51 of 52 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


(109) Genova, A.; Pavanello, M. Exploiting the locality of periodic subsystem densityfunctional theory: efficient sampling of the Brillouin zone. J. Phys. Condens. Matter 2015, 27, 495501. (110) Pulay, P. Convergence acceleration of iterative sequences. the case of scf iteration. Chem. Phys. Lett. 1980, 73, 393–398. (111) Genova, A.; Ceresoli, D.; Krishtal, A.; Andreussi, O.; DiStasio, R. A.; Pavanello, M. eQE: An open-source density functional embedding theory code for the condensed phase. Int. J. Quantum Chem. 2017, 117 . (112) Nakai, H.; Katouda, M.; Kawamura, Y. Energy density analysis of cluster size dependence of surface-molecule interactions: H 2, C 2 H 2, C 2 H 4, and CO adsorption onto Si (100)-(2× 1) surface. J. Chem. Phys. 2004, 121, 4893–4900.

51



Graphical TOC Entry periodic DFT-in-DFT cluster approximation WF-in-periodic DFT

52


Page 52 of 52

Projection-Based Correlated Wave Function in Density Functional

Recommend Documents