Transcript lecture9 - WordPress.com

Ultrafast processes in molecules
IX – Surface hopping implementation
Mario Barbatti
[email protected]
Energy
Wave packet propagation
Surface hopping propagation
Reaction coordinate
• Tully, Preston, JCP 55, 562 (1971)
2
the semi-classical propagation
3
 r, R ,t  
c

c j t Φ j r; R
c
t 
 Φ  |
k
Φk Φl
r
  kl 
j



i

H
 (r , R , t )  0
e 

 t

dc k
i

dt
2
dt
 Φ  | H
k
 Φ  | H
d
k
2
kj
kj
c
c
c
 cj  0
Φk
j
c
d Rm
   H kj  i Fkj  v
c ,m
F kj
c

Fm
0
Mm
t
c
 F kj  v
c
Φj
r
 Φk  Rm Φ j
A  AR
 V k  kj 
 W kj , F kj  0

c
c
r
t 
Adiabatic basis

Diabatic basis
4
2
c
d Rm
dt
2
c

Fm
Fm   R m H ll
c
0
Mm
c
Any standard method can be used in the integration of the Newton’s
equations.
A good one is the Velocity Verlet
For each nucleus m: R cm ( t   t )  R cm ( t )  v cm ( t )  t 
1
2
a m (t )  t
c
2
t 
1 c

c
v t 
  v m (t )  a m (t )  t
2
2


c
m
a m (t )  
c
1
Mm
c
 R E  R m ( t   t ) 
t  1 c
c
c 
v m (t   t )  v m  t 
  a m (t   t )  t
2  2

• Swope et al. J. Chem. Phys. 76, 637 (1982)
5
• Schlick, Barth and Mandziuk, Annu. Rev. Biophys. Struct. 26, 181 (1997)
6
Time step should not be larger than 1 fs (1/10v).
Dt = 0.5 fs assures a good level of conservation of energy.
Exceptions:
• Dynamics close to the conical intersection may require 0.25 fs
• Dissociation processes may require even smaller time steps
7
i
dc k
dt

 H
c
kj
 i Fkj  v
c
c
c
j
0
SC-TDSE
j
The SC-TDSE is solved with standard methods
(Unitary Propagator, Adams Moulton 6th-order, Butcher 5th-oder)
0.7
t = 0.5 fs; ms = 1
Runge-Kutta
Adams Moulton
U. Propagator
Butcher
0.6
|c0|
2
0.5
0.4
0.3
0.2
0.1
0.0
14
16
18
20
22
24
Time (fs)
26
28
30
8
t/ms
...
h(t)
h(t+t)


n 
n
h  t + t 1 =
h
(
t
)+
( h ( t +  t )- h ( t ))

ms
 ms  

 |c |
2
i
( n = 1 ..m s - 1)
t = 0.5 fs
i
3
ms = 1
2
1.04
1.03
ms = 10
1.02
1.01
ms = 20
1.00
0
10
20
30
40
50
60
70
80
90
100
Time (fs)
9
1.0
Occupation
Ad. pop.
0.8
0.6
0.4
0.2
0.0
0
20
40
60
80
100
Time (fs)
Uncorrected
H 2C
C
C
Butatriene cation
CH2
10
TDSE
E
Q
SC-TDSE
E
Because in the SC-TDSE the
“wave-packet” splitting among
the several electronic surfaces
is kept correlated by the
coordinate Rc, the time
propagation is fully coherent.
Q
• Schwartz, Bittner, Prezhdo, Rossky, J. Chem. Phys. 104, 5942 (1996)
• Zhu, Jasper, Truhlar, J. Chem. Phys. 120, 5543 (2004)
• Granucci, Persico, Zoccante, J. Chem. Phys. 133, 134111 (2010)
11
Decoherence is introduced ad hoc by correcting the time dependent
coefficients:
c k  c k exp    t /  ki 
'

2 
ci  ci   i  1 


 ki

'
2
ck
k i



1/ 2

C 

1 

Ek  Ei 
E kin 
C  0.1 hartree
• Granucci and Persico, J. Chem. Phys. 126, 134114 (2007)
12
1.0
1.0
Occupation
Ad. pop.
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0
20
40
60
80
100
Time (fs)
C
C
Butatriene cation
0.0
0
20
40
60
80
100
Time (fs)
Uncorrected
H 2C
Occupation
Ad. pop.
Corrected
CH2
13
    F12 v

H 12 
 exp   2 

 

2
P *   *
 1
H 11  H 22   t
0.57
Energy (au)
-224.65
*
cs
-224.70
*
-224.75
0.43
-224.80
cs
-224.85
0
2
4
6
8
10
12
14
Time (fs)
14
Pl  k 
P opulation increm ent in k due to flux fro m l during  t
P opulation of l
 lk  t   c l c k
*
Pl  k
 2t
 m ax  0,
 ll


1
Im   kl  H
• Tully, J Chem Phys 93, 1061 (1990)
c
lk
 R e   kl 

F  v 

c
kl
c
15
A hopping will take place if two conditions are satisfied:
1) A uniformly selected random number rt in the [0, 1] interval is such that
k 1
k
 P t   r   P t 
l n
n 1
l n
t
n 1
2) The energy gap between the final and initial states satisfies
Vk  R
c
 t    Vl  R c  t   
 N at c
c ,m 
  v m  F kl 
 m

N at
2 M
1
m
2
F 
c ,m
kc
2
m
16
E
Forbidden hop
Total energy
R
How to treat such situations:
 Reject all classically forbidden hop and keep the momentum
 Reject all classically forbidden hop and invert the momentum
 Use the time uncertainty principle to search for a point where the hop is allowed
• Jasper, Stechmann, Truhlar, J. Chem. Phys. 116, 5424 (2002)
17
E
Total energy
KN(t)
KN(t+t)
R
After hop, what are the new nuclear velocities?
 Redistribute the energy excess equally among all degrees
 Adjust velocities components in the direction of the nonadiabatic coupling h12
 Adjust velocities components in the direction of the difference gradient vector g12
• Pechukas, Phys. Rev. 181, 174 (1969)
• Fabiano, Keal, Thiel, Chem. Phys. 349, 334 (2008)
18
1.
Wigner Distribution for harmonic oscillator
PW  Q , P  
3N 6
  
1
i 1
2.
i2
  i  Hi O Q i 2

P
exp  
 i i
  R , P


HO


For each Rn in {R,P}
H eΦ k  r ; R n   V k Φ k  r ; R n   V k  R n  , f1k  R n 
3.
Compute cross section
 E 
4.
e
2
2 m c 0
N
fs

l2
 1

 N p
Np

n

f 1 l  R n  g  E   V1 l  R n  ,   

Compare to experiments
  10 ln 10 
3

NA
• Barbatti, Aquino, Lischka, PCCP 12, 4959 (2010)
• Barbatti, PCCP 13, 4686 (2011)
19
Wavelength (nm)
300
Ura
200
simulated
{R,P}
measured
-1
Cross section (Å .molecule )
0.4
250
2
0.3
0.2
*+nR3s
0.1
0.0
*
n*
*
4
5
6
7
8
Energy (eV)
20
surface hopping: panoramic view
21
1. Solve SE for Rc
H eΦ k  r ; R
c
  V Φ  r; R   V
c
k
k
k
,  V k , Fkj  Φ k  Φ j
2. Solve Newton’s equation on one surface
2
c
d Rm
dt
2

 mV l
c
0
Step 1 is the
computational
bottleneck
Mm
3. Integrate the SC-TDSE
i
dc k

dt
  V
c
k
 kj  i Fkj  v
c
c
c
j
0
j
4. Compute transition probability
Pl  k
 2t

*
c
c
 m ax  0,
R e c k c l F kl  v 
2
cl




• If it hops, then adjust momentum
• Repeat procedure until the end of
the trajectory
• Compute many trajectories
5. Decide surface for next time step
k 1
k
 P  t   r   P  t ,
l n
n 1
l n
t
rt  random  0,1
n 1
22
Simple implementation
• Propagation in Cartesian coordinates
• Trivial connection to different quantum chemical methods,
including QM/MM (regarding these methods can provide excited
state energies, energy gradients and nonadiabatic couplings)
• Independent trajectories: trivial parallelization
Local approximation
• No need of precomputing multidimensional potential energy
surfaces
• Straightforward on-the-fly implementation
• All nuclear degrees of freedom are propagated
23
Simple implementation
• High computational costs in the on-the fly approach
Local approximation
• Inconsistent treatment of zero point energy
• No treatment of tunneling effects
• Wrong coherence between states
24
Comparison to other methods
• Cattaneo and Persico, J. Phys. Chem. A 101, 3454 (1997)
• Worth, Hunt, Robb, J. Phys. Chem. A 107, 621 (2003)
Comparison between hopping algorithms
• Zhu, A. W. Jasper, and D. G. Truhlar, JCTC 1, 527 (2005)
• Fabiano, Groenhof, Thiel, Chem. Phys. 351, 111 (2008)
Conceptual background
•
•
•
•
Herman, J. Chem. Phys. 103, 8081 (1995)
Schwartz, Bittner, Prezhdo, Rossky, J. Chem. Phys. 104, 5942 (1996)
Tully, Faraday Discuss. 110, 407 (1998)
Schmidt, Parandekar, Tully, J. Chem. Phys. 129, 044104 (2008)
Surface hopping reviews
• Doltsinis, NIC series, 2002
• Barbatti, WIREs: Comp. Mol. Sci. 1, 620 (2011)
25
electrons
• classical
treatment
• quantum
treatment
nuclei
Atomic and molecular collision reactions
Molecular photochemistry and photophysics
Condensed phase dynamics
See References in:
• Barbatti, WIREs: Comp. Mol. Sci. 1, 620 (2011)
26
• classical
treatment
Slow degrees
of freedom
Fast degrees
of freedom
• quantum
treatment
• Tully, Faraday Discuss. 110, 407 (1998)
• Herman, J. Chem. Phys. 103, 8081 (1995)
Singlet-triplet transitions
• Carbogno, Behler, Reuter, Gross, Phys. Rev. B 81, 035410 (2010)
Electric field interactions
• Mitric, Petersen, Bonacic-Koutecky, Phys. Rev. A 79, 053416 (2009)
Solvent-induced vibrational relaxation
• Hammes-Schiffer, Tully, J. Chem. Phys. 101, 4657 (1994)
27
the newton-x program
28
Barbatti, Granucci, Ruckenbauer, Plasser, Crespo-Otero, Pittner, Persico, Lischka
NEWTON-X: A Package for Newtonian Dynamics Close to the Crossing Seam (2007-2013)
Interface
Method
Spectrum
Dynamics
Adiabatic
Nonadiabatic (TSH)
NAC
COLUMBUS
TURBOMOLE
GAUSSIAN
GAMESS
ANALYTICAL
DFTB
DFTB+
TINKER
DFT-MRCI
MRCI, MCSCF
TDDFT
CC2, ADC(2)
CASSCF
TDDFT
MCSCF
User defined
TD-DFTB
DFTB
MM
DFT-MRCI
CIO
+MM
+MM
+MM
+MM
+MM
• Barbatti, Granucci, Persico, Ruckenbauer, Vazdar, Eckert-Maksic, Lischka,
J Photochem Photobiol A 190, 228 (2007)
LD
Easy and practical of using:
• just make the inputs and start the simulations; monitor partial
results on-the-fly; get relevant summary of results at the end
Robust:
• if the input is right, the job will run: in case of error, messages must
guide the user to fix the problem
Flexible:
• some different case to study or new method to implement? It
should be easy to change the code
Open source:
• NX is the first MQCD-oriented program freely available and open to
the community
30
-----------------------------------------NEWTON-X
Newton dynamics close to the crossing seam
-----------------------------------------MAIN MENU
1. GENERATE INITIAL CONDITIONS
2. SET BASIC INPUT
3. SET GENERAL OPTIONS
4. SET NONADIABATIC DYNAMICS
5. GENERATE TRAJECTORIES
6. SET STATISTICAL ANALYSIS
7. EXIT
Select one option (1-7):
31
-----------------------------------------NEWTON-X
Newton dynamics close to the crossing seam
-----------------------------------------SET BASIC OPTIONS
nat: Number of atoms.
There is no value attributed to nat
Enter the value of nat : 6
Setting nat = 6
nstat: Number of states.
The current value of nstat is: 2
Enter the new value of nstat : 3
Setting nstat = 3
nstatdyn: Initial state (1 - ground state).
The current value of nstatdyn is: 2
Enter the new value of nstatdyn : 2
Setting nstatdyn = 2
prog:
Quantum chemistry program and method
0
- ANALYTICAL MODEL
1
- COLUMBUS
2.0 - TURBOMOLE RI-CC2
2.1 - TURBOMOLE TD-DFT
The current value of prog is: 1
Enter the new value of prog : 1
32
FSSH needs coefficients c:
Pl  k
 2t

*
 m ax  0,
R e c k c l Fkj  v 
2
cl


1) Through:

dc k

i
dt

V k c k    Fkj  v  c j
j
A) First order nonadiabatic couplings (NAC)
Explicit evaluation of Fkj
Only few methods (MCSCF, MRCI)
B) CI overlap (OVL)
Numerical evaluation of Fkj  v
F jk  v  t  
S jk  t  
1
4t
 3S  t   3S  t   S  t   t   S  t   t  
jk
kj
jk
kj
j t  t  k t 
Any method allowing to get a CI-like wavefunction (TDDFT, ADC(2), …)
• Hammes-Schiffer and Tully, J Chem Phys 101, 4657 (1994)
• Pittner, Lischka, MB, Chem Phys 356, 147 (2009)
34
FSSH needs coefficients c:
 2t

*
Pl  k  m ax  0,
R e c k c l Fkj  v 
2
cl


1


V t   TV t  t  T
1
1
t  c t 
2) Through: c  t   t   T exp   i ħ
2




C) Local diabatization (LD)
Löwdin orthogonalization of CI overlap matrix S
T  SO 
 1/ 2
O
t
S SO  O 
t
Any method allowing to get a CI-like wavefunction (TDDFT, ADC(2), …)
• Granucci, Toniolo, Persico, J Chem Phys 114, 10608 (2001)
• Plasser, Granucci, Pittner, MB, Persico, Lischka, J Chem Phys 137, 22A314 (2012)
35
Percent relative error
30
NAC
OVL
LD
20
10
Landau-Zener model system:
Relative error for the asymptotic
diabatic transition probability x t
0
-10
LD is especially adequate to weak
nonadiabatic coupling regimes
-20
-30
0.01
0.1
1
Timestep length (fs)
• Plasser, Granucci, Pittner, MB, Persico, Lischka, J Chem Phys 137, 22A314 (2012)
36
The nuclear Ensemble Method
OH
4
6
5
HN
3
2
N
4
9
67 7
hn
8
1
trans
6
5
HN
7
3
O
10
2
N
8
1
cis
HO
9
• Review: Gibbs, Tye, Norval, Photochem Photobiol Sci 7, 655 (2008)
O
10
OH
4
6
5
HN
2
N
hn
8
67 7
3
4
9
6
O
2
10
cis
8
HO
9
O
10
UVA
Ftrans→cis (289 nm) = 0.08
2
-1
Cross section (Å .molecule )
UVB
N
1
trans
UVC
7
3
1
1.0
5
HN
0.5
Ftrans→cis (302 nm) = 0.31
Ftrans→cis (313 nm) = 0.49
0.0
260
280
300
320
Wavelength (nm)
340
N
NH
OH
HN
OH
N
O
O
2
Cross section (Å )
1.0
Exp.
0.5
Theor.
0.0
280
320
Wavelength (nm)
• Barbatti, PCCP 13, 4686 (2011)
First order of the time-dependent perturbation theory
 E 

3 c 0 n r E

 E 00 , nk  R   00  R   0   r , R   n
 nk  R     E 00 , nk  R   E  d R
2
*
r
n ,k
 E 00 , nk  E nk  E 00   E 0 , n
The problem can be recast in the time domain:
 E 
1
6
2
c 0 n r

E
n
R e   E 0 , n  R  M 0 n  R     00  R   n  R , t  e

M 0n  R   0  e  r; R  n
 n t   e
 iH n t /
2
r
2
*
 E 0 , n   E 00 , nk
i  E  E 00  t /
dt  d R

k
 00
• Sakurai (1994) Modern Quantum Mechanics
• Tannor, Heller, J Chem Phys 77, 202 (1982)
41
Relation between cross section (cm2) and extinction coefficient (M-1cm-1)
  10 ln 10 
3

NA
42
The core of the method is to compute the overlap
  R , t    00  R   n  R , t 
*
needed to integrate
 E 
1
6
2
c 0 n r

Re  E  R  M  R      R    R ,t  e


E
2
0 ,n
2
0n
*
00
n
i E  E 00  t /
n
dt  d R

Tannor and Heller proposed a analytical solution based on harmonic oscillator
 00  n  t 


   2nj
i j t  i
 i j t
a
 exp   
1 e

   E 0 , n t 
 2
2 
 j 



43
But we want to go beyond the Condon approximation.
Starting from Tannor-Heller equation:
 00  n  t 


   2nj
i j t  i
 i j t
a
 exp   
1 e

   E 0 , n t 
 2
2
 j 




The expansion to second order is


 2
1
2
a
2
2


1



t





  nj j  t 

nj
j
0n 
2 j
2 j

 
1
 00  n  t   exp   i 

This motivates to introduce the following functional:

  n  R , t    00  R  exp  
*
00
2

i
 E 0 ,n  R  t 
i
 nt 
i
E 00 t 
1

nt 
2
8

• Crespo-Otero, Barbatti, Theor Chem Acc 131, 1237 (2012)
2 2
44

 00  n  R , t    00  R  exp  
2
*
i
 E 0 ,n  R  t 

 E 
 E 
1
6
2
c 0 n r
 e
i
 nt 
i
E 00 t 
1
2 2

Re  E  R  M  R      R    R ,t  e


E
2
2 m c 0 n r E
2
0 ,n
2
0n
*
00
n


nt 
2
8

i E  E 00  t /
n
dt  d R

 00  R   E 0 , n  R  f 0 n  R  g G auss  E   E 0 , n  R    n  R  ,  n  d R
2
n
g G auss  E   E 0 , n  R    n  R  ,  n  
1
 2  
n
/ 2
2

1/ 2
   E  E  R    2
0 ,n
n
exp 
2

2  n / 2 





45
 E 
 e
2
2 m c 0 n r E

 00  R   E 0 , n  R  f 0 n  R  g  E   E 0 , n  R    n  R  ,  n  d R
2
n
If we have a ground state distribution of points:  00  R l 
2
We can integrate the cross section by Monte-Carlo and get:
 E 
e
2
2 m c 0 n r E
N
fs

n
1
Np
Np
 E  R  f  R  g  E  E  R  , 
0 ,n
l
0n
l
l
n
l
46
Using a Wigner distribution to a harmonic oscillator, we can sample the ground
state:
 00  q  
2
3N 6

j 1
  j j 


 

1/ 2
exp    j  j q j /
2

Another way of sampling is to run a very long trajectory in the ground state
and pick points from it.
But be careful with the right temperature!
47
Pros:
• Easy to use
• Clear conceptual basis
• Absolute heights
• Absolute widths
• Post Condon Approximation
• Dark vibronic bands
• Implemented in Newton-X
Cons:
• No vibrational resolution
• No non-adiabatic info
• No info on excited state wavefunction
• One arbitrary parameter
48
0.3
3
0.010
0.2
2
0.005
0.1
1
Simulated
Expt. Bolovinos et al. 1982
2
-1
Cross section (Å .molecule )
0.015
0.000
4.5
5.0
Energy (eV)
5.5
0.0
5.6
5.8
6.0
6.2
Energy (eV)
6.4
0
4.5
5.0
5.5
6.0
6.5
7.0
Energy (eV)
TD-CAM-B3LYP/TZVP
49
2
Cross section (Å )
NE
LVC
Expt.
0.002
H 3C
N
N
CH3
0.001
azomethane
0.000
2.5
3.0
3.5
4.0
4.5
5.0
Energy (eV)
• Szalay, Aquino, Barbatti, Lischka, Chem Phys 380, 9 (2011)
tra n s -2
tra n s -1
7
2
1.0
5
-1
-1
 (10 M cm )
0.8
e
trans-1
trans-2
cis-A1
Total
Expt.
d
e'
0.6
d'
N
H
21
N
24
22
23
NH2
N
N
9
H
N
H
N
NH2
N
H
N
12
17
c is -A 2
c is -A 1
b
0.4
c
b'
a
c'
0.2
N
H
N
H
N
N
NH2
N
N
H
N
H
N
NH2
a'
0.0
300
400
500
600
700
800
900
c is -B 2
c is -B 1
Wavelength (nm)
N
N
• Lan, Nonell, Barbatti, J Phys Chem A 116, 3366 (2012)
N
H
NH2
N
H
H
N
H
N
Porphycene
N
N
NH2
• Homem, Lopez-Castillo, Barbatti, Rosa, Iza, Cavasso-Filho, Farenzena, Lee, Iga,
J Chem Phys 137, 184305 (2012)
52
Challenges in computational
simulations of nucleobases dynamics
E2
1H
6
MRCIS-2n
6H
1
2E
6
E2
1H
6
OM2/MRCI
N1
6E
2E
2
E2
C2
H
E6
Ab initio and semiempirical
MRCI produce
divergent results
54
TD-PBE0
TD-BHLYP
E2
6
2
E
H1
6
E
6
6
H1
2
E
2E
E
E2
E2
1H
6
E6
55
Ground state population % at 1 ps
C2 puckering
NH2
C6 puckering
H elimination
C
Experimental
S0 population at 1 ps (%)
85 8
80
68 2
60
57 12
59 8
N
Right for
different reasons
C
CH
HC
C
N
40
5 8
10 16
0
wrong
Ex
pe
r
im
en
ta
AD l
C
(2
M )
R
O
C
I
M
2/ S
M
R
TD CI
TD -P
- BE
B9
TD 7X
-B D
3L
TD TD YP
-C -P
AM BE
-B 0
3
TD LY
-B P
TD HLY
-M P
06
-H
F
0
N
H
25 16
20 21
20 15
12 8
20
N
• MB, Lan, Crespo-Otero, Szymczak, Lischka, Thiel, J Chem Phys 137, 22A503
(2012)
56
Pump-probe signal depends
on changes in the ionization
cross section and in the
ionization energy
Variations in the ionization
energy are not usually taken
into account neither in
simulations nor in
experimental analysis
• Barbatti and Ullrich, PCCP 13, 15492 (2011)
57
58
It is not only an excited-state
problem. Ground state is also bad.
59
Ground state simulations: 1- 2 kcal/mol accuracy
Excited state simulations: 5-10 kcal/mol accuracy
Reason:
• Excited states spectral region has high density of states
• Small variation in geometry leads to a change of electronic
character
• No affordable method can describe all characters at the same
level
60
• Merchan and Serrano-Andres, JACS 125, 8108 (2003)
61
Next lecture:
• Case studies of dynamics simulations
Contact
[email protected]