Dual basis set approach for density functional and wave function

Subscriber access provided by University of Sussex Library

Quantum Electronic Structure

Dual basis set approach for density functional and wave function embedding schemes Bence Hégely, Péter R. Nagy, and Mihaly Kallay J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00350 • Publication Date (Web): 26 Jul 2018 Downloaded from http://pubs.acs.org on July 27, 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 54 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

Dual basis set approach for density functional and wave function embedding schemes Bence Hégely, Péter R. Nagy, and Mihály Kállay∗ MTA-BME Lendület Quantum Chemistry Research Group, Department of Physical Chemistry and Materials Science, Budapest University of Technology and Economics, H-1521 Budapest, P.O.Box 91, Hungary E-mail: [email protected]

Abstract A dual-basis (DB) approach is proposed which is suitable for the reduction of the computational expenses of the Hartree–Fock, Kohn–Sham, and wave-function-based correlation methods. The approach is closely related to the DB approximation of HeadGordon and co-workers [J. Chem. Phys. 125, 074108 (2006)] but specifically designed for embedding calculations. The new approach is applied to our variant of the projectorbased embedding theory utilizing the Huzinaga-equation, multi-level local correlation methods, and combined density functional-multi-level local correlation approximations. The performance of the resulting DB density functional and wave function embedding methods are evaluated in extensive benchmark calculations and also compared to that of the corresponding embedding schemes exploiting the mixed-basis approximation. Our results show that, with an appropriate combination of basis sets, the DB approach significantly speeds up the embedding calculations, and, for chemical processes where the electronic structure considerably changes, it is clearly superior to the mixed-basis ∗

To whom correspondence should be addressed

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

approximation. The results also demonstrate that the DB approach, if integrated with the mixed-basis approximation, efficiently eliminates the major weakness of the latter, and the combination of the DB and mixed-basis schemes is the most efficient strategy to accelerate embedding calculations.

1

Introduction

The computational expenses of electronic structure calculations steeply grow with the size of the studied systems. 1 The computation time of the Hartree–Fock (HF) and Kohn–Sham (KS) density functional theory (DFT) calculations using a hybrid functional scale formally as the fourth-power of the systems size (N ), and the scaling of KS calculations with a “pure”, i.e., non-hybrid functional is still N 3 . The more advanced ab initio correlation methods, such as the second-order Møller–Plesset (MP2) approach 2 or the coupled-cluster (CC) models, 3 e.g., CC with single and double excitations (CCSD) 4 and CCSD with perturbative triples [CCSD(T)], 5 exhibit a scaling of at least N 5 . Consequently, numerous approaches have been developed over the years to break down the expenses of quantum chemical calculations. Since the computation time for any electronic structure calculation strongly depends on the size of the atomic orbital (AO) basis set applied, the most plausible strategy to decrease the computational costs is the reduction of the size of the AO basis using, for instance, the mixed-basis set approach. In this simplistic approximation the atoms of the chemically relevant part of the system are described using a large AO basis set, whereas a smaller basis set is employed for the remaining atoms. 6–14 Although this approach is conceptually simple and trivial to implement, it is not elegant from the theoretical point of view. The use of basis sets of different size gives rise to the artificial polarization of the electronic structure since the electron density grows in the region where more basis functions can be found. These artifacts considerably worsen the efficiency of the approach for chemical processes where the electronic structure significantly changes, such as charge-manifestation processes. 15 In addition, the use of mixed basis sets might cause convergence problems for self-consistent 2


Page 2 of 54

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


field (SCF) algorithms and also increases the probability that the SCF procedure converges to an unstable solution (see Sect. 3.4). The so-called dual basis (DB) set approach, which is free from the aforementioned problems, is another related strategy to decrease the basis set size. In a DB HF or KS SCF calculation the SCF equations are first solved in a relatively small basis set, the resulting electron density or molecular orbitals (MOs) are some way expanded in a bigger AO basis, and usually an approximate SCF energy is also evaluated with the latter basis set, hereafter referred to as the large basis. For the computation of the correlation energy in a DB framework the simplest approach is to use the larger basis set together with the MOs constructed in the preceding DB HF calculation, but more complicated correlated DB methods also exist, in which two basis sets are employed also at the correlation calculation. The DB approximation was first employed in 1983 at the HF level by Havriliak and King, 16 who studied the Rydberg states of the ammonium radical. Since HF calculations with large basis sets including diffuse functions were not feasible at that time, the authors solved the SCF equations with a smaller basis set. Then an appropriate virtual MO space was constructed by orthonormalizing a bigger AO basis containing the required diffuse functions to the occupied space and diagonalizing the Fock operator in the new basis. For the first time, DB sets were applied in correlated calculations by Jurgens-Lutovsky and Almlöf. 17 In their DB MP2 method, a DB HF calculation was performed following Havriliak and King, and the MP2 energy was evaluated as usual but utilizing that occupied orbitals were expanded in the small basis, and consequently only integrals with two indices in the larger basis are required. This approach was later also adapted in the DB MP2 method of Wolinski and Pulay. 18,19 Significant progress has been made on DB methods by Head-Gordon and co-workers. 20–29 In their DB SCF approach an SCF calculation is converged in a small basis set, the density matrix is projected into the large basis set, a single Roothaan step is taken, and an energy correction is evaluated. 22,23 Alternatively, the expression for the energy lowering can be reformulated to avoid the diagonalization of the Fock matrix, which makes the method

3



suitable for a linear-scaling implementation. 22 Based on the new DB HF approach a DB MP2 method was also developed utilizing the density-fitting (DF) MP2 technique. 23 Later analytic gradients were also implemented for the DB SCF 24 and the DF-MP2 26 methods. The utility of these approaches was demonstrated in ab initio molecular dynamics simulations 28 and for non-covalent interactions. 27,30 Further promising DB HF approximations were introduced by Gill et al. deriving perturbative corrections that estimate the complete basis set HF energy by diagonalizing the Fock potential from a small basis in a large one. 31,32 New DB MP2 methods were also proposed by these authors, where the MOs are obtained from the large basis HF calculation, and the DB approximation is invoked at the integral transformation. 33,34 An interesting dual-level DFT approach was published by Nakajima and Hirao, where, in the primary SCF calculation, not only a smaller basis set but also a less sophisticated functional was used. 35 This approach was also extended to time-dependent DFT 36 and relativistic DFT, 37 while Gill and co-workers generalized the method by also applying a cruder quadrature grid at the first SCF run. 38 The DB approximation was used in reduced-scaling correlation methods by Kobayashi and Nakai, 39 Friedrich and Dolg, 40–43 and very recently Røeggen and Gao proposed a scheme closely related to the DB approach for the calculation of the correlation energy of molecules and solids. 44 We also note that the complete auxiliary basis set singles correction employed frequently for improving the HF reference energy in explicitly correlated methods 45–48 can also be regarded as a special DB approximation. The other large class of reduced-cost methods that are relevant for our study includes quantum embedding approaches. 49–57 In these methods the system is divided into an active subsystem comprising the chemically important part and its environment. The latter is described at a lower-level of theory, such as semi-empirical or DFT, while a more accurate approximation, such as a more advanced DFT approach or a wave function theory (WFT) method is used for the active subsystem. The arising DFT-in-DFT and WFT-inDFT schemes face the important problem of the coupling of the different approaches. In an

4


Page 4 of 54

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


orbital free representation of the electron density or if the orbitals representing the embedded and embedding electron density are not orthogonal, the sum of the subsystem electron densities does not agree with the total density. Consequently, the kinetic energy is also not additive, which can be remedied by the use of approximate nonadditive functionals 58,59 or with optimized effective potential methods. 60–66 In the alternative projector-based embedding theories developed by Rajchel et al., 67 Manby and co-workers, 68–71 and Hoffmann et al. 72 projectors are employed to impose the orthonormality of the active and the environmental orbitals. Another new approach, which resolves the above issue in an elegant way, is the embedded mean-field theory (EMFT) developed by Miller and Manby and their co-workers. 73,74 In EMFT the system is partitioned at the basis set level, the two subsystems are treated by a low- and a high-level DFT approach, and the subsystem densities are simultaneously iterated in a single SCF run. Unfortunately, EMFT does not allow for WFTin-DFT-type embedding calculations. The embedding of a WFT method into a lower-level WFT approach does not seem to be a trivial task either. Nevertheless, simple WFT-in-WFT type embedding schemes can be developed if the locality of the electron correlation is exploited, and the electrons localized on different parts of a molecule are treated at different levels of theory. 75–81 Recently, we put forward a simple extension of the projector-based embedding of Manby and co-workers, 80 which, instead of a somewhat arbitrary level-shift projector, uses the Huzinaga-equation. 82 In the Huzinaga-equation-based embedding, hereafter shortly referred to as Huzinaga-embedding, the orthonormality of the active and inactive orbitals is strictly enforced. The potential of the Huzinaga-equation was subsequently also realized by Chulhai and Goodpaster for embedding schemes where the subsystem’s density or wave function is expanded in terms of basis functions associated with the subsystem’s atoms. 83 A simple alternative to the Huzinaga-equation-based scheme was proposed by Hammes-Schiffer et al. utilizing the orthogonality constrained basis set expansion procedure, which also allows both the active and the environmental densities to respond to each other. 84

5



The mixed-basis approach has also been combined with quantum embedding schemes, where the use of mixed basis sets seems to be a natural way to reduce the computational costs because the system is anyway partitioned into high- and low-level regions. 73,74 Of course, the major disadvantage of the mixed-basis methods, the imbalanced description of the high- and low-level densities is also inherited by the combined approaches and can be even amplified for embedding theories where the environmental densities are not frozen, such as EMFT, 15 especially when using a minimal basis set for the environment. For the latter method it was recently demonstrated that the electron density artifacts can be largely eliminated by the semi-empirical Fock-corrected DFT approach, 85 where the energy evaluated with an inexpensive functional and a minimal basis is corrected to recover the KS energy obtained with an accurate functional and a larger basis set. Nevertheless, to our knowledge, embedding schemes have not yet been combined with the genuine DB approaches presented above. In this paper we propose a slightly modified version of the DB method of Head-Gordon and co-workers. 23 The new approach might be interesting on its own right but is particularly suitable for combination with embedding theories. We apply the new approach to our Huzinaga-embedding method. The performance of the resulting DB embedding theories are evaluated in extensive benchmark calculations and also compared to that of the corresponding mixed-basis embedding schemes.

2 2.1

Theory Dual basis set approach

Let us assume that we have a system of n electrons and two basis sets, a small and a large one; the former is not necessarily a subset of the latter. We are interested in the HF or KS energy evaluated with the large basis set. The corresponding energy functional will simply be denoted by E[D], where D is a one-particle density matrix. The Fock-matrix required in

6


Page 6 of 54

Page 7 of 54 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 SCF procedure is defined as F=

∂E . ∂D

(1)

If not obvious, notation F[D] will be used for the Fock-matrix computed from density matrix D. The DB calculation starts with solving the SCF equations in the small basis set. The resulting occupied MOs are projected into the large basis and orthonormalized. Since the method is invariant with respect to the orthogonalization algorithm, any one can be chosen. From the orthonormal MOs a density matrix, D0 , is computed. In the case the initial SCF is converged with a density-based algorithm, i.e., no MOs are available, we can directly project the density matrix onto the large basis and use some purification technique to construct an n-representable D0 . If MOs are needed, they can be obtained as the natural orbitals of D0 transformed to the orthogonalized AO basis, that is, by diagonalizing matrix S1/2 D0 S1/2 , where S stands for the AO overlap matrix of the large basis. Subsequently, one SCF iteration is performed with D0 in the large basis set. The Fockmatrix F[D0 ] is constructed, the E[D0 ] energy is calculated, the Fock-matrix is diagonalized, and an improved density matrix, D1 , is computed from the new MOs. To obtain a good approximation to the converged SCF energy with the large basis set, a first-order estimate for E[D1 ], E [1] [D1 ], is derived. If we expand E[D] into Taylor series around D0 , drop the quadratic and higher-order terms, and utilize Eq. (1), we arrive at the

E [1] [D1 ] = E[D0 ] + Tr [F[D0 ] (D1 − D0 )]

(2)

approximation, which also defines our DB SCF energy. If a DB correlation calculation is carried out, the DB SCF energy is evaluated as described above, and the correlation energy is calculated with Fock-matrix F[D0 ] and the MOs obtained by the diagonalization of F[D0 ]. Obviously, both the DB SCF and the DB correlation approaches are exact in the sense that, if the small basis and large basis sets are identical, the SCF or the correlation energy

7



computed with the given basis set is retrieved. The DB approach outlined here differs at two important points from that presented in Ref. 23. First, in our energy expression, Eq. (2), the zeroth-order term, E[D0 ], is an approximate SCF energy evaluated in the large basis, while Head-Gordon and co-workers use the converged small basis SCF energy instead. In this respect our method is analogous to the approach of Gill et al. 31 Second, in contrast to Ref. 23 and most of the other DB schemes our density matrix D0 is n-representable. On the one hand, these changes makes the theory somewhat more accurate and theoretically more satisfactory. On the other hand, the calculation of energy gradients is more complicated with our approach, and our purified D0 is less sparse than the density matrix obtained simply by projecting the small basis density matrix into the large basis, which worsens the efficiency of the Fock-matrix construction in the large basis. Anyway, the n-representability of the density matrix or, equivalently, the orthonormality of the orbitals used to evaluate the large-basis Fock-matrix is a must for embedding theories, which are of primary concern in this study. Furthermore, as it will be obvious soon, the generalization of Eq. (2) to our embedding method would be problematic if E[D0 ] was replaced by the converged small basis SCF energy.

2.2

Huzinaga-equation-based embedding

Here we briefly review our Huzinaga-embedding, 80 which is a modified version of the projectorbased embedding theory of Manby and co-workers. 68,69 As customary in embedding theories the system is divided into two subsystems defined by the corresponding atoms, an embedded (active) part A and its environment B. In the case of DFT-in-DFT embedding the two subsystems are described by a high- and a low-level DFT method, defined by energy functionals E1 [D] and E2 [D], respectively. The corresponding Fock-matrices will be denoted by F1 [D] and F2 [D], respectively. In WFT-in-DFT calculations, a WFT method is used for the active subsystem. Then a HF-in-DFT calculation is carried out first, and E1 [D] will stand for the HF energy. In the initial step of the embedding 8


Page 8 of 54

Page 9 of 54 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


calculation the SCF equations of the low-level method defined by E2 [D] are converged for the entire system. The occupied MOs are localized, and for each localized MO (LMO) the atoms are identified on which the orbital is localized invoking the Boughton–Pulay (BP) algorithm 86 or an algorithm based on Mulliken populations (see the supplementary material for more details). The LMOs which are localized on the atoms of the embedded subsysB tem and the environment will be denoted, respectively, by φA i and φi . The selection of

the orbitals also determines the number of the electrons in the embedded and embedding B subsystems, for which symbols nA and nB will stand, respectively. From orbitals φA i and φi

the density matrices of the two subsystems, DA and DB , respectively, can be computed, and the density matrix of the entire system, DAB , is obviously their sum. In the second step of the embedding calculation a further SCF calculation is performed for the active subsystem of nA electrons keeping the φB i orbitals frozen. The DFT-in-DFT (HF-in-DFT) energy is evaluated as

DFT ˜ A E12 [D ; DA , DB ]

˜ A ; DA , DB ] + Tr = E12 [D

A A ∂E12 ˜ D −D , ∂DA

(3)

where ˜ A ; DA , DB ] = E2 [DAB ] + E1 [D ˜ A ] − E2 [DA ], E12 [D

(4)

˜ A is the self-consistent density matrix of the second SCF run computed from orbitals and D φÃ i . Note that E12 is the SCF energy of the entire system calculated with the low-level method replacing the contribution of the density of the active subsystem, DA , by that of the ˜ A , and the second term on the right-hand-side of Eq. (3) higher-level self-consistent density D ˜ A; D ˜ A , DB ] and E12 [D ˜ A ; DA , DB ]. is a first-order correction to the difference between E12 [D B To ensure the orthogonality of orbitals φÃ i to the frozen φi orbitals, the

A B ˜ Ã Ã FDFT − SDB FDFT − FDFT 12 12 12 D S C = SC E

9


(5)


Page 10 of 54

Ã A B Huzinaga-equation is solved, where Fock-matrix FDFT 12 [D ; D , D ] is defined as FDFT = 12

∂E12 ∂E12 + , ˜ A ∂DA ∂D

(6)

˜ A denotes the matrix of MO coefficients for orbitals φÃ , and E ˜ A is a diagonal matrix with C i the corresponding orbital energies on its diagonal. Here and in the following, for brevity, we omit the arguments of E12 and FDFT 12 . By default, these will be deemed to be constructed ˜ A , DA , and DB , and it will be pointed out if the corresponding approximate of matrices D density matrices are used instead of the latter. It is easy to see that the above embedding is exact, which means that the energy of the entire system evaluated with a mean-field method is recovered if that particular method is embedded into itself. In the case of WFT-in-DFT embedding, after the HF-in-DFT calculation, a WFT calculation is performed with the above Fock-matrix and frozen φB i orbitals. The WFT-in-DFT energy is evaluated as

WFT ˜ A , DA , DB ] E12 [ΨA ; D

AB

= E2 [D

˜ A , DA , DB ]−E2 [DA ]−Tr ]+E1WFT [ΨA ; D

∂E12 DA ∂DA

, (7)

where E1WFT is the total energy computed with the WFT method using Fock-matrix FDFT 12 , and ΨA is the corresponding wave function.

2.3

Huzinaga-embedding with dual basis set

We seek the generalization of the DB approach presented in Sect. 2.1 for the Huzinagaembedding outlined in Sect. 2.2, which is not unambiguous. Of the several possibilities, we explored two ones. In Ansatz 1 (A1), which is probably the most straightforward DB variant of the HuzinagaDFT embedding, E12 is computed with the small basis set as described above. The occupied B B φA i and φi orbitals are projected into the large basis set. The projected φi orbitals are

orthonormalized, and the resulting orbitals, φB 0i , are used to compute the nB -representable 10


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


B B density matrix DB 0 . We suppose that the φ0i orbitals and the D0 density matrix are good

enough approximations for the corresponding large basis quantities, or at least their quality is less important compared to that of the orbitals and densities of subsystem A. Hence, we do not endeavor to improve them, and they will be kept fixed in the following. The B projected φA i orbitals are also orthonormalized, but beforehand orbitals φ0i are projected out

by multiplying the corresponding MO coefficient matrix with the projector Q = 1 − DB 0 S to ensure the orthogonality of the MOs of the active subsystem and the environment. From the A φA 0i orbitals constructed in this way, the nA -representable density matrix D0 is computed, and B Ã = DA we also calculate the n-representable DAB 0 + D0 density matrix. The φi orbitals are 0 B processed similarly to φA i : they are projected onto the large basis, orbitals φ0i are projected

out of them using Q, the resulting MOs are orthonormalized to construct orbitals φÃ 0i , which ˜ A . Notice again that φÃ and φB are orthogonal are used to evaluate the nA -representable D 0i 0i 0 to each other. Utilizing the new density matrices estimates for the large-basis values of the three terms on the right-hand-side of Eq. (4) are calculated relying on our DB approach. First, the Fockand diagonalized. A matrix corresponding to the low-level method is constructed with DAB 0 AB ] new approximation to DAB , DAB 1 , is computed from the MOs, and an estimate for E2 [D

is evaluated as [1] AB AB AB AB E2 [DAB . 1 ] = E2 [D0 ] + Tr F2 [D0 ] D1 − D0

(8)

Second, to calculate an approximate E2 [DA ] in a similar way, an appropriate DA 1 is required. For that purpose one possibility would be to localize the orbitals obtained in the first step and assign the LMOs to the subsystems A and B as described in Sect. 2.2. In this case, however, B orbitals φB 0i and density matrix D0 would be redefined, and we should avoid this to preserve

the consistency with the other contribution. Instead, we project out the φB 0i orbitals from the occupied MOs calculated in the first step using projector Q. This results in nA pieces of orbitals similar to φA 0i and nB MOs of close to zero norm. These orbitals are orthonormalized employing canonical orthogonalization dropping nB eigenvectors of the overlap matrix with 11



Page 12 of 54

small eigenvalues. The remaining orbitals are the updated MOs of the active subsystem, φA 1i , A which are used to compute DA 1 . The F2 [D0 ] Fock-matrix is constructed, and the approximate

E2 [DA ] contribution is evaluated as [1] A A A A E2 [DA . 1 ] = E2 [D0 ] + Tr F2 [D0 ] D1 − D0

(9)

A Note that the diagonalization of F2 [DA 0 ] is avoided here as D1 is calculated in an indi-

˜ A ] Fock-matrix is computed, which is then used to construct rect way. Third, the F1 [D 0 DFT DFT B Ã A B FDFT − SDB − FDFT 12 [D0 ; D0 , D0 ] and the F12 0 F12 12 D0 S Huzinaga-matrix. The latter is B diagonalized to get orbitals φÃ 1i , which are still orthogonal to orbitals φ0i . The updated den-

˜ A , is evaluated using orbitals φÃ . The approximate sity matrix of the active subsystem, D 1i 1 ˜ A ] contribution is calculated similarly to the other ones as large basis E1 [D h i [1] ˜ A A A A A ˜ ˜ ˜ ˜ ] = E [ D ] + Tr F [ D ] D − D . E 1 [D 1 1 1 0 0 1 0

(10)

After introducing the above three contributions the large basis estimate for E12 can be defined by the [1] ˜ A [1] [1] ˜ A [1] A B AB A E12 [D 1 ; D1 , D0 ] = E2 [D1 ] + E1 [D1 ] − E2 [D1 ]

(11)

expression. To derive our large basis embedding energy we should consider the last term of Eq. (3). This is a first-order correction for the relaxation of the density of the active subsystem, thus we do not attempt to add another first-order correction for this term for the basis˜ A , DA , and E12 [D ˜ A ; DA , DB ] set extension, but it is simply evaluated with the available D 1 1 0 0 0 quantities. Consequently, our final DB DFT-in-DFT (HF-in-DFT) energy of A1 reads as DFT[1] ˜ A B E12 [D1 ; DA 1 , D0 ]

=

[1] ˜ A A B E12 [D 1 ; D1 , D0 ]

+ Tr

Ã D 1

−

DA 1

∂E 12 ∂DA 0

,

(12)

In the case of a WFT-in-DFT approach the WFT calculation is carried out with the B Ã Ã A B FDFT 12 [D0 ; D0 , D0 ] Fock-matrix in the basis of the φ1i MOs keeping orbitals φ0i frozen, and

12


Page 13 of 54 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 A1 DB WFT-in-DFT energy is evaluated according to the WFT[1] ˜ A , DA , DB ] E12 [ΨA ; D 1 1 0

=

[1] E2 [DAB 1 ]

+

˜ A , DA , DB ] E1WFT [ΨA ; D 1 1 0

−

[1] E2 [DA 1]

− Tr

∂E12 DA 1 ∂DA 0

(13) equation. It is easy to see that Eq. (12) reduces to Eq. (2) in the limit that the entire system is active, and A1 is exact in the sense that the DB DFT energy evaluated for the entire system with the high-level DFT method is recovered in this limit. A1 also inherits the exactness of the DB approach, that is, independently of the partitioning of the system, the DFT-in-DFT energy, Eq. (3), calculated with the particular basis set is reproduced if the small basis and large basis sets are the same. These also imply that similar statements hold for the DB WFT-in-DFT energy obtained with A1. We have also designed a somewhat simpler DB embedding scheme, Ansatz 2 (A2). Here the SCF equations are solved with the small basis set and the low-level method for the entire system, the orbitals are localized and assigned to the subsystems, but the further steps of the embedding calculation with the small basis are not taken. With the large basis AB B A AB B A A B density set the φA 0i , φ0i , and φ1i orbitals as well as the D0 , D0 , D0 , D1 , D1 , and D1 [1]

[1]

A matrices are constructed and the low-level energy contributions, E2 [DAB 1 ] and E2 [D1 ],

are calculated as in the case of A1. In the third step, however, the Huzinaga-equations DFT DFT Ã A B are iterated to convergence with the FDFT − SDB − 12 [D ; D0 , D0 ] Fock- and the F12 0 F12 B ˜ A is the self-consistent density matrix formed from FDFT 12 D0 S Huzinaga-matrix, where D

orbitals φÃ i . The DFT-in-DFT energy is calculated according to Eqs. (11) and (12) with ˜ A by D ˜ A and E [1] [D ˜ A ] by E1 [D ˜ A ]. If a WFT-in-DFT calculation is carried out, replacing D 1 1 1 B Ã A B the correlation energy is computed using the Fock-matrix FDFT 12 [D ; D0 , D0 ] with the φ0i WFT[1] ˜ A , DA , DB ]. and φÃ [ΨA ; D i orbitals, and the energy is evaluated as E12 1 0

One can simply prove that the A2 DFT-in-DFT (WFT-in-DFT) approach, in contrast to A1, is exact in the sense that the high-level DFT (WFT) energy of the whole system is

13



retrieved if subsystem A incorporates the entire system. In turn, it is also true for A2 that the DFT-in-DFT (WFT-in-DFT) energy of the original Huzinaga-embedding is reproduced with identical small and large basis sets. Concerning the computational expenses of our DB DFT-in-DFT embedding scheme, we should take in account that, in the small basis, both the low- and the high-level SCF equations must be solved with A1, whereas only the low-level SCF iterations have to be converged with A2. With the large basis set, however, altogether three Fock-matrix builds—one with n and two with nA electrons—are necessary for A1, while a Fock-matrix of nA electrons has to be constructed in each SCF cycle with A2. This suggests that A2 is only competitive if the embedded subsystem is small with respect to the entire system but turns expensive when enlarging the active subsystem, even though the energy converges to a more accurate value. Concerning our WFT-in-DFT method, the DB approximation only reduces the costs of the SCF part but those of the subsequent correlation calculation are not affected. Thus, the DB approximation can efficiently decrease the expenses of the WFT-in-DFT embedding calculation if the SCF run is a bottleneck, that is, the size of subsystem A is small with respect to that of the entire system. As mentioned above, there are further possibilities to introduce DB embedding Ansätze which are compatible with our DB approach and the Huzinaga-embedding. For example, one could directly project the small-basis DAB density onto the large basis, build the F2 Fock matrix from the projected and purified density matrix, and diagonalize it. The resulting B orbitals could be localized to obtain orbitals φA 1i and φ1i , and then the calculations could

proceed along the route similar to either A1 or A2. Alternatively, the small-basis orbitals projected into the large basis could also be localized before constructing the corresponding density matrices. We did not explore these possibilities since these are all more complicated than the two Ansätze presented above, which seem the simplest combinations of the Huzinaga-embedding and the DB approximations. Finally, we also note that the DB embedding schemes introduced can be combined with

14


Page 14 of 54

Page 15 of 54 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 mixed-basis approach. In the resulting dual mixed basis (DMB) scheme the small basis set of the DB embedding calculation is also the small basis of the mixed-basis approximation, while the large basis set is a mixed basis with the same small basis for the environment and a bigger one for the active subsystem. The major advantage of this scheme is that the initial SCF calculation is carried out with a homogeneous basis, and the density obtained will not be artificially polarized due to the mixing of basis sets. This density is frozen in the subsequent steps of the embedding calculation and expected to generate a more balanced embedding potential than the ones computed with mixed basis sets. Moreover, these steps of the embedding calculation will be faster because of the decreased size of the basis set.

2.4

Further embedding approaches

In our previous publication 80 we also proposed a simple WFT-in-WFT multi-level local correlation approach based on our local correlation methods. 78,87–90 In a local correlation calculation the MOs are localized, and a domain of LMOs is assembled for each occupied LMO. The approximate contribution of the latter to the correlation energy is evaluated in this domain. Additionally, pair domains are constructed for each pair of occupied LMOs, and the long-range correlation is estimated as the sum of the approximate pair correlation energies calculated within these domains. Our multi-level local correlation scheme utilizes this partitioning of the correlation energy: the occupied LMOs are classified as active or inactive as described in Sect. 2.2, and the correlation contributions are evaluated at different levels of theory depending on the type of the corresponding LMOs. For example, in a LNO-CCSD(T)-in-LMP2 embedding calculation, the active subsystem is treated at the CCSD(T) level using local natural orbitals [LNO-CCSD(T)], 87,88,90 while our local MP2 (LMP2) model 89 is applied to the environment. This WFT-in-WFT approach can also be combined with our Huzinaga-embedding outlined in Sect. 2.2 to define WFT-in-WFT-inDFT embedding schemes, which will be referenced as the Huzinaga multi-level local correlation approximation. In these three-layer approaches, first, a HF-in-DFT calculation is 15



performed, and the WFT-in-WFT energy is calculated in the basis of the HF orbitals keeping the KS orbitals of the environment frozen. We also note that the embedding approaches discussed herein can also be integrated with the quantum mechanics/molecular mechanics (QM/MM) method 91 to introduce further three- or four-level embedding schemes, such as DFT-in-DFT-in-MM, WFT-in-DFT-in-MM, or WFT-in-WFT-in-DFT-in-MM. The DB approximation can also accelerate the above embedding methods and the multilevel local correlation schemes as well. In the case of WFT-in-WFT approaches the original DB scheme of Sect. 2.1 can be used for the HF calculation. This is only efficient if the HF is the rate-determining step, e.g., if a LNO-CCSD(T)-in-LMP2 calculation is carried out, and the embedded subsystem is relatively small. To speed up the three- or four-layer embedding calculations the DB method of Sect. 2.3 is applicable. Again, for WFT-in-WFT-inDFT, WFT-in-DFT-in-MM, or WFT-in-WFT-in-DFT-in-MM embedding considerable gain in speed can only be expected if the subsystem treated at the highest level is compact. Finally, the DMB approach can also be combined with our multi-layer methods. Its major advantage is that, in the case of WFT-in-WFT or WFT-in-WFT-in-DFT embeddings, it speeds up not only the SCF iterations but also the correlation calculations.

3

Benchmark calculations

In this section, after the technical details, first, the accuracy of our DB scheme is discussed briefly with an emphasis on choosing the small basis set. Second, the results of benchmark DFT embedding calculations are presented using our two DB set approaches, namely DB-A1 and DB-A2, and their DMB set alternatives, DMB-A1 and DMB-A2. Finally, the performance of the WFT-in-DFT, WFT-in-WFT, and WFT-in-WFT-in-DFT methods is assessed employing the corresponding DB approaches. The accuracy and efficiency of these techniques are compared to the high-level methods using the large basis set, as well as to the original Huzinaga-embedding employing the large basis and mixed-basis sets. We note that

16


Page 16 of 54

Page 17 of 54 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 efficiency of the various techniques will be characterized here by the size of the active subsystem required for a particular accuracy. The efficiency of the methods will also be quantified by presenting factual computation times for a larger example in Sect. 4.

3.1

Computational details

All the schemes presented here have been implemented in the Mrcc program suite 92 and are available in its current release. In our calculations Pople’s double- (6-31G, 6-31G*, 6-31G**, 6-31+G*, 6-31+G**, 631++G**) and triple-ζ (6-311G, 6-311G*, 6-311+G*, 6-311G**, 6-311+G**, 6-311++G**) basis sets 93,94 as well as Dunning’s (augmented) correlation consistent polarized double-, triple-, and quadruple-ζ bases [(aug-)cc-pVDZ, (aug-)cc-pVTZ, and (aug-)cc-pVQZ] 95–97 were used. Both the SCF and the correlation calculations employed the density fitting approximation using the auxiliary bases developed by Weigend and co-workers. 98,99 In all calculations spherical harmonic Gaussians were employed. The DFT-in-DFT, WFT-in-DFT, and WFT-in-WFT-in-DFT embedding methods were tested with three different exchange-correlation functionals, which were chosen in order to represent the various rungs on the Jacob’s ladder of density functionals. The simplest functional of them is the local (spin) density approximation (LDA), 100,101 which only depends on the electron density itself. An improved functional compared with LDA is the Perdew– Burke–Ernzerhof (PBE) 102 functional of the general gradient approximation (GGA) class, which also depends on the gradient of density. The most complete functional in our tests was the hybrid GGA functional B3LYP (Becke’s three-parameter hybrid functional including the correlation functional of Lee, Yang, and Parr). 103,104 As the exchange-correlation functional becomes more sophisticated, the accuracy, in general, also improves. However, this comes with increased computational expenses: B3LYP is much more expensive because of the exact exchange term compared with the PBE and LDA functionals, which have more or less similar costs. Since we have found that the convergence of the errors with the subsystem size is 17



similar for these functionals with the Huzinaga-embedding scheme, only the B3LYP-in-LDA, WFT-in-PBE, and WFT-in-WFT-in-PBE results are shown in the main text. The rest of the results can be found in the supplementary material. The WFT-in-DFT, WFT-in-WFT, and WFT-in-WFT-in-DFT computations used the LNO-CCSD(T) 87,88,90 method as the high-level approach with tight truncation thresholds except for aromatic systems, where the default thresholds were applied. In the multi-level local correlation calculations our LMP2 approximation 89 was invoked as the low- or mediumlevel method. The Pipek–Mezey 105 localization scheme was employed for the occupied orbitals throughout. Unless otherwise noted, the aMul orbital selection technique (see the supplementary material) was utilized with the default = 0.30 population threshold. The improved πMul algorithm was used for the aromatic systems with the κ = 0.95 completeness criterion and the default population threshold for the σ orbitals. The benchmark systems were taken from Ref. 73: the substitution reaction of 1-chlorodecane with a hydroxide ion to yield 1-decanol and a chloride ion, the deprotonation of decanoic acid, the Diels–Alder reaction between the conjugated octadecanonaene and 1,3butadiene (see Fig. 2d of Ref. 73), and the hydrogenation of pentacene (see Fig. 3d of Ref. 73). These test systems will shortly be referenced as a, b, c, and d, respectively. The corresponding geometries and the detailed description of the definition of the embedded subsystems can be found in our previous paper. 80 The reference reaction energies for our benchmark tests were calculated using the highlevel methods with the large basis sets. The reference energies are subtracted from the reaction energies computed by the embedding techniques, and the errors are plotted as a function of the embedded carbon atoms, except for the Diels–Alder reaction where the abscissa represents the number of embedded carbon atoms of octadecanonaene. In the limit of zero embedded carbon atoms, the convention of Ref. 73 was adopted. This means that the presented errors are the errors of the reaction energies calculated with the low-level method

18


Page 18 of 54

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


using the basis set compatible with the embedding scheme, namely, the small basis for the DMB and mixed-basis variants, the DB with the DB Ansätze, and the large basis for the original Huzinaga-embedding or multi-level local correlation approaches. For the WFT-inWFT-in-DFT methods, the HF and DFT regions at the initial HF-in-DFT calculations were kept fixed, and only the border of the two correlation approaches was varied. At DMB and mixed-basis WFT-in-WFT-in-DFT calculations, the large basis set was only used for the region treated by the LNO-CCSD(T) method, while the small basis set was employed for the rest of the molecule.

3.2

Dual basis set approximation

To assess the performance of our DB approach and to select the possible small basis and large basis pairings we carried out benchmark calculations for the test reactions with the B3LYP, PBE, and LDA DFT approaches and the LNO-CCSD(T) method without any embedding approximation. Since in our DB scheme the small basis is not necessarily the subspace of the large basis set, one can choose an appropriate small basis from various basis set families. For the DFT theories various combinations of the Pople basis sets were tested, while for the LNOCCSD(T) method Dunning’s (aug-)cc-pVTZ and aug-cc-pVQZ basis sets were employed as the large basis set combined with several Pople and Dunning basis sets as the small basis to investigate the effects of the DB approximation on the correlation energies. The basis sets augmented with diffuse functions were used as the large basis set if the phenomenon considered requires the presence of such functions, that is, for reactions a and b, where one of the reactants or products is an anion. Our findings are in line with the conclusions of Head-Gordon et al., 23,27 who found that accurate reaction energies can be achieved by eliminating the functions with high angular momenta from the large basis set. As an example, the errors for the B3LYP and LNOCCSD(T) approaches with the DB approximation for the deprotonation of decanoic acid are shown in Fig. 1, while the other tests can be found in the supplementary material. 19



3.0 Error [kcal/mol]

B3LYP / 6-311++G**

Error of KS energy (product)

2.0

Error of KS energy (reactant) Error of reaction energy

1.0 0.0 -1.0

G -31

6

** * * ** * * ** G* G G* 1+G* 11+G 11G* 1++G 11+G G 1 + G 1 1 1 1 1 1 6-3 6-3 6-3 6-3 6-3 6-3 6-3 6-3 6-3 6-3

14.0 12.0

Error of HF energy (product)

LNO-CCSD(T) / aug-cc-pVTZ

Error of HF energy (reactant)


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

Page 20 of 54

Error of correlation energy (product)

8.0

Error of correlation energy (reactant)

6.0

Error of reaction energy

4.0 2.0 0.0 -2.0 -4.0

* 31G

6-

6-3

1

* +G

6-

** 31G

p cc-

Z VD

6-3

* 11G

6-3

* +G

1

* 6-3

1

G ++

** 6-3

* 11G

Z

VD

*

c-p g-c

au

Figure 1: Error of B3LYP (top) and LNO-CCSD(T) (bottom) reaction energies calculated with the DB approximation using various small basis sets. The KS, HF, and correlation energy errors are decomposed for the educt and the product. The 6-311++G** and augcc-pVTZ large basis sets were used for the B3LYP and the LNO-CCSD(T) calculations, respectively. The large basis B3LYP and LNO-CCSD(T) reference reaction energies are 350.63 and 359.03 kcal/mol, respectively.

20


Page 21 of 54 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 results for B3LYP reveal that the DB approach can reproduce the reaction energies within a half kcal/mol when at least a double-ζ basis set with polarization functions on heavy atoms is employed. Moreover, the errors can be reduced to a few tenths of a kcal/mol with the addition of diffuse functions on heavy atoms, and the accuracy can be further improved to a couple of hundredths of a kcal/mol when the aforementioned functions are also placed on the hydrogens. These findings suggest that diffuse functions in the small basis are less important for DB DFT methods. However, we should keep in mind that the electronic structure of the decanoate anion is less diffuse as the negative charge is delocalized on the carboxylate group. For small anions with rather diffuse electron distribution, such as the chloride and hydroxide ions, diffuse functions are essential in the small basis, which can be seen from the results for test reaction a (see the supplementary material). The LNO-CCSD(T) DB calculations require at least a double-ζ small basis set with polarization and diffuse functions on heavy atoms for an accuracy of less than a half kcal/mol in the reaction energies. Similarly to the B3LYP results, the addition of the aforementioned functions to hydrogens reduces the reaction energy errors to a couple of hundredths of a kcal/mol. This trend holds on for other reactions as well, which means that the correlated calculations seem to be more sensitive to diffuse functions. If the reactants or products of the reaction require diffuse functions, the small basis must also include them at least on the heavy atoms for a tolerable error. The decomposition of the reaction energy errors also reveals that individual HF energies may be more accurately approximated using a smaller basis set (for example, 6-31+G**) than a bigger one (aug-cc-pVDZ), however, the error cancellation is more beneficial when the small basis and large basis are from the same basis set family. It is also interesting that the error of the correlation energy is less dependent on the small basis when polarization functions are placed at least on the heavy atoms, which can be rationalized by the fact that the correlation energy is more sensitive to the underlying AO basis than the quality of the occupied orbitals. The conclusions for the DB LNO-CCSD(T) results with the aug-cc-pVQZ large basis set

21



are similar to the triple-ζ case (see the supplementary material) since the smallest errors are obtained when aug-cc-pVTZ is used as the small basis set. The calculations employing the aug-cc-pVDZ basis as the small basis set are also remarkably accurate, and the reaction energy errors are still a kcal/mol or smaller even with the 6-31+G* small basis set. It is also instructive to compare the results of the DB calculations and the simple calculations which are only performed with the corresponding small basis set without any DB approximation. Our results suggest (see the supplementary material) that the DB calculations are roughly by an order of magnitude more accurate than the corresponding small basis only calculations if reliable small bases are used. In the light of these findings we choose the double-ζ counterparts of the triple-ζ large basis sets as small basis in the following to maximize the accuracy of our DB embedding ansatz.

3.3

DFT embedding

Our benchmark calculations were designed with the purpose to explore the convergence behavior of the embedding approaches with respect to the size of the embedded subsystem. Please note that the considered reactions are more and more difficult for the embedding theories. The simplest test system is reaction a, where the subsystems are separated across a singe σ bond. This is also the case for reaction b, which is, however, more difficult because the error cancellation is less significant due to the differently charged species. Dividing a system through conjugated bonds is modeled by reaction c, and the most difficult one is reaction d, where an aromatic system is broken. The errors for the B3LYP-in-LDA calculations are presented in Fig. 2. For better visualization, the results calculated with ansatz A2 are not shown because they are very close to those obtained by A1. The A2 results together with the details of the DFT embedding calculations employing other functionals can be found in the supplementary material. The DB error for the high-level method, which is equivalent to the error of the DB 22


Page 22 of 54

Page 23 of 54

10.0 a) Substitution reaction of chlorodecane


Error [kcal/mol]

2.0 1.0 0.0 -1.0

Error [kcal/mol]

0.0

-2.0

2.5 0.0

0 1 2 3 4 5 6 7 8 9 10 Number of carbon atoms in the embedded subsystem 5.0 d) Hydrogenation of pentacene 2.5

c) Diels-Alder reaction of octadecanonane

1.0

-1.0

5.0

-5.0

0 1 2 3 4 5 6 7 8 9 10 Number of carbon atoms in the embedded subsystem 2.0

b) Deprotonation of decanoic acid

-2.5

-2.0

Error [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


Huzinaga-LB Huzinaga-MB Huzinaga-DB-A1 Huzinaga-DMB-A1

0.0 -2.5 -5.0 -7.5

0 2 4 6 8 10 12 14 16 18 Number of carbon atoms in the embedded subsystem

-10.0

0 2 4 6 8 10 12 14 16 18 20 22 Number of carbon atoms in the embedded subsystem

Figure 2: Error of the reaction energies using various embedding approaches as a function of the number of carbon atoms included in the embedded subsystem. The LB and MB acronyms stand for large basis and mixed basis, respectively. The 6-311++G** (a), 6-311++G** (b), 6-311G** (c), and 6-311G** (d) large basis sets were used in the calculations, while the corresponding small basis sets in the mixed-basis, DB, and DMB calculations were 631+G**, 6-31+G**, 6-31G*, and 6-31G*, respectively. The large basis reference reaction energies are −56.71, 350.63, −16.05, and 34.76 kcal/mol, respectively.

23



Huzinaga-embedding at the maximum number of carbon atoms, is less than 0.05 kcal/mol for all the reactions, which suggests that the selection of the small basis sets is appropriate. The results show that the errors of the DB, DMB, mixed-basis, and large basis schemes are similar, except for reactions a and b at smaller embedded subsystems. For the former reaction the DB and DMB variants produce errors of 1.0 to 1.5 kcal/mol larger compared to the parent method for the smaller subsystems. The opposite trends can be observed for test reaction b, where the DB and DMB approaches yield 1.5 to 2.0 kcal/mol smaller errors compared with the large basis approach, but the difference of the errors vanishes again as the embedded subsystem enlarges. These findings show that the error cancellation may have an important role in the DB calculations with small active regions but is not consistent. Nevertheless, the absolute errors are still moderate even at the small embedded subsystems, and the various DB schemes perform similarly well when the subsystem boundary is 5 to 6 Å away from the reaction center. Note that this distance approximately coincides with the distance where the error fluctuations of all the embedding variants disappear or reduces to a low level. Interestingly, for the complicated test systems, reactions c and d, the DB Huzinaga-embedding is practically indistinguishable from the parent Huzinaga-embedding method for any reasonable embedded subsystem. The simple mixed-basis approach performs remarkably well for the considered reactions, especially for reactions c and d, but its slower convergence is conspicuous for reaction b, where a charged species is treated with the mixedbasis. The performance of the DMB approach is as good as that of the pure DB scheme. The errors obtained with the DB-A1 and DB-A2 methods are very similar, and this trend can also be observed for the performances of the DMB-A1 and DMB-A2 techniques (see the supplementary material). Taking into account that the A2 ansatz is more expensive because a high-level SCF calculation is carried out with the large basis, the A1 ansatz is clearly more beneficial. All in all, considering also the computational costs, all the approximate embedding schemes are more efficient than the original Huzinaga-embedding. Using larger embedded

24


Page 24 of 54

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


subsystems, which are anyway required for accurate results, the DB approach is at least as accurate as the mixed-basis scheme. However, increasing the size of the embedded subsystem the costs of the latter are getting more expensive, and the DB approach tends to be more efficient. The most efficient variant seems to be the DMB approach, which inherits the accuracy of the DB technique but is less expensive than both the DB and the mixed-basis methods.

3.4

WFT embedding

In this subsection, the results of LNO-CCSD(T)-in-PBE, LNO-CCSD(T)-in-LMP2, and LNO-CCSD(T)-in-LMP2-in-PBE calculations using the large basis as well as the mixedbasis, DB, and DMB approximations are presented for our benchmark systems. The results of further embedding calculations with the LDA and B3LYP functionals can be found in the supplementary material. The convergence of the Huzinaga multi-level local correlation method was tested with a fixed DFT and HF partitioning, and only the LNO-CCSD(T) and LMP2 regions were varied. The HF and DFT subsystems were chosen so that the size of the HF region will be minimal while also keeping the error of the LNO-CCSD(T)-in-DFT method using the same HF region as small as possible. In the case of Huzinaga multi-level local correlation with DMB or mixed-basis calculations, the large basis and small basis boundary also varies, which means that the large basis set is only placed on the LNO-CCSD(T) region, whereas the small basis set is used for the rest of the system. The results of the various embedding schemes for the benchmark systems a, b, c, and d are presented in Figs. 3, 4, 5, and 6, respectively. Note that the accuracy of the DB approximation for the LNO-CCSD(T) method without any embedding (see the caption of the figures) somewhat decreases compared to the DFT results for test reactions a and d. The increased error of the DB ansatz is likely due to the greater sensitivity of the HF method to the quality of the small basis in the case of challenging electronic structures, such as aromatic systems and anions with diffuse charge distributions (cf. Fig. 1). Nevertheless, these errors 25



are still tolerable and could be consistently decreased with a less aggressive truncation of the large basis set. 23,27 Despite the worse performance of the DB approximation in some cases, the errors of the Huzinaga embedding with the DB and DMB schemes follow similar trends as we have seen for the DFT embedding. The DB and DMB approximations are capable of reproducing the reaction energies computed with the large basis scheme within 1 kcal/mol for test systems a and c, while even better convergence with respect to the active subsystem size can be seen for reaction b. For reaction d, the DB (DMB) technique performs worse by 1.0 (3.0 to 7.0) kcal/mol compared with the large basis scheme if the active subsystem is small, but the differences vanish (or significantly decrease) when the embedded subsystem consists of at least 3 carbon rings. This active subsystem size is also necessary for the original method to produce meaningful reaction energies. In contrary to the DFT results, the mixed-basis technique is unreliable for reactions a and b because large error fluctuations can be observed even at extended embedded subsystems. This is a consequence of the fact that an unstable HF-in-PBE reference is found in particular points with the default settings of the SCF algorithm. This phenomenon disappears for the most difficult test systems, reactions c and d, where no charged species are included in the reaction, but also less accurate results are produced compared to the large basis method. The multi-level local correlation approach is not only more accurate than the Huzinaga embedding but also more expensive, which makes the application of DB approaches even more desirable. The results show that the accuracy of the multi-level local correlation approaches is practically unaffected for any test reaction by the use of the DB approximation. The reaction energies of the DB and DMB schemes are within 0.2 kcal/mol of the large basis and mixed-basis approaches, respectively, independently of the size of the embedded subsystem. The combination of the DFT embedding and the multi-level local correlation approach has even more potential for accurate and efficient calculations. The results of the Huzinaga

26


Page 26 of 54

Page 27 of 54

Error [kcal/mol]

3.0

LNO-CCSD(T)-in-PBE


2.0 1.0 0.0 -1.0 -2.0 -3.0

0

1

2

3

4

5

6

7

8

9

10

Error [kcal/mol]

3.0

LNO-CCSD(T)-in-LMP2

2.0 1.0 0.0

MLC-LB MLC-MB MLC-DB MLC-DMB

-1.0 -2.0 -3.0

0

1

2

3

4

5

6

7

8

9

10

3.0

Error [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


LNO-CCSD(T)-in-LMP2-in-PBE

2.0 1.0 0.0

HMLC-LB HMLC-MB HMLC-DB-A1 HMLC-DMB-A1

-1.0 -2.0 -3.0

0

1

2

3

4

5

Number of carbon atoms in the embedded subsystem Figure 3: Error of the reaction energies for the substitution reaction of 1-chlorodecane using various embedding approaches as a function of the number of carbon atoms included in the embedded subsystem. The MLC abbreviation refers to the multilevel local correlation approach, while HMLC means the Huzinaga-embedding based MLC scheme. See Fig. 2 for further abbreviations. The aug-cc-pVTZ and aug-cc-pVDZ basis sets were used as large basis and small basis sets, respectively. In the Huzinaga-embedding based multilevel local correlation calculations the boundary of DFT and HF regions was fixed at the 5th carbon atom. The reference LNO-CCSD(T)/aug-cc-pVTZ energy is −52.12 kcal/mol, and the error of the DB approximation is −0.18 kcal/mol.

27


Error [kcal/mol] Error [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

Error [kcal/mol]


4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 -4.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 -4.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 -4.0

Page 28 of 54

LNO-CCSD(T)-in-PBE


0

1

2

3

4

5

6

7

8

9

10

LNO-CCSD(T)-in-LMP2


0

1

2

3

4

5

6

7

8

9

10



0

1

2

3

4

5

6

Number of carbon atoms in the embedded subsystem Figure 4: Error of the reaction energies for the deprotonation of decanoic acid using various embedding approaches as a function of the number of carbon atoms included in the embedded subsystem. See Figs. 2 and 3 for the abbreviations. The aug-cc-pVTZ and aug-cc-pVDZ basis sets were used as large basis and small basis sets, respectively. In the Huzinaga-embedding based multilevel local correlation calculations the boundary of DFT and HF regions was fixed at the 6th carbon atom. The reference LNO-CCSD(T)/aug-ccpVTZ energy is 359.03 kcal/mol, and the error of the DB approximation is −0.01 kcal/mol.

28


Page 29 of 54

Error [kcal/mol]

6.0

LNO-CCSD(T)-in-PBE


4.0 2.0 0.0 -2.0 -4.0 -6.0

0

2

4

6

8

10

12

14

16

18

Error [kcal/mol]

6.0

LNO-CCSD(T)-in-LMP2

4.0 2.0 0.0


-2.0 -4.0 -6.0

0

2

4

6

8

10

12

14

16

18

6.0

Error [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



4.0 2.0 0.0


-2.0 -4.0 -6.0

0

2

4

6

8

10

12

Number of carbon atoms in the embedded subsystem Figure 5: Error of the reaction energies for the Diels–Alder reaction of octadecanonaene and 1,3-butadiene using various embedding approaches as a function of the number of carbon atoms included in the embedded subsystem. See Figs. 2 and 3 for the abbreviations. The ccpVTZ and cc-pVDZ basis sets were used as large basis and small basis sets, respectively. In the Huzinaga-embedding based multilevel local correlation calculations the boundary of DFT and HF regions was fixed at the 4th carbon atom from the reaction center (that is, 12 carbon atoms were included in the embedded subsystem). The reference LNO-CCSD(T)/cc-pVTZ energy is −35.55 kcal/mol, and the error of the DB approximation is −0.05 kcal/mol.

29


Error [kcal/mol] Error [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

Error [kcal/mol]


15.0 10.0 5.0 0.0 -5.0 -10.0 -15.0 -20.0 -25.0 15.0 10.0 5.0 0.0 -5.0 -10.0 -15.0 -20.0 -25.0 15.0 10.0 5.0 0.0 -5.0 -10.0 -15.0 -20.0 -25.0

Page 30 of 54

LNO-CCSD(T)-in-PBE


0

2

4

6

8

10

12

14

16

18

20

22

LNO-CCSD(T)-in-LMP2


0

2

4

6

8

10

12

14

16

18

20

22



0

2

4

6

8

10

12

14

Number of carbon atoms in the embedded subsystem Figure 6: Error of the reaction energies for the hydrogenation of pentacene using various embedding approaches as a function of the number of carbon atoms included in the embedded subsystem. See Figs. 2 and 3 for the abbreviations. The cc-pVTZ and cc-pVDZ basis sets were used as large basis and small basis sets, respectively. In the Huzinaga-embedding based multilevel local correlation calculations the boundary of DFT and HF regions was fixed at the 6th carbon atom from the reaction center (that is, 14 carbon atoms of 3 aromatic rings were included in the embedded subsystem). The reference LNO-CCSD(T)/cc-pVTZ energy is 26.02 kcal/mol, and the error of the DB approximation is −0.36 kcal/mol.

30


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


multi-level local correlation technique show that, if the HF and DFT subsystems are carefully chosen, then for any basis-set combination the LNO-CCSD(T) region has to cover only 2 to 3 bonds from the reaction center to reproduce the corresponding LNO-CCSD(T)-in-DFT reaction energy within 0.5 kcal/mol. The only exception is reaction d, where larger LNOCCSD(T) layers are needed. The convergence pattern of the errors with the size of the LNO-CCSD(T) region is reminiscent of that for the multi-level local correlation approach, excepting the mixed-basis scheme, which inherits the large error fluctuations from the parent Huzinaga-embedding. These tests indicate that it is unnecessary to correlate all the occupied orbitals of the HF region at the LNO-CCSD(T) level as the LMP2 subsystem serves as a cost-effective intermediate layer. In conclusion, the WFT embedding test calculations revealed that the DB and DMB schemes, either with the A1 or the A2 approaches, are capable of accelerating the Huzinaga calculations over the original large basis schemes, while the mixed-basis technique shows unacceptable inaccuracies in specific cases. Again, the A2 approach is less favorable over ansatz A1 because no improvement is provided in the accuracy in exchange for the increased computational cost. All the basis set reduction approaches considered here effectively accelerate the multi-level local correlation calculations, however, the DMB method is the clear winner because it saves computation time not only for the HF step but also for the correlation calculation. The most promising technique is the combination of the Huzinaga DFT embedding Ansatz A1 with the DMB and multi-level local correlation, which effectively breaks down the cost of all computational steps while maintaining high accuracy.

4

Application: Catalytic reaction of Catechol-O-Methyltransferase

To check if the conclusions drawn for the small, quasi-one-dimensional molecules are valid for real-life 3D systems and to demonstrate the applicability of our multi-level approaches, 31



in this section, the performance of the methods is tested for an enzyme reaction of biological relevance. In these calculations our methods are also combined with the QM/MM technique. The projector-based embedding approach of Manby and co-workers was already applied to such systems, 106,107 however, this is the first time when the accuracy of the WFT-in-WFT-inMM and WFT-in-WFT-in-DFT-in-MM approximations as well as multi-layer DB techniques are explored with respect to the active subsystem size. Here, we choose the catechol-Omethyltransferase (COMT) enzyme as our model system, which has an important role in neurotransmitter regulation by transferring a methyl group from S-adenosyl-L-methionine (SAM) to deprotonated catechol (CAT) amines. An important feature of the COMT system is that the reacting species, the SAM and the CAT molecules, are noncovalently bound to the protein, thus the boundary effects of the QM/MM technique are small. In the upcoming subsections, the computational details are discussed. Then, the accuracy of our most efficient methods, namely the Huzinaga-A1, the multi-level local correlation method, and the combination of the latter with DMB, are compared to the reference large basis and DB LNO-CCSD(T) method. Finally, timings for the most relevant calculations are presented.

4.1

Computational details

The reactant and the product structure of the COMT system were adopted from the work of Martínez et al. 108 The Antechamber 109 and the Tleap utilities of the AmberTools17 program package 110 were invoked to set up the force field parameters for the MM calculations. The ff14SB 111 parameters were used for the protein, while the TIP3P 112 model was employed for the water and the ions. For the CAT and SAM molecules, the General Amber Force Field 2 (GAFF2) 113 was utilized with the AM1-BCC charge scheme. 114,115 The non-bonded energy terms were not cut off, and only those bonded energy terms were retained that included at least one MM or MM host atom. The residual charges, which are the results of zeroing the charges of the QM and MM host atoms due to electrostatic embedding, were equally 32


Page 32 of 54

Page 33 of 54 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


distributed among all MM atoms. The QM/MM calculations were also carried out with AmberTools17 invoking the Amber-Mrcc interface 91 and using the aforementioned force fields for the MM region. The QM calculations were performed with the LNO-CCSD(T) 87,88,90 and LMP2 89 local correlation methods, and the PBE 102 density functional employing Grimme’s D3 dispersion correction 116,117 together with Dunning’s cc-pVTZ and cc-pVDZ basis sets 95 as the large basis and small basis sets, respectively. Again, the QM calculations employed the density fitting approximation with the auxiliary basis sets of Weigend and co-workers. 98,99 The exchange contributions in the HF calculations were approximated using local fitting domains 89,118 with the default domain settings proposed in Ref. 89. All the calculations used the Boys localization technique, 119 whereas the Huzinaga-embedding, the multi-level local correlation approach, and the Huzinaga multi-level local correlation method employed the aMul orbital selection scheme. The reference reaction energy was computed at the LNO-CCSD(T)-in-MM level with the cc-pVTZ basis set. The presented errors were calculated by subtracting the reference reaction energy from the reaction energies computed with the various approaches. Other conventions that were mentioned in Sect. 3.1 are also adopted here. The selected QM region and QM subsystems are illustrated in Fig. 7 and Fig. 8, respectively, while the details about the selection of the various subsystems can be found in the supplementary material.

4.2

Results

The errors of the various schemes are presented in Fig. 9. The DB approximation proves very accurate for this system as the DB error is only a few hundredths of a kcal/mol. The results obtained with the Huzinaga-embedding follow similar trends as we have seen in Sect. 3.4. The errors are less than 1.5 kcal/mol even if small embedded subsystems are used. Some error fluctuation can be observed in this region for the DMB scheme, but, for most subsystems, the errors are smaller compared to the parent large basis scheme. Interestingly, 33



Figure 7: The COMT protein including 3420 atoms. The QM region (highlighted in red) consists of 36 amino acid residues, 3 water molecules, an Mg2+ ion, the SAM, and the CAT molecules, which results in 571 QM atoms including 18 link atoms. See the supplementary material for more details regarding our selection approach.

34


Page 34 of 54

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


’reac_01_sat.jpg’ binary filetype=jpg



A) 50 atoms (2.0 Å)

B) 64 atoms (2.0 Å)

C) 81 atoms (2.5 Å)

’reac_03b_sat.jpg’ binary filetype=jpg



D) 95 atoms (2.5 Å)

E) 140 atoms (4.5 Å)

F) 224 atoms (5.5 Å)

Figure 8: Visualization of the selected QM subsystems, which are shortly referenced as the corresponding capital letters. The MM, PBE-D3 (LMP2), and LNO-CCSD(T) methods are applied to the regions highlighted in green, blue, and red, respectively, in the case of the Huzinaga (multi-level local correlation) approach. As a reference point, the reacting methyl group is highlighted in purple. The region treated with LNO-CCSD(T) is gradually expanded from the reaction center towards the full QM region (see the supplementary material for details).

35



the errors of the DMB scheme are also much smaller in the case of subsystems A and C compared to the B and D partitionings, which suggests that the error cancellation is more fortunate if the embedded subsystem is more compact. Nevertheless, the errors decrease to a low level if the subsystem border is at least 5.5 Å away from the reaction center, which is in line with our previous findings. It is also important to note that the neglect of the D3 correction (not shown) increases the errors by about 1 kcal/mol for the small embedded regions, which proves the importance of the inclusion of dispersion interactions for embedding calculations even at the low-level method. The performance of the multi-level local correlation scheme with the large basis set is excellent, its maximum error is very low, 0.3 kcal/mol, for any embedded subsystem, and the errors monotonically decrease beyond subsystem C. On the other hand, moderate error fluctuations occur in the case of the multi-level local correlation technique with DMB, but the approach is still very accurate as the maximum error does not exceed 0.5 kcal/mol. These attributes are also inherited by the Huzinaga multi-level local correlation method because their errors approximately equal those of the corresponding multi-level local correlation schemes shifted by the error of the Huzinaga-embedding with active subsystem F. The errors of the Huzinaga multi-level local correlation methods are still within 0.5 kcal/mol compared to the reference reaction energy, which shows that accurate reaction energies can be achieved by applying the highest-level technique only to a small region around the reaction center. These results, in accordance with the conclusions for the small molecules, reveal that all the considered methods have tolerable accuracy even with relatively small embedded subsystems. Note also that the subsystem boundaries cross only one or two bonds for the small molecules, while here the subsystems are separated by much more bonds, but the error does not increase with the number of such bonds. The DFT/HF border is recommended to be at least 5 Å away from the reaction center in the case of the Huzinaga embedding to minimize the ambiguity stemming from the error cancellation effects and to produce close to converged results. The errors of the Huzinaga multi-level local correlation method converge faster with

36


Page 36 of 54

Page 37 of 54

Error [kcal/mol]

1.5

LNO-CCSD(T)-in-PBE-D3-in-MM

1.0 0.5 0.0 -0.5

Huzinaga-LB Huzinaga-DMB

-1.0 -1.5

0

50 100 150 200 250 300 350 400 450 500 550 600

Error [kcal/mol]

1.5

LNO-CCSD(T)-in-LMP2-in-MM

1.0 0.5 0.0 -0.5

MLC-LB MLC-DMB

-1.0 -1.5

0

50 100 150 200 250 300 350 400 450 500 550 600

1.5

Error [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


LNO-CCSD(T)-in-LMP2-in-PBE-D3-in-MM

1.0 0.5 0.0 -0.5

HMLC-LB HMLC-DMB

-1.0 -1.5

0

25

50

75

100 125 150 175 200 225 250

Number of QM atoms in the embedded subsystem Figure 9: Error of the reaction energies for the methyl transfer reaction catalyzed by the COMT enzyme using various multi-level approaches as a function of the number of QM atoms included in the embedded subsystem. See Figs. 2 and 3 for the abbreviations. The cc-pVTZ and cc-pVDZ basis sets were used as large basis and small basis sets, respectively. In the Huzinaga-embedding based multilevel local correlation calculations the boundary of the DFT and HF regions was fixed, and the latter included the atoms of subsystem F. The reference LNO-CCSD(T)/cc-pVTZ reaction energy is −23.18 kcal/mol, and the error of the DB approximation is −0.06 kcal/mol.

37



respect to the subsystem size, but a domain of radius of 2.5 Å around the reaction center still has to be treated with the high-level method for accurate results. Taking into account these findings, to demonstrate what we can gain with the various approximations, computation times were measured for the Huzinaga and multi-level local correlation approaches using the F and C partitionings, respectively. The timings for the Huzinaga multi-level local correlation approach were determined employing the LMP2 and LNO-CCSD(T) methods for subsystems F and C, respectively. The timings of the various schemes are presented in Table 1. The results show that a speedup of 1.6 can be achieved for the reference HF calculation by using the DB approximation without embedding approaches. At first glance this can be surprising because the cc-pVTZ basis set consists of more than twice as many AOs as the cc-pVDZ basis. The reason for the relatively low efficiency of the DB approximation is that small local fitting domains are used during the SCF iterations, and finally, to get an accurate HF energy, one iteration is performed with much larger fitting domains. 89,118 This final step usually requires computation time comparable to that for the preceeding SCF iterations altogether. In a DB HF calculation employing the local fitting approximation the single iteration in the large basis set must be carried out with the large fitting domains, thus only the expenses of the SCF iterations performed with small fitting domains can be saved. Nevertheless, the speedup for the HF step can be increased by the DMB ansatz to 2.1, while the speedup for the SCF part is 3.2 if the DMB HF approach is also combined with DFT via the Huzinaga-embedding. Concerning the overall speedup of the local correlation calculation, we can see that the DB approximation alone does not significantly affect the costs of the calculations, though it only results in a negligible loss of accuracy. The speedups are more pronounced for the multi-level schemes at the expense of 0.3 kcal/mol loss in accuracy. The DMB approach roughly doubles the speedup with respect to the corresponding multi-level scheme applying the large basis set. The most efficient but reliable method is clearly the Huzinaga multi-level local correlation technique, followed by the Huzinaga and the multi-level local correlation approach.

38


Page 38 of 54

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


We note that the potential speedups for the aforementioned methods can be even greater if larger errors can be afforded and a smaller domain is selected for the highest-level method. Table 1: Representative timings for the various multi-level approaches. See Figs. 2 and 3 for the abbreviations. The reported computation times are wall-clock times measured on a 6-core 3.5 GHz Intel Xeon E5-1650 processor. The total computation times are decomposed into the DFT, HF, and correlation contributions. The errors and speedups are given with respect to the LNO-CCSD(T)-in-MM method using the cc-pVTZ basis set. See Figs. 2 and 3 for the abbreviations. Method Basis set approach HF region CCSD(T) region DFT [day] HF [day] Correlation [day] Total [day] Speedup Error [kcal/mol]

5

LNO-CCSD(T) LB full full − 2.64 3.65 6.29 − −

LNO-CCSD(T) DB full full − 1.61 3.65 5.26 1.20 −0.06

Huzinaga LB F F 0.43 1.32 1.08 2.83 2.23 0.15

Huzinaga DMB F F 0.16 0.66 0.90 1.71 3.68 −0.29

MLC LB full C − 2.64 1.16 3.81 1.65 −0.30

MLC DMB full C − 1.24 0.61 1.85 3.39 0.00

HMLC LB F C 0.43 1.32 0.67 2.41 2.61 −0.11

HMLC DMB F C 0.16 0.66 0.39 1.21 5.20 −0.21

Conclusions

In this paper, we have presented an alternative dual basis ansatz, which performs similarly to the previous approaches when applied to conventional electronic structure methods but is suitable for the generalization to embedding schemes. Relying on the new ansatz two possible DB approximations have been proposed for the projector-based embedding method using the Huzinaga-equation. We have demonstrated that the DB approximation considerably increases the efficiency of embedding calculations in exchange for negligible loss of accuracy and also outperforms the notorious mixed-basis set approach. However, if the DB and mixed-basis approximations are combined, the resulting DMB scheme inherits the advantageous features of both and further improves the efficiency for both the WFT and the DFT embeddings. The DB and DMB schemes have also been tested for our multi-level correlation approach, where various local correlation methods are used for the different regions of the system. Our results show that the accuracy of the multi-level local correlation method is practically unaffected by the basis set approximations. The combination of the

39



Huzinaga-embedding and the multi-level local correlation schemes has also been introduced, namely the Huzinaga multi-level local correlation approach. The benchmark calculations show that it is sufficient to only treat a small subsystem of the molecule at the CCSD(T) level if an MP2 layer is also inserted between the CCSD(T) and the DFT regions. For the CCSD(T) layer, if no aromatic system is broken, it is sufficient to include atoms separated by 2-3 bonds from the reactive center of the molecule, while the MP2 region is recommended to consist of further atoms that are separated by 5-6 bonds from the bonds affected in the reaction. In addition, the DMB approach can also be safely employed together with the Huzinaga multi-level local correlation method using the large basis set only for the atoms of the CCSD(T) layer. The proposed approaches have also been combined with QM/MM techniques and tested for the biologically relevant COMT system. The trends observed for the enzyme reaction verify our conclusions drawn for the small benchmark systems. This application also demonstrates that the proposed methods are useful for the accurate treatment of large molecular systems including thousands of atoms.

Supplementary material See supplementary material for molecular structures, detailed results of DB calculations, DB embedding calculations with various functionals. This information is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgments The authors are grateful for the financial support from the National Research, Development, and Innovation Office (NKFIH, Grant No. KKP126451). This work was also supported by the BME-Biotechnology FIKP grant of EMMI (BME FIKP-BIO). The work of PRN is supported through the New National Excellence Program of the Ministry of Human Capacities, ID: ÚNKP-17-4-BME-55. The computing time granted on the Hungarian HPC 40


Page 40 of 54

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


Infrastructure at NIIF Institute, Hungary, is gratefully acknowledged.

References (1) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Electronic Structure Theory; Wiley: Chichester, 2000. (2) Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618. (3) Gauss, J. In Encyclopedia of Computational Chemistry; Schleyer, P. R., Jorgensen, W. L., Schaefer III, H. F., Schreiner, P. R., Thiel, W., Eds.; Wiley: New York, 1998, p. 615. (4) Purvis III, G. D.; Bartlett, R. J. A full coupled-cluster singles and doubles model: The inclusion of disconnected triples. J. Chem. Phys. 1982, 76, 1910. (5) 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. (6) Bauschlicher Jr., C. W.; Bagus, P. S. Mixed basis set calculations for atomic hydrogen on beryllium(0001). Chem. Phys. Lett. 1982, 90, 355. (7) Chesnut, D. B.; Moore, K. D. Locally dense basis sets for chemical shift calculations. J. Comput. Chem. 1989, 10, 648. (8) Orozco, M.; Luque, F. J. On the use of mixed basis sets to compute accurate molecular electrostatic potentials. Chem. Phys. Lett. 1989, 160, 305. (9) Hinton, J. F.; Guthrie, P.; Pulay, P.; Wolinski, K. Ab initio quantum mechanical calculation of the chemical shift anisotropy of the hydrogen atom in the (H2 O)17 water cluster. J. Am. Chem. Soc. 1992, 114, 1604. 41



(10) Kirby, R. A.; Hansen, A. E. Study of locally dense and locally saturated basis sets in localized molecular orbital calculations of nuclear shielding: Ab initio LORG calculations for

13

C and

17

O in norbornenone. Int. J. Quantum Chem. 1996, 57, 199.

(11) DiLabio, G. A.; Wright, J. S. Calculation of bond dissociation energies for large molecules using locally dense basis sets. Chem. Phys. Lett. 1998, 297, 181. (12) DiLabio, G. A. Using Locally Dense Basis Sets for the Determination of Molecular Properties. J. Phys. Chem. A 1999, 103, 11414. (13) Moon, S.; Case, D. A. A comparison of quantum chemical models for calculating NMR shielding parameters in peptides: Mixed basis set and ONIOM methods combined with a complete basis set extrapolation. J. Comput. Chem. 2006, 27, 825. (14) Reid, D. M.; Kobayashi, R.; Collins, M. A. Systematic Study of Locally Dense Basis Sets for NMR Shielding Constants. J. Chem. Theory Comput. 2014, 10, 146. (15) Lee, S. J. R.; Miyamoto, K.; Ding, F.; Manby, F. R.; Miller III, T. F. Density-based errors in mixed-basis mean-field electronic structure, with implications for embedding and QM/MM methods. Chem. Phys. Lett. 2017, 683, 375. (16) Havriliak, S.; King, H. F. Rydberg Radicals. 1. Frozen-Core Model for Rydberg Levels of the Ammonium Radical. J. Am. Chem. Soc. 1983, 105, 4. (17) Jurgens-Lutovsky, R.; Almlöf, J. Dual basis sets in calculations of electron correlation. Chem. Phys. Lett. 1991, 178, 451. (18) Wolinski, K.; Pulay, P. Second-order Møller–Plesset calculations with dual basis sets. J. Chem. Phys. 2003, 118, 9497. (19) Ksiazek, A.; Wolinski, K. Molecular properties with dual basis set methods. Mol. Phys. 2008, 106, 769.

42


Page 42 of 54

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


(20) Lee, M. S.; Head-Gordon, M. Polarized atomic orbitals for self-consistent field electronic structure calculations. J. Chem. Phys. 1997, 107, 9085. (21) Lee, M. S.; Head-Gordon, M. Absolute and relative energies from polarized atomic orbital self-consistent field calculations and a second order correction. Convergence with size and composition of the secondary basis. Comput. Chem. 2000, 24, 295. (22) Liang, W.; Head-Gordon, M. Approaching the Basis Set Limit in Density Functional Theory Calculations Using Dual Basis Sets without Diagonalization. J. Phys. Chem. A 2004, 108, 3206. (23) Steele, R. P.; DiStasio Jr., R. A.; Shao, Y.; Kong, J.; Head-Gordon, M. Dual-basis second-order Møller–Plesset perturbation theory: A reduced-cost reference for correlation calculations. J. Chem. Phys. 2006, 125, 074108. (24) Steele, R. P.; Shao, Y.; DiStasio, R. A.; Head-Gordon, M. Dual-Basis Analytic Gradients. 1. Self-Consistent Field Theory. J. Phys. Chem. A 2006, 110, 13915. (25) Steele, R. P.; Head-Gordon, M. Dual-basis self-consistent field methods: 6-31G* calculations with a minimal 6-4G primary basis. Mol. Phys. 2007, 105, 2455. (26) Distasio Jr., R. A.; Steele, R. P.; Head-Gordon, M. The analytical gradient of dualbasis resolution-of-the-identity second-order Møller–Plesset perturbation theory. Mol. Phys. 2007, 105, 2731. (27) Steele, R. P.; DiStasio Jr., R. A.; Head-Gordon, M. Non-Covalent Interactions with Dual-Basis Methods: Pairings for Augmented Basis Sets. J. Chem. Theory Comput. 2009, 5, 1560. (28) Steele, R. P.; Head-Gordon, M.; Tully, J. C. Ab Initio Molecular Dynamics with Dual Basis Set Methods. J. Phys. Chem. A 2010, 114, 11853.

43



(29) Mao, Y.; Horn, P. R.; Mardirossian, N.; Head-Gordon, T.; Skylaris, C.-K.; Head-Gordon, M. Approaching the basis set limit for DFT calculations using an environment-adapted minimal basis with perturbation theory: Formulation, proof of concept, and a pilot implementation. J. Chem. Phys. 2016, 145, 044109. (30) Marshall, M. S.; Sears, J. S.; Burns, L. A.; Brédas, J.-L.; Sherrill, C. D. An Error and Efficiency Analysis of Approximations to Møller–Plesset Perturbation Theory. J. Chem. Theory Comput. 2010, 6, 3681. (31) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. Approaching the Hartree–Fock limit by perturbative methods. J. Chem. Phys. 2009, 130, 231101. (32) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. Hartree–Fock perturbative corrections for total and reaction energies. J. Chem. Phys. 2010, 133, 044116. (33) Deng, J.; Gill, P. M. W. A new approach to dual-basis second-order Møller–Plesset calculations. J. Chem. Phys. 2011, 134, 081103. (34) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. MP2[V] – A Simple Approximation to Second-Order Møller–Plesset Perturbation Theory. J. Chem. Theory Comput. 2015, 11, 1639. (35) Nakajima, T.; Hirao, K. A dual-level approach to density-functional theory. J. Chem. Phys. 2006, 124, 184108. (36) Tokura, S.; Sato, T.; Tsuneda, T.; Nakajima, T.; Hirao, K. A dual-level state-specific time-dependent density-functional theory. J. Comput. Chem. 2008, 29, 1187. (37) Mizukami, W.; Nakajima, T.; Hirao, K.; Yanai, T. A dual-level approach to fourcomponent relativistic density-functional theory. Chem. Phys. Lett. 2011, 508, 177. (38) Deng, J.; Gilbert, A. T. B.; Gill, P. M. W. Density functional triple jumping. Phys. Chem. Chem. Phys. 2010, 12, 10759. 44


Page 44 of 54

Page 45 of 54 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) Kobayashi, M.; Nakai, H. Dual-level hierarchical scheme for linear-scaling divide-andconquer correlation theory. Int. J. Quantum Chem. 2009, 109, 2227. (40) Friedrich, J.; Dolg, M. Fully Automated Incremental Evaluation of MP2 and CCSD(T) Energies: Application to Water Clusters. J. Chem. Theory Comput. 2009, 5, 287. (41) Zhang, J.; Dolg, M. Approaching the complete basis set limit of CCSD(T) for large systems by the third-order incremental dual-basis set zero-buffer F12 method. J. Chem. Phys. 2014, 140, 044114. (42) Zhang, J.; Dolg, M. Third-Order Incremental Dual-Basis Set Zero-Buffer Approach for Large High-Spin Open-Shell Systems. J. Chem. Theory Comput. 2015, 11, 962. (43) Anacker, T.; Tew, D. P.; Friedrich, J. First UHF Implementation of the Incremental Scheme for Open-Shell Systems. J. Chem. Theory Comput. 2016, 12, 65. (44) Røeggen, I.; Gao, B. Combination of large and small basis sets in electronic structure calculations on large systems. J. Chem. Phys. 2018, 148, 134118. (45) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (46) Noga, J.; Kedžuch, S.; Šimunek, J. Second order explicitly correlated R12 theory revisited: A second quantization framework for treatment of the operators’ partitionings. J. Chem. Phys. 2007, 127, 034106. (47) Köhn, A.; Tew, D. P. Towards the Hartree–Fock and coupled-cluster singles and doubles basis set limit: A study of various models that employ single excitations into a complementary auxiliary basis set. J. Chem. Phys. 2010, 132, 024101. (48) Kong, L.; Valeev, E. F. Perturbative correction for the basis set incompleteness error of complete-active-space self-consistent field. J. Chem. Phys. 2010, 133, 174126.

45



(49) Warshel, A.; Karplus, M. Calculation of ground and excited state potential surfaces of conjugated molecules. I. Formulation and parametrization. J. Am. Chem. Soc. 1972, 94, 5612. (50) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227. (51) Maseras, F.; Morokuma, K. IMOMM – A new integrated ab-initio plus molecular mechanics geometry optimization scheme of equilibrium structures and transitionstates. J. Comput. Chem. 1995, 16, 1170. (52) Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. Angew. Chem. Int. Ed. 2009, 48, 1198. (53) Gomes, A. S. P.; Jacob, C. R. Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2012, 108, 222. (54) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation methods: A route to accurate calculations on large systems. Chem. Rev. 2012, 112, 632. (55) Libisch, F.; Huang, C.; Carter, E. A. Embedded correlated wavefunction schemes: Theory and applications. Acc. Chem. Res. 2014, 47, 2768. (56) Wesolowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891. (57) Knizia, G.; Chan, G. K.-L. Density Matrix Embedding: A Strong-Coupling Quantum Embedding Theory. J. Chem. Theory Comput. 2013, 9, 1428. (58) Cortona, P. Self-consistently determined properties of solids without band-structure calculations. Phys. Rev. B 1991, 44, 8454. 46


Page 46 of 54

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


(59) Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050. (60) Goodpaster, J. D.; Ananth, N.; Manby, F. R.; Miller III, T. F. Exact nonadditive kinetic potentials for embedded density functional theory. J. Chem. Phys. 2010, 133, 084103. (61) Elliott, P.; Cohen, M. H.; Wasserman, A.; Burke, K. Density Functional Partition Theory with Fractional Occupations. J. Chem. Theory Comput. 2009, 5, 827. (62) Elliott, P.; Burke, K.; Cohen, M. H.; Wasserman, A. Partition density-functional theory. Phys. Rev. A 2010, 82, 024501. (63) 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. (64) Huang, C.; Pavone, M.; Carter, E. A. Quantum mechanical embedding theory based on a unique embedding potential. J. Chem. Phys. 2011, 134, 154110. (65) Goodpaster, J. D.; Barnes, T. A.; Miller III, T. F. Embedded density functional theory for covalently bonded and strongly interacting subsystems. J. Chem. Phys. 2011, 134, 164108. (66) Nafziger, J.; Wu, Q.; Wasserman, A. Molecular binding energies from partition density functional theory. J. Chem. Phys. 2011, 135, 234101. (67) Rajchel, L.; Żuchowski, P. S.; Szcz¸eśniak, M. M.; Chałasiński, G. Derivation of the supermolecular interaction energy from the monomer densities in the density functional theory. Chem. Phys. Lett. 2010, 486, 160. (68) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller III, T. F. A Simple, Exact DensityFunctional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564. 47



(69) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller III, T. F. Accurate and systematically improvable density functional theory embedding for correlated wavefunctions. J. Chem. Phys. 2014, 140, 18A507. (70) Barnes, T. A.; Goodpaster, J. D.; Manby, F. R.; Miller III, T. F. Accurate basis set truncation for wavefunction embedding. J. Chem. Phys. 2013, 139, 024103. (71) Bennie, S. J.; Stella, M.; Miller III, T. F.; Manby, F. R. Accelerating wavefunction in density-functional-theory embedding by truncating the active basis set. J. Chem. Phys. 2015, 143, 024105. (72) 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. (73) Fornace, M. E.; Lee, J.; Miyamoto, K.; Manby, F. R.; Miller III, T. F. Embedded Mean-Field Theory. J. Chem. Theory Comput. 2015, 11, 568. (74) Ding, F.; Manby, F. R.; Miller III, T. F. Embedded Mean-Field Theory with BlockOrthogonalized Partitioning. J. Chem. Theory Comput. 2017, 13, 1605. (75) Li, S.; Shen, J.; Li, W.; Jiang, Y. An efficient implementation of the “cluster-inmolecule” approach for local electron correlation calculations. J. Chem. Phys. 2006, 125, 074109. (76) Mata, R. A.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106. (77) Li, W.; Piecuch, P. Multilevel Extension of the Cluster-in-Molecule Local Correlation Methodology: Merging Coupled-Cluster and Møller–Plesset Perturbation Theories. J. Phys. Chem. A 2010, 114, 6721. (78) Rolik, Z.; Kállay, M. A general-order local coupled-cluster method based on the clusterin-molecule approach. J. Chem. Phys. 2011, 135, 104111. 48


Page 48 of 54

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


(79) Myhre, R. H.; Sánchez de Merás, A. M. J.; Koch, H. Multi-level coupled cluster theory. J. Chem. Phys. 2014, 141, 224105. (80) 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. (81) Sparta, M.; Retegan, M.; Pinski, P.; Riplinger, C.; Becker, U.; Neese, F. Multilevel Approaches within the Local Pair Natural Orbital Framework. J. Chem. Theory Comput. 2017, 13, 3198. (82) Huzinaga, S.; Cantu, A. A. Theory of Separability of Many Electron Systems. J. Chem. Phys. 1971, 55, 5543. (83) Chulhai, D. V.; Goodpaster, J. D. Improved Accuracy and Efficiency in Quantum Embedding through Absolute Localization. J. Chem. Theory Comput. 2017, 13, 1503. (84) Culpitt, T.; Brorsen, K. R.; Hammes-Schiffer, S. Density functional theory embedding with the orthogonality constrained basis set expansion procedure. J. Chem. Phys. 2017, 146, 211101. (85) Miyamoto, K.; Miller III, T. F.; Manby, F. R. Fock-Matrix Corrections in Density Functional Theory and Use in Embedded Mean-Field Theory. J. Chem. Theory Comput. 2016, 12, 5811. (86) Boughton, J. W.; Pulay, P. Comparison of the Boys and Pipek–Mezey Localizations in the Local Correlation Approach and Automatic Virtual Basis Selection. J. Comput. Chem. 1993, 14, 736. (87) Rolik, Z.; Szegedy, L.; Ladjánszki, I.; Ladóczki, B.; Kállay, M. An efficient linearscaling CCSD(T) method based on local natural orbitals. J. Chem. Phys. 2013, 139, 094105. 49



(88) Kállay, M. Linear-scaling implementation of the direct random-phase approximation. J. Chem. Phys. 2015, 142, 204105. (89) Nagy, P. R.; Samu, G.; Kállay, M. An integral-direct linear-scaling second-order Møller–Plesset approach. J. Chem. Theory Comput. 2016, 12, 4897. (90) Nagy, P. R.; Kállay, M. Optimization of the linear-scaling local natural orbital CCSD(T) method: Redundancy-free triples correction using Laplace transform. J. Chem. Phys. 2017, 146, 214106. (91) Hégely, B.; Bogár, F.; Ferenczy, G. G.; Kállay, M. A QM/MM program for calculations with frozen localized orbitals based on the Huzinaga equation. Theor. Chem. Acc. 2015, 134, 132. (92) Mrcc, a quantum chemical program suite written by M. Kállay, Z. Rolik, J. Csontos, P. Nagy, G. Samu, D. Mester, J. Csóka, B. Szabó, I. Ladjánszki, L. Szegedy, B. Ladóczki, K. Petrov, M. Farkas, P. D. Mezei, and B. Hégely. See also Ref. 87 as well as http://www.mrcc.hu/ (accessed June 15, 2018). (93) Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta 1973, 28, 213. (94) 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. (95) Dunning Jr., 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. (96) Kendall, R. A.; Dunning Jr., T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796. 50


Page 50 of 54

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


(97) Woon, D. E.; Dunning Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358. (98) Weigend, F. Hartree–Fock Exchange Fitting Basis Sets for H to Rn. J. Comput. Chem. 2008, 29, 167. (99) Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116, 3175. (100) Dirac, P. A. M. Quantum mechanics of many-electron systems. Proc. Roy. Soc. (London) A 1929, 123, 714. (101) Slater, J. C. A simplification of the Hartree–Fock method. Phys. Rev. 1951, 81, 385. (102) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (103) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic-behavior. Phys. Rev. A 1988, 38, 3098. (104) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648. (105) Pipek, J.; Mezey, P. A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys. 1989, 90, 4916. (106) Bennie, S. J.; van der Kamp, M. W.; Pennifold, R. C. R.; Stella, M.; Manby, F. R.; Mulholland, A. J. A Projector-Embedding Approach for Multiscale Coupled-Cluster Calculations Applied to Citrate Synthase. J. Chem. Theory Comput. 2016, 12, 2689. (107) Zhang, X.; Bennie, S. J.; van der Kamp, M. W.; Glowacki, D. R.; Manby, F. R.; Mulholland, A. J. Multiscale analysis of enantioselectivity in enzyme-catalysed ‘lethal synthesis’ using projector-based embedding. Royal Soc. Open Sci 2018, 5, 171390. 51



(108) Kulik, H. J.; Zhang, J.; Klinman, J. P.; Martínez, T. J. How Large Should the QM Region Be in QM/MM Calculations? The Case of Catechol O-Methyltransferase. J. Phys. Chem. B 2016, 120, 11381. (109) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 2006, 25, 247. (110) Case, D. A.; Cerutti, D. S.; Cheatham III, T. E.; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Götz, A. W.; Greene, D.; Homeyer, N.; Izadi, S.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Liu, J.; Luchko, T.; Luo, R.; Mermelstein, D.; Merz, K. M.; Monard, G.; Nguyen, H.; Omelyan, I.; Onufriev, A.; Pan, F.; Qi, R.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; BotelloSmith, W. M.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Xiao, L.; York, D. M.; Kollman, P. A. Amber 2017. University of California, San Francisco, 2017. (111) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696. (112) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (113) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general Amber force field. J. Comput. Chem. 2004, 25, 1157. (114) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of highquality atomic charges. AM1-BCC model: I. Method. J. Comput. Chem. 2000, 21, 132.

52


Page 52 of 54

Page 53 of 54 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


(115) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem. 2002, 23, 1623. (116) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (117) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456. (118) Polly, R.; Werner, H.-J.; Manby, F. R.; Knowles, P. J. Fast Hartree–Fock theory using local fitting approximations. Mol. Phys. 2004, 102, 2311. (119) Foster, J. M.; Boys, S. F. Canonical Configurational Interaction Procedure. Rev. Mod. Phys. 1960, 32, 300.

53



Graphical TOC Entry

54


Page 54 of 54

Dual basis set approach for density functional and wave function

Recommend Documents