- Email: [email protected]
Three-dimensional wave-packet dynamics on vibronically coupled dissociative potential energy surfaces
Volume 178, number I
CHEMICAL PHYSICS LETTERS
15 March 1991
Three-dimensional wave-packet dynamics on vibronically coupled dissociative potential energy surfaces U. Manthe and H. Kiippel Theoretirche Chemie, Physikalisch-Chemisches Institut, UniversitiitHeidelberg, W-6900Heidelberg, Germany Received IONovember 1990;in final form 17 December 1990
A scheme for calculating threedimensional wave-packet dynamics on conically intersecting Z- and H-surfaces of a triatomic molecule is presented. General potentials are employed and ali three nuclear degrees of freedom remaining after separation of the rotational motion (I=O) are included. Difficulties in the treatment of the angular coordinate due to the different centrifugal potentials on both surfaces are removed by a suitable transformation of the Hamiltonian. The scheme is applied to study numerically dissociation dynamics starting from an excited repulsive B electronic state which is vibronically coupled to a bound ITstate.
1.Introduction
Dissociation processes following photoexcitation or photoionisation are the subject of a considerable amount of experimental and theoretical work #’ [ 28]. Many theoretical techniques have been developed to describe these processes: classical trajectory and semiclassical methods [ 21, coupled-channel time-independent scattering methods [ 3 ] and timedependent wave-packet dynamics [ 4-81. In this Letter, a method for treating a dissociation process on coupled surfaces is presented. Three degrees of freedom are included in the calculation: two stretching modes and one bending mode. The majority of wave-packet-dynamics studies treats only one electronic surface [4,5 1. While in most calculations, this surface is two-dimensional, recently also calculations including three nuclear coordinates have been published [ 51. However, in many systems more than one electronic state contributes to the system dynamics. Time-dependent calculations of vibronically coupled systems employing structurally simple Hamiltonians have been performed for C2H: [ 91, pyrazine [ lo], CsH6+ [ 111 and NO, [ 121. Up to seven vibrational degrees of freedom are included in these calculations, but they are restricted to model type of Hamiltonians which *’ A recent review is given, for example, by ref. [ 11, 36
imply bound nuclear motion. Calculations on general coupled surfaces including dissociation were done for example for Hz [ 61, H2S and methyl iodide [ 71. All this latter work is restricted to two-dimensional surfaces. Only motion in two stretchingtype coordinates is studied, all other nuclear degrees of freedom are frozen. Motion in an angular coordinate has been investigated only for spin-orbit coupled systems with vanishing electronic angular momentum in each electronic state [ 8 1. The neglect of angular motion can yield only a crude description of the system dynamics. Especially, the coupling between the different electronic states typically results in a strong redistribution of energy between the different vibrational modes [ 9 I. Interacting Z and H states of linear molecules cannot be described without the angular coordinate anyway, because only the bending motion causes the vibronic coupling in such systems. In this work, a scheme for full three-dimensional wave-packet dynamics on conically intersecting Z and Il surfaces of linear triatomics is presented and a numerical example is given. To the best of our knowledge, this is the first three-dimensional calculation on general coupled surfaces. The inclusion of a bending coordinate allows to account for vibronic coupling effects by well established concepts [ 131. The importance of all three vibrational modes is shown in a model example: A linear triatomic molecule ABC
0009-2614/91/$ 03.50 0 1991 - Elsevier Science Publishers B.V. (North-Holland)
Volume 178, number 1
CHEMICALPHYSICSLETTERS
dissociates after photoionisation into the fragments A+ and BC. Only the J=O rotational state is investigated which allows the exact separation of rotational motion. The treatment is presented as follows: Section 2 gives a short description of the Hamiltonian. The numerical methodology is described in section 3. The treatment of the bending motion is discussed in some detail, since a new numerical scheme has to be developed due to the different angular momenta on the C-and II-surfaces. Finally, section 4 shows the initial numerical results.
2. The Hamiltonian In this section, we briefly describe the Hamiltonian of the triatomic which forms the basis for our investigation. A more complete derivation and discussion is found elsewhere [ 141. We consider a C and a II electronic state and a total (electronic plus nuclear) angular momentum J=O. The two components of the l-l state are chosen in such a way that one is symmetric with respect to the molecular plane of the triatomic (called II, in the following) and the other one is antisymmetric (IL). The antisymmetric component is coupled to the E and l-I+ states only via Renner-Teller coupling. Since this coupling is weak for J= 0 [ 13,151, the II_ state decouples to a good approximation and is dropped in the further treatment. Thus, the problem is reduced to two electronic states: C and II,. Following well established concepts [ 131, we use a diabatic electronic basis and neglect the dependence of the electronic wavefunctions on the nuclear coordinates. Our formulae are given in a vector notation for the wavefunction: the upper component represents the C state, the lower one the II+ state. The Hamiltonian is written by employing Jacobi coordinates for the positions of the nuclei in the bodyfixed frame: r, (v for vibrational coordinate) is the distance between the atoms B and C, rd (d for dissociative coordinate) the dis@ncebetween A and the center of mass of the diatomic BC and Bis the angle between the two corresponding axes. In this way, the overall rotational motion can be separated from the internal motion for J= 0 [ 16J. The conventional scheme for a single electronic surface can rather eas-
15March 1991
ily be generalized to the present situation of two surfaces with different electronic angular momenta [ 17] _ Neglecting one small term of minor importance [ 141, we arrive at the Hamiltonian, 1 a2 ---_r i a= H= ---r 2,&r, a($ d 2&r, ar: v 1 - --sinB-
a
a
2b sin 8 ae
ae (I)
The reduced masses and the moment of inertia are defined by 1 -= -&t&, P-V 1
-=Pe
l+l Pvrf
1 -=
L,‘,
fld
MA
pdri
’
b
(2)
Vx and Vn are the diabatic potential surfaces, Vzn is the non-diabatic or potential coupling.. An important difference between the Hamiltonian ( 1) and analogous Hamiltonians for systems with only one electronic surface is the appearance of the centrifugal potential,
(3) It is a consequence of the addition of nuclear and electronic angular momenta. For a E electronic state, a vanishing total angular momentum J implies also a vanishing nuclear angular momentum. But for a II state, the electronic angular momentum differs from zero; therefore, J=O results in a non-zero nuclear angular momentum and a corresponding centrifugal potential_
3. Collocation scheme for the angularcoordinate The time evolution of the system is studied by direct numerical integration of the time-dependent Schrijdinger equation. The propagation of the wave packet is done by a Ianczos time integrator [ 18] and a pseudospectral or collocation scheme [ 19,201 for the representation of the Hamiltonian. In the present section, we focus on the question how to apply the 37
Volume 178,number 1
collocation scheme to the angular coordinate. By the familiar transformation +r”‘rd’yl,
H=r,*r,H(r,T&’
,
(4)
the Hamiltonian ( 1) can be written in a form more convenient for the application of pseudospectral schemes, I?_
-
1a2
~2
2~ ar;
2~
+&%(:
i
ar:
--sin 2hsin
;)+(:n
2).
a.
eae
uation equals that for a single surface. But the tentrifugal potential V_,,,,+rueal (eq. (3) ) contributes an additional term for the II-surface. This term diverges for 0= 0. Its matrix elements in the Legendre basis, *
s .I
p,(c0s
0
a
0%
15)
In the reaction ABC+BC+A, the dissociation occurs in the r,j direction, while an oscillatory motion is found in the r, direction, We represent the Hamiltonian by employing basis sets (and corresponding grids) which are best adapted to the situation in the dissociation limit. The F’FT scheme of Kosloff and co-workers [ 191 is used for the rd coordinate. Numerical efficiency of the fast Fourier transform in case of large grid sizes and the plane wave basis set make this method the best choice for dissociative degrees of freedom. In the r, coordinate, we employ a harmonic oscillator discrete variable (DVR) scheme [ 201 with a frequency corresponding to the diatomic fragment BC, This method is attractive for this mode, because there the zeroth order description given by the harmonic oscillator basis is not too bad. Compared with an PFI scheme, a smaller number of basis functions and corresponding grid points can be used. Up to now, we have employed established methodology of wavefunction representation. Problems occur in the treatment of the bending angle. Standard (DVR-) methodology [ 52 11, used in the single-surface situation, represents the angular degree of freedom in a basis set spanned by Legendre polynomials P,(cos 8). Corresponding grid points are connected by Gaussian quadrature to these basis functions. In this basis, the angular differential operator is diagonal, namely,
(6) Let us now consider our Hamiltonian fl. This treatment is appropriate for the Z-surface, where the sit38
1SMarch 1991
CHEMICAL PHYSICS LETTERS
1 Pr ( cos e) sin 8 d0 , 2~~ sin%
e) -
also diverge. Therefore, the Legendre basis (and the corresponding grid) is unable to describe the dynamics on the II-surface. Numerical calculations show the lack of convergencein this basis. The problem can be resolved on the II-surface by employing a basis of spherical harmonics Y,,(0, +O) which diagonalizes the differential operator including the II-surface centrifugal term,
(
_-- ’ a sin 0 za sin 9 ae =1(/t
+ &ww
l)Y[l(B, 0) .
(7)
Using this basis, one encounters the same problem as above for the Z-surface. Therefore, different basis sets would be needed for the X- and II-surface. Since the grids corresponding to the &(cos 0) and the Yj,(6, 0) basis are different, a transformation from one grid to the other is complicated. But for the calculation of the off-diagonal elements of the potential matrix in the Hamiltonian (5), such a transformation would be necessary. We emphasize that these problems occur in every basis set and also in FFT schemes. They result from the inability to represent sin 8 in a basis of co03 functions, and vice versa. This problem can be overcome by transforming the wavefunction according to
-JO
Y-
(
1 I
0 sin-l@ ”
(8)
(In this and all following formulae, sin -“8 means (sin@-“.) No singularities in 9 result from the transformation, because the II-component of Q is, due to symmetry, proportional to sin 8. The transformed Hamiltonians read
_----1 a2 2bar:
1 iTI2
=-zS
The centrifugal potential is now absorbed in the angular kinetic energy term. Next, we turn to the representation of the wavefunction in the 8 direction. Because p is the 2x periodic, it can be expanded in a Fourier series,
1 2h
N-l p(fl,
sin-‘8 xK0
15March 1991
CHEMICALPHYSICSLE?TERs
Volume178,number 1
0 sin-*8)~sin8$(~
0
rd,
8)~ lim N-m
)
c
sine
*
(9)
Following standard FFT techniques [ 191, the series is truncated at a finite value of N. An equidistant grid (8,) is connected to the plane-wave basis via a FFT transformation:
n=-N
xexp(im8,/2x),
Due to this definition of the scalar product, the matrix elements of the centrifugal potential do not diverge any longer. All other terms in the Hamiltonian (9) also have regular matrix elements. Therefore, the transformation (8) removes the above problem: It avoids the representation of wavefunction components proportional to sin 0. Finally, by using the identity
- ’ +m(m+l) 38 sin”
, (Ill
the Hamiltonian (9) is simplified to w=
+
1 a2 2&z--@--
i
a2
V&n 8 V&/sin 8 V, *
V,
i
(12)
(14a)
N-l
xexp( -im6,/27c) ,
sin”+‘8 ae
.
(13)
We mention that also V&sin 8 has no singularity for 8= 0, because V,, is also proportional to sin 8 for symmetry reasons. The transformations (4) and (8) have changed the definition of the Hilbert-space scalar product. It is now given by
i a =- -_sin2m+l8?
&,(ry, ra) exp( -in8/2x)
,I=-N
e,=x(nti)/N.
(14b)
(14c)
The wavefunction is represented on the grid (8,; n=O, .... N- 1). As a result of the even symmetry of q, p(r,, rd, 8) =c(rv, rd, -8), it is not necessary to calculate the value of the wavefunction at negative 8 values. The multiplication with the potential matrix and the coordinate-dependent matrices in the kinetic energy expression of the Hamiltonian ( 11) is done in this grid representation. The two a/a8 op erators appearing in the Hamiltonian are evaluated employing the FFT transformations (14): with eq. ( 14a), the wavefunction is transformed to the planewave basis. The points in the interval ( -x, 0) needed in the formula arc given by the symmetry relation I( r, , rd, e) = Q( r,, rd, - 8). In the plane-wave basis, the action of the a/a8 operator is calculated easily. After its evaluation, the wavefunction is transformed back to the grid by eq. ( 14b). Again, this operation is done only for positive 8 values. In the above procedure, the whole interval of angles, [ -R, x], is covered by the grid. In many cases, this is not necessary. The wavefunction takes nonnegligible values only in an interval [ -it&,,, t9,,,,,] with 8,,,msx
Volume 178, number 1
CHEMICAL PHYSICS LEl’TERS
is done by replacing all ICby f3,, in eqs. ( 13) and ( 14). Such a procedure allows one to diminish N while keeping the spacing of the grid points constant.
4. Numerical example In section 3, a scheme for propagation on coupled general potential surfaces has been developed. NoiK, in this section, a first exemplary calculation is presented. Dissociation on vibronically coupled surfaces is quite a complex process. To get a first insight, we restrict ourselves to a physically appealing model potential rather than presenting a particular example. The potentials and potential couplings of our model read in terms of the internuclear distances r&.B,rBc and the bond angle a, T/n=tluvo~(r~c-Rn)2+f~~dW~(r*B--R0)2
v,,=Asin(a)
exp[-&S(rAB-&)l
.
(15a)
The ground-state surface of the neutral molecule ABC is
t&w:(
1-cos a) .
(15b)
All parameters are displayed in table 1. These potentials model important aspects of the HCN system. The symbols A, B, C denote the atoms H, C,N, respectively as indicated by their masses given in table I. The potential V, is taken to present HCN in its ground state, while Vz and Vn describe the ion
15 March 1991
with a hole in the C-H o-bond and in the C-N ICbond, respectively. A more complete discussion is given elsewhere [ 141. Our numerical method is based on Jacobi coordinates; therefore, the potentials ( 15) have to be trat%fOrtWd from the bond coordinates rAa, r,, (Y to the Jacobi coordinates r,, r& 8. This transformation yields complicated formulae for Vz, Vnand V,,. We will not give these cumbersome expressions explicitly, but we stress the fact that there is no special form of the potentials which allows a simplification of the numerical algorithm, Propagation on coupled general non-harmonic potentials has to be calculated. Some technical details of the calculation should be mentioned now. The initially prepared state is assumed to be the ground state of the neutral molecule ABC vertically shifted to the E-surface of the ion [ABC]+. It is numerically calculated by employing the ground-state potential surface [ 141. The time integration is done by a Lanczos-type integrator with variable step size [ 141. In order to obtain converged results, it proved necessary to employ grid sizes of up to 96 and 48 in the FFT scheme for rd and 0, respectively, and 20 in the harmonic oscillator basis employed for r,. Including the electronic degree of freedom, this amounts to a length of the state vector of 184320. Figs. l-3 show important characteristics of the evolution of the wave packet. In fig. 1, the probability density to find the system at a defined value of r,, II [ tyl 2sin 8rt dr, de, is plotted versus r,. The potentials of the E and ll state are indicated for 8=0 and rv=RZ. During the first few femtoseconds, the initial wave packet moves downwards on the C-surface where it is originally located, e.g. the A-B bond length increases. Due to the form of the potential energy surface, the wave packet broadens in parallel. At around 4 fs, the conical intersection is reached. There
Table 1 System parameters for the example shown in figs. l-3 m,=l u R~=1_156~10-~~m ~=4.102~10’~s-’ Emin=- 1 ev o OI =8433~1O’~s-’ .
mB=12u Rn=1.216x10-10m ~=3.494x10i4s-’ Erun= 1 eV Jihs=4,OOOX 10’Om-’
J~0[~‘(R~+R~~/m~)-~+~~~~*]-‘=1.916X10~~~~m~
40
mc= 14u I&= 1.064x lo-‘” m od=6.048x lOI s-’ I=leV
Volume 178. number 1
-11
1.0
CHEMICAL PHYSICS LETTERS
1.5
2.0
2.5
3.0
I5 March 1991
1 oo” Fig. 3. The probability density to End tbe nondissociating part (see text ) of the wave packet at a defined value of fl is displayed at 10 fs (full line), 20 fs (dashed line) and 30 fs (dotted line).
2
1
3
5
4
Fig. 1. Probability density to End tbe system at a defined value of rd ( rdin lB-‘” m). Dashed tines show the wave packet at 0 fs (A) and 10 fs (B), dotted at 5 fs (A) and 15 fs (B). Thepotentials of the Z and TIstates are displayed asfultlines.Thescaleof the y axis indicatesthe energies (in eV) for the potentials. Tbe units for the probability density are arbitrary but are of the same scale relative to the energies in both panels.
1. I
0.9
I
1.0
.
“““I
1.1
,
1.2
1.3
1.4
Fig. 2. Probability density to find the non-dissociating part (see text) of the wave packet at a defined value of r, (rv in lo-l0 m). The wave packet is displayed at 10 fs (full line), 20 fs (dashed line) and 30 fs (dotted line).
the main part of the wave packet stays in the C state. A minor part of around 6% interconverts to the IIsurface. The part of the wave packet located on the C-surface undergoes a direct dissociation. Its angular distribution is only slightly changed by the vibronic coupling with the II-surface, and its other behavior remains virtually unchanged. An explicit comparison and further details will be given in following work. More interesting is the part of the wave packet which has interconverted to the electronic II state. It does not dissociate directly (see fig. 1B) but moves back to smaller values of the dissociative coordinate rd. This leads to a clear separation of the wave packet into a directly dissociating and a (so far) non-dissociating part. Strong multimode effects are found in the non-dissociating part. The equilibrium distance of the BC bond is larger in the II than in the 1 state. Therefore, the BC stretching vibration is excited. Oscillatory motion of the wave packet is found in this direction (fig. 2 ) . Since the wave packet moves now on the lower adiabatic potential surface, the doubleminimum form of this surface leads to a strong broadening of the wave packet in the 8 direction (tig. 3). Such a behavior is typical for the motion in the coupling coordinate on conically intersecting surfaces [221. This latter part of the wave packet dissociates quite slowly. Describing the situation in diabatic electronic states, the packet has to change from the II to the C state in order to dissociate. This delayed dis41
Volume 178, number 1
CHEMICAL PHYSICS LETTERS
sociation process occurs on a timescale of more than a few ten femtoseconds.
5. Conclusions A method has been developed which allows one to calculate wave-packet dynamics on coupled general potential energy surfaces. Contrary to known numerical schemes, different centrifugal potentials for each surface can be treated. A linear triatomic molecule has been investigated numerically, with all three vibrational modes being included in the calculation. This is, to the best of our knowledge, the first threedimensional calculation with coupled electronic states and general potentials. The results show the importance of all three vibrational modes for the system dynamics. Pronounced motion in the angular coordinate is found. This motion was completely neglected by all earlier studies which were restricted to two dimensions. Only the very first results are shown in this Letter; a more complete presentation of results and an extended discussion will be given in a subsequent article [ 14 1.
Acknowledgement The authors are indebted to L.S. Cederbaum and H.-D. Meyer for stimulating discussions and their interest in this work. One of us (UM) acknowledges financial support by the Verband der Chemischen Industrie.
References [I] S. Mukamel, Ann. Rev. Phys. Chem. 41 (1990) 647; J.P. Simons, J. Phys. Chem. 90 (1987) 5378. [Z] E.J. Heller, J. Chem. Phys. 62 (i975) 1544; M. Desouter-Lecomte, C. Galley, J.C. Lorquet and M. Vaz Pires, J. Chem. Phys. 71 (1979) 3661; R. Schinke, J. Phys. Chem. 90 (1986) 1742; N.E. Hendriksen and E.J. Heller, J. Chem. Phys. 91 ( 1989) 4700.
42
15 March 1991
[ 31 K.C. Kulander and J.C. Lit, J. Chem. Phys. 73 ( 1980) 4337; DC. Clary, J. Chem. Phys. 84 (1986) 4288. [4] K. Weide and R. Schinke, J. Chem. Phys. 90 (1989) 7150; V. Engel, R. Schinke, S. Hennig and H. Metiu, J. Chem. Phys. 92 (1990) 1; SE. Bradforth, A. Weaver, D.W. Arnold, R.B. Metz and D.M. Neumark, J. Chem. Phys. 92 (1990) 7205. [5] SK. Gray and C.E. Wozny, J. Chem. Phys. 91 (1989) 7671; F. Le Querc and C. Leforcstier, J. Chem. Phys. 92 (1990) 247. 16lA.E. Ore1 and K.C. Kulander, Chem. Phys. Letters 146 (1988) 428; X.P. Jiang, R. Heather and H. Metiu, J. Chem. Phys. 90 (1989) 2555. [7] K. Wcide, V. StimmIer and R. Schinke, J. Chem. Phys. 93 (1990) 861; H. Guo and G.C. Schatz, J. Chem. Phys. 93 ( 1990) 393. [ 81 H. Guo and G.C. Schatz, J. Chem. Phys. 92 ( 1990) 1634. [9] H. Kiippel, Chem. Phys. 77 (1983) 359. [IO] R. Schneider and W. Domcke, Chem. Phys. Letters 150 (1988) 235; R. Schneider, W. Domcke and H. Kiippel, J. Chem. Phys. 92 (1990) 1045. [ 111H. Kiippel, L.S. Cederbaum and W. Domcke, J. Cbem. Phys. 88 (1988) 2023. [ 121U. Mantheand H. Ktippel, J. Chem. Phys. 93 (1990) 1658. [ 131H. Kiippel, W. Domcke and L.S. Cederbaum, Advan. Cbem. Phys. 57 (1984) 59. [ 141U. Manthe, H. Kiippeland L.S. Cederbaum, to be published. [IS] Ch. Jungen and A.J. Mew, Mol. Phys. 40 (1984) I; M. Peric, S.D, Peyerimhoff and R.J. Buenker, Intern. Rev. Phys. Chem. 4 (1985) 85; L.J. Dunne, H. Guo and J.N. Murrell, Mol. Phys. 62 (1987) 283. [ 161C.F. Curtiss, J.O. Hirschfelder and ET. Adler, J. Chem. Phys. 18 (1950) 1638. [ 171C. Petrongolo, J. Chem. Phys. 89 (1988) 1297; B.T. Sutcliffe, in: Studies in physical and theoretical chemistry, Vol. 2 1. Current aspects of quantum chemistry, ed. R. Garbo (Elsevier, Amsterdam, 1981) p. 99. [18]T.J.ParkandJ.C.Light,J.Chem.Phys.85 (1986)5870. [ 191 D. Kosloff and R. Kosloff, J. Comput. Phys. 52 (1983) 35. [20] D.O. Harris, G.G. Engerholm and W.D. Gwinn, J. C&em. Phys.43 (1965) 1515; A.S. Dickinson andP.R. Certain, J. Chem. Phys. 49 ( 1968) 4209; J.V. Lill, G.A. Parker and J.C. Light, J. Chem. Phys. 82 (1985) 1400. [21] Z. BaEiCand J.C. Light, J. Chem. Phys. 85 (1986) 4594. [22] U. Manthe and H. Kiippel, J. Chem. Phys. 93 (1990) 345.
CHEMICAL PHYSICS LETTERS
15 March 1991
Three-dimensional wave-packet dynamics on vibronically coupled dissociative potential energy surfaces U. Manthe and H. Kiippel Theoretirche Chemie, Physikalisch-Chemisches Institut, UniversitiitHeidelberg, W-6900Heidelberg, Germany Received IONovember 1990;in final form 17 December 1990
A scheme for calculating threedimensional wave-packet dynamics on conically intersecting Z- and H-surfaces of a triatomic molecule is presented. General potentials are employed and ali three nuclear degrees of freedom remaining after separation of the rotational motion (I=O) are included. Difficulties in the treatment of the angular coordinate due to the different centrifugal potentials on both surfaces are removed by a suitable transformation of the Hamiltonian. The scheme is applied to study numerically dissociation dynamics starting from an excited repulsive B electronic state which is vibronically coupled to a bound ITstate.
1.Introduction
Dissociation processes following photoexcitation or photoionisation are the subject of a considerable amount of experimental and theoretical work #’ [ 28]. Many theoretical techniques have been developed to describe these processes: classical trajectory and semiclassical methods [ 21, coupled-channel time-independent scattering methods [ 3 ] and timedependent wave-packet dynamics [ 4-81. In this Letter, a method for treating a dissociation process on coupled surfaces is presented. Three degrees of freedom are included in the calculation: two stretching modes and one bending mode. The majority of wave-packet-dynamics studies treats only one electronic surface [4,5 1. While in most calculations, this surface is two-dimensional, recently also calculations including three nuclear coordinates have been published [ 51. However, in many systems more than one electronic state contributes to the system dynamics. Time-dependent calculations of vibronically coupled systems employing structurally simple Hamiltonians have been performed for C2H: [ 91, pyrazine [ lo], CsH6+ [ 111 and NO, [ 121. Up to seven vibrational degrees of freedom are included in these calculations, but they are restricted to model type of Hamiltonians which *’ A recent review is given, for example, by ref. [ 11, 36
imply bound nuclear motion. Calculations on general coupled surfaces including dissociation were done for example for Hz [ 61, H2S and methyl iodide [ 71. All this latter work is restricted to two-dimensional surfaces. Only motion in two stretchingtype coordinates is studied, all other nuclear degrees of freedom are frozen. Motion in an angular coordinate has been investigated only for spin-orbit coupled systems with vanishing electronic angular momentum in each electronic state [ 8 1. The neglect of angular motion can yield only a crude description of the system dynamics. Especially, the coupling between the different electronic states typically results in a strong redistribution of energy between the different vibrational modes [ 9 I. Interacting Z and H states of linear molecules cannot be described without the angular coordinate anyway, because only the bending motion causes the vibronic coupling in such systems. In this work, a scheme for full three-dimensional wave-packet dynamics on conically intersecting Z and Il surfaces of linear triatomics is presented and a numerical example is given. To the best of our knowledge, this is the first three-dimensional calculation on general coupled surfaces. The inclusion of a bending coordinate allows to account for vibronic coupling effects by well established concepts [ 131. The importance of all three vibrational modes is shown in a model example: A linear triatomic molecule ABC
0009-2614/91/$ 03.50 0 1991 - Elsevier Science Publishers B.V. (North-Holland)
Volume 178, number 1
CHEMICALPHYSICSLETTERS
dissociates after photoionisation into the fragments A+ and BC. Only the J=O rotational state is investigated which allows the exact separation of rotational motion. The treatment is presented as follows: Section 2 gives a short description of the Hamiltonian. The numerical methodology is described in section 3. The treatment of the bending motion is discussed in some detail, since a new numerical scheme has to be developed due to the different angular momenta on the C-and II-surfaces. Finally, section 4 shows the initial numerical results.
2. The Hamiltonian In this section, we briefly describe the Hamiltonian of the triatomic which forms the basis for our investigation. A more complete derivation and discussion is found elsewhere [ 141. We consider a C and a II electronic state and a total (electronic plus nuclear) angular momentum J=O. The two components of the l-l state are chosen in such a way that one is symmetric with respect to the molecular plane of the triatomic (called II, in the following) and the other one is antisymmetric (IL). The antisymmetric component is coupled to the E and l-I+ states only via Renner-Teller coupling. Since this coupling is weak for J= 0 [ 13,151, the II_ state decouples to a good approximation and is dropped in the further treatment. Thus, the problem is reduced to two electronic states: C and II,. Following well established concepts [ 131, we use a diabatic electronic basis and neglect the dependence of the electronic wavefunctions on the nuclear coordinates. Our formulae are given in a vector notation for the wavefunction: the upper component represents the C state, the lower one the II+ state. The Hamiltonian is written by employing Jacobi coordinates for the positions of the nuclei in the bodyfixed frame: r, (v for vibrational coordinate) is the distance between the atoms B and C, rd (d for dissociative coordinate) the dis@ncebetween A and the center of mass of the diatomic BC and Bis the angle between the two corresponding axes. In this way, the overall rotational motion can be separated from the internal motion for J= 0 [ 16J. The conventional scheme for a single electronic surface can rather eas-
15March 1991
ily be generalized to the present situation of two surfaces with different electronic angular momenta [ 17] _ Neglecting one small term of minor importance [ 141, we arrive at the Hamiltonian, 1 a2 ---_r i a= H= ---r 2,&r, a($ d 2&r, ar: v 1 - --sinB-
a
a
2b sin 8 ae
ae (I)
The reduced masses and the moment of inertia are defined by 1 -= -&t&, P-V 1
-=Pe
l+l Pvrf
1 -=
L,‘,
fld
MA
pdri
’
b
(2)
Vx and Vn are the diabatic potential surfaces, Vzn is the non-diabatic or potential coupling.. An important difference between the Hamiltonian ( 1) and analogous Hamiltonians for systems with only one electronic surface is the appearance of the centrifugal potential,
(3) It is a consequence of the addition of nuclear and electronic angular momenta. For a E electronic state, a vanishing total angular momentum J implies also a vanishing nuclear angular momentum. But for a II state, the electronic angular momentum differs from zero; therefore, J=O results in a non-zero nuclear angular momentum and a corresponding centrifugal potential_
3. Collocation scheme for the angularcoordinate The time evolution of the system is studied by direct numerical integration of the time-dependent Schrijdinger equation. The propagation of the wave packet is done by a Ianczos time integrator [ 18] and a pseudospectral or collocation scheme [ 19,201 for the representation of the Hamiltonian. In the present section, we focus on the question how to apply the 37
Volume 178,number 1
collocation scheme to the angular coordinate. By the familiar transformation +r”‘rd’yl,
H=r,*r,H(r,T&’
,
(4)
the Hamiltonian ( 1) can be written in a form more convenient for the application of pseudospectral schemes, I?_
-
1a2
~2
2~ ar;
2~
+&%(:
i
ar:
--sin 2hsin
;)+(:n
2).
a.
eae
uation equals that for a single surface. But the tentrifugal potential V_,,,,+rueal (eq. (3) ) contributes an additional term for the II-surface. This term diverges for 0= 0. Its matrix elements in the Legendre basis, *
s .I
p,(c0s
0
a
0%
15)
In the reaction ABC+BC+A, the dissociation occurs in the r,j direction, while an oscillatory motion is found in the r, direction, We represent the Hamiltonian by employing basis sets (and corresponding grids) which are best adapted to the situation in the dissociation limit. The F’FT scheme of Kosloff and co-workers [ 191 is used for the rd coordinate. Numerical efficiency of the fast Fourier transform in case of large grid sizes and the plane wave basis set make this method the best choice for dissociative degrees of freedom. In the r, coordinate, we employ a harmonic oscillator discrete variable (DVR) scheme [ 201 with a frequency corresponding to the diatomic fragment BC, This method is attractive for this mode, because there the zeroth order description given by the harmonic oscillator basis is not too bad. Compared with an PFI scheme, a smaller number of basis functions and corresponding grid points can be used. Up to now, we have employed established methodology of wavefunction representation. Problems occur in the treatment of the bending angle. Standard (DVR-) methodology [ 52 11, used in the single-surface situation, represents the angular degree of freedom in a basis set spanned by Legendre polynomials P,(cos 8). Corresponding grid points are connected by Gaussian quadrature to these basis functions. In this basis, the angular differential operator is diagonal, namely,
(6) Let us now consider our Hamiltonian fl. This treatment is appropriate for the Z-surface, where the sit38
1SMarch 1991
CHEMICAL PHYSICS LETTERS
1 Pr ( cos e) sin 8 d0 , 2~~ sin%
e) -
also diverge. Therefore, the Legendre basis (and the corresponding grid) is unable to describe the dynamics on the II-surface. Numerical calculations show the lack of convergencein this basis. The problem can be resolved on the II-surface by employing a basis of spherical harmonics Y,,(0, +O) which diagonalizes the differential operator including the II-surface centrifugal term,
(
_-- ’ a sin 0 za sin 9 ae =1(/t
+ &ww
l)Y[l(B, 0) .
(7)
Using this basis, one encounters the same problem as above for the Z-surface. Therefore, different basis sets would be needed for the X- and II-surface. Since the grids corresponding to the &(cos 0) and the Yj,(6, 0) basis are different, a transformation from one grid to the other is complicated. But for the calculation of the off-diagonal elements of the potential matrix in the Hamiltonian (5), such a transformation would be necessary. We emphasize that these problems occur in every basis set and also in FFT schemes. They result from the inability to represent sin 8 in a basis of co03 functions, and vice versa. This problem can be overcome by transforming the wavefunction according to
-JO
Y-
(
1 I
0 sin-l@ ”
(8)
(In this and all following formulae, sin -“8 means (sin@-“.) No singularities in 9 result from the transformation, because the II-component of Q is, due to symmetry, proportional to sin 8. The transformed Hamiltonians read
_----1 a2 2bar:
1 iTI2
=-zS
The centrifugal potential is now absorbed in the angular kinetic energy term. Next, we turn to the representation of the wavefunction in the 8 direction. Because p is the 2x periodic, it can be expanded in a Fourier series,
1 2h
N-l p(fl,
sin-‘8 xK0
15March 1991
CHEMICALPHYSICSLE?TERs
Volume178,number 1
0 sin-*8)~sin8$(~
0
rd,
8)~ lim N-m
)
c
sine
*
(9)
Following standard FFT techniques [ 191, the series is truncated at a finite value of N. An equidistant grid (8,) is connected to the plane-wave basis via a FFT transformation:
n=-N
xexp(im8,/2x),
Due to this definition of the scalar product, the matrix elements of the centrifugal potential do not diverge any longer. All other terms in the Hamiltonian (9) also have regular matrix elements. Therefore, the transformation (8) removes the above problem: It avoids the representation of wavefunction components proportional to sin 0. Finally, by using the identity
- ’ +m(m+l) 38 sin”
, (Ill
the Hamiltonian (9) is simplified to w=
+
1 a2 2&z--@--
i
a2
V&n 8 V&/sin 8 V, *
V,
i
(12)
(14a)
N-l
xexp( -im6,/27c) ,
sin”+‘8 ae
.
(13)
We mention that also V&sin 8 has no singularity for 8= 0, because V,, is also proportional to sin 8 for symmetry reasons. The transformations (4) and (8) have changed the definition of the Hilbert-space scalar product. It is now given by
i a =- -_sin2m+l8?
&,(ry, ra) exp( -in8/2x)
,I=-N
e,=x(nti)/N.
(14b)
(14c)
The wavefunction is represented on the grid (8,; n=O, .... N- 1). As a result of the even symmetry of q, p(r,, rd, 8) =c(rv, rd, -8), it is not necessary to calculate the value of the wavefunction at negative 8 values. The multiplication with the potential matrix and the coordinate-dependent matrices in the kinetic energy expression of the Hamiltonian ( 11) is done in this grid representation. The two a/a8 op erators appearing in the Hamiltonian are evaluated employing the FFT transformations (14): with eq. ( 14a), the wavefunction is transformed to the planewave basis. The points in the interval ( -x, 0) needed in the formula arc given by the symmetry relation I( r, , rd, e) = Q( r,, rd, - 8). In the plane-wave basis, the action of the a/a8 operator is calculated easily. After its evaluation, the wavefunction is transformed back to the grid by eq. ( 14b). Again, this operation is done only for positive 8 values. In the above procedure, the whole interval of angles, [ -R, x], is covered by the grid. In many cases, this is not necessary. The wavefunction takes nonnegligible values only in an interval [ -it&,,, t9,,,,,] with 8,,,msx
Volume 178, number 1
CHEMICAL PHYSICS LEl’TERS
is done by replacing all ICby f3,, in eqs. ( 13) and ( 14). Such a procedure allows one to diminish N while keeping the spacing of the grid points constant.
4. Numerical example In section 3, a scheme for propagation on coupled general potential surfaces has been developed. NoiK, in this section, a first exemplary calculation is presented. Dissociation on vibronically coupled surfaces is quite a complex process. To get a first insight, we restrict ourselves to a physically appealing model potential rather than presenting a particular example. The potentials and potential couplings of our model read in terms of the internuclear distances r&.B,rBc and the bond angle a, T/n=tluvo~(r~c-Rn)2+f~~dW~(r*B--R0)2
v,,=Asin(a)
exp[-&S(rAB-&)l
.
(15a)
The ground-state surface of the neutral molecule ABC is
t&w:(
1-cos a) .
(15b)
All parameters are displayed in table 1. These potentials model important aspects of the HCN system. The symbols A, B, C denote the atoms H, C,N, respectively as indicated by their masses given in table I. The potential V, is taken to present HCN in its ground state, while Vz and Vn describe the ion
15 March 1991
with a hole in the C-H o-bond and in the C-N ICbond, respectively. A more complete discussion is given elsewhere [ 141. Our numerical method is based on Jacobi coordinates; therefore, the potentials ( 15) have to be trat%fOrtWd from the bond coordinates rAa, r,, (Y to the Jacobi coordinates r,, r& 8. This transformation yields complicated formulae for Vz, Vnand V,,. We will not give these cumbersome expressions explicitly, but we stress the fact that there is no special form of the potentials which allows a simplification of the numerical algorithm, Propagation on coupled general non-harmonic potentials has to be calculated. Some technical details of the calculation should be mentioned now. The initially prepared state is assumed to be the ground state of the neutral molecule ABC vertically shifted to the E-surface of the ion [ABC]+. It is numerically calculated by employing the ground-state potential surface [ 141. The time integration is done by a Lanczos-type integrator with variable step size [ 141. In order to obtain converged results, it proved necessary to employ grid sizes of up to 96 and 48 in the FFT scheme for rd and 0, respectively, and 20 in the harmonic oscillator basis employed for r,. Including the electronic degree of freedom, this amounts to a length of the state vector of 184320. Figs. l-3 show important characteristics of the evolution of the wave packet. In fig. 1, the probability density to find the system at a defined value of r,, II [ tyl 2sin 8rt dr, de, is plotted versus r,. The potentials of the E and ll state are indicated for 8=0 and rv=RZ. During the first few femtoseconds, the initial wave packet moves downwards on the C-surface where it is originally located, e.g. the A-B bond length increases. Due to the form of the potential energy surface, the wave packet broadens in parallel. At around 4 fs, the conical intersection is reached. There
Table 1 System parameters for the example shown in figs. l-3 m,=l u R~=1_156~10-~~m ~=4.102~10’~s-’ Emin=- 1 ev o OI =8433~1O’~s-’ .
mB=12u Rn=1.216x10-10m ~=3.494x10i4s-’ Erun= 1 eV Jihs=4,OOOX 10’Om-’
J~0[~‘(R~+R~~/m~)-~+~~~~*]-‘=1.916X10~~~~m~
40
mc= 14u I&= 1.064x lo-‘” m od=6.048x lOI s-’ I=leV
Volume 178. number 1
-11
1.0
CHEMICAL PHYSICS LETTERS
1.5
2.0
2.5
3.0
I5 March 1991
1 oo” Fig. 3. The probability density to End tbe nondissociating part (see text ) of the wave packet at a defined value of fl is displayed at 10 fs (full line), 20 fs (dashed line) and 30 fs (dotted line).
2
1
3
5
4
Fig. 1. Probability density to End tbe system at a defined value of rd ( rdin lB-‘” m). Dashed tines show the wave packet at 0 fs (A) and 10 fs (B), dotted at 5 fs (A) and 15 fs (B). Thepotentials of the Z and TIstates are displayed asfultlines.Thescaleof the y axis indicatesthe energies (in eV) for the potentials. Tbe units for the probability density are arbitrary but are of the same scale relative to the energies in both panels.
1. I
0.9
I
1.0
.
“““I
1.1
,
1.2
1.3
1.4
Fig. 2. Probability density to find the non-dissociating part (see text) of the wave packet at a defined value of r, (rv in lo-l0 m). The wave packet is displayed at 10 fs (full line), 20 fs (dashed line) and 30 fs (dotted line).
the main part of the wave packet stays in the C state. A minor part of around 6% interconverts to the IIsurface. The part of the wave packet located on the C-surface undergoes a direct dissociation. Its angular distribution is only slightly changed by the vibronic coupling with the II-surface, and its other behavior remains virtually unchanged. An explicit comparison and further details will be given in following work. More interesting is the part of the wave packet which has interconverted to the electronic II state. It does not dissociate directly (see fig. 1B) but moves back to smaller values of the dissociative coordinate rd. This leads to a clear separation of the wave packet into a directly dissociating and a (so far) non-dissociating part. Strong multimode effects are found in the non-dissociating part. The equilibrium distance of the BC bond is larger in the II than in the 1 state. Therefore, the BC stretching vibration is excited. Oscillatory motion of the wave packet is found in this direction (fig. 2 ) . Since the wave packet moves now on the lower adiabatic potential surface, the doubleminimum form of this surface leads to a strong broadening of the wave packet in the 8 direction (tig. 3). Such a behavior is typical for the motion in the coupling coordinate on conically intersecting surfaces [221. This latter part of the wave packet dissociates quite slowly. Describing the situation in diabatic electronic states, the packet has to change from the II to the C state in order to dissociate. This delayed dis41
Volume 178, number 1
CHEMICAL PHYSICS LETTERS
sociation process occurs on a timescale of more than a few ten femtoseconds.
5. Conclusions A method has been developed which allows one to calculate wave-packet dynamics on coupled general potential energy surfaces. Contrary to known numerical schemes, different centrifugal potentials for each surface can be treated. A linear triatomic molecule has been investigated numerically, with all three vibrational modes being included in the calculation. This is, to the best of our knowledge, the first threedimensional calculation with coupled electronic states and general potentials. The results show the importance of all three vibrational modes for the system dynamics. Pronounced motion in the angular coordinate is found. This motion was completely neglected by all earlier studies which were restricted to two dimensions. Only the very first results are shown in this Letter; a more complete presentation of results and an extended discussion will be given in a subsequent article [ 14 1.
Acknowledgement The authors are indebted to L.S. Cederbaum and H.-D. Meyer for stimulating discussions and their interest in this work. One of us (UM) acknowledges financial support by the Verband der Chemischen Industrie.
References [I] S. Mukamel, Ann. Rev. Phys. Chem. 41 (1990) 647; J.P. Simons, J. Phys. Chem. 90 (1987) 5378. [Z] E.J. Heller, J. Chem. Phys. 62 (i975) 1544; M. Desouter-Lecomte, C. Galley, J.C. Lorquet and M. Vaz Pires, J. Chem. Phys. 71 (1979) 3661; R. Schinke, J. Phys. Chem. 90 (1986) 1742; N.E. Hendriksen and E.J. Heller, J. Chem. Phys. 91 ( 1989) 4700.
42
15 March 1991
[ 31 K.C. Kulander and J.C. Lit, J. Chem. Phys. 73 ( 1980) 4337; DC. Clary, J. Chem. Phys. 84 (1986) 4288. [4] K. Weide and R. Schinke, J. Chem. Phys. 90 (1989) 7150; V. Engel, R. Schinke, S. Hennig and H. Metiu, J. Chem. Phys. 92 (1990) 1; SE. Bradforth, A. Weaver, D.W. Arnold, R.B. Metz and D.M. Neumark, J. Chem. Phys. 92 (1990) 7205. [5] SK. Gray and C.E. Wozny, J. Chem. Phys. 91 (1989) 7671; F. Le Querc and C. Leforcstier, J. Chem. Phys. 92 (1990) 247. 16lA.E. Ore1 and K.C. Kulander, Chem. Phys. Letters 146 (1988) 428; X.P. Jiang, R. Heather and H. Metiu, J. Chem. Phys. 90 (1989) 2555. [7] K. Wcide, V. StimmIer and R. Schinke, J. Chem. Phys. 93 (1990) 861; H. Guo and G.C. Schatz, J. Chem. Phys. 93 ( 1990) 393. [ 81 H. Guo and G.C. Schatz, J. Chem. Phys. 92 ( 1990) 1634. [9] H. Kiippel, Chem. Phys. 77 (1983) 359. [IO] R. Schneider and W. Domcke, Chem. Phys. Letters 150 (1988) 235; R. Schneider, W. Domcke and H. Kiippel, J. Chem. Phys. 92 (1990) 1045. [ 111H. Kiippel, L.S. Cederbaum and W. Domcke, J. Cbem. Phys. 88 (1988) 2023. [ 121U. Mantheand H. Ktippel, J. Chem. Phys. 93 (1990) 1658. [ 131H. Kiippel, W. Domcke and L.S. Cederbaum, Advan. Cbem. Phys. 57 (1984) 59. [ 141U. Manthe, H. Kiippeland L.S. Cederbaum, to be published. [IS] Ch. Jungen and A.J. Mew, Mol. Phys. 40 (1984) I; M. Peric, S.D, Peyerimhoff and R.J. Buenker, Intern. Rev. Phys. Chem. 4 (1985) 85; L.J. Dunne, H. Guo and J.N. Murrell, Mol. Phys. 62 (1987) 283. [ 161C.F. Curtiss, J.O. Hirschfelder and ET. Adler, J. Chem. Phys. 18 (1950) 1638. [ 171C. Petrongolo, J. Chem. Phys. 89 (1988) 1297; B.T. Sutcliffe, in: Studies in physical and theoretical chemistry, Vol. 2 1. Current aspects of quantum chemistry, ed. R. Garbo (Elsevier, Amsterdam, 1981) p. 99. [18]T.J.ParkandJ.C.Light,J.Chem.Phys.85 (1986)5870. [ 191 D. Kosloff and R. Kosloff, J. Comput. Phys. 52 (1983) 35. [20] D.O. Harris, G.G. Engerholm and W.D. Gwinn, J. C&em. Phys.43 (1965) 1515; A.S. Dickinson andP.R. Certain, J. Chem. Phys. 49 ( 1968) 4209; J.V. Lill, G.A. Parker and J.C. Light, J. Chem. Phys. 82 (1985) 1400. [21] Z. BaEiCand J.C. Light, J. Chem. Phys. 85 (1986) 4594. [22] U. Manthe and H. Kiippel, J. Chem. Phys. 93 (1990) 345.