Corrected Polarizable Embedding: Improving the Induction

Subscriber access provided by READING UNIV

Article

Corrected Polarizable Embedding: Improving the Induction Contribution to Perichromism for Linear Response Theory Heiner Schröder, and Tobias Schwabe J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01033 • Publication Date (Web): 15 Dec 2017 Downloaded from http://pubs.acs.org on December 15, 2017

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 free 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 accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

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

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

Corrected Polarizable Embedding: Improving the Induction Contribution to Perichromism for Linear Response Theory Heiner Schröder and Tobias Schwabe∗ University of Hamburg, ZBH – Center for Bioinformatics and Institute of Physical Chemistry, Bundesstraße 43, 20146 Hamburg, Germany E-mail: [email protected]

Abstract An extension of the Polarizable Embedding (PE) approach for the computation of perichromatic shifts within linear response theory, termed corrected PE, is presented. It covers the change in induction effects in addition to contributions from electrostatics and non-resonant excitonic coupling and thereby presents a combination of the corrected Linear Response and the PE method. Using this method, we analyzed the individual contributions for six different excitations from four molecules in different solvents to clarify the question, which effects should be accounted for by a polarizable solvation model. The (vertical) reference excitation energies are evaluated by the means of full quantum-mechanical computations of large solute-solvent clusters. Excellent agreement is achieved, when both the shift due to the change in induction and non-resonant excitonic coupling in addition to the shift due to electrostatics are accounted for.

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

1

Introduction

The understanding of perichromism, the spectral tuning of a chromophore by its environment, is of great importance for all kinds of spectroscopy. 1,2 Usually spectroscopic experiments are conducted in a condensed phase, where the properties of a molecule (including the optical one) can differ significantly from those in vacuum. Results from spectroscopic experiments are often hard to analyze and the spectroscopist needs to seek aid from theoretical methods or, more precisely, quantum mechanical (QM) methods. The most straightforward way to model environmental effects is to include the whole environment into a QM computation. Unfortunately, this is computationally prohibitive for highly accurate QM methods like coupled cluster (CC) or even density functional theory (DFT), since at least several hundreds of atoms have to be included into the QM region for convergence with respect to the system size. 3 Many electronic transitions are however local in nature. Therefore, the whole system may be separated into a chromophore part, where the (local) excitation takes place, and the environment, which has only an indirect effect and may thus be treated on a lower level of accuracy. In this work, we will deal with the computationally efficient quantum mechanics/molecular mechanics (QM/MM) approach. In its simplest version, it just covers the electrostatic effects of the environment, but a mutual interaction of both subsystem has proven to be crucial for accurate excitation energies. 4–9 This is realized through a polarizable MM region. Two main strategies to compute electronic excitations via QM have emerged. For the first one, the Schrödinger equation is solved explicitly for higher solutions, which yields excited states, and for the second, excited state properties are determined from the response function of the ground state. For exact electronic structure theories, the results of both approaches are identical. 10 In practice, different strategies have to be applied for different electronic structure theories. A consequence of the difference between both approaches becomes apparent, when a non2


Page 2 of 33

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


linear Hamiltonian which depends on its own solutions is used. An example are polarizable solvation models where different results can be found even in the limit of exact states. 11,12 Because some solvent models apply a state-dependent Hamiltonian, they are called statespecific approaches. The reason for their apparently contradictory outcome in contrast to response theory is that a straight forward extension of the pure QM methods to QM/MM methods covers different physical effects on the solute. In general, in a supra-system of non-overlapping charge densities, the perichromatic shift ∆∆Eci0 e0 of a chromophore going from the ground-state c0 to an excited state ci in an environment, which is in its electronic ground state e0 , is the difference between the excitation energy within the environment ∆Eci0 e0 and the one in the free chromophore ∆Eci0 :

∆∆Eci0 e0 = ∆Eci0 e0 − ∆Eci0

(1)

One of us analyzed the interactions of both subsystems using long-range perturbation theory and showed that the perichromatic shift up to second order may be split up into four distinct contributions, 13

disp + Ecexcoupl , + ∆Ecind ∆∆Eci0 e0 = ∆EcCoulomb i e + ∆E i ie i e0 c e0 0 0 0

0

(2)

0

the shift due to (i) the change in the Coulomb interaction with the electrostatic environment, (ii) the change in mutual induction effects, (iii) the change in the London dispersion interaction, and (iv) a non-resonant excitonic coupling. Note that non-overlapping charge densities between the chromophore and the environment were assumed for the derivation and thus no Pauli repulsion is considered, even though it might be crucial in some cases. 14 All solvation models including electrostatic interactions account for contribution (i). The non-resonant excitonic coupling, contribution (iv), is the coupling of the transition (dipole) moment to the environmental response and is—together with contribution (iii)—part of the dispersion interaction of subsystems of which at least one is in an excited state. 12,13,15,16

3



As the transition moment is also related to the oscillator strength, this contribution is of particular importance for experimentally easily observable states. 2,17,18 It is accounted for by solvent models based on response theory. 12 On the other hand, contribution (ii), the (classical) induction effect, is fully taken into account by state-specific approaches. Linear response methods include the mutual polarization of subsystems only in the ground state and do not account for the change of the polarization due to excitation. Instead, the induced dipoles of the ground state enter these methods as an additional electrostatic contribution. Although both approaches describe independent physical effects, linear response was sometimes seen as a less accurate but computationally more efficient alternative to a statespecific approach. 11,19 Corni et al. argued that the neglect of the non-resonant excitonic coupling is justified, since considering dispersion interactions goes beyond the intent of an electrostatic model. 12 By similar arguments, the corrected Linear Response (cLR) method has been developed. 20 Here, the excited state density is computed from response theory and the induction contribution is corrected for in a non-iterative (perturbation like) method based on the change in the density. The cLR method was originally developed for the PCM model and has been recently adopted for an explicit solvent model. 21,22 Such a perturbation scheme was also introduced for the state-specific CIS- 23 and CIS(D)-method 24 in combination with the effective fragment potential (EFP). Both methods have been investigated based on their differences without finding a preference for either. 25,26 In contrast to that, Köhn et al. advocated for an inclusion of both effects based on a general discussion of the physical effects and their different relevance for different systems. They developed a state-specific ADC(2)/COSMO method 27 and added the non-resonant excitonic coupling by including all additional contributions coming from response theory into their working equations. They concluded by comparison with experimental solvent shifts that both contributions should be accounted for. Likewise, Daday et al. found in an thorough study of the green fluorescent protein system that they had to include induction

4


Page 4 of 33

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


and non-resonant excitonic coupling effects fully to reproduce the experimental data. For that, they combined state-specific CASPT2 and linear-response TDDFT results to estimate the final result. 9 Albeit the fact that a general conclusion about the importance of both contributions based on theoretical and numerical considerations is now slowly emerging, the situation is still not satisfying. The evaluation by means of experimental data is prone to error cancellation for implicit as well as explicit solvent models. Implicit solvent models lack any specific solute-solvent interactions like H-bonds or π − π-stacking. Furthermore, all computations are usually done on the basis of a single configuration, thus any statistical effects are missing. These effects may, however, enter explicit solvent models, but those rely on a representative collection of configurations from a molecular dynamics (MD) simulation such that the real system with all its degrees of freedom is described correctly. Car-Parinello MD 28 or semiempirical approaches like density functional tight binding 29 offer a way to gain configurations of QM quality, but are computationally very demanding. In most cases, one has to rely on molecular mechanics simulations using classical force fields. However, the chosen MD setup may significantly affect the results. 30 A better strategy would be a direct comparison to full QM results with a method, which allows to address both relevant contributions to the perichromic shift of the excitation energy. For this, a correction to the induction contribution is derived for the Polarizable Embedding (PE) method by Kongsted and coworkers, 31–33 which is equivalent to the energy correction of the cLR method for the explicit solvent model. 21 In contrast to cLR the nonresonant excitonic coupling is included and the new method is consequently labeled corrected Polarizable Embedding cPE. As in cLR, the shift due to the change in London dispersion (iii) is not accounted for, since it is expected to have a negligible effect for valence transitions. Thus the perichromatic shift for cLR (∆∆EcLR ) is given by

∆∆E cLR = ∆E Coulomb + ∆E ind

5


(3)


Page 6 of 33

while the shift for cPE (∆∆EcP E ) additionally includes the non-resonant excitonic coupling contribution

∆∆E cPE = ∆E Coulomb + ∆E ind + E excoupl .

(4)

Note that all terms are now related to a QM/MM description. Therefore, referencing the electronic state of the environment is obsolete and ∆E X referes to the contribution of effect X to the electronic excitation energy of the chromophore. Further, the polarizabilities used for this study are not frequency-dependent, which is an approximation to the problem. Frequency dependence of the polarizabilities would make the model more complicated without providing substantial numerical advantages, 34 while it might even cause problems for finding the correct poles of the subsystem response function. 35 Another benefit is that isotropic frequency-independent polarizabilities can be obtained easily not only from quantum chemical approaches 36 but also from C6 dispersion coefficients. 37 For reasons given above, we explicitly refrain from using any experimental data as a reference and compare the used methods to a full quantum mechanical treatment of an explicit solute-solvent system only. Of course the final goal is to reproduce or even predict experimental spectra, but this is only meaningful if the error of all approximations can be quantified and one does not rely on error cancellation. This way we can especially circumvent any error related to an inaccurate sampling of the configurational space. Our method is combined with DFT in this work, because it allows for the computation of a significant number of reference data points on structures with hundreds of atoms. Only this should enable a solid assessment of the the individual contributions to the perichromic shift. Several representative chromophores with prototypic electronic excitations in different solvents have been chosen for this study. A very similar correction has been introduced in the effective fragment potential framework but was tested only for one electronic excitation and without comparison to full QM results. 38 Therefore, the objective of the present work is a more thorough assessment of such an approach. After a more comprehensive presentation 6


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


of the underlying theory, results will be discussed in more detail before a final conclusion is given.

2

Theory

To be able to include the full induction contribution to the perichromic shift, the missing terms in linear response theory have to be identified first. For that, remember that the reference state (usually the ground state) already includes the effect of the mutual polarization of the subsystems. In a QM/MM setup in which polarization is described by induced dipoles this leads to induced dipoles due to the multipoles of the MM region, the nuclei in the QM region, and the electrons:

µind = µMM,ind + µnuc,ind + µel,ind

(5)

During an electronic excitation, the former two components do not change (assuming a vertical excitation and no relaxation of the nuclei) and their field affects the excitation energy in the same way as a static dipole field. Response theory accounts for this correctly but does not include the change of the induced dipoles due to the change of the electric field in the excited state. Instead, the field of the ground state induced dipoles enters the response theory formulation. Both factors, the missing change in polarization and a the wrong field from the ground state contribution, leads to an erroneous description of the perichromic shift except for the special case where the change in the electron electric field due to excitation is exactly zero. Both factors can be addressed by correcting the perichromic shift with a polarization energy correction based on the difference density of the excited state (see the Appendix for details). The difference density is the natural quantity to describe first order properties of excited states within response theory. Its derivation for TDDFT has been derived not only for the vacuum case 39 but also in the context of polarizable embedding models. 38,40–42 Here, 7



Page 8 of 33

the necessary equations to obtain the difference density for PE-TDDFT are given. For this purpose, the PE method is briefly recapitulated first in the following section.

2.1

The Polarizable Embedding method

The general concepts and theory of the PE method for DFT in the ground state as well as linear and quadratic response were presented before in second quantization formalism. 31 In the PE method, the dynamic response of the environment to a change in the charge density is modeled by induced dipoles. Isotropic or anisotropic dipole-dipole-polarizabilities are therefore assigned to U MM sites. The energy contribution due to induced dipoles is given as: 1 E ind = − F † RF 2

(6)

where F is a vector containing the electric field components, due to the multipoles, nuclei and electrons, at all polarizable sites,

F † = (Fx1 Fy1 Fz1 Fx2 · · · FzU )

(7)

and R is a classical response matrix in super matrix form: 

−1

 (α1 )   −T(2) 21  R= ..  .   (2) −TU 1

(2) −T12

...

(2) −T1U

(α2 )−1 . . . .. ... .

(2) −T2U

(2)

−TU 2

.. .

. . . (αU )−1

−1        

(8)

where the diagonal elements are the inverse of the polarizability tensors αu at site u and the off-diagonal elements are the dipole-dipole-interaction tensors between site u and v,

2 T(2) uv = ∇u

1 , |ru − rv |

8


(9)

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


with ∇, the Nabla operator, and ru the Cartesian coordinate vector of site u. Consequently, we have for the induced dipole vector:

µind = RF .

(10)

As mentioned, only the electric field due to the electrons is of interest in the following, since the field due to the multipoles and nuclei is constant during an electronic excitation. The electric field and thus the dipoles induced by it are functions of the electron density P, 1 E ind = − F [P]† µind [P] 2

(11)

This equation is central to this work and forms the point of origin for the following derivations. The electric field of the electrons is a first order property which can be obtained as an one-electron expectation value. Applying the respective operator εˆ(i) for electron i and exploiting the symmetry in Eq. (6), one can derive an effective one-electron operator which has to be added to the Kohn–Sham- (or Fock-) operator to determine the self-consistent field ground state:

PE ˆ G[P](i) =−

XX u

εˆαu (i)µ[P]ind αu ,

(12)

α

where a labels a Cartesian component. Again, the density dependency of the effective operator as well as of the induced dipoles is given explicitly. This information is sufficient to turn now to the derivation of the difference density matrix for the PE method.

2.2

The time-dependent PE-Kohn–Sham eigenvalue problem

In the following sections, we follow the notation of Furche and Ahlrichs. 39 The contributions due to the polarizable environment were incorporated in analogy to the work of Scalmani et al. for PCM-TDDFT. 40 Indices i, j, . . . label occupied molecular orbitals (MOs), a, b, . . . virtual and p, q, . . . general MOs. The time-dependent Kohn-Sham eigenvalue problem reads 9



(Λ − Ω∆) |X, Yi = 0.

Page 10 of 33

(13)

The eigenvalues Ω are the excitation energies and the eigenvectors |X, Yi the corresponding transitionen density matrices. Λ and ∆ are operators in supermatrix form, 



A B Λ= , ∗ ∗ B A

  1 0  ∆= , 0 −1

(14)

with matrices A and B, which sometimes are called orbital rotational Hessians. Solution of Eq. (13) requires the evaluation of matrix-vectors products |U, Vi:

|U, Vi = Λ |X, Yi ,

(15)

which is most efficiently performed as:

(U ± V) = (A ± B)(X ± Y),

(16)

where the PE contribution ν PE only enters the (A + B) part,

xc (A + B)iajb = (ǫa − ǫi ) δij δab + 2 (ia|jb) + 2fiajb PE −cx [(ja|ib) + (ab|ij)] + 2νiajb .

(17)

Here, ǫp is the energy of the SCF orbital p, (pq|rs) denotes two-electron integrals in Mulliken notation, cx is the well-known scaling parameter for inclusion of exact Fock exchange and f xc is the linear response contribution due to the exchange-correlation functional. The PE potential is given by

PE νiajb

XXXX ∂GPE ia =− Rαu βv hi| εˆβv |ai hj| εˆαu |bi . = ∂Pjb u α v β 10


(18)

Page 11 of 33 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 the spirit of the PECC-implentation, 32 we define a contraction of the four-index potential matrix with a two-index density matrix-like vector:

GPE ia [X + Y ] =

X jb

PE νiajb (X + Y )ia .

(19)

Note that this is identical to the construction of the effective operator for the SCF part but with a different input density. This allows to reformulate Eq. (16) as:

(U + V) = (A + B)′ (X + Y) + 2GPE [X + Y],

(20)

where the prime just indicates the removal of 2ν P E from the (A + B)-matrix. Solving equation Eq. (13) under the constraint

hX, Y| ∆ |X, Yi − 1 = 0

(21)

yields the excitation energies and corresponding eigenvectors vectors, sometimes also called reponse vectors. In the following section, we will use the eigenvectors to obtain the one-particle difference density matrix ∆P needed to compute the correction of the induction contribution.

2.3

The difference density

The unrelaxed difference density matrix T is directly determined from the response vectors,

1X [(X + Y )ia (X + Y )ib + (X − Y )ia (X + Y )ib ] 2 i i 1 Xh = − (X + Y )ia (X + Y )ja + (X − Y )ia (X + Y )ja 2 a

Tab = Tij

Tia = Tai = 0,

11


(22)


Page 12 of 33

and contributions due to orbital relaxation effects are determined from the modified Z-vector equation: X (A + B)′iajb′ Zjb + 2GPE ia [Z] = −Ria .

(23)

jb′

Before giving the right-hand side R, two additional contractions for a Fock-matrix like object are defined:

+ Hpq [V] =

X

xc 2 (pq|rs) + 2fpqrs − cx [(ps|rq) + (pr|sq)] Vrs

rs

− Hpq [V ] =

X rs

!

+ 2GPE pq [V],

cx [(ps|rq) − (pr|sq)] Vrs .

(24)

(25)

The right-hand side Ria is then given as:

Ria =

X b

X j

+ − (X + Y )ib Hab [X + Y] + (X − Y )ib Hab [X − Y] −

(X + Y )ja Hji+ [X + Y] + (X − Y )ia Hji− [X − Y] +

+ Hia [T] + 2

X

xc giajbkc (X + Y )jb (X + Y )kc .

(26)

jbkc

where gxc denotes the third order functional derivative of the exchange-correlation energy. Solving Eq. (23) and adding the result to the unrelaxed difference density yields the relaxed difference density matrix ∆P:

∆P = T + Z.

(27)

Once ∆P is known, it can be applied to compute the change in the electric field and the correction of the excitation energies Ω based on Eq. (11):

12


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


1 ∆E corr,ind = − F [∆P]† µind [∆P] 2

(28)

The correction is then simply added to the excitation energy. Since the electric field vector F [∆P] enters Eq. (11) quadratically (see Eq. (6)), the correction is always negative and thus lowers the excitation energy.

3

Technical Details

Structures of the solute/solvent-complexes were obtained as snapshots from a molecular dynamics (MD) simulation. For all MD simulations, the solute and solvents were treated as rigid bodies. Snapshots for acrolein in water, 43 pyridazine in water and acetonitrile, 43 and p-nitroaniline (PNA) in water and dichlormethan 30 were taken from previous studies. For octatetraene and acrolein in acetonitrile, an MD simulation was made following the same protocol. For acrolein the same structure as for water in ref. 43 was used. The geometry for trans-octatetraene was optimized at the B3LYP 44,45 -D3(BJ) 46,47 /def2TZVPP 48 /COSMO 49 (ǫ = 2.23) level of theory. The subsequent steps followed the protocol detailed in ref. 30. 120 relevant snapshots have been obtained from each MD simulation. The final clusters for our computations were obtained from the MD snapshots by cutting out a drop around the solute molecule. A solvent molecule was included into the final cluster, if any of its atoms was within a distance of 5˚ A of any of the solute’s atoms. The distance was chosen such that a large number of full QM computations were feasible, while one to two solvation shells were included. All full quantum chemical computations on these clusters were done using a local development version of TURBOMOLE 7.1. 50 Excitation energies were obtained with the ESCF module 51 following a ground state optimization with the DSCF module. 52 We employed the BHLYP functional 53–55 in order to reduce self-interaction errors. This functional has been successfully applied in previous studies on electronic excitation energies and offers a 13



good compromise when range-separated functionals are not available (as is the case in TURBOMOLE). 56 We rely on the TURBOMOLE program because it is highly optimized and allowed us to conduct several hundreds of computations on big solute-solvent clusters within an acceptable timeframe. Although common density functional approximations are not always reliable to compute solvatochromic shifts, 57 it is valid for the present study, because the direct comparison of QM/MM vs full QM results is not affected by general shortcomings in the description of the electronic excitation. The numerical integration grid m4 58 was used for the ground state optimization and the numerical grid 4 for the response calculations. The cc-pVTZ 59 basis set was used for all solute molecules and the SVP 60 basis set for the solvent molecules to reduce the computational cost. Test computations showed that the excitation energies varied only slightly if a larger basis set was used for the solvent (See the SI for details). For the PE-calculations, the solvent atoms were replaced by an embedding potential represented by a multipole expansion up to quadrupoles and anisotropic polarizabilities centered at the nuclei positions, while the solute was treated with the same basis set as for the full QM computations. All embedding potentials were set up using the WHIRPOOL program. 61 The static multipoles and anisotropic dipole-dipole polarizabilities for all solvents have been obtained via the LoProp 36 method and were taken from previous studies. 30,43 All ground states were then optimized with the PE potential included. The missing PE contributions for the relaxed difference density were implemented into the EGRAD module 39,51,62 of our development version of TURBOMOLE, which is used to compute the correction to electronic excitation energies.

4

Results and Discussion

Our analysis was conducted on four different model systems, which exhibit low lying valence excitation with varying proportions of the effects of interest. The model systems are trans-

14


Page 14 of 33

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


octatetraene, trans-arcolein, pyridazine and PNA. A list of all solute/solvent combinations with critical properties can be found in Tab. 1. All properties are obtained from a PETDDFT calculation and averaged over 120 snapshots from an MD simulation. The critical properties in the context of this work, which allowed us to estimate the magnitude of the shift due to non-resonant excitonic coupling and induction are the oscillator strength and the change in dipole moment upon excitation, respectively. The change of the dipole moment is represented by the change of the norm between the ground and excited state and the angle between both dipole vectors. All investigated excitation are either classified as n − π ∗ - or π − π ∗ -excitations. The n − π ∗ -excitations show characteristically low oscillator strenghts, but a significant decrease in the dipole moment between −0.78 and −1.00 a.u.. The π − π ∗ -excitation of octatetraene in tetrachlormethane was chosen, because of its huge oscillator strength but almost no dipole moment in the ground- and excited state. Consequently, the comparatively high angle between both dipole vectors is a numerical artifact, since the difference between both vectors is almost zero. The excitation is thus excellently suited to study the effect of the non-resonant excitonic coupling on the solvent shift. The remaining π − π ∗ -excitations have significant oscillator strengths of around 0.5 a.u. and an increase in the dipole moment, which is around three times higher for PNA compared to acrolein. In these cases both effects might have a significant contribution.

4.1

Analysis of solvent contributions to the excitation energies

In the following sections, vertical excitation energies are shown and the energy shifts due to the different contributions are analyzed. All excitation energies shown in the following, except for full QM energies, are computed on a reference state optimized with PE. For the first energy E0 , all PE-contributions in the linear response equations are neglected and no induction correction is computed, it corresponds to a vacuum LR calculation on a PE-optimized ground state. If one adds the induction correction to this energy, one arrives at the cLR 15



Page 16 of 33

Table 1: Investigated excitations of four molecules given with solvent, excitation character, the change of the Euclidean norm of the dipole moment upon excitation (∆|µ|), the angle between ground- and excited state dipole vector (6 (~µGS , ~µES )) and the oscillator strength in length representation (osc.l.). All numerical properties were obtained from a PEDFT computation and averaged over 120 snapshots from a MD simulation. solute octatetraene acrolein acrolein acrolein acrolein pyridazine pyridazine pyridazine pyridazine PNA PNA

solvent CCl4 H2 O H2 O MeCN MeCN H2 O H2 O MeCN MeCN H2 O DCM

excitation 1 1 Bu 1 1 A′′ 1 1 A′ 1 1 A′′ 1 1 A′ 1 1 B2 1 1 A2 1 1 B2 1 1 A2 1 1 A′ 1 1 A′

character π − π∗ n − π∗ π − π∗ n − π∗ π − π∗ n − π∗ n − π∗ n − π∗ n − π∗ π − π∗ π − π∗

∆|µ| / a.u. 0.02 ± 0.00 −0.91 ± 0.00 0.54 ± 0.00 −1.00 ± 0.00 0.43 ± 0.01 −0.76 ± 0.00 −0.79 ± 0.00 −0.83 ± 0.00 −0.87 ± 0.00 1.42 ± 0.00 1.75 ± 0.03

(~µGS , ~µES ) 19.1 ± 0.7 9.1 ± 0.4 9.6 ± 0.1 15.9 ± 0.2 10.9 ± 0.0 4.3 ± 0.2 3.8 ± 0.0 3.5 ± 0.0 4.0 ± 0.0 1.4 ± 0.1 1.1 ± 0.0 6

osc. l. / a.u. 1.70 ± 0.00 0.00 ± 0.00 0.50 ± 0.00 0.00 ± 0.00 0.51 ± 0.00 0.01 ± 0.00 0.00 ± 0.00 0.01 ± 0.00 0.00 ± 0.00 0.54 ± 0.00 0.50 ± 0.00

excitation energy, which accounts for the shift due the change in electrostatics and induction. The standard PE-energy ∆EPE takes the shift due to non-resonant excitonic coupling into account, but does not provide the correct induction term for the excited state. Using the same induction correction for this energy, gives the cPE exctiation energies ∆EcPE , which should include all relevant contributions.

4.1.1

n − π ∗ -excitations

Vertical excitation energies of the n − π ∗ -excitations are shown in Tab. 2. With no significant oscillator strenght present, and therefore no significant transition dipole moment, non-resonant excitonic coupling is negligible here. Consequently, practically no difference between excitation energies ∆EcLR and ∆EcPE can be observed for any excitation. On the other hand, as expected, the correction of the induction contribution, which enters the cLR and the cPE method, is important to reproduce the results of the full QM computations. Overall, the magnitude of the correction for the n − π ∗ -excitations is in the 16


Page 17 of 33 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 energy and standard error of the n−π ∗ -excitations of acrolein and pyridazine in vacuum, water and acetonitrile. Excitation energies in solution were averaged over 120 snapshots from an MD simulation. All values in eV. solvent water acetonitrile water acetonitrile water acetonitrile

∆E0

∆EcLR ∆EPE ∆EcPE acrolein – (∆Evac = 3.96eV ) 4.38 ± 0.02 4.28 ± 0.02 4.38 ± 0.02 4.28 ± 0.02 4.14 ± 0.01 4.08 ± 0.01 4.14 ± 0.01 4.08 ± 0.01 pyridazine (1st exc.) – (∆Evac = 4.09eV ) 4.51 ± 0.01 4.44 ± 0.01 4.50 ± 0.01 4.43 ± 0.01 4.28 ± 0.01 4.23 ± 0.01 4.27 ± 0.01 4.22 ± 0.01 pyridazine (2nd exc.) – (∆Evac = 4.95eV ) 5.53 ± 0.02 5.42 ± 0.01 5.52 ± 0.02 5.42 ± 0.01 5.18 ± 0.01 5.11 ± 0.01 5.18 ± 0.01 5.11 ± 0.01

∆EFQM 4.29 ± 0.01 4.13 ± 0.01 4.42 ± 0.01 4.25 ± 0.01 5.38 ± 0.02 5.15 ± 0.01

range of 0.06–0.11 eV and is always slightly larger in water than in acetonitrile. The corrected excitation energies show an excellent agreement for the excitation energies in water, but yield a blue shift which is slightly too small for acetonitrile. In all cases, the excitation energies without any direct contributions from the polarizable environment to the excitation energies (∆E0 ) overestimate the blue shift. Consequently, the induction correction always goes into the right direction, albeit it overshots to some extend for acetonitrile. Anyhow, the error in comparison to the full QM computations is not larger than 0.05 eV on average for all systems. 4.1.2

π − π ∗ -excitations

As stated before, the π − π ∗ -excitation of octatetraene in an unpolar solvent is perfectly suited to study the effect of non-resonant excitonic coupling, since there is practically no change in the dipole moment upon excitation, but an exceptionally high oscillator strength. Our computed data Tab. 3 confirms that assumption. While there is a considerable shift of 0.16 eV due to the non-resonant excitonic coupling, no contribution due to the change in induction can be observed. Because cLR is blind to this effect, it does not improve on the description. But PE and cPE both match the full QM reference perfectly. It should also be noted that the effective solvatochromic shift in comparison is Zero here, but this due to 17



Page 18 of 33

cancellation effects from different contributions. So what might seem as an innocent solvent in experiment may just have opposing solvatochromic effects on a choromphore. Table 3: Vertical excitation energies and standard errors of low-lying π − π ∗ -excitation of the considered solute/solvent systems. Excitation energies in solution were averaged over 120 snapshots. All values in eV. solvent

∆E0 ∆EcLR ∆EPE octatetraene – (∆Evac = 4.31eV ) tetrachlormethane 4.47 ± 0.00 4.47 ± 0.00 4.31 ± 0.00 acrolein – (∆Evac = 6.52eV ) water 6.38 ± 0.01 6.33 ± 0.01 6.27 ± 0.01 acetonitrile 6.53 ± 0.00 6.51 ± 0.01 6.47 ± 0.00 para-nitroaniline – (∆Evac = 4.38eV ) water 4.04 ± 0.01 3.97 ± 0.01 3.95 ± 0.02 dichlormethane 4.30 ± 0.01 4.19 ± 0.01 4.22 ± 0.01 para-nitroaniline with 2.2 ˚ A QM water – (∆Evac water 4.03 ± 0.01 3.97 ± 0.01 3.97 ± 0.01

∆EcPE

∆EFQM

4.31 ± 0.00

4.31 ± 0.00

6.22 ± 0.01 6.44 ± 0.01

6.18 ± 0.01 6.38 ± 0.01

3.87 ± 0.02 4.12 ± 0.01 = 4.38eV ) 3.91 ± 0.01

3.93 ± 0.01 4.15 ± 0.01 3.93 ± 0.01

After investigating cases, where only either one of the effects was significant we now turn to two cases, which might combine contributions from both effects, as it is for the 1 1 A′ excitation of acrolein in water. Here, ∆E0 yields a result which is 0.2 eV too large. The induction correction lowers this by 0.05 eV but only the combination of both contributions with cPE gives a result close to the reference: 6.22 vs. 6.18 eV. The picture is simlar for the acetonitrile solution but again, the overall effects are smaller. So far, only statistical averages were presented. To verify that the correction is systematically improving the standard PE model, we have picked the excitation of acrolein in water and plotted the individual deviations of the excitation energies ∆E0 , ∆EP E and ∆EcP E to the full QM excitation energies for all 120 snapshots in Fig. 1. All excitation energies ∆E0 , that just cover the electrostatic contributions, overestimate the excitation energy. Inclusion of non-resonant excitonic coupling (∆EP E ) improves the energies towards the full QM for every configuration, but only the additional inclusion of the induction energy brings almost all of the remaining (one third) data points into the 0.1 eV error range. Notable exceptions are for example configuration 92, for which we could observe resonant exitonic coupling be18


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


tween the 1 1 A′ -excitation of acrolein and an excitation in the surrounding water. Such a coupling may only be modeled within a QM approach, but is currently out of the scope of a classical solvation model. Consequently, the error for all methods is over 0.3 eV for this configuration. Other outliers might suffer from very close contacts of the chromophore and the environment as discussed for PNA later.

Figure 1: Individual errors for all 120 snapshots of acrolein in water. Errors are the difference between the respective energy and the reference. Errors are shown for a) including electrostatic contributions ∆E0 ,b) non-resonant excitonic coupling contributions in addition to a) ∆EP E and c) induction contributions in addition to b) ∆EcP E . The last excitation investigated here, is the 1 A′ excitation of PNA. The induction correction and non-resonant excitonic coupling contribution are again of equal magnitude (0.07– 0.09 eV) in water. Both on their own correct the ∆E0 energy closely towards to the reference, but overcorrect the energy by 0.06 eV together. In contrast, results for dichlormethane as a solvent are very satisfying when cPE is considered. Again, both effects add up for a better result than cLR or pure PE alone. The most probable reason for an overcorrection in the case of water is the neglect of Pauli repulsion, which has also been observed by others 14 and is supported by the fact that the

19



effect gets worse, if a more diffuse basis set is employed for the solute (see SI). Since Pauli repulsion is not considered within the PE approach, the easiest way to include it, is to add a few of the closest solvent molecules to the QM region. Using our distance based scheme described earlier, we included all solvent molecules within 2.2 ˚ A into the QM region. The number of added molecules ranges from 0 to 7 with an average of 3.2 per configuration. While the ∆E0 energy is hardly effected by this, the PE excitation energy which includes the non-resonant excitonic coupling is raised by 0.02 eV while the induction correction is reduced by the same amount. In total, this reduces the error of the fully corrected excitation energy to one third of the former error, 0.02 eV vs. 0.06 eV. We have plotted the individual errors of the cPE energies ∆EcP E without and with a few QM solvents in Fig. 2. The inclusion of a few QM water is most beneficial for the configurations, which have the highest errors for the pure PE environment. Since Pauli repulsion is most significant, when there are very close contacts between the chromophore and the environment, a partial inclusion of the environment into the QM region can significantly improve the results. For the same excitation in dichlormethane, the initial error for ∆E0 is 0.15 eV. Again the non-resonant excitonic coupling (∼ 0.08 eV) and the induction correction (∼ 0.07 eV) have significant contributions of similar magnitude. The respective corrected energies are already in good agreement with the reference, but a combination of all contributions yields a perfect match with the reference. The contribution due to non-resonant excitonic coupling proved to be significant for all π−π ∗ -excitations here, ranging from 0.06-0.16 eV. The best results with very low errors below 0.05 eV could be obtained, when both the non-resonant excitonic coupling and the induction contribution were accounted for simultaneously. The magnitude of both contributions was almost equal for all π−π ∗ excitations considered here, thus no contribution might be regarded as more significant.

20


Page 20 of 33

Page 21 of 33 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 2: Individual errors for all 120 snapshots of PNA in water. Errors are the difference between the excitation energies including electrostatic, induction and non-resonant excitonic coupling contributions and the reference. Errors are shown for a pure PE-environment and with QM water included (2.2 ˚ A cutoff). An additional line, parallel to the x-axis, is included to mark the 0.1 eV error range.

21



5

Conclusion

We have derived and implemented a model to compute the contribution to the perichromatic shift due to induction for TDDFT based on the Polarziable Embedding approach. Using this method, we analyzed the energy contributions to the solvatochromatic shift for six different excitations from four molecules in different solvents. To ensure direct comparability of our data with the reference, we compared it with full QM-calculations of the solute-solvent cluster of the same size. We could show that both—the contribution due to induction and due to non-resonant excitonic coupling—may have a non-negligible effect on vertical excitation energies depending on the nature of the considered excitation. The order of magnitude of both contributions may be in the order of magnitude of the desired accuracy (∼ 0.1 eV), which is the usually accepted ”chemical accuracy” for electronic excitation energies. The best results with very low errors below 0.05 eV could be obtained, when both the non-resonant excitonic coupling and the induction contribution were accounted for simultaneously. The magnitude of both contributions was almost equal for all π−π ∗ excitations considered here, thus no contribution might be regarded as more significant. Excluding contributions due to the non-resonant excitonic coupling on principal, as done for the cLR method, not only gives inferior results, but also produces computational overhead, since the TDKS eigenvalue problem has to be solved twice. Thus the computational cost for cLR is at least around 50% higher compared to cPE. For n − π ∗ -excitations, which have characteristically very low oscillator strengths, the excitonic coupling contribution is almost zero and its computation might be omitted, but the computational advantage gained is negligible. However, the computation of the induction contribution involves the solution of the Z-vector equation, which is about as expensive as the computation of an additional excitation. Unfortunately, except for special cases, like octatetraene here, one cannot predict, when the computation is dispensable. When the chromophore and the environment get in very close contact, Pauli repulsion 22


Page 22 of 33

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


starts to play a significant role. We could show that inclusion of the closest environment into the QM reduces the error. Another approach to the problem might be an effective operator that accounts for Pauli repulsion. 14 Contributions due to the change of the London dispersion interaction between the groundand excited state are not considered by our approach, since it is not expected to add a significant contribution. Tkatchenko presented a method to compute atomic C6 dispersion coefficients from a partition of the full density. 63 Using this scheme with the relaxed excited state density from this work may give rise to an approximate dispersion correction of the excitation energies to proof this assumption. We are aware that our method presents a perturbation correction to linear response methods which is only justified based on arguments from long-range perturbation theory. Nevertheless, we are confident that it should be possible to derive a method fully consistent with response theory. For that, one has to include the induction energy term for the excited state into the excited state Lagrangian and start to derive all working equations from there. Of course, this would couple the equations for the response vectors and the Z-vectors, thereby increasing the size of the numerical problem. Whether this would lead to substancially improved results compared to our corrected PE approach has to be seen. Nevertheless, our approach seems to be more consistent to us than state-specific variants, where the reference (ground) state is also iteratively updated. This is in contradiction to the general idea of response theory that all excited state properties can be derived directly from the reference state. Our findings here should in principal be valid for all polarizable solvation models, including implicit solvent models, where a direct comparison to full QM computations cannot be conducted and assessment of these methods has to be more indirect. Further, combining the corrected PE approach with a high level QM method, like CC2, gets us in a position to assess the quality of MD simulations for perichromatic shifts by direct comparison with experimental data.

23



Page 24 of 33

Acknowledgement We are thankful to Giovanni Scalmani for helpful feedback about the TDDFT/PCM gradient implementation and to Christian Reichardt for suggesting the octatetraene test case.

6

Appendix

In contrast to Eq. (2), from linear response we get for the perichromic shift of the standard PE method:

excoupl es LR + ∆E(0)ind = ∆Ei0 ∆∆Ei0 i0 + Ei

(29)

Here, ∆E(0)ind i0 emphasizes that the change in the interaction energy is taken with respect disp to the (frozen) ground state induced dipoles. As ∆Ei0 is assumed to be negligible for most

cases, this seems to be the largest difference between both solutions. For reasons which become obvious in the following, we split the induced dipole interaction in to contributions from non-electronic and electronic induced dipoles:

noel,ind el,ind ind ∆Ei0 = ∆Ei0 + ∆Ei0

(30)

Note that for a given state n: 1 noel,ind noel,ind µel,ind + Fel ) = −Fel Ennoel,ind = − (Fnoel n n µn n µn 2 n 1 Enel,ind = − Fel µel,ind 2 n n

(31) (32)

The difference between the induction effect from LR and a state-specific approach is (Note that for ∆E(0)ind i0 no dipoles are induced through the excitation and therefore, no

24


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


work against the electric field occurs and hence, the factor

1 2

is removed):

1 el el,ind Fi µi − F0el µel,ind 0 2 1 el el,ind el el,ind el el,ind + + F µ −Fiel µel,ind F µ − F µ 0 0 0 0 i 0 i 2 1 −2Fiel µ0el,ind + F0el µel,ind + Fiel µel,ind 0 i 2 1 el el,ind el el,ind el el,ind −Fiel µel,ind − F µ + F µ + F µ 0 0 0 0 i i i 2 1 Fiel − F0el µel,ind − µel,ind 0 i 2 1 el,ind ∆Fi0el ∆µi0 2

el,ind ∆E(0)el,ind − ∆Ei0 = −∆Fi0el µel,ind + 0 i0

= = = = =

(33) (34) (35) (36) (37) (38)

This is a result of the symmetry in Eq. (6) so that Fiel µel,ind = F0el µel,ind . Once access 0 i el,ind to ∆Pi0 is available (via standard LR techniques), a correction ∆∆Ei0 can be computed

such that:

el,ind el,ind ∆Ei0 = ∆E(0)el,ind + ∆∆Ei0 i0

(39)

Supporting Information Available Basis set study for full QM cluster computations and all individual excitation energies with cc-pVTZ and (aug)-cc-pVDZ basis sets.

This material is available free of charge via the

Internet at http://pubs.acs.org/.

References (1) Canuto, S., Ed. Solvation effects on molecules and biomolecules; Springer Science + Business Media, 2008. (2) Reichardt, C.; Welton, T. Solvents and solvent effects in organic chemistry; WileyBlackwell, 2010. 25



(3) Milanese, J. M.; Provorse, M. R.; Alameda, E.; Isborn, C. M. Convergence of Computed Aqueous Absorption Spectra with Explicit Quantum Mechanical Solvent. J. Chem. Theory Comput. 2017, 13, 2159–2171. (4) Sneskov, K.; Schwabe, T.; Christiansen, O.; Kongsted, J. Scrutinizing the effects of polarization in QM / MM excited state calculations. Phys. Chem. Chem. Phys. 2011, 13, 18551–18560. (5) Söderhjelm, P.; Husberg, C.; Strambi, A.; Olivucci, M.; Ryde, U. Protein Influence on Electronic Spectra Modeled by Multipoles and Polarizabilities. J. Chem. Theory Comput. 2009, 5, 649–658. (6) Schwabe, T.; Beerepoot, M. T. P.; Olsen, J. M. H.; Kongsted, J. Analysis of computational models for an accurate study of electronic excitations in GFP. Phys. Chem. Chem. Phys. 2015, 17, 2582–2588. (7) Lewandowski, J. R.; Dumez, J.-N. N.; Akbey, U.; Lange, S.; Emsley, L.; Oschkinat, H. Enhanced Resolution and Coherence Lifetimes in the Solid-State NMR Spectroscopy of Perdeuterated Proteins under Ultrafast Magic-Angle Spinning. J. Phys. Chem. Lett. 2011, 2, 2205–2211. (8) Schwabe, T.; Olsen, J. M. H.; Sneskov, K.; Kongsted, J.; Christiansen, O. Solvation Effects on Electronic Transitions: Exploring the Performance of Advanced Solvent Potentials in Polarizable Embedding Calculations. J. Chem. Theory Comput. 2011, 7, 2209–2217. (9) Daday, C.; Curutchet, C.; Sinicropi, A.; Mennucci, B.; Filippi, C. ChromophoreProtein Coupling beyond Nonpolarizable Models: Understanding Absorption in Green Fluorescent Protein. J. Chem. Theory Comput. 2015, 11, 4825–4839. (10) Christiansen, O.; Jørgensen, P.; Hättig, C. Response Functions from Fourier Compo-

26


Page 26 of 33

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


nent Variational Perturbation Theory Applied to a Time-Averaged Quasienergy. Int. J. Quantum Chem. 1998, 68, 1–52. (11) Cammi, R.; Corni, S.; 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, 122, 104513. (12) 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. (13) Schwabe, T. General theory for environmental effects on (vertical) electronic excitation energies. J. Chem. Phys. 2016, 145, 154105. (14) N˚ abo, L. J.; Olsen, J. M. H.; List, N. H.; Solanko, L. M.; W¨ ustner, D.; Kongsted, J. Embedding beyond electrostatics—The role of wave function confinement. J. Chem. Phys. 2016, 145, 104102. (15) Zhu, C.; Dalgarno, A.; Porsev, S. G.; Derevianko, A. Dipole polarizabilities of excited alkali-metal atoms and long-range interactions of ground- and excited-state alkali-metal atoms with helium atoms. Phys. Rev. A 2004, 70, 032722. (16) Buhmann, S. Y. Springer Tracts in Modern Physics; Springer Science + Business Media, 2012; pp 113–147. (17) Suppan, P. Invited review solvatochromic shifts: The influence of the medium on the energy of electronic states. J. Photochem. Photobiol. A 1990, 50, 293–330. (18) Catalán, J.; Valle, J. C. D. A spectroscopic rule from the solvatochromism of aromatic solutes in nonpolar solvents. J. Phys. Chem. B 2014, 118, 5168–5176. 27



(19) Caricato, M. A comparison between state-specific and linear-response formalisms for the calculation of vertical electronic transition energy in solution with the CCSD-PCM method. J. Chem. Phys. 2013, 139, 044116. (20) Caricato, M.; Mennucci, B.; Tomasi, J.; Ingrosso, F.; Cammi, R.; Corni, S.; Scalmani, G. Formation and relaxation of excited states in solution: A new time dependent polarizable continuum model based on time dependent density functional theory. J. Chem. Phys. 2006, 124, 124520. ´ Caprasecca, S.; Lagardère, L.; Lipparini, F.; Piquemal, J. P.; (21) Loco, D.; Polack, E.; Mennucci, B. A QM/MM Approach Using the AMOEBA Polarizable Embedding: From Ground State Energies to Electronic Excitations. J. Chem. Theory Comput. 2016, 12, 3654–3661. (22) Loco, D.; Lagardére, L.; Caprasecca, S.; Lipparini, F.; Mennucci, B.; Piquemal, J.P. Hybrid QM/MM Molecular Dynamics with AMOEBA Polarizable Embedding. J. Chem. Theory Comput. 2017, 13, 4025–4033, PMID: 28759205. (23) Arora, P.; Slipchenko, L. V.; Webb, S. P.; DeFusco, A.; Gordon, M. S. Solvent-Induced Frequency Shifts: Configuration Interaction Singles Combined with the Effective Fragment Potential Method. J. Phys. Chem. A 2010, 114, 6742–6750. (24) Kosenkov, D.; Slipchenko, L. V. Solvent effects on the electronic transitions of pnitroaniline: A QM/EFP study. J. Phys. Chem. A 2011, 115, 392–401. (25) Guido, C. A.; Jacquemin, D.; Adamo, C.; Mennucci, B. Electronic excitations in solution: The interplay between state specific approaches and a time-dependent density functional theory description. J. Chem. Theory Comput. 2015, 11, 5782–5790. (26) Guareschi, R.; Zulfikri, H.; Daday, C.; Floris, F. M.; Amovilli, C.; Mennucci, B.; Filippi, C. Introducing QMC/MMpol: Quantum Monte Carlo in polarizable force fields for excited states. J. Chem. Theory Comput. 2016, 12, 1674–1683. 28


Page 28 of 33

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


(27) Lunkenheimer, B.; Köhn, A. Solvent Effects on Electronically Excited States Using the Conductor- Like Screening Model and the Second-Order Correlated Method. J. Chem. Theory Comput. 2013, 9, 977–994. (28) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and DensityFunctional Theory. Phys. Rev. Lett. 1985, 55, 2471–2474. (29) Seifert, G.; Joswig, J. O. Density-functional tight binding-an approximate densityfunctional theory method. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 456–465. (30) Schwabe, T. Assessing Molecular Dynamics Simulations with Solvatochromism Modeling. J. Phys. Chem. B 2015, 119, 10693–10700. (31) Olsen, J. M.; Aidas, K.; Kongsted, J. Excited states in solution through polarizable embedding. J. Chem. Theory Comput. 2010, 6, 3721–3734. (32) Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. The polarizable embedding coupled cluster method. J. Chem. Phys. 2011, 134, 104108. (33) Schwabe, T.; Sneskov, K.; Haugaard Olsen, J. M.; Kongsted, J.; Christiansen, O.; Hättig, C. PERI-CC2: A polarizable embedded RI-CC2 method. J. Chem. Theory Comput. 2012, 8, 3274–3283. (34) Nørby, M. S.; Vahtras, O.; Norman, P.; Kongsted, J. Assessing frequency-dependent site polarisabilities in linear response polarisable embedding. Molecular Physics 2017, 115, 39–47. (35) List, N. H.; Norman, P.; Kongsted, J.; Jensen, H. J. A. A quantum-mechanical perspective on linear response theory within polarizable embedding. J. Chem. Phys. 2017, 146, 234101. (36) Gagliardi, L.; Lindh, R.; Karlström, G. Local properties of quantum chemical systems

29



: The LoProp approach Local properties of quantum chemical systems : The LoProp approach. J. Chem. Phys. 2004, 121, 4494–4500. (37) Schröder, H.; Schwabe, T. Efficient determination of accurate atomic polarizabilities for polarizeable embedding calculations. J. Comput. Chem. 2016, 37, 2052–2059. (38) Minezawa, N.; Silva, N. D.; Zahariev, F.; Gordon, M. S. Implementation of the analytic energy gradient for the combined time-dependent density functional theory/effective fragment potential method: Application to excited-state molecular dynamics simulations. J. Chem. Phys. 2011, 134, 054111. (39) Furche, F.; Ahlrichs, R. Adiabatic time-dependent density functional methods for excited state properties. J. Chem. Phys. 2002, 117, 7433–7447. (40) Scalmani, G.; Frisch, M. J.; Mennucci, B.; Tomasi, J.; Cammi, R.; Barone, V. Geometries and properties of excited states in the gas phase and in solution : Theory and application of a time-dependent density functional theory polarizable continuum model Geometries and properties of excited states in the gas phase and in solution :. J. Chem. Phys. 2006, 124, 094107. (41) Si, D.; Li, H. Analytic energy gradient in combined time-dependent density functional theory and polarizable force field calculation. J. Chem. Phys. 2010, 133, 144112. (42) Menger, M. F. S. J.; Caprasecca, S.; Mennucci, B. Excited-State Gradients in Polarizable QM/MM Models: An Induced Dipole Formulation. J. Chem. Theory Comput. 2017, 13, 3778–3786. (43) Schwabe, T.; Magnus, J.; Olsen, H.; Sneskov, K.; Kongsted, J.; Christiansen, O. Solvation Effects on Electronic Transitions: Exploring the Performance of Advanced Solvent Potentials in Polarizable Embedding Calculations. J. Chem. Theory Comput. 2011, 7, 2209–2217.

30


Page 30 of 33

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


(44) Becke, A. D. Density Functional Thermochemistry III The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. (45) Stephens, P.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab-Initio Calculation of Vibrational Absorption and Circular-Dichroism Spectra Using Density-Functional Force-Fields. J. Phys. Chem. 1994, 98, 11623–11627. (46) 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. (47) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (48) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297. (49) Klamt, A. Calculation of UV / Vis Spectra in Solution. J. Phys. Chem. 1996, 100, 3349–3353. (50) TURBOMOLE, V7.1, 2016, developed by University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH since 2017;available from http://www.turbomole.com. (51) Furche, F.; Rappoport, D. Comput. Photochem.; 2005; Vol. 16; pp 93–128. (52) Häser, M.; Ahlrichs, R. Improvements on the direct SCF method. J. Comput. Chem. 1989, 10, 104–111. (53) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Physical Review A 1988, 38, 3098–3100.

31



(54) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (55) Becke, A. D. A new mixing of HartreeFock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372. (56) Dierksen, M.; Grimme, S. The vibronic structure of electronic absorption spectra of large molecules: A time-dependent density functional study on the influence of ”Exact” hartree-fock exchange. J. Phys. Chem. A 2004, 108, 10225–10237. (57) Eriksen, J. J.; Sauer, S. P.; Mikkelsen, K. V.; Christiansen, O.; Jensen, H. J. A.; Kongsted, J. Failures of TDDFT in describing the lowest intramolecular charge-transfer excitation in para-nitroaniline. Mol. Phys. 2013, 111, 1235–1248. (58) Treutler, O.;

Ahlrichs, R. Efficient molecular numerical integration schemes.

J. Chem. Phys. 1995, 102, 346–354. (59) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (60) Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to K. J. Chem. Phys. 1992, 97, 2571–2577. (61) Aidas, K. Whirlpool, a QM/MM analysis program, 1.0. 2010. (62) Rappoport, D.; Furche, F. Analytical time-dependent density functional derivative methods within the RI-J approximation, an approach to excited states of large molecules. The Journal of Chemical Physics 2005, 122, 064105. (63) Tkatchenko, A.; Scheffler, M. Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. Phys. Rev. Lett. 2009, 102, 073005.

32


Page 32 of 33

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

33


Corrected Polarizable Embedding: Improving the Induction

Recommend Documents