Polarizable Density Embedding Coupled Cluster Method - Journal of

Subscriber access provided by The University of British Columbia Library

Article

The Polarizable Density Embedding Coupled Cluster Method Dalibor Hršak, Jógvan Magnus Haugaard Olsen, and Jacob Kongsted J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01153 • Publication Date (Web): 14 Feb 2018 Downloaded from http://pubs.acs.org on February 15, 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.

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 31 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

The Polarizable Density Embedding Coupled Cluster Method Dalibor Hrˇsak, Jógvan Magnus Haugaard Olsen, and Jacob Kongsted∗ Department of Physics, Chemistry and Pharmacy, University of Southern Denmark Campusvej 55, 5230 Odense M, Denmark E-mail: [email protected] Phone: +45 65502304 Abstract We present the theory and implementation of the polarizable density embedding (PDE) model in combination with coupled cluster (CC) theory (PDE-CC). This model has been implemented in the Dalton quantum chemistry program by adapting the CC code to the polarizable embedding library (PElib). In the PDE-CC method, the smaller, but chemically important core region is described with a high-level CC method. The environment surrounding the core region can be separated into two levels of description: an inner and an outer region. The effect of the inner region on the core region is described by an embedding potential consisting of a set of fragment densities obtained from calculations on isolated fragments with a quantum-chemistry method such as Hartree–Fock (HF) or Kohn–Sham density functional theory (KS-DFT) and dressed with a set of atom-centered anisotropic dipole-dipole polarizabilities. The outer region consists of distributed atom-centered multipoles and polarizabilities, i.e. in the same way as in the polarizable embedding (PE) model. The PDE-CC method contains embedding potential operators that account for the electrostatic and polarization interactions between the core region and the environment, as well as for non-electrostatic

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

(also known as Pauli and exchange) repulsion. All environmental effects are included through one-electron operators and account very efficiently for the response of the environment due to the change in the electron density of the core region, e.g., upon an electronic transition.

1

Introduction

One of the important topics of theoretical and computational chemistry is an accurate description of spectroscopic properties of large molecular systems, such as proteins or solutesolvent systems. Hybrid quantum mechanics/molecular mechanics (QM/MM) methods 1–5 offer efficient calculations on such systems, where interactions between the chromophore and its environment modulate the molecular properties. The main advantage of such methods is the focused effort for accurate description of the chemically important part, while still retaining an atomistic description of the environment and thereby the capability of describing explicit interactions of the molecule with its nearest environment. Another, less refined, but computationally less expensive environment model is the polarizable continuum model (PCM), which emulates the bulk properties of the surrounding, but lacks the explicit interactions between the two subsystems. 6,7 Use of polarizable force fields enables an advanced and flexible modeling of solvents and biochemical environments. 8–10 One of the models able to describe polarization effects from a discrete environment is the polarizable embedding (PE) model, 11,12 which has been developed in recent years. The PE model has a convenient mathematical structure that allows for addition of simple one-electron effective operators to the vacuum Hamiltonian and thereby accounts for electrostatic interactions, as well as mutual polarization between the quantum region and its environment in an efficient fashion. This is particularly amenable for the description of molecular response properties, where the environment needs to be able to adapt to changes in the electron density of the molecule in the quantum region, e.g., upon an electronic excitation. The environment in the PE model is described by atom-centered mul2


Page 2 of 31

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


tipoles and polarizabilities up to a given order. Such embedding potentials are here denoted as MxPy, where x is an integer denoting the order of multipoles (0: charges, 1: charges and dipoles, etc.) and y denotes the type of polarizabilities (0: no polarizabilities, 1: isotropic polarizabilities, 2: anisotropic polarizabilities, etc.). The multipoles and polarizabilities are calculated for isolated fragment molecules (e.g. solvent molecules or amino acids) using a given quantum chemistry method. The mathematical structure of the PE operators stems from classical electrostatic expressions that describe permanent electrostatics and polarization. 13 In addition to Hartree–Fock (HF) and Kohn–Sham density functional theory (KS-DFT), 11 PE has also been combined with coupled cluster (CC), 14,15 multi-configurational self-consistent field (MCSCF), 16 as well as a number of other theoretical methods 17–19 and response theory approaches. 20,21 These approaches have been shown to be able to model the effects of the environment on various spectroscopic properties in different types of systems and in general achieve good agreement with experiments. Examples of such applications include solvent effects on one- and two-photon absorption spectra, 22 phosphorescence lifetimes of chromophores in solutions, 23 spectroscopic properties of proteins, 24–26 as well as UV spectra of molecules embedded into DNA environments. 27 The promising performance of the PE model calls for further advances in the field, one of which being a more refined description of the surrounding environment in the nearest vicinity of the target molecule. This effort has recently resulted in the development of a three-layered QM/QM/MM model termed polarizable density embedding (PDE). 28 The total system in this model can be separated into three regions; 1) the core region consisting of the target molecule (solute or chromophore) fully described by a quantum-chemistry method, 2) an inner region where each fragment (solvent molecule) is described semi-classically and in isolation, and 3) an outer region consisting of the rest of the environment where each fragment is described classically in the same fashion as in the PE model. The polarization interactions in the inner region retains the same form as in the PE model. The embedding

3



potential that describes the inner region is here denoted as FDPyR, where FD denotes the description of electrostatic interactions with fragment electron densities, y denotes the type of the polarizabilities and R denotes the inclusion of the non-electrostatic repulsion in the model. This latter component is brought about by another layer of complexity needed for an accurate description of explicit interactions in cases where the electron density of the chromophore extends into the environment and overlaps with the electron densities of the environment molecules. 29,30 The PDE model has recently been successfully applied in a case where solvent confinement effect was the main solvent contribution. 31 Until now, the PDE model has only been combined with HF and KS-DFT. In this work, we report the theory and implementation of the PDE model combined with CC methods (PDE-CC). Conventional CC is known to be a very accurate method, but characterized by a steep scaling in computational cost with respect to the size of the system, and this feature prohibits a full CC treatment of larger molecular systems. Here, we will harness both the high accuracy of CC and the ability of PDE to accurately represent intermolecular interactions using relatively simple one-electron interaction operators. A crucial feature of the PDE model is to include polarization effects classically via polarizabilites just as in the PE model. This allows us to use the same density-driven strategy as outlined in Ref. 14. This polarization plays an important role when addressing molecular response-property calculations. Coupled cluster has previously been combined with other embedding type models for polarizable force fields. One such implementation is the combination of equation-of-motion (EOM) CC and the effective fragment potential 32 (EFP) model. 33–36 Embedding using EFP treats the environment in a similar way as PE and the parameters for the embedding potential are also obtained through first-principles calculations on isolated fragments. Equation-ofmotion CCSD has also been combined with a two-level fluctuating charges / polarizable continuum model 37 (EOM-CCSD-FQ-PCM), where polarization is accounted for through the redistribution of the point charges according to the electronegativity equalization principle.

4


Page 4 of 31

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


A multilevel QM/QM analogue of the frozen density embedding 38 (FDE) is the CC-inDFT approach, 39 that has been implemented for the calculation of the molecular response properties. The latter approach bears resemblances to PDE-CC due to the CC-in-DFT formalism, while the difference lies in the treatment of the polarization effects, where CC-inDFT applies the freeze-and-thaw procedure while PDE-CC introduces polarization through classical electrostatics.

2

Theory

2.1

The PDE model

The first step in the derivation of the PDE-CC method is to introduce the relevant interactions between the core part and its environment. The contributions specific for the PDE model are detailed in the original work (at the SCF level of theory), 28 but will be described here along with contributions that are shared with the PE model. For an in-depth description of PE-HF/DFT or PE-CC, we refer the reader to Refs. 11, 12, and 14. The polarizable environment is divided into two layers that surround the central molecule: an inner region that consists of the fragments in the nearest vicinity and an outer region that consists of the rest of the environment. The embedding potential of the inner region is defined by a set of fragment electronic densities, point nuclei, and a set of atom-centered polarizabilities, all obtained through calculations on isolated fragments. The embedding potential of the outer region is defined by fragment-based atom-centered multipoles and polarizabilities. The expansion points used for the atom-centered multipoles and polarizabilities will be referred to as the classical sites. The interactions with the environment are described through one-electron interaction operators accounting for electrostatic interaction of the outer-region multipoles, electrostatic interaction of the inner-region fragment densities, non-electrostatic repulsion and polarization terms. The operator describing the electrostatic interaction of the core wave function 5



Page 6 of 31

with the outer-region multipoles, Vˆ mul , is equivalent to the electrostatic interaction operator in the PE model. This operator has the form

Vˆ mul =

S X Ks X (−1)|k| (k) X (k) ˆ Ms vs,pq Epq k! pq s=1

(1)

|k|=0

where the summation over s runs over the S sites in the classical region, k is a multi-index, 40 (k)

Ks is the truncation level of the multipole expansion, and Ms

is a component of a multipole

(k)

moment located at site s, and vs,pq is a one-electron interaction integral which involves partial derivatives of the

1 r

operator where r is the distance between an electron and a multipole

site. The superscript (k) notation is used to indicate order and Cartesian component. For example, in the case of multipoles, k = (0, 0, 0) specifies a monopole (i.e. charge), k = (0, 1, 1) corresponds to the yz-component of a quadrupole, and so on. The electrostatic interactions between the core and the fragments in the inner region are described by the electrostatic fragment density operator

Vˆ fd = −

Na X X

a ˆ vpq Epq

+

Na X X X

a a ˆ vpq,ij Dij Epq

(2)

a=1 pq ij∈a

a=1 pq

a is a nuclear where the summations over a run over the Na fragments in the inner region, vpq

attraction integral involving orbitals p and q belonging to the core and the nuclei in the a inner-region fragment a, vpq,ij is an inter-fragment two-electron repulsion integral involving

orbitals p and q belonging to the core region and orbitals i and j belonging to inner-region a fragment a, and Dij is an element of the density matrix of the inner-region fragment a.

This set of operators describe the intermolecular interactions under the assumption that the wave functions of the two interacting fragments do not overlap. To model non-electrostatic exchange repulsion, another operator is introduced that is based on the Huzinaga-Cantu

6


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


projection operator method 41 and has the following form in the MO basis

Vˆ rep = −

Na N occ X X X a=1 i∈a

εi Spi Siq Eˆpq

(3)

pq

Here, Spi and Siq are elements of the intermolecular overlap matrix involving the occupied orbitals of the inner-region fragments and εi is a fragment MO energy. Transformation of this operator from the MO to the AO basis gives rise to the presence of an energy-weighted density matrix as given in Ref. 28. The projection operator method simulates orthogonality between the wave functions of the inner-region fragments and the core-region wave function. In order to derive the equations for optimization of the cluster amplitudes and multipliers, we proceed by adding the above one-electron operators to the vacuum Hamiltonian. This results in a time-independent PDE-CC Lagrangian which given as

L

PDE−CC

S X S E 1X D ˆ vac ˆ mul ˆ fd ˆ rep ¯ Es (D)| Bss0 Es0 (D) (t, t) = Λ H + V + V + V CC − 2 s=1 s0 =1

(4)

based on the usual definitions of the CC states. 42 The last part of Eq. (4) is the energy associated with polarization of the environment. The Bss0 matrix is a 3 × 3 subblock of the classical linear response matrix (also known as the relay matrix) in which the diagonal subblocks (s = s0 ) contain the inverse of the polarizabilities and the off-diagonal subblocks (s 6= s0 ) contain dipole–dipole interaction tensors. The electric field at a site s, i.e. Es (D), contains the field from the multipoles and fragment densities in the environment, and in addition the fields from the nuclei and electrons in the core region which thus depends on the CC density. By enforcing stationarity of the Lagrangian with respect to both amplitudes and multipliers, we obtain the zeroth-order CC amplitudes and multipliers. The resulting PDE-CC optimization equations exhibit the exact same structure as the corresponding PE-CC opti-

7



Page 8 of 31

mization equations ∂LPDE-CC D ˆ vac ˆ emb E = µ ¯ i H + V CC = 0 ∂ t¯µi i E ∂LPDE-CC D h ˆ vac ˆ emb = Λ H + V , τˆµi CC = 0 ∂tµi

(5a) (5b)

The embedding operator Vˆ emb is defined in a density-driven polarization fashion through the equation Vˆ emb = Vˆ mul + Vˆ fd + Vˆ rep −

S X

| ˆ el µind s (D) Es ,

(6)

s=1

ˆ el is the electric-field operator, due to the electrons at site s, and µind is the induced where E s s dipole at the site s. As for the case of PE-CC, the PDE-CC amplitudes and multipliers equations are coupled because the embedding operator contains both sets of parameters (amplitudes and multipliers), which calls for an iterative optimization procedure. This problem is solved as outlined in Ref. 43, i.e. by performing the so-called inner and outer iterations and rely on a full solution of the coupled set of equations.

2.2

Polarizable Density Embedding CC response theory

To describe molecular response properties and excited states, we apply response theory as defined in Ref. 44 and formulated through Fourier component perturbation theory applied to the quasi-energy in Ref. 42. We refer to these papers for more details on the formalism. Frequency-dependent molecular response properties are calculated from response functions, defined as the derivatives of the time-averaged quasi-energy Lagrangian with respect to the perturbation strength parameters. For this work, the linear and quadratic response functions are considered and their final expressions after applying the chain-rule for partial derivatives

8


Page 9 of 31 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 the 2n + 1 rule 45 for the amplitudes and multipliers are

hhX; Y iiωY hhX; Y, ZiiωY ,ωZ

1 ˆ ±ω ˆ XY X Y 1 ¯X 1 X Y Y = C P (7a) η t (ωY ) + Ft (ωX )t (ωY ) − t (ωX )¯t (ωY )J 2 2 2 1 1 ˆ ±ω ˆ XY 1 X Y F t (ωY ) + GtX (ωX )tY (ωY ) = C P 2 2 6 1 ¯X 1 ¯X Y X Y Y ¯ ¯ + t (ωX )A + t (ωX )Bt (ωY ) + t (ωX )t (ωY )K tZ (ωZ ) (7b) 2 2

where ωX , ωY and ωZ are the frequencies associated with the external perturbations, Cˆ ±ω is the frequency symmetrizer, and Pˆ is the permutation operator. The response matrices and vectors entering into the expression for the linear and quadratic response functions are defined in Table 1. Enforcing the stationarity of the second-order time-averaged quasi-energy Table 1: General expressions for the matrices and vectors in the PDE-CC response theory. Quantity ηµXi Fµi νj Jµi νj

Vacuum contribution D h i E ˆ Λ X, τˆµi CC D hh i i E ˆ Λ H ˆµi , τˆνj CC 0, τ

PDE contribution

Bµi νj σk

D h i E ˆ µ ¯i H ˆνj CC 0, τ D E µ ¯i Yˆ CC D hh i i E ˆ τˆµi , τˆνj CC Λ X, D hhh i i i E ˆ 0 , τˆµ , τˆν , τˆσ CC Λ H i j k D h i E µ ¯i Yˆ , τˆνj CC D hh i i E ˆ µ ¯i H ˆνj , τˆσk CC 0, τ

Kµi νj σk

-

Aµi νj ξµYi FµXi νj Gµi νj σk AYµi νj

i E h D hh i i Λ Vˆ emb , τˆµi , τˆνj + 12 Pˆ µi νj Vˆ emb,µi , τˆνj CC E D µ ¯i νj Vˆ emb CC E D h i µ ¯i Vˆ emb , τˆνj + Vˆ emb,νj CC D hhh i hh i i E i i Λ Vˆ emb , τˆµi , τˆνj , τˆσk + 12 Pˆ µi νj σk Vˆ emb,µi , τˆνj , τˆσk CC D hh i i h i E µ ¯i Vˆ emb , τˆνj , τˆσk + Pˆ νj σk Vˆ emb,νj , τˆσk + 21 Vˆ emb,νj σk CC D h i E Pˆ µi νj µ ¯i νj Vˆ emb , τˆσk CC

Lagrangian with respect to the first-order amplitudes and multipliers results in the first-order

9



Page 10 of 31

CC equations, which enable us to obtain the first-order parameters

ξ Y + (A − ωY I)tY (ωY ) + ¯tY (ωY )J = 0

(8a)

η Y + FtY (ωY ) + ¯tY (ωY )(A + ωY I) = 0

(8b)

These equations show the same structure as for the case of PE-CC and, exactly as for the case of PE-CC, these equations are coupled in the amplitudes and multipliers and are therefore solved based on a similar strategy as the one outlined for the ground state parameters. As seen from Table 1, the form of the PDE response matrices are similar in structure to the form of the PE response matrices, with the only difference being the new contributions to the effective environment operator entering in the definition of the Vˆ emb operator. The single and double partial derivatives of the effective environment operator with respect to the νj amplitudes (Vˆ emb,µi and Vˆ emb,µi νj ) and with respect to the multipliers ( Vˆ emb ) have exactly

the same form as in the PE-CC theory because only the polarization part depends on the cluster amplitudes and Lagrange multipliers, while the new contributions specific for PDE do not. The exact forms of these operators can be found in Ref. 14.

3

Implementation Aspects

The PDE-CC method has been implemented in a development version of Dalton 2016 46 by coupling the existing CC routines to the polarizable embedding library, PElib. 47 The available CC wave-function methods in the new PDE-CC implementation are CC2, 48 CCSD and CCSDR(3). 49 In addition to the PDE model, PElib provides a number of new features and functionalities not present in the previous PE-CC implementation, such as for example inclusion of multipole moments up to fifth order (M5). In order to efficiently solve the response equations, the matrices are transformed with left or right trial vectors. 50 The expressions for

10


Page 11 of 31 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 transformed matrices are listed in the following equation APDE c =

h i E i E X D h ˆ emb,C µ ¯i Vˆ emb , τˆνj + Vˆ emb,νj CC cνj = µ ¯i Vˆ emb , C CC + ξ V

(9a)

D h i E X D h i E C ¯ Vˆ emb , τˆν CC + η Vˆ emb cµi µ ¯i Vˆ emb , τˆνj + Vˆ emb,νj CC = C j

(9b)

D ν E C ˆ emb cµi µ ¯i j Vˆ emb CC = ξ V

(9c)

XD µi νj

cAPDE =

X µi νj

cJ =

X µi νj

FPDE c =

µi

νj

i E i i E i i h X D hh X D hh ˆ emb,C Λ Vˆ emb , τˆµi , C CC + η V Λ Vˆ emb , τˆµi , τˆνj + Vˆ emb,µi , τˆνj CC cνj =

µi νj

GPDE bc =

(9d)

µi

i i E i i hh i i hh X D hh Λ Vˆ emb,µi , τˆνj , τˆσk + Vˆ emb,νj , τˆµi , τˆσk + Vˆ emb,σk , τˆµi , τˆνj CC bµi cνj

µi νj σk

ˆ emb,C

ˆ emb,B

ˆ emb,BC

c + ηV b + FV = FV i h i h i i E X D hh σ ¯k Vˆ emb , τˆµi , τˆνj + Vˆ emb,µi , τˆνj + Vˆ emb,νj , τˆµi + Vˆ emb,µi νj CC bµi cνj BPDE bc =

(9e)

µi νj σk

=

XD σk

bcK =

hh E i i ˆ emb,B ˆ emb,C ˆ emb,BC c + AV b + ξV σ ¯k Vˆ emb , B , C CC + AV

X

bµi cνj

µi νj σk

D

h E i E D h µ i B ˆ emb C ˆ emb ν = Cη V + Bη V µ ¯i j Vˆ emb , τˆσk CC + ν¯j i Vˆ emb , τˆσk CC

(9f) (9g)

Matrices J and K represent relatively small contributions and they can be neglected when calculating the residues of the response functions needed for the transition moments. This approximation facilitates the calculations and produces only small errors as previously observed. 51,52 The matrices and vectors in the above equation have been constructed in the density-driven response polarization fashion, in which first the transformed densities (DB , DC , B D, C D and DBC ) are calculated and used to construct the effective environment operators. The difference from the previous PE-CC implementation is in the response calculation, whereby PElib allows to pass multiple densities simultaneously for the construction of the Vˆ emb operators instead of passing only one density at a time. This has the potential to save computational time when calculating the response as it avoids repeated calculations of electric-field integrals.

11



4

Calculations

4.1

Formamide in water

To illustrate the performance of the PDE-CC method – and compare the results to PE-CC – we first consider the case of formamide in water solution. The electronic spectrum of formamide is characterized by a number of Rydberg transitions in the region around the most intense π → π ∗ transition, 53–56 and such Rydberg states are destabilized by the solvent through non-electrostatic repulsion. 57,58 The spectra of aqueous formamide have already been a target of PE-CC studies, 14,59,60 but the aforementioned destabilization of the Rydberg states has not been clearly analyzed. For our analysis, 120 configurations (snapshots) of aqueous formamide were extracted from an MD simulation as previously reported. 59 Each snapshot includes approximately 345 water molecules corresponding to a buffer size of 12 ˚ A from the solute. The water solvent was described with one of two types of embedding potentials; for PDE-CC the FDP2R potential was used for all solvent molecules and for PE-CC the potential was M2P2. 61 Local properties (atom-centered multipole moments and polarizabilities) were calculated based on B3LYP 62–65 /aug-cc-pVDZ. 66–69 Before the local property calculations the basis set has been recontracted to atomic natural orbital (ANO) form. The local property calculations were performed in Molcas 70 using the LoProp formalism, 71 whereas the fragment densities were calculated in the Dalton program 46 also with B3LYP/aug-cc-pVDZ. The non-electrostatic repulsion was scaled by a factor of 0.7 to remedy the overestimation of short-range effects that arises from the fact that the fragment densities are optimized for isolated solvent molecules. 72 For the vacuum calculations the formamide molecule was optimized at the level of B3LYP/aug-cc-pVTZ using the Gaussian 09 program. 73 The 10 lowest singlet excited states were calculated with either CC2, CCSD or CCSDR(3) in combination with PDE or PE and with the aug-cc-pVDZ basis set. From experiments in vacuum a series of transitions have been observed consisting of one n → π ∗ transition,

12


Page 12 of 31

0.4 0.3 0.2 0.1

CCSD

0.006 0.004 0.002

S1 S2 S3 S4 S4 S6 S7 S8 S9 S10

0.006 0.004 0.002 5.8

6.2

6.6

7.0

0.3

0.006 0.004 0.002 5.8

0.006 0.004 0.002

0.4

CCSDR(3)

6.2

6.6

7.0

0.006 0.004 0.002 5.8

6.2

6.6

7.0

5.8

6.2

6.6

7.0

PDE

CC2 oscillator strength oscillator strength

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


5.8

6.2

6.6

7.0

PE

Page 13 of 31

0.006 0.004 0.002 5.8

6.2

6.6

7.0

0.2 0.1

6.0

7.0

8.0

9.0

E / eV

10.0

11.0

6.0

7.0

8.0

9.0

10.0

11.0

E / eV

6.0

7.0

8.0

9.0

10.0

11.0

E / eV

Figure 1: Distribution of the 10 lowest-lying singlet excited states across the 120 snapshots of formamide solvated in water, calculated with PDE-CC or PE-CC methods. The insets show in more detail the region where the lowest singlet excited state (S1 ) can be found (the units are same as in the main plots). one Rydberg transition, an intense π → π ∗ transition, followed by a Rydberg and a Q-type band. 53,54,57 Here, the calculated electronic spectrum of formamide in isolation is characterized by a weak n → π ∗ transition at ∼5.75 eV for the three CC models, followed by a series of weak Rydberg-type excitations, an intense π → π ∗ transition, and another set of weak transition. For CC2 the intense π → π ∗ transition corresponds to the seventh excited state located at 7.62 eV. For CCSD and CCSDR(3) it is the sixth excited state located at 7.73 eV and 7.66 eV, respectively. The distribution – based on the 120 configurations – of the excited states of formamide in aqueous solution for the given embedding models are shown in Fig. 1. The differences between PDE and PE for all three CC methods are clearly seen. The PE-CC results exhibit a substantial mixture of Rydberg and valence excited states, whereas for PDE-CC valence states and Rydberg states are clearly separated, where all of the configurations have the intense transition with clear π → π ∗ character as the second excited state, and a weak transition with clear n → π ∗ character as the lowest excited state. We have further investigated the non-electrostatic repulsion effect on Rydberg states by changing the repulsion factor. The stick spectra are shown in Fig. S1 in the Supplementary Information. We see that the positions and intensities of n → π ∗ and π → π ∗ are not visibly affected by 13



either reduction or increase in non-electrostatic repulsion. However, the Rydberg states are more affected. For the low repulsion factor (0.5), the mixing with the valence states as seen in PE-CC start to occur (especially for CC2), and for the high repulsion factor (0.9), the Rydberg states are shifted towards higher energies than for the optimal repulsion factor of 0.7. This confirms that the Rydberg states are strongly affected by non-electrostatic repulsion, which is expected because these states are diffuse in nature.

4.2

Acrolein in water

The next example is the case of acrolein in a water solution. These calculations are based on 120 configurations produced by an MD simulation already reported in a previous publication. 74 The embedding potentials were obtained in the same fashion as for the case of formamide and they can be found in Ref. 75. The target molecular property was the UV/vis spectrum, or more precisely, the dependence of the excitation energies and solvent shifts of n → π ∗ and π → π ∗ transitions on the embedding model. The acrolein molecule was solvated in a sphere consisting of approximately 330 water molecules reaching up to 12.0 ˚ A from the solute. Absorption spectra of acrolein in water were based on the 10 lowest singlet excited states, with either PDE or PE in combination with CCSDR(3), CCSD and CC2 methods. The spectra were convoluted using a Lorentzian line shape function with full-width-at-half-maximum (FWHM) set at 0.2 eV (for further details on how to produce UV/vis absorption spectra see Ref. 22) and averaged over the snapshots. The resulting spectra can be seen in Fig. 2 The intense peaks located at ∼6.5 eV for the different methods correspond to the π → π ∗ electronic excitation. While there are no significant differences between the PDE- or PEbased results in that part of the spectrum, this is not the case for the higher transitions. For PDE-CC, the peaks corresponding to the higher excited states are well separated from the intense peak and located at ∼8.5 – 9 eV, while for PE-CC the higher transitions lie closer to the intense peak, namely at ∼7.5 – ∼8.5 eV. These higher transitions in aqueous acrolein 14


Page 14 of 31

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


10-4 × ε / (M-1cm-1)

Page 15 of 31

PDE-CCSDR(3) PE-CCSDR(3) PDE-CCSD PE-CCSD PDE-CC2 PE-CC2

28 26 24 22 20 18 16 14 12 10 8 6 4 2 5.5

6

6.5

7

7.5

8

8.5

9

9.5

10

E / eV Figure 2: Absorption spectra of a water-solvated acrolein molecule described by the PDE or PE model in combination with CCSDR(3), CCSD or CC2.

15



Page 16 of 31

are not readily available with the experiment due to the overlap with the solvent excitations, so it is not clear which results are more accurate in this case. However, these excited states show similar behavior as the Rydberg states in formamide. The results for the averaged excitation energies (based directly on the values for the calculated transition energies and not extracted from the positions of the peaks in the absorption spectrum) and solvatochromic shifts of the two lowest transitions, as well as the comparison with the experiments are given in Table 2. The given standard errors of the mean √ (SEM) are calculated as σ/ n, where σ is the variance of the excitation energies through the snapshots and n is the number of snapshots. Excitation energies as a function of snapshot are plotted in Fig. S3. The difference in the performance between PDE-CC and PE-CC is Table 2: Average excitation energies (∆E) and solvatochromic shifts (∆∆E) given in eV of the two lowest transitions in aqueous acrolein, calculated with PDE or PE in combination with CCSDR(3), CCSD or CC2. Experimental results from Ref. 74 are given for orientation. n → π∗ ∆E ∆∆E PDE 4.184 ± 0.013 PE 4.212 ± 0.016 PDE 4.207 ± 0.013 CCSD PE 4.247 ± 0.016 PDE 4.119 ± 0.013 CCSDR(3) PE 4.161 ± 0.016 experiment 3.94 CC2

0.285 0.313 0.275 0.315 0.280 0.323 0.25

π → π∗ ∆E ∆∆E 6.394 ± 0.013 6.248 ± 0.012 6.537 ± 0.013 6.427 ± 0.014 6.375 ± 0.012 6.272 ± 0.014 5.90

-0.303 -0.449 -0.698 -0.808 -0.581 -0.683 -0.52

in general small compared to other errors involved, when describing the excitation energies of both n → π ∗ and π → π ∗ transitions. It has previously been observed 28,72 that for some systems PE performs well in case of valence excitations based on cancellation of errors due to lack of non-electrostatic repulsion causing charge-leaks into the environment. PDE, on the contrary, provides a physically more correct picture of the short-range interactions, but this may not necessarily be reflected in a significantly better agreement with the experimental data. When comparing calculated solvent shifts with experimental solvent shifts, we see, however, that PDE-CC performs marginally better for all three CC models, except in case 16


Page 17 of 31 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 π → π ∗ transition calculated with CC2.

4.3

4-nitrophenolate in water

The last system considered is the anionic 4-nitrophenolate (p-nitrophenolate or PNP− ) molecule solvated in water. We used 120 configurations consisting of one PNP− molecule, one Na+ ion, and 400 water molecules, that were extracted from a QM/MM Born-Oppenheimer molecular dynamics (BOMD) trajectory. 76 Table 3: Average excitation energies (in eV) of the S1 excited state of aqueous PNP− calculated with PDE or PE in combination with CCSDR(3), CCSD or CC2. CC2

CCSD

CCSDR(3)

PDE 2.916 ± 0.018 3.147 ± 0.019 2.996 ± 0.019 PE 2.851 ± 0.020 3.089 ± 0.021 2.946 ± 0.022

The FDP2R and M2P2 embedding potentials 77 for the PDE- and PE-CC calculations were parametrized in the same manner as in the case of formamide and acrolein, but this time using the 6-31+G* basis set 78–81 also recontracted to ANO-type basis representation as required by the LoProp approach. The same scaling factor of 0.7 was included for the non-electrostatic repulsion. We calculated the excitation energies and oscillator strengths of the two lowest singlet excited states with either PDE or PE in combination with CCSDR(3), CCSD and CC2 methods using the 6-31+G* basis set. The intense band in the spectrum (Fig. 3) corresponds to the HOMO → LUMO transition occurring upon excitation to the lowest singlet excited state (S1 ). The corresponding experimental peak for PNP− in water has been observed at 3.09 eV (i.e. 401 nm). 82 The peaks for PDE are located at somewhat higher energies than the PE peaks and closer to experimental peaks, however the differences between the embedding models are most likely not as significant as errors associated with other approximations involved in the overall computational procedure such as, e.g., the quality of the geometrical parameters of the solute and configurational sampling. We also observe that there is no systematic convergence towards 17



the experimental value when including more electron correlation, i.e. the PDE-CCSD result is very close to the experimental result, while PDE-CCSDR(3) is closer to PDE-CC2.

30

PDE-CCSDR(3) PE-CCSDR(3) PDE-CCSD PE-CCSD PDE-CC2 PE-CC2

25

10-4 × ε / (M-1cm-1)

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

Page 18 of 31

20

15

10

5

2.5

2.6

2.7

2.8

2.9

3.0

3.1

3.2

3.3

3.4

3.5

3.6

E / eV Figure 3: Absorption spectra of a water solvated 4-nitrophenolate described using CCSDR(3), CCSD and CC2 in combination with either PDE or PE model. The averaged excitation energies are given along with SEM in Table 3, and excitation energies as function of snapshot is plotted in Fig. S3. The lower excitation energies from PE-CC can be explained by the lack of non-electrostatic repulsion in PE-CC which leads to over-stabilization of the excited state. 83 Because of the more diffuse character of the excited state it will be stabilized more than the ground state, thus, leading to a lower excitation energy. It is expected that such phenomena are more important for anionic solutes, where a more significant contribution of the diffuse functions can be found, which makes them prone to electron spill-out effects.

18


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


5

Conclusion

We have presented the theory and implementation of the PDE-CC method designed for calculation of molecular response properties of molecules in solutions. Specifically, the PDE model is combined with the CC2, CCSD, and CCSDR(3) methods. The PDE model is an extension of the previous PE model. In PDE the electrostatic effect of each solvent molecule on the core quantum part is modeled through the electron density and point nuclei in contrast to PE which relies on a multipole expansion. This provides a more accurate description of the electrostatic interactions between the solute and the solvent molecules in its nearest environment. Short-range non-electrostatic repulsion effects are described through a projection operator, which is not present in the PE model, and which effectively prevents electron spill-out and constitutes a significant improvement for cases where confinement of the solute wave function by the solvent plays an important role. Another advantage of PDE is that it retains the same efficient procedure for inclusion of environment polarization effects as the PE model. The performance of the PDE-CC method was exemplified on three different molecules, namely aqueous solutions of formamide, acrolein and 4-nitrophenolate, and compared to the corresponding PE-CC results. The case of formamide has clearly shown the advantages of PDE-CC over PE-CC where we have demonstrated that non-electrostatic repulsion plays a crucial role for a correct description of the aqueous UV spectrum of formamide. For acrolein, PDE-CC is somewhat better at predicting the solvent shifts of the n → π ∗ and π → π ∗ transitions when comparing to the experiment. The overall spectrum of acrolein is not heavily affected by non-electrostatic repulsion in the commonly observed part (where the first two excited states are located), but the higher transitions are more strongly affected. The PDE-CC method predicts higher excitation energies for the π → π ∗ transition in 4nitrophenolate compared to PE-CC which is explained by the difference in the stabilization of the involved states. The PE model is prone to over-stabilization effects due to the lack of non-electrostatic repulsion which plays an important role especially for anionic solutes. 19



In conclusion, we find PDE-CC to be a valuable extension of the PE-CC model that introduces significant refinements in the description of especially short-range interactions in systems with a CC core and a polarizable environment.

Acknowledgement The authors thank Erik Rosendahl Kjellgren for providing the structures of 4-nitrophenolate in water solution. Computation/simulation for the work described in this paper was supported by the DeIC National HPC Centre, SDU. J. M. H. O. acknowledges financial support from the Carlsberg Foundation (grant id. CF15-0823). J.K. acknowledges financial support from the Danish e-Infrastructure Cooperation (DeIC), the Villum Foundation, and the Danish Council for Independent Research (DFF).

Supporting Information Available This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227–249. (2) Lin, H.; Truhlar, D. G. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2006, 117, 185–199. (3) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. 2009, 48, 1198–1229.

20


Page 20 of 31

Page 21 of 31 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) Tu, Y.; Laaksonen, A. Implementing Quantum Mechanics into Molecular Mechanics– Combined QM/MM Modeling Methods. Adv. Quantum Chem. 2010, 59, 1–15. (5) Brunk, E.; Röthlisberger, U. Mixed Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217–6263. (6) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (7) Kongsted, J.; Mennucci, B. How to Model Solvent Effects on Molecular Properties Using Quantum Chemistry? Insights from Polarizable Discrete or Continuum Solvation Models. J. Phys. Chem. A 2007, 111, 9890–9900. (8) Vesely, F. J. N-particle dynamics of polarizable Stockmayer-type molecules. J. Comp. Phys. 1977, 24, 361–371. (9) Halgren, T. A.; Damm, W. Polarizable force fields. Curr. Opin. Struct. Biol. 2001, 11, 236–242. (10) Baker, C. M. Polarizable force fields for molecular dynamics simulations of biomolecules. WIREs Comput. Mol. Sci. 2015, 5, 241–254. (11) Olsen, J. M. H.; Aidas, K.; Kongsted, J. Excited States in Solution through Polarizable Embedding. J. Chem. Theory Comput. 2010, 6, 3721–3734. (12) Olsen, J. M. H.; Kongsted, J. Molecular Properties through Polarizable Embedding. Adv. Quantum Chem. 2011, 107–143. (13) Stone, A. J. The Theory of Intermolecular Forces, 1st ed.; Oxford University Press, 1997. (14) Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. The polarizable embedding coupled cluster method. J. Chem. Phys. 2011, 134, 104108. 21



(15) Schwabe, T.; Sneskov, K.; Olsen, J. M. H.; Kongsted, J.; Christiansen, O.; Hättig, C. PERI–CC2: A Polarizable Embedded RI-CC2 Method. J. Chem. Theory Comput. 2012, 8, 3274–3283. (16) Hedeg˚ ard, E. D.; List, N. H.; Jensen, H. J. Aa..; Kongsted, J. The multi-configuration self-consistent field method within a polarizable embedded framework. J. Chem. Phys. 2013, 139, 044101–1. (17) Eriksen, J. J.; Sauer, S. P. A.; Mikkelsen, K. V.; Jensen, H. J. Aa..; Kongsted, J. On the importance of excited state dynamic response electron correlation in polarizable embedding methods. J. Comput. Chem. 2012, 33, 2012–2022. (18) Hedeg˚ ard, E. D.; Reiher, M. Polarizable Embedding Density Matrix Renormalization Group. J. Chem. Theory Comput. 2016, 12, 4242–4253. (19) Hedeg˚ ard, E. D.; Olsen, J. M. H.; Knecht, S.; Kongsted, J.; Jensen, H. J. Aa.. Polarizable embedding with a multiconfiguration short-range density functional theory linear response method. J. Chem. Phys. 2015, 142, 114113. (20) Pedersen, M. N.; Hedeg˚ ard, E. D.; Olsen, J. M. H.; Kauczor, J.; Norman, P.; Kongsted, J. Damped Response Theory in Combination with Polarizable Environments: The Polarizable Embedding Complex Polarization Propagator Method. J. Chem. Theory Comput. 2014, 10, 1164–1171. (21) Steindal, A. H.; Beerepoot, M. T. P.; Ringholm, M.; List, N. H.; Ruud, K.; Kongsted, J.; Olsen, J. M. H. Open-ended response theory with polarizable embedding: multiphoton absorption in biomolecular systems. Phys. Chem. Chem. Phys. 2016, 18, 28339–28352. (22) Hrˇsak, D.; Holmegaard, L.; Poulsen, A. S.; List, N. H.; Kongsted, J.; Denofrio, M. P.; Erra-Balsells, R.; Cabrerizo, F. M.; Christiansen, O.; Ogilby, P. R. Experimental and computational study of solvent effects on one- and two-photon absorption spectra of chlorinated harmines. Phys. Chem. Chem. Phys. 2015, 17, 12090–12099. 22


Page 22 of 31

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


(23) Krause, K.; Bauer, M.; Klopper, W. Approaching Phosphorescence Lifetimes in Solution: The Two-Component Polarizable-Embedding Approximate Coupled-Cluster Method. J. Chem. Theory Comput. 2016, 12, 2853–2860. (24) Beerepoot, M. T. P.; Steindal, A. H.; Kongsted, J.; Brandsdal, B. O.; Frediani, L.; Ruud, K.; Olsen, J. M. H. A polarizable embedding DFT study of one-photon absorption in fluorescent proteins. Phys. Chem. Chem. Phys. 2013, 15, 4735–4743. (25) Sneskov, K.; Olsen, J. M. H.; Schwabe, T.; Hättig, C.; Christiansen, O.; Kongsted, J. Computational screening of one- and two-photon spectrally tuned channelrhodopsin mutants. Phys. Chem. Chem. Phys. 2013, 15, 7567–7576. (26) List, N. H.; Pimenta, F. M.; Holmegaard, L.; Jensen, R. L.; Etzerodt, M.; Schwabe, T.; Kongsted, J.; Ogilby, P. R.; Christiansen, O. Effect of chromophore encapsulation on linear and nonlinear optical properties: the case of ”miniSOG”, a protein-encased flavin. Phys. Chem. Chem. Phys. 2014, 16, 9950–9959. (27) Nørby, M. S.; Steinmann, C.; Olsen, J. M. H.; Li, H.; Kongsted, J. Computational Approach for Studying Optical Properties of DNA Systems in Solution. J. Chem. Theory Comput. 2016, 12, 5050–5057. (28) Olsen, J. M. H.; Steinmann, C.; Ruud, K.; Kongsted, J. Polarizable Density Embedding: A New QM/QM/MM-Based Computational Strategy. J. Phys. Chem. A 2015, 119, 5344–5355. ´ (29) Surján, P.; Angy´ an, J. G. The reliability of the point charge model representing intermolecular effects in ab initio calculations. Chem. Phys. Lett. 1994, 225, 258–264. (30) Fradelos, G.; Wesolowski, T. A. Importance of the Intermolecular Pauli Repulsion in Embedding Calculations for Molecular Properties: The Case of Excitation Energies for a Chromophore in Hydrogen-Bonded Environments. J. Phys. Chem. A 2011, 115, 10018–10026. 23



(31) 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. (32) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The Effective Fragment Potential Method: A QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105, 293– 307. (33) Slipchenko, L. V. Solvation of the Excited States of Chromophores in Polarizable Environment: Orbital Relaxation versus Polarization. J. Phys. Chem. A 2010, 114, 8824– 8830. (34) Ghosh, D. Perturbative approximation to hybrid equation of motion coupled cluster/effective fragment potential method. J. Chem. Phys. 2014, 140, 094101. (35) Ghosh, D. Hybrid Equation-of-Motion Coupled-Cluster/Effective Fragment Potential Method: A Route toward Understanding Photoprocesses in the Condensed Phase. J. Phys. Chem. A 2017, 121, 741–752. (36) 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. (37) Caricato, M.; Lipparini, F.; Scalmani, G.; Cappelli, C.; Barone, V. Vertical Electronic Excitations in Solution with the EOM-CCSD Method Combined with a Polarizable Explicit/Implicit Solvent Model. J. Chem. Theory Comput. 2013, 9, 3035–3042. (38) Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050–8053. (39) Höfener, S.; Gomes, A. S. P.; Visscher, L. Solvatochromic shifts from coupled-cluster theory embedded in density functional theory. J. Chem. Phys. 2013, 139, 104106. 24


Page 24 of 31

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


(40) Saint Raymond, X. Elementary Introduction to the Theory of Pseudodifferential Operators; CRC Press, 1991; Vol. 3. (41) Huzinaga, S.; Cantu, A. A. Theory of Separability of Many-Electron Systems. J. Chem. Phys. 1971, 55, 5543–5549. (42) 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. (43) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Coupled Cluster/Molecular Mechanics Method: Implementation and Application to Liquid Water. J. Phys. Chem. A 2003, 107, 2578–2588. (44) Olsen, J.; Jørgensen, P. Linear and nonlinear response functions for an exact state and for an MCSCF state. J. Chem. Phys. 1985, 82, 3235–3264. (45) Handy, N. C.; Schaefer III, H. F. On the evaluation of analytic energy derivatives for correlated wave functions. J. Chem. Phys. 1984, 81, 5031–5033. (46) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Hst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jans´ık, B.; Jensen, H. J. Aa..; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; ˇ Ruden, T. A.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkeviˇcius, Z.; Ruud, K.; Rybkin, V. V.; Salek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; 25



Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; SylvesterHvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; ˚ Agren, H. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269–284. (47) Olsen, J. M. H. PElib: The Polarizable Embedding library (version 1.0.8). 2014; https://gitlab.com/pe-software/pelib-public. (48) Christiansen, O.; Koch, H.; Jørgensen, P. The second-order approximate coupled cluster singles and doubles model CC2. Chem. Phys. Lett. 1995, 243, 409–418. (49) Christiansen, O.; Koch, H.; Jørgensen, P. Perturbative triple excitation corrections to coupled cluster singles and doubles excitation energies. J. Chem. Phys. 1996, 105, 1451–1459. (50) Olsen, J.; Jensen, H. J. Aa..; Jørgensen, P. Solution of the large matrix equations which occur in response theory. J. Comp. Phys. 1988, 74, 265–282. (51) List, N. H.; Coriani, S.; Christiansen, O.; Kongsted, J. Identifying the Hamiltonian structure in linear response theory. J. Chem. Phys. 2014, 140, 224103. (52) List, N. H.; Coriani, S.; Kongsted, J.; Christiansen, O. Lanczos-driven coupled–cluster damped linear response theory for molecules in polarizable environments. J. Chem. Phys. 2014, 141, 244107. (53) Hunt, H. D.; Simpson, W. T. Spectra of Simple Amides in the Vacuum Ultraviolet. J. Am. Chem. Soc. 1953, 75, 4540–4543. (54) Kaya, K.; Nagakura, S. Vacuum ultraviolet absorption spectra of simple amides. Theor. Chem. Acc. 1967, 7, 117–123.

26


Page 26 of 31

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


(55) Basch, H.; Robin, M. B.; Kuebler, N. A. Electronic States of the Amide Group. J. Chem. Phys. 1967, 47, 1201–1210. (56) Serrano-Andrés, L.; F¨ ulscher, M. P. Theoretical Study of the Electronic Spectroscopy of Peptides. 1. The Peptidic Bond: Primary, Secondary, and Tertiary Amides. J. Am. Chem. Soc. 1996, 118, 12190–12199. (57) Besley, N. A.; Hirst, J. D. Ab Initio Study of the Effect of Solvation on the Electronic Spectra of Formamide and N-Methylacetamide. J. Phys. Chem. A 1998, 102, 10791– 10797. (58) Besley, N. A.; Hirst, J. D. Ab Initio Study of the Electronic Spectrum of Formamide with Explicit Solvent. J. Am. Chem. Soc. 1999, 121, 8559–8566. (59) Sneskov, K.; Matito, E.; Kongsted, J.; Christiansen, O. Approximate Inclusion of Triple Excitations in Combined Coupled Cluster/Molecular Mechanics: Calculations of Electronic Excitation Energies in Solution for Acrolein, Water, Formamide, and N-Methylacetamide. J. Chem. Theory Comput. 2010, 6, 839–850. (60) Hrˇsak, D.; Marefat Khah, A.; Christiansen, O.; Hättig, C. Polarizable Embedded RI-CC2 Method for Two-Photon Absorption Calculations. J. Chem. Theory Comput. 2015, 11, 3669–3678. (61) Hrˇsak, D.; Olsen, J. M. H.; Kongsted, J. Embedding potentials of aqueous formamide. 2018; https://doi.org/10.5281/zenodo.1149019. (62) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (63) 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.

27



(64) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can. J. Phys. 1980, 58, 1200–1211. (65) Stephens, P. J.; 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. (66) 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. (67) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358–1371. (68) Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3d elements ScZn. J. Chem. Phys. 2005, 123, 064107. (69) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (70) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P.-A.; Neogrády, P.; Pedersen, T. B.; Pitoˇ nák, M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. MOLCAS 7: The Next Generation. J. Comput. Chem. 2010, 31, 224–247. (71) Gagliardi, L.; Lindh, R.; Karlström, G. Local properties of quantum chemical systems: The LoProp approach. J. Chem. Phys. 2004, 121, 4494–4500. (72) Hrˇsak, D.; Olsen, J. M. H.; Kongsted, J. Optimization and transferability of non-

28


Page 28 of 31

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


electrostatic repulsion in the polarizable density embedding model. J. Comput. Chem. 2017, 38, 2108–2117. (73) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannen¨ Foresman, J. B.; Ortiz, J. V.; berg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision D.01. Gaussian Inc., Wallingford CT, 2013. (74) Aidas, K.; Møgelhøj, A.; Nilsson, E. J. K.; Johnson, M. S.; Mikkelsen, K. V.; Christiansen, O.; Söderhjelm, P.; Kongsted, J. On the performance of quantum chemical methods to predict solvatochromic effects: The case of acrolein in aqueous solution. J. Chem. Phys. 2008, 128, 194503. (75) Hrˇsak, D.; Olsen, J. M. H.; Kongsted, J. Embedding potentials of aqueous acrolein. 2018; https://doi.org/10.5281/zenodo.1149061. (76) Kjellgren, E. MM-MD and QM/MM-MD trajectories of para-nitrophenolate derivatives in aqueous solution. 2017; https://doi.org/10.5281/zenodo.1043593.

29



(77) Hrˇsak, D.; Olsen, J. M. H.; Kongsted, J. Embedding potentials of aqueous 4nitrophenolate. 2018; https://doi.org/10.5281/zenodo.1149072. (78) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li–F. J. Comput. Chem. 1983, 4, 294–301. (79) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. J. Chem. Phys. 1982, 77, 3654–3665. (80) Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chem. Acc. 1973, 28, 213–222. (81) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self—Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian—Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257–2261. (82) Bowers, G. N.; McComb, R. B.; Christensen, R. G.; Schaffer, R. High-purity 4nitrophenol: purification, characterization, and specifications for use as a spectrophotometric reference material. Clin. Chem. 1980, 26, 724–729. (83) Reinholdt, P.; Kongsted, J.; Olsen, J. M. H. Polarizable Density Embedding: A Solution to the Electron Spill-Out Problem in Multiscale Modeling. The Journal of Physical Chemistry Letters 2017, 8, 5949–5958.

30


Page 30 of 31

Graphical TOC Entry UV/vis spectrum of formamide(aq) π→π∗

Polarizable density embedding Vemb=Vmul+Vfd+Vrep+Vind

oscillator strength

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


Rydberg states

n→π∗ π→π∗

Polarizable embedding Vemb=Vmul+Vind

oscillator strength

Page 31 of 31

Rydberg states

n→π∗ energy

31


Polarizable Density Embedding Coupled Cluster Method - Journal of

Recommend Documents