Analytical Nuclear Excited-State Gradients for the Second-Order

Subscriber access provided by MT ROYAL COLLEGE

Quantum Electronic Structure

Analytical nuclear excited-state gradients for the secondorder approximate coupled-cluster singles and doubles (CC2) method employing uncoupled frozen-density embedding Johannes Heuser, and Sebastian Hoefener J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00369 • Publication Date (Web): 07 Aug 2018 Downloaded from http://pubs.acs.org on August 12, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

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

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

Journal of Chemical Theory and Computation

Analytical nuclear excited-state gradients for the second-order approximate coupled-cluster singles and doubles (CC2) method employing uncoupled frozen-density embedding Johannes Heuser and Sebastian Höfener∗ Institute of Physical Chemistry, Karlsruhe Institute of Technology (KIT), P.O. Box 6980, D-76049 Karlsruhe, Germany E-mail: [email protected] Phone: +49 721 608-43321. Fax: +49 721 608-47225 Abstract We report the derivation and implementation of analytical orbital-relaxed properties and nuclear gradients for excited states using the second-order approximate coupled-cluster singles and doubles (CC2) model combined with uncoupled frozendensity embedding (FDEu). An implementation of the algebraic diagrammatic construction through second order ADC(2), which arises from simplification of RICC2 FDEu, is also presented. In order to ensure a RICC2 FDEu Lagrange functional that is linear in the Lagrange multipliers the Hartree-Fock density is employed for the target subsystem in the embedding contributions. The accuracy of the new scheme is assessed using the carbon monoxide molecule, 4-aminophthalimide, and a benzonitrile dimer, revealing that the obtained errors are below the method error of RICC2. Using density functional theory (DFT) for the environment, the efficiency of the new method is

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

illustrated by computing the perturbed excited-state dipole moment of a chromophore in a biological environment. For this system, comprising 32 molecules consisting of 366 atoms in total, the computation requires only a couple of days on a standard compute node. RICC2 FDEu thus enables large-scale calculations of ab-initio wave functions for molecules in complex environments as routine applications.

1

Introduction

The accurate description of potential energy surfaces of excited states is still a major challenge for quantum chemistry. While a large variety of problems, ranging from spectroscopy to photochemistry or from measurements in vacuum to biological applications, is taking place on excited states’ potential energy surfaces, current state-of-the art quantum chemical methods are often not applicable to the complex environments of interest. Density functional theory for the electronic ground state often provides a good costaccuracy ratio and can thus be applied to large systems nowadays in particular when efficient approximations such as the resolution-of-the-identity and multipole-accelerated techniques are employed. 1 However, for excited states DFT-based methods exhibit an increased scaling compared to the ground-state, thus restricting the potential applicability to medium-sized systems. One possible solution to the problem is to employ subsystem schemes which partition the entire supermolecule into smaller subsystems. Among various schemes, frozendensity embedding (FDE) has proven to be a valuable tool for a significant reduction in computation time. 2–9 This is achieved by dividing a system into subunits, allowing also for a convenient subsystem analysis and chemical interpretation. 10 For DFT FDE, analytical ground-state gradients 11–15 as well as excited-state gradients have been reported recently. 16,17 While being rather accurate for molecular properties in the electronic ground state, timedependent DFT (TDDFT) methods can yield a rather poor description for excited states, in particular for charge-transfer excitations. The inclusion of exact Hartree-Fock exchange reduces the numerical errors to some extend but in many cases unsystematic errors remain. 2


Page 2 of 41

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


In order to overcome this problem, a number of range-separated functionals were developed, for example the Coulomb-attenuation method (CAM), 18 which can lead to more robust descriptions of charge-transfer excitations for some molecules. 19,20 However, it does not provide accurate results in all cases and, for example, the CAMB3LYP method is not able to provide accurate results for rhodopsin models. 21 The skillful inclusion of exact Hartree-Fock exchange can thus reduce the errors in selected cases but in general a careful analysis of the outcome of TDDFT calculations is necessary. Ab-initio methods such as coupled cluster do not suffer from increased errors for chargetransfer states 22 and are thus often able to cover a broader range of excited states exhibiting a single-reference character. With respect to excited state gradient theory, pioneering work has been carried out by Stanton and Gauss for the coupled-cluster singles and doubles (CCSD) model which exhibits sixth order scaling computational costs. 23–26 To reduce the scaling, in 1995 the second-order approximate coupled-cluster singles and doubles (CC2) was reported which exhibits fifth order scaling computational costs while the energies are correct to second order. 27 The prefactor of the scaling can be reduced significantly when density fitting (DF) and a strictly doubles-direct formulation is used, leading to applicability to systems with more than 1000 basis functions routinely. 28 Due to the efficiency, for this method excited-state gradients and properties can be computed for a wide range of mediumsized molecules. 29–32 However, even for very efficient implementations the fifth order scaling computational costs of conventional RICC2 prevent the application to molecules in complex environments and additional approximations such as the Laplace-transformation are invoked to reduce the scaling. 33 In the present work, we combine the ab-initio RICC2 method with FDE to compute analytical excited-state properties including nuclear gradients, enabling us to study locally excited states in complex environments. The work is organized as follows. We first provide the derivation analytical nuclear gradients of an efficient variant of the RICC2-FDE for excited states, followed by the assessment of the new method using four test cases. The

3



Page 4 of 41

efficiency is discussed using a chromophore in a biological environment.

2

Theory

In this section, the derivation of analytical nuclear gradients is presented for the excitedstate uncoupled FDE (FDEu) 5,34 Lagrangian as implemented in the KOALA program. In FDEu no orbital information is shared among the subsystems. The equations as reported are exact if the environment densities are not relaxed with respect to their surrounding as illustrated using a numerical example in the results section. The method thus corresponds to approximate gradients in case of one or more freeze-thaw (FT) iterations, but only small numerical deviations are introduced as also illustrated using a numerical example in the results section.

2.1

Frozen-Density Embedding

In FDE, the total electron density is constructed as the sum of the subsystem electron densities:

ρ(r) =

X

ρz (r) .

(1)

z

Each of those subsystem densities ρz is obtained from subsystem calculations for which any desired electronic structure method can be used provided the electron density is available. From this partitioning, the energy can then be expressed as the sum of the subsystem energies and an interaction contribution. For two subsystems I and II the total energy reads:

Etot [ρ(r)] = Eint [ρI (r), ρII (r)] + EI [ρI (r)] + EII [ρII (r)] ,

4


(2)

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


where the interaction between the subsystems is given as: 35

Eint [ρI , ρII ] =

Z +

II ρI (r)vnuc (r)dr

+

Z

I I,II ρII (r)vnuc (r)dr + Enuc

ρI (r)ρII (r′ ) nadd dr dr′ + Eemb:xck [ρI , ρII ] . |r − r′ |

Z Z

(3)

Independent of the methods used for the subsystems, the non-additive exchange-correlation and the kinetic energy contributions are always expressed using density-functional theory (DFT) and read: nadd [ρI , ρII ] = Exc [ρI + ρII ] − Eemb:xck

X

Exc [ρz ] + Ts [ρI + ρII ] −

z=I,II

X

Ts [ρz ] .

(4)

z=I,II

Finally, the embedding potential is defined as the functional derivative of the interaction energy, e.g., for subsystem I while keeping ρII fixed: 36 I vemb (r)

2.2

=

δEint [ρ] . δρ(r) ρ=ρI

(5)

The FDE procedure for coupled-cluster methods

Based on the partitioning in Eq. (1), for FDE coupled-cluster methods in general the corresponding (unrelaxed) coupled-cluster density is used in embedding calculations: CC Etot [ρ(r)] = Eint [ρCC I (r), ρII (r)] + ECC [ρI (r)] + EII [ρII (r)] .

(6)

The coupled-cluster density is obtained from projections of the Schrödinger equation which can be combined in one (orbital-relaxed ground-state) coupled-cluster Lagrange functional,

LCC,grd = hHF| e

−TN

He

TN

|HFi +

N X X i=1

5

t¯µi Ωµi +

µi


X µ0

, κ ¯ µ0 Fµvac 0

(7)


where Ωµi are the ground-state amplitude equations, t¯µi are the ground-state Lagrange multipliers, the cluster operator T contains the coupled cluster amplitudes, and N denotes the truncation level, e.g., 2 in case of coupled cluster with singles and doubles amplitudes. 32 In the present work, µ0 is used as a general index pair (pq), µ1 denote single excitations (ai), µ2 denote double excitations (ai, bj). Fµvac denotes the Fock matrix and κ ¯ µ0 are the Lagrangian 0 multipliers for orbital rotations. In the coupled cluster case, the FDE procedure is as follows. First the electron density of all environment subsystems is calculated without any FDE contributions. In the next step, the usual Hartree-Fock and coupled-cluster iterations are carried out to obtain orbital energies and cluster amplitudes, denoted micro iterations. Then the unrelaxed CC2 density has to be constructed for which the Lagrangian multipliers have to be obtained. In the last step the embedding potential is constructed and the interaction energy is calculated. This procedure is repeated until convergence, denoted macro iterations. These macro iterations are subsequently carried out for all environment subsystems. Finally, the target subsystem is relaxed using macro iterations with respect to the environment densities. One cycle to update all molecules once is called one FT cycle. After convergence or a given number of FT cycles, excitation energies or molecular properties such as gradients can be calculated for the target subsystem. Whereas in supermolecular approaches the response of the entire system is included, in embedding schemes the excited-state polarization of the environment arises as a separate perturbation additional to the ground-state polarization. One possible ansatz to account for this contribution is given by state-specific embedding approaches. 37–41 While state-specific FDE would increase the complexity significantly as it requires macro iterations over excitation energies and excited-state Lagrangian multipliers to construct the target-state unrelaxed excited-state CC2 density, followed by freeze-thaw iterations to polarize the environment with the excited-state contributions, the present article is concerned with efficient approximations to reduce the computation time.

6


Page 6 of 41

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


Instead of employing the coupled-cluster density, an efficient approximation is to use the Hartree-Fock density for the embedding treatment in a perturbation to the energy (PTE) manner for the target subsystem I, 42 HF Etot [ρ(r)] = Eint [ρHF I (r), ρII (r)] + ECC [ρI (r)] + E[ρII (r)] .

(8)

If only the Hartree-Fock density is computed, no cluster amplitudes or Lagrangian multipliers have to be calculated during the macro iterations and the FT cycles, resulting in a significant speedup. Furthermore, this approximate treatment provides a Lagrangian that is linear in the coupled cluster multipliers, since no coupled-cluster amplitudes or multipliers arise in the non-additive contributions in the interaction energy: ∂Eint [ρHF ∂Eint [ρHF I , ρII ] I , ρII ] = =0 ∂t ∂¯t

(9)

∂ρHF (r) ∂ρHF I (r) = I¯ = 0. ∂t ∂t

(10)

since

The accuracy of this method in particular for analytical nuclear gradients is discussed in the following.

2.3

The excited-state CC2 FDEu Lagrange function

The conventional orbital-relaxed CC2 Lagrange function of an excited state f for the vacuum Hamiltonian H is given as: 27,30 ˆ + [H, ˆ T2 ] |HFi + LfCC2,exc,vac = hHF| H

2 X X i=1

+¯ ωf

1−

2 X X i=1

E¯µfi Eµfi

µi

7

!

+

t¯fµi Ωµi +

µi

X

2 X X

E¯µfi Aµi ,νj Eνfj

i,j=1 µi ,νi

κ ¯ fµ0 Fvac µ0

µ0


(11)


Page 8 of 41

¯ and E are the left and right excitation vectors, where ω ¯ f denotes the excitation energy, E respectively, A is the CC2 Jacobian matrix, 28

Aµν =

∂Ωµ , ∂tν

(12)

ˆ = e−T1 HeT1 . The Lagrange multipliers now and T1 -transformed quantities are defined as H include a superscript

f

because the excited-state contributions are now included compared

to Eq. (7). For our purpose it is more convenient to rewrite the Lagrangian using effective one-particle and two-particle densities: 29 LfCC2,exc,vac = hHF| H |HFi +

X pq

vac F,ex,f κ Dpq + Dpq Fpq

1 X ˆnsep,ξ,f ˆnsep,A,f ˆ (pq|rs)DF d + dpqrs + 2 pqrs pqrs

F,ex,f F,ξ,f F,A,f Dpq = Dpq + Dpq ,

(13) (14)

where it is assumed that the left and right eigenvectors form a biorthonormal set. The usual indices for the orbital subsets are used, i.e. (p, q, . . . ) for general orbitals, (a, b, . . . ) to indicate virtual orbitals, (I, J, . . . ) for all occupied orbitals, (I ⋆ , J ⋆ , . . . ) for frozen occupied orbitals and (i, j, . . . ) to indicate active occupied orbitals. The definition of the effective densities DF,ξ,f and DF,A,f is given in Ref. 29, for example. Dκ is a matrix containing the Lagrange multipliers for orbital rotations. For FDE, only a few additional terms need to be included compared to the vacuum Lagrangian. For example, the subsystem Fock matrices contain the embedding potential also, F = Fvac + vemb ,

(15)

so that only the matrix representation of the Fock operator including the embedding potential

8


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


is diagonal: emb εp δpq = Fvac pq + vpq .

(16)

In principle, the interaction energy Eint as well as the environment Lagrangian LII are added to Eq. (14). However, as the inclusion of LII is problematic if the derivative of the Lagrangian with respect to orbital rotations of the active subsystem has to be calculated since it requires information about orbitals to be shared between the subsystems. We therefore neglect the contributions of LII in the spirit of FDEu. 14,17 The FDEu CC2 Lagrange functional can then be written as: LfCC2,exc,FDEu = hHF| H |HFi +

X pq

F,ex,f κ Dpq + Dpq Fpq

1 X ˆnsep,ξ,f ˆnsep,A,f ˆ + (pq|rs)DF + Eint . d + dpqrs 2 pqrs pqrs

2.4

(17)

Lagrange multipliers

For excited states, the Lagrange multipliers ¯tf consist of a ground-state and an excited-state contribution, ¯f t¯fµ = t¯(0) µ + Nµ ,

(18)

which arise from the individual terms when enforcing stationarity of the Lagrangian, 30,43 X X ∂LfCC2,exc ! = ην(0) + t¯fµ Aµν + E¯µf Bµγν Eγf = 0 , ∂tν µ µγ

9


(19)


Page 10 of 41

with ˆ + [H, ˆ T2 ] |HFi ∂ hHF| H , ∂tν ∂Aµν . = ∂tγ

ην(0) = Bµνγ

(20) (21)

The excited-state contribution of the multiplier response can be separated, yielding the ¯ f : 30 response equation for the excited-state Lagrangian multipliers N X

¯µf Aµν = − N

µ

X

E¯µf Bµγν Eγf .

(22)

µγ

In the case of FDEu and using the Hartree-Fock density for the embedding contribution, FDE contributions enter via orbital energies only, so that no further derivative contributions have to be implemented for this equation, cf. Eq. (10).

2.5

Orbital rotations

In order to ensure stationarity of the Lagrangian with respect to orbital rotations, Eq. (17) is differentiated as follows: X ∂LfCC2,exc ∂Fµ0 X F,ex,f ∂Fpq = + κ ¯ fµ0 Dpq ∂κν ∂κν ∂κν µ pq 0

ˆ 1 X ˆnsep,ξ,f ˆnsep,A,f ∂ (pq|rs) DF ! + = 0, dpqrs + dpqrs 2 pqrs ∂κν

(23)

Rearranging Eq. (23) and using the general pair index pq, cf. Eq. (7), yields the Z-vector equation, X pq

κ ¯ fpq

∂Fpq RHS = −ηrs ∂κrs

vac ∂Frs = ∂κpq

(24)

vac vac vac vac δqr Fsp − δps Fqr − δqs Frp + δpr Fqs + (nr − ns )A(κ) pqrs , 10


(25)

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


from which the orbital rotation multipliers κ ¯ are obtained. ηνRHS is the right-hand side (RHS) which is defined by comparison to Eq. (23), np are the occupation numbers, and the supermatrix A(κ) reads: A(κ) pqrs = 4 (pq|rs) − (pr|sq) − (ps|rq) .

(26)

Differentiating the embedding potential results in embedding kernel contributions, 34 I wemb:xck (r′ , r′′ )

=

nadd nadd δEemb:xck [ρ] δEemb:xck [ρ] − , δρ(r′ )δρ(r′′ ) ρ=ρI +ρII δρ(r′ )δρ(r′′ ) ρ=ρI

emb:xck I wpqrs = (pq|wemb:xck |rs) .

(27) (28)

The frozen–occupied orbital rotation multipliers are obtained in a first step since the virtual–occupied orbital rotation multipliers do not contribute. For frozen–occupied rotations, the following Z-vector equation is obtained: X

emb:xck κ ¯ fkL⋆ δik δJ ⋆ L⋆ (ǫi − ǫJ ⋆ ) + 2wiJ ⋆ kL⋆

kL⋆

RHS = ηiJ ⋆ ,

(29)

where RHS ηiJ = ⋆

ex

HiJ ⋆ − ex ΩiJ ⋆ − ex HJ ⋆ i + ex ΩJ ⋆ i .

(30)

Here, the frozen–occupied block of DF,ex,f is zero. The intermediates ex H and ex Ω are defined as: 29 ex

Hmn − ex Ωmn − ex Hnm + ex Ωnm o X n ˆ ˆ ˆnsep,ex,f + dˆnsep,ex,f (mp|rs) ˆnsep,ex,f (pn|rs) . (31) d − dˆnsep,ex,f + d = nprs pnrs pmrs mprs DF DF prs

Details how to obtain an efficient implementation of these intermediates is given in Ref. 29.

11



Page 12 of 41

The Z-vector equation for the virtual–occupied rotations reads: X JB

(κ) RHS emb:xck = ηIA , κ ¯ fJB δAB δIJ (ǫB − ǫJ ) + AIAJB + 2wIAJB

(32)

for which the right-hand side is given as: RHS ηIA =

2.6

(κ) 1 X F,ex F,ex F,ex emb:xck F,ex + DAI DIA (ǫI − ǫA ) − ApqIA + 2wpqIA Dpq + Dpq 2 pq 1 X κ (κ) emb:xck − + (ex HAI − ex ΩAI − ex HIA + ex ΩIA ) . (33) DkL⋆ AkL⋆ IA + 2wkL ⋆ IA 2 kL⋆

Analytical nuclear gradients

Once the molecular orbital (MO) coefficients, coupled cluster amplitudes, multipliers and excitation vectors, as well as the Lagrange multipliers for the orbital rotations are computed, the gradient can be obtained as the partial derivative of the total Lagrangian, Eq. (17), with ~ respect to displacement of the nuclear coordinates of the active subsystem X: dLCC2,exc,FDEu ∂LCC2,exc X eff,emb [X] ~ = − Fpq Spq ~ ~ dX ∂X pq Z ∂ I,II ∂ I + + ρII (r) E v (r) dr ~ nuc ~ nuc ∂X ∂X X X ~ ~ HF I HF [X] I + Dpq hp|vemb:coul |qi[X] + 2 Dpq hp |vemb:xck |qi pq

+

X

pq

I F,ex κ |qi (Dpq + Dpq )hp|vemb

~ [X]

.

(34)

pq

The analytical nuclear excited-state gradients for conventional CC2,

∂LfCC2,exc , ~ ∂X

is given in

Ref. 29. The embedding contributions to the effective Fock matrix to account for the metric change upon geometric perturbations of the underlying Gaussian basis functions are as

12


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


follows: ← Feff,emb Ip

X

(35)

F,ex F,ex κ emb:xck (Drs + Dsr + Drs ) wrsIp .

rs

2.7

Related wavefunction models

Starting from a CC2 implementation, the excitation energies and properties for the methods CIS(D∞ ) 44 and the (strict) algebraic diagrammatic construction through second order, ADC(2), 45 can be obtained using simple recipes. 31 The CIS(D∞ ) model uses the ground-state second-order Møller-Plesset perturbation theory (MP2) amplitudes in combination with the CIS(D∞ ) Jacobian, 31 

 ACIS(D∞ ) = 

hai |H

−

EHF |ck i

+

h|[[H, T2MP2 ], τkc ]|i

c hab ij |H|k i

hai |H|cd kl i cd hab ij |(F − E0 )|kl i



 ,

(36)

where EHF = hHF|H|HFi and E0 = hHF|F |HFi. The secular matrix for (strict) ADC(2) is obtained from the CIS(D∞ ) Jacobi matrix by symmetrization: 46 AADC(2) =

1 ACIS(D∞ ) + (ACIS(D∞ ) )† . 2

(37)

For canonical orbitals only the singles–singles block needs to be symmetrized because all other blocks are already Hermitian. The implementation of the methods CIS(D∞ ) and (strict) ADC(2) for excitation energies, response properties, and analytical nuclear gradients is thus straightforward if a CC2 implementation is available. This holds not only for the vacuum case but also in particular for the FDEu method when employing the Hartree-Fock density.

13



3

Computational details

All FDE calculations were carried out using the KOALA program. 47,48 The embedding density of the target subsystem, ρI , was constructed from the Hartree-Fock density if not stated otherwise. The embedding density of the environment subsystems, ρII , was obtained using the density corresponding to the method. For example, in case of RIMP2 the unrelaxed RIMP2 density was used to construct the embedding potential. The embedding exchangecorrelation functional used was PBE, 49 the embedding kinetic energy functional used was PW91k. 50 To account for dispersion effects between the subsystems, Grimme’s D3 correction is employed. 51–53 The KOALA program uses grids based on the TURBOMOLE 54,55 grids. 56 For FDE geometry optimizations, grid 2 was used. For subsequent single-point calculations to obtain excitation energies as well as dipole moments, grid 5 was used. Geometry optimizations were considered converged when all of the following criteria were fulfilled: Change of (a) energy < 1 · 10−6 , (b) RMS displacement < 5 · 10−4 , (c) RMS gradient < 5·10−4 , (d) max. displacement < 1·10−3 and (e) max. gradient < 1·10−3 . The coordinates of the environment subsystems were kept fixed during the optimization. Supermolecule calculations were carried out with TURBOMOLE. 54,55 For better comparison, the coordinates of the atoms that belong to the environment subsystems during the FDE runs were also kept frozen in the supermolecule optimization. Unless stated otherwise, the environment density was not relaxed with respect to the other subsystems. Numerical gradients were calculated by displacing the carbon atom by ±0.0005 Bohr along the three Cartesian coordinates at a given structure. Vibrational corrections are computed at RICC2/def2-TZVP level of theory using semi-numerical second derivatives (obtained from analytical gradients). The same vibrational corrections are used for all methods in the present work.

14


Page 14 of 41

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


4

Results

In this section, results as obtained using the new implementation are presented. First, the use of the Hartree-Fock density is assessed and the accuracy compared to numerical gradients is investigated. Second, the implementation is used to investigate local excited states in two molecular complexes, for which excitation and de-excitation energies are studied. Third, the new method is used to study the influence of the environment density upon the computation of orbital-relaxed excited-state dipole moments of photosynthetic chlorophyll clusters. Finally, the new implementation is used to compute vertical excitation energies and dipole moments in a biological rhodopsin model.

4.1

Assessing the Hartree-Fock density approximation

In order to assess the accuracy of the new approach and to verify our implementation, the lowest singlet excited state (S1 ) of carbon monoxide (CO) in the presence of one water molecule is investigated, see Fig. 1. For the CO molecule, it is known that the HartreeFock method fails to predict the correct sign of the (ground-state) dipole moment while the correlated density provides the correct sign of the dipole moment. Therefore, this molecule seems to be a well-suited candidate to assess the accuracy of the presented method because the environment molecule is polarized with an electron density that is shifted in the wrong direction. However, in the following it is shown that the wrong sign of the (ground-state)

Figure 1: The geometry of the CO–H2 O complex as used to assess the accuracy of the new method. Oxygen atoms are red, carbon atoms grey and hydrogen atoms white. 15



Table 1: Fluorescence energies obtained as vertical de-excitation energies (VDE) from the lowest singlet excited state (S1 ) to the ground-state in the CO–H2 O complex at a C–H distance of 2.53 Å in eV. The supermolecular RICC2 result is 7.31 eV employing the same geometry.

b

FTa HFb CC2c 0 7.331 7.331 1 7.332 7.334 2 7.332 7.334 3 7.332 7.334 a Number of freeze-thaw iterations. Hartree-Fock density employed for FDE. c CC2 density employed for FDE.

dipole moment has only a small influence on the excitation energies. For this system, the entire complex was optimized using RICC2/def2-TZVPP in the ground state, followed by a re-optimization of the CO bond distance in the lowest singlet excited state using RICC2–in–RICC2/def2-TZVPPD and grid 2 with a frozen water geometry. In the following, the subsystem geometries are kept fixed and only the intermolecular

Figure 2: Interaction energies of the lowest singlet state in the CO–H2 O complex, cf. Fig. 1, using the Hartree-Fock (HF) and CC2 density in milli-Hartree, computed as the difference between the total energy obtained from a FDE calculation and the sum of the ground-state energy of the water molecule and the lowest singlet excited state of the CO molecule.

16


Page 16 of 41

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


distance is probed. All calculations were carried out using the RICC2–in–RICC2 method and the def2-TZVP basis set. A large reference grid was used to compare analytical and numerical gradients. In Fig. 2 potential energy surfaces of the lowest excited singlet state (S1 ) are displayed for the CO–H2 O complex using the Hartree-Fock and the coupled-cluster density in the embedding contribution. For FT [0], the curves obtained with Hartree-Fock and coupled-cluster densities only differ marginally. Such a small error is obtained because the environment water molecule is not polarized and is thus not affected by the wrong sign of the CO dipole moment. Polarization of the water molecule via FT iterations leads to increased errors. In case of FT [1], the PES of the Hartree-Fock density differs by about 0.05 milli-Hartree with respect to the FT[0] curve in the minimum region, while a shift of about 0.15 milli-Hartree is observed in case of the coupled-cluster density. The PES obtained from FT [2] coincide with the PES obtained with FT [1]. We thus conclude that the errors of FT [0] as introduced using the Hartree-Fock density even for such close intermolecular distances of about 2.5 Å are negligible compared to the method error. A comparison of analytical and numerical gradients is given in Fig. 3. The figure shows that the curves as obtained for FT [0] lie on top of each other. The figure thus shows the correct implementation and agreement of analytical and numerical excited-state RICC2 FDE gradients for FT [0]. In Tab. 1, vertical de-excitation energies at an intermolecular distance of 2.53 Å are shown for different FT cycles employing the Hartree-Fock and coupled-cluster density. Relaxation of the environment increases the difference between Hartree-Fock and coupled-cluster density slightly due to the different polarization stemming from the the wrong sign of the dipole moment in case of the Hartree-Fock density. However, even for FT [3], the resulting transition energies differ only by about 0.002 eV, which is significantly smaller than the method error of RICC2–in–RICC2 embedding. 42 Furthermore, the de-excitation energy obtained from a supermolecule RICC2 calculation is computed to be 7.31 eV, yielding an error of about 0.02

17



Figure 3: Gradient norm of the carbon atom in the CO–H2 O complex computed with analytical and numerical gradients as obtained with RICC2 FDE FT [0] employing the Hartree-Fock density for the embedding contributions. eV due to the FDE approximation, which is well below the method error of RICC2. 57 In order to judge these differences it is helpful to compare also the total environment shift. The isolated CO molecule shows a VDE of 7.234 eV when re-optimized in the excited state, yielding a solvatochromic shift of about 0.07 eV in case of the supermolecule result in Tab. 1. Due to the difference of 0.02 eV, the FDE methods using both the HF and CC2 density yield solvatochromic shifts of about 0.10 eV.

4.2

4-Aminophthalimide-water complex

The emission spectrum of 4-aminophthalimide (4AP) is known to be sensitive to the environment, particularly if hydrogen bonding is possible. For the isolated 4AP molecule, the experimentally observed 0-0 transition of the electronic ground state and lowest singlet excited state in a supersonic jet was measured to be 3.59 eV. 58,59 Chen et. al. have reported a 4AP–(H2 O)2 cluster in which the 0-0 experiences a red shift of 0.31 eV compared to the 18


Page 18 of 41

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


Table 2: Vertical excitation energies (VEE), vertical de-excitation energies (VDE), and adiabatic excitation energies (AEE), all in eV, as well as the dipole moments in Debye of the 4-aminophthalimide monomer. The experimentally observed 0-0 transition in a supersonic jet was measured to be 3.59 eV. 58,59 Method ADC(2)

d

Statea Basis VEE VDE AEEb v-AEEc Dipoled S1 def2-SVP 4.259 3.611 3.924 3.808 9.6459 S1 def2-SVPD 4.083 3.439 3.744 3.628 9.9102 S1 def2-TZVP 4.149 3.448 3.784 3.668 9.8176 RICC2 S1 def2-SVP 4.235 3.592 3.903 3.787 9.5772 S1 def2-SVPD 4.053 3.415 3.715 3.599 9.9204 S1 def2-TZVP 4.118 3.424 3.759 3.643 9.7413 e MC-QDPT2 S1 10.6 a State label corresponds to the ground-state geometry. b No vibrational correction included. c Vibrational correction of −0.116 eV included. Norm of the (orbital-relaxed) dipole moments of the 4AP molecule in the optimized excited state in Debye. e Taken from Ref. 59.

isolated 4AP molecule. 60 In order to compute the solvatochromic shift upon complexation, absorption and emission transitions need to be determined for the isolated molecule and the complex. In Tab. 2 results are displayed for the isolated 4AP monomer. For example, an adiabatic excitation energy (AEE) without vibrational corrections of 3.759 eV at RICC2/def2-TZVP level of theory is obtained corresponding to an error of about 0.2 eV compared to the experimental value. For this case, the zero-point vibrational energies are 131.50 mEh and 127.23 mEh for the ground state and excited state, respectively, yielding a vibrational correction of −0.116 eV. Including this correction leads to an AEE of 3.643 eV which exhibits an error of about 0.05 eV to the experimental value. AEE including the same correction of −0.116 eV are displayed in the column v-AEE in Tab. 2. In Tab. 3 results are collected for the 4AP–water complex. The table reveals that the supermolecular calculations are in general in good agreement with the experimental value. In particular when including a vibrational correction of −0.082 eV obtained at RICC2/def219



Page 20 of 41

Table 3: Vertical excitation energies (VEE), vertical de-excitation energies (VDE), adiabatic excitation energies (AEE), and solvatochromic shift (∆) in eV and dipole moments in Debye of the 4-aminophthalimide complex with two water molecules [4AP–(H2 O)2 ]. The experimentally observed electronic origin transition for the two-water complex is 3.28 eV, yielding a red shift of 0.31 eV compared to the isolated molecule. 60 Method ADC(2)

Basis VEE VDE AEEa ∆a,c v-AEEb v-∆b,c 4AP dipoled def2-SVP 3.778 3.088 3.424 −0.50 3.342 −0.47 12.24e def2-SVPD 3.621 2.971 3.281 −0.46 3.199 −0.43 13.40e def2-TZVP 3.691 2.987 3.355 −0.43 3.273 −0.40 13.69e RICC2 def2-SVP 3.766 3.104 3.430 −0.47 3.348 −0.44 11.96e def2-SVPD 3.603 2.978 3.279 −0.44 3.197 −0.40 13.18e def2-TZVP 3.669 2.993 3.318 −0.44 3.236 −0.41 13.40e ADC(2)–in–PBE def2-SVP 3.936 3.305 3.609 −0.32 3.527 −0.28 11.12 ADC(2)–in–PBE def2-SVPD 3.754 3.148 3.438 −0.31 3.356 −0.27 11.48 ADC(2)–in–PBE def2-TZVP 3.795 3.138 3.470 −0.31 3.388 −0.28 11.56 RICC2–in–PBE def2-SVP 3.915 3.322 3.608 −0.30 3.526 −0.26 11.08 RICC2–in–PBE def2-SVPD 3.715 3.141 3.410 −0.31 3.328 −0.27 11.49 RICC2–in–PBE def2-TZVP 3.762 3.137 3.394 −0.36 3.313 −0.33 11.47 a No vibrational correction included. b Vibrational correction of −0.081 eV included. c Solvatochromic shift in eV obtained as AEE difference with respect to the monomer result in the corresponding basis, cf. Tab. 2. d Norm of the (orbital-relaxed) dipole moments of the 4AP molecule in the optimized excited state in Debye. e Obtained from a supermolecular ESP fit of the (orbital-relaxed) excited-state density as implemented in the KOALA program.

TZVP level of theory for all supermolecular calculations, both ADC(2) and RICC2 exhibit small errors of about 0.05 eV or less. However, the errors in the solvatochromic shifts ∆ and v-∆ are about 0.1 eV which results from the sum of errors for monomer and complex. Using the FDE method, absolute AEE deviations are increased slightly to about 0.1 eV compared to the supermolecular results, while the solvatochromic shifts ∆ and v-∆ exhibit a decreased deviation due to fortunate error compensation. Tab. 3 also shows dipole moments for the 4AP molecule obtained from both supermolecular and FDE calculations. Using FDE, molecular properties such as dipole moments are obtained as a by-product due to the subsystem partitioning, while no unambiguous method 20


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


exists to extract general subsystem properties from supermolecular calculations. In order to obtain an estimate of the subsystem dipole moment of the 4AP molecule in the 4AP– water complex, an electro-static potential (ESP) fit 61 was carried out for the orbital-relaxed excited-state dipole moment from which subsequently only the atoms of the 4AP molecule have been used to determine the 4AP molecular dipole moment. The table shows a strong basis-set dependence for the estimated dipole moments from supermolecular calculations. For example, the difference between the def2-SVP and def2-SVPD basis in case of ADC(2) is more than 1 Debye, while for the monomer and FDE calculations only a weak basis-set dependence with less than 0.5 Debye is obtained. We thus conclude that the FDE approximation is a viable tool to study subsystem properties of molecules in complex environments while subsystem properties estimated from supermolecular calculations might need further analysis. The geometries for the ground-state and the first singlet excited states as obtained from a RICC2–in–PBE/def2-TZVP optimization were the subject of further investigations, for which RICC2–in–PBE and RICC2–in–MP2 single-point calculations were carried out with the def2-TZVPPD basis set, up to two freeze-thaw iterations and both the SCF and CC2 density for the target subsystem. Results are compiled in Tab. 4. Carrying out one FT iteration improves the vertical transition energies by about 0.08 eV, while the second FT iteration yields a further change of about 0.01 eV. In contrast, the choice of the density for the active system has negligible influence on the result. To give an example, a maximum change of 0.002 eV was observed when using the SCF density instead of the CC2 density, while the choice for the environment method yields differences of about 0.01 eV. Here, the vertical transition energies are overestimated in all FDE calculations compared to the supermolecule RICC2 result. The table also reveals that increasing the number of FT iterations can help to reduce the errors and the FT[2] results are the closest to the supermolecular calculations. For RICC2–in–RIMP2 the difference amounts up to 0.026 eV, whereas for the RICC2–in–PBE case the error is slightly larger with up to 0.036 eV. In case of adiabatic transition energies,

21



Table 4: Comparison of vertical excitation energies (VEE), vertical de-excitation energies (VDE) and adiabatic excitation energies (AEE) obtained using CC2 and SCF density (def2TZVPPD, grid 3) for the 4AP water complex optimized at RICC2–in–PBE/def2-TZVP level of theory in eV. The experimentally obtained AEE is measured to be 3.28 eV. 60 Method RICC2–in–PBE RICC2–in–PBE RICC2–in–PBE RICC2–in–PBE RICC2–in–PBE RICC2–in–PBE RICC2–in–RIMP2 RICC2–in–RIMP2 RICC2–in–RIMP2 RICC2–in–RIMP2 RICC2–in–RIMP2 RICC2–in–RIMP2 RICC2e

e

FTa Densityb VAE VEE AAEc v-AAEd FT[0] SCF 3.748 3.125 3.423 3.342 FT[1] SCF 3.664 3.045 3.356 3.275 FT[2] SCF 3.655 3.036 3.348 3.267 FT[0] CC2 3.748 3.125 3.422 3.341 FT[1] CC2 3.667 3.047 3.357 3.276 FT[2] CC2 3.657 3.038 3.350 3.269 FT[0] SCF 3.732 3.109 3.444 3.363 FT[1] SCF 3.653 3.034 3.379 3.298 FT[2] SCF 3.645 3.026 3.371 3.290 FT[0] CC2 3.732 3.109 3.443 3.362 FT[1] CC2 3.655 3.036 3.379 3.298 FT[2] CC2 3.646 3.028 3.371 3.290 3.621 3.002 3.358 3.277 a Number of freeze-thaw iterations. b The density used for the active subsystem. c Without ZPVE. d With ZPVE of −0.081 eV. Using the same geometry optimized at RICC2–in–PBE/def2-TZVP level of theory. Calculation carried out with the KOALA program.

the conventional RICC2 calculations including ZPVE show reasonable agreement with the experiment with a deviation of about 0.005 eV. Including ZPVE, RICC2–in–PBE FT[2] and RICC2–in–RIMP2 FT[2] deviate from the experimental adiabatic transition energy by about 0.01 eV only.

4.3

The benzonitrile dimer

Investigating the potential energy surface (PES) of the benzonitrile dimer poses a further challenge for the FDE ansatz. Despite the fact that subsystems are computed at coupledcluster level of theory in the present work, all inter-subsystem interactions are described based on DFT. The deficiencies of FDE to accurately describe PESs are well known and can 22


Page 22 of 41

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


Table 5: Vertical excitation energies (VEE) and vertical de-excitation energies (VDE) in eV, as well as dipole moments in Debye of the benzonitrile dimer as obtained in supermolecular calculations without embedding.

a

Method ADC(2) ADC(2) ADC(2) ADC(2) RICC2 RICC2 RICC2 RICC2 State label

Basis Statea VEE def2-SVP S1 5.094 def2-TZVP S2 5.059 def2-TZVPPD S2 5.031 aug-cc-pVTZ S2 5.021 def2-SVP S1 5.071 def2-TZVP S2 5.034 def2-TZVPPD S2 5.007 aug-cc-pVTZ S2 4.995 corresponds to the ground-state

VDE 4.738 4.703 4.679 4.670 4.713 4.675 4.652 4.642 geometry.

often traced back to the quality of the kinetic energy which is computed from approximate functionals. 14,17,62 In Tabs. 5 and 6 results are collected for supermolecule and FDEu vertical excitation energies at the optimized ground-state geometry and vertical de-excitation energies at the optimized excited-state geometries. The results reveal that for such a weakly interacting system the FDE approximation yields rather accurate excitation energies with maximum errors of about 0.02 eV as the excitation are purely localized on one subsystem. The overall error is thus consistent with earlier TDA–in–DFT excited state geometry optimizations. For example, the benzonitrile dimer has been investigated using TDA–in–DFT at B3LYP–in– B3LYP and PBE–in–PBE level of theory, for which errors of 0.01 eV or less have been obtained with respect to supermolecular results. 17 Despite the small increase of +0.01 eV in the FDE error for CC2 and ADC(2) compared to DFT FDE, errors of 0.01 eV and 0.02 eV are nevertheless one order of magnitude smaller than the method errors of TDA and CC2, respectively. Furthermore, the difference of the def2-TZVP and the aug-cc-pVTZ results is about 0.04 eV and thus approximately twice as large as the error introduced by FDE in this case. In Tabs. 7 and 8 inter- and intramolecular distances are displayed for geometries ob23



Table 6: Vertical excitation energies (VEE) and vertical de-excitation energies (VDE) in eV, as well as dipole moments in Debye of the benzonitrile dimer.

b

Method Basis Statea VEE VDE Dipoleb ADC(2)–in–PBE def2-SVP S1 5.105 4.753 4.547 ADC(2)–in–PBE def2-TZVP S1 5.070 4.718 4.972 ADC(2)–in–PBE def2-TZVPPD S1 5.045 4.697 4.964 ADC(2)–in–PBE aug-cc-pVTZ S1 5.032 4.686 4.989 ADC(2)–in–RIMP2 def2-SVP S1 5.101 4.745 4.580 ADC(2)–in–RIMP2 def2-TZVP S1 5.069 4.716 4.992 ADC(2)–in–RIMP2 def2-TZVPPD S1 5.042 4.691 4.992 ADC(2)–in–RIMP2 aug-cc-pVTZ S1 5.031 4.681 5.086 RICC2–in–PBE def2-SVP S1 5.055 4.728 4.578 RICC2–in–PBE def2-TZVP S1 5.045 4.690 4.989 RICC2–in–PBE def2-TZVPPD S1 5.020 4.670 4.987 RICC2–in–PBE aug-cc-pVTZ S1 5.007 4.659 5.015 RICC2–in–RIMP2 def2-SVP S1 5.081 4.720 4.629 RICC2–in–RIMP2 def2-TZVP S1 5.045 4.690 5.020 RICC2–in–RIMP2 def2-TZVPPD S1 5.018 4.665 5.026 RICC2–in–RIMP2 aug-cc-pVTZ S1 5.007 4.654 5.100 a State label corresponds to the ground-state geometry. Norm of the (orbital-relaxed) dipole moments of subsystem I in the optimized excited state in Debye.

Table 7: Selected intermolecular and intramolecular bond lengths in Å in the state S1 of the benzonitrile dimer as obtained from supermolecular calculations.

a

Method Basis N(I) – H(II) H(I) – N(II) C1(I) – C2(I)a ADC(2) def2-TZVP 2.349 2.349 1.439 ADC(2) def2-TZVPPD 2.308 2.315 1.439 ADC(2) aug-cc-pVTZ 2.284 2.305 1.439 RICC2 def2-TZVP 2.329 2.334 1.442 RICC2 def2-TZVPPD 2.289 2.299 1.442 RICC2 aug-cc-pVTZ 2.281 2.283 1.442 Intramolecular bond length with the largest change upon excitation, cf. Tab. 5 of Ref. 63.

24


Page 24 of 41

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


Table 8: Selected intermolecular and intramolecular bond lengths in Å in the state S1 of the benzonitrile dimer. The deviation with respect to the supermolecule optimization is given in parentheses, cf. Tab. 7.

a

Method Basis N(I) – H(II) H(I) – N(II) C1(I) – C2(I)a ADC(2)–in–PBE def2-TZVP 2.180 (–0.169) 2.228 (–0.121) 1.438 (–0.001) ADC(2)–in–PBE def2-TZVPPD 2.179 (–0.129) 2.230 (–0.085) 1.438 (–0.001) ADC(2)–in–PBE aug-cc-pVTZ 2.180 (–0.104) 2.229 (–0.076) 1.438 (–0.001) ADC(2)–in–MP2 def2-TZVP 2.289 (–0.060) 2.342 (–0.007) 1.438 (–0.001) ADC(2)–in–MP2 def2-TZVPPD 2.215 (–0.093) 2.310 (–0.005) 1.439 (+0.000) ADC(2)–in–MP2 aug-cc-pVTZ 2.239 (–0.045) 2.126 (–0.179) 1.438 (–0.001) RICC2–in–PBE def2-TZVP 2.175 (–0.154) 2.224 (–0.110) 1.441 (–0.001) RICC2–in–PBE def2-TZVPPD 2.172 (–0.117) 2.227 (–0.072) 1.441 (–0.001) RICC2–in–PBE aug-cc-pVTZ 2.172 (–0.109) 2.226 (–0.057) 1.441 (–0.001) RICC2 –in–MP2 def2-TZVP 2.155 (–0.174) 2.312 (–0.022) 1.441 (–0.001) RICC2 –in–MP2 def2-TZVPPD 2.146 (–0.143) 2.317 (–0.018) 1.442 (+0.000) RICC2 –in–MP2 aug-cc-pVTZ 2.249 (–0.032) 2.179 (-0.104) 1.441 (–0.001) Intramolecular bond length with the largest change upon excitation, cf. Tab. 5 of Ref. 63.

tained using supermolecular calculations and FDE FT [0], respectively. It can be seen that a PBE environment underestimates the repulsion between the monomers yielding too short intermolecular distances. The overall error in the geometries compared to the supermolecule calculations is less than 0.2 Å which is also consistent with earlier RICC2–in–RICC2 and DFT–in–DFT ground state and excited state geometry optimizations. 14,17 Using basis sets with diffuse functions such as Dunning’s aug-cc-pVTZ basis, errors of 0.01 Å are obtained for both ADC(2) FDE and CC2 FDE with respect to the supermolecular geometry optimizations. However, it must be emphasized that a strict comparison of supermolecular and FDE optimizations is not straight forward because both methods as employed here suffer from different intrinsic drawbacks. On one hand, in FDE approximate kinetic and exchange-correlation energy functionals are applied leading to an approximate treatment of intermolecular interactions. On the other hand, the supermolecular optimizations suffer from basis-set superposition error (BSSE), which are avoided in the FDE calculations since a monomer expansion is used in the present work.

25



4.4

Chlorophyll P700 dimer

In order to illustrate the effect of the chosen environment method upon molecular properties, the P700 type pigment dimer is presented. 33 Results for different methods obtained using the def2-TZVP basis are collected in Tab. 9. The table reveals a difference of about 0.2 eV for RICC2 and ADC(2), which corresponds to the difference in the monomers for which 2.17 eV and 1.98 eV are obtained using RICC2 and ADC(2), respectively. 33 While for the monomer different methods can be used, Suomivuori et al. only employed the semi-empirically scaled reduced virtual space (RVS) LT-SOS-RVS-ADC(2)/def2-TZVP method for the dimer (and tetramer). In the present work the RICC2 FDE method is used to study the chlorophyll subsystems, allowing to avoid the shift of 0.2 eV in the monomer. However, while the absolute values of RICC2 and ADC(2) differ, the resulting splitting of about 0.03–0.04 eV is similar in all cases and agrees with the LT-SOS-RVS-ADC(2) shift of 0.05 eV. 33 Tab. 9 reports also excited-state subsystem dipole moments for the RICC2 FDE and ADC(2) FDE method. In these calculations, additionally an ESP fit has been carried out after the relaxed excited-state density was computed. The resulting excited-state atomic charges are displayed in Fig. 4, showing positive and negative atomic charges in blue and red, respectively. Light colors indicate a small charge and intense colors a strong charge, revealing that the choice of PBE or MP2 as the environment method has only little influence. Table 9: Vertical excitation energies (VEE) and relaxed dipole moments after vertical excitation employing different embedding schemes for the investigated chlorophyll cluster, cf. Fig. 4.

b

Subsystem I Subsystem II Method Exc.a Dip.b Exc.a Dip.b ADC(2)–in–MP2 2.01 2.507 2.05 2.390 ADC(2)–in–PBE 2.01 2.482 2.04 2.375 CC2–in–MP2 2.21 2.696 2.24 2.635 CC2–in–PBE 2.21 2.673 2.24 2.612 a Vertical excitation energy in eV. Relaxed dipole moment in atomic units after vertical excitation to S1 .

26


Page 26 of 41

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


Figure 4: Atomic charges as obtained from an ESP fit for subsystem I (top) and subsystem II (bottom) in the chlorophyll cluster, cf. Tab. 9. Blue and red indicate positive and negative atomic charges, respectively.

4.5

Rhodopsin model

In case of biological environments, molecules in close proximity to the chromophore can introduce major perturbations leading to significantly altered properties of the biologically active compound. Here, a rhodopsin model is studied in terms of excitation energies and relaxed dipole moments. 21,42 The geometry, comprising 32 molecules with 366 atoms in total, has been taken from Ref. 21 and no further optimization was carried out. In the following, the def2-TZVP basis is used for the chromophore and the def2-SVP basis for the environment unless stated otherwise. 64–66 Results for the rhodopsin model are collected in Tab. 10. For example, using ADC(2)– in–PBE yields an absorption energy of 2.50 eV of the chromophore inside the biological model environment. Removing the environment molecules yields an absorption energy of 1.97 eV employing ADC(2), thus yielding a shift of +0.53 eV beyond the pure geometry 27



Table 10: Vertical excitation energies (VEE) and relaxed dipole moments after vertical excitation for rhodopsin, cf. Fig. 5. Monomer values are computed at the supermolecular geometry from which the environment molecules have been removed. Complexes are computed using FDE with a PBE environment, e.g., ADC(2)–in–PBE.

b

VEEa Dipoleb Method Monomer Complex shift Monomer Complex shift CIS(D∞ ) 1.97 +0.53 15.377 +2.511 ADC(2) 1.97 +0.53 15.380 +2.509 CC2 2.10 +0.51 15.107 +2.562 a Vertical excitation energy in eV. Relaxed dipole moment in atomic units after vertical excitation to S1 .

effects which are not investigated in the present work. In order to analyze the electrostatic effects, the relaxed dipole moment after (vertical) excitation was computed. Tab. 10 shows that the environment shift is about 2.5 atomic units corresponding to about 6.3 Debye. From the relaxed density, further analysis can be carried out. An ESP fit, for example, provides atomic charges from which regions can be identified that experience a significant change of electron density upon excitation. The atomic charges as obtained from an ESP fit using the ADC(2)–in–PBE method using the conductor-like screening model (COSMO) are displayed in Fig. 5. 67,68 However, as briefly discussed for the 4AP molecule, atomic charges require further analysis in general. A discussion concerning the physical meaning of the atomic charges for this particular system is beyond the scope of the current manuscript. While conventional RICC2 exhibits fifth order scaling computational costs with respect to the total system size, the method as reported in the present work exhibits fifth order scaling computational costs only with respect to the size of the individual subsystems but a linear scaling with respect to the number of subsystems due to the FDE partitioning which can be illustrated unambiguously using computation times. For example, even with the very efficient TURBOMOLE program, a supermolecular RICC2 calculation (without further approximations) of this biological model system containing 366 atoms and 3601 basis functions was not feasible on computers available to us, cf. Tab. 11. Send et al. have

28


Page 28 of 41

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


Figure 5: Atomic charges of the chromophore in a biological rhodopsin model as computed using ADC(2)–in–PBE with COSMO for the first excited state after vertical excitation. The colors of the atoms displayed using their van der Waals radii denote the charge (re-) distribution upon excitation: Blue indicates a charge accumulation, red indicates a charge reduction. reported large supermolecule RICC2 excitation energy calculations (without properties) of a rhodopsin model system in the def2-TZVP basis containing 165 atoms (2565 contracted spherical basis functions) which required 13.6 days wall time on 64 cores. 69,70 These timings show on one hand that selected case studies are possible but, on the other hand, that these are beyond routine application on standard computer clusters in particular when the system sizes are further increased, e.g. to 366 atoms in the present model, or a number of different configurations needs to be studied. Based on the computation time of the isolated chromophore at RICC2/def2-TZVP level of theory using TURBOMOLE and assuming an ideal quintic scaling we estimated a wall time of about 100 days using 5 cores of an Intel® Xeon® E5-2630 v4 @ 2.20GHz processor. Using FDE, this demanding supermolecular calculation is replaced with 2 × 32 significantly smaller subsystem calculations in the present case of rhodopsin and one freeze-thaw iteration. For example, the RICC2–in–PBE calculation on the rhodopsin model using the Hartree-Fock

29



Page 30 of 41

Table 11: Number of contracted spherical atomic orbitals (AOs) for the different computational setups for the rhodopsin model. Employing the FDE approach an explicit computation of the supermolecule is avoided. Supermolecule Chromophore a Method def2-SVP/def2-TZVP def2-TZVP def2-TZVPP def2-TZVP def2-TZVPP Hartree-Fock 3601 5814 7605 1003 1470 RI-JK (universal) 14989 14989 14989 2561 2561 Coulomb fitting 10800 14385 17745 2470 3040 a def2-TZVP basis for chromophore, def2-SVP basis for environment molecules.

density takes about 4.3 days on 5 cores on the same Intel® Xeon® E5-2630 v4 @ 2.20GHz processor including all (repeated) Kohn-Sham, Hartree-Fock and coupled-cluster iterations including excitation energies and the calculation of relaxed excited-state properties for the chromophore. Using ADC(2)–in–PBE, the computation time for this system, consisting of 366 atoms in total, is even reduced to 2.7 days on 5 cores on the same machine. Additionally, the partitioning of the supermolecule also reduces disk and memory requirements significantly. The full strength of the new method becomes evident in particular when the basis is further increased. Employing the def2-TZVP basis for the environment molecules leads to a supermolecular calculation with 5814 AOs, cf. Tab. 11. On the same Intel® Xeon® E5-2630 v4 @ 2.20GHz processor, the SCF procedure takes a factor of about 4.7 longer compared to the calculation time needed employing the def2-SVP basis for the environment molecules, both obtained using TURBOMOLE. This corresponds to an effective scaling of about O(N 3.2 ) in this particular case, where N is a measure of the system size. Assuming an ideal quintic scaling for conventional RICC2, the coupled-cluster part would take an additional factor of 10 longer compared to the calculation which employs a def2-SVP basis for the environment molecules. We want to emphasize that this increase in computation time is not related to the efficiency of the implementation but arises solely due to the quintic scaling of conventional RICC2. However, in case of frozen density embedding, which avoids

30


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


the bottleneck of a quintic scaling with the total system size, the computation time for the chromophore remains constant and only the time for each environment subsystem calculation increases slightly when the def2-TZVP basis is used for all molecules, leading effectively to a linear scaling with respect to the number of subsystems. This advantage becomes crucial if even larger basis sets are required. For example, employing the def2-TZVPP basis for all molecules would result in a supermolecular calculation with 7605 basis functions, cf. Tab. 11. In the case of FDE, the computation of the chromophore, being the largest subsystem in this model, includes 1470 basis functions in this case, leading to a total computation time of 4.2 days on 5 cores for excited-state properties using ADC(2)–in–PBE. The new method can thus be used to efficiently study excited-state properties of medium-sized molecules in complex environments using ab-initio methods as routine applications within a couple of days.

5

Summary

In the present work, we have presented the derivation and implementation of analytical nuclear gradients for excited states using the second-order approximate coupled-cluster singles and doubles model combined with uncoupled frozen-density embedding using density fitting. As a by-product, methods which are related to CC2 such as ADC(2) or CIS(D∞ ) are also available. In order to provide exact analytical gradients without approximations the Hartree-Fock density is employed for the target subsystem in the embedding contributions. The method is assessed using three critical test cases. In case of the CO–H2 O complex the Hartree-Fock dipole moment exhibits a wrong sign but the results show that nevertheless accurate results are obtained. In case of the 4AP–water complex solvatochromic shifts for the adiabatic excitation energy is computed with an accuracy of about 0.1 eV. In case of the weakly-interacting benzonitrile dimer, a maximum error for the vertical transition energies of about 0.02 eV for all methods and basis sets investigated is observed. For example in

31



case of the 4AP–water complex it was found that the adiabatic excitation energies can be improved further when the additional freeze-thaw cycles are carried out at the converged geometry. When comparing FDE with supermolecular calculations two key differences must be kept in mind. First, in case of FDE the intermolecular contributions are computed using approximate kinetic and exchange-correlation functionals independently of the methods used for the subsystems. Second, while no counterpoise correction was used to correct the BSSE in the supermolecular calculations, a monomer expansion is used for the FDE calculations avoiding BSSE. The application of the new method to study subsystem properties of molecules in complex environments reveals two key aspects. First, due to the FDE approximation it becomes straightforward to study local properties of target molecules such as particular chromophores in complex biological environments, allowing for further analysis in terms of chemicallyoriented concepts. Second, the method as reported provides a significant speed-up due to the Lagrangian which is linear in the coupled-cluster multipliers, reducing the computational complexity significantly. Additionally, the reduced dimensions increase numerical stability and significantly reduce disk and memory requirements. The new method can thus be used to efficiently study excited-state properties of medium-sized molecules in complex environments using ab-initio methods as routine applications within a couple of days.

Acknowledgement S. H. gratefully acknowledges financial support from the Deutsche Forschungsgemeinschaft DFG (HO-4605/2-1).

References (1) Sierka, M.; Hogekamp, A.; Ahlrichs, R. Fast evaluation of the Coulomb potential for electron densities using multipole accelerated resolution of identity approximation. J. 32


Page 32 of 41

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


Chem. Phys. 2003, 118, 9136–9148. (2) Cortona, P. Self-consistently Determined Properties of Solids Without Band-structure Calculations. Phys. Rev. B 1991, 44, 8454–8458. (3) Cortona, P. Direct Determination of Self-consistent Total Energies and Charge-densities of Solids - A Study of the Cohesive Properties of the Alkali-halides. Phys. Rev. B 1992, 46, 2008–2014. (4) Wesolowski, T. A.; Warshel, A. Frozen Density-functional Approach For Ab-initio Calculations of Solvated Molecules. J. Phys. Chem. 1993, 97, 8050–8053. (5) Neugebauer, J. Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory. J. Chem. Phys. 2007, 126, 134116. (6) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesolowski, T. A. The merits of the frozen-density embedding scheme to model solvatochromic shifts. J. Chem. Phys. 2005, 122, 094115. (7) Neugebauer, J.; Jacob, C. R.; Wesolowski, T. A.; Baerends, E. J. An explicit quantum chemical method for modeling large solvation shells applied to aminocoumarin C151. J. Phys. Chem. A 2005, 109, 7805–7814. (8) Neugebauer, J.; Baerends, E. J. Exploring the ability of frozen-density embedding to model induced circular dichroism. J. Phys. Chem. A 2006, 110, 8786–8796. (9) S. P. Gomes, A.; Jacob, C. R. Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems. Annu. Rep. Prog. Chem., Sect. C 2012, 108, 222–277. (10) Neugebauer, J. Chromophore-specific theoretical spectroscopy: From subsystem density functional theory to mode-specific vibrational spectroscopy. Phys. Reports 2010, 489, 1 – 87. 33



(11) Wesolowski, T. A.; Goursot, A.; Weber, J. Properties of CO adsorbed in ZSM5 zeolite: Density functional theory study using the embedding scheme based on electron density partitioning. J. Chem. Phys. 2001, 115, 4791–4797. (12) Dułak, M.; Kamiński, J. W.; Wesołowski, T. A. Equilibrium Geometries of Noncovalently Bound Intermolecular Complexes Derived from Subsystem Formulation of Density Fu nctional Theory. J. Chem. Theory Comput. 2007, 3, 735–745. (13) Iannuzzi, M.; Kirchner, B.; Hutter, J. Density functional embedding for molecular systems. Chem. Phys. Lett. 2006, 421, 16 – 20. (14) Heuser, J.; Höfener, S. Wave-function frozen-density embedding: Analytical nuclear ground-state gradients. J. Comput. Chem. 2016, 37, 1092–1101. (15) Schlüns, D.; Franchini, M.; Götz, A. W.; Neugebauer, J.; Jacob, C. R.; Visscher, L. Analytical gradients for subsystem density functional theory within the slater-functionbased amsterdam density functional program. J. Comput. Chem. 2016, 38, 238–249. (16) Kovyrshin, A.; Neugebauer, J. Analytical gradients for excitation energies from frozendensity embedding. Phys. Chem. Chem. Phys. 2016, 18, 20955–20975. (17) Heuser, J.; Höfener, S. Analytical nuclear excited-state gradients for the Tamm Dancoff approximation using uncoupled frozen-density embedding. J. Comput. Chem. 2017, 38, 2316–2325. (18) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. (19) Peach, M. J. G.; Helgaker, T.; Salek, P.; Keal, T. W.; Lutnaes, O. B.; Tozer, D. J.; Handy, N. C. Assessment of a Coulomb-attenuated exchange-correlation energy functional. Phys. Chem. Chem. Phys. 2006, 8, 558–562. 34


Page 34 of 41

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


(20) Höfener, S.; Trumm, M.; Koke, C.; Heuser, J.; Ekström, U.; Skerencak-Frech, A.; Schimmelpfennig, B.; Panak, P. J. Computing UV/vis spectra using a combined molecular dynamics and quantum chemistry approach: bis-triazin-pyridine (BTP) ligands studied in solution. Phys. Chem. Chem. Phys. 2016, 18, 7728. (21) Zhou, X.; Sundholm, D.; Wesolowski, T. A.; Kaila, V. R. I. Spectral Tuning of Rhodopsin and Visual Cone Pigments. J. Amer. Chem. Soc. 2014, 136, 2723–2726. (22) Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. Chem. Rev. 2005, 105, 4009–4037. (23) Stanton, J. F. Many-body methods for excited state potential energy surfaces. I. General theory of energy gradients for the equation of motion coupled cluster method. J. Chem. Phys. 1993, 99, 8840–8847. (24) Stanton, J. F.; Gauss, J. Many-body methods for excited state potential energy surfaces. II. Analytic second derivatives for excited state energies in the equation of motion coupled cluster method. J. Chem. Phys. 1995, 103, 8931–8943. (25) Stanton, J. F.; Gauss, J. Analytic energy derivatives for the equation of motion coupledcluster method: Algebraic expressions, implementation and application to the S1 state of HFCO. Theor. Chem. Acc. 1995, 91, 267–289. (26) Stanton, J. F.; Gauss, J. The first excited singlet state of s-tetrazine: A theoretical analysis of some outstanding questions. J. Chem. Phys. 1996, 104, 9859–9869. (27) Christiansen, O.; Koch, H.; Jørgensen, P. The 2nd-order Approximate Coupled-cluster Singles and Doubles Model CC2. Chem. Phys. Lett. 1995, 243, 409–418. (28) Hättig, C.; Weigend, F. CC2 excitation energy calculations on large molecules using the resolution of the identity approximation. J. Chem. Phys. 2000, 113, 5154–5161.

35



(29) Köhn, A.; Hättig, C. Analytic gradients for excited states in the coupled-cluster model CC2 employing the resolution-of-the-identity approximation. J. Chem. Phys. 2003, 119, 5021–5036. (30) Hättig, C.; Köhn, A. Transition moments and excited-state first-order properties in the coupled-cluster model CC2 using the resolution-of-the-identity approximation. J. Chem. Phys. 2002, 117, 6939–6951. (31) Hättig, C. In Response Theory and Molecular Properties (A Tribute to Jan Linderberg and Poul Jørgensen); Jensen, H., Ed.; Advances in Quantum Chemistry; Academic Press, 2005; Vol. 50; pp 37 – 60. (32) Hättig, C. Geometry optimizations with the coupled-cluster model CC2 using the resolution-of-the-identity approximation. J. Chem. Phys. 2003, 118, 7751–7761. (33) Suomivuori, C.-M.; Winter, N. O. C.; Hättig, C.; Sundholm, D.; Kaila, V. R. I. Exploring the Light-Capturing Properties of Photosynthetic Chlorophyll Clusters Using Large-Scale Correlated Calculations. Journal of Chemical Theory and Computation 2016, 12, 2644–2651. (34) Wesołowski, T. A. Hydrogen-Bonding-Induced Shifts of the Excitation Energies in Nucleic Acid Bases: An Interplay between Electrostatic and Electron Density Overlap Effects. J. Amer. Chem. Soc. 2004, 126, 11444–11445. (35) Gomes, A. S. P.; Jacob, C. R.; Visscher, L. Calculation of local excitations in large systems by embedding wave-function theory in density-functional theory. Phys. Chem. Chem. Phys. 2008, 10, 5353–5362. (36) Höfener, S.; Gomes, A. S. P.; Visscher, L. Molecular properties via a subsystem density functional theory formulation: A common framework for electronic embedding. J. Chem. Phys. 2012, 136, 044104.

36


Page 36 of 41

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


(37) Corni, S.; Cammi, R.; Mennucci, B.; Tomasi, J. Electronic excitation energies of molecules in solution within continuum solvation models: Investigating the discrepancy between state-specific and linear-response methods. J. Chem. Phys. 2005, 123, 134512. (38) Cammi, R.; Corni, S.; Mennucci, B.; Tomasi, J. Electronic excitation energies of molecules in solution: State specific and linear response methods for nonequilibrium continuum solvation models. J. Chem. Phys. 2005, 122, 104513. (39) Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A state-specific polarizable continuum model time dependent density functional theory method for excited state calculations in solution. J. Chem. Phys. 2006, 125, 054103. (40) Daday, C.; König, C.; Valsson, O.; Neugebauer, J.; Filippi, C. State-Specific Embedding Potentials for Excitation-Energy Calculations. J. Chem. Theory Comput. 2013, 9, 2355–2367. (41) Lunkenheimer, B.; Köhn, A. Solvent Effects on Electronically Excited States Using the Conductor-Like Screening Model and the Second-Order Correlated Method ADC(2). J. Chem. Theory Comput. 2013, 9, 977–994. (42) Heuser, J.; Höfener, S. Communication: Biological applications of coupled-cluster frozen-density embedding. J. Chem. Phys. 2018, 148, 141101. (43) Christiansen, O.; Jørgensen, P.; Hättig, C. Response functions from Fourier component variational perturbation theory applied to a time-averaged quasienergy. Int. J. Quant. Chem. 1998, 68, 1–52. (44) Head-Gordon, M.; Oumi, M.; Maurice, D. Quasidegenerate second-order perturbation corrections to single-excitation configuration interaction. Mol. Phys. 1999, 96, 593–602.

37



(45) Schirmer, J. Beyond the random-phase approximation: A new approximation scheme for the polarization propagator. Phys. Rev. A 1982, 26, 2395–2416. (46) Winter, N. O. C.; Hättig, C. Scaled opposite-spin CC2 for ground and excited states with fourth order scaling computational costs. J. Chem. Phys. 2011, 134, 184101. (47) KOALA, an ab-initio electronic structure program, written by S. Höfener, with contributions from A.-S. Hehn, J. Heuser and N. Schieschke. (48) Höfener, S. Coupled-cluster frozen-density embedding using resolution of the identity methods. J. Comput. Chem. 2014, 35, 1716–1724. (49) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (50) Lembarki, A.; Chermette, H. Obtaining a gradient-corrected kinetic-energy functional from the Perdew-Wang exchange functional. Phys. Rev. A 1994, 50, 5328–5331. (51) 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. (52) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (53) Briggs, E. A.; Besley, N. A. Modelling excited states of weakly bound complexes with density functional theory. Phys. Chem. Chem. Phys. 2014, 16, 14455–14462. (54) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165 – 169. (55) TURBOMOLE V7.2 2017,

a development of University of Karlsruhe and

Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; 38


Page 38 of 41

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


available from http://www.turbomole.com. (56) Treutler, O.; Ahlrichs, R. Efficient molecular numerical integration schemes. J. Chem. Phys. 1995, 102, 346–354. (57) Kánnár, D.; Szalay, P. G. Benchmarking Coupled Cluster Methods on Valence Singlet Excited States. J. Chem. Theory Comput. 2014, 10, 3757–3765. (58) Pryor, B. A.; Palmer, P. M.; Andrews, P. M.; Berger, M. B.; Troxler, T.; Topp, M. R. Spectroscopy of jet-cooled polar complexes of aminophthalimides. Chem. Phys. Lett. 1997, 271, 19 – 26. (59) Sajadi, M.; Obernhuber, T.; Kovalenko, S. A.; Mosquera, M.; Dick, B.; Ernsting, N. P. Dynamic Polar Solvation Is Reported by Fluorescing 4-Aminophthalimide Faithfully Despite H-Bonding. J. Phys. Chem. A 2009, 113, 44–55. (60) Chen, Y.; Topp, M. R. Infrared-optical double-resonance measurements of hydrogenbonding interactions in clusters involving aminophthalimides. Chem. Phys. 2002, 283, 249 – 268. (61) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 5, 129–145. (62) Schlüns, D.; Klahr, K.; Muck-Lichtenfeld, C.; Visscher, L.; Neugebauer, J. SubsystemDFT potential-energy curves for weakly interacting systems. Phys. Chem. Chem. Phys. 2015, 17, 14323–14341. (63) Schmitt, M.; Böhm, M.; Ratzer, C.; Siegert, S.; van Beek, M.; Meerts, W. L. Electronic excitation in the benzonitrile dimer: The intermolecular structure in the {S0} and {S1} state determined by rotationally resolved electronic spectroscopy. J. Mol. Struct. 2006, 795, 234 – 241. 39



(64) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571–2577. (65) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829–5835. (66) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials. Theor. Chem. Acc. 1997, 97, 119–124. (67) Klamt, A.; Schüürmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799–805. (68) Schieschke, N.; Di Remigio, R.; Frediani, L.; Heuser, J.; Höfener, S. Combining frozendensity embedding with the conductor-like screening model using Lagrangian techniques for response properties. J. Comput. Chem. 2017, 38, 1693–1703. (69) Send, R.; Kaila, V. R. I.; Sundholm, D. Reduction of the virtual space for coupledcluster excitation energies of large molecules and embedded systems. J. Chem. Phys. 2011, 134, 214114. (70) Robert Send, personal communcation, 2015.

40


Page 40 of 41

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


Graphical TOC Entry

41


Analytical Nuclear Excited-State Gradients for the Second-Order

Recommend Documents