Relativistic Polarizable Embedding - Journal of Chemical Theory and

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Relativistic Polarizable Embedding Erik Donovan Hedegård, Radovan Bast, Jacob Kongsted, Jógvan Magnus Haugaard Olsen, and Hans Jørgen Aagaard Jensen J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00162 • Publication Date (Web): 11 May 2017 Downloaded from http://pubs.acs.org on May 25, 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 35 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

Relativistic Polarizable Embedding Erik Donovan Hedeg˚ ard,∗,† Radovan Bast,‡ Jacob Kongsted ,¶ Jógvan Magnus Haugaard Olsen,¶ and Hans Jørgen Aagaard Jensen∗,¶ †Department of Theoretical Chemistry, Lund University, Lund, Sweden ‡High Performance Computing Group, UiT The Arctic University of Norway, Tromsø, Norway ¶Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Odense M, Denmark E-mail: [email protected]; [email protected]

Abstract Most chemistry, including chemistry where relativistic effects are important, occurs in an environment, and in many cases this environment has a significant effect on the chemistry. In non-relativistic quantum chemistry a lot of progress has been achieved with respect to including environments such as a solvent or a protein in the calculations, and now is the time to extend the possibilities for doing that also in relativistic quantum chemistry. The polarizable embedding (PE) model efficiently incorporates electrostatic effects of the environment by describing it as a collection of localized electric multipoles and polarizabilities obtained through quantum chemical calculations. In this paper, we present the theory and implementation of four- or exact two-components Hamiltonians within a PE framework. We denote the methods the PE-4c-DFT and PE-X2C-DFT models. The models include a linear response formalism to calculate time-dependent (TD) properties; PE-TD-4c-DFT and PE-TD-X2C-DFT. With this first implementation, we calculate the PE-TD-4c-PBE0 excitation energies of the TcO4– and ReO4– ions

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

in an explicit water solvent. This initial investigation focuses on the relative size of relativistic and solvent contributions to the excitation energies. The solvent effect is divided into an indirect solvent effect due to the structural perturbation of the XO4– ion, and a direct, electrostatic effect. The relativistic effects as well as both types of solvent effects are found to contribute to a shift of the excitation energies, but to different extents depending on the ion and the electronic transition in question.

Introduction The theory of special relativity affects chemistry with what usually is called “relativistic effects”. 1–3 Most of the chemistry involving only lighter elements from the periodic table can be described adequately without considering relativistic effects, but not all. For instance, relativistic effects are essential for Electron Spin Resonance spectroscopy and phosphorescence. For molecules with heavy elements, inclusion of at least some relativistic effects is always essential for useful computational results. The chemistry of the heavier elements remains relevant; they are increasingly exploited industrially, e.g. as catalysts or in batteries, 4,5 and a better understanding of actinides in solution is vital for safer handling of nuclear waste. 6 Consequently, many non-relativistic electronic structure methods for atoms and molecules have been extended to the relativistic regime. 7,8 However, only a small fraction of the chemistry where relativistic effects are important occurs in the gas phase, and we therefore see an obvious need for further development of efficient, relativistic computational methods that can include an environment (e.g., a solvent or a protein) to simulate the actual experimental conditions. From a number of non-relativistic quantum mechanical calculations it is known that ignoring the presence of such environments can lead to erroneous results. 9–14 To include an environment in relativistic calculations, we describe here the theory and first implementation of the polarizable embedding (PE) model 15,16 within four-component (“4c”) and eXact two-component (“X2C’) relativistic frameworks. We present the theory and implementation of both the optimization of a 4c/X2C Hartree-Fock (HF) or Kohn-Sham 2


Page 2 of 35

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


(KS) determinant with PE and its subsequent use in linear response calculations, here concentrating on electronic excitation energies and intensities. The PE model is different from implicit solvent models such as the polarizable continuum model (PCM), 17 whose implementation in combination with time dependent (TD) 4c-DFT and X2C-DFT was described recently. 18 Implicit solvent models lack explicit interactions, meaning that they are less accurate for solvents with extensive hydrogen bonding and/or heterogeneous charge distributions. Further, they cannot handle highly inhomogeneous environments such as chromophores inside proteins (see e.g. Ref. 11). Methods that include explicit interactions are for instance the hybrid models based on quantum mechanics and molecular mechanics (QM/MM). 19,20 In many QM/MM models, the electrostatic interactions with the environment are described through point charges and explicit polarization of the environment is ignored. However, some chemical phenomena require higher accuracy: This is for instance the case for the description of electronic excitations and excited states. The accuracy of explicit environment models can be enhanced by replacing the point charges by higher-order multipoles. For even greater accuracy, the mutual polarization of the quantum mechanical system and the environment must be accounted for. We have over the last years developed the PE model to achieve this goal by representing the environment in terms of atom-centered multipoles and polarizabilities. Both multipoles and polarizabilities are obtained by dividing the environment into chemically meaningful fragments and performing quantum-chemical calculations on each of these fragments in vacuum. The PE model has today been combined with several quantum-chemical methods. 15,16,21–29 From the literature, we know of only two other explicit models that have addressed some of the issues of environment effects in relativistic quantum chemistry: frozen density embedding, either as 4c DFT-in-DFT or as 4c coupled-cluster (CC) and DFT (CC-in-DFT). 30,31 In addition to the work presented here, the PE model has formed the basis for a very recent X2C-PERI-CC2 model. 22,32,33 The paper is organized as follows. In the next two sections, we describe the theory,

3



and comment on a few aspects of the implementation, which is based on the polarizable embedding library (PElib), 34 and the general four-component self-consistent-field (4c-SCF) and response modules developed within the framework of the DIRAC program. 35–37 The implementation has been done in a development version of the DIRAC16 release 38 and is expected to be part of the DIRAC17 release. The theory presented here focuses on 4c Kohn-Sham density functional theory (KS-DFT), but our implementation also covers the 4c-HF scheme and X2C exact two-component Hamiltonians. After the first two sections, we illustrate one of the possibilities with the new code by investigation of both relativistic effects and solvent effects on the electronic spectra of TcO4– and ReO4– in a snapshot of water molecules, and we finish with some conclusions and perspectives.

Theory In the following, the Born-Oppenheimer approximation is always assumed and all Hamiltonians cover exclusively electronic degrees of freedom within a scalar potential from the nuclei. In the first two subsections we summarize and introduce notation for the building blocks: the existing Dirac KS and polarizable embedding methodologies. We then derive the extensions needed to combine the two models. In the third subsection we combine relativistic Dirac-KS response theory with polarizable embedding, and in the final subsection we summarize external effective field effects as described by the polarizable theory.

The Vacuum Dirac Kohn-Sham Equations Both HF theory 35,39–41 and the KS formulation of DFT 42 have been extended to relativistic frameworks. Here we summarize the latter using the formulation by Saue and Helgaker, 36 except we will in some places use zero as a subscript to designate the vacuum (no environment) case. Following Saue and Helgaker, we employ an exponential parameterization of the

4


Page 4 of 35

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


KS reference state, |0(κ)i = exp(−ˆ κ)|0i,

(1)

where κ ˆ is defined in terms of creation (ˆ a†p ) and annihilation (ˆ aq ) operators

κ ˆ=

X

κpq a ˆ†p a ˆq ,

(2)

pq

and κpq are elements of the anti-Hermitian complex spinor rotation matrix. The |0i wavefunction is represented as a single determinant, comprised of spinors, {φp }. We use p, q, r, s for general spinors; i, j, k, l for occupied spinors; and a, b, c, d for unoccupied spinors. The energy functional can be written

E0 [ρ(D(κ))] =

X

hpq Dpq (κ) +

pq

X

jpq [D(κ)]Dpq (κ) + Exc [ρ(D(κ))] + Enn ,

(3)

pq

where D(κ) is the one-electron reduced density matrix (1-RDM). The terms in Eq. (3) are to a large degree formally similar to non-relativistic theory; a discussion of the terms and their differences to their non-relativistic counterparts is given with a KS-DFT framework in Ref. 36. We refer to this paper for the original vacuum implementation. Here, we only note that the main difference lies in the first term of Eq. (3) which contains integrals over the one-electron Hamiltonian  ˆ=h ˆ D + Vêxt =  h 







ˆ )   Vext I2 c(σ · p 02  + . 2 ˆ ) −2c I2 c(σ · p 02 Vext I2 02

(4)

We again refer to Ref. 36 for the definition of constants and operators contained in Eq. (4). The wave function is optimized under the requirement that the electronic gradient with respect to the wave function parameters vanishes

[1] Epq =

∂E0 [ρ(D(κ))] = 0. ∂κ∗pq 5


(5)


Page 6 of 35

[1]

The elements of the gradient, Epq , can be shown to be equal to the negative off-diagonal elements in the occupied-unoccupied block of the 4c-KS matrix 36 Eia = −h0|[fˆ0KS , a ˆ†i a â ]|0i = −FiaKS . [1]

(6)

The zero-field KS-operator (no environment), corresponding to the KS-matrix in Eq. (6), can generally be written as

fˆ0KS =

X

fpq a ˆ†p a ˆq =

pq

X

hD,pq + jpq [D(κ)] + vxc,pq [ρ(D(κ))] a ˆ†p a ˆq .

(7)

pq

ˆ D must deal Note that variational theories that employ an SCF scheme including the 4c h with unbound solutions. In the DIRAC program, the 4c optimization is carried out using the min-max principle. 43,44

Polarizable Embedding We now summarize the calculation of the PE energy as defined in the PE model. 15,16 The PE energy is associated with the interaction between an embedded system and its environment, resulting in both permanent electrostatic and induction (polarization) interactions, E es and E ind , respectively, E PE = E es + E ind .

(8)

The electrostatic energy is described through the interactions between the nuclei and electrons in the quantum region and the permanent charges, qs0 , dipole moments, µs0 , quadrupole moments, Qs0 , etc., defined for each site, s0 , in the environment, es . E es = h0(κ)|ˆ v es |0(κ)i + Enes + Em

(9)

The operator vês in the first term of Eq. (9) describes the electron-multipole interactions,

6


Page 7 of 35 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 is defined as, vês =

X

es † ˆq , a ˆp a vpq

(10)

pq

so that the expectation value becomes

h0(κ)|ˆ v es |0(κ)i =

X

es vpq Dpq (κ).

(11)

pq

The second term in Eq. (9), Enes , is the nuclear-multipole interaction energy, and the third es and final term, Em , is the multipole-multipole interaction energy.

The induction energy, which is due to polarization of the environment, is defined as 1 1 E ind = − E tot [D(κ)]T T E tot [D(κ)] = − µind [D(κ)]T E tot [D(κ)], 2 2

(12)

where T is the inverse of a (super) matrix with inverse dipole-dipole polarizabilities, α−1 s0 , on (2)

the block-diagonal and (negative) dipole-dipole interaction tensors, −Ts0 s00 , as off-diagonal blocks. The E tot [D(κ)] and µind [D(κ)] terms are (super) vectors that contain the electric fields and induced dipoles from the electrons, nuclei, and multipoles at all polarizable sites. A Cartesian component of the field at a (polarizable) site can thus be written as a sum

EXtot [D(κ)](rs0 ) = h0(κ)|EˆXe (rs0 )|0(κ)i + EXn (rs0 ) + EXm (rs0 ),

(13)

where the electronic electric field operator is defined as

EˆXe (rs0 ) =

X

e Epq,X (rs0 )ˆ a†p a ˆq .

(14)

pq

The expectation value of Eê can be written in terms of the 1-RDM as

h0(κ)|EˆXe (rs0 )|0(κ)i =

X

e Epq,X (rs0 )Dpq (κ),

pq

7


(15)


Page 8 of 35

e where Epq,X (rs0 ) is a Cartesian component of an electric field integral,

e Epq,X (rs0 ) = −hφp |

rX − rX,s0 |φq i. |r − rs0 |3

(16)

We can now combine Eqs. (9), (12), and (3) to yield the total energy-functional 1 es . E = E0 [ρ(D(κ))] + h0(κ)|ˆ v es |0(κ)i − µind [D(κ)]T E tot [D(κ)] + Enes + Em 2

(17)

Before we continue with the gradient contributions from the environment, we note that the induced dipole moment can be split into the individual contributions according to the fields given in Eq. (13)

µind [D(κ)] = µeind [D(κ)] + µnind + µm ind ,

(18)

and we will use the notation described above for the different contribution in the next subsection. In this subsection we continue to derive the electronic gradient with respect to the wave function parameters [1] Eai

∂E ∂E0 [ρ(D(κ))] = = ∂κ∗ia κ=0 ∂κ∗ia κ=0 ∂h0(κ)|ˆ v es |0(κ)i 1 ∂µind [D(κ)]T E tot [D(κ)] + − . ∂κ∗ia 2 ∂κ∗ia κ=0 κ=0

(19)

Inserting the expectation values from Eq. (11) and (15) in Eq. (19) and employing a BakerCampbell-Hausdorff (BCH) expansion of the 1-RDM gives [1] Eai = −h0|[fˆ0KS , a ˆ†a a î ]|0i − h0|[ˆ v PE , a ˆ†a a î ]|0i,

8


(20)

Page 9 of 35 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 e

e

ˆ = vês − µind [D(κ)] E ˆ . vˆPE = vês − E tot [D(κ)]T T E

(21)

Evaluating the integrals gives the modified gradient expression [1]

Eia = −FiaKS − FiaPE .

(22)

Relativistic Kohn-Sham response theory with polarizable embedding For the calculation of molecular properties we shall employ time-dependent perturbation theory, following Refs. 37,45–47. The theory of 4c-DFT response theory and its implementation in the DIRAC program were originally described by Bast et al. 48 Here we combine this theory with the PE model. The interaction operator, Vˆ (t), is a periodic, time-dependent perturbation defined as Vˆ (t) =

X

exp(−iωk t)

X

ˆX , EXωk H

(23)

X

k

with the Fourier transform components defined in terms of the field strength components, ˆ X . Note that in EXωk , of the external perturbation described by the property operator, H the following, fields with origin in the operator in Eq. (23) will be labeled with frequencies, for example ωk , ωa , or ωb . For the components of the field from the quantum region and the static multipoles at the polarizable sites in the environment, we use the same labels as in previous section. Thus, EXe (rs0 ), EXn (rs0 ), and EXm (rs0 ) denote the X component of the field at a polarizable site, rs0 , from the quantum-region electrons, quantum-region nuclei, and multipoles in the environment, respectively, while EXωk denotes component X of the external field. The contribution from EXωk to the field on the polarizable sites, EXωk (rs0 ), is not considered at first, but will be taken into account in the next subsection. Molecular properties can be obtained from response functions that in turn can be defined from an 9



Page 10 of 35

expansion of the time-dependent expectation value

ˆ X |˜0i = h0|H ˆ X |0i + h˜0|H

X

exp(−iωk t)

X

ˆ X ; Vˆ ωk iiω hhH k Y

Y

k

X 1X ˆ X ; Vˆ ωk ; Vˆ ωk0 iiω ,ω 0 + · · · . exp(−i(ωk + ωk0 )t) hhH + k k Y Z 2 k,k0 Y,Z

(24)

ˆ X ; Vˆ ωk iiω , and hhH ˆ X ; Vˆ ωk ; Vˆ ωk0 iiω ,ω 0 The linear and quadratic response functions are hhH k k k Y Y Z and |˜0i is the time-dependent wave function. Within a DFT framework, |˜0i is obtained by extending the parameterization of the KS reference determinant in Eq. (1) to the timedependent form

|˜0i = exp(−ˆ κ(t))|0i,

(25)

where κ ˆ (t) can be expanded in the perturbation series. Since we here are only concerned with linear response we terminate the perturbation series at first order and write

κ ˆ (t) =

X X k

† −ωk ∗ + κai qâi exp(−iωk t). κωaik qâi

(26)

ai

The quantities qîa = a ˆ†i a â include both rotations between occupied and unoccupied electronic spinors, and occupied electronic spinors and positronic spinors (the latter are always unoccupied). We now assign the periodic perturbation to have a period of T , and define the R T /2 time-averaged quasi energy, {Q}T = −T /2 dt Q(t), within a KS framework 47,49 (see also Ref. 37 for a derivation within an HF framework)

¯ T+ {Q}T = {Q}

XX k

10

EXωk {QX }T ,

X


(27)

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


where we have introduced

¯ T = Q0 + {QPE }T . {Q}

(28)

The time-averaged quasi-energy concerned with the time-derivative and the KS operator, fˆ0KS , is contained in Q0 , and explicit expressions for this term as well as for {QX }T can be found in Refs. 37,47,49. Here we are mainly concerned with {QPE }T which is defined as

{Q

PE

1 es T ˜ ˜ }T = h˜0|ˆ v |˜0i − µind [D(κ)] E[D(κ)] 2 T X 1 m T m es ˜ n T n e T e ˜ ˜ = vpq Dpq (κ) − µ E + µind E + µind [D(κ)] E pq Dpq (κ) 2 ind T pq o Xn T ˜ pq (κ) + E n + µn T E e D ˜ µm E epq D , (29) − ind ind pq pq (κ) T

pq

˜ ˜ pq (κ) is the time-dependent 1-RDM and µe [D(κ)] where D is the time-dependent analogue ind of µeind [D(κ)] in Eq. (18). It can be shown that the linear response function can be obtained as 50,51 ¯ T d2 {Q} ˆ X ; Vˆ ωk iiω = hhH Y dEXω dEYωk

ω + ωk = 0,

(30)

which by carrying out the differentiation becomes 2 ¯ ¯ T d2 {Q} ∂ {Q}T ∂κωpqk ∂κωrsk ∂{QA }T ∂κωpqk ∂{QB }T ∂κωrsk = + + δ(ωa + ωb ). dEAωa dEBωb E=0 ∂κωpqk ∂κωrsk ∂EBωb ∂EAωa ∂κpq ∂EBωb ∂κωrsk ∂EAωa E=0 (31) The amplitudes,

∂κωk , ∂EYω

can be obtained from

X ∂ 2 {Q} ¯ T ¯ T ∂κωk ∂{QX }T d ∂{Q} rs = + = 0. ωk ωk ωk ω dEXω ∂κωpqk E=0 ∂κ ∂κ ∂E ∂κ pq rs pq X rs

11


(32)


Page 12 of 35

Here we note that the difference in Eqs. (30)–(32) from corresponding vacuum equations is ¯ T in Eq. (28) rather than Q0 . Eqs. (31) and (32) are evaluated by expanding the use of {Q} the time-dependent 1-RDM through a BCH expansion, and inserting it into the expression for Q0 , {QX }T and {QPE }T . The resulting expression for the amplitudes is in matrix notation given as −1 [1] ∂κωk EX δ(ω + ωk ), = − E[2] − ωk S[2] ω ∂EX

(33)

and the linear response function becomes ˆ X ; Vˆ ωk iiω = −E[1] E[2] − ωk S[2] hhH k Y X

−1

[1]

EY .

(34)

The elements of the property-gradient are [1] ˆ X ]|0i, Eai,X = h0|[ˆ qai , H

(35)

and the E[2] matrix is given as 



 A B  E[2] =  , B ∗ A∗

(36)

† † Aai;bj = h0|[ˆ qai , [ˆ qbj , fˆ0 + vˆPE ]]|0i + h0|[ˆ qai , vînd ]]|0i

(37)

Bai;bj = h0|[ˆ qai , [ˆ qbj , fˆ0 + vˆPE ]]|0i + h0|[ˆ qai , vînd ]|0i.

(38)

with the matrix elements

The part that contains fˆ0 + vˆPE is, just as the metric, S[2] , equivalent to the vacuum expressions 37,45–47,49 and we refer to the literature for these expressions. However, the E[2] part also has an additional contribution (the last term of Eqs. 37–38) that is an unique (linear

12


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


response) term of polarizable embedding models. 15,16,52,53 The vînd operator of this term is given by

vînd =

X

e

† ˆ ]|0iT TE e a ˆj , h0|[ˆ qbj , E bj ˆb a

(39)

bj

and accounts for the change of the environment interaction due to the alternation of the molecular electron density under influence of the perturbing (electromagnetic) field.

External Effective Field Effects In the previous subsection, the contributions from the environment were exclusively included in the electronic Hessian, E[2] , in Eq. (36). However, the introduction of an external field has an additional consequence that we have not yet accounted for. The external field will, even in a fictitious system without the quantum mechanical part, induce a polarization (and thereby a field) within the environment. The effect is well known from continuum models 54–57 and has also been described within the discrete reaction field model. 58 We have previously denoted this effect the effective external field (EEF) effect and it was recently described within a non-relativistic PE scheme. 59 To include the effect we extend the total field in Eq. (13) with ˜ ˜ the perturbing field, E(t), to obtain E tot (t) = hEê i + E n + E m + E(t). The total field can thus be obtained as an expectation value over the time-dependent wave function ˜ tot (t) = h˜0|E ˆ tot (t)|˜0i = E ˜ e (t) + E n + E m + E(t), ˜ E

(40)

˜ e is the time-dependent form of Eq. (15). By the modification of Eq. (13), we also where E obtain a modification of the time-averaged quasi energy in Eq. (27), which now becomes

¯ T+ {Qeef }T = {Q}

XX k

13

¯ X }T . EXωk {Q

X


(41)


Page 14 of 35

In analogy with Eq. (27), we have in Eq. (41) used the notation

¯ X }T = {QX }T + {QPE {Q X }T ,

(42)

with ( {QPE X }T =

−

X s00 s0

Ts00 s0

1 X ωk E n (rs0 ) + E m (rs0 ) + h˜0|Eê (rs0 )|˜0i + EY (rs0 ) exp(−iωk t) 2 Y

) . T

(43) Replacing Eq. (27) with Eq. (41) naturally also leads to modification of Eqs. (31) and (32) which for the former case changes to 2 ¯ ∂ {Q}T ∂κωpqk ∂κωrsk d2 {Qeef }T = δ(ωa + ωb ) dEAωa dEBωb E=0 ∂κωpqk ∂κωrsk ∂EBωb ∂EAωa E=0 ¯ ¯ B }T ∂κωrsk ∂{QA }T ∂κωpqk ∂{Q + + δ(ωa + ωb ), ∂κωpqk ∂EBωb ∂κωrsk ∂EBωb E=0

(44)

and a similar modification, i.e., extension with {QPE X }T , is obtained for the amplitudes in ¯ X }T terms out Eq. (32). By writing the {Q 2 ¯ d2 {Qeef }T ∂ {Q}T ∂κωpqk ∂κωrsk ∂QA ∂κωpqk ∂QB ∂κωrsk = + + δ(ωa + ωb ) dEAωa dEBωb E=0 ∂κωpqk ∂κωrsk ∂EBωb ∂EAωa ∂κωpqk ∂EBωb ∂κωrsk ∂EAωa E=0 PE ωk ωk ∂QA ∂κpq ∂QPE B ∂κrs δ(ωa + ωb ), (45) + + ∂κωpqk ∂EBωb ∂κωrsk ∂EBωb E=0 and comparing with Eq. (31), we can identify the last term of Eq. (45) as the EEF correction [1],PE

term, EX

, that can be added to the property-gradient in Eq. (35). Thus ¯ [1] = E[1] + E[1],PE , E X X X

(46)

replaces the property-gradient terms in Eqs. (33) and (34). The elements of the EEF correction gradient-term can be written as the derivative over an induced dipole from the external 14


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


field [1],PE

Eai,X =

dµext,ω ind,X h0|[ˆ qai , Eê ]|0i, dEXω

where we have used the notation µext,ω ind,X =

P

s00 s0

(47)

Ts00 s0 EXω (rs0 ). We finally note that Eq. (47)

does not include terms of higher order in the external field, although these enters in Eq. (43). Leaving these terms out is consistent with a linear response formalism.

Implementation The elements of the gradient that belong to the vacuum part in Eq. (22) are just the (negative) off-diagonal elements of the regular KS matrix. 36 The optimization condition [1]

Eia = 0

∀ ai,

(48)

therefore corresponds to a conventional diagonalization of the KS matrix. By means of Eq. (22), the environment is thus included in the SCF step. In a 4c formulation the individual elements

FiaKS = hi|fˆ0 |ai,

(49)

can be formulated as two-component spinors, i.e.,

|ii =

NL X κ





L  |χκ i  L   Cκ + 0

Ns X λ



 0

 

|χSλ i

 S  Cλ

15

 ;




X  |χiα i  |χX i =  . i X |χiβ i

(50)


Page 16 of 35

The primitive basis functions |χLκ i and |χSλ i are themselves two component spinors, comprised of individual (spin) components of the Kramers pairs   |χLm i = 

 |χLα i |χLβ i

  |χSm i = 

 

 |χSα i |χSβ i

 ,

(51)

where X = L, S. As such, the KS matrix in Eq. (49) can be composed in the following 2 × 2 block 

 LL

LS

F   F FKS , 0 =  FSL FSS

(52)

XY ˆKS Y Fκλ = hχX κ |f0 |χλ i.

(53)

with the generic matrix blocks

As also noted in the context of PCM, 18 the static electrostatic and polarization contributions within the PE operator are diagonal 4c operators and thus do not couple large and small component wave functions, i.e., from Eq. (16) we obtain

E eκγ (rs0 )

=

−hχX κ |

r − rs0 |r − rs0 |3

|χYλ i

Z = δXY

r − rs0 Ωκλ (r)dr, |r − rs0 |3

(54)

es which is only non-zero for X = Y (the same argument can be put forth for vκλ ). Note

that the last integral does not reduce to δκλ as the primitive orbital basis is not necessarily orthogonal for molecular calculations. Accordingly, the structure of the implemented KS matrix is 



LL PE,LL FLS  F +F  FKS =  . SL SS PE,SS F F +F

16


(55)

Page 17 of 35 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 note should be added regarding the use of X2C Hamiltonians. In this case, the one-electron part of the Fock matrix is transformed (see e.g. Ref. 7 for details). Strictly, all external fields (including a field from a environment) should also be be transformed. Currently, our PEX2C-DFT implementation only involves transformation of the vacuum terms. Thus, the PE operator is left untransformed, and the untransformed operator is carried over to the linear response formalism described below (following the general implementation in the DIRAC program, all property-operators are always transformed). For the linear response part, the response eigenvalue equations yielding predicted excitation energies

E[2] Xk = ωk S[2] Xk ,

(56)

are solved for the requested number of roots using a Davidson diagonalization scheme with trial vectors, 48 {bi }. What now remains is to specify the changes in the algorithm for the σ-vectors

σ = E[2] bi ,

(57)

arising from PE compared to the previously published vacuum expressions. The final computational expression becomes

σai =

X

ˆ e ]|0i, ˜ eind T h0|[−ˆ (Aai;bj bbj + Bai;bj bjb ) = h0|[−ˆ qai , f˜0KS + v˜PE ]|0i + µ qai , E

(58)

bj

where it was used that † h0|[qai ,[

X

† bbj qˆbj , fˆ0KS + vˆPE ]]|0i = h0|[qai , f˜0KS + v˜PE ]|0i,

bj

17


(59)


Page 18 of 35

and

˜ eind = Th0|[ µ

X

e

ˆ ]|0i. bbj qˆbj , E

(60)

bj

˜ are so-called one-index transformed operators (see e.g. Ref. 37 The operators of the form X for the vacuum expressions and Refs. 15 or 24 for the PE variants). We have in Eq. (58) followed Ref. 37 and used that the inactive-inactive and secondary-secondary blocks are zero, which means that we can we can use a free summation in the one-index transformed operators. Evaluating Eq. (58) yields

KS PE σia = f˜ia + v˜ia +µ ˜ eind T E eai .

(61)

Computational Details The implementation described in the previous Section was done in the DIRAC program 38 using PElib 34 for the PE contributions, and has been used to calculate the excitation energies of TcO4– and ReO4– in an explicit water solvent.

Optimized Structures For both TcO4– and ReO4– , the structures were optimized with the Turbomole program, 60 using a def2-SV(P) basis set 61,62 and the BP86 functional. 63,64 For the solvated structures, the XO4– (where X is Tc or Re) were immersed in an environment comprised of 516 water molecules. The initial structures for these calculations were constructed from a single snapshot taken from an ab-initio molecular dynamics simulation of MnO4– in a water droplet (see Ref. 65 for details). The Mn atom was replaced by a Tc or Re atom and all water molecules more than 4 ˚ A away were removed. The structures of the systems with XO4– and water molecules within 4 ˚ A were optimized with the water molecules frozen to their initial posi-

18


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


tions. The optimized structures of XO4– were then inserted back into the full environment of 516 water molecules. Using the same functional and basis set we also optimized symmetric (Td ) structures in vacuum of both ions.

Excitation Energies From the optimized structures, we carried out excitation energy calculations with XO4– as the quantum region, employing the PBE0 66 functional. For the large component, we employed the aug-cc-pVDZ 67,68 basis set for oxygen and the Dyall double-zeta basis set 69–71 for the metals (all uncontracted). The small component basis sets were generated by employing the kinetic balance relation. The calculations of excitation energies were carried out in a step-wise manner. Thus, we initially carried out calculations on symmetric (Td ) structures in vacuum with a non-relativistic Hamiltonian (using the “.NONREL” keyword in DIRAC). Next, we carried out KS-4c Dirac-Coulomb calculations for the symmetric structures. We then performed the calculation on the solvated structure, but without including the solvent in the calculation. Comparing these excitation energies to the excitation energies calculated for a symmetric structure will allow us to analyze the indirect effect of the solvent environment, i.e., the perturbation of the structure due to the interactions with explicit solvent molecules. Finally, we included the solvent potential, employing the PE model, which gives the direct (electrostatic) solvent effect on the electron density of the quantum region. In these PE calculations, the solvent water molecules are represented by a potential comprised of multipoles up to (and including) 2nd order and anisotropic dipole-dipole polarizabilities. We denote this potential “M2P2” in accordance with previous papers. 72 Both multipoles and polarizabilities are calculated with DFT (B3LYP 64,73,74 ) using the LoProp 75 procedure as implemented in MOLCAS. 76 The LoProp procedure requires specially constructed basis sets 75 and, accordingly, we employed here a recontracted aug-cc-pVTZ 67,68 basis set. At this point we note that the XO4– ions are usually assumed to have Td symmetry, and the lowest-lying transition is thus triply-degenerate. In accordance with experimental stud19



ies, 77 we use here the non-relativistic point-group labels for these transitions (either of T1 or T2 ). However, a structure obtained from a MD trajectory has lower symmetry leading to a (sometimes significant) splitting of the T1 and T2 bands, which complicates the spectra. Additionally, the symmetry-lowering also leads to that excitations of T1 symmetry gain intensity, adding a further complication, when comparing to a vacuum calculation (excitations between the ground-state and he T1 states are symmetry-forbidden in Td symmetry). The analysis of the results is here carried out on the basis of a spectrum made by a convolution of Gaussian-broadened peaks that are based on the calculated vertical excitation energies and using a broadening parameter (related to the full width at half maximum) of 0.15 eV. The broadening is not selected to mimic experimental conditions (where 0.3–0.4 eV is normally employed). Rather it is chosen to obtain a sufficient resolution to discuss the composition of the individual bands, but without discussing each individual transition.

− TcO− 4 and ReO4 ions in water In this section we investigate the effect of an aqueous environment on the electronic spectra of pertechnate, TcO4– and perrhenate, ReO4– . The lighter permanganate (MnO4– ) has been the target of a plethora of studies (see e.g. Ref. 78 for an overview), but its heavier congeners, TcO4– and ReO4– , have been much less studied. TcO4– was recently studied with timedependent DFT (TD-DFT), equation-of-motion coupled cluster (EOM-CC) and restricted active space second-order perturbation theory (RASPT2). 79 Earlier studies also exist with symmetry adapted cluster configuration interaction (SAC-CI) 80 and Xα 81 methods. For ReO4– , we are only aware of one study, which employed TD-X2C-DFT. 82 Neither TcO4– , nor ReO4– have been studied within an explicit solvent environment. On the experimental side, a number of studies on the UV-vis spectra of TcO4– and ReO4– have been reported. 77,83,84 Before discussing the two molecules and their spectra, we should emphasize that we here focus on a detailed analysis, based on one specific snapshot. While we can thus analyze

20


Page 20 of 35

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


indirect and direct solvent effect in great detail, we lack the effect of dynamics, and the analysis should be viewed in this light. Indeed, general conclusions as well as comparison with experimental values can only be done after averaging over a statistically relevant sample of snapshots and – in case of the XO4– ions – possibly also after inclusion of vibronic coupling such as done in Ref. 85 for MnO4– . We start by discussing TcO4– . The calculated absorption spectra are shown in Figure 1 (a)–(d), and we denote the peaks with numbers 1–4. We also report the band-maxima in Table S1 in the Supporting Information (SI). Only peaks 2 and 4 are intense for the Td symmetric systems in Figure 1 (a) and (b). By comparing these two figures, the change upon introducing a relativistic Hamiltonian can be inferred. This effect is mainly a shift of the two bands in the spectrum. Analysis of the transitions under these two bands shows that they are the symmetry-allowed 1T2 and 2T2 states. The shifts due to relativity can be read from columns 1 and 2 of Table S1 to be around 0.2 eV. We note that the first band (1T2 ) also includes a low-intensity triplet state with an excitation energy of 35900 cm−1 (4.45 eV). Next, we investigate the spectrum in Figure 1 (c), which is calculated as a vacuum calculation, but employing the solvated structure. The reduction of symmetry leads to the introduction of a new band, 1, at 33000 cm−1 (4.09 eV). The new band has its origin in the 1T1 states, which gain intensity as they are no longer strictly symmetry-forbidden. The same lowering of symmetry leads to splitting of the 1T2 states and the spectrum now has two maxima (band 2) at 36300 and 38800 cm−1 (i.e. 4.50 and 4.81 eV) with origin in these states. This splitting also occurs for the 2T2 state (peak 4), but here it is somewhat smaller. The electrostatic effect from the multipoles and polarizabilities in the environment can be seen by comparison of Figure 1 (c) and (d). The effect is rather small for the bands 1 and 2 with origin in 1T1 and 1T2 (for the numerical values, see column 3 and 4 for bands 1 and 2 in Table S1). However, we observe an effect on the intensity which changes the appearance of the spectrum. For the second intense band, 4, with origin in 2T2 , a larger direct solvent effect is found. Both the splitting and the relative intensities between the three 2T2 states

21



changes. Due to the former, the spectrum now has two peaks with maxima at 43600 and 44200 cm−1 (5.41 and 5.48 eV) and thus, the solvent shifts the energy with up to ∼ 0.2 eV. We can now combine the total shift due to relativity and solvent effects: Here we obtain for the first intense band (2) in Figure 1 a total shift of about 0.4 eV, of which approximately 0.2 eV is due to relativistic effects, and 0.2 eV is an indirect solvent shift, caused by a change in the TcO4– structure. The electrostatic effects (i.e. direct solvent shifts) are in this case negligible. For the second intense band (4), the direct solvent shift is larger and the total shift is ∼ 0.4 eV from 5.22 eV in vacuum (non-relativistic) to 5.41–5.58 eV in water (relativistic). The relativistic effect is approximately 0.2 eV, whereas the indirect solvent shift in this case is negligible, but the direct solvent effects contribute with about 0.2 eV. We now move on to discuss ReO4– . The spectra that we have obtained for ReO4– are depicted in Figure 1 (e)–(h) in the same manner as for TcO4– , and the numerical values for the peak-maxima are given in the SI (Table S2). With the heavier Re nucleus, we expect larger relativistic effects. Indeed, comparing the first two vacuum calculations (nonrelativistic and relativistic) in Figure 1 (e) and (f), we see that the two bands (2 and 4) are shifted significantly (∼0.3–0.5 eV). We also note that two new transitions appear, splitting the bands. Analyzing the states under each band shows that both bands 2 and 4 are of T2 parentage. In band 2, the first peak can be ascribed to 13 T2 whereas the second peak has origin in 11 T2 , although both peaks show pronounced singlet-triplet mixing (which was not seen to the same degree for TcO4– ). This mixing also explains the increased intensity of the bands with origin in 13 T2 . The transitions within band 4 show a similar pattern, cf. Figure 1 (b). Using the structure obtained within the environment, leads (as for TcO4– ) to the appearance of a small band immediately before the large band from the 11 T2 and 13 T2 states. Also in ReO4– , the band has origin in transitions of 11 T1 parentage that becomes partially allowed in the reduced symmetry. Unlike TcO4– the band has significant contributions from

22


Page 22 of 35

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


0.060

3.7

eV 5.0

4.3

5.6

6.2

(a)

Vacuum (Sym, nonrel) 0.050

6.2

6.8

(e)

4 0.040

0.030

2

2

0.020

0.010

21T2

11T2

0.030

1 1T 2

eV

0.050

21T2

0.020

5.6

Vacuum (Sym, nonrel)

4

0.040

5.0

0.060

0.010

25000

30000

35000

3.7

0.050

40000 cm −1

45000

50000

5.6

6.2

eV 5.0

4.3

30000

0.020

(b)

Vacuum (sym)

35000

4.3

40000

5.0

45000 cm −1

50000

eV 5.6

6.2

6.8

(f)

Vacuum (sym)

4

0.040

0.015

21T2

21T2

11T2

13T2

0.030

23T2

0.010

2

11T2

0.020

25000

30000

3.1

3.7

4

2

33T2

0.005

13T2

0.010

0.020

55000

35000

40000 cm −1 eV

4.3

Vacuum

45000

50000

5.6

6.2

5.0

30000

0.020

(c)

1

2 T2

40000

4.3

5.0

eV

0.010

6.8

4 2 1T 2

1 13T2 1 T2

3 3 T2

31T1

1 T2

0.005

3 11T1

1

30000

35000

3.7

4.3

0.020

6.2

5.6

2

1

25000

55000

0.015

13T2

0.005

50000

(g)

2 0.010

45000 cm −1

Vacuum

4

0.015

35000

40000 cm −1 eV 5.0

45000

50000

5.6

6.2

30000

0.020

(d)

M2P2

1 1 T1

1

35000

40000

4.3

5.0

23T2

3

45000 cm −1

50000

eV 5.6

6.2

55000

6.8

(h)

M2P2

1

2 T2 0.015

4

0.015

2 2

0.010

4 0.010

13T2 1 T2

0.005 1

1 T1 25000

30000

21T2

1

1 T2 1 3 T2

1

31T1 3

0.005

3

1 35000

40000 cm −1

11T1 45000

50000

30000

23T2/33T2

35000

1 40000

45000 cm −1

50000

55000

Figure 1: Calculated spectra of TcO4– [left column: (a)–(d)] and of ReO4– [right column: (e)–(h)]. 23



the triplet, 13 T2 , that splits sufficiently to contribute to both band 1 and 2. We now consider the indirect solvent effect, which is obtained by comparing Figure 1 (f) with Figure 1 (g). The band systems 2 and 4 are not much affected (for the numerical values, compare column 2 and 3 in Table S2). It should be noted that a third band (3) appears between these two bands (cf. Figure 1 g and h), which from analysis of the response vectors is shown to have origin in 31 T1 (the 21 T2 transitions are hidden under the 11 T2 /13 T1 transitions). The same band is seen for TcO4– , cf. Figures 1 (c) and (d), where it appears as a shoulder to band 4. In the case of direct solvent effect, we note that it is most pronounced for band 4, where the effect amounts to ∼ 0.1 eV, cf. Figures 1 (f) and (g). As for TcO4– we can now estimate the total effect of the relativistic and solvent effects to be ∼ 0.6 eV for band 2 and ∼ 0.7 eV for band 4. For the former, both indirect and direct solvent effects are negligible. In the latter case, the direct solvent effect account for about 0.1 eV. The calculated intensities are also remarkably affected by the solvent, which significantly alter the features of the spectrum.

External Effective Field Effects We finally probe the effect of including the EEF effects described in the last part of the theory section. The calculated spectra are shown in Figure 2 in black (for easier comparison, the spectra obtained without involving EEF are shown in red). It should be noted that the excitation energies are not affected by including the EEF, although the peak maxima (which also depend on the intensities) can be affected. For TcO4– , both bands 2 and 4 (11 T2 and 21 T2 ) were significantly split without considering EEF effects. The EEF reduces the intensity of the highest of the transitions of 11 T2 parentage, which was the main contributor to the second peak within band 2. At the same time, including EEF effects increases the intensity related to some of the transitions of 13 T2 parentage. The same scenario repeats for band 4, for which the two peaks also reduces to one as can be seen by comparing the red 24


Page 24 of 35

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


0.020

3.7

eV 5.0

4.3

5.6

2

4.3

eV 5.6

5.0

6.2

4

0.015

2 1T 2

21T2

2

2 1T 2

0.010

1

1 T2

6.8

(b)

M2P2 M2P2(EEF)

4

0.015

0.010

0.020

(a)

M2P2 M2P2(EEF)

13T2

6.2

23T2 0.005

0.005

25000

30000

35000

40000 cm −1

23T2

3

3

1 45000

50000

30000

35000

40000

45000 cm −1

50000

55000

Figure 2: Calculated spectra of (a) TcO4– and (b) ReO4– (b). Only peaks with intensities very affected by inclusion of EEF effects are assigned. and black spectrum in Figure 2 (a). In ReO4– , the largest EEF effect is for band 4, where the intensity is re-distributed among the peaks of 21 T2 parentage, and the same occurs for the peaks of 23 T2 parentage. The final result is a small (red) shift of the peak maxima of band 4. At this point it should be emphasized that the peak is only red-shifted in the spectrum, due a change of intensity in the underlying transitions. The EEF model does not alter the transition frequencies, which are identical in the spectra. As previously stated, the results presented here are only valid for the specific snapshot, and will most likely change for other snapshots. The above discussion thus only serves to show that the EEF effects can alter the appearance of the spectrum significantly, by decreasing or increasing the intensities of certain transitions.

Discussion and Conclusion We have in this paper shown the formal equations and presented the first implementation of the PE model within fully relativistic, 4c-DFT and TD-4c-DFT frameworks. As an illustration, the new PE-4c-TD-DFT method was employed to estimate solvent and relativistic effects on the absorption spectra of TcO4– and ReO4– . By this, we have been able to compare both direct electrostatic effects due to inclusion of a solvent potential, indirect structural ef25



Page 26 of 35

fects from structural distortions caused by the solvent, and effects due to the relativistic treatment. 0.020

4.3

5.0

eV

5.6

6.2

0.030

(a)

M2P2 M2P2(nonrel)

0.025

4.3

5.0

eV 5.6

6.2

6.8

(b)

M2P2 M2P2(nonrel)

0.015 0.020 0.010

0.015 0.010

0.005 0.005

25000

30000

35000

40000 cm −1

45000

50000

30000

35000

40000

45000 cm −1

50000

55000

Figure 3: Effect of the 4c Dirac-Coulomb method versus a non-relativistic treatment for the calculated spectra of (a) TcO4– and (b) ReO4– in aqueous solvent (M2P2). For excitation energies, we find that both indirect and direct solvent effects can contribute to a shift, both up to around 0.1–0.2 eV in TcO4– and ReO4– , each. In the investigated cases, the indirect and direct solvent effect both lead to blue shifts, and therefore do not cancel out. Furthermore, the magnitudes of both indirect and direct solvent effects differ between the individual transitions, even within the same system. Accordingly, a uniform shift of the spectrum is not sufficient to account for these effects. The use of a relativistic framework also shifts the energies; the amount of the shift is (as expected) strongly correlated with the nature of the metal. In TcO4– , the relativistic effect amounts to a shift of up to 0.2 eV, which is comparable to the effects from the solvent, whereas for ReO4– , the relativistic effects are much more substantial (up to 0.6–0.7 eV). In both cases the effects on the calculated spectra is thus significant. In particular it should be emphasized that the impact of spinorbit coupling is significant for ReO4– , where singlet and triplet states are significantly mixed. This effect is not obtained with scalar relativistic approaches or if the relativistic effects are described exclusively through a basis set with a pseudopotential. Spectra calculated within the PE model for both non-relativistic and relativistic models are shown in Figure 3.

26


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


At this point we can also high-light a few differences between this work and the interesting X2C-PERI-CC2 study by Krause et al. 33 Our current implementation includes all the X2C Hamiltonians implemented in DIRAC, although we have not employed any X2C Hamiltonian in this paper. While the computational cost of X2C models is considerably lower than for 4c models, the scaling with system size is the same, and the theory and implementation presented here enables the choice of a 4c Hamiltonian when desired. One cannot directly compare the two studies, as they employ different electronic structure methods, and have different target systems. Ref. 33 targets small organic systems, where the solvated molecule shows large solvent effects, but rather benign relativistic effects (unlike the systems studied here). Further, the descriptions of the environment are very different. Ref. 33 studies an average solvent shift over a number of snapshots from a classical MD simulation. The inclusion of the dynamics of the solvent-solute interactions is indeed an interesting aspect of Ref. 33, but has not been done in our study since MD simulations for transition metal systems are much more involved than for small organic systems. The embedding potential in Ref. 33 employs a rather simplified electrostatic description of the environment, where we have instead chosen to focus on a systematic study of different types of solvent effects with accurate solvent potentials. Finally, we will make a few comments on our results in comparison to the experimental spectra, although the comparison is naturally hampered by the fact that we only include one structure. The experimental spectra have rather broad transitions with marked vibrational structure. The (overlapping) bands are in the regions 32300–38700 cm−1 (∼ 4.0–4.8 eV) and 37900–43500 (∼ 4.7–5.4 eV) for TcO4– . 77,84 For ReO4– , the corresponding transitions are 40300–46000 (∼ 5.0–5.7 eV) and 46000–54000 (∼ 5.7–6.7 eV). 77,84 We see that our results are in reasonable agreement with these numbers. Furthermore, it has been suggested 77 that triplet excitations contribute significantly to the two bands in ReO4– , and this we can confirm.

27



Page 28 of 35

Acknowledgement E.D.H. (grant id. CF15-0208,) and J.M.H.O. (grant id. CF15-0823) thank the Carlsberg foundation. J.M.H.O. and J.K. thank the Danish Council for Independent Research (Sapere Aude research career program, grant id. DFF–1325-00091B for J.M.H.O. and support within the Sapere Aude program for J.K.). E.D.H and J.K. thank the Villum Foundation for financial support. The authors thank DeIC for computational resources.

Supporting Information Available Peak-maxima for bands 1–4 for TcO4– and ReO4– are provided in Table S1 and S2.

This

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

References (1) Kupratakuln, S. J. Phys. C 1970, 3, S109–S119. (2) Pyykkö, P.; Desclaux, J. P. Acc. Chem. Res. 1979, 12, 276–281. (3) Pyykkö, P. Chem. Rev. 1988, 88, 563–594. (4) Gorin, D. J.; Toste, F. D. Nature 2007, 446, 395–403. (5) Ahuja, R.; Blomqvist, A.; Larsson, P.; Pyykkö, P. Phys. Rev. Lett. 2011, 106, 018301. (6) Sieffert, N.; Wipff, G. Dalton Trans. 2015, 44, 2623–2638. (7) Saue, T. ChemPhysChem 2011, 12, 3077–3094. (8) Reiher, M.; Wolf, A. The Theory of Intermolecular Forces; Wiley-VCH, 2009. (9) Hu, L.; Söderhjelm, P.; Ryde, U. J. Chem. Theory Comput. 2013, 9, 640–649. (10) Finkelmann, A. R.; Senn, H. M.; Reiher, M. Chem. Sci. 2014, 5, 4474–4482. 28


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


(11) Hedeg˚ ard, E. D.; Kongsted, J.; Ryde, U. Angew. Chem. Int. Ed. 2015, 54, 6246–6250. (12) Crescenzi, O.; Pavone, M.; de Angelis, F.; Barone, V. J. Phys. Chem. B 2005, 109, 445–453. (13) Pedone, A. J. Chem. Theory Comput. 2013, 9, 4087–4096. (14) Thallmair, S.; Zauleck, J. P. P.; de Vivie-Riedle, R. J. Chem. Theory Comput. 2015, 11, 1987–1995. (15) Olsen, J. M.; Aidas, K.; Kongsted, J. J. Chem. Theory and Comput. 2010, 6, 3721– 3734. (16) Olsen, J. M.; Kongsted, J. Adv. Quant. Chem. 2011, 61, 107–143. (17) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999–3093. (18) Di Remigio, R.; Bast, R.; Frediani, L.; Saue, T. J. Phys. Chem. A 119, 119, 5061–5077. (19) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227–249. (20) Senn, H. M.; Thiel, W. Angew. Chem. Int. Ed. 2009, 48, 1198–1229. (21) Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. J. Chem. Phys. 2011, 134, 104108. (22) Schwabe, T.; Sneskov, K.; Olsen, J. M. H.; Kongsted, J.; Christiansen, O.; Hättig, C. J. Chem. Theory Comput. 2012, 8, 3274–3283. (23) Eriksen, J. J.; Sauer, S. P. A.; Mikkelsen, K. V.; Jensen, H. J. A.; Kongsted, J. J. Comput. Chem. 2012, 33, 2012–2022. (24) Hedeg˚ ard, E. D.; List, N. H.; Jensen, H. J. Aa.; Kongsted, J. J. Chem. Phys. 2013, 139, 044101.

29



Page 30 of 35

(25) Pedersen, M. N.; Hedeg˚ ard, E. D.; Olsen, J. M. H.; Kauczor, J.; Norman, P.; Kongsted, J. J. Chem. Theory Comput. 2014, 10, 1164–1171. (26) Steinmann, C.; Olsen, J. M. H.; Kongsted, J. J. Chem. Theory Comput. 2014, 10, 981–988. (27) Hedeg˚ ard, E. D.; Olsen, J. M. H.; Knecht, S.; Kongsted, J.; Jensen, H. J. Aa. J. Chem. Phys. 2015, 114, 1102–1107. (28) Steindal, A. H.; Beerepoot, M. T. P.; Ringholm, M.; List, N. H.; Ruud, K.; Kongsted, J.; Olsen, J. M. H. Phys. Chem. Chem. Phys. 2016, 18, 28339–28352. (29) List, N. H.; Olsen, J. M. H.; Kongsted, J. Phys. Chem. Chem. Phys. 2016, 18, 20234– 20250. (30) Gomes, A. S. P.; Jacob, C. R.; Visscher, L. Phys. Chem. Chem. Phys 2008, 10, 5353– 5362. (31) Gomes, A. S. P.; Christoph, C. R.; Réal, F.; Visscher, L.; Vallet, V. Phys. Chem. Chem. Phys. 2013, 15, 15153–15162. (32) Krause, K.; Klopper, W. J. Chem. Phys. 2016, 144, 041101. (33) Krause, K.; Bauer, M.; Klopper, W. J. Chem. Theory Comput. 2016, 12, 2853–2860. (34) Olsen, J. M. H. PElib,

The Polarizable Embedding library,

https://gitlab.com/pe-software/pelib-public, 2016. (35) Saue, T.; Jensen, H. J. Aa. J. Chem. Phys. 1999, 111, 6211–6222. (36) Saue, T.; Helgaker, T. J. Comput. Chem. 2002, 23, 814–823. (37) Saue, T.; Jensen, H. J. Aa. J. Chem. Phys. 2003, 118, 522–536.

30


version 1.2.4.

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


(38) DIRAC, a relativistic ab initio electronic structure program, Release DIRAC16 (2016), written by H. J. Aa. Jensen, R. Bast, T. Saue, and L. Visscher, with contributions from V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, A. S. P. Gomes, T. Helgaker, J. Henriksson, M. Iliaˇs, Ch. R. Jacob, S. Knecht, S. Komorovsk´ y, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, H. S. Nataraj, M. K. Nayak, P. Norman, G. Olejniczak, J. Olsen, Y. C. Park, J. K. Pedersen, M. Pernpointner, R. di Remigio, K. Ruud, P. Salek, B. Schimmelpfennig, J. Sikkema, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, S. Villaume, O. Visser, T. Winther, and S. Yamamoto (see http://www.diracprogram.org). (39) Kagawa, T. Phys. Rev. A 1975, 12, 2245–2256. (40) Malli, G.; Oreg, J. J. Chem. Phys. 1975, 63, 830–841. (41) Aoyama, T.; Matsuoka, O. J. Chem. Phys. 1980, 73, 1329–1333. (42) Rajagopal, A. K.; Callaway, J. Phys. Rev. B 1973, 7, 1912–1919. (43) Talman, J. D. Phys. Rev. Lett. 1986, 57, 1091. (44) LaJohn, L.; Talman, J. D. Chem. Phys. Lett. 1992, 189, 383. (45) Olsen, J.; Jørgensen, P. J. Chem. Phys. 1985, 82, 3235–3264. (46) Olsen, J.; Jørgensen, P. In Modern Electronic Structure Theory; Yarkony, D. R., Ed.; World Scientific, 1995; Vol. 2; pp 857–990. (47) Salek, P.; Helgaker, T.; Saue, T. Chem. Phys. 2005, 311, 187–201. (48) Bast, R.; Jensen, H. J. Aa.; Saue, T. Int. J. Quantum Chem. 2009, 109, 981–988. (49) Salek, P.; Vahtras, O.; Helgaker, T.; ˚ Agren, H. J. Chem. Phys. 2002, 117, 9630–9645. (50) Christiansen, O.; Jørgensen, P.; Hättig, C. Int. J. Quantum Chem. 1998, 68, 1–52.

31



(51) Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristiansen, K.; Olsen, J.; Ruud, K. Chem. Rev. 2012, 112, 543–631. (52) Mikkelsen, K. V.; Jørgensen, P.; Jensen, H. J. Aa. J. Chem. Phys. 1994, 100, 6597– 6607. (53) Cammi, R.; Mennucci, B. J. Chem. Phys. 1999, 110, 9877–9877. (54) Cammi, R.; Cossi, M.; Tomasi, J. J. Chem. Phys. 1996, 104, 4611–4620. (55) Cammi, R.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 1998, 102, 870–875. (56) Wortmann, R.; Bishop, D. M. J. Chem. Phys. 1998, 108, 1001–1007. (57) Pipolo, S.; Corni, S.; Cammi, R. J. Chem. Phys. 2014, 140, 164114. (58) Jensen, L.; Swart, M.; van Duijnen, P. T. J. Chem. Phys. 2005, 122, 034103. (59) List, N. H.; Jensen, H. J. Aa.; Kongsted, J. Phys. Chem. Chem. Phys. 2016, 10070– 10080. (60) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165–169. (61) Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571–2577. (62) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119–124. (63) Perdew, J. P. Phys. Rev. B 1986, 33, 8822. (64) Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100. (65) Olsen, J. M. H.; Steinmann, C.; Ruud, K.; Kongsted, J. J. Phys. Chem. A 2015, 119, 5344–5355.

32


Page 32 of 35

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


(66) Adamo, C.; Scuseria, G. E.; Barone, V. J. Chem. Phys. 1999, 111, 2889–2899. (67) Dunning Jr., T. H. J. Chem. Phys. 1989, 90, 1007–1023. (68) Dunning Jr., T. H. J. Chem. Phys. 1971, 55, 716–723. (69) Dyall, K. G. Theor. Chem. Acc. 2007, 117, 483–489. (70) Dyall, K. G. Theor. Chem. Acc. 2004, 112, 403–409. (71) Dyall, K. G.; Gomes, A. S. P. Theor. Chem. Acc. 2009, 125, 97. (72) Söderhjelm, P.; Husberg, C.; Strambi, A.; Olivucci, M.; Ryde, U. J. Chem. Theory Comput. 2009, 5, 649–658. (73) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–789. (74) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (75) Gagliardi, L.; Lindh, R.; Karlström, G. J. Chem. Phys. 2004, 121, 4494–4500. (76) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P., P. ˚ A. Neogrády; Pedersen, T. B.; Pitonak, M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. J. Comput. Chem. 2010, 31, 224–247. (77) G¨ udel, H. U.; Ballhausen, C. J. Theor. Chem. Acc. 1972, 25, 331–337. (78) Zielger, T. Struct. Bond. 2012, 143, 1–38. (79) Su, J.; Xu, W.-H.; Xu, C.-F.; Schwarz, W. H. E.; Li, J. Inorg. Chem. 2013, 52, 9867– 9874. (80) Hasegawa, J.; Toyota, K.; Hada, M.; Nakai, H.; Nakatsuji, H. Theor. Chim. Acc. 1995, 92, 351–359. (81) Ziegler, T.; Rauk, A.; Baerends, E. J. Chem. Phys. 1976, 16, 209–217. 33



(82) Xu, W.; Ma, J.; Peng, D.; Zou, W.; Liu, W.; Staemmler, V. Chem. Phys. 2009, 356, 219–228. (83) Bailey, N.; Carrington, A.; Lott, K. A. K.; Symons, M. C. R. J. Chem. Soc. 1958, 290–297. (84) Mullen, P.; Schwochau, K.; Jørgensen, C. K. Chem. Phys. Lett. 1969, 3, 49–51. (85) Neugebauer, J.; Baerends, E. J.; Nooijen, M. J. Phys. Chem. A 2005, 109, 1168–1179.

34


Page 34 of 35

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

4.3

5.0

eV

5.6

6.2

0.030

TcO4–

M2P2 M2P2(nonrel)

0.025

4.3

5.0

eV 5.6

6.2

6.8

ReO4–

M2P2 M2P2(nonrel)

0.015 0.020 0.010

0.015 0.010

0.005 0.005

25000

30000

35000

40000 cm −1

45000

50000

30000

35

35000

40000

45000 cm −1


50000

55000

Relativistic Polarizable Embedding - Journal of Chemical Theory and

Recommend Documents