Mario Barbatti Institute for Theoretical Chemistry University of Vienna

1. Photochemistry: adiabatic and nonadiabatic molecular dynamics with  multireference ab initio methods  Mario Barbatti Institute for Theoretical Chemistry University of Vienna COLUMBUS in BANGKOK (3-TS2C2) Apr. 2 - 5, 2006 Burapha University, Bang Saen, Thailand

2. Outline • First Lecture: An introduction to molecular dynamics • Dynamics, why? • Overview of the available approaches • Second Lecture: Towards an implementation of surface hopping dynamics • The NEWTON-X program • Practical aspects to be adressed • Third Lecture: Some applications: theory and experiment • On the ambiguity of the experimental raw data • On how the initial surface can make difference • Intersection? Which of them? • Readressing the DNA/RNA bases problem

3. Part IITowards an implementation of surface hopping dynamics Cândido Portinari, O lavrador de café, 1934

4. Quick overview of the method

5. where Time derivative Nonadiabatic coupling vector Tully, JCP 93, 1061 (1990); Ferretti et al. JCP 104, 5517 (1996) Population: • Two electronic states are coupled only by the nonadiabatic coupling vector hij (adiabatic representation). Surface hopping approach: adiabatic population

6. Hammes-Schiffer and Tully, JCP 101, 4657 (1994) Surface hopping approach: fewest switches Tully, JCP 93, 1061 (1990)

7. E • The transition probability gkj between two electronic states is calculated at each time step of the classical trajectory. • A random event decides whether the system hops to other adiabatic state. t Surface hopping approach: semiclassical trajectories • Nuclear motion is obtained by integrating the Newton eqs. • At each time, the dynamics is performed on one unique adiabatic state, Ei = Hii. • In the adiabatic representation, Ei(R), Ei, and hji are obtained with traditional quantum chemistry methods. • aji is obtained by integrating

8. Conventional Preparation of surface dynamics On-the-fly dynamics Total time (Dis)advantages of the on-the-fly approach • Advantages: • It is not need to get the complete surface. Only that regions spanned during the dynamics • It dispenses interpolation, extrapolation and fitting schemes • Disadvantages: • Time-expensive dynamics • No non-local effects (tunneling)

9. Practical aspects to be addressed

10. E Sn S0 R,P How to prepare initial conditions • Microcanonical distributions: • Classical harmonic distribution • R-P uncorrelated quantum harmonic distrib. (Wigner) • R-P correlated quantum harmonic ditrib. • Canonical distribution: • Boltzmann

11. Classical dynamics: integrator Any standard method can be used in the integration of the Newton equations. A good one is the Velocity Verlet (Swope et al. JCP 76, 637 (1982)): For each nucleus I Quantum chemistry calculation

12. Time-step for the classical equations Schlick, Barth and Mandziuk, Annu. Rev. Biophys. Struct. 26, 181 (1997).

13. Time-step for the classical equations • Time step should not be larger than 1 fs (1/10v). • Dt = 0.5 fs assures a good level of conservation of energy most of time. • Exceptions: • Dynamics close to the conical intersection may require 0.25 fs • Dissociation processes may require even smaller time steps

14. TDSE: integrator Numerical Recipes in Fortran Constant time step: • Fourth-order Runge-Kutta (RK4) Adaptive time step: • Bulirsh-Stoer Adaptive works better than constant time step

15. TDSE: integrators • Some step-constant integrators available in NEWTON-X: • Polynomial, 3rd order • Runge-Kutta, 4th order • Adams Moulton predictor-corrector, 5th order • Adams Moulton predictor-corrector, 6th order • Unitary propagator • Butcher, 5th order

16. Dt = 0.5 fs Time-step for the quantum equations Dt/ms . . . h(t) h(t+Dt)

17. 1 Probability of hopping P2→1 0 Fewest switches: two states Population in S2: Trajectories in S2: Minimum number of hoppings that keeps the correct number of trajectories:

18. 1 Example: Three states P3→2 +P3→1 Only the fraction of derivative connected to the particular transition P3→2 0 Fewest switches: several states Tully, JCP 93, 1061 (1990) Hammer-Schiffer and Tully, JCP 101, 4657 (1994)

19. Forbidden hop Total energy E R Forbidden hops • Forbidden hop makes the classical statistical distributions deviate from the quantum populations. • How to treat them: • Reject all classically forbidden hop and keep the momentum. • Reject all classically forbidden hop and invert the momentum. • Use the time uncertainty to search for a point in which the hop is allowed (Jasper et al. 116 5424 (2002)).

20. E R Adjustment of momentum after hopping Total energy KN(t) KN(t+Dt) 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 vector h12 • Adjust velocities components in the direction of the difference gradient vector g12 • Adjust velocities in the direction

21. Phase control CNH4+: MRCI/CAS(4,3)/6-31G* Compare h(t) and h(t+ Dt)

22. Phase control CNH4+: MRCI/CAS(4,3)/6-31G* Compare h(t) and h(t+ Dt)

23. 12.00 fs Abrupt changes control 11.75 fs

24. The routine also gives the linear parameters: Orthogonalization g-h space orthogonalization [Yarkony, JCP112, 2111 (2000)]

25. When surface hopping fails • SH is supposed to reproduce quantum distributions, in the sense that • fraction of trajectories (t) = adiabatic population (t) Eq. (1) • This statement should be true for: • Number of forbidden hops →0 • Number of trajectories → infinity • Granucci and Persico have shown that for some cases, even if these conditions are satisfied, Eq. (1) may be not true. • SH, as any trajectory-independent semiclassical method, cannot account for quantum interference effects and quantization of vibrational and rotational motions. • It is unclear how good the fewest switches approach in the proximity of conical intersections is.

26. NEWTON-X: a package for Newtonian dynamics close to the crossing seamM. Barbatti, G. Granucci, H. Lischka and M. Ruckenbauer (2005-2006)

27. NX aims • 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: in the future, NX should be opened to the community.

28. NX input facility: nxinp ------------------------------------------ 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):

29. NX input facility: nxinp ------------------------------------------ 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

30. R(t), v(t) t+t, R(t+t), v(t+t/2) Tables and graphics provide: Ek(t+t), hkl(t+t) akk, Pkl(t+t) v(t+t) Statistical analysis NX modular design • Fortran 90 routines • Perl controller Initial condition generation

31. Adiabatic dynamics R(t), v(t) t+t, R(t+t), v(t+t/2) provide: Ek(t+t) v(t+t)

32. R(t), v(t) t+t, R(t+t), v(t+t/2) provide: Ek(t+t), hkl(t+t) akk, Pkl(t+t) v(t+t) Methods available • Presently: • COLUMBUS [(non)adiabatic dynamics] • MCSCF • MRCI • TURBOMOLE [adiabadic dynamics] • TD-DFT • RI-CC2 • Analytical models [user provided] • Being implemented: • COLUMBUS + TINKER • QM/MM [(non)adiabatic dynamics] • To be implemented: • ACES II • EOM-CC [(non)adiabatic dynamics]

33. So many choices… What method should I use?

34. Present situation of quantum chemistry methods

35. RI-CC2 TD-DFT CASSCF Comparison among methods CNH4+: MRCI/CAS(4,3)/6-31G* Adiabatic dynamics can be used to find out the most relevant relaxation paths. But be careful with the limitation of each method (CT states in TD-DFT for example).

36. A basic protocol • Use TD-DFT for large systems (> 10 heavy atoms) with one single configuration dominating the region of the phase space spanned by the dynamics. Test against CASSCF and RI-CC2. • Use RI-CC2 for medium systems (6-10 heavy atoms) under the same conditions as in the previous point. • Use CASSCF for medium systems with strong multireference character in all phase space. Test against MRCI. • Use MRCI for small systems (< 6 heavy atoms). • In all cases, when the number of relevant internal coordinates is small (2-4) and they can easily be determined, test against wave-packet dynamics.

37. Conclusions Correlation Method Full-CI … MRCI CASSCF HF Basis set DZ TZ … BS limit

38. < 6 heavy atoms 6-10 heavy atoms Static calculations Adiabatic dynamics Surf. Hoppingand Mean Field Multiple Spawning Wavepacket dynamics Dynamics method Conclusions Correlation Method Full-CI … MRCI CASSCF HF Basis set DZ TZ … BS limit

39. This lecture: • Surface hopping is one of the most popular methods available for nonadiabatic dynamics • Its implementation is direct and it can be used with any quantum chemistry method that can provide analytical excited-state gradients and analytical nonadiabatic coupling vectors • These requirements are fulfilled by only a few methods such as CASSCF, MRCI and (partially) EOM-CC • Next lecture: • Adiabatic and nonadiabatic dynamics methods will be used in the investigation of some examples of photoexcited systems