Advanced techniques: Quantum treatment of nuclei and non-adiabatic - - PowerPoint PPT Presentation

1

Advanced techniques: Quantum treatment of nuclei and non-adiabatic (surface hopping) approaches

Mauro Boero

Institut de Physique et Chimie des Matériaux de Strasbourg University of Strasbourg - CNRS, F-67034 Strasbourg, France and @Institute of Materials and Systems for Sustainability, Nagoya University - Oshiyama Group, Nagoya Japan

About Excited States

• Photoactive molecules are also the target of potential

technological applications in molecular optoelectronics, photocatalysis and photo-biochemistry

• They involve electron excitations
• Time-dependent DFT (TDDFT) has been proposed as a

way to include electron excitation (see M. E. Casida, Recent Advances in Density Functional Methods, Vol. 1,

• ed. by Chong, D.P., World Scientific, Singapore, 1995)
• Although TDDFT is computationally expensive…

2

• Generally organic photoreactions involve mainly the first

excited singlet state (S1) and the lowest triplet state (T1). Other excited states have a too short lifetime to be of real practical interest and can generally be neglected

• see N. J. Turro, Modern Molecular Photochemistry,

University Science Books, Mill Valley1991

• A. Zewail* Femtochemistry: Ultrafast Dynamics of the

Chemical Bond, World Scientific Series in 20th Century Chemistry, Vol. 3, World Scientific, Singapore 1994 * 1999 Chemistry Nobel Prize

3

About Excited States

The “minimal” excited states

• If a single valence electron is excited from the highest
• ccupied (ground state) orbital a to the lowest unoccupied
• rbital b, four different determinants can be obtained

according to the Pauli’s principle

4

The “minimal” excited states: Wavefunctions

5

The “minimal” excited states: Energies

6

ROKS: Excited states@KS

• Instead of separate wavefunctions for the t and m states, it

has been shown that it is possible to determine a single set

• f spin-restricted single-particle orbitals yi(x) for the states

i = 1,…, N+1 in such a way that

• A new DFT functional, the restricted open shell Kohn-

Sham (ROKS) functional, can be written as

7

) ( ) ( ) ( ) ( ) ( x x x x x

t t m m    

        

     

 

 

 

   

1 1 , * 3 KS KS ROKS

) ( ) ( 2 )} ( {

N j i ij j i ij t m i

x d E E H  y y    y x x x

• I. Frank, J. Hutter, D. Marx, M. Parrinello, J. Chem. Phys. 108, 4060 (1998)

Total energy functionals for excited states

• The functionals with the superscript KS are Kohn-Sham total

energy functionals with the difference reduced only to the exchange-correlation term

8

       

II eI m m xc H i k i m

E E E E E E          y y

 

] , [ } { } {

KS

       

II eI t t xc H i k i t

E E E E E E          y y

 

] , [ } { } {

KS

…and associated equations to solve (I)

• The minimization of the functional HROKS[yi(x)] with respect

to the orbitals leads to two sets of Schrödinger-like equations,

• ne for the doubly occupied orbitals…

9

   

        



m m xc m m xc eI

v v V y d

     

     , , ) ( ) ( 2 1

3 2

x y

• x

y

   



 

     

1 1

) ( ) ( , 2 1 , 2 1

N j j ij i t t xc t t xc

v v x x y  y    

     

t m t m t m xc xc t m t m t m xc xc

E v E v

, , , , , ,

] , [ , ] , [

       

           

• …and one for the singly occupied a and b states

10

   

) ( , 2 1 , ) ( ) ( 2 1 2 1

3 2

x x y

• x

y

a t t xc m m xc eI

v v V y d y     

     

                     





 



1 1

) (

N j j aj

x y 

   

) ( , 2 1 , ) ( ) ( 2 1 2 1

3 2

x x y

• x

y

a t t xc m m xc eI

v v V y d y     

     

                     

 

 



1 1

) (

N j j aj

x y 

…and associated equations to solve (II)

11

…and does it actually work?

12

…well, maybe yes !

Where can it be used ?

• This approach has been used to study the isomerization and

energy changes of the rhodopsin chromophore.

• This is the photosensitive protein in the rod cells of the retina
• f vertebrates and the process of vision.
• It, involves the photoisomerization

as a response to the absorption of photons (in about 200 fs) and triggers a cascade of slower reactions that produce a specific biological signal (C. Molteni et al. JACS 121, 12177 (1999)).

13

ћw

14

Rhodopsin: Light-sensitive receptor protein responsible for vision

15

16

Rhodopsin: ROKS simulated isomerization

Doing CPMD-like dynamics with more than one PES

17

• If the ground state S0 and the ROKS excited S1 surface are two

accessible states (e.g. photochemistry) it is possible to adopt a Tully scheme (J. C. Tully, J. Chem. Phys. 93, 1061 (1990); ibid. 55, 562 (1971))

• The electronic wavefunction for the whole system reads



 

   

excit j

N j dt E i j j

e a

Adiabatic state on Sj Energy expectation value on Sj

Doing CPMD-like dynamics with more than one PES

18

S1 S0

(adiabatic) surfaces between which we want to hop

F1 F0

j j j

H ˆ E   

..and aj?

Doing CPMD-like dynamics with more than one PES

19

• aj are determined by the solution of the time-dependent

Schrödinger equation

• For the closed shell KS ground state S0
• And for the excited ROKS S1
• n = half the (even) number of electrons

     t i H ˆ

1 1 * n n *

      

 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

2 1

 

          

n * n * * n n *

 

Doing CPMD-like dynamics with more than one PES

20

• These 0 and 1 are normalized on S0 and S1, respectively but

they are not orthogonal to each other

• It is possible to define the quantities

Non-adiabatic coupling matrix → easy to do: wfs velocities are directly available in CPMD

ij j i ij

S     

1

10 01

  

ii

S S S S

j i j i ij

t D         

  

ii ji ij

D D D

Doing CPMD-like dynamics with more than one PES

21

• Solving for gives

       t i H ˆ

* *

 

 

   



dt E i a

j j j j 1

exp

 

           S D a p p D a E E S p p ia S a

10 1 01 1 1 1 1 2

1 1 

 

         

1 2 1 01 1 1 10 2 1

1 1 E E S ia S D a p p D a S a 

Doltsinis & Marx, Phys. Rev. Lett. 88, 166402 (2002)

Doing CPMD-like dynamics with more than one PES

22

• Note that
• If the wavefunctions were eigenfunctions of the KS Hamiltonian,

then |a0|2 and |a1|2 would be occupation numbers

• …but they are not. So what ?

 



dt E i j

j

e p

Doltsinis & Marx, Phys. Rev. Lett. 88, 166402 (2002)

S E H H H H ˆ E

jj j j j 10 01

     

Doing CPMD-like dynamics with more than one PES

23

• Expand on an orthonormal auxiliary set of wfs ’j

j j j

p a b b b d d          

1 1 1 1

1

2 1 2

  d d

true state population

1 1

    

j j j

c c

                 

11 01 1 10 00

c c c c c c

Eigenvectors of H ci = Ei S ci

Doing CPMD-like dynamics with more than one PES

24

• Hence, we get
• and the orthonormal auxiliary wavefunctions and occupations

> E1 if E0 < E1

                  1 1

1

S c c

2 2 1 1

1 S E S E E E E      

2 1 1

1 S S          

 

1 2 1 2 2 2

2 b b Re S b S b d

*

  

 

2 1 2 2 1

1 b S d  

Doing CPMD-like dynamics with more than one PES

25

• Transitions from one surface to the other is done according to the

fewest switches criterion of Tully: select a random number rnd[]=z. If z > pi then

• jump from i to j if Ej-Ei < Ti = kinetic energy in the i state
• and do not hop if Ti is insufficient to compensate for the jump

Ej - Ei → Tully’s forbidden transitions

• Warning: The accessible time scale is of the order of fs (as in TD-

DFT, Ehrenfest MD, etc.) Dt = CPMD time step

2 2 i i i

d dt / d d t p D  

ћw

Example of application: Photo-isomerization of formaldimine

26

Example of application: Photo-isomerization of formaldimine

27

Doltsinis & Marx, Phys. Rev. Lett. 88, 166402 (2002)

Beyond classical nuclei: Path Integral MD

28

Sometimes nuclei cannot be approximated as classical point-like

• particles. Textbook examples are:
• Tunneling processes (mainly protons)
• Quantum broadening of nuclei comparable to (or larger than)

thermal fluctuations Proposed solution: Consider also (light) nuclei as quantum objects First Principles Path Integral Molecular Dynamics (PIMD)

Path Integrals ?

29

A quantum particle travelling from an initial point x0 at time t0 can reach a point x’ at time t’ along any possible path joining the two points… …each weighted with its own probability exp (-iS), being S the action of the system.

(R.P. Feynman, Statistical Mechanics, Addison Wesley, Massachusetts, 1972)

x x0 x’ t  3 2 x1 x2 x3 R.P. Feynman 1918 - 1988

) , ( x x L dt S 





Well, yes… Path Integrals !

The probability amplitude of finding at (x’, t’) a particle started at (x0, t0) is then a sum of all possible weighted paths      

 

        



    



n n n

L dt i L dt i L dt i n n H t t i

e e e dx dx dx t x e t x

) 1 ( 1 2 2 1

... ... ' , ' ,

1 2 1 ˆ ) ' (

lim

where ne = t’, x0= x(t0) and x’ = x(t’) = x(ne. Formally we write

1 2 1 1 1

... lim lim ) (

      

 



n n n i i n

dx dx dx dx t x D

 

) ( ), ( ˆ ) ' (

) ( ' , ' ,

t x t x iS x x H t t i

e t x t x e t x

    



 D

30

Path Integral formulation of the partition function

The partition function Z() of a system having a Hamiltonian H = K + V , i.e. a kinetic term K= mv2/2 (depending on the velocities) and a potential term V (depending on positions) at a temperature T ( = 1/kBT) can be written as

 

1 1 1 1 1 1

lim Tr ) ( x x dx x e x dx e Z

P P V K H

   

 

      



31

where we make use of the Trotter-Suzuki factorization

 

P P P V P K P V P P V K

e e e e         

        

lim lim

2 2    

H.F. Trotter Proc. Amer. Math. Soc. 10 545 (1959)

• M. 鈴木, 物性研究 4, 415 (1965)

Path Integral formulation of the partition function

Using the completeness relation we get

1 1 1 1 1 3 2 2 1 2 1

lim ... ... lim ) (

       

     

 

P i i P i i P P P P

x x dx x x x x x x dx dx dx Z 

32

x x dx



 1

and each matrix element can be computed as

 

             

 

          

) ( ) ( 2 ) ( 2 1 2 2 1

1 2 1

2

i i i i

x V x V P x x mP i V P K P V P i i i

e mP x e e e x x x

    



Path Integral formulation of the partition function

33

Hence, the partition function reads (wP

2=P/ 2)

 

 

) ( 1 ) ( 1 ) ( 2 1 2 / ) ( ) ( 2 1 2 /

1 1 1 2 1 2 1 1 1 2 1 1 1

) ( 2 lim ... 2 lim

 w   

   

x S x x x V P x x m x x P i i P P x V P x x mP x x P P P

e x dx e dx mP e dx dx mP Z

P i i i i P P P i i i i P

                         

    

                    

     

D

map each QM particle onto an effective classical system of P beads coupled by a harmonic potential plus V(x)

Path Integral formulation in a Car-Parrinello Lagrangean scheme

34

 

 

 

   

     

 

                   

P I I I P I I I P i j i ij j i ij I i DFT i

M M E P

1 2 2 2 1 , 2

2 2 , ) ( 1

          

w  y y y y  R R R x   L

• M and M’ are different in order to keep on the same time

scale all the modes involved in the QM treatment of the nuclei

• A Nosé-Hoover thermostat is added to ensure a canonical

sampling and to control the adiabaticity

• D. Marx and M. Parrinello,
• Zeit. Phys. B 95, 143 (1994)

The Euler-Lagrange equations of motion:

35

) ( ) ( 1 1 ) ( 1

1 *

x x x

     

y   y y  y 

i j j ij i DFT i

P P E P P         



   

  y  

2 1 ) 1 ( 2 1 ) 1 (

) ( 2     

e i kin e i e

Q E Q         



x

 

 

M e e e e

Q Q Q

           

         

   

1 ] [

1 ) ( 1 2 1 ) 1 ( ) (

    

with  = 2,…,M and for the QM nuclei we have

 

        

w



I I I I I I P DFT I I

s M M E P M

I

R R R R R

R

   

1 1 1 2

2 1         

 

and related thermostats

Example of application of PI-CPMD: proton propagation in water

36

The proton is no longer a point-like object RH(t) but a (quantum) probability distribution (RH). Eigen/Zundel transition between oxygen Oa and Ob

• M. Tuckerman et. al. Science 275, 817 (1997)
• D. Marx et al. Nature 397, 601 (1999)

ROa ROb

Example of application of PI-CPMD: Proton propagation in water

37

Probability distributions of proton position in terms of H+ displacement  = |ROa-ROb| relative to the Oa-H-Ob midpoint. Quantum Classical

Proton propagation in water

Free energy profile as a function of : no single dominant structure

38

• 1.0 0.0 +1.0

 = |ROa-ROb| DF (kcal/mol) 0 1.0 2.0

300 K Quantum Classical

39

Acknowledgements

• Michele Parrinello, ETHZ-USI
• Roberto Car, Princeton University
• Kiyoyuki Terakura, NIMS
• Atsushi Oshiyama, Nagoya University
• Michiel Sprik, Cambridge University
• Pier Luigi Silvestrelli, Padova University
• Alessandro Laio, SISSA
• Marcella Iannuzzi, PSI
• Wanda Andreoni, Swiss Federal Institute of Technology (EPFL)
• Carla Molteni, King’s College
• Kenichi Koizumi, IMS-Okazaki

Special thanks to Katsumasa Kamiya (Kanagawa Institute of Technology), Ikeda Takashi (JAEA) Carlo Massobrio & Guido Ori (IPCMS) for a critical reading of the original draft and for correcting several mistakes