FULL PAPER
WWW.C-CHEM.ORG
Wave-Function Frozen-Density Embedding: Approximate Analytical Nuclear Ground-State Gradients €fener* Johannes Heuser and Sebastian Ho We report the derivation of approximate analytical nuclear ground-state uncoupled frozen density embedding (FDEu) gradients for the resolution of identity (RI) variant of the secondorder approximate coupled cluster singles and doubles (RICC2) as well as density functional theory (DFT), and an efficient implementation thereof in the KOALA program. In order to guarantee a computationally efficient treatment, those gradient terms are neglected which would require the exchange of orbital information. This approach allows for geometry optimizations of single molecules surrounded by numerous molecules with fixed nuclei at RICC2-in-RICC2, RICC2-in-DFT, and DFT-in-DFT FDE level of theory using a dispersion correction, required due to the DFT-based treatment of the interaction in FDE theory. Accuracy and applicability are assessed by the
example of two case studies: (a) the Watson-Crick pair adenine-thymine, for which the optimized structures exhibit a ˚ for our best scheme commaximum error of about 0.08 A pared to supermolecular reference calculations, (b) carbon monoxide on a magnesium oxide surface model, for which the ˚ for our best scheme. Efficiency is error amount up to 0.1 A demonstrated by successively including environment molecules and comparing to an optimized conventional supermolecular implementation, showing that the method is able to outperform conventional RICC2 schemes already with a rather small number of environment molecules, gaining significant C 2016 Wiley Periodicals, Inc. speed up in computation time. V
Introduction
wave-function methods while the treatment of the subsystem interaction (and its response) is based on density-functional theory (DFT).[20–23] For many applications the ground-state density and potential is sufficient, and the embedding of a wave-function based scheme in a DFT environment like RICC2in-DFT is adequate, in particular when solute molecules are surrounded by a solvent. Embedding a wave-function in a wave-function based environment treatment becomes important when similar subsystems are to be investigated and all subsystems thus need to be treated on the same footing. One of the most fundamental molecular properties is a molecule’s geometry and its associated potential-energy surface (PES), as it defines a molecule’s chemistry, that is, possible reaction pathways, transition states, and bond-breaking behavior, to name only a few. Therefore, the calculation of analytical nuclear gradients and geometry optimizations has become standard practice in computational chemistry to locate minima and transition states. In this context, low-scaling methods are of special interest, and among the most common ones are DFT schemes using the resolution of the identity (RI) and further approximations, see for example Ref. [24], as well as RIbased second-order Møller-Plesset perturbation theory (RIMP2), see for example Ref. [25]. The main purpose of RICC2 in the context of ground-state optimizations can be seen as a
Coupled-cluster methods are powerful tools as they allow for a systematic improvement of accuracy when solving the €dinger equation up to the full configuration interaction Schro (FCI) limit. However, this mightiness is limited by the scaling of the computation time with system size, in particular when explicit molecular interactions need to be taken into account. Thus only rather small molecular arrangements can be treated in practice with high accuracy, reducing the applicability range significantly. This is not only the case for ground-state energies, but in particular the calculation of molecular properties is even more limited due to increased prefactors, leading to a further decrease in the molecular size which can be investigated. A solution to the problem of a nonlinear scaling of computation time with increasing molecular size is given by reducedscaling methods that introduce small numerical approximations while decreasing the computation time significantly. A reduced scaling can be achieved in different ways and, alongside supermolecular (linear-scaling) methods, one ansatz that has become popular are so-called subsystem or embedding methods, which split the entire system into subsystem fragments.[1–4] One method for ab initio schemes is frozen-density embedding (FDE),[5–19] because it offers the possibility to combine any quantum-chemistry methods which offer electron densities since the interaction of the subsystems can be formulated in terms of electron densities only and large numbers of subsystems are easy to treat. In particular, for molecular response properties, FDE has been proven to be an efficient approach to divide a complex consisting of several molecules into subsystems, with all subsystems computed using ab initio 1092
Journal of Computational Chemistry 2016, 37, 1092–1101
DOI: 10.1002/jcc.24301
J. Heuser, S. H€ ofener Institut f€ ur Physikalische Chemie, Fakult€ at f€ ur Chemie und Biowissenschaften, Karlsruher Institut f€ ur Technologie (KIT), D-76131 Karlsruhe E-mail: [email protected] Contract grant sponsor: Deutsche Forschungsgemeinschaft DFG; Contract grant number: HO-4605/2-1 (SH) C 2016 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
needed precursor step for excited-state geometry optimizations. Analytical nuclear gradients have been explored in a pioneering work for DFT-in-DFT embedding,[26] and more implementations followed,[27,28] but to the best of our knowledge no analytical gradients for wave-function embedding have been presented yet. The purpose of this article is to fill the gap of missing FDE gradients for general computation setups, providing not only analytical nuclear gradients for DFT-in-DFT but also RICC2-inRICC2 as well as RICC2-in-DFT. The method reported in this work together with a geometry optimizer based on a levelshifted Newton procedure[29] is implemented in the program package KOALA.[21,30] The article is organized as follows. We start with a presentation of the method developed, followed by a discussion of the results of the new implementation. The article closes with a summary and a brief outlook.
Methodology This work uses three key features to ensure computational efficiency. 1. FDE techniques are used to reduce the overall scaling. The starting point for the FDE treatment is to partition the electron density q into N subsystem densities, N X q5 qi ;
(1)
i
where i runs over all subsystems. Then, the supermolecular Lagrangian L for the energy expression can be divided, for example, in case of two subsystems into three parts, L5LI 1Eint 1LII ;
(2)
2. The CC2 method is used to ensure a rather low scaling but wave-function based treatment. The conventional energy Lagrange functional for the CC2 method using singles (T1) and doubles (T2) excitations reads:[31] ^ LCC2 5hSCFjexp ð2T2 ÞHexp ðT2 ÞjSCFi X ^ ^ t l1 hl1 jH1½H; T2 jSCFi 1 l1
(5)
X ^ t l2 hl2 jH1½F; 1 T2 jSCFi l2
where H^ denotes the Hamilton operator transformed with the T1 operator, F is the Fock operator (vide infra), SCF is an abbreviation for self-consistent field, l1 indicates all single excitations, l2 indicates all double excitations, and t are the Lagrange multipliers for the cluster amplitudes. 3. The resolution of the identity is used to guarantee a scaling with a small prefactor. Using the Coulomb metric, one obtains for the integrals (written in ”chemists’ notation”):[32,33] X 1 1 jjkÞ ðlmj jjkÞRI 5 ðlmj GPlm GPjk ; (6) jr2r0 j jr2r0 j P X ðPj GPlm 5 Q
ðQjlmÞ5ðQj
1 jQÞ21=2 ðQjlmÞ jr2r0 j
(7)
1 jlmÞ ; jr2r0 j
(8)
21=2 1 denotes the matrix square root and P, Q where ðPj jr2r 0 j jQÞ are indices of an auxiliary basis and l; m; j; k belong to the conventional orbital basis. This allows to formulate efficient intermediates and avoids the need to store four-index quantities. In the following, these ingredients are used to derive approximate analytical nuclear gradients for the FDE RICC2 method to allow for efficient geometry optimizations of molecules in large molecular surroundings.
where Eint is the interaction contribution and LI and LII are the two subsystem energy functionals. For each subsystem arbitrary quantum-chemistry methods can be chosen, provided the electron density is available, while the interaction energy is expressed in terms of subsystem densities using DFT to capture also nonadditive (nadd) contributions[19]: ð ð II I I;II ðrÞdr1 qII ðrÞvnuc ðrÞdr 1Enuc Eint ½qI ; qII 5 qI ðrÞvnuc
Unrelaxed RICC2 FDE Ground-State Lagrangian
(3)
(11)
ðð 1
qI ðrÞqII ðr0 Þ nadd ½qI ; qII ; dr dr0 1Eemb:xck jr2r0 j
where “nuc” indicates potential of the nuclei of the subsystem I;II is the repulsion energy specified in the superscript index. Enuc of the nuclei of subsystems I and II. The nonadditive contribution is given as:
N X 1Ekin ½qI 1qII 2 Ekin ½qi : i
DRICC2 LRICC2 unrelax 5ESCF ½H1Lunrelax ½H1Eint 1LII ;
ESCF ½H5h SCF jHj SCF i ; X 1 X K;nsep ^ LDRICC2 DFpq Fpq 1 d ðpqjrsÞRI ; unrelax ½H5 2 pqrs pqrs pq
(4)
(9) (10)
where the subscript RI indicates that the two-electron integrals are approximated using the resolution of the identity,[34] and ðpq^jrsÞRI denotes T1-transformed two-electron integrals calculated using density fitting. The one-electron density DF reads: DFpq 5hSCFj½Epq ; T1 jSCFi X t l1 hl1 je2T1 Epq eT1 1½Epq ; T2 jSCFi 1
N X Exc ½qi
nadd Eemb:xck ½qI ; qII 5Exc ½qI 1qII 2
i
The conventional RICC2 Lagrangian can be found elsewhere.[31] For FDE, only the interaction energy and the environment contribution is added:
l1
(12)
X t l2 hl2 j½Epq ; T2 jSCFi ; 1 l2
Journal of Computational Chemistry 2016, 37, 1092–1101
1093
FULL PAPER
WWW.C-CHEM.ORG
and the nonseparable two-electron density is given as:
qðrÞCC 5hKje
2T
! X t l hlj e2T DðrÞeT jSCFi ; DðrÞe jSCFi5 hSCFj1 T
l
K;nsep 5hSCFj½½epqrs ; T1 ; T1 1½epqrs ; T2 jSCFi dpqrs
(21)
X t l1 hl1 j½epqrs ; T1 1½epqrs ; T2 jSCFi 1
(13)
l1
X t l2 hl2 jepqrs jSCFi ; 1 l2
where Epq and epqrs are the singlet one-electron and twoelectron excitation operators, respectively.[35] “[H]” denotes the dependence on the vacuum Hamilton operator only, that is, without embedding contributions. Note that, for the sake of clarity, all contributions without an explicit subsystem label belong to subsystem I and no working equations are presented for terms which arise also in conventional vacuum RICC2 gradients because they can be found in Ref. [34], for instance. It is important to note that due to the FDE ansatz, there is no explicit embedding-potential contribution in the Hamiltonian, for example, in eq. (10), because these contributions arise from the interaction energy expression. In analogy to the Hamilton operator, the Fock matrix F is defined without the embedding potential, X Fpq 5hpjhjqi1 ðpqj k
1 1X 1 ðpkj jkkÞ2 jqkÞ ; 0 jr2r j 2 k jr2r0 j
(16)
;
where “xck” indicates that an exchange-correlation and a kinetic contribution are included. To give an example, the embedding potential for subsystem I is calculated according to:
where DSCF pq denotes the SCF density. Compared to a conventional vacuum treatment, the Lagrangian terms are not zero since the amplitudes are optimized in the presence of the embedding potential. However, due to the stationary condition of the coupled-cluster amplitudes which include the embedding potential,[36] the energy thus reads emb emb 2LDRICC2 1Eint 1LII : (24) Eunrelax 5ESCF ½H1LDRICC2 unrelax ½H1v unrelax ½v emb , the Lagrangian contriNote that in case of LDRICC2 unrelax ½H1v butions vanish and a conventional energy expression can be emb used, but not in case of LDRICC2 : unrelax ½v
X X hpq DSCF Gpq ½DSCF DSCF ESCF ½H5 pq 1 pq 1Enuc ; pq
emb DRICC2 5Eunrelax ½H1v emb ; LDRICC2 unrelax ½H1v emb 5 LDRICC2 unrelax ½v
X
(25)
pq
emb DFpq vpq :
(26) (27)
pq
ð q ðr0 Þ II vIemb:coul ðrÞ5vnuc 1 II 0 dr0 jr2r j
(17)
vIemb:xck ðrÞ5vIemb:xc ðrÞ1vIemb:kin ðrÞ
(18) (19) I
dEkin ½q dEkin ½q 2 (20) vIemb:kin ðrÞ5 dq q5qI 1qII dq q5q I
Using coupled-cluster methods, the interaction energy is obtained from (unrelaxed) coupled-cluster densities, Journal of Computational Chemistry 2016, 37, 1092–1101
(23)
(15)
where “½H1vemb ” denotes the dependence on both the vacuum Hamilton operator as well as the embedding potential vemb obtained from the first derivative of eq. (3) with respect to the density of the active subsystem. The embedding potential can be split into two parts:
1094
X F emb:coul jqi1Eemb:xck ½q; qenv ; Eint 5 ðDSCF pq 1Dpq Þhpjv
(22)
(14)
emb dpq ep ½H1v emb 5Fpq 1vpq ;
dExc ½q dExc ½q vIemb:xc ðrÞ5 2 dq q5qI 1qII dq q5q
DRICC2 LRICC2 unrelax 5ESCF ½H1Lunrelax ½H1Eint 1LII ;
pq
where p, q denote general molecular orbitals (MOs), k denote occupied MOs, and h the sum of the one-electron operators. Therefore, the following relation between Fock matrix and orbital energies is obtained:
emb emb:coul emb:xck vpq 5vpq 1vpq
where DðrÞ is the density operator at point r. Note that, as a consequence, the energy functional L does not depend linearly on the Lagrangian multipliers due to the exchange-correlation and kinetic contributions. Environment effects in the subsystem coupled-cluster calculation are incorporated as effective embedding potentials, which are obtained from an expansion of the interaction energy, eq. (3).[36] Although the embedding contributions are written similarly to a conventional quantum operator, for instance, the energy cannot be calculated as an expectation value of the potential representation times the electron density. Therefore, the embedding contribution is thus denoted embedding potential (instead of embedding operator) analogously to the exchange-correlation potential in DFT. The full Lagrangian can thus explicitly be formulated as:
The exact treatment of the Lagrangian constraints and multipliers is discussed in prior publications.[36] The entire frozendensity embedding procedure can be summarized as follows. In a first step, the energies and densities of all molecules are calculated in vacuum, since all densities are required to calculate the embedding potential of any subsystem. In the next step, an embedding potential for the target subsystem is constructed and, since the embedding potential contains nonlinear contributions also from the “active” density, the embedding potential for the target subsystem is updated until convergence, denoted macroiterations. Conventional iterations to obtain a density or wave-function, needed for solving the WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Kohn-Sham or coupled-cluster equations, for instance, are referred to as micro iterations here. After this procedure, the electronic structure of the target molecule is relaxed with respect to the environment, but the environment molecules have not yet been relaxed to any surrounding. Subsequently, the target molecule is frozen and one environment molecule after the other is thawed, updated using macroiterations, and frozen again. In the last step, the target molecule is thawed and updated again using macroiterations, required due to the relaxed environment. One “super-iteration” of which the last step is always the calculation of the target system, is denoted freeze–thaw iteration. These freeze–thaw iterations are carried out until convergence or a predefined number of cycles is reached, followed by a single calculation of the molecular properties of the target system, for example, in this work, the orbital relaxation multipliers and the gradient itself. Up to this point, the full supermolecular Lagrangian including the environment energy LII has been taken into account. In general, all three contributions, LI ; Eint ; LII , depend explicitly or implicitly on the electron density of the active subsystem and will change upon a geometrical distortion, for instance. In order to ensure computational efficiency, we neglect orbitaldependent contributions of the environment subsystems in the gradient calculation of the active subsystem, which is often denoted “uncoupled” FDE. Applying this approximate uncoupled scheme, efficient gradients are available independently of the quantum-chemistry methods used for the subsystems. We therefore neglect the contributions of LII in the following for LRICC2 relax . A comparison to numerical FDE geometry optimizations to investigate the errors introduced by the proposed approximations is contained in the Results and Discussion section.
Relaxed RICC2 FDE Ground-State Lagrangian In order to obtain analytical nuclear gradients, orbital-relaxed densities are needed to account for the orbital change at the perturbed geometry. The SCF determinant is therefore expressed as j SCF i5exp
X
! †
jrs ðErs 2Esr Þ j SCF 0 i ;
of all dependencies would lead to additional relaxation cycles and since the coupling elements are expected to be small, all of these contributions are neglected in the following. The approximate Lagrangian used in this work is explicitly given as: DRICC2 LRICC2 relax 5ESCF ½H1Lunrelax ½H1Eint X X jI? FI? j 1vIemb aI FIa 1vIaemb ; 1 j j 1 ?j I? j
(30)
Ia
X F emb:coul jqi1Eemb:xck ½q; qenv ; Eint 5 ðDSCF pq 1Dpq Þhpjv
(31)
pq
where the indices I? ; J? refer to frozen occupied orbitals, i, j refer to nonfrozen occupied orbitals, I, J refer to both frozen and active occupied orbitals, a, b refer to virtual orbitals, and p, q refer to all occupied and virtual orbitals. In this work, the occupied–occupied Lagrangian multipliers need to be calculated explicitly due to additional terms arising in the embedding treatment and cannot be calculated on-thefly in contrast to standard vacuum implementations.[34] The orbital relaxation coefficients are determined by differentiating with respect to orbital rotations. Note that since the Hamilton operator is defined without an embedding potential but, for example, the SCF iterations are solved in presence of the embedding potential, both contributions are not zero by themselves but cancel each other when added. This leaves the task to solve: !
05
@ DRICC2 @ X @ X jI? FI? j 1vIemb aI FIa 1vIaemb : Lunrelax ½H1 j j 1 ?j @jIJ @jIJ I? j @jIJ Ia
(32) Since the occupied–occupied rotations do not depend on the occupied–virtual rotations, one set of equations is solved for j, gLHS 52gRHS ;
(33)
with a standard Davidson solver[37–39] with the left-hand side (LHS):
(28)
X emb:xck KL dIK dJL ðeJ 2eI Þ12fKLIJ gLHS j ; IJ 5
rs
(34)
KL
and j SCF 0 i is chosen such that for the unperturbed system jrs 5 0. This leads to the introduction of orbital relaxation Lagrangian multipliers for the active subsystem only due to the neglect of the contributions arising from LII. For FDE, this yields: X DRICC2 l0 Fl0 1vlemb j 1Eint : LRICC2 relax 5ESCF ½H1Lunrelax ½H1 0
(29)
l0
where j denote Lagrangian multipliers for orbital rotations and the dimension of l0 depends on the choice of the constraint. A rigorous derivation for coupled cluster amplitudes and multipliers would therefore include also couplings with the orbital relaxation Lagrangian multipliers. Since a rigorous accounting
and a conventional RICC2 right-hand side (RHS)[34]: X
gRHS IJ 5
ðxÞ
ðxÞ
ðxÞ
ðxÞ
fHIJ 2XIJ 1HJI 2XJI g :
(35)
x5a;b;c
femb:xck denotes kernel contributions which arise as the second functional derivative of the exchange-correlation and kinetic interaction energy. The nonseparable contribution is given as: X
ðxÞ ðxÞ ðxÞ fHmn 2XðxÞ mn 1Hnm 2Xnm g
x5a;b;c
i Xh nsep nsep nsep nsep 51 ðdpmrs 1dmprs ÞðpnjrsÞRI 2ðdpnrs 1dnprs ÞðmpjrsÞRI
(36)
prs
Journal of Computational Chemistry 2016, 37, 1092–1101
1095
FULL PAPER
WWW.C-CHEM.ORG
X
nsep dpqrs 5
kppp0 khqq0 kprr0 khss0 dpK;nsep 0 q0 r 0 s 0
(37)
;
p0 q0 r0 s0 p
!
@ DRICC2 @ X @ X jI? FI? j 1vIemb aI FIa 1vIaemb ; j j Lunrelax ½H1 1 ?j @jbJ @jbJ I? j @jbJ Ia
1
X
emb:xck;½x
Dj;ao lm vlm
lm
vac-RICC2 is the convenwhere [x] denotes the partial derivative, Lrelax [34] tional vacuum RICC2 Lagrangian, S is the overlap matrix, Dj is the matrix of the Lagrange multipliers for orbital rotations, and a superscript ao denotes backtransformed densities:
X ao 5 Xpq Clp Cmq : Xlm
(39)
(44)
pq
which can be solved using a standard Davidson solver with
Ia
X eff;emb;ao ½x ½DSCF 1DF 2 Flm Slm ;
lm
(38)
X X emb:xck Ia dab dIJ ðea 2eI Þ1AaIbJ 12faIbJ j ; gLHS aI 5
(43)
lm
h
with the transformation matrices k 512t1 and k 511t1 using the single excitation t1 amplitudes.[34] In case of the occupied-virtual orbital relaxation multipliers, again the contributions from SCF and the interaction energy cancel each other by virtue of eq. (15). This leads to the equation: 05
dLRICC2 @Lvac-RICC2 X eff;ao emb:coul;½x relax Dlm vlm 5 relax 1 dx @x lm X ½x emb:xck 1 ðDSCF;ao 1DF;ao jmi lm lm Þ2hl jv
The effective density is given here as:
bJ
SCF F j pq 1 j qp Þ=2 : Deff pq 5Dpq 1Dpq 1ð
(45)
and 1 F 1X F F ðD 1DFqp ÞApqAI gRHS aI 5 ðDAI 1DIA ÞðeI 2eA Þ2 2 4 pq pq X emb:xck 1 DjKL ðAKLAI 12fKLAI Þ
One should highlight the integral derivative in the last term in eq. (43), which is given for LDA-type functionals as
(40)
@Eemb:xck ðlðrÞmðrÞÞ½x dr @qa r ! ðX @ 2 Eemb:xck @ 2 Eemb:xck Djk 1 1 @q2a r @qa @qb r jk
emb:xck;½x vlm ½D5
KL
1
X
ðxÞ
ðxÞ
ðxÞ
ðxÞ
fHAI 2XAI 1HIA 2XIA g :
x5a;b;c
(46)
3lðrÞmðrÞðjðrÞkðrÞÞ½x dr ;
The matrix A, Apqrs 54ðpqjrsÞ2ðprjsqÞ2ðpsjrqÞ ;
(41)
arises from the differentiation of the Fock matrix elements with respect to the orbital rotations,[40] ^ 2j X @F½ej He pq 2 ^ 5 hSCFj½a1 pr ; ½aqr ; ½Ers ; H1 jSCFi @jrs x5x0 r
where qa, qb denote spin densities, since it corresponds to kernel contributions which only arise in a wave-function based embedding treatment due to the orbital rotation multipliers for gradients, while in conventional methods these terms arise in the calculation of analytical second nuclear derivatives. Such a kernel contribution does not only enter the gradient itself, but also the effective Fock matrix, X eff;emb emb:xck Fpq 5 Djrs frspq ;
5dqr Fsp 2dps Fqr 2dqs Frp 1dpr Fqs 1ðmr 2ms ÞApqrs ; (42) where m are Hartree-Fock occupation numbers of orbitals r and s, respectively.
Analytical Nuclear Gradients for FDE CoupledCluster Once all wave-function parameters such as orbitals, cluster amplitudes, and Lagrangian parameters have been determined, the nuclear gradients can be calculated. In quantumchemistry methods, an elegant way to express any analytic nuclear gradient is to formulate it in terms of effective densities which are contracted with differentiated integrals. Note that the gradient expression for the active subsystem is identical for RICC2-in-RICC2 and RICC2-in-DFT using the present work’s scheme. Using RICC2 embedded in a (DFT or coupledcluster) FDE environment, this leads to the total gradient expression: 1096
ð
Journal of Computational Chemistry 2016, 37, 1092–1101
(47)
rs
which is contracted with differentiated overlap integrals. All other contributions are conventional RICC2 gradient contributions.[34]
Analytical Nuclear Gradients for FDE DFT For the sake of completeness, we also report the equations for analytical nuclear gradients for DFT-in-DFT embedding for an implementation using RI and Gaussian basis functions: dLKS @Lvac-KS 5 dx @x X emb:coul;½x ½x emb:xck DSCF;ao v 12hl jv jmi ; 1 lm lm
(48)
lm
where[x] denotes the partial derivative and Lvac-KS then conventional vacuum Kohn-Sham Lagrangian.[41,42] It can be seen that WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
changing the role of “active” and “environment” molecules successively, this was not carried out for this study. Therefore, all geometry optimizations using embedding refer to the geometry optimization of one subsystem, while the other subsystem’s geometry was kept fixed.
The adenine–thymine dimer Figure 1. The adenine–thymine dimer in its Watson–Crick pairing.
eq. (48) is significantly simplified compared to the RICC2 case in eq. (43). The most important fact is, however, that no kernel contributions have to be calculated due to the absence of orbital rotation multipliers.
Results and Discussion Computational details Reference supermolecular calculations have been carried out using Turbomole 6.6[43] with default convergence criteria for both energy calculations as well as geometry optimizations. All embedding calculations have been carried out using the KOALA program using PBE as the exchange-correlation functional and PW91k as the kinetic energy functional.[21,30] Dispersion corrections are included using Grimme’s D3 program.[44,45] The number of times the environment subsystems were relaxed is denoted “Freeze–Thaw [m],” with “Freeze–Thaw [0]” indicating that the environment was not relaxed and “Freeze– Thaw [1]” that it has been relaxed once.[21] The geometry optimizer consists of a level-shifted quasi Newton procedure,[29] using a Broyden–Fletcher–Shanno–Goldfarb (BFGS) update for the inverse Hessian.[46–50] Geometry optimizations were considered converged when all of the following criteria were fulfilled: Change of (a) energy
WWW.C-CHEM.ORG
Wave-Function Frozen-Density Embedding: Approximate Analytical Nuclear Ground-State Gradients €fener* Johannes Heuser and Sebastian Ho We report the derivation of approximate analytical nuclear ground-state uncoupled frozen density embedding (FDEu) gradients for the resolution of identity (RI) variant of the secondorder approximate coupled cluster singles and doubles (RICC2) as well as density functional theory (DFT), and an efficient implementation thereof in the KOALA program. In order to guarantee a computationally efficient treatment, those gradient terms are neglected which would require the exchange of orbital information. This approach allows for geometry optimizations of single molecules surrounded by numerous molecules with fixed nuclei at RICC2-in-RICC2, RICC2-in-DFT, and DFT-in-DFT FDE level of theory using a dispersion correction, required due to the DFT-based treatment of the interaction in FDE theory. Accuracy and applicability are assessed by the
example of two case studies: (a) the Watson-Crick pair adenine-thymine, for which the optimized structures exhibit a ˚ for our best scheme commaximum error of about 0.08 A pared to supermolecular reference calculations, (b) carbon monoxide on a magnesium oxide surface model, for which the ˚ for our best scheme. Efficiency is error amount up to 0.1 A demonstrated by successively including environment molecules and comparing to an optimized conventional supermolecular implementation, showing that the method is able to outperform conventional RICC2 schemes already with a rather small number of environment molecules, gaining significant C 2016 Wiley Periodicals, Inc. speed up in computation time. V
Introduction
wave-function methods while the treatment of the subsystem interaction (and its response) is based on density-functional theory (DFT).[20–23] For many applications the ground-state density and potential is sufficient, and the embedding of a wave-function based scheme in a DFT environment like RICC2in-DFT is adequate, in particular when solute molecules are surrounded by a solvent. Embedding a wave-function in a wave-function based environment treatment becomes important when similar subsystems are to be investigated and all subsystems thus need to be treated on the same footing. One of the most fundamental molecular properties is a molecule’s geometry and its associated potential-energy surface (PES), as it defines a molecule’s chemistry, that is, possible reaction pathways, transition states, and bond-breaking behavior, to name only a few. Therefore, the calculation of analytical nuclear gradients and geometry optimizations has become standard practice in computational chemistry to locate minima and transition states. In this context, low-scaling methods are of special interest, and among the most common ones are DFT schemes using the resolution of the identity (RI) and further approximations, see for example Ref. [24], as well as RIbased second-order Møller-Plesset perturbation theory (RIMP2), see for example Ref. [25]. The main purpose of RICC2 in the context of ground-state optimizations can be seen as a
Coupled-cluster methods are powerful tools as they allow for a systematic improvement of accuracy when solving the €dinger equation up to the full configuration interaction Schro (FCI) limit. However, this mightiness is limited by the scaling of the computation time with system size, in particular when explicit molecular interactions need to be taken into account. Thus only rather small molecular arrangements can be treated in practice with high accuracy, reducing the applicability range significantly. This is not only the case for ground-state energies, but in particular the calculation of molecular properties is even more limited due to increased prefactors, leading to a further decrease in the molecular size which can be investigated. A solution to the problem of a nonlinear scaling of computation time with increasing molecular size is given by reducedscaling methods that introduce small numerical approximations while decreasing the computation time significantly. A reduced scaling can be achieved in different ways and, alongside supermolecular (linear-scaling) methods, one ansatz that has become popular are so-called subsystem or embedding methods, which split the entire system into subsystem fragments.[1–4] One method for ab initio schemes is frozen-density embedding (FDE),[5–19] because it offers the possibility to combine any quantum-chemistry methods which offer electron densities since the interaction of the subsystems can be formulated in terms of electron densities only and large numbers of subsystems are easy to treat. In particular, for molecular response properties, FDE has been proven to be an efficient approach to divide a complex consisting of several molecules into subsystems, with all subsystems computed using ab initio 1092
Journal of Computational Chemistry 2016, 37, 1092–1101
DOI: 10.1002/jcc.24301
J. Heuser, S. H€ ofener Institut f€ ur Physikalische Chemie, Fakult€ at f€ ur Chemie und Biowissenschaften, Karlsruher Institut f€ ur Technologie (KIT), D-76131 Karlsruhe E-mail: [email protected] Contract grant sponsor: Deutsche Forschungsgemeinschaft DFG; Contract grant number: HO-4605/2-1 (SH) C 2016 Wiley Periodicals, Inc. V
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
needed precursor step for excited-state geometry optimizations. Analytical nuclear gradients have been explored in a pioneering work for DFT-in-DFT embedding,[26] and more implementations followed,[27,28] but to the best of our knowledge no analytical gradients for wave-function embedding have been presented yet. The purpose of this article is to fill the gap of missing FDE gradients for general computation setups, providing not only analytical nuclear gradients for DFT-in-DFT but also RICC2-inRICC2 as well as RICC2-in-DFT. The method reported in this work together with a geometry optimizer based on a levelshifted Newton procedure[29] is implemented in the program package KOALA.[21,30] The article is organized as follows. We start with a presentation of the method developed, followed by a discussion of the results of the new implementation. The article closes with a summary and a brief outlook.
Methodology This work uses three key features to ensure computational efficiency. 1. FDE techniques are used to reduce the overall scaling. The starting point for the FDE treatment is to partition the electron density q into N subsystem densities, N X q5 qi ;
(1)
i
where i runs over all subsystems. Then, the supermolecular Lagrangian L for the energy expression can be divided, for example, in case of two subsystems into three parts, L5LI 1Eint 1LII ;
(2)
2. The CC2 method is used to ensure a rather low scaling but wave-function based treatment. The conventional energy Lagrange functional for the CC2 method using singles (T1) and doubles (T2) excitations reads:[31] ^ LCC2 5hSCFjexp ð2T2 ÞHexp ðT2 ÞjSCFi X ^ ^ t l1 hl1 jH1½H; T2 jSCFi 1 l1
(5)
X ^ t l2 hl2 jH1½F; 1 T2 jSCFi l2
where H^ denotes the Hamilton operator transformed with the T1 operator, F is the Fock operator (vide infra), SCF is an abbreviation for self-consistent field, l1 indicates all single excitations, l2 indicates all double excitations, and t are the Lagrange multipliers for the cluster amplitudes. 3. The resolution of the identity is used to guarantee a scaling with a small prefactor. Using the Coulomb metric, one obtains for the integrals (written in ”chemists’ notation”):[32,33] X 1 1 jjkÞ ðlmj jjkÞRI 5 ðlmj GPlm GPjk ; (6) jr2r0 j jr2r0 j P X ðPj GPlm 5 Q
ðQjlmÞ5ðQj
1 jQÞ21=2 ðQjlmÞ jr2r0 j
(7)
1 jlmÞ ; jr2r0 j
(8)
21=2 1 denotes the matrix square root and P, Q where ðPj jr2r 0 j jQÞ are indices of an auxiliary basis and l; m; j; k belong to the conventional orbital basis. This allows to formulate efficient intermediates and avoids the need to store four-index quantities. In the following, these ingredients are used to derive approximate analytical nuclear gradients for the FDE RICC2 method to allow for efficient geometry optimizations of molecules in large molecular surroundings.
where Eint is the interaction contribution and LI and LII are the two subsystem energy functionals. For each subsystem arbitrary quantum-chemistry methods can be chosen, provided the electron density is available, while the interaction energy is expressed in terms of subsystem densities using DFT to capture also nonadditive (nadd) contributions[19]: ð ð II I I;II ðrÞdr1 qII ðrÞvnuc ðrÞdr 1Enuc Eint ½qI ; qII 5 qI ðrÞvnuc
Unrelaxed RICC2 FDE Ground-State Lagrangian
(3)
(11)
ðð 1
qI ðrÞqII ðr0 Þ nadd ½qI ; qII ; dr dr0 1Eemb:xck jr2r0 j
where “nuc” indicates potential of the nuclei of the subsystem I;II is the repulsion energy specified in the superscript index. Enuc of the nuclei of subsystems I and II. The nonadditive contribution is given as:
N X 1Ekin ½qI 1qII 2 Ekin ½qi : i
DRICC2 LRICC2 unrelax 5ESCF ½H1Lunrelax ½H1Eint 1LII ;
ESCF ½H5h SCF jHj SCF i ; X 1 X K;nsep ^ LDRICC2 DFpq Fpq 1 d ðpqjrsÞRI ; unrelax ½H5 2 pqrs pqrs pq
(4)
(9) (10)
where the subscript RI indicates that the two-electron integrals are approximated using the resolution of the identity,[34] and ðpq^jrsÞRI denotes T1-transformed two-electron integrals calculated using density fitting. The one-electron density DF reads: DFpq 5hSCFj½Epq ; T1 jSCFi X t l1 hl1 je2T1 Epq eT1 1½Epq ; T2 jSCFi 1
N X Exc ½qi
nadd Eemb:xck ½qI ; qII 5Exc ½qI 1qII 2
i
The conventional RICC2 Lagrangian can be found elsewhere.[31] For FDE, only the interaction energy and the environment contribution is added:
l1
(12)
X t l2 hl2 j½Epq ; T2 jSCFi ; 1 l2
Journal of Computational Chemistry 2016, 37, 1092–1101
1093
FULL PAPER
WWW.C-CHEM.ORG
and the nonseparable two-electron density is given as:
qðrÞCC 5hKje
2T
! X t l hlj e2T DðrÞeT jSCFi ; DðrÞe jSCFi5 hSCFj1 T
l
K;nsep 5hSCFj½½epqrs ; T1 ; T1 1½epqrs ; T2 jSCFi dpqrs
(21)
X t l1 hl1 j½epqrs ; T1 1½epqrs ; T2 jSCFi 1
(13)
l1
X t l2 hl2 jepqrs jSCFi ; 1 l2
where Epq and epqrs are the singlet one-electron and twoelectron excitation operators, respectively.[35] “[H]” denotes the dependence on the vacuum Hamilton operator only, that is, without embedding contributions. Note that, for the sake of clarity, all contributions without an explicit subsystem label belong to subsystem I and no working equations are presented for terms which arise also in conventional vacuum RICC2 gradients because they can be found in Ref. [34], for instance. It is important to note that due to the FDE ansatz, there is no explicit embedding-potential contribution in the Hamiltonian, for example, in eq. (10), because these contributions arise from the interaction energy expression. In analogy to the Hamilton operator, the Fock matrix F is defined without the embedding potential, X Fpq 5hpjhjqi1 ðpqj k
1 1X 1 ðpkj jkkÞ2 jqkÞ ; 0 jr2r j 2 k jr2r0 j
(16)
;
where “xck” indicates that an exchange-correlation and a kinetic contribution are included. To give an example, the embedding potential for subsystem I is calculated according to:
where DSCF pq denotes the SCF density. Compared to a conventional vacuum treatment, the Lagrangian terms are not zero since the amplitudes are optimized in the presence of the embedding potential. However, due to the stationary condition of the coupled-cluster amplitudes which include the embedding potential,[36] the energy thus reads emb emb 2LDRICC2 1Eint 1LII : (24) Eunrelax 5ESCF ½H1LDRICC2 unrelax ½H1v unrelax ½v emb , the Lagrangian contriNote that in case of LDRICC2 unrelax ½H1v butions vanish and a conventional energy expression can be emb used, but not in case of LDRICC2 : unrelax ½v
X X hpq DSCF Gpq ½DSCF DSCF ESCF ½H5 pq 1 pq 1Enuc ; pq
emb DRICC2 5Eunrelax ½H1v emb ; LDRICC2 unrelax ½H1v emb 5 LDRICC2 unrelax ½v
X
(25)
pq
emb DFpq vpq :
(26) (27)
pq
ð q ðr0 Þ II vIemb:coul ðrÞ5vnuc 1 II 0 dr0 jr2r j
(17)
vIemb:xck ðrÞ5vIemb:xc ðrÞ1vIemb:kin ðrÞ
(18) (19) I
dEkin ½q dEkin ½q 2 (20) vIemb:kin ðrÞ5 dq q5qI 1qII dq q5q I
Using coupled-cluster methods, the interaction energy is obtained from (unrelaxed) coupled-cluster densities, Journal of Computational Chemistry 2016, 37, 1092–1101
(23)
(15)
where “½H1vemb ” denotes the dependence on both the vacuum Hamilton operator as well as the embedding potential vemb obtained from the first derivative of eq. (3) with respect to the density of the active subsystem. The embedding potential can be split into two parts:
1094
X F emb:coul jqi1Eemb:xck ½q; qenv ; Eint 5 ðDSCF pq 1Dpq Þhpjv
(22)
(14)
emb dpq ep ½H1v emb 5Fpq 1vpq ;
dExc ½q dExc ½q vIemb:xc ðrÞ5 2 dq q5qI 1qII dq q5q
DRICC2 LRICC2 unrelax 5ESCF ½H1Lunrelax ½H1Eint 1LII ;
pq
where p, q denote general molecular orbitals (MOs), k denote occupied MOs, and h the sum of the one-electron operators. Therefore, the following relation between Fock matrix and orbital energies is obtained:
emb emb:coul emb:xck vpq 5vpq 1vpq
where DðrÞ is the density operator at point r. Note that, as a consequence, the energy functional L does not depend linearly on the Lagrangian multipliers due to the exchange-correlation and kinetic contributions. Environment effects in the subsystem coupled-cluster calculation are incorporated as effective embedding potentials, which are obtained from an expansion of the interaction energy, eq. (3).[36] Although the embedding contributions are written similarly to a conventional quantum operator, for instance, the energy cannot be calculated as an expectation value of the potential representation times the electron density. Therefore, the embedding contribution is thus denoted embedding potential (instead of embedding operator) analogously to the exchange-correlation potential in DFT. The full Lagrangian can thus explicitly be formulated as:
The exact treatment of the Lagrangian constraints and multipliers is discussed in prior publications.[36] The entire frozendensity embedding procedure can be summarized as follows. In a first step, the energies and densities of all molecules are calculated in vacuum, since all densities are required to calculate the embedding potential of any subsystem. In the next step, an embedding potential for the target subsystem is constructed and, since the embedding potential contains nonlinear contributions also from the “active” density, the embedding potential for the target subsystem is updated until convergence, denoted macroiterations. Conventional iterations to obtain a density or wave-function, needed for solving the WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
Kohn-Sham or coupled-cluster equations, for instance, are referred to as micro iterations here. After this procedure, the electronic structure of the target molecule is relaxed with respect to the environment, but the environment molecules have not yet been relaxed to any surrounding. Subsequently, the target molecule is frozen and one environment molecule after the other is thawed, updated using macroiterations, and frozen again. In the last step, the target molecule is thawed and updated again using macroiterations, required due to the relaxed environment. One “super-iteration” of which the last step is always the calculation of the target system, is denoted freeze–thaw iteration. These freeze–thaw iterations are carried out until convergence or a predefined number of cycles is reached, followed by a single calculation of the molecular properties of the target system, for example, in this work, the orbital relaxation multipliers and the gradient itself. Up to this point, the full supermolecular Lagrangian including the environment energy LII has been taken into account. In general, all three contributions, LI ; Eint ; LII , depend explicitly or implicitly on the electron density of the active subsystem and will change upon a geometrical distortion, for instance. In order to ensure computational efficiency, we neglect orbitaldependent contributions of the environment subsystems in the gradient calculation of the active subsystem, which is often denoted “uncoupled” FDE. Applying this approximate uncoupled scheme, efficient gradients are available independently of the quantum-chemistry methods used for the subsystems. We therefore neglect the contributions of LII in the following for LRICC2 relax . A comparison to numerical FDE geometry optimizations to investigate the errors introduced by the proposed approximations is contained in the Results and Discussion section.
Relaxed RICC2 FDE Ground-State Lagrangian In order to obtain analytical nuclear gradients, orbital-relaxed densities are needed to account for the orbital change at the perturbed geometry. The SCF determinant is therefore expressed as j SCF i5exp
X
! †
jrs ðErs 2Esr Þ j SCF 0 i ;
of all dependencies would lead to additional relaxation cycles and since the coupling elements are expected to be small, all of these contributions are neglected in the following. The approximate Lagrangian used in this work is explicitly given as: DRICC2 LRICC2 relax 5ESCF ½H1Lunrelax ½H1Eint X X jI? FI? j 1vIemb aI FIa 1vIaemb ; 1 j j 1 ?j I? j
(30)
Ia
X F emb:coul jqi1Eemb:xck ½q; qenv ; Eint 5 ðDSCF pq 1Dpq Þhpjv
(31)
pq
where the indices I? ; J? refer to frozen occupied orbitals, i, j refer to nonfrozen occupied orbitals, I, J refer to both frozen and active occupied orbitals, a, b refer to virtual orbitals, and p, q refer to all occupied and virtual orbitals. In this work, the occupied–occupied Lagrangian multipliers need to be calculated explicitly due to additional terms arising in the embedding treatment and cannot be calculated on-thefly in contrast to standard vacuum implementations.[34] The orbital relaxation coefficients are determined by differentiating with respect to orbital rotations. Note that since the Hamilton operator is defined without an embedding potential but, for example, the SCF iterations are solved in presence of the embedding potential, both contributions are not zero by themselves but cancel each other when added. This leaves the task to solve: !
05
@ DRICC2 @ X @ X jI? FI? j 1vIemb aI FIa 1vIaemb : Lunrelax ½H1 j j 1 ?j @jIJ @jIJ I? j @jIJ Ia
(32) Since the occupied–occupied rotations do not depend on the occupied–virtual rotations, one set of equations is solved for j, gLHS 52gRHS ;
(33)
with a standard Davidson solver[37–39] with the left-hand side (LHS):
(28)
X emb:xck KL dIK dJL ðeJ 2eI Þ12fKLIJ gLHS j ; IJ 5
rs
(34)
KL
and j SCF 0 i is chosen such that for the unperturbed system jrs 5 0. This leads to the introduction of orbital relaxation Lagrangian multipliers for the active subsystem only due to the neglect of the contributions arising from LII. For FDE, this yields: X DRICC2 l0 Fl0 1vlemb j 1Eint : LRICC2 relax 5ESCF ½H1Lunrelax ½H1 0
(29)
l0
where j denote Lagrangian multipliers for orbital rotations and the dimension of l0 depends on the choice of the constraint. A rigorous derivation for coupled cluster amplitudes and multipliers would therefore include also couplings with the orbital relaxation Lagrangian multipliers. Since a rigorous accounting
and a conventional RICC2 right-hand side (RHS)[34]: X
gRHS IJ 5
ðxÞ
ðxÞ
ðxÞ
ðxÞ
fHIJ 2XIJ 1HJI 2XJI g :
(35)
x5a;b;c
femb:xck denotes kernel contributions which arise as the second functional derivative of the exchange-correlation and kinetic interaction energy. The nonseparable contribution is given as: X
ðxÞ ðxÞ ðxÞ fHmn 2XðxÞ mn 1Hnm 2Xnm g
x5a;b;c
i Xh nsep nsep nsep nsep 51 ðdpmrs 1dmprs ÞðpnjrsÞRI 2ðdpnrs 1dnprs ÞðmpjrsÞRI
(36)
prs
Journal of Computational Chemistry 2016, 37, 1092–1101
1095
FULL PAPER
WWW.C-CHEM.ORG
X
nsep dpqrs 5
kppp0 khqq0 kprr0 khss0 dpK;nsep 0 q0 r 0 s 0
(37)
;
p0 q0 r0 s0 p
!
@ DRICC2 @ X @ X jI? FI? j 1vIemb aI FIa 1vIaemb ; j j Lunrelax ½H1 1 ?j @jbJ @jbJ I? j @jbJ Ia
1
X
emb:xck;½x
Dj;ao lm vlm
lm
vac-RICC2 is the convenwhere [x] denotes the partial derivative, Lrelax [34] tional vacuum RICC2 Lagrangian, S is the overlap matrix, Dj is the matrix of the Lagrange multipliers for orbital rotations, and a superscript ao denotes backtransformed densities:
X ao 5 Xpq Clp Cmq : Xlm
(39)
(44)
pq
which can be solved using a standard Davidson solver with
Ia
X eff;emb;ao ½x ½DSCF 1DF 2 Flm Slm ;
lm
(38)
X X emb:xck Ia dab dIJ ðea 2eI Þ1AaIbJ 12faIbJ j ; gLHS aI 5
(43)
lm
h
with the transformation matrices k 512t1 and k 511t1 using the single excitation t1 amplitudes.[34] In case of the occupied-virtual orbital relaxation multipliers, again the contributions from SCF and the interaction energy cancel each other by virtue of eq. (15). This leads to the equation: 05
dLRICC2 @Lvac-RICC2 X eff;ao emb:coul;½x relax Dlm vlm 5 relax 1 dx @x lm X ½x emb:xck 1 ðDSCF;ao 1DF;ao jmi lm lm Þ2hl jv
The effective density is given here as:
bJ
SCF F j pq 1 j qp Þ=2 : Deff pq 5Dpq 1Dpq 1ð
(45)
and 1 F 1X F F ðD 1DFqp ÞApqAI gRHS aI 5 ðDAI 1DIA ÞðeI 2eA Þ2 2 4 pq pq X emb:xck 1 DjKL ðAKLAI 12fKLAI Þ
One should highlight the integral derivative in the last term in eq. (43), which is given for LDA-type functionals as
(40)
@Eemb:xck ðlðrÞmðrÞÞ½x dr @qa r ! ðX @ 2 Eemb:xck @ 2 Eemb:xck Djk 1 1 @q2a r @qa @qb r jk
emb:xck;½x vlm ½D5
KL
1
X
ðxÞ
ðxÞ
ðxÞ
ðxÞ
fHAI 2XAI 1HIA 2XIA g :
x5a;b;c
(46)
3lðrÞmðrÞðjðrÞkðrÞÞ½x dr ;
The matrix A, Apqrs 54ðpqjrsÞ2ðprjsqÞ2ðpsjrqÞ ;
(41)
arises from the differentiation of the Fock matrix elements with respect to the orbital rotations,[40] ^ 2j X @F½ej He pq 2 ^ 5 hSCFj½a1 pr ; ½aqr ; ½Ers ; H1 jSCFi @jrs x5x0 r
where qa, qb denote spin densities, since it corresponds to kernel contributions which only arise in a wave-function based embedding treatment due to the orbital rotation multipliers for gradients, while in conventional methods these terms arise in the calculation of analytical second nuclear derivatives. Such a kernel contribution does not only enter the gradient itself, but also the effective Fock matrix, X eff;emb emb:xck Fpq 5 Djrs frspq ;
5dqr Fsp 2dps Fqr 2dqs Frp 1dpr Fqs 1ðmr 2ms ÞApqrs ; (42) where m are Hartree-Fock occupation numbers of orbitals r and s, respectively.
Analytical Nuclear Gradients for FDE CoupledCluster Once all wave-function parameters such as orbitals, cluster amplitudes, and Lagrangian parameters have been determined, the nuclear gradients can be calculated. In quantumchemistry methods, an elegant way to express any analytic nuclear gradient is to formulate it in terms of effective densities which are contracted with differentiated integrals. Note that the gradient expression for the active subsystem is identical for RICC2-in-RICC2 and RICC2-in-DFT using the present work’s scheme. Using RICC2 embedded in a (DFT or coupledcluster) FDE environment, this leads to the total gradient expression: 1096
ð
Journal of Computational Chemistry 2016, 37, 1092–1101
(47)
rs
which is contracted with differentiated overlap integrals. All other contributions are conventional RICC2 gradient contributions.[34]
Analytical Nuclear Gradients for FDE DFT For the sake of completeness, we also report the equations for analytical nuclear gradients for DFT-in-DFT embedding for an implementation using RI and Gaussian basis functions: dLKS @Lvac-KS 5 dx @x X emb:coul;½x ½x emb:xck DSCF;ao v 12hl jv jmi ; 1 lm lm
(48)
lm
where[x] denotes the partial derivative and Lvac-KS then conventional vacuum Kohn-Sham Lagrangian.[41,42] It can be seen that WWW.CHEMISTRYVIEWS.COM
FULL PAPER
WWW.C-CHEM.ORG
changing the role of “active” and “environment” molecules successively, this was not carried out for this study. Therefore, all geometry optimizations using embedding refer to the geometry optimization of one subsystem, while the other subsystem’s geometry was kept fixed.
The adenine–thymine dimer Figure 1. The adenine–thymine dimer in its Watson–Crick pairing.
eq. (48) is significantly simplified compared to the RICC2 case in eq. (43). The most important fact is, however, that no kernel contributions have to be calculated due to the absence of orbital rotation multipliers.
Results and Discussion Computational details Reference supermolecular calculations have been carried out using Turbomole 6.6[43] with default convergence criteria for both energy calculations as well as geometry optimizations. All embedding calculations have been carried out using the KOALA program using PBE as the exchange-correlation functional and PW91k as the kinetic energy functional.[21,30] Dispersion corrections are included using Grimme’s D3 program.[44,45] The number of times the environment subsystems were relaxed is denoted “Freeze–Thaw [m],” with “Freeze–Thaw [0]” indicating that the environment was not relaxed and “Freeze– Thaw [1]” that it has been relaxed once.[21] The geometry optimizer consists of a level-shifted quasi Newton procedure,[29] using a Broyden–Fletcher–Shanno–Goldfarb (BFGS) update for the inverse Hessian.[46–50] Geometry optimizations were considered converged when all of the following criteria were fulfilled: Change of (a) energy
