Faraday Discuss. Chem. SOC., 1991, 91, 239-247 Non-adiabatic Reactive Routes and the Applicability of Multiconfiguration Time-dependent Self-consistent Field Approximations Ronnie Kosloff" and Audrey Dell Hammerich Department of Physical Chemistry and The Fritz Haber Research Center for Molecular Dynamics, The Hebrew University of Jerusalem, Jerusalem 91 904, Israel A time-dependent non-adiabatic formulation is considered for chemical reactions. Basic algorithms are presented for the photodissociation of CH31 under the influence of strong short laser pulses. Detailed insight into the dynamics of this process is demonstrated. A time-dependent self-consistent field approach is suggested when many degrees of freedom have to be considered simultaneously. The derivation is based on a Liouville space description for which quantum and classical mechanics are treated on equal grounds.A multiconfiguration approach is formulated for explicitly includ- ing correlation due to splitting of the probability density and for non- adiabatic motion. 1. Introduction It is apparent that the naive adiabatic picture of a chemical reaction proceeding on only one Born-Oppenheimer potential-energy surface is an oversimplification. The nature of a chemical reaction, in which large changes in the electronic configuration occur within small distances, leads to a non-adiabatic formulation. Photochemical reactions commonly involve more than one potential-energy surface. A typical example is the photoisomerization of rhodopsin, for which at least three potential surfaces can be identified.' A similar phenomenon arises in the interaction of molecules with metal surfaces.For the dissociation of nitrogen on catalytic surfaces, an important step in ammonia synthesis, at least two stable chemical species have been found to participate in the reactiom2 Also for the dissociation of oxygen on silver the different electronic species 02, 0, and 0- have been ide~~tified.~ The adiabatic picture fails completely for hyperthermal surface ionization processes in which a continuum of electronic states has to be considered, owing to the proximity of the ionization process to the s ~ r f a c e . ~ In this paper a strategy is proposed based on a time-dependent quantum-dynamical framework in which the molecular dynamics are studied. The traditional approach to this problem starts with the work of Landau and Zener.' Their approach preserves the Born-Oppenheimer potential-energy picture and corrects for isolated breakdown of the approximation.This is basically a one-dimensional scheme, suitable for isolated curve- crossing events. An extensive review of this approach and related developments can be found in the work of Nikitim6 The one-dimensional model for non-adiabatic processes, however, is not adequate for multidimensional cases. In particular, the LZ approximation fails when the motion is parallel to the crossing seam between the electronic surfaces. For large-scale realistic simulations, practical procedures have to be developed to overcome the above difficulty. The semiclassical method of Tully and Preston' is such a procedure.Basically it is a classical trajectory method which includes surface hopping calculated by the LZ formula or by some other improved non-adiabatic method. Its main disadvantage is that phase coherence is lost on each crossing event. Full quantum- mechanical calculations of multidimensional non-adiabatic scattering events have been 239240 Non-adiabatic Reactive Routes performed, as in the work of Baer on the HT system and the H l + H2 reaction.' The problem becomes far more complicated when the dimensionality increases further, as is found in full 3D calculations or for chemical reactions in a solvent or on a surface for which infinitely many degrees of freedom participate. In the next sections a procedure based on a time-dependent quantum-mechanical approach is presented.Section 2 presents the basic procedures of a non-adiabatic process by the example of the photodissociation of CH31. Section 3 summarizes the time- dependent self-consistent field (TDSCF) approximation which is the main tool used to reduce the dimensionality of the problem. Section 4 considers correlated motion by a multiconfiguration approach in Liouville space. 2. Time-dependent Quantum-mechanical Approach It is best to review the basic time-dependent methods through a specific example. The photodissociation of CHJ will serve this purpose. It is a system well studied both experimentally and theoretically.' In this process two electronically different iodine atoms are produced. This means that the photodissociation proceeds on at least three electronic surfaces.The time-dependent Schrodinger equation of the system has the form: where t,bulv is the projection of the w$vefunction onto the upper surfaces. +, isjhe Rrojeztion onto the ground surface. Hi is the upper/lower surf%ce Hamiltonian: Hi = T, + V, where T, = P2/2p the kiTetic energy operator and V, is the upper/lower potential-energy surface. V,, and Vvg are the interaction potentials: Vig = -pi,E( t ) , where pig is the dipole operator to each syrface, and E ( t ) represents the amplitude of a semiclassical electromagnetic field. The Vu, term represents the non-adiabatic coupling between the two upper surfaces. It is customary to treat the system within a limited dimensionality configuration space, including the dissociative I -CH3 coordinate and the umbrella CH3 vibrational motion.The time-dependent approach solves the dynamics, including the radiative coupling, starting at t = 0 with the system in the ground state of the lower electronic surface." The calculation is discretized on a spatial grid with evenly distributed sampling points. Each electronic surface has its own grid. It should be noted that the electromagnetic field is included directly in a time-dependent Hamiltonian operator allowing the simula- tion of strong field effects and pulse shaping." The kinetic energy term in the Hamiltonian is calculated by transforming to momentum space by the fast Fourier transform (FFT) algorithm, multiplying by ( k 2 / 2 ) p and then back-transforming by an inverse FFT to configuration space.Using this method the Hamiltonian is represented with very high accuracy owing to exponential convergence with respect to grid spacing. The computa- tional effort scales as N log N, where N is the total number of grid points needed to represent the problem." The initial ground-state wavefunction is foynd by propagating the Schrodinger equation in imaginary time with the Hamiltonian H, until the ground state is ~onverged.'~ The dynamics of the dissociation are obtained by solving eqn. (2.1) subject to the initial conditions. Because the Hamiltonian is time dependent, the total propagation time is divided into small segments Af in which the variation of the Hamiltonian is small. A short-time propagator is used to advance + one interval of time.+ ( t + A t ) =eo+(t) (2.2)R. Koslofl and A. D. Hammerich where the operator d is obtained from the Magnus series: 241 t+At s H ( s ) d s - l I t 2h2 t [ I ; l ( s ) , I j ( s ’ ) ] d s d s ’ + . . . The propagator e0 in eqn. ( 2 . 2 ) is expanded as a polynomial allowing a recursive algorithm. It has been found14 that an optimal way to represent this polynomial is by the Newton interpolation formula: e 0 = a,i + a, ( O - z o i ) + a,( O - z& d - z , i > + - * where zi are sampling points located within the energy spectrum of the operator 6 and ai are the divided difference coefficients. For example a, =f( z,) = ezo, a , = (ezl - e’o)/(z, - zo). The sampling points zi are optimized so that there will be a mjnimum residuum from the approximation with respect to the operator function f(0).The familiar Chebyshev expan~ion’~ can be obtained by choosing the sampling points z, as zzros of the Chebyshev polynomial of degree M + 1 which spans the energy range of 0. This time propagation method has been shown to have accuracy matching the convergence characteristics of the Hamiltonian representation. l4 A typical encounter of the photodissociation of CHJ induced by a short pulse of 10 fs is displayed in Fig. 1. The field intensity is 1 x 1014 W cm-*, leading to 85% excitation. The three panels exhibit the wavefunction density on the ground, radiatively coupled and non-adiabatically coupled surfaces with irradiation by a 10 fs resonant Gaussian pulse. The x axis is the dissociative coordinate and the y axis the CH3 umbrella vibration.Stimulated emission is evidenced by the loss of symmetry and structure apparent in the ground-state wavefunction. The dynamics are highly non-linear as the 17 7~ pulse employed cycles density many times between the two strongly coupled electronic surfaces. The non-adiabatic coupling is also strong resulting in a branching ratio of 3/1 between the two exit channels. The product in the I channel resulting from the curve crossing is vibrationally hotter than that in the I* channel, attributing to the greater width along the y axis (umbrella motion) in the uppermost panel of the figure. A more detailed examination of the dynamics reveals that the non-adiabatic coupling cycles density between the coupled surfaces even far into the dissociative channel.The richness of the photodissociation dynamics can be better illuminated by compar- ing the diabatic potential-energy surfaces with their ‘adiabatic’ counterparts formed by diagonalizing the potential matrix which includes the off-diagonal radiative and non-adiabatic coupling terms. As these potentials include the field they are time depen- dent. Fig. 2 shows the minimum-energy path (reaction coordinate) along the diabatic curves and their ‘adiabatic’ analogues, the latter at 30fs, the peak field intensity. The upper field dressed potential exhibits a shallow well capable of supporting a metastable state that can undergo several cycles of absorption and emission. The lowering of the ground-state potential and raising of the excited-state potential (for instance around 8 a.u.) implies that the energy gap between the two surfaces is closer to the laser central frequency over a long distance, offering an efficient mechanism for the transfer of amplitude between the two surfaces.Recent calculations by Morokuma and co-workers’6 have shown that the validity of the limited dimensionality calculation is questionable due to breaking of the C,, sym- metry. As a result at least a four-dimensional calculation has to be performed in order to include the symmetry-breaking degree of freedom. Although a 4D calculation is feasible on modern computers, a 9D calculation which includes the C-H stretch coordinates or a simulation of the photochemical reaction in a condensed medium is beyond the performance of current computers and therefore requires a new strategy.The next section will outline the TDSCF approach.242 Non-adiabatic Reactive Routes Fig. 1 Simulation of the photodissociation of CH31 induced by a resonant 1Ofs Gaussian laser pulse. From bottom to top are displayed the square modulus of the wavefunction at 44 fs on the ground, radiatively coupled and non-adiabatically coupled electronic potential-energy surfaces. The x axis represents the dissociative I-CH3 coordinate and the y axis the CH3 umbrella vibrational motion 3. TDSCF Formulation The goal of this section is to outline a consistent procedure for modelling a quantum system in contact with a larger system which can serve as a heat bath. The physical system under investigation is very demanding. The main problem relates to its multiple bodied aspect so that a full quantum-mechanical description of all degrees of freedom is beyond current computational capabilities.A practical modelling of the processes, therefore, requires a hierarchical procedure in which only the relevant degrees of freedom are treated in a full quantum-mechanical fashion. Other degrees of freedom are described by a classical or semiclassical approximation. Mixed quantum/classical calculations have been tried by many authors.” Owing to the difference in formulation of the quantum and classical mechanics, mixing the theories often leads to conflicts with questionable interpretation. A consistent treatment has to be formulated on a common denominator. A density description in Liouville space can serve such a purpose.This approach overcomes the difficulty that a single trajectory has no direct analogue in quantum me ch an i cs.R. Koslofland A. D. Hammerich Oa20 n 243 0.03 - I_- Fig. 2 Reaction coordinates on the diabatic CH31 potential-energy surfaces and their ‘adiabatic’ analogues which include the radiative and non-adiabatic couplings. The upper panel portrays the diabatic surfaces’ with the two horizontal lines denoting the energy of the ground state before and after full resonant excitation. In the lower panel are the corresponding ‘adiabatic’ curves at 30 fs. Since the electric field is time dependent the ‘adiabatic’ potentials are also The correct many-bodied solution is formally solved by the Liouville-von Neumann equation for the density operator: dP&& at where the Liouville operator can be separated into individual contributions: 9 = zs + 9 b + z s b (3.2) where the subscript s indicates the primary system and b the bath.A numerical scheme to solve the Liouville-von Neumann equation under dissipative conditions which result from contact with a heat bath has been developed.** The method is based on a direct expansion of the density operator on a grid and a propagation scheme able to incorporate complex eigenvalues of the evolution operator for dissipative systems. Despite the appealing aspect of a direct solution to eqn. (3.1), the procedure is numerically expensive and therefore limited to model systems. For realistic simulations, further approximations are required, the main one being the use of a mean-field approach which expresses the full density operator as a tensor product of the density operators of the primary system s and of the bath degree of freedom b: p^‘p^s@p^b (3.3)244 Non-adiabatic Reactive Routes Such a separation replaces the full multiple bodied correlated motion with a mean-field interaction with the equations of motion: (3.4a) (3.4b) Once a partition of the density has been made, each of the two parts can be treated in a different fashion.For example, the nuclear motion of the bath can be described classically or semiclassically. A further approximation is to specify that the initial state be a pure state (zero temperature). In the mean-field approach of eqn. (3.4) this reduces eqn. (3.4~2) to a Schrodinger-type equation: (3.5a) where and is is the primary system Hamiltonian, psb(&) is the system-bath interaction potential, and R b is the bath nuclear coordinate.The effective potentials in eqn. ( 3 . 5 ~ ) assume a semiclassical approximation in which the averaged potential over the nuclear wavefunction can be replaced by the potential of the mean position. Similarly eqn. (3.4b) is replaced by Hamilton’s equations for the mean position and momentum: a: - (a::)s (3.5b) 6 s = l ~ s ) ( ~ s l -- - dPRb= -(..EJs at where H is the averaged classical Hamiltonian function and (-)s represents an average over the primary system’s degrees of freedom. A generalized Langevin (GLE) approach’’ can be used to describe the influence of the rest of the heat bath. Explicitly, the primary part exerts a force on the ith bath nuclei which is calculated by the matrix element: Mixed quantum/ classical calculations of this type have been applied successfully to many but in view of the non-adiabatic reactions considered in this work this approximation is bound to fail due to the omission of important correlations in eqn.(3.2). First, the bath subsystem has to know on which electronic surface the primary system resides in order to feel the correct force. Secondly, in reactive scattering, an important correlation is built up which is the result of the density splitting into the reactive and non-reactive parts.2’ These important correlations imply a modification of the simple mean-field approach. The method described below, a multiconfiguration approach, addresses this problem.4. MC-TDSCF Formulation For simplicity the formulation is developed for three configurations. Extension to more configurations follows the same lines. Define a projection operator P in Liouville space: P2 = P. (4.1)R. Koslof and A. L3. Hammerich 245 The projection operator P is defined in the s space and separates the system into two parts 1 and 2. These two parts for example could be the two electronic surfaces: PP: = P: Pp,2 = 0. Using the projection operator P in the s space the projection Q is defined for the other surface: P + Q + W = I (4.3) where the projection operator W selects the parts that are not included in the projections P and Q. These parts of the density operator are off-diagonal and represent correlation between the 1 and 2 parts.The tensor product density operator (3.3) is replaced by a correlated density operator: The bath density operator Pb is chosen to be normalized: tr {p:} = tr {p:} = tr {p;} = I b b b The total normalization leads to: tr {p: + p:} = 1 tr {p:} = o S S (4.5) because, by construction, the projection W selects off-diagonal terms of the density operator. The time derivative of eqn. (4.4) becomes: (4.7) now: PQ = 0 then: p = p,'Op:+p:Op;+ p,3op;+p;Op:+p,2op:+ p:op: Pdp=P2(P+Q+ W ) p at dp- QX(P+ Q+ W ) p Qdt- (4.8) Wdp= WX(P+Q+ W)p Combining eqn. (4.8) with the definition of the projection operators, and taking partial traces the following equation is obtained: at (4-9) Similar equations are derived for the other parts.246 Non-adiabatic Reactive Routes (4.10) which represents couple$ motio$ on two electronic surf5ces.Combining eqn. (4.9) with (4.10) and using PX = I1)( 11x1 1)( 1 I for any operator X one obtains: (4.1 1 ) where the [ , ]+ symbol represents the anticommutator, v: 1( r,) = tr{ Vll( r, , Rb)pi} and v;2( r,) = tr{ V12( r,, Rb)p:} are potentials averaged over the bath density operators p: and p:. The equation of motion for the bath mode becomes: (4.12) where v:l(Rb) =tr{Vll(rs, Rb)p:} and v;2(Rb) =tr{V12(rs, Rb)p?} A similar equation is obtained for p2. It should be noted that there is no direct coupling between the equations of motion of pi and pf. Owing to coupling to the bath and the fact that each surface bath moves under different conditions, the equations of motion are able to create statistical mixing in the individual components p:.Examining the structure of eqn. (4.12) it can be found that a semiclassical approximation for the bath will contain three coupled bath coordinates. Each of these coordinates will move under a different average force. 5. Discussion Time-dependent quantum-mechanical methods for molecular dynamics have evolved considerably. Calculations of realistic molecular encounters are able to shed light on fundamental processes. The purpose of this work is to gain insight into non-adiabatic molecular processes. The photodissociation of CHJ was brought as a demonstration of how one can elucidate a molecular process induced by a short and strong laser pulse. This type of intimate understanding, in particular the dynamical picture displayed by the instantaneous dipole, can lead to control of the chemical outcome.The coherent evolution is a key factor in achieving such control. Considering reactions with many degrees of freedom such as in condensed phases the main issue is whether the phase coherence is maintained long enough in order to achieve chemical control of the process. The surface-hopping approximation of Preston and Tully represents the other extreme view in which phase coherence is only maintained for a very short period. The density formulation presented in sections 3 and 4 has been constructed with this issue in mind. Simulation of ionization processes in strong fields22 has shown that phase coherence is preserved for very long times. It was found that extremely accurate integration schemes are therefore required in order to follow long-time interference effects.It is therefore expected that coherent chemical manipulations are possible in condensed phases. The time-dependent quantum-mechanical formulation, in particular eqn. (4.9)- (4.12) is a rich source of physical approximations. The possibilities range from a full quantum-mechanical calculation to a mixed quantum and semiclassical approach.R. Koslofl and A. D. Hammerich 247 Multiconfi uration time-dependent self-consistent field approximations have been tried before.*l* ' 3-26 . Their derivation was based on a wavefunction representation. The work of Kotler et u Z . * ~ directly addresses multisurface dynamics. For such cases, the choice of the projection operator is obvious. Other work2' emphasizes the correlation produced in reactive scattering.The work of Meyer et al.26 is a generalization to many configurations. The density operator formulation described here differs in major aspects from previous work. In quantum mechanics it is possible that a projection of a combined pure state on a local one-particle description of state leads to a mixed state. The influence of the bath degrees of freedom in eqn. (4.12) will lead to this effect. This research was supported by a grant from the G.I.F., the German-Israeli Foundation for Scientific Research and Development. The Fritz Haber Research Center for Molecular Dynamics is supported by the Minerva Gesellschaft fur die Forschung, GmbH Munchen, Germany. References 1 I. Ohmine, personal communication. 2 G.Haase, M. Asscher and R. Kosloff, J. Chem. Phys., 1989, 90, 3346. 3 M. E. M. Spruit, PhD Thesis, University of Amsterdam, 1989; M. Rocca and U. Valbusa, Phys. Rev. 4 A. Dannon, E. Kolodney and A. Amirav, Surf: Sci., 1988, 193, 132; A. Amirav and M. J. Cardillo, 5 L. D. Landau, Phys. 2. Sow., 1932, 2, 46; C. Zener, Roc. R. SOC. London, Ser. A, 1932, 137, 696. 6 E. E. Nikitin, Chemische Elementarprozesse, ed. Hermann Hartman, Springer-Verlag, Berlin, 1968. 7 R. K. Preston and J. C. Tully, J. Chem. Phys., 1971, 54, 4297. 8 Z. Top and M. Baer, Chem. Phys., 1977, 25, 1; M. Baer and C. Y. Ng, J. Chem. Phys., 1990,93, 7787. 9 H. Guo and G. C. Schatz, J. Chem. Phys., 1990, 93, 393; M. Shapiro, J. Phys. Chem., 1986, 90, 3644. Lett., 1990, 63, 2398. Phys. Rev. Lett., 1986, 57, 2299. 10 A. D. Hammerich, R. Kosloff, H. Guo and M. Ratner, in preparation. 1 1 R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni and D. Tannor, Chem. Phys., 1989, 139, 201. 12 R. Kosloff, J. Phys. Chem., 1988, 92, 2087. 13 R. Koslop and H. Tal-Ezer, Chem. Phys. Lett., 1986, 127, 223. 14 H. Tal-Ezer, R. Kosloff and C. Cerjan, J. Comput. Phys., to be published. 15 H. Tal-Ezer and R. Kosloff, J. Chem. Phys., 1984,81, 2967. 16 Y. Amatatsu, K. Morokuma and S. Yabushita, J. Chem. Phys., in the press. 17 See e.g. G. D. Billing, Comput. Phys. Rep., 1990, 12, 383. 18 M. Berman, R. Kosloff and H. Tal-Ezer, submitted. 19 See e.g. B. J. Garrison and S. A. Adelman, Sur$ Sci., 1977, 66, 253; J. C. Tully, J. Chem. Phys., 1980, 20 Y. Zeiri and R. Kosloff, J. Chem. Phys., 1990, 93, 6890. 21 A. D. Hammerich, R. Kosloff and M. A. Ratner, Chem. Phys. Lett., 1990, 171, 97. 22 C. Cerjan, R. Kosloff, N. Bar-Tal, in preparation. 23 R. Kosloff, A. D. Hammerich and M. A. Ratner, in Large Finife Systems, 20th Jerusalem Symposium, 24 Z. Kotler, A. Nitzan and R. Kosloff, Chem. Phys. Lett., 1988, 153, 483. 25 N. Makri and W. H. Miller, J. Chem. Phys., 1987, 87, 5781. 26 H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73. 73, 1975; R. R. Lucchese and J. C. Tully, Surf: Sci., 1983, 137, 570. ed. J. Jortner, A. Pullman and B. Pullman, Reidel, Dordrecht, 1987, p. 53. Paper 1/00263E; Received 18th January, 1991.