- Email: [email protected]
Classical and quantum dynamics on the collinear potential energy surface for the reaction of Li with H2
Chemical Physics 233 Ž1998. 9–27
Classical and quantum dynamics on the collinear potential energy surface for the reaction of Li with H 2 Nick J. Clarke a,1, Maurizio Sironi a , Mario Raimondi a , Sanjay Kumar b, Franco A. Gianturco b,2 , Erasmo Buonomo c , David L. Cooper d a
Dipartimento di Chimica Fisica ed Elettrochimica, UniÕersita´ di Milano, e Centro CNR-CSRSRC, Via Golgi 19, 20133 Milano, Italy b Department of Chemistry, The UniÕersity of Rome, Citta´ UniÕersitaria, 00185 Rome, Italy c Center for Supercomputing Applications to Research (CASPUR), Citta´ UniÕersitaria, 00185 Rome, Italy d Department of Chemistry, UniÕersity of LiÕerpool, P.O. Box 147, LiÕerpool L69 7ZD, UK Received 22 December 1997
Abstract Extensive Valence Bond ŽVB. calculations using non-orthogonal, spin-coupled wavefunctions with optimized orbitals have been employed to analyse the reactive region for the title process. Both the exothermic channel and the reverse, endothermic process are of strong astrophysical interest although no previous calculations have been available on both the reactive dynamics and the interaction energy surface. The specific features of the potential are analysed for some indicative configurations and classical trajectory calculations are carried out for the special collinear arrangement. In the latter instance, quantum time-dependent wavepacket calculations have also been performed and the two sets of results are found to be in rather good accord with each other. Some consequences of these exploratory calculations are discussed in detail. q 1998 Elsevier Science B.V. All rights reserved.
1. Introduction In the general scenario of what is known as ‘‘Early Universe’’ chemistry, the formation and depletion of LiH and LiHq is known to play a rather important role w1–3x since such species, together with H, Hq, H 2 , Hq 2 and others, can induce in the primordial gas spectral distortions and spatial anisotropies of the Cosmic Background Radiation ŽCBR.. Thus, the possible presence of molecules in such region can act as a heat pump, absorbing energy
in their rotovibrational levels and re-emitting this energy at lower frequencies by undergoing rotational cooling which will thus cause deviations of the CBR spectrum from the pure Planck spectrum w4x. It is usually assumed that the lithium chemistry was initiated, as the universe expanded, by the recombination of lithium ions Liqq ey™ Li q hn
Ž 1.
or by mutual neutralization Liqq Hy™ Li q H
Ž 2.
By radiative association 1
Present address: School of Chemistry, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK. 2 Corresponding author. E-mail: [email protected]
Liqq H ™ LiHqq hn
Ž 3.
the ionic molecule could then enter the cycle while
0301-0104r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 Ž 9 8 . 0 0 1 3 1 - 1
10
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the corresponding neutral could be produced by the exchange reaction LiHqq H ™ LiH q Hq Ž 4. or by either radiative association Li q H ™ LiH q hn Ž 5. or by associative detachment Li q Hy™ LiH q ey Ž 6. The neutral molecule thus formed could now undergo photodissociation LiH q hn ™ Li q H Ž 7. or charge transfer LiH q Hq™ LiHqq H Ž 8. q thus producing further LiH , while the exchange reaction LiH q H ™ Li q H 2 Ž 9. could additionally deplete the neutral molecule concentration. While several of the rates for the suggested processes of above have been obtained from theoretical calculations w5x, very little is known, either theoretically or experimentally, on the interaction forces or the reaction rates of Eq. Ž9. and they were simply ‘‘estimated’’ as a result of educated guesses w6x. This means that one of the important elements for the development of a global modelling of the early Universe and the protogalactic clouds is still fraught with uncertainties and that it is of interest to try to find out, as accurately as possible, both the features of the interaction guiding reaction Ž9. and the actual values of the reaction rates in the temperature regimes of interest. One of the recent studies on the LiH 2 system w7x has analyzed the vertical and adiabatic ionization potentials using second-order Møller-Plesset ŽMP2. calculations with triple-zeta plus polarization and diffuse basis functions. The aim was to study the y LiHy transition and the correspond2 ™ LiH 2 q e ing structures. It was found there that LiHy 2 has a stable linear D`h structure. The detachment of the extra electron leads in turn to the 2 Sq state of u neutral LiH 2 with an imaginary bending frequency value that prevents it from being a local minimum. The authors therefore found for the neutral molecule another optimal geometry of C2 Õ ,2 B2 symmetry that corresponds to a local minimum and a further C2 Õ Ž 2A1 . bound structure even lower than the previous one but still higher in energy than the Li q H 2
fragments. The latter geometry, therefore leads to . asymptotes without any barrier LiŽ 2 S. qH 2 Ž 1 Sq g and therefore it can dissociate. Since it also undergoes a nonadiabatic crossing with the C2 Õ Ž 2 B2 . state only 2.4 kcalrmol above the minimum of the latter, then one can say that also the 2 B2 geometry is not stable and it quickly dissociates into the same asymptotes via a nonadiabatic crossing with the C2 Õ Ž 2A1 . surface w7x. The above considerations are of importance in deciding on the possible existence of a stable transition state for this system and on its consequences in carrying out the reactive depletion of LiH molecule via process Ž9.. We have therefore completed a rather extensive calculation for the most important features of the lowest reactive surface of process Ž9. and have begun an exploratory analysis of its reactive behaviour in the special two-dimensional Ž2D. collinear arrangements involving Li-H-H and H-Li-H structures. The following Section discusses the details of the present computational method for the potential energy while Section 3 reports the most important features of the potential energy surface for the reaction ŽPESR.. Section 4 will in turn present and analyse the classical reactive calculations in the collinear arrangements while Section 5 will discuss the details of the wavepacket analysis. The general conclusions and a detailed comparison of methods will be presented in Section 6. 2. The valence-bond (VB) calculations The adiabatic LiH 2 potential energy surface has been calculated using the spin-coupled valence bond ŽSCVB. non-orthogonal configuration interaction ŽCI. method w8–11x. The expansion length of the CI wave function has been substantially reduced from that associated with the traditional implementation of the SCVB method by utilizing fully optimized nonorthogonal virtual orbitals. The orbital optimization was achieved using the highly efficient scheme developed recently by Clarke et al. w12x. The wave functions constructed from optimized virtual orbitals are denoted SCVB ) to differentiate them from the traditional SCVB wave functions, constructed using virtual orbitals produced by diagonalizing one-electron effective operators. The resulting SCVB ) wave
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
functions are even more compact than their SCVB counterparts, already in general smaller than energetically equivalent orthogonal MO-CI expansions w9,8,11,12x, but maintain the extremely high accuracy associated with the SCVB method. It is essential that effective methods of optimizing virtual orbitals are available to enable large, many-electron systems to be studied using sufficiently large basis sets, without the size of the Hamiltonian representation matrix becoming impossibly large. In the SCVB ) method, a set of non-orthogonal virtual orbitals is optimized within the framework of a multiconfiguration spin-coupled wave function. The reference function, from which the excited configurations are generated, is the single configuration, Nelectron spin-coupled wave function ŽSC. f SN
C0 s A
ž
f 10f 20 ... f N0
Ý cSk0 QSNM ; k ks1
/
Ž 10 .
where the f i0 4 are the N non-orthogonal occupied spin-coupled orbitals, f SN is the spin-degeneracy and the index k labels each of the f SN spin eigenfunctions Q SNM ; k of Sˆ2 and Sˆz . The f i0 4 are expanded as usual in a basis of M atomic orbitals M
f i0 s
Ý d i0a xa
Ž 11 .
as1
and the SC wave function is freely optimized with respect to both the orbital expansion coefficient, d i0a , 0 and the spin-coupling coefficients, c Sk . The f i0 4 are then divided, if appropriate, into a set containing the orbitals that describe the motion of the core electrons, and a set that describe the motion of the valence electrons. To each of the N Õ valence orbitals, f i0 4 Õ , we now associate a virtual orbital fq j which is orthogonal to its associated SC orbital, f j0 , but non-orthogonal to all other SC orbitals and virtual orbitals f k0 and fq k . From this expanded orbital set, we construct a multiconfiguration SC wave function consisting of the SC configuration, plus all possible configurations corresponding to a double replacement of SC valence orbitals with their associated virtual orbitals, N
C s c 0C 0 q
N
Ý Ý c i jCi j is1 j)i
Ž 12 .
11
where fSN
ž
q 0 Ci j s A f 10f 20 ... fq i ... f j ... f N
Ý cSki j QSNM ; k ks1
/
Ž 13 .
This results in a total of Ž 2N . excited configurations. The core set of SC orbitals, f i0 4 c, are not assigned virtual orbitals and so remain present in each of the excited configurations. As a result the core electrons are not correlated in this method beyond the level that already exists in the SC configuration. We have found that correlating them further leads to what is essentially just a uniform lowering of the potential surface at all points w12x. 4 The set of virtual orbitals fq i , is then optimized within the framework of this multiconfiguration wave function. To make the optimization as efficient as possible, we have used two simplifying approximations. Firstly, the excited configuration, spin-couij pling coefficients c Sk are fixed to the values for the 0 SC reference configuration c Sk and, secondly, the orbitals are optimized relative to the second-order perturbation energy w13,14x, not the variational energy. The major saving introduced by the second approximation is that it becomes necessary only to evaluate the diagonal and first row elements of the Hamiltonian and overlap matrices. The second-order perturbation energy expression is N Ž2.
E s H00 q
N
ÝÝ is1 j)i
HŽ0,i j. y H00 SŽ0 ,i j. H00 SŽ i j,i j. y HŽ i j,i j.
2
Ž 14 .
where H00 s ²C 0 < Hˆ
12
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
mate form for the Hessian, with the use of the second-order perturbation expression for the energy, results in an overall orbital optimization strategy that scales extremely favourably with both the number of active electrons and the number of basis functions. As a result, the SCVB ) method should now make the study of larger systems, not previously accessible to SCVB methods, possible. In the calculations for the current work the number of electrons, N, is five, the lowest energy state is a doublet, so f SN s 5. The lithium 1s 2 electrons were considered to be core since they have minimal influence on bond formation w16,17x, leaving a total of three valence spin-coupled orbitals and therefore three doubly-excited configurations. The results published here are a subset of a complete study of the full H-Li-H reactive surface. To maintain absolute consistency in our calculations, as the spatial arrangement of the three atoms changed symmetry, we have evaluated separately four discrete sets of virtual orbitals. Each set of virtual orbitals were expanded only in a subset of the atomic orbital ŽAO. basis relating to a certain class of functions. The four sets were created by dividing the full AO basis into one group containing just s functions and three groups each containing just one of the components of the p orbitals. The components of the Cartesian d functions were assigned to the appropriate group of the same symmetry. Each of the four sets of virtual orbitals were optimized separately, giving a total of twelve symmetry pure virtual orbitals. From these orbitals, it was possible to create a SCVB ) wave function which was equally flexible over the entire potential surface and did not change in character when the spatial arrangement of atoms changed symmetry. Without this technique, problems would have been experienced, particularly when the potential passes through the collinear geometries considered in this work. The separation of the virtual orbitals into separate symmetry pure groups had the added benefit of removing any linear dependence problems from which the SCVB ) method has a tendency to suffer. The final SCVB ) wave function included the DC reference configuration and all possible doubly-excited configurations which involved only ‘ vertical’ excitations of SC occupied orbitals into one of their four associated virtual orbitals. To this set, was
Table 1 Exponents a i and contraction K i for the GTO basis set Lithium
Hydrogen
Type
ki
ai
s
0.001367 0.010425 0.049859 0.160701 0.344604 0.425197 1.0 1.0 1.0 1.0 0.038770 0.236257 1.0 1.0 1.0 1.0
921.3 138.7 31.94 9.353 3.158 1.157 0.4446 0.07666 0.02864 0.00676 1.488 0.2667 0.07201 0.02370 0.3 0.1
s s s s p p p d d
Type
ki
ai
s s s s p p d
0.00255 0.01938 0.09280 1.0 1.0 1.0 1.0 1.0 1.0 1.0
68.160 10.2465 2.34648 0.67332 0.22466 0.082217 0.04 1.2 0.3 0.3
added all singly ionic configurations where a single virtual orbital was doubly occupied. This results in a total of just 73 configurations Ž293 VB structures. making it an extremely compact wave function. The a to m ic b a sis se t u se d c o n sid e re d o f Ž10s4p2dr7s2p1d. Cartesian GTOs contracted to Ž5s3p2dr5s2p1d. and is reported in Table 1. Ab initio points were calculated for the geometries RŽLi-H 1 . s 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.04 and RŽLi-H 2 . s RŽLi-H 1 . q 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 10.04 , where all the distances are in atomic units. In addition to these, a further set of points at RŽLi-H 1 . s 1.4, 2.1 and RŽLi-H 2 . s RŽLi-H 1 . q 0.7, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 10.04 were calculated to describe accurately the repulsive part of the potential. For all geometries, /ŽH 1 LiH 2 . s 0. These points were then interpolated by a two-dimensional spline function in order to be able to treat the two variables within the dynamical codes discussed in the following Sections, where the specific interpolated reactive surface will be employed. 3. Features of the collinear surface Fig. 1 reports in detail the general behaviour of the collinear arrangement. The plot on the left shows
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
13
Fig. 1. Computed reactive potential energy surface for the collinear reaction. Left: energy level diagrams for the two coordinates. Right: three-dimensional view of the reaction channels with the asymptotic profiles.
the potential energy contours associated to the relevant collinear coordinates, while the three dimensional presentation on the right shows clearly the presence of the large drop in energy on the side of
Fig. 2. Reaction coordinate profile as function of the two relevant internuclear distances along the minimum energy path.
the exothermic reaction Ž9.. The strongly repulsive nature of the interaction at short distances within the system Li-H-H is also seen in both presentations of the 2D surface and will be further discussed in the following sections. Another interesting feature of the 2D surface discussed here is presented in Fig. 2, where the minimum energy path is shown from the lowest electronic state of the Li q H 2 system Žleft. up to the energetically higher, lowest electronic state of the LiH q H system Žright.. One sees immediately that the LiH depletion reaction Ž9. is strongly exothermic with an energy gain on nearly 2.0 eV, while the opposite formation reaction Li q H 2 ™ LiH q H Ž 15 . is obviously strongly endothermic and requires highly excited H 2 molecular partners Žwith the vibrational index ÕX G 4. before it can occur. It is also of interest to notice the slight barrier which appears along the path for the depletion reaction when the H atom is close to the target molecule and the D R s R 2 y R 1 difference becomes negative. This feature will also be discussed in the following Sections. The five occupied orbitals from the single configuration SC calculation included a clearly defined lithium 1s 2 pair, which were overwhelmingly singlet coupled, at all points on the collinear potential sur-
14
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
face. These were easily identified as the orbitals containing the core electrons and there was never any uncertainty as to which SC orbitals should not be correlated further in the SCVB ) stage. The three
remaining occupied orbitals, free to delocalize over all of the atomic centers, changed nature according to the position on the potential surface. At some points they were heavily localized on a particular
Fig. 3. Comparison of the two asymptotic molecular partners as computed from the present VB correlated calculations. See text for definition of the acronyms.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
atomic center, for instance when there was one isolated hydrogen atom, whilst at others they were delocalized over more than one center, for instance when all three atoms were in close proximity. This is a very important feature of the spin-coupled method, since there is no external predetermination of the degree of localization experienced by the orbitals: it results as a natural consequence of the method. The spin-coupling of the electrons in the three valence occupied orbitals also varied with geometry, as the nature of the orbitals changed. At some points a perfect pairing scheme dominated, as two valence electrons became predominantly singlet coupled, whilst at others a more complicated situation was present and more than one mode of spin-coupling was active. This illustrates the importance, especially in open shell systems of this kind, of including the full spin space in the calculation and not just assuming that a perfect-pairing scheme will be a reasonable approximation. A good test of the quality of the calculated points is to consider the asymptotic regions of the potential that correspond effectively to the potential curves for the isolated diatomic species ŽLiH when RŽH 1-H 2 .. is very large and H 2 when RŽLi-H 1 . is very large..
15
Fig. 3 shows the asymptotic LiH potential curve calculated in this work at the single configuration SC and multiconfiguration SCVB ) levels, and compares it with the isolated LiH potential presented in our SCVB ) benchmark calculations w12x. The SC wave function recovers 76% of the true binding energy and overestimates the position of the minimum by ˚ The SCVB ) wave function in this work 0.04 A. recovers a greatly improved 93% of the binding energy, and overestimates the minimum by only ˚ which compares extremely well with the 0.015 A, ˚ respectively. benchmark values of 95% and 0.015 A, The SCVB ) binding energy is slightly lower in this work than for the benchmark as a smaller LiH basis set was used Ž43 contracted functions compared to 87 in the benchmark.. It is clear that including the extra configurations in going from the SC to the SCVB ) wave function has a considerable quantitative effect on the quality of the potential. Fig. 3 also shows that essentially the same story applies to the H 2 asymptotic potential curve. The SC and SCVB ) potentials from this work are shown compared to the position of the experimental minimum w18x. The SC wave function recovers 87% of the H 2 binding energy, whilst the SCVB ) wave
Fig. 4. Computed H-H distance dependence of the interaction using two different levels of correlated calculations as described in the main text.
16
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
function recovers 97%, and both reproduce the position of the minimum almost exactly. The results of Fig. 3 therefore show that the asymptotic regions of the two reaction channels, calculated at the SCVB ) level, are reproduced with very high degree of accuracy and thus indicate that the rest of the potential should also be of a very high quality, as our procedure in no way favours one region of the potential over another. In the non-asymptotic regions of the potential, it is interesting to note that one of the main effects of introducing the excited configurations was to lower the height of a barrier that blocked the path to a potential well as RŽH 1-H 2 . was reduced. This is illustrated in Fig. 4 which plots the energy as a function of RŽH 1-H 2 . with RŽLi-H 1 . fixed at 3 au, a geometry at which this effect is most pronounced. The asymptotic, large RŽH 1-H 2 ., energy has been rescaled to zero for both the SC and SCVB ) curves to facilitate the comparison of the potentials. In reality, the asymptotic part of the SCVB ) curve would sit well below that of the SC curve as there is still an interacting LiH system present. It is plainly seen that a barrier of about 10 mH Ž0.27 eV. at the SC level is reduced to one of just 2 mH Ž0.05 eV. at the SCVB ) level. This is likely to have a significant impact on the low energy dynamics for this reaction and emphasises the importance of including the additional dynamical correlation that the SCVB ) procedure provides.
4. The reactive scattering: a classical treatment We have carried out collinear classical dynamics for the depletion reaction H q LiHŽ Õ . ™ Li q H 2 as well as for the formation reaction Li q H 2 Ž Õ . ™ LiH q H reaction over a range of relative translational energies Ž Etrans . between 0–1 eV. The index Õ represents here the vibrational state of the diatom. For a given Õ and Etrans values the vibrational phase Ž f i . of the diatom was varied systematically in the range between 0–2p . Typically, 1000 trajectories were run to compute the reaction probability Ž P .. The initial separation of the reactants was taken to be in the asymptotic region at 10 a o . The Hamilton’s equations were solved in mass scaled and skewed coordinates w19x using the 4th order Runge-Kutta
method with a step size of 0.21548 fs. The accuracy of the integration was checked by employing the Gear integration along the evolution. In the present study, considering that the earlier calculation for a bent LiH 2 , species w7x suggested its inherent instability, we felt that it was sufficient to analyse the collinear reaction without the effects, which should therefore be minimal, of a bent geometry which dissociates quickly without forming an intermediate structure. 4.1. The depletion reaction: H q LiH(Õ) ™ Li q H2 Since the depletion reaction is exothermic and has practically no barrier, it follows that this reaction goes downhill over the whole reactive path. However, there appears to be a rather small barrier of 0.036 eV in the interaction zone Žas can be seen in Fig. 2.. The reaction is exothermic by 2.0238 eV considering the zero-point energies of the diatoms and therefore one expects to see only a small influence, and then only near threshold, of the above barrier. The calculated probability Ž P depletion . for the reaction is shown in Fig. 5 as a function of Etrans for different Õ-states of LiH. The corresponding probabilities obtained from the collinear time-dependent wave packet ŽTDWP. calculations, which will be described in Section 5, are also displayed in the figure. The classical results are marked with a solid line connecting the circles. The representative statistical error bars are also shown at Etrans s 1.0 eV for all Õ-states. We note the following: Ža. For all Õ states P depletion increases sharply at the low energy, reaching its maximum around Etrans s 0.1 eV. It then drops down to its minimum around 0.2 eV. Žb. With further increase in Etrans the corresponding reaction probability shows a slow monotonic increase, which however remains smaller for higher Õ states. Žc. Below 0.1 eV the classical probability shows some oscillations, which are also well marked when one looks at the wave packet calculations in the same energy region. To analyse the unusual energy dependence of P depletion between Etrans s 0.1–0.2 eV, we have plotted in Fig. 6 a family of Žten. trajectories given in the usual mass-scaled and skewed coordinates Ž x, y . w19x and which also sample uniformly f i of LiH between
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
0 and 2p for different Etrans . Here x represents the relative motion of the atom with respect to the diatom and y denotes the internal motion of the diatom. At Etrans s 0.1 eV, most of the trajectories are seen to end into the reactive channel. At the other energies examined, however, only a few of them find a way to go into the reactive channel, while the majority of them bounces back into the entrance channel after dropping downhill and hitting the repulsive core of the interaction region. Despite the fact that a trajectory has to travel all the way down the reaction path Žsee Fig. 2. for a reaction to occur, the dynamical constraints force it to get back into the entrance channel whenever an excess of translational energy is available. In order to get into the exit channel, therefore, they have to reach the interaction
17
zone with the right phase which then would allow them go into the exit channel, a condition which appears to be achieved only at smaller Etrans . Increasing Etrans does not help to go into the reactive region as a result the nonreactive trajectories which bounce off the repulsive wall at a narrow angle and come back to yield an increasingly higher vibrational excitation content of LiH, as seen, for example, for Etrans s 0.9 eV. To further analyse this special behaviour we have plotted in Fig. 7 the final classical vibrational action Ž n f . of the diatoms as a function of f i of LiH for both reactive and nonreactive cases systematically sampling 100 trajectories between 0–2p for a range of Etrans values. The left panel displays the Etrans values around 0.1 eV where the P depletion shows in
Fig. 5. Computed depletion reaction probabilities as function of collision energy, Etrans , and of the vibrational content of the LiH partner. The classical results are given by the open circles connected by solid lines. The quantum wavepacket results are given by the dotted lines.
18
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
Fig. 5 its first maximum and then a dip, while the right panel shows the results for Etrans ) 0.1 eV. The crosses and the circles indicate the nonreactive ŽNR. and reactive ŽR. trajectories, repectively. We see that, in general, the reactivity is a banded function of f i , that is, within certain range of f i values the reactive events occur while the NR events dominate within another range. The width of the reactive band starts decreasing with the increase of Etrans above 0.1
eV. Even when only a small increase occurs, the R phases dramatically shrink around 0.13 eV and remain roughly the same up to 0.2 eV. With the further increase of Etrans its range then starts to slowly increase. This explains in some way the slow rise in the reaction probability for Etrans ) 0.2 eV. We note that the width of the R band is quite sensitive for collisions at low energies. At Etrans s 0.13, in fact, we observe a marked fractured behaviour between
Fig. 6. Classical calculations of reactive and non reactive trajectories, as function of translational energy, for the depletion reaction. See text for the meaning of symbols.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27 Fig. 7. Computed classical vibrational action Ž n f . of the LiH diatom as function of the phase f i and for different values of translational energy. The circles describe the reactive trajectories while the crosses refer to the nonreactive trajectories.
19
20
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the switch over regions of the R and NR bands. Such switch over regions are known w20x to have fractal characteristics and correspond to the onset of chaotic behaviour for the trajectories. What it means is that they show the existence of a ‘‘smooth’’ banded and a switchover ‘‘chattering’’ region buried within a finer resolution of f i . The trajectories in the chattering regions show multiple collisions in the interaction zone before finally going into either the entrance or the exit channels. We believe that the small oscillations in the classical probability arise presumably because of the existence of nonlinearity in the dynamics at low Etrans , which however becomes less and less so as the collision energy is increased. It is
interesting to note, as we shall further discuss below that the quantum mechanical ŽQM. results also show oscillations at low energies. Finally, it is worth pointing out that in the present 2D reactive process for the depletion reaction the product H 2 molecules appear to be formed in vibrationally excited states. 4.2. The formation reaction: Li q H2 (Õ) ™ LiH q H The formation reaction is an endothermic process which requires an additional energy of 2.06 eV. In order for the reaction to occur, a given trajectory has to then go up along the uphill path Žsee Figs. 1 and 2. and the barrier could be regarded as located in the
Fig. 8. Computed reactive probabilities for the endothermic formation reaction: Li q H 2 ™ LiH q H as function of H 2 vibrational content and of relative translational energy. The classical calculations are given by open circles connected by solid lines while the quantum calculations are given by dotted lines.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
exit channel, i.e., to be of the late uphill type. This suggests that by putting energy into the vibrational mode of H 2 one would more efficiently help the formation reaction than when putting the equivalent amount of energy into the relative translational energy, Etrans . Since we are interested in the low-energy collisions beween 0–1.0 eV, it turns out that H 2 should be at least vibrationally excited in its Õ s 4 vibrational level for the reaction to occur in the above energy range. The computed classical probabilities Ž P formation . for the formation reaction are shown as a function of Etrans , in Fig. 8 for Õ s 4,5,6 and 7. The corresponding probabilities obtained from the quantum mechanical computations are also displayed in that figure. As before, the classical results are marked with a solid line connecting the circles. For Õ s 4 the classical probability rises from the energetic threshold, becomes flat up to 0.7 eV, and then rises with a further increase of Etrans . For Õ s 5,6 and 7 the probability rises sharply at low collision energies, follows a marked decrease and then increases slowly and monotonically with the further rise of Etrans . When comparing the classical results with the QM probabilities, one sees that the agreement is in general quite good, particularly for Õ s 5–7 and the predicted probabilites agree quite well with each other both at low and at high collision energies. For Õ s 4, although the general trend is maintained in the classical results, there seems to be some quantitative discrepancies at the higher values of Etrans . It is interesting to note that the probability of formation is larger at low collision energies for Õ s 5–7. However, for Õ s 4 it becomes larger when the collision energy is increased. Understandably, this effect arises because of the fact that the energy content for Õ G 5 is already nearly equal or above the reaction threshold, and little Etrans can already help a given trajectory to reach the interaction zone with the right ‘‘phase’’ for the reaction to occur efficiently, while it does not seem to be effective when the collision energy is further increased. In that case the trajectories get back into the reactant channel after hitting the repulsive core of the interaction zone as previously discussed. For Õ s 4, an increase in Etrans from the energetic threshold helps the occurence of the reaction since part of the available energy is used up to further climb over the energy barrier.
21
5. The reactive scattering: a quantal treatment 5.1. Computational method The Time Dependent Wave Packet ŽTDWP. method has been widely used in the recent years for reactive scattering processes w21x. Such method is based on the discretisation of the wavefunction on a finite grid, which is propagated in time according to the time-dependent Schrodinger equation. Since the ¨ method is well known, in this section we will only describe the details of the present implementation relevant for the subsequent discussion. In the following the hamiltonian of the system A q BC in its collinear configuration is represented using the mass scaled and skewed coordinates Ž x, y . introduced before w21x and it is given by "2 Hsy
2m
ž
E2 E x2
E2 q
E y2
/
qV Ž x, y.
Ž 16 .
with m s Ž mA m B m C rmA q m B q m C .Ž1r2.. The initial wavefunction C Ž x, y,t s 0. is chosen as a product of a minimum uncertainty Gaussian wavepacket ŽGWP., F Ž x ., which describes the translational part of the wavefunction while the vibration of the diatomic molecule is represented by a Morse oscillator eigenfunction, x Õ Ž y ., corresponding to the vibrational state Õ of the diatomic fragment with energy eÕ
C Ž x , y,t s 0 . s F Ž x . x Õ Ž y .
Ž 17 .
where FŽ x. s
1
Ž 2pd . 1r4 Ž
.
ž
exp y
Ž x y x0 . 4d 2
2
y ik 0 x
/ Ž 18 .
Here x 0 is the center of the wavepacket, placed in the reactant region, while k 0 is the momentum determined from the translational energy of the wavepacket, k 0 s w2 m Ž Etot y e Õ . r" x where Etot is the total energy, corresponding to the center of the wavepacket in the momentum space. The width parameter d controls the spread of the wavepacket in both the coordinate and momentum space: optimal values of the parameter are those which lead to a wavepacket sufficiently narrow in the coordinate
(
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
22
space to be included in the asymptotic reactant region as well as a distribution in the momentum space sufficiently large to allow for the determination of the reaction probabilities for a wide range of translational energies. The time evolution of the wavepacket is carried out using the well-known Split Operator method w22x, where the evolution operator associated with a time independent hamiltonian is approximated as
ž
exp y
iHD t "
/
ž
iTD t
/ ž ž /
f exp y
exp y
2"
=exp y
iVD t "
/
iTD t 2"
Ž 19 .
where T and V are the kinetic and potential energy operator, respectively, and D t is the time step used in the propagation. The action of the kinetic energy operator is evaluated in its local representation by using the Fast Fourier Transform ŽFFT. algorithm. The reaction probability is evaluated using the time-energy mapping of the reactive flux w23x: by fixing a dividing line at y s yI in the product region, the reaction probability is then given by PÕ s
ž
"2 k
mm x < F Ž k . < 2
ž¦
Im C Ž x , yI , E .
/
EC Ž x , yI , E .
;/
Ey
Ž 20 .
where m x s w mAŽ m B q m C .rŽ mA q m B q m C .x, k s 2 m x Ž E y e Õ . r", F Ž k . is the representation in the momentum space of F Ž x . and C Ž x, yI , E . is the Fourier transform of the wavepacket, C Ž x, yI ,t ., calculated along the dividing line y s yI
(
C Ž x , yI , E . 1
s
ž / 2p
1r2
q`
Hy` C Ž x , y ,t . exp Ž iEtr" . d t I
Ž 21 . After the flux analysis, the wavepacket is absorbed by a negative imaginary potential ŽNIP. placed near the edge of the grid in both reactant and product channels. A linear type of NIP, as suggested by
Neuhauser and Baer w24x is used in the present calculation
Vabs
°yiV s~ ¢0
0
R y R abs D R abs
R abs - R F R abs q D R abs R F R abs
Ž 22 . where V0 is the height of the potential, D R abs is its width and R abs is its location, which usually takes values far out from the interaction region. The time propagation is carried out until a sufficiently large fraction of the wavepacket has reached the boundaries of the grid, where it is absorbed by the NIP. 5.2. Details of the calculation In the present work, we have calculated the reaction probabilitites of the Ži. depletion and Žii. formation reactions of LiH in collision with H, in the collinear configuration, as a function of the energy, for different initial vibrational states of the target molecule. Using the formalism introduced above, those reactions correspond to the following way of labelling the three atoms: Ž1. A s Li, B s C s H and Ž2. A s H, B s Li, C s H. In both calculations, we have used a grid in Ž x, y . space composed by Ž128 = 128. points; due to the different masses of the target molecules, on which the scaled and skewed coordinates depend, the grid parameters are different for the depletion and formation reaction. In the former case, the grid starts from ˚ y s 0.74 A˚ with step sizes D x s 0.060 x s 1.02 A, ˚A and D y s 0.035A, ˚ while in the latter case the ˚ y s 0.28 A˚ with step reaction starts at x s 1.17 A, ˚ ˚ sizes of D x s 0.047 A and D y s 0.055 A. The calculations have been carried out for a range of translational energies between 0.2 eV and 1.1 eV, with the initial vibrational level of LiH ranging from Õ s 0 to Õ s 4 in the depletion reaction, while, due to the topography of the surface as discussed before, we have included only excited vibrational levels of H 2 in the formation reaction, from Õ s 4 to Õ s 7. The initial wavepackets have been centered at x 0 f 5 ˚ using different widths in the two reactions, to take A, into account the different distributions of the transla-
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
tional energy: in the depletion reaction, d has been ˚ while we have used d s 0.13 in the fixed at 0.2 A formation reaction. As a consequence of this choice, we have been able to calculate reaction probabilities for an energy range of 0.3 eV with a single propagation for the two reactions under study.
23
The dividing lines for the calculation of the flux ˚ for the depletion have been fixed at yI s 4.0 A reaction and at yI s 4.62 for the formation reaction. Together with the main dividing line, other two lines have been placed 10 steps before and after the main dividing line in order to check, from the values of
Fig. 9. Snapshots of the propagation of a wavepacket for the depletion reaction of LiH, with LiH in its initial vibrational level Õ s 1 and an initial mean translational energy of about 0.3 eV.
24
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the reaction probabilities, that the position of yI has been placed well inside the product asymptotic region. The location of the absorbing potentials in the product channels have been chosen to be 15 steps from the dividing line used to calculate the flux; in
˚ for the reactant channels, we have used x abs s 6.0 A ˚ for the the depletion reaction and x abs s 5.32 A reactant channels with a height of about 1 eV in both cases. Finally, a time step of 0.04 fs was chosen for the depletion reaction, reduced to 0.02 fs in the forma-
Fig. 10. Snapshots of the propagation of a wavepacket for the reaction H 2 q Li, with H 2 in its Õ s 4 level and an initial translational energy of about 0.9 eV.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
tion reaction: the latter value has been taken to be smaller than the former to take into account the larger energy content of the wavepackets used in the formation reaction. It is worth mentioning that such numerical parameters have been the results of a rather extensive number of test calculations for numerical convergence. The propagation has been carried out until the 99.9% of the wavepacket has been absorbed; such large fraction ensures the convergence of the reaction probabilities on the number of chosen time steps. We have found convenient to interface the program of the time dependent propagation with a graphical package that allows us to watch the propagation of the bidimensional wavepacket during the execution of the program. In order to have some insight on the physical processes under investigation, we have reported in Fig. 9 a series of snapshots of the propagation of a specific wavepacket describing the depletion reaction, with initial vibrational state Õ s 1 of LiH and a small initial translational energy, at equally spaced instants of time. From the reaction probabilities reported in Fig. 5 such wavepacket contains energy components placed around the minimum of the curve. We can see that, due to the nature of the potential, the wavepacket is not captured in the interaction region therefore the reaction appears to be a direct process. Furthermore, we can observe from the shape of the wavepacket in the product channel that H 2 is formed with quite a large vibrational energy content. At any rate, we can also see that, after the collisional event, the wavepacket in the product region has a different shape from that at t s 0. It is very probable that all the available vibrational levels of LiH are populated after the collision so that there is a reduction of the average translational energy and the wavepacket moves now moves slowly towards the asymptote of the reactant region. In such a case a very large number of time steps is needed to reach the end of the propagation. The convergence of the reaction probabilities in the cases where the total energy is close to the threshold of the reaction needs to be carefully examined: to this aim, we have carried out a test calculation using a grid of Ž256 = 256. points, for a wavepacket with the same initial translational energy as the one reported in Fig. 9 and with the lowest vibrational state, Õ s 0, in order to have the minimum total energy: we have found that
25
the reaction probabilities are in good agreement with those obtained with the the less dense grid. We have also reported further snapshots for the formation reaction in Fig. 10, for a wavepacket with H 2 in the vibrational state Õ s 4 and average translational energy of about 0.9 eV. Even in this case the reaction is direct; furthermore, we can observe that the reactive flux is due to the outer part Žlarge y . of the initial wavepacket, while the inner part remains almost unaltered by the collision, the opposite of what we have found in the reverse, depletion reaction discussed above. It is also interesting to note that the corresponding depletion and formation probabilities reported in Figs. 5 and 8 appear to be in very good agreement with the classical results. In other words, we see from the depletion data of Fig. 5 that the quantum calculations also suggest a direct reaction where the excess energy is not always driven into the translational component: little additional relative energy content, in fact, is able to make the reaction to occur with large probabilities no matter what the vibrational level of the LiH is initially chosen to be. As the collision energy is increased, on the other hand, the efficiency of the reaction decreases rather markedly and only picks up rather slowly especially for initially excited, vibrationally ’hot’ partner molecules. The situation is shown to be similar for the formation reaction of Fig. 8 where, however, we see larger discrepancies between quantum and classical calculations at low energies possibly indicating the more significant role of vibrational coupling, in the quantum case, between Etrans and H 2 energy content.
6. Concluding remarks In the present work we have approached, for the first time, the study of the depletion and formation reactions of astrophysical interest which involve LiH molecules. Starting from first principles we therefore have computed over a very extensive grid of points a highly correlated reactive surface for the two-dimensional collinear processes relative to the above two reactions. After accurate numerical interpolation of the initial points, we have carried out classical trajectory calculations for the reactive case and also time-de-
26
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
pendent, wavepacket propagation studies for the quantum reactive process. The results from the two approaches, and a comparison of their computed probabilities, shows the following: Ži. The depletion reaction is strongly exothermic and appears to be, in the 2D case, a direct reaction without any transition state region. Žii. The shape of the reactive potential energy surface shows that low translational energy content markedly favours the depletion reaction and that its increase leads initially to very small reaction probabilities as the reagents ‘bounce back’ into the entrance channel due to the special repulsive features of the skewed coordinate system for the partner masses. Žiii. The initial vibrational content of LiH in the depletion reaction plays a rather minor role Žsee Fig. 5. in the sense that the final reaction probability depends very little on it. Živ. Both classical and quantum treatments at low and intermediate energies agree remarkably well with each other and show the time-dependent approach to be only adding a more marked interference structure in the very-low energy regions. Žv. The formation reaction requires highly excited H 2 molecules for it to occur at low collision energies. Both classical and quantum results indicate, in fact, that the H 2 vibrational levels Õ s 4,5,6,7 provide very efficient starting partners for the endothermic reaction to take place. Žvi. In the formation reaction the calculations show once more the rather inefficient coupling of vibrational and translational forms of energy already observed before for the depletion reaction: low Etrans values markedly favour reaction even when excess energy would be available from Evib . On the whole, this preliminary study gives us some useful indications on what we should expect from a more extensive calculation of the full 3D reaction probability: the quantum effects are not seen to be important for the depletion process, a fairly direct reactive event without complex formation. Furthermore, we see that internal energy content of LiH is not an helping factor for the process to occur and therefore one should expect that ‘cold’ LiH molecules can be depleted by reactive encounters with H atom at fairly low values of translational energy. With the same token, the formation of LiH
from vibrationally excited H 2 molecules appears to lead to the formation of vibrationally excited LiH, as seen from the time-dependent snapshots described before. Clearly one further needs to analyse the behaviour of the full 3D potential for the reaction rates and therefore the rather delicate process of fitting the available array of Born-Oppenheimer points is currently under way in our group. The preliminary collinear results, on the other hand, are already able to tell us a great deal on the specific features of the reaction and on its likely success to be realistically described by classical calculations. As a matter of fact, and as mentioned earlier, the computed instabilities of the bent structures w7x and the absence of complex formation in the 2D study conjure up a physical picture where head-on collisions greatly favour reactions, and therefore one should not expect to find that 3D calculations strongly ‘push away’ flux from the collinear path.
Acknowledgements The financial support of the Italian National Research Council ŽCNR. and of the Italian Ministry for University and Research ŽMURST. is gratefully acknowledged. N.J.C. thanks the European union for a training grant under the TMR project. S.K. wishes to thank the Max-Planck-Gesellshaft for supporting his stay at the University of Rome through the funding of a Max-Planck-Research-Prize to F.A.G. and E.B. also wishes to thank Dr. Pietro Lanucara for his generous help with graphic interfaces.
References w1x V.K. Dubrovich, Astron.&Astrophys. Trans. 5 Ž1994. 57. w2x R. Maoli, F. Melchiorri, D. Tosti, Astrophys. J. 425 Ž1994. 372. w3x S. Lepp, J.M. Shull, Astrophys. J. 280 Ž1984. 465. w4x V.K. Dubrovich, A.A. Lipovka, Astron. &Astrophys. 296 Ž1995. 301. w5x P.C. Stancil, S. Lepp, A. Dalgarno, Astrophys. J. 548 Ž1996. 401. w6x E. Bougleux, D. Galli, private comm. w7x A.I. Boldyrev, J. Simons, J. Chem. Phys. 99 Ž1993. 4628. w8x J. Gerrat, M. Raimondi, Proc. Poy. Soc. London A 371 Ž1980. 525.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27 w9x D.L. Cooper, J. Gerrat, M. Raimondi, Int. Rev. Phys. Chem. 7 Ž1985. 59. w10x J. Gerrat, D.L. Cooper, M. Raimondi, in: D.J. Klein, N. ´ ŽEds.., Elsevier, Amsterdam, 1990, p. 287. Trinastic w11x D.L. Cooper, J. Gerrat, M. Raimondi, Chem. Rev. 91 Ž1991. 929. w12x N.J. Clarke, M. Raimondi, M. Sironi, J. Gerratt, D.L. Cooper, Theor. Chim. Acta 95 Ž1997. 1328. w13x P.O. Lowdin, J. Chem. Phys. 19 Ž1951. 1396. ¨ w14x R. McWeeny, Methods of Molecular Quantum Mechanics, Academic Press, London, 2nd edition, 1992. w15x S.M. Goldfield, R.E. Quandt, H.F. Trotter, Research Memorandum No. 95, Economic Research Program, Princeton University, 1968. w16x N.C. Pyper, J. Gerratt, Proc. Roy. Soc. Lond. A 355 Ž1977. 406.
27
w17x H. Partridge, C.W.B.Jr., P.E.M. Siegbahn, Chem. Phys. Lett. 97 Ž1983. 198. w18x G. Herzberg, J. Mol. Spectry 33 Ž1970. 147. w19x P.J. Kuntz, in: W.H. Miller ŽEd.., Dynamics of Molecular Collisions Plenum, 1976. w20x V. Balasubramanian, B.K. Mishea, A. Bahel, S. Kumar, N. Sathyamurthy, J. Chem. Phys. 95 Ž1991. 4160. w21x N. Balakrishnan, C. Kalyanaraman, N. Sathyamurthy, Phys. Rep. 280 Ž1976. 129. w22x J.A. Fleck Jr., M.D. Feit, J. Chem. Phys. 78 Ž1982. 301; M.D. Feit, J.A. Fleck Jr., A. Steiger, J. Comp. Phys. 47 Ž1982. 412; D. Kosloff, R. Kosloff, J. Comp. Phys. 52 Ž1983. 35. w23x N. Balakrishnan, N. Sathyamurthy, Chem. Phys. Lett. 240 Ž1995. 119. w24x D. Neuhauser, M. Baer, J. Chem. Phys. 90 Ž1989. 4351.
Classical and quantum dynamics on the collinear potential energy surface for the reaction of Li with H 2 Nick J. Clarke a,1, Maurizio Sironi a , Mario Raimondi a , Sanjay Kumar b, Franco A. Gianturco b,2 , Erasmo Buonomo c , David L. Cooper d a
Dipartimento di Chimica Fisica ed Elettrochimica, UniÕersita´ di Milano, e Centro CNR-CSRSRC, Via Golgi 19, 20133 Milano, Italy b Department of Chemistry, The UniÕersity of Rome, Citta´ UniÕersitaria, 00185 Rome, Italy c Center for Supercomputing Applications to Research (CASPUR), Citta´ UniÕersitaria, 00185 Rome, Italy d Department of Chemistry, UniÕersity of LiÕerpool, P.O. Box 147, LiÕerpool L69 7ZD, UK Received 22 December 1997
Abstract Extensive Valence Bond ŽVB. calculations using non-orthogonal, spin-coupled wavefunctions with optimized orbitals have been employed to analyse the reactive region for the title process. Both the exothermic channel and the reverse, endothermic process are of strong astrophysical interest although no previous calculations have been available on both the reactive dynamics and the interaction energy surface. The specific features of the potential are analysed for some indicative configurations and classical trajectory calculations are carried out for the special collinear arrangement. In the latter instance, quantum time-dependent wavepacket calculations have also been performed and the two sets of results are found to be in rather good accord with each other. Some consequences of these exploratory calculations are discussed in detail. q 1998 Elsevier Science B.V. All rights reserved.
1. Introduction In the general scenario of what is known as ‘‘Early Universe’’ chemistry, the formation and depletion of LiH and LiHq is known to play a rather important role w1–3x since such species, together with H, Hq, H 2 , Hq 2 and others, can induce in the primordial gas spectral distortions and spatial anisotropies of the Cosmic Background Radiation ŽCBR.. Thus, the possible presence of molecules in such region can act as a heat pump, absorbing energy
in their rotovibrational levels and re-emitting this energy at lower frequencies by undergoing rotational cooling which will thus cause deviations of the CBR spectrum from the pure Planck spectrum w4x. It is usually assumed that the lithium chemistry was initiated, as the universe expanded, by the recombination of lithium ions Liqq ey™ Li q hn
Ž 1.
or by mutual neutralization Liqq Hy™ Li q H
Ž 2.
By radiative association 1
Present address: School of Chemistry, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK. 2 Corresponding author. E-mail: [email protected]
Liqq H ™ LiHqq hn
Ž 3.
the ionic molecule could then enter the cycle while
0301-0104r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 Ž 9 8 . 0 0 1 3 1 - 1
10
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the corresponding neutral could be produced by the exchange reaction LiHqq H ™ LiH q Hq Ž 4. or by either radiative association Li q H ™ LiH q hn Ž 5. or by associative detachment Li q Hy™ LiH q ey Ž 6. The neutral molecule thus formed could now undergo photodissociation LiH q hn ™ Li q H Ž 7. or charge transfer LiH q Hq™ LiHqq H Ž 8. q thus producing further LiH , while the exchange reaction LiH q H ™ Li q H 2 Ž 9. could additionally deplete the neutral molecule concentration. While several of the rates for the suggested processes of above have been obtained from theoretical calculations w5x, very little is known, either theoretically or experimentally, on the interaction forces or the reaction rates of Eq. Ž9. and they were simply ‘‘estimated’’ as a result of educated guesses w6x. This means that one of the important elements for the development of a global modelling of the early Universe and the protogalactic clouds is still fraught with uncertainties and that it is of interest to try to find out, as accurately as possible, both the features of the interaction guiding reaction Ž9. and the actual values of the reaction rates in the temperature regimes of interest. One of the recent studies on the LiH 2 system w7x has analyzed the vertical and adiabatic ionization potentials using second-order Møller-Plesset ŽMP2. calculations with triple-zeta plus polarization and diffuse basis functions. The aim was to study the y LiHy transition and the correspond2 ™ LiH 2 q e ing structures. It was found there that LiHy 2 has a stable linear D`h structure. The detachment of the extra electron leads in turn to the 2 Sq state of u neutral LiH 2 with an imaginary bending frequency value that prevents it from being a local minimum. The authors therefore found for the neutral molecule another optimal geometry of C2 Õ ,2 B2 symmetry that corresponds to a local minimum and a further C2 Õ Ž 2A1 . bound structure even lower than the previous one but still higher in energy than the Li q H 2
fragments. The latter geometry, therefore leads to . asymptotes without any barrier LiŽ 2 S. qH 2 Ž 1 Sq g and therefore it can dissociate. Since it also undergoes a nonadiabatic crossing with the C2 Õ Ž 2 B2 . state only 2.4 kcalrmol above the minimum of the latter, then one can say that also the 2 B2 geometry is not stable and it quickly dissociates into the same asymptotes via a nonadiabatic crossing with the C2 Õ Ž 2A1 . surface w7x. The above considerations are of importance in deciding on the possible existence of a stable transition state for this system and on its consequences in carrying out the reactive depletion of LiH molecule via process Ž9.. We have therefore completed a rather extensive calculation for the most important features of the lowest reactive surface of process Ž9. and have begun an exploratory analysis of its reactive behaviour in the special two-dimensional Ž2D. collinear arrangements involving Li-H-H and H-Li-H structures. The following Section discusses the details of the present computational method for the potential energy while Section 3 reports the most important features of the potential energy surface for the reaction ŽPESR.. Section 4 will in turn present and analyse the classical reactive calculations in the collinear arrangements while Section 5 will discuss the details of the wavepacket analysis. The general conclusions and a detailed comparison of methods will be presented in Section 6. 2. The valence-bond (VB) calculations The adiabatic LiH 2 potential energy surface has been calculated using the spin-coupled valence bond ŽSCVB. non-orthogonal configuration interaction ŽCI. method w8–11x. The expansion length of the CI wave function has been substantially reduced from that associated with the traditional implementation of the SCVB method by utilizing fully optimized nonorthogonal virtual orbitals. The orbital optimization was achieved using the highly efficient scheme developed recently by Clarke et al. w12x. The wave functions constructed from optimized virtual orbitals are denoted SCVB ) to differentiate them from the traditional SCVB wave functions, constructed using virtual orbitals produced by diagonalizing one-electron effective operators. The resulting SCVB ) wave
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
functions are even more compact than their SCVB counterparts, already in general smaller than energetically equivalent orthogonal MO-CI expansions w9,8,11,12x, but maintain the extremely high accuracy associated with the SCVB method. It is essential that effective methods of optimizing virtual orbitals are available to enable large, many-electron systems to be studied using sufficiently large basis sets, without the size of the Hamiltonian representation matrix becoming impossibly large. In the SCVB ) method, a set of non-orthogonal virtual orbitals is optimized within the framework of a multiconfiguration spin-coupled wave function. The reference function, from which the excited configurations are generated, is the single configuration, Nelectron spin-coupled wave function ŽSC. f SN
C0 s A
ž
f 10f 20 ... f N0
Ý cSk0 QSNM ; k ks1
/
Ž 10 .
where the f i0 4 are the N non-orthogonal occupied spin-coupled orbitals, f SN is the spin-degeneracy and the index k labels each of the f SN spin eigenfunctions Q SNM ; k of Sˆ2 and Sˆz . The f i0 4 are expanded as usual in a basis of M atomic orbitals M
f i0 s
Ý d i0a xa
Ž 11 .
as1
and the SC wave function is freely optimized with respect to both the orbital expansion coefficient, d i0a , 0 and the spin-coupling coefficients, c Sk . The f i0 4 are then divided, if appropriate, into a set containing the orbitals that describe the motion of the core electrons, and a set that describe the motion of the valence electrons. To each of the N Õ valence orbitals, f i0 4 Õ , we now associate a virtual orbital fq j which is orthogonal to its associated SC orbital, f j0 , but non-orthogonal to all other SC orbitals and virtual orbitals f k0 and fq k . From this expanded orbital set, we construct a multiconfiguration SC wave function consisting of the SC configuration, plus all possible configurations corresponding to a double replacement of SC valence orbitals with their associated virtual orbitals, N
C s c 0C 0 q
N
Ý Ý c i jCi j is1 j)i
Ž 12 .
11
where fSN
ž
q 0 Ci j s A f 10f 20 ... fq i ... f j ... f N
Ý cSki j QSNM ; k ks1
/
Ž 13 .
This results in a total of Ž 2N . excited configurations. The core set of SC orbitals, f i0 4 c, are not assigned virtual orbitals and so remain present in each of the excited configurations. As a result the core electrons are not correlated in this method beyond the level that already exists in the SC configuration. We have found that correlating them further leads to what is essentially just a uniform lowering of the potential surface at all points w12x. 4 The set of virtual orbitals fq i , is then optimized within the framework of this multiconfiguration wave function. To make the optimization as efficient as possible, we have used two simplifying approximations. Firstly, the excited configuration, spin-couij pling coefficients c Sk are fixed to the values for the 0 SC reference configuration c Sk and, secondly, the orbitals are optimized relative to the second-order perturbation energy w13,14x, not the variational energy. The major saving introduced by the second approximation is that it becomes necessary only to evaluate the diagonal and first row elements of the Hamiltonian and overlap matrices. The second-order perturbation energy expression is N Ž2.
E s H00 q
N
ÝÝ is1 j)i
HŽ0,i j. y H00 SŽ0 ,i j. H00 SŽ i j,i j. y HŽ i j,i j.
2
Ž 14 .
where H00 s ²C 0 < Hˆ
12
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
mate form for the Hessian, with the use of the second-order perturbation expression for the energy, results in an overall orbital optimization strategy that scales extremely favourably with both the number of active electrons and the number of basis functions. As a result, the SCVB ) method should now make the study of larger systems, not previously accessible to SCVB methods, possible. In the calculations for the current work the number of electrons, N, is five, the lowest energy state is a doublet, so f SN s 5. The lithium 1s 2 electrons were considered to be core since they have minimal influence on bond formation w16,17x, leaving a total of three valence spin-coupled orbitals and therefore three doubly-excited configurations. The results published here are a subset of a complete study of the full H-Li-H reactive surface. To maintain absolute consistency in our calculations, as the spatial arrangement of the three atoms changed symmetry, we have evaluated separately four discrete sets of virtual orbitals. Each set of virtual orbitals were expanded only in a subset of the atomic orbital ŽAO. basis relating to a certain class of functions. The four sets were created by dividing the full AO basis into one group containing just s functions and three groups each containing just one of the components of the p orbitals. The components of the Cartesian d functions were assigned to the appropriate group of the same symmetry. Each of the four sets of virtual orbitals were optimized separately, giving a total of twelve symmetry pure virtual orbitals. From these orbitals, it was possible to create a SCVB ) wave function which was equally flexible over the entire potential surface and did not change in character when the spatial arrangement of atoms changed symmetry. Without this technique, problems would have been experienced, particularly when the potential passes through the collinear geometries considered in this work. The separation of the virtual orbitals into separate symmetry pure groups had the added benefit of removing any linear dependence problems from which the SCVB ) method has a tendency to suffer. The final SCVB ) wave function included the DC reference configuration and all possible doubly-excited configurations which involved only ‘ vertical’ excitations of SC occupied orbitals into one of their four associated virtual orbitals. To this set, was
Table 1 Exponents a i and contraction K i for the GTO basis set Lithium
Hydrogen
Type
ki
ai
s
0.001367 0.010425 0.049859 0.160701 0.344604 0.425197 1.0 1.0 1.0 1.0 0.038770 0.236257 1.0 1.0 1.0 1.0
921.3 138.7 31.94 9.353 3.158 1.157 0.4446 0.07666 0.02864 0.00676 1.488 0.2667 0.07201 0.02370 0.3 0.1
s s s s p p p d d
Type
ki
ai
s s s s p p d
0.00255 0.01938 0.09280 1.0 1.0 1.0 1.0 1.0 1.0 1.0
68.160 10.2465 2.34648 0.67332 0.22466 0.082217 0.04 1.2 0.3 0.3
added all singly ionic configurations where a single virtual orbital was doubly occupied. This results in a total of just 73 configurations Ž293 VB structures. making it an extremely compact wave function. The a to m ic b a sis se t u se d c o n sid e re d o f Ž10s4p2dr7s2p1d. Cartesian GTOs contracted to Ž5s3p2dr5s2p1d. and is reported in Table 1. Ab initio points were calculated for the geometries RŽLi-H 1 . s 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.04 and RŽLi-H 2 . s RŽLi-H 1 . q 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 10.04 , where all the distances are in atomic units. In addition to these, a further set of points at RŽLi-H 1 . s 1.4, 2.1 and RŽLi-H 2 . s RŽLi-H 1 . q 0.7, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 10.04 were calculated to describe accurately the repulsive part of the potential. For all geometries, /ŽH 1 LiH 2 . s 0. These points were then interpolated by a two-dimensional spline function in order to be able to treat the two variables within the dynamical codes discussed in the following Sections, where the specific interpolated reactive surface will be employed. 3. Features of the collinear surface Fig. 1 reports in detail the general behaviour of the collinear arrangement. The plot on the left shows
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
13
Fig. 1. Computed reactive potential energy surface for the collinear reaction. Left: energy level diagrams for the two coordinates. Right: three-dimensional view of the reaction channels with the asymptotic profiles.
the potential energy contours associated to the relevant collinear coordinates, while the three dimensional presentation on the right shows clearly the presence of the large drop in energy on the side of
Fig. 2. Reaction coordinate profile as function of the two relevant internuclear distances along the minimum energy path.
the exothermic reaction Ž9.. The strongly repulsive nature of the interaction at short distances within the system Li-H-H is also seen in both presentations of the 2D surface and will be further discussed in the following sections. Another interesting feature of the 2D surface discussed here is presented in Fig. 2, where the minimum energy path is shown from the lowest electronic state of the Li q H 2 system Žleft. up to the energetically higher, lowest electronic state of the LiH q H system Žright.. One sees immediately that the LiH depletion reaction Ž9. is strongly exothermic with an energy gain on nearly 2.0 eV, while the opposite formation reaction Li q H 2 ™ LiH q H Ž 15 . is obviously strongly endothermic and requires highly excited H 2 molecular partners Žwith the vibrational index ÕX G 4. before it can occur. It is also of interest to notice the slight barrier which appears along the path for the depletion reaction when the H atom is close to the target molecule and the D R s R 2 y R 1 difference becomes negative. This feature will also be discussed in the following Sections. The five occupied orbitals from the single configuration SC calculation included a clearly defined lithium 1s 2 pair, which were overwhelmingly singlet coupled, at all points on the collinear potential sur-
14
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
face. These were easily identified as the orbitals containing the core electrons and there was never any uncertainty as to which SC orbitals should not be correlated further in the SCVB ) stage. The three
remaining occupied orbitals, free to delocalize over all of the atomic centers, changed nature according to the position on the potential surface. At some points they were heavily localized on a particular
Fig. 3. Comparison of the two asymptotic molecular partners as computed from the present VB correlated calculations. See text for definition of the acronyms.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
atomic center, for instance when there was one isolated hydrogen atom, whilst at others they were delocalized over more than one center, for instance when all three atoms were in close proximity. This is a very important feature of the spin-coupled method, since there is no external predetermination of the degree of localization experienced by the orbitals: it results as a natural consequence of the method. The spin-coupling of the electrons in the three valence occupied orbitals also varied with geometry, as the nature of the orbitals changed. At some points a perfect pairing scheme dominated, as two valence electrons became predominantly singlet coupled, whilst at others a more complicated situation was present and more than one mode of spin-coupling was active. This illustrates the importance, especially in open shell systems of this kind, of including the full spin space in the calculation and not just assuming that a perfect-pairing scheme will be a reasonable approximation. A good test of the quality of the calculated points is to consider the asymptotic regions of the potential that correspond effectively to the potential curves for the isolated diatomic species ŽLiH when RŽH 1-H 2 .. is very large and H 2 when RŽLi-H 1 . is very large..
15
Fig. 3 shows the asymptotic LiH potential curve calculated in this work at the single configuration SC and multiconfiguration SCVB ) levels, and compares it with the isolated LiH potential presented in our SCVB ) benchmark calculations w12x. The SC wave function recovers 76% of the true binding energy and overestimates the position of the minimum by ˚ The SCVB ) wave function in this work 0.04 A. recovers a greatly improved 93% of the binding energy, and overestimates the minimum by only ˚ which compares extremely well with the 0.015 A, ˚ respectively. benchmark values of 95% and 0.015 A, The SCVB ) binding energy is slightly lower in this work than for the benchmark as a smaller LiH basis set was used Ž43 contracted functions compared to 87 in the benchmark.. It is clear that including the extra configurations in going from the SC to the SCVB ) wave function has a considerable quantitative effect on the quality of the potential. Fig. 3 also shows that essentially the same story applies to the H 2 asymptotic potential curve. The SC and SCVB ) potentials from this work are shown compared to the position of the experimental minimum w18x. The SC wave function recovers 87% of the H 2 binding energy, whilst the SCVB ) wave
Fig. 4. Computed H-H distance dependence of the interaction using two different levels of correlated calculations as described in the main text.
16
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
function recovers 97%, and both reproduce the position of the minimum almost exactly. The results of Fig. 3 therefore show that the asymptotic regions of the two reaction channels, calculated at the SCVB ) level, are reproduced with very high degree of accuracy and thus indicate that the rest of the potential should also be of a very high quality, as our procedure in no way favours one region of the potential over another. In the non-asymptotic regions of the potential, it is interesting to note that one of the main effects of introducing the excited configurations was to lower the height of a barrier that blocked the path to a potential well as RŽH 1-H 2 . was reduced. This is illustrated in Fig. 4 which plots the energy as a function of RŽH 1-H 2 . with RŽLi-H 1 . fixed at 3 au, a geometry at which this effect is most pronounced. The asymptotic, large RŽH 1-H 2 ., energy has been rescaled to zero for both the SC and SCVB ) curves to facilitate the comparison of the potentials. In reality, the asymptotic part of the SCVB ) curve would sit well below that of the SC curve as there is still an interacting LiH system present. It is plainly seen that a barrier of about 10 mH Ž0.27 eV. at the SC level is reduced to one of just 2 mH Ž0.05 eV. at the SCVB ) level. This is likely to have a significant impact on the low energy dynamics for this reaction and emphasises the importance of including the additional dynamical correlation that the SCVB ) procedure provides.
4. The reactive scattering: a classical treatment We have carried out collinear classical dynamics for the depletion reaction H q LiHŽ Õ . ™ Li q H 2 as well as for the formation reaction Li q H 2 Ž Õ . ™ LiH q H reaction over a range of relative translational energies Ž Etrans . between 0–1 eV. The index Õ represents here the vibrational state of the diatom. For a given Õ and Etrans values the vibrational phase Ž f i . of the diatom was varied systematically in the range between 0–2p . Typically, 1000 trajectories were run to compute the reaction probability Ž P .. The initial separation of the reactants was taken to be in the asymptotic region at 10 a o . The Hamilton’s equations were solved in mass scaled and skewed coordinates w19x using the 4th order Runge-Kutta
method with a step size of 0.21548 fs. The accuracy of the integration was checked by employing the Gear integration along the evolution. In the present study, considering that the earlier calculation for a bent LiH 2 , species w7x suggested its inherent instability, we felt that it was sufficient to analyse the collinear reaction without the effects, which should therefore be minimal, of a bent geometry which dissociates quickly without forming an intermediate structure. 4.1. The depletion reaction: H q LiH(Õ) ™ Li q H2 Since the depletion reaction is exothermic and has practically no barrier, it follows that this reaction goes downhill over the whole reactive path. However, there appears to be a rather small barrier of 0.036 eV in the interaction zone Žas can be seen in Fig. 2.. The reaction is exothermic by 2.0238 eV considering the zero-point energies of the diatoms and therefore one expects to see only a small influence, and then only near threshold, of the above barrier. The calculated probability Ž P depletion . for the reaction is shown in Fig. 5 as a function of Etrans for different Õ-states of LiH. The corresponding probabilities obtained from the collinear time-dependent wave packet ŽTDWP. calculations, which will be described in Section 5, are also displayed in the figure. The classical results are marked with a solid line connecting the circles. The representative statistical error bars are also shown at Etrans s 1.0 eV for all Õ-states. We note the following: Ža. For all Õ states P depletion increases sharply at the low energy, reaching its maximum around Etrans s 0.1 eV. It then drops down to its minimum around 0.2 eV. Žb. With further increase in Etrans the corresponding reaction probability shows a slow monotonic increase, which however remains smaller for higher Õ states. Žc. Below 0.1 eV the classical probability shows some oscillations, which are also well marked when one looks at the wave packet calculations in the same energy region. To analyse the unusual energy dependence of P depletion between Etrans s 0.1–0.2 eV, we have plotted in Fig. 6 a family of Žten. trajectories given in the usual mass-scaled and skewed coordinates Ž x, y . w19x and which also sample uniformly f i of LiH between
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
0 and 2p for different Etrans . Here x represents the relative motion of the atom with respect to the diatom and y denotes the internal motion of the diatom. At Etrans s 0.1 eV, most of the trajectories are seen to end into the reactive channel. At the other energies examined, however, only a few of them find a way to go into the reactive channel, while the majority of them bounces back into the entrance channel after dropping downhill and hitting the repulsive core of the interaction region. Despite the fact that a trajectory has to travel all the way down the reaction path Žsee Fig. 2. for a reaction to occur, the dynamical constraints force it to get back into the entrance channel whenever an excess of translational energy is available. In order to get into the exit channel, therefore, they have to reach the interaction
17
zone with the right phase which then would allow them go into the exit channel, a condition which appears to be achieved only at smaller Etrans . Increasing Etrans does not help to go into the reactive region as a result the nonreactive trajectories which bounce off the repulsive wall at a narrow angle and come back to yield an increasingly higher vibrational excitation content of LiH, as seen, for example, for Etrans s 0.9 eV. To further analyse this special behaviour we have plotted in Fig. 7 the final classical vibrational action Ž n f . of the diatoms as a function of f i of LiH for both reactive and nonreactive cases systematically sampling 100 trajectories between 0–2p for a range of Etrans values. The left panel displays the Etrans values around 0.1 eV where the P depletion shows in
Fig. 5. Computed depletion reaction probabilities as function of collision energy, Etrans , and of the vibrational content of the LiH partner. The classical results are given by the open circles connected by solid lines. The quantum wavepacket results are given by the dotted lines.
18
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
Fig. 5 its first maximum and then a dip, while the right panel shows the results for Etrans ) 0.1 eV. The crosses and the circles indicate the nonreactive ŽNR. and reactive ŽR. trajectories, repectively. We see that, in general, the reactivity is a banded function of f i , that is, within certain range of f i values the reactive events occur while the NR events dominate within another range. The width of the reactive band starts decreasing with the increase of Etrans above 0.1
eV. Even when only a small increase occurs, the R phases dramatically shrink around 0.13 eV and remain roughly the same up to 0.2 eV. With the further increase of Etrans its range then starts to slowly increase. This explains in some way the slow rise in the reaction probability for Etrans ) 0.2 eV. We note that the width of the R band is quite sensitive for collisions at low energies. At Etrans s 0.13, in fact, we observe a marked fractured behaviour between
Fig. 6. Classical calculations of reactive and non reactive trajectories, as function of translational energy, for the depletion reaction. See text for the meaning of symbols.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27 Fig. 7. Computed classical vibrational action Ž n f . of the LiH diatom as function of the phase f i and for different values of translational energy. The circles describe the reactive trajectories while the crosses refer to the nonreactive trajectories.
19
20
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the switch over regions of the R and NR bands. Such switch over regions are known w20x to have fractal characteristics and correspond to the onset of chaotic behaviour for the trajectories. What it means is that they show the existence of a ‘‘smooth’’ banded and a switchover ‘‘chattering’’ region buried within a finer resolution of f i . The trajectories in the chattering regions show multiple collisions in the interaction zone before finally going into either the entrance or the exit channels. We believe that the small oscillations in the classical probability arise presumably because of the existence of nonlinearity in the dynamics at low Etrans , which however becomes less and less so as the collision energy is increased. It is
interesting to note, as we shall further discuss below that the quantum mechanical ŽQM. results also show oscillations at low energies. Finally, it is worth pointing out that in the present 2D reactive process for the depletion reaction the product H 2 molecules appear to be formed in vibrationally excited states. 4.2. The formation reaction: Li q H2 (Õ) ™ LiH q H The formation reaction is an endothermic process which requires an additional energy of 2.06 eV. In order for the reaction to occur, a given trajectory has to then go up along the uphill path Žsee Figs. 1 and 2. and the barrier could be regarded as located in the
Fig. 8. Computed reactive probabilities for the endothermic formation reaction: Li q H 2 ™ LiH q H as function of H 2 vibrational content and of relative translational energy. The classical calculations are given by open circles connected by solid lines while the quantum calculations are given by dotted lines.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
exit channel, i.e., to be of the late uphill type. This suggests that by putting energy into the vibrational mode of H 2 one would more efficiently help the formation reaction than when putting the equivalent amount of energy into the relative translational energy, Etrans . Since we are interested in the low-energy collisions beween 0–1.0 eV, it turns out that H 2 should be at least vibrationally excited in its Õ s 4 vibrational level for the reaction to occur in the above energy range. The computed classical probabilities Ž P formation . for the formation reaction are shown as a function of Etrans , in Fig. 8 for Õ s 4,5,6 and 7. The corresponding probabilities obtained from the quantum mechanical computations are also displayed in that figure. As before, the classical results are marked with a solid line connecting the circles. For Õ s 4 the classical probability rises from the energetic threshold, becomes flat up to 0.7 eV, and then rises with a further increase of Etrans . For Õ s 5,6 and 7 the probability rises sharply at low collision energies, follows a marked decrease and then increases slowly and monotonically with the further rise of Etrans . When comparing the classical results with the QM probabilities, one sees that the agreement is in general quite good, particularly for Õ s 5–7 and the predicted probabilites agree quite well with each other both at low and at high collision energies. For Õ s 4, although the general trend is maintained in the classical results, there seems to be some quantitative discrepancies at the higher values of Etrans . It is interesting to note that the probability of formation is larger at low collision energies for Õ s 5–7. However, for Õ s 4 it becomes larger when the collision energy is increased. Understandably, this effect arises because of the fact that the energy content for Õ G 5 is already nearly equal or above the reaction threshold, and little Etrans can already help a given trajectory to reach the interaction zone with the right ‘‘phase’’ for the reaction to occur efficiently, while it does not seem to be effective when the collision energy is further increased. In that case the trajectories get back into the reactant channel after hitting the repulsive core of the interaction zone as previously discussed. For Õ s 4, an increase in Etrans from the energetic threshold helps the occurence of the reaction since part of the available energy is used up to further climb over the energy barrier.
21
5. The reactive scattering: a quantal treatment 5.1. Computational method The Time Dependent Wave Packet ŽTDWP. method has been widely used in the recent years for reactive scattering processes w21x. Such method is based on the discretisation of the wavefunction on a finite grid, which is propagated in time according to the time-dependent Schrodinger equation. Since the ¨ method is well known, in this section we will only describe the details of the present implementation relevant for the subsequent discussion. In the following the hamiltonian of the system A q BC in its collinear configuration is represented using the mass scaled and skewed coordinates Ž x, y . introduced before w21x and it is given by "2 Hsy
2m
ž
E2 E x2
E2 q
E y2
/
qV Ž x, y.
Ž 16 .
with m s Ž mA m B m C rmA q m B q m C .Ž1r2.. The initial wavefunction C Ž x, y,t s 0. is chosen as a product of a minimum uncertainty Gaussian wavepacket ŽGWP., F Ž x ., which describes the translational part of the wavefunction while the vibration of the diatomic molecule is represented by a Morse oscillator eigenfunction, x Õ Ž y ., corresponding to the vibrational state Õ of the diatomic fragment with energy eÕ
C Ž x , y,t s 0 . s F Ž x . x Õ Ž y .
Ž 17 .
where FŽ x. s
1
Ž 2pd . 1r4 Ž
.
ž
exp y
Ž x y x0 . 4d 2
2
y ik 0 x
/ Ž 18 .
Here x 0 is the center of the wavepacket, placed in the reactant region, while k 0 is the momentum determined from the translational energy of the wavepacket, k 0 s w2 m Ž Etot y e Õ . r" x where Etot is the total energy, corresponding to the center of the wavepacket in the momentum space. The width parameter d controls the spread of the wavepacket in both the coordinate and momentum space: optimal values of the parameter are those which lead to a wavepacket sufficiently narrow in the coordinate
(
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
22
space to be included in the asymptotic reactant region as well as a distribution in the momentum space sufficiently large to allow for the determination of the reaction probabilities for a wide range of translational energies. The time evolution of the wavepacket is carried out using the well-known Split Operator method w22x, where the evolution operator associated with a time independent hamiltonian is approximated as
ž
exp y
iHD t "
/
ž
iTD t
/ ž ž /
f exp y
exp y
2"
=exp y
iVD t "
/
iTD t 2"
Ž 19 .
where T and V are the kinetic and potential energy operator, respectively, and D t is the time step used in the propagation. The action of the kinetic energy operator is evaluated in its local representation by using the Fast Fourier Transform ŽFFT. algorithm. The reaction probability is evaluated using the time-energy mapping of the reactive flux w23x: by fixing a dividing line at y s yI in the product region, the reaction probability is then given by PÕ s
ž
"2 k
mm x < F Ž k . < 2
ž¦
Im C Ž x , yI , E .
/
EC Ž x , yI , E .
;/
Ey
Ž 20 .
where m x s w mAŽ m B q m C .rŽ mA q m B q m C .x, k s 2 m x Ž E y e Õ . r", F Ž k . is the representation in the momentum space of F Ž x . and C Ž x, yI , E . is the Fourier transform of the wavepacket, C Ž x, yI ,t ., calculated along the dividing line y s yI
(
C Ž x , yI , E . 1
s
ž / 2p
1r2
q`
Hy` C Ž x , y ,t . exp Ž iEtr" . d t I
Ž 21 . After the flux analysis, the wavepacket is absorbed by a negative imaginary potential ŽNIP. placed near the edge of the grid in both reactant and product channels. A linear type of NIP, as suggested by
Neuhauser and Baer w24x is used in the present calculation
Vabs
°yiV s~ ¢0
0
R y R abs D R abs
R abs - R F R abs q D R abs R F R abs
Ž 22 . where V0 is the height of the potential, D R abs is its width and R abs is its location, which usually takes values far out from the interaction region. The time propagation is carried out until a sufficiently large fraction of the wavepacket has reached the boundaries of the grid, where it is absorbed by the NIP. 5.2. Details of the calculation In the present work, we have calculated the reaction probabilitites of the Ži. depletion and Žii. formation reactions of LiH in collision with H, in the collinear configuration, as a function of the energy, for different initial vibrational states of the target molecule. Using the formalism introduced above, those reactions correspond to the following way of labelling the three atoms: Ž1. A s Li, B s C s H and Ž2. A s H, B s Li, C s H. In both calculations, we have used a grid in Ž x, y . space composed by Ž128 = 128. points; due to the different masses of the target molecules, on which the scaled and skewed coordinates depend, the grid parameters are different for the depletion and formation reaction. In the former case, the grid starts from ˚ y s 0.74 A˚ with step sizes D x s 0.060 x s 1.02 A, ˚A and D y s 0.035A, ˚ while in the latter case the ˚ y s 0.28 A˚ with step reaction starts at x s 1.17 A, ˚ ˚ sizes of D x s 0.047 A and D y s 0.055 A. The calculations have been carried out for a range of translational energies between 0.2 eV and 1.1 eV, with the initial vibrational level of LiH ranging from Õ s 0 to Õ s 4 in the depletion reaction, while, due to the topography of the surface as discussed before, we have included only excited vibrational levels of H 2 in the formation reaction, from Õ s 4 to Õ s 7. The initial wavepackets have been centered at x 0 f 5 ˚ using different widths in the two reactions, to take A, into account the different distributions of the transla-
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
tional energy: in the depletion reaction, d has been ˚ while we have used d s 0.13 in the fixed at 0.2 A formation reaction. As a consequence of this choice, we have been able to calculate reaction probabilities for an energy range of 0.3 eV with a single propagation for the two reactions under study.
23
The dividing lines for the calculation of the flux ˚ for the depletion have been fixed at yI s 4.0 A reaction and at yI s 4.62 for the formation reaction. Together with the main dividing line, other two lines have been placed 10 steps before and after the main dividing line in order to check, from the values of
Fig. 9. Snapshots of the propagation of a wavepacket for the depletion reaction of LiH, with LiH in its initial vibrational level Õ s 1 and an initial mean translational energy of about 0.3 eV.
24
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
the reaction probabilities, that the position of yI has been placed well inside the product asymptotic region. The location of the absorbing potentials in the product channels have been chosen to be 15 steps from the dividing line used to calculate the flux; in
˚ for the reactant channels, we have used x abs s 6.0 A ˚ for the the depletion reaction and x abs s 5.32 A reactant channels with a height of about 1 eV in both cases. Finally, a time step of 0.04 fs was chosen for the depletion reaction, reduced to 0.02 fs in the forma-
Fig. 10. Snapshots of the propagation of a wavepacket for the reaction H 2 q Li, with H 2 in its Õ s 4 level and an initial translational energy of about 0.9 eV.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
tion reaction: the latter value has been taken to be smaller than the former to take into account the larger energy content of the wavepackets used in the formation reaction. It is worth mentioning that such numerical parameters have been the results of a rather extensive number of test calculations for numerical convergence. The propagation has been carried out until the 99.9% of the wavepacket has been absorbed; such large fraction ensures the convergence of the reaction probabilities on the number of chosen time steps. We have found convenient to interface the program of the time dependent propagation with a graphical package that allows us to watch the propagation of the bidimensional wavepacket during the execution of the program. In order to have some insight on the physical processes under investigation, we have reported in Fig. 9 a series of snapshots of the propagation of a specific wavepacket describing the depletion reaction, with initial vibrational state Õ s 1 of LiH and a small initial translational energy, at equally spaced instants of time. From the reaction probabilities reported in Fig. 5 such wavepacket contains energy components placed around the minimum of the curve. We can see that, due to the nature of the potential, the wavepacket is not captured in the interaction region therefore the reaction appears to be a direct process. Furthermore, we can observe from the shape of the wavepacket in the product channel that H 2 is formed with quite a large vibrational energy content. At any rate, we can also see that, after the collisional event, the wavepacket in the product region has a different shape from that at t s 0. It is very probable that all the available vibrational levels of LiH are populated after the collision so that there is a reduction of the average translational energy and the wavepacket moves now moves slowly towards the asymptote of the reactant region. In such a case a very large number of time steps is needed to reach the end of the propagation. The convergence of the reaction probabilities in the cases where the total energy is close to the threshold of the reaction needs to be carefully examined: to this aim, we have carried out a test calculation using a grid of Ž256 = 256. points, for a wavepacket with the same initial translational energy as the one reported in Fig. 9 and with the lowest vibrational state, Õ s 0, in order to have the minimum total energy: we have found that
25
the reaction probabilities are in good agreement with those obtained with the the less dense grid. We have also reported further snapshots for the formation reaction in Fig. 10, for a wavepacket with H 2 in the vibrational state Õ s 4 and average translational energy of about 0.9 eV. Even in this case the reaction is direct; furthermore, we can observe that the reactive flux is due to the outer part Žlarge y . of the initial wavepacket, while the inner part remains almost unaltered by the collision, the opposite of what we have found in the reverse, depletion reaction discussed above. It is also interesting to note that the corresponding depletion and formation probabilities reported in Figs. 5 and 8 appear to be in very good agreement with the classical results. In other words, we see from the depletion data of Fig. 5 that the quantum calculations also suggest a direct reaction where the excess energy is not always driven into the translational component: little additional relative energy content, in fact, is able to make the reaction to occur with large probabilities no matter what the vibrational level of the LiH is initially chosen to be. As the collision energy is increased, on the other hand, the efficiency of the reaction decreases rather markedly and only picks up rather slowly especially for initially excited, vibrationally ’hot’ partner molecules. The situation is shown to be similar for the formation reaction of Fig. 8 where, however, we see larger discrepancies between quantum and classical calculations at low energies possibly indicating the more significant role of vibrational coupling, in the quantum case, between Etrans and H 2 energy content.
6. Concluding remarks In the present work we have approached, for the first time, the study of the depletion and formation reactions of astrophysical interest which involve LiH molecules. Starting from first principles we therefore have computed over a very extensive grid of points a highly correlated reactive surface for the two-dimensional collinear processes relative to the above two reactions. After accurate numerical interpolation of the initial points, we have carried out classical trajectory calculations for the reactive case and also time-de-
26
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27
pendent, wavepacket propagation studies for the quantum reactive process. The results from the two approaches, and a comparison of their computed probabilities, shows the following: Ži. The depletion reaction is strongly exothermic and appears to be, in the 2D case, a direct reaction without any transition state region. Žii. The shape of the reactive potential energy surface shows that low translational energy content markedly favours the depletion reaction and that its increase leads initially to very small reaction probabilities as the reagents ‘bounce back’ into the entrance channel due to the special repulsive features of the skewed coordinate system for the partner masses. Žiii. The initial vibrational content of LiH in the depletion reaction plays a rather minor role Žsee Fig. 5. in the sense that the final reaction probability depends very little on it. Živ. Both classical and quantum treatments at low and intermediate energies agree remarkably well with each other and show the time-dependent approach to be only adding a more marked interference structure in the very-low energy regions. Žv. The formation reaction requires highly excited H 2 molecules for it to occur at low collision energies. Both classical and quantum results indicate, in fact, that the H 2 vibrational levels Õ s 4,5,6,7 provide very efficient starting partners for the endothermic reaction to take place. Žvi. In the formation reaction the calculations show once more the rather inefficient coupling of vibrational and translational forms of energy already observed before for the depletion reaction: low Etrans values markedly favour reaction even when excess energy would be available from Evib . On the whole, this preliminary study gives us some useful indications on what we should expect from a more extensive calculation of the full 3D reaction probability: the quantum effects are not seen to be important for the depletion process, a fairly direct reactive event without complex formation. Furthermore, we see that internal energy content of LiH is not an helping factor for the process to occur and therefore one should expect that ‘cold’ LiH molecules can be depleted by reactive encounters with H atom at fairly low values of translational energy. With the same token, the formation of LiH
from vibrationally excited H 2 molecules appears to lead to the formation of vibrationally excited LiH, as seen from the time-dependent snapshots described before. Clearly one further needs to analyse the behaviour of the full 3D potential for the reaction rates and therefore the rather delicate process of fitting the available array of Born-Oppenheimer points is currently under way in our group. The preliminary collinear results, on the other hand, are already able to tell us a great deal on the specific features of the reaction and on its likely success to be realistically described by classical calculations. As a matter of fact, and as mentioned earlier, the computed instabilities of the bent structures w7x and the absence of complex formation in the 2D study conjure up a physical picture where head-on collisions greatly favour reactions, and therefore one should not expect to find that 3D calculations strongly ‘push away’ flux from the collinear path.
Acknowledgements The financial support of the Italian National Research Council ŽCNR. and of the Italian Ministry for University and Research ŽMURST. is gratefully acknowledged. N.J.C. thanks the European union for a training grant under the TMR project. S.K. wishes to thank the Max-Planck-Gesellshaft for supporting his stay at the University of Rome through the funding of a Max-Planck-Research-Prize to F.A.G. and E.B. also wishes to thank Dr. Pietro Lanucara for his generous help with graphic interfaces.
References w1x V.K. Dubrovich, Astron.&Astrophys. Trans. 5 Ž1994. 57. w2x R. Maoli, F. Melchiorri, D. Tosti, Astrophys. J. 425 Ž1994. 372. w3x S. Lepp, J.M. Shull, Astrophys. J. 280 Ž1984. 465. w4x V.K. Dubrovich, A.A. Lipovka, Astron. &Astrophys. 296 Ž1995. 301. w5x P.C. Stancil, S. Lepp, A. Dalgarno, Astrophys. J. 548 Ž1996. 401. w6x E. Bougleux, D. Galli, private comm. w7x A.I. Boldyrev, J. Simons, J. Chem. Phys. 99 Ž1993. 4628. w8x J. Gerrat, M. Raimondi, Proc. Poy. Soc. London A 371 Ž1980. 525.
N.J. Clarke et al.r Chemical Physics 233 (1998) 9–27 w9x D.L. Cooper, J. Gerrat, M. Raimondi, Int. Rev. Phys. Chem. 7 Ž1985. 59. w10x J. Gerrat, D.L. Cooper, M. Raimondi, in: D.J. Klein, N. ´ ŽEds.., Elsevier, Amsterdam, 1990, p. 287. Trinastic w11x D.L. Cooper, J. Gerrat, M. Raimondi, Chem. Rev. 91 Ž1991. 929. w12x N.J. Clarke, M. Raimondi, M. Sironi, J. Gerratt, D.L. Cooper, Theor. Chim. Acta 95 Ž1997. 1328. w13x P.O. Lowdin, J. Chem. Phys. 19 Ž1951. 1396. ¨ w14x R. McWeeny, Methods of Molecular Quantum Mechanics, Academic Press, London, 2nd edition, 1992. w15x S.M. Goldfield, R.E. Quandt, H.F. Trotter, Research Memorandum No. 95, Economic Research Program, Princeton University, 1968. w16x N.C. Pyper, J. Gerratt, Proc. Roy. Soc. Lond. A 355 Ž1977. 406.
27
w17x H. Partridge, C.W.B.Jr., P.E.M. Siegbahn, Chem. Phys. Lett. 97 Ž1983. 198. w18x G. Herzberg, J. Mol. Spectry 33 Ž1970. 147. w19x P.J. Kuntz, in: W.H. Miller ŽEd.., Dynamics of Molecular Collisions Plenum, 1976. w20x V. Balasubramanian, B.K. Mishea, A. Bahel, S. Kumar, N. Sathyamurthy, J. Chem. Phys. 95 Ž1991. 4160. w21x N. Balakrishnan, C. Kalyanaraman, N. Sathyamurthy, Phys. Rep. 280 Ž1976. 129. w22x J.A. Fleck Jr., M.D. Feit, J. Chem. Phys. 78 Ž1982. 301; M.D. Feit, J.A. Fleck Jr., A. Steiger, J. Comp. Phys. 47 Ž1982. 412; D. Kosloff, R. Kosloff, J. Comp. Phys. 52 Ž1983. 35. w23x N. Balakrishnan, N. Sathyamurthy, Chem. Phys. Lett. 240 Ž1995. 119. w24x D. Neuhauser, M. Baer, J. Chem. Phys. 90 Ž1989. 4351.