Subspace density matrix functional embedding theory: Theory

Subscriber access provided by University of South Dakota

Quantum Electronic Structure

Subspace density matrix functional embedding theory: Theory, implementation, and applications to molecular systems Xing Zhang, and Emily A. Carter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00990 • Publication Date (Web): 21 Dec 2018 Downloaded from http://pubs.acs.org on December 22, 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.

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

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

Journal of Chemical Theory and Computation

Subspace density matrix functional embedding theory: Theory, implementation, and applications to molecular systems Xing Zhang1 and Emily A. Carter2, * 1

Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ

08544-5263 2

School of Engineering and Applied Science, Princeton University, Princeton, NJ 08544-5263

We introduce the subspace density matrix functional embedding theory (sDMFET), in which optimization of the nonlocal embedding potential and subsequent embedded correlated wavefunction calculations are carried out within a truncated subspace determined by a Schmidt decomposition. As compared to the original density matrix functional embedding theory [K. Yu and E. A. Carter, Proceedings of the National Academy of Sciences, 114, E10861 (2017)], the computational cost of sDMFET is significantly reduced while the accuracy is preserved. We perform test calculations for both covalently and non-covalently bound molecular systems to demonstrate the feasibility of our theory.

ACS Paragon Plus Environment

1

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

Page 2 of 41

1 Introduction Kohn-Sham density functional theory (KS DFT) is a powerful method for determining the electronic structures of large molecules and materials. However, the commonly used approximate density functionals are often insufficient to provide accurate results in situations such as electronic excitation, charge transfer, dispersion interactions, etc. Correlated wavefunction (CW) methods, in contrast to DFT approximations, are more accurate, although their high computational cost limits their application to small numbers of atoms. In quantum embedding theories, especially embedded CW (ECW) methods,1–3 one tries to combine the efficiency of DFT and the accuracy of CW methods. This is achieved by partitioning a system into a small region of chemical/physical interest and its surrounding environment, with the former and the latter treated at CW and DFT levels, respectively. A variety of quantum embedding methods exist, e.g., frozen density embedding4–6 with kinetic energy density functionals7–12 or with KS inversions,13–17 projectionbased embedding,18–22 and density matrix embedding theory.23,24 Besides these, a series of KS inversion embedding schemes based on a unique embedding potential shared by each subsystem have been developed in our group,25,26 namely, density functional embedding theory (DFET),27 potential functional embedding theory (PFET),28–31 and density matrix functional embedding theory (DMFET).32 In our methods, the density/density matrix partitioning is uniquely determined, which is not usually the case in other embedding theories. In this work, we focus on improving DMFET, in which a non-local embedding potential is introduced to describe the interactions between the different subsystems. Such a formalism is found to be particularly useful for the situations where subsystem boundaries cut across covalent bonds. Although DMFET is accurate in predicting chemical reaction energies, it is computationally expensive, because both the embedding potential and the subsequent energy


2

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


calculations are performed in the total system basis. However, the high computational cost of DMFET can be avoided by exploiting the following two facts. First, the density matrices are built from the occupied orbitals, so the embedding potential can be computed without any loss of exactness in a truncated subspace spanned by only the occupied orbitals. Second, for chemical reactions that take place in a spatially confined region, a large portion of the system is chemically inert such that the electron density associated with it can be treated as frozen, which may further reduce the size of the subspace. In order to construct such a subspace, we perform a Schmidt decomposition33 for the KS determinant of the total system. The original DMFET is then reformulated within the subspace. This newly developed method, which we call subspace DMFET (sDMFET), is a simple approximation to DMFET in which the embedding potential optimization and local refinement of the electronic structure are carried out in the subspace. Although sDMFET generally yields different density partitions, embedding potentials and energy corrections than DMFET, we will show that sDMFET is as accurate as DMFET while also being much more efficient. The remainder of the paper is organized as follows. In Section 2, we briefly introduce the formalism of DMFET, followed by providing the definition of the subspace in sDMFET based on a Schmidt decomposition. In Section 3, we examine the performance of sDMFET by means of benchmark calculations for systems ranging from hydrogen-bonded molecular clusters to those that are covalently bonded. We also apply sDMFET to study a surface defect in a large Si nanocrystal. Thereafter, we compare computational costs of sDMFET and DMFET. In Section 4, we state our conclusions.

2 Theory


3


Page 4 of 41

2.1 DMFET In DMFET,32 the one-electron density matrix of the total system 𝛾𝑡𝑜𝑡, computed at the DFT level, typically is divided into two pieces, 𝛾𝐴 and 𝛾𝐵. (Although it is straightforward to generalize DMFET to multiple subsystems, we only consider two subsystems for simplicity.) These pieces then are assigned to two subsystems A (embedded region) and B (environment), respectively: 𝛾𝑡𝑜𝑡(𝒓,𝒓′) = 𝛾𝐴(𝒓,𝒓′) + 𝛾𝐵(𝒓,𝒓′).

(1)

Here, the subsystem density matrices are obtained by solving the embedded KS equations:

( ― 12∇2 + 𝑣𝑋𝑒𝑥𝑡 + 𝑣𝑋𝐻𝑋𝐶[𝜌𝑋] + 𝑣𝑒𝑚𝑏)𝜓𝑋𝑖(𝒓) = 𝜀𝑋𝑖𝜓𝑋𝑖(𝒓),

(2)

where X labels subsystems A and B, 𝑣𝑋𝐻𝑋𝐶 denotes the Hartree and exchange-correlation (XC) potentials, the 𝜓’s and 𝜀’s are the KS orbitals and their respective energies, and the subsystem density is defined as 𝜌𝑋(𝒓) ≡ 𝛾𝑋(𝒓,𝒓′)|𝒓 = 𝒓′ = ∑𝑖𝜓𝑋𝑖 ∗ (𝒓)𝜓𝑋𝑖(𝒓′)|𝒓 = 𝒓′.

(3)

In Eq. (2), the embedding potential is assumed to have a non-local form 𝑣𝑒𝑚𝑏(𝒓,𝒓′) and is the same for both subsystems.27,32 It can be proven that unique solutions for the embedding potential and the density matrix partitioning can be obtained32 once the electron number and the external potential 𝑣𝑋𝑒𝑥𝑡(𝒓) for each subsystem are fixed. To solve for the embedding potential, one can perform a constrained minimization for the sum of the subsystem energies 𝐸𝐴 and 𝐸𝐵, with the constraint shown in Eq. (1). In practice, the constrained minimization is recast into an unconstrained maximization,27,32,34 with the objective functional defined as 𝑊[𝑣𝑒𝑚𝑏] = 𝐸𝐴[𝜌𝐴;𝑣𝑒𝑚𝑏] + 𝐸𝐵[𝜌𝐵;𝑣𝑒𝑚𝑏] + ∬𝑑𝒓𝑑𝒓′𝑣𝑒𝑚𝑏 × (𝛾𝐴 + 𝛾𝐵 ― 𝛾𝑡𝑜𝑡),

(4)

where the embedding potential appears as a Lagrange multiplier. DMFET is currently implemented in atom-centered basis sets, since it is prohibitively expensive to optimize the nonlocal embedding potential in real space or using planewave basis.32


4

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


After one obtains the embedding potential, a CW calculation may be carried out to refine the electronic structure of the embedded region. For example, the total-system DFT energy can be improved by a first-order local energy correction such that27 𝐷𝐹𝑇 𝑒𝐶𝑊 ), 𝐸𝐶𝑊/𝐷𝐹𝑇 ― 𝐸𝑒𝐷𝐹𝑇 𝐷𝑀𝐹𝐸𝑇 = 𝐸𝑡𝑜𝑡 + (𝐸𝐴 𝐴

(5)

where the letter “e” in the superscripts indicates that the energies are obtained in the presence of the embedding potential.27 It should be noted that the local electronic structure refinement is performed only for the electrons in the embedded region. In other words, the number of electrons contained in the embedded region is set a priori.

2.2 sDMFET and Schmidt decomposition The computational cost of the original DMFET is one-to-two orders of magnitude higher than that of the regular DFT calculation on the total system,32 because maximization of the functional W in Eq. (4) is performed in the total-system basis set. Furthermore, the subsequent ECW calculation is carried out in the same basis set, which prevents it from being applied to large systems. Nevertheless, as we will show next, the high computational cost can be avoided by performing the DMFET calculation within a truncated space. Effectively, DMFET can be viewed as a molecular orbital (MO) localization method that rotates the occupied MOs of the total system to two new sets of MOs localized on each of the subsystems, respectively.32 It is clear that the virtual MOs do not contribute to the functional W and the concomitant optimized embedding potential. As such, it is equivalent to performing the DMFET calculation within a subspace spanned by only the occupied MOs of the total system. Such a choice of subspace has two drawbacks, however. First, the ECW calculation requires a set of virtual orbitals that must be determined separately. Second, the possibility that the electron density far


5


Page 6 of 41

away from the embedded region can be kept frozen is not exploited. We overcome these drawbacks here by instead constructing the subspace by the so-called Schmidt decomposition.35 Although other orbital localization schemes such as Boys36 and Pipek−Mezey37 may also be applied to construct the subspace, we find using Schmidt decomposition more robust in our case. First, with Schmidt decomposition, there is no need to localize the occupied and virtual orbitals separately. Second, by following the strategy presented next, the subspace is determined unambiguously with its size limited by the size of the embedded region. The Schmidt decomposition of a Slater determinant can be obtained in the following way.38 First, project the canonical MOs (denoted as 𝜓’s) of the total system onto subsystem A:

|𝜓𝐴𝑖⟩ = ∑𝜇 ∈ 𝐴|𝜒𝜇⟩⟨𝜒𝜇│𝜓𝑖⟩,

(6)

where the 𝜒’s denote a set of orthonormal atomic orbitals (AOs). Next, a new set of natural orbitals localized on subsystem A can be obtained as

|𝐴𝑖⟩ =

1 𝑁𝑜𝑐𝑐 ∑ 𝑑𝑖 𝑘

|𝜓𝐴𝑘⟩𝑈𝑘𝑖.

(7)

Here, 𝑁𝑜𝑐𝑐 is the number of occupied MOs, d and U satisfy the eigenvalue equations 𝐌𝐴𝐔 = 𝐔𝐝, and the overlap matrix 𝐌𝐴 has elements 𝑀𝐴𝑖𝑗 = ⟨𝜓𝐴𝑖│𝜓𝐴𝑗⟩. Similarly, the natural orbitals localized on subsystem B can be derived as 1 𝑁 ∑ 𝑜𝑐𝑐 1 ― 𝑑𝑖 𝑘

(8)

|𝜓𝐵𝑖⟩ = (1 ― ∑𝜇 ∈ 𝐴|𝜒𝜇⟩⟨𝜒𝜇|)|𝜓𝑖⟩.

(9)

|𝐵𝑖⟩ =

|𝜓𝐵𝑘⟩𝑈𝑘𝑖,

where

Then, the Slater determinant can be expanded in the basis of natural orbitals |𝐴𝑖⟩ and |𝐵𝑖⟩:

|𝛷0⟩ = ∏𝑁𝑖 𝑜𝑐𝑐( 𝑑𝑖𝑎 † (𝐴𝑖) + 1 ― 𝑑𝑖𝑎 † (𝐵𝑖))|0⟩,

(10)


6

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


where 𝑎 † (𝐴𝑖) and 𝑎 † (𝐵𝑖) are the creation operators associated with orbitals |𝐴𝑖⟩ and |𝐵𝑖⟩, respectively, and |0⟩ denotes the vacuum state. It then becomes obvious from Eq. (10) that the probability of having an electron in orbitals |𝐴𝑖⟩ and |𝐵𝑖⟩ is 𝑑𝑖 and (1 ― 𝑑𝑖), respectively. In practice, it is more straightforward to construct the natural orbitals mentioned above by diagonalizing the reduced density matrices.24,33 To illustrate this, one can write the one-electron density matrix for the Slater determinant in Eq. (10) as 𝑁

𝛾 = ∑𝑖 𝑜𝑐𝑐|𝜓𝑖⟩⟨𝜓𝑖|,

(11)

|𝜓𝑖⟩ = 𝑑𝑖|1𝐴𝑖0𝐵𝑖⟩ + 1 ― 𝑑𝑖|0𝐴𝑖1𝐵𝑖⟩.

(12)

where

In Eq. (12), |0𝐴𝑖⟩ and |1𝐴𝑖⟩ represent empty and fully-occupied |𝐴𝑖⟩ states, respectively, and a similar notation is used for states |𝐵𝑖⟩. By taking the partial trace over states |𝐵𝑖⟩, we obtain the one-electron reduced density matrix on subsystem A as 𝑁

𝛾𝐴 = ∑𝑖 𝑜𝑐𝑐[𝑑𝑖|1𝐴𝑖⟩⟨1𝐴𝑖| + (1 ― 𝑑𝑖)|0𝐴𝑖⟩⟨0𝐴𝑖|].

(13)

As the empty states do not contribute to the density matrix, Eq. (13) is equivalent to 𝑁

𝛾𝐴 = ∑𝑖 𝑜𝑐𝑐𝑑𝑖|𝐴𝑖⟩⟨𝐴𝑖|.

(14)

By using Eqs. (6) and (7), it is easy to show that 𝑁

𝛾𝐴 = ∑𝜇𝜈 ∈ 𝐴∑𝑖 𝑜𝑐𝑐𝐶𝜇𝑖𝐶𝑖𝜈† |𝜒𝜇⟩⟨𝜒𝜈| , = ∑𝜇𝜈 ∈ 𝐴𝑃𝜇𝜈|𝜒𝜇⟩⟨𝜒𝜈|

(15)

where C is the canonical MO coefficient matrix and P is the one-electron density matrix in the AO basis. Because Eq. (15) is just a different representation of Eq. (14), 𝐏𝐴, with its elements defined as 𝑃𝐴𝜇𝜈 = 𝑃𝜇𝜈 ∈ 𝐴, should have the same eigenvalues as d. This is true regardless of whether the dimension


7


Page 8 of 41

of 𝐏𝐴, labeled as 𝑁𝐴, is larger than 𝑁𝑜𝑐𝑐 or not. If 𝑁𝐴 is larger than 𝑁𝑜𝑐𝑐, then 𝐏𝐴 will have

(𝑁𝐴 ― 𝑁𝑜𝑐𝑐) zero eigenvalues in addition to the eigenvalues of d. Otherwise, the upper limit of the summation in Eq. (14) is replaced by 𝑁𝐴. It then is possible to show that the eigenvectors of 𝐏𝐴 associated with eigenvalues 𝑑𝑖 are exactly those natural orbitals |𝐴𝑖⟩. Consider the overlap matrix 𝐌𝐴, which can be rewritten as 𝐌𝐴 = 𝐂𝐴† 𝐂𝐴 = 𝐔𝐝𝐔 † ,

(16)

where 𝐂𝐴 has its elements defined as 𝐶𝐴𝜇𝑝 = 𝐶𝜇 ∈ 𝐴,𝑝 ∈ 𝑜𝑐𝑐. By multiplying 𝐂𝐴 on the left and 𝐂𝐴† on the right of Eq. (16), we have 𝐂𝐴𝐌𝐴𝐂𝐴† = 𝐂𝐴𝐂𝐴† 𝐂𝐴𝐂𝐴† = 𝐂𝐴𝐔𝐝𝐔 † 𝐂𝐴† , = 𝐏𝐴𝐏𝐴 = 𝐕𝐝𝐝𝐕 †

(17)

where V is the eigenvector matrix of 𝐏𝐴. Moreover, Eq. (17) indicates that 𝐕 = 𝐂𝐴𝐔𝐝 ―1/2 ,

(18)

and, combined with Eqs. (6) and (7), we conclude that ∑𝜇 ∈ 𝐴|𝜒𝜇⟩𝑉𝜇𝑖 = |𝐴𝑖⟩.

(19)

Similarly, diagonalizing 𝐏𝐵 gives |𝐵𝑖⟩ with corresponding eigenvalues (1 ― 𝑑𝑖). In sDMFET, the subspace is spanned by the natural orbitals |𝐴𝑖⟩ and |𝐵𝑖⟩ with the occupation numbers 0 < 𝑑𝑖 ≤ 1 and 0 < 1 ― 𝑑𝑖 < 1, respectively. In practical calculations, we identify pairs of entangled |𝐴𝑖⟩ and |𝐵𝑖⟩ orbitals by examining their occupation numbers, where the sum should be 1. (The accuracy of such an identification procedure depends on the error introduced by the eigenvalue solver, which in our case is on the order of machine precision.) The rest of the natural orbitals are unentangled and remain either empty or fully occupied. Those pairs of entangled orbitals then are sorted in descending order by the extent of entanglement (the smaller the difference between 𝑑𝑖 and 1 ― 𝑑𝑖 is, the larger the entanglement is). If there is any |𝐴𝑖⟩ orbital with


8

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


occupation number 𝑑𝑖~1 that cannot find a pairing |𝐵𝑖⟩ orbital (i.e., |𝐴𝑖⟩ is fully occupied and unentangled), then this usually means that the embedded region is larger than the environment region and that all the entangled orbitals are important for the subsequent ECW calculation on the embedded region. In this case, we include all pairs of entangled orbitals in the subspace. Otherwise, only the most entangled pairs of |𝐴𝑖⟩ and |𝐵𝑖⟩ orbitals are included in the subspace, truncated at the orbital pair where 𝑑𝑖 has its maximum value. This pair of orbitals, according to Eq. (12), corresponds to a total-system occupied MO that has the smallest overlap with the environment. The remaining entangled orbitals |𝐴𝑖⟩ (|𝐵𝑖⟩), treated approximately as empty (fully occupied), are not included in the subspace, because they have occupation numbers close to 0 (1), which represent the total-system occupied MOs that are almost entirely localized in the environment. Such an approximation is sufficiently accurate for the systems studied in this work: the entangled orbitals

|𝐴𝑖⟩ that are excluded from the subspace have occupation numbers less than 10 ―6. We also include the fully occupied |𝐴𝑖⟩ orbitals in the subspace for simplicity, but treat those fully occupied |𝐵𝑖⟩ orbitals as if they are frozen. Finally, for the ECW calculation, the subspace is expanded further by adding in the empty |𝐴𝑖⟩ orbitals, and the embedding potential must be projected onto this enlarged space. To summarize, the size of the subspace for embedding potential optimization is the number of occupied |𝐴𝑖⟩ orbitals plus the number of partially occupied |𝐵𝑖⟩ orbitals; while for ECW calculations, it is the number of basis functions in the embedded region plus the number of partially occupied |𝐵𝑖⟩ orbitals.

2.3 Orthogonality between subsystem MOs After the functional W in Eq. (4) is maximized, one obtains the optimized embedding potential along with the subsystem density matrices 𝛾𝐴 and 𝛾𝐵 that satisfy the constraint in Eq. (1). It can be


9


Page 10 of 41

proven that the occupied subsystem MOs used to construct these density matrices are mutually orthogonal (see Appendix). However, the virtual MOs of one subsystem are in general nonorthogonal to the MOs of the other subsystem. Since the virtual orbitals of the embedded region have contributions from the environment occupied orbitals, the subsequent ECW calculations on the embedded region will introduce effectively the electron excitations to the environment occupied orbitals, thereby violating the Pauli Exclusion Principle. To solve this problem, we project out the contributions of the environment occupied orbitals from the virtual orbital space of the embedded region. We do this by introducing a level shifting operator, so that the modified embedding potential applied to the ECW calculation now reads 𝑣𝑒𝑚𝑏 = 𝑣𝑒𝑚𝑏 +𝜆∑𝑖 ∈ 𝑜𝑐𝑐|𝜓𝐵𝑖⟩⟨𝜓𝐵𝑖|,

(20)

where 𝜓𝐵 denotes the environment canonical MOs (represented in the basis of natural orbitals introduced in Section 2.2), and 𝜆 is set to be a large number (e.g., 105). The modification to the embedding potential only affects the virtual orbital space of the embedded region; hence, it should not alter the density partitioning that satisfies Eq. (1). Note the discussion here applies to both sDMFET and DMFET.

3 Results and discussions In this section, we examine the performance of sDMFET and compare it with that of the original DMFET. All of the calculations were performed using a modified version of the PySCF program.39 A Python script was implemented to determine the subspace and to optimize the embedding potential using the algorithm given in Ref. [40]. We also apply Fermi-Dirac smearing, with a smearing width of 0.005 hartree to stabilize the embedded DFT calculations during optimization.40 ―5 The optimizations were considered to be converged when max |𝐷𝐴𝑖𝑗 + 𝐷𝐵𝑖𝑗 ― 𝐷𝑡𝑜𝑡 , where 𝑖𝑗 | < 10 𝑖𝑗


10

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


the D’s are the one-electron density matrices computed within the subspace. If the boundaries between the embedded region and the environment cut across covalent bonds, then we assign the bonding electrons to the embedded region and follow the procedure introduced in Ref. [32] to cap the subsystems with point charges. Since DMFET is not very sensitive to the exact capping sites,32 these point charges are simply placed at the positions of the boundary atoms in the surroundings. We performed benchmark calculations classified into three categories, depending on the different types of boundaries between the embedded region and the environment. The first three examples partition the total system across hydrogen bonds, including a proton transfer reaction of the solvated (FHF) ― anion, binding energies of water hexamers, and low-lying excited states of solvated uracil. Next, we studied deprotonation of decanoic acid, where the boundary between the two subsystems cuts across a single covalent bond. Finally, we examined an example of partitioning conjugated systems, namely, for the hydrogenation of pentacene. For comparison, we also carried out ONIOM41 calculations for some of the tests. The dangling bonds in the ONIOM calculations were capped by hydrogen atoms positioned at a distance of 1.087 Å, which roughly matches the C—H bond lengths in methane and ethylene. In the remainder of this article, we also will treat ONIOM as an embedding method with a zero embedding potential. For all of the embedding calculations, the notation “method 1/method 2” implies that the embedded region and the environment are treated by methods 1 and 2, respectively. Furthermore, we performed nonembedded calculations on bare embedded regions, denoted with a “cap-” prefix followed by the level of theory (e.g., cap-CCSD), using the same capping strategy as that in the ONIOM calculations. Finally, we applied sDMFET to investigate a proposed defect-induced conical intersection of a Si nanocrystal, which illustrates the efficacy of sDMFET for studying properties


11


Page 12 of 41

of isolated defects in semiconductors. The computational cost of sDMFET also will be discussed briefly.

3.1 Proton transfer reaction in solvated (FHF) ― In this section, we investigate the proton transfer reaction in the model system (FHF) ― (H2O)4 (Fig. 1). The geometry of the system is obtained from Ref. [26], with the two F ― anions separated by 3.4 Å and the water molecules relaxed at the MP2/6-31G*42 level. The distance between the proton and the left F ― anion (𝑅F ― H in Fig. 1) is chosen to be the reaction coordinate, and the proton moves along the straight line connecting the two F ― anions. We carried out potential energy surface (PES) scans along the reaction pathway at different levels of theory using the cc-pVDZ43 basis set. The reference PES was computed by coupledcluster singles and doubles plus perturbative triples44 [CCSD(T)] on the entire system. We compare three embedding methods, namely, sDMFET, DMFET, and ONIOM, at the CCSD(T)/Hartree-Fock (HF) and CCSD(T)/DFT-B3LYP45,46 levels, for which the results are plotted in Figs. 2 and 3, respectively. Examination of the energy errors displayed in Figs. 2(b) and 3(b) reveals that sDMFET and DMFET give very similar PESs with overall errors of ~1 kcal/mol compared to the reference CCSD(T) results. Moreover, the use of different mean-field theories (i.e., HF and DFT-B3LYP) to treat the environment has little effect on the quality of the computed PESs. Unlike sDMFET and DMFET, ONIOM is sensitive to the low-level theories employed in the calculations. When the environment is treated at the HF level, ONIOM overestimates the reaction barrier and exothermicity by ~5 kcal/mol [Fig. 2(b)]. However, the reaction energy is underestimated by ~2 kcal/mol when DFT-B3LYP is used [Fig. 3(b)].


12

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


H O

H

F RF-H H F embedded region

Figure 1. Molecular structure of the (FHF) ― (H2O)4 complex. (FHF) ― is treated as the embedded region. The distance between the proton and the left F ― anion (𝑅F ― H) is the reaction coordinate.

Figure 2. PESs (a) and relative errors (b) along the proton transfer pathway for the (FHF) ― (H2O)4 complex. All curves have been shifted to match at 𝑅F ― H = 1.0 Å, which corresponds to a local minimum geometry. CCSD(T)/cc-pVDZ energies of the entire system provide the benchmarks for comparison (black). For the embedding calculations, CCSD(T) and HF are employed to treat the embedded region and the environment, respectively. While ONIOM (blue) overestimates the reaction barrier and exothermicity, sDMFET (green) and DMFET (red) energies agree with each other very well and match the benchmarks within ~1 kcal/mol.


13


Page 14 of 41

Figure 3. As with Fig. 2, except that DFT-B3LYP is employed to treat the environment in the embedding calculations, and here ONIOM underestimates the exothermicity.

3.2 Binding energies of water hexamers Following Manby et al.,18 we test sDMFET for its ability to capture many-body effects through embedding. As an example, we compute the binding energies of water hexamer isomers with the embedded many-body expansion truncated at the two-body level (EMBE2).18 The structures of those (H2O)6 isomers, optimized at the MP2/aug′-cc-pVTZ level,47 were obtained from Ref. [47]. Since our sDMFET method requires a total-system energy calculation at the mean-field level, the many-body expansion (MBE) can be performed only for the correlation energy:18 𝑒𝐶𝑜𝑟𝑟 ― (𝑁 ― 2)∑𝑖𝐸𝑒𝐶𝑜𝑟𝑟 = 𝐸𝐻𝐹 𝐸𝐸𝑀𝐵𝐸2 ], 𝑖 𝑡𝑜𝑡 𝑡𝑜𝑡 + [∑𝑖 < 𝑗𝐸𝑖𝑗

(21)

where 𝐸𝑒𝐶𝑜𝑟𝑟 and 𝐸𝑒𝐶𝑜𝑟𝑟 are the correlation energies of the embedded monomers and dimers, 𝑖 𝑖𝑗 respectively, and N is the number of monomers. For comparison, the energy expression of the nonembedded many-body expansion at the two-body level (MBE2) reads:48 𝐶𝑊 = ∑𝑖 < 𝑗𝐸𝐶𝑊 𝐸𝑀𝐵𝐸2 𝑖𝑗 ― (𝑁 ― 2)∑𝑖𝐸𝑖 , 𝑡𝑜𝑡

(22)


14

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


where 𝐸𝐶𝑊 and 𝐸𝐶𝑊 𝑖 𝑖𝑗 label the monomer and dimer energies, respectively, calculated by the CW methods. In principle, the total-system HF energy (𝐸𝐻𝐹 𝑡𝑜𝑡) in Eq. (21) also should be expanded at the two-body level in EMBE2,48 i.e., 𝑒𝐻𝐹 𝑒𝐻𝐹 . 𝐸𝐻𝐹 𝑡𝑜𝑡 = ∑𝑖 < 𝑗𝐸𝑖𝑗 ― (𝑁 ― 2)∑𝑖𝐸𝑖

(23)

However, only the local correlation energy (𝐸𝑒𝐶𝑊 ― 𝐸𝑒𝐻𝐹 𝐴 𝐴 ) can be evaluated unambiguously due to the fact that the embedding potential determined by sDMFET contains an arbitrary constant shift. As such, we cannot perform the embedded MBE at the mean-field level [as shown in Eq. (23)] within the framework of sDMFET. Fig. 4 displays the binding energies of (H2O)6 isomers computed by MP2 on the entire system, by EMBE2 with the embedding described by sDMFET, by HF-Δ12,18,49 and by MBE2. In the HFΔ12 method, one defines the total-system energy exactly as in Eq. (21), except that a zero embedding potential is used.18,49 Thus the difference between HF-Δ12and MBE2 is the latter’s partitioning of the HF total system energy as in Eq. (23), while the difference between HF-Δ12 and our EMBE2 method is the presence of the nonlocal embedding potential. The aug-cc-pVDZ50 basis set was employed for all calculations. For the calculations on monomers and dimers, the totalsystem basis was applied to avoid basis set superposition error (BSSE). This also requires projection of the embedding potential determined by sDMFET onto the total-system basis for the subsequent ECW calculations. Shown on a finer scale are the errors in our EMBE2 method and the HF-Δ12 method relative to the supersystem MP2 results. As the binding energies computed by HF-Δ12 and MBE2 (green and blue bars in Fig. 4) reveal, the majority of the error in MBE2 comes from the approximation for the total-system HF energies [i.e., Eq. (23)], while using MBE2 (as in HF-Δ12) for the correlation energies alone is more accurate. However, the total HF energies for the (H2O)6 isomers studied here are ~455 hartrees


15


Page 16 of 41

whereas the total MP2 correlation energies are only ~1.3 hartrees. As such, using MBE2 for the correlation energies actually produces higher percentage errors than for the HF energies. Our EMBE2 method, by contrast, agrees with the supersystem MP2 almost perfectly for most of the

(H2O)6 isomers, because the many-body effects are captured correctly through the embedding potentials. The binding energy errors given by our EMBE2 method are also significantly smaller than those from the HF-Δ12 calculations (lower panel of Fig. 4).

Figure 4. Binding energies of (H2O)6 isomers computed by supersystem MP2, EMBE2, HF-Δ12, and MBE2 (upper panel). All calculations employ the aug-cc-pVDZ basis set. Errors in the binding energies of EMBE2 and HF-Δ12 calculations relative to the supersystem MP2 results are shown in the lower panel.

3.3 Solvatochromic shifts in low-lying excited states of uracil


16

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


Solvatochromism refers to changes in solute absorption spectra observed when solutes are immersed in varying solvents.51 An explicit quantum-mechanical treatment of the solvent molecules is usually necessary when “specific” solvent effects are important, e.g., strong hydrogen bonding.52 While it is computationally intractable to include all of the solvent molecules in highlevel CW calculations, our sDMFET method provides a scheme where only the solute molecule is treated at the CW level whereas the solvent effects are captured through the embedding potential computed efficiently at the mean-field level. In this section, we test sDMFET for its ability to describe solvatochromic shifts in uracil, which were shown to have a significant impact on excited-state deactivation mechanisms.53 The model system studied here is a uracil molecule microhydrated by four water molecules [(uracil)(H2O)4]. The ground-state equilibrium geometry is obtained from Ref. [53], which was optimized at the spin-flip time-dependent DFT54 level using the BH&HLYP hybrid exchange-correlation functional (50% HF exchange plus 50% Becke exchange55 with LYP correlation46) and the 631+G**56 basis set. We compute the vertical excitation energies for the first two singlet excited states of both the gas-phase uracil and the (uracil)(H2O)4 complex, using the equation-of-motion coupled-cluster singles and doubles57 (EOM-CCSD) method and the 6-31+G** basis. In Table 1, a solvent-induced blue-shift of ~0.2 eV in the excitation energy for the S1 state, which has a 𝑛𝜋 ∗ character, is predicted due to the hydrogen bonding between uracil and water. The same microhydration environment leads to a red-shift of ~0.2 eV in the excitation energy for the S2 state, which has a 𝜋𝜋 ∗ character. Note that the blue-shift for the 𝑛𝜋 ∗ state is slightly smaller than the previously reported value of ~0.5 eV,53 possibly because long-range bulk polarization by the solvent, which was treated by implicit solvent models in the previous studies,53 are not considered here.


17


Page 18 of 41

In the sDMFET calculations, uracil is taken as the embedded region and the water molecules are defined as the environment. The embedding potentials were computed at either the HF or DFTB3LYP levels. However, the excitation energies, calculated at the embedded EOM-CCSD level, are insensitive to the choice of the mean-field theories used to determine the embedding potentials (Table 1). Our sDMFET method is able to reproduce the total-system EOM-CCSD results for the excitation energies of the 𝑛𝜋 ∗ and 𝜋𝜋 ∗ states with errors of ~0.1 eV and ~0.04 eV, respectively. We also present predicted probability amplitudes of the important orbital transitions for each state. The results from the sDMFET and the total-system EOM-CCSD calculations agree with each other qualitatively, implying that sDMFET can correctly describe transition properties, as well. For the 𝑛𝜋 ∗ state, sDMFET yields excitation energies with larger errors, which can be understood by examining the frontier MOs that have major contributions to the transitions (Fig. 5). From the supersystem calculation for the (uracil)(H2O)4 complex, we find that the nonbonding orbitals involved in the 𝑛→𝜋 ∗ transitions are delocalized over the oxygen atoms of uracil and those of the neighboring water molecules [Figs. 5(a)-(b)]. However, because the number of electrons included in each subsystem is fixed a priori, any delocalization across the subsystems will not be captured in our (s)DMFET method. This is evident by observing the corresponding nonbonding orbitals from the embedded calculation [Figs. 5(e)-(f)], which are localized on uracil. On the other hand, the 𝜋 and 𝜋 ∗ orbitals are well localized on uracil, even in the total-system calculation [Figs. 5(c)(d)]. As such, sDMFET, in the current example, is less accurate for describing the 𝑛𝜋 ∗ state than the 𝜋𝜋 ∗ state. Nevertheless, such problems can be solved easily by including few water molecules that are closest to the solute in the embedded region.


18

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


Table 1. Vertical excitation energies and transition characters for the S1 and S2 states of gas-phase uracil and the (uracil)(H2O)4 complex, computed at the EOM-CCSD level on the entire system, and at the sDMFET level with the embedded region (uracil) treated by EOM-CCSD and the environment (water molecules) treated by either HF or DFT-B3LYP. The 6-31+G** basis set was employed in all calculations. Excitation energies (eV) sDMFETc B3LY Gasa b States phase Hydrated HF P S1 5.30 5.49 5.59 5.61 a

S2 5.80 5.61 5.65 5.65 Gas-phase uracil at the EOM-CCSD level.

b

(uracil)(H2O)4 at the EOM-CCSD level.

c

(uracil)(H2O)4 at the sDMFET level.

Probability amplitude sDMFETc Orbital Gastransitions phasea Hydratedb HF B3LYP n1→π* 0.190 0.349 0.298 0.265 n2→π* 0.662 0.608 0.635 0.561 π→π* 0.747 0.740 0.754 0.667


19


Page 20 of 41

Figure 5. (a)-(d) Frontier MOs of the (uracil)(H2O)4 complex computed at the HF/6-31+G** level. (e)-(h) Frontier MOs of the uracil molecule computed at the embedded HF/6-31+G** level for which the embedding potential is obtained from the sDMFET calculation for the (uracil)(H2O)4 complex at the same level of theory. The isosurface values are set to 0.05 bohr ―3.

3.4 Deprotonation of decanoic acid In this section, we further test the performance of sDMFET by studying the deprotonation energy of decanoic acid, which serves as an example where the boundary between the embedded region and the environment cuts across a single covalent bond. The reference reaction energy (i.e., the difference in energy between the protonated and deprotonated forms in the gas phase) was computed at the total-system CCSD/cc-pVDZ level of theory, using structures optimized at the DFT-B3LYP/cc-pVDZ level. In Fig. 6(a), we plot the errors in the deprotonation energies obtained from different embedding methods and the cap-CCSD method (i.e., capping the dangling bond with a hydrogen atom) as a function of the embedded regions’ size. Figure 6(b) presents the structures of decanoic acid and its deprotonated form, as well as the partitioning of the total system. As expected, the cap-CCSD results exhibit the largest errors. However, since the reaction is localized in the terminal carboxylic acid group, even cap-CCSD can predict the reaction energy with sub-kcal/mol accuracy when the embedded region is extended to include five C atoms. The three embedding methods improves the cap-CCSD results significantly, while our sDMFET and DMFET methods perform better than the ONIOM method. This is most evident in the case where only the carboxylic acid group is treated as the embedded region (i.e., one C atom in the embedded region); both sDMFET and DMFET already are reasonably accurate with reaction energy errors less than 2 kcal/mol, whereas ONIOM overestimates the protonation energy by ~4 kcal/mol.


20

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


Again, we find that sDMFET agrees with DMFET for all of selections of the embedded regions, thereby validating our subspace formalism.

Figure 6. Embedding errors for the deprotonation of decanoic acid, where the boundary between the embedded region and its environment cuts across a single covalent bond. (a) Errors in the reaction energy plotted as a function of the size of the embedded region. (b) The structures of decanoic acid and its deprotonated form, calculated at the DFT-B3LYP level. The cc-pVDZ basis set was employed in all calculations. The reference deprotonation reaction energy, 367.5 kcal/mol, was obtained at the CCSD level. The embedding calculations were performed at the CCSD/HF level. The areas enclosed by the blue and red boxes are examples where the embedded regions contain one and two carbon atoms, respectively. Along the horizontal axis of the plot in (a), the


21


Page 22 of 41

embedded region is expanded by adding one “–CH2–” group at a time from the left to the right ends of the molecule.

3.5 Hydrogenation of pentacene As a more challenging case for embedding methods, we consider a terminal hydrogenation of pentacene (Scheme 1). The equilibrium geometries of pentacene and its hydrogenated form were optimized with DFT-B3LYP. We employed the 6-31G*42 basis set for all of the calculations in this section. The reference reaction energy, computed at the total-system MP2 level, is 33.7 kcal/mol. Figure 7 presents the errors in the reaction energy obtained from three embedding methods; these errors are plotted as functions of the number of C atoms inside the embedded region, starting with the two terminal C atoms involved in the hydrogenation and then adding two C atoms at a time. It is clear that the reaction energy predicted by DMFET converges quickly to the correct value as we increase the size of the embedded region. In contrast, ONIOM performs poorly when the number of C atoms in the embedded region is in multiples of four. In these cases, the hydrogenated pentacene divides into a diradical (embedded region, in which all broken sigma bonds are capped by hydrogens) and a polycyclic aromatic hydrocarbon (environment with same capping of broken sigma bonds), as shown in Scheme 2(a). Such partitions break the 𝜋-conjugation in the total system, leading to localized unpaired electrons in the embedded region. ONIOM is unable to recover the correct electronic structure of the total system due to the lack of a proper embedding potential, which explains its failure. For the other cases where the embedded region contains 2, 6, 10, … carbon atoms, the double bonds in the hydrogenated pentacene are preserved in the two subsystems [Scheme 2(b)], indicating that the local electronic structure of the total system is not


22

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


perturbed significantly by the partition. Hence, ONIOM performs comparably to DMFET. Note that the partition of pentacene always leads to two conjugated systems (Scheme 3); thus, it is likely not a major source of error in the ONIOM calculations. Finally, sDMFET agrees with DMFET in most cases. In sDMFET, however, the virtual orbital space is truncated, and electron excitations from the occupied orbitals in the embedded region to the virtual orbitals in the environment are forbidden. This causes sDMFET to be less accurate than DMFET if the electrons are very delocalized. For example, when pentacene is divided at the middle ring (corresponding to 10 or 12 carbon atoms in the embedded region), the virtual orbitals in both the embedded region and the environment may contribute equally importantly to the correlation energies from the ECW calculations. Therefore, in such cases, sDMFET gives larger errors than DMFET, as can be seen in Fig. 7. Nevertheless, the overall performance of sDMFET is satisfactory.

Scheme 1. The hydrogenation reaction of pentacene. The areas enclosed by the blue and red boxes are examples of embedded regions that contain two and four C atoms, respectively.

(a) (b)


23


Page 24 of 41

Scheme 2. Illustration of dividing the hydrogenated pentacene into two subsystems: (a) a diradical (embedded region) and a tetracene, and (b) a hydrogenated naphthalene (embedded region) and a conjugated system. Both of the two subsystems are capped with two hydrogen atoms.

(a) (b) Scheme 3. Illustration of dividing pentacene into two subsystems: (a) a butadiene (embedded region) and a tetracene, and (b) a naphthalene (embedded region) and a conjugated system. Both of the two subsystems are capped with two hydrogen atoms.

Figure 7. Embedding calculations at the MP2/HF/6-31G* level for hydrogenation of pentacene. All errors are plotted as functions of the number of carbon atoms in the embedded region, starting with the two terminal carbon atoms that are hydrogenated (displayed inside the blue boxes in Scheme 1). The reference reaction energy was computed at the total-system MP2/6-31G* level to be 33.7 kcal/mol.


24

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


3.6 Defect-induced conical intersection in a Si nanocrystal Semiconductor nanocrystals can undergo nonradiative recombination via defect-induced conical intersections,58,59 which is an unwanted process in optoelectronics and phosphors as it lowers the energy conversion efficiency. As these intersections involve electronic excitations that are localized at specific defects in the nanocrystal,59 sDMFET may serve as an ideal technique to study such processes due to its efficiency and accuracy for treating localized excited states (Section 3.3). In this work, we apply sDMFET to investigate a potential conical intersection in a silicon nanocrystal with a surface silanol defect (Si44H45OH), which has been studied previously60 by the configuration interaction singles natural orbital complete active space configuration interaction (CISNOr-CASCI) method.61 The ground-state equilibrium and the S0/S1 minimum energy conical intersection (MECI) geometries were taken from Ref. [60], with the former optimized at the DFTCAM-B3LYP62/LANL2DZ63,64 level and the latter optimized at the CISNOr-CASCI/LANL2DZ level using a small active space of two electrons in three orbitals (2e, 3o). The embedded region in the sDMFET calculations was chosen to be the silanol group plus the neighboring functional groups to which the silanol group connects, as shown in Fig. 8. The appropriateness of such a partition can be assessed by comparing the highest-occupied MO (HOMO) and the lowestunoccupied MO (LUMO) computed on the total system with those from the embedded calculation, as the S1 state is associated with excitations from the HOMO to the LUMO. At the S0/S1 MECI geometry found earlier, we find that the Si nanocrystal has its HOMO and LUMO localized at the defect site [Figs. 9(a)-(b)], identical to the results from the embedded calculation [Fig. 9(c)], which validates our definition of the embedded region. It should be noted, however, that the frontier MOs of the Si nanocrystal are delocalized at the ground-state equilibrium structure;59 hence, it is


25


Page 26 of 41

inappropriate to study the excited states in the Franck-Condon region with the current sDMFET method. However, the ground-state energy difference between the two structures might still be evaluated accurately by sDMFET (as in pentacene in Section 3.5, where the electrons are delocalized over the entire system as well). In Table 2, we compare the energies of the S0 and S1 states computed at different levels of theory at the S0/S1 MECI structure. In the sDMFET calculations, the embedded region was treated by the EOM-CCSD method as well as its spin-flip variant65 (EOM-SF-CCSD), with the latter being known to perform better for open-shell multiplets.65 We also tested different mean-field theories to compute the embedding potentials, namely, HF and DFT-B3LYP, although both of these provided similar energetics. As can be seen from Table 2, sDMFET predicts the MECI to be ~4.6 eV above the ground-state minimum energy, which is slightly lower than the CISNOr-CASCI (2e, 3o) result at 4.88 eV. We therefore can confirm that this conical intersection, if it exists, is unlikely to cause nonradiative recombination due to its high energy. However, we also observe large energy gaps between the S0 and S1 states: ~0.7 eV by embedded EOM-CCSD and ~0.3 eV by embedded EOM-SF-CCSD. This, together with the fact that a very small active space was employed in the earlier CASCI calculations, indicates that the MECI studied here may not be a true conical intersection. As such, further research is required to locate the potential conical intersection at higher levels of theory.


26

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


Figure 8. (a) Ground-state equilibrium and (b) S0/S1 MECI structures of the Si nanocrystal adopted from Ref. [60]. The functional groups enclosed by the blue boxes are treated as the embedded region in the sDMFET calculations.

Figure 9. Comparison of the total-system and embedded calculations for the HOMO and LUMO of the Si nanocrystal at its S0/S1 MECI geometry. (a)-(b) HOMO and LUMO of the total system computed at the HF/LANL2DZ level. (c) HOMO and LUMO of the embedded region computed at the embedded HF/LANL2DZ level. (d) Density partition determined by sDMFET: the blue and red isosurfaces represent the electron densities of the embedded region and the environment, respectively. The isosurface values are set to 0.05 bohr ―3. The electron density that corresponds to the frozen orbitals is not shown.


27


Page 28 of 41

Table 2. Energies (in eV) of the S0 and S1 states for the Si nanocrystal at the previously predicted S0/S1 MECI geometry,60 relative to the ground-state minimum energy structure. The LANL2DZ basis set and effective core potential were employed in all calculations. Methods CISNOr-CASCI (2e, 3o)a EOM-CCSD/HFb EOM-CCSD/B3LYPb EOM-SF-CCSD/HFb EOM-SF-CCSD/B3LYPb a From Ref. [60]. b

S0 S1 S1-S0 gap 4.88 ~4.88 ~0.0 4.65 5.33 0.68 4.59 5.32 0.73 4.68 4.94 0.26 4.61 4.98 0.37

sDMFET results from this work.

3.7 Computational cost of sDMFET The main reason for developing sDMFET is to reduce the high computational cost of DMFET. In this section, we explore the efficiency of sDMFET through a numerical example. For coronene (C24H12), we use sDMFET and DMFET, respectively, to compute the embedding potentials between the C6 ring in the center and its surrounding environment. Figure 10 presents the wall times for calculating the embedding potential, along with the wall time for the total-system calculation for the ground-state energy. All of the calculations were performed at the HF/cc-pVDZ level. From the timing data, we see that sDMFET is two orders of magnitude faster than DMFET. This speedup is because the computational complexity of sDMFET is roughly 𝑂(𝐶 × 𝑛3), where n is the size of the subspace and C denotes the number of iterations required to converge the embedding potential, while the complexity of DMFET is 𝑂(𝐶 × 𝑁3) with N being the number of basis functions of the total system. Since n is usually much smaller than N, as the size of the embedded region is small (for example, n/N in this case is about 1/4), we expect sDMFET to be much more efficient than DMFET. Note that the wall times for the embedding calculations in Fig.


28

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


10 also include the time for a total-system energy calculation at the mean-field level. As such, the overall cost of sDMFET to determine the embedding potential is roughly the same as the cost of a supersystem mean-field calculation.

Figure 10. Computer time required to perform the total-system ground-state energy calculation, and sDMFET and DMFET optimizations of the embedding potentials on coronene at the HF/ccpVDZ level. The red and yellow bars also include the time to perform the total system HF calculation (black). The density partition obtained from the sDMFET calculation is also shown, with the red and blue isosurfaces represent the electron densities of the embedded region (C6 ring in the center) and the environment, respectively. The isosurface values are set to 0.05 bohr ―3.

4 Conclusions In this work, we improved the efficiency of the original DMFET by reformulating the original density matrix embedding problem within a truncated subspace determined by a Schmidt decomposition. The newly developed method, called sDMFET, is formally equivalent to DMFET for computing the unique embedding potential and density matrix partition between the


29


Page 30 of 41

subsystems. The only approximation introduced in sDMFET, as compared to DMFET, is that the subsequent ECW calculation does not allow electron transitions between the occupied orbital space of one subsystem and the virtual orbital space of the other subsystem. However, such transitions are usually less important when the phenomenon of interest take place in a small region of the total system. We carried out benchmark calculations including examples where systems are partitioned across hydrogen bonds, single covalent bonds, and even double bonds. In all cases, sDMFET agreed with DMFET very well in predicting accurate energetics. The example of hydrated uracil also proved the potential of sDMFET for describing localized electronic transition properties. The application of sDMFET to the study of a potential defect-induced conical intersection in a Si nanocrystal further demonstrated the method’s promise for investigating defects in condensed matter. Finally, periodic boundary conditions (PBCs) can be adopted for sDMFET in a straightforward way by use of the crystalline Gaussian-based atomic orbitals66 as the basis. This variant of sDMFET has been implemented by us and will be reported elsewhere. A more appealing way to incorporate PBCs is to represent KS orbitals derived from more accurate plane-wave basis DFT calculations by a smaller basis so that the size of the density matrix becomes tractable. These challenges form the basis of our ongoing work.

Appendix: Orthogonality between the occupied orbitals from two subsystems in DMFET If the total system is treated at the DFT level, the corresponding one-electron density matrix is idempotent: 𝐃𝐃 = 𝐃.

(A1)


30

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


At convergence, we have that the sum of subsystem density matrices is equal to the total density matrix: 𝐃A + 𝐃B = 𝐃.

(A2)

From the above two equations and the fact that subsystem density matrices are also idempotent, we therefore get 𝐃 A𝐃 B + 𝐃 B 𝐃 A = 𝟎 .

(A3)

Now, if we express the subsystem density matrices by the occupied MO coefficients, 𝐨A and 𝐨B, then Eq. (A3) becomes 𝐨A𝐨A† 𝐨B𝐨B† + 𝐨B𝐨B† 𝐨A𝐨A† = 𝟎.

(A4)

By multiplying 𝐨A† from the left and 𝐨B from the right of the equation above, we obtain 𝐨A† 𝐨B + 𝐨A† 𝐨B𝐨B† 𝐨A𝐨A† 𝐨B = 𝟎,

(A5)

𝐨A† 𝐨B(𝟏 + 𝐨B† 𝐃A𝐨B) = 𝟎.

(A6)

and

Since 𝐃A as a density matrix is positive semi-definite, we have ‖𝟏 + 𝐨B† 𝐃A𝐨B‖ ≥ 1. Therefore, 𝐨A† 𝐨B = 𝟎, which means that the occupied orbitals from the two subsystems are mutually orthogonal.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] ORCID Xing Zhang: 0000-0002-1892-1380


31


Page 32 of 41

Emily A. Carter: 0000-0001-7330-7554 Notes The authors declare no competing financial interest. ACKNOWLEDGMENT This material is based upon work supported by the National Science Foundation (Grant No. CHE-1265700), and the U.S. Department of Energy, Office of Science, Offices of Basic Energy Sciences and Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program. We thank Ms. Nari Baughman for careful editing of this manuscript.

REFERENCES (1)

Huang, P.; Carter, E. A. Advances in Correlated Electronic Structure Methods for Solids,

Surfaces, and Nanostructures. Annu. Rev. Phys. Chem. 2008, 59, 261–290. (2)

Libisch, F.; Huang, C.; Carter, E. A. Embedded Correlated Wavefunction Schemes: Theory

and Applications. Acc. Chem. Res. 2014, 47, 2768–2775. (3)

Jacob, C. R.; Neugebauer, J. Subsystem Density-Functional Theory. Wiley Interdiscip.

Rev.: Comput. Mol. Sci. 2014, 4, 325–362. (4)

Cortona, P. Self-Consistently Determined Properties of Solids without Band-Structure

Calculations. Phys. Rev. B 1991, 44, 8454–8458.


32

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


(5)

Wesolowski, T. A.; Warshel, A. Frozen Density Functional Approach for Ab Initio

Calculations of Solvated Molecules. J. Phys. Chem. 1993, 97, 8050–8053. (6)

Wesolowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for

Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891–5928. (7)

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. (8)

Govind, N.; Wang, Y. A.; Carter, E. A. Electronic-Structure Calculations by First-

Principles Density-Based Embedding of Explicitly Correlated Systems. J. Chem. Phys. 1999, 110, 7677–7688. (9)

Klüner, T.; Govind, N.; Wang, Y. A.; Carter, E. A. Periodic Density Functional Embedding

Theory for Complete Active Space Self-Consistent Field and Configuration Interaction Calculations: Ground and Excited States. J. Chem. Phys. 2002, 116, 42–54. (10) 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. (11) Jacob, C. R.; Visscher, L. A Subsystem Density-Functional Theory Approach for the Quantum Chemical Treatment of Proteins. J. Chem. Phys. 2008, 128, 155102. (12) 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.


33


Page 34 of 41

(13) Roncero, O.; de Lara-Castells, M. P.; Villarreal, P.; Flores, F.; Ortega, J.; Paniagua, M.; Aguado, A. An Inversion Technique for the Calculation of Embedding Potentials. J. Chem. Phys. 2008, 129, 184104. (14) Roncero, O.; Zanchet, A.; Villarreal, P.; Aguado, A. A Density-Division Embedding Potential Inversion Technique. J. Chem. Phys. 2009, 131, 234110. (15) Goodpaster, J. D.; Ananth, N.; Manby, F. R.; Miller, T. F. Exact Nonadditive Kinetic Potentials for Embedded Density Functional Theory. J. Chem. Phys. 2010, 133, 084103. (16) 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. (17) Goodpaster, J. D.; Barnes, T. A.; Miller, T. F. Embedded Density Functional Theory for Covalently Bonded and Strongly Interacting Subsystems. J. Chem. Phys. 2011, 134, 164108. (18) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F. A Simple, Exact DensityFunctional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564–2568. (19) 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. (20) Libisch, F.; Marsman, M.; Burgdörfer, J.; Kresse, G. Embedding for Bulk Systems Using Localized Atomic Orbitals. J. Chem. Phys. 2017, 147, 034110. (21) 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.


34

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


(22) Chulhai, D. V.; Goodpaster, J. D. Projection-Based Correlated Wave Function in Density Functional Theory Embedding for Periodic Systems. J. Chem. Theory Comput. 2018, 14, 1928– 1942. (23) Knizia, G.; Chan, G. K. Density Matrix Embedding: A Strong-Coupling Quantum Embedding Theory. J. Chem. Theory Comput. 2013, 9, 1428–1432. (24) Wouters, S.; Jiménez-Hoyos, C. A.; Sun, Q.; Chan, G. K. L. A Practical Guide to Density Matrix Embedding Theory in Quantum Chemistry. J. Chem. Theory Comput. 2016, 12, 2706– 2719. (25) Yu, K.; Krauter, C. M.; Dieterich, J. M.; Carter, E. A. Density and Potential Functional Embedding: Theory and Practice. In Fragmentation: Toward Accurate Calculations on Complex Molecular Systems; Gordon, M. S., Ed.; John Wiley & Sons, Ltd: Chichester, West-Sussex, United Kingdom, 2017; pp 81–117. (26) Zhang, X.; Carter, E. A. Kohn-Sham Potentials from Electron Densities Using a Matrix Representation within Finite Atomic Orbital Basis Sets. J. Chem. Phys. 2018, 148, 034105. (27) Huang, C.; Pavone, M.; Carter, E. A. Quantum Mechanical Embedding Theory Based on a Unique Embedding Potential. J. Chem. Phys. 2011, 134, 154110. (28) Huang, C.; Carter, E. A. Potential-Functional Embedding Theory for Molecules and Materials. J. Chem. Phys. 2011, 135, 194104. (29) 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.


35


Page 36 of 41

(30) 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. (31) Ou, Q.; Carter, E. A. Potential Functional Embedding Theory with an Improved Kohn– Sham Inversion Algorithm. J. Chem. Theory Comput. 2018, 14, 5680–5689. (32) Yu, K.; Carter, E. A. Extending Density Functional Embedding Theory for Covalently Bonded Systems. Proc. Natl. Acad. Sci. 2017, 114, E10861–E10870. (33) Peschel, I. Special Review: Entanglement in Solvable Many-Particle Models. Braz. J. Phys. 2012, 42, 267–291. (34) 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. (35) Schmidt, E. Zur Theorie Der Linearen Und Nichtlinearen Integralgleichungen. Math. Ann. 1907, 63, 433–476. (36) Foster, J. M.; Boys, S. F. Canonical Configurational Interaction Procedure. Rev. Mod. Phys. 1960, 32, 300–302. (37) Pipek, J.; Mezey, P. G. A Fast Intrinsic Localization Procedure Applicable for Ab Initio and Semiempirical Linear Combination of Atomic Orbital Wave Functions. J. Chem. Phys. 1989, 90, 4916–4926. (38) Klich, I. Lower Entropy Bounds and Particle Number Fluctuations in a Fermi Sea. J. Phys. A: Math. Gen. 2006, 39, L85–L91.


36

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


(39) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K.-L. PySCF: The Python-Based Simulations of Chemistry Framework. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2018, 8, e1340. (40) Huang, C. Extending the Density Functional Embedding Theory to Finite Temperature and an Efficient Iterative Method for Solving for Embedding Potentials. J. Chem. Phys. 2016, 144, 124106. (41) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H. B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678–5796. (42) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213–222. (43) 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. (44) 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. (45) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (46) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789.


37


Page 38 of 41

(47) Dahlke, E. E.; Olson, R. M.; Leverentz, H. R.; Truhlar, D. G. Assessment of the Accuracy of Density Functionals for Prediction of Relative Energies and Geometries of Low-Lying Isomers of Water Hexamers. J. Phys. Chem. A 2008, 112, 3976–3984. (48) Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-Body Expansion for Large Systems, with Applications to Water Clusters. J. Chem. Theory Comput. 2007, 3, 46–53. (49) Taylor, C. R.; Bygrave, P. J.; Hart, J. N.; Allan, N. L.; Manby, F. R. Improving Density Functional Theory for Crystal Polymorph Energetics. Phys. Chem. Chem. Phys. 2012, 14, 7739– 7743. (50) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. III. The Atoms Aluminum through Argon. J. Chem. Phys. 1993, 98, 1358–1371. (51) Liptay, W. Electrochromism and Solvatochromism. Angew. Chem., Int. Ed. Engl. 1969, 8, 177–188. (52) Mewes, J. M.; You, Z. Q.; Wormit, M.; Kriesche, T.; Herbert, J. M.; Dreuw, A. Experimental Benchmark Data and Systematic Evaluation of Two a Posteriori, PolarizableContinuum Corrections for Vertical Excitation Energies in Solution. J. Phys. Chem. A 2015, 119, 5446–5464. (53) Zhang, X.; Herbert, J. M. Excited-State Deactivation Pathways in Uracil versus Hydrated Uracil: Solvatochromatic Shift in the 1nπ ∗ State Is the Key. J. Phys. Chem. B 2014, 118, 7806– 7817.


38

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


(54) Shao, Y.; Head-Gordon, M.; Krylov, A. I. The Spin–Flip Approach within TimeDependent Density Functional Theory: Theory and Applications to Diradicals. J. Chem. Phys. 2003, 118, 4807–4818. (55) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098–3100. (56) Hehre, W. J.; Ditchfield, K.; 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. (57) Stanton, J. F.; Bartlett, R. J. The Equation of Motion Coupled-cluster Method. A Systematic Biorthogonal Approach to Molecular Excitation Energies, Transition Probabilities, and Excited State Properties. J. Chem. Phys. 1993, 98, 7029–7039. (58) Shu, Y.; Fales, B. S.; Levine, B. G. Defect-Induced Conical Intersections Promote Nonradiative Recombination. Nano Lett. 2015, 15, 6247–6253. (59) Shu, Y.; Fales, B. S.; Peng, W.-T.; Levine, B. G. Understanding Nonradiative Recombination through Defect-Induced Conical Intersections. J. Phys. Chem. Lett. 2017, 8, 4091– 4099. (60) Shu, Y.; Levine, B. G. First-Principles Study of Nonradiative Recombination in Silicon Nanocrystals: The Role of Surface Silanol. J. Phys. Chem. C 2016, 120, 23246–23253. (61) Shu, Y.; Hohenstein, E. G.; Levine, B. G. Configuration Interaction Singles Natural Orbitals: An Orbital Basis for an Efficient and Size Intensive Multireference Description of Electronic Excited States. J. Chem. Phys. 2015, 142, 024102.


39


Page 40 of 41

(62) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange–Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. (63) Dunning, T. H.; Hay, P. J. Gaussian Basis Sets for Molecular Calculations. In Methods of Electronic Structure Theory. Modern Theoretical Chemistry, vol 3; Schaefer, H. F., Ed.; Springer US: Boston, MA, 1977; pp 1–27. (64) Wadt, W. R.; Hay, P. J. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for Main Group Elements Na to Bi. J. Chem. Phys. 1985, 82, 284–298. (65) Krylov, A. I. Size-Consistent Wave Functions for Bond-Breaking: The Equation-ofMotion Spin-Flip Model. Chem. Phys. Lett. 2001, 338, 375–384. (66) McClain, J.; Sun, Q.; Chan, G. K.-L.; Berkelbach, T. C. Gaussian-Based Coupled-Cluster Theory for the Ground-State and Band Structure of Solids. J. Chem. Theory Comput. 2017, 13, 1209–1218.


40

Page 41 of 41

6000

10 5

5000

0 -5

Wall time (s)

∆E-∆EMP2 (kcal/mol)

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


-10 -15

H2

-20 -25

sDMFET DMFET ONIOM

-30 -35 -40

2

4

6

8

10

12

14

16

18

20

No. of C atoms in the embedded region

4000 3000 2000 100

22

0

sDMFET

DMFET


41

Subspace density matrix functional embedding theory: Theory

Recommend Documents