A bond-order LiFH potential energy surface for 3D quantum-mechanical calculations

Volume 143, number 2

CHEMICAL PHYSICS LETTERS

8 January 1988

A BOND-ORDER LiFH POTENTIAL ENERGY SURFACE FOR 3D QUANTUM-MECHANICAL CALCULATIONS

A. LAGANA Dipartimento di Chimica deK!Jniversit~, 06100 Perugia, Italy

0. GERVASI Centro di Calcolo dell’Universitri,06100 Perugia_Italy

and E. GARCIA Departamento de Quimica Fisica de I’Universidad, Bilbao, Spain Received 15 October 1987

A new tit to ab initio estimates of the LiFH potential energy has been obtained for use in a full three-dimensional quantummechanical calculation of the reactivity of this system.

1. Introduction

We have recently discussed methods of fitting [ l-41 ab initio potential energy values. Although

most of the proposed schemes are applicable to general atom-diatom reactions, our main concern has been for systems having a heavy central atom and only one reactive channel open. As a particular application, checks on the accuracy of surfaces fitted using the abovementioned procedures have focused on some alkali/alkaline earth+ hydrogen halide systems. In the particular case of the LiFH reaction, the quality of the different fits to the ab initio points of Chen and Schaefer [ 51 has been investigated using classical trajectories. As a result of these tests [ 21 a piece-wise (CSGL) potential energy surface (PES) has been found to lead to the best agreement between calculated and measured reactive properties

161. The CSGL PES is designed mainly for use in those reduced dimensionality quantum methods [ 7- 121 based on the assumption that the angle of approach 174

can be held constant when integrating scattering equations. For this reason, the GSCL potential functional is formulated as a function of two internuclear distances (HF and LiF) and of the induced angle 19 (defined as the smallest LiFH angle. It has to be pointed out here that for LiFH this angle and the approaching angle are almost coincident). The choice of these coordinates is motivated by the fact that the LiH channel is closed at energies of interest and that HF and LiF are the only two bonds which can be broken or formed during the reaction. The funo tional of the CSGL PES partitions the plane defined by the LiF and HF distances TL~F and THFat a fixed value of B into four regions. For all regions diatomic potentials ( Vz( rr )) are expressed as a polynomial of the third order in the displacement of the diatom from its equilibrium distance damped by an exponential in the same variable [ 131. In contrast in each region of the fixed 8 plane the three-body term ( V3(r,, r,, r3)) of the multi-body expansion of the potential [ 141 is given a different functional representation according to the local character of the interaction. In the inner (strong in-

Cl009-26 14/88/S 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

Volume 143, number 2

8 January 1988

CHEMICAL PHYSICS LETTERS

teraction) region VXis expressed as a polynomial in the inverse power of the internuclear distances. In the outer (long-range interaction) regions the threebody term is represented by either a LiF of a HF (or a mixture of both) scaled diatomic functionals depending on the closest asymptotic fragment. This choice, inevitably, builds into the functional interpolator the impossibility of properly describing the formation of the LiH product and the insertion of Li into HF. It is obvious that such a bias of the potential representation is unimportant for the usual quantum reduced dimensionality approaches in which the collision angle is held constant (i.e. when rotations are considered as decoupled from other degrees of freedom). In addition, as already commented above, the CSGL potential has shown itself to be adequate for trajectory calculations. This happens both because the formation of LiH is energetically forbidden at energies of interest and because, due to the bent geometry of the transition state, the lithium atom can attack directly its tinal partner F even when approaching from the H side. On the contrary, as discussed in more detail elsewhere [ 15 1, for a full 3D quantum scattering calculation [ 16-221 two main difficulties arise when using this piece-wise potential functional. The first is due to the fact that the potential energy cross section constructed at constant value of the propagation coordinate sometimes shows spurious holes. These spurious features, although classically unaccessible, demand large computational resources at the expense of the part of the potential energy surface important for reactive dynamics. A second difficulty originates from the fact the CSGL PES allows an appropriate asymptotic expansion for the product channel only at large values of the propagation variable (both these features particularly affect hyperspherical coordinate treatments). For these reasons during the preliminary stages of a full 3D quantum mechanical investigation of the Li+ HF reactivity we had to abandon the CSGL PES in favour of a global representation of the potential energy surface [ 151. In this paper the litting procedure and the graphical analysis of the new PES is reported. Its functional form and parameters are also reported for use in exact and approximate benchmark calculations.

2. The fitting procedure The procedure followed for fitting the new PES is similar to that used for previous interpolations of the potential energy surfaces of other alkali and alkaline earth+ hydrogen fluoride systems [ l-4,23,24]. Accordingly, ab initio values have been modified to make asymptotic limits reproduce dissociation energies of diatomic fragments and energetics of the reactive process. A mrther correction of the calculated points has been introduced to obtain a realistic value of the barrier to reaction. At the same time an extensive investigation of the accuracy obtainable from this type of CI ab initio evaluation of the barrier to reaction has been carried out for the Li + HCl system [ 231. It has shown that at the level of accuracy used for evaluating the LiFH potential energies [ 51, the geometry of the transition state can be considered already converged. In contrast, the same calculations have shown that related energy values are still quite sensitive to an extension of both the orbital basis set and the CI expansion. Therefore, as also suggested by the value of the measured threshold energy of this reaction, empirical corrections are needed to lower the transition state barrier to a realistic height. A correction of 3 kcal/mole for the transition state energy of this system was adopted by Murrell and coworkers for their fit [ 251. In our trajectory calculations we have also tested surfaces generated by applying an additional correction of 3 kcal/mole (the CSGL surface is one of them). As already mentioned, calculated quasiclassical cross sections showed a better agreement with experimental tindings when estimated on potential energy surfaces having the lower barrier to reaction. However, more recent approximate quantum calculations [ 261 performed according to a light-heavy-light 10s scheme [ 9,271 have shown that the reactivity of this system is higher than indicated by trajectory calculations because tunneling effects are quite important. Therefore, for the transition state of the new fit we. have adopted the same correction as ref. [ 251. As a second step, diatomic terms were fitted using a fourth-order expansion in bond order variables V2(ri)=-Dei

c ajmnT ) M

(1)

where ni is exp[ - bi(ri-rci)] and the parameters Dci, 175

Volume 143, number 2

CHEMICAL PHYSICS LETTERS

8 January 1988

Table 1 Parameters for the diatomic BO terms. Energies in kcal/mole, distances in 8,

LiF FH LiH

137.59 141.20 58.00

bi

r,,

0.9708

1.5639 0.1968 1.5997

2.1942 1.1709

ai,,, and bi are determined by solving a system of equations obtained by forcing the reproduction of spectroscopic constants [ 28 1. Values of these parameters for LiFH are reported in table 1. As a final step V3was calculated for all geometries considered by subtracting the value of the three dia-

2.3044 2.0781 2.1512

ai

ai3

ai

-2.3017 - 1.2567 - 1.3822

1.6903 0.2791 0.3107

- 0.6930 - 1.0005 -0.0797

tomic terms from the modified ab intio points. The resulting values were again fitted using a polynomial expansion in the bond order variables, V,(r1~rZ!~r3)=‘C C C j

k

cjkldn%ni

(2)

*

I

Table 2 Parameters for the three-body BO term (in kcal/mole) j kl

j kl

c Jkf

O.l046546E+04

122

-0.1007101E+03

011 210 120

0.6041740Et03 0.4203070E t 03 -O.l349205E+04 -0.3248449Et03

032 203 113

-0.1738752Et03 -0.5171799E+03

201 1 1 1 021

-0.10973268+04 -0.7902607EtO3 -O.l262961E+04

102 012 310

-0.5112102Et03 -0.1446229Et03 0.6956992E+03

220 130 301 211

0.1587417EtO4 -0.1675520Bt04 0.8526413E+03 0.5228272E t 03 0.7727870E+03 0.2121133E+04

110 101

121 031 202 112 022 103 013 410 320 230 140 401 311 221 131 041 302 212

176

c,

0.1090760E t 04 0.6315115E+O2 O.l851443E+03 0.4547308Et02 0.1704061Et03 -O.l625348E+03 -0.3074558E t 03 -0.73838618$03 0.1750511Et04 -0.2350221Et03 -0.7683703EtO3 -0.4342706Et03 -0.2941487E+03 -0.1569613Et04 -0,81622398+03 0.9306138Et03

023 104 014 510 420 330 240 150 501 411 321 231 141 051 402 312 222 132 042 303 213 123 033 204 114 024 105 015

-0.1514158Et03

-0.1152083Et02 0.4385045E t 02 -0.5360226E+Ol 0.1823785EtOl -O.l471057E+03 0.1173729Et03 O.l115220E+03 -0.4845248Et03 0.1004951Et01 0.4434655E t 03 0.5675508Et03 -0.4708619Et02 0.329444OE t 02 0.40126OOEtO3 0.2076121Et03 -0.8350344Et03 -0.4439566Et03 0.17636028 t03 0.4160904Et02 0.72979608 t02 0.5640212EtO3 0.9578297E+O2 -0.4226320Et02 -0.15889948+00 -0.10909348+03 0.1695968Et02 -O.l066076E+O2 O.l237160E+02

_

CHEMICAL PHYSICS LETTERS

Volume 143, number 2 Table 3 RMS deviations of the BO PES (in kcal/mole) Energy

rms

O-20 20-40 40-60 60-80 80-100 >lOO

1.02 2.10 3.85 2.23 3.51 4.91

Best fit values of the c,k/paratIIeterS are reported in table 2. The accuracy of the fitted surface is illustrated in table 3 by reporting rms deviations of the fitted potential from modified ab initio values for different energy bands. As apparent from the table, the lower energy regions of the BO PES give the smaller contributions to the overall root-mean-square deviation.

3. An analysis of the BO surface 4 peculiarity of the BO functional is the exponential nature of the BO variables. This characteristic ensures a natural falloff of the three-body interaction as expected from physical arguments when moving from equilibrium to large distances. As a consequence there is no need for damping or switching functions, whose use is critical for introducing spurious structures in polynomial fits. This makes the 30 PES gccurately fit the potential near the minimum energy path (MEP) and connect smoothly asymptotes and internal regions. However, if no constraints are set for the shortrange regiqp the value of the BO PES may oscillate at very close distances. To prevent this kind of unrealistic behaviour one can either force analytically a smooth evolution of the interpolating functional or impose additional constraints at short distances by adding in regions of interest a set of values empirically derived. Although the addition of further points can easily be done by means of a careful local extrapolation of calculated points, this procedure has the disadvantage of lowering the accuracy of the fit in those regions of the potential energy surface which are important for reactive dynamics. A different approach to forcing a smooth behaviour into the short-

8 January 1988

range potential consists in acting directly on the functional form to smooth it down. This can be achieved either by adopting for the expansion of the potential a different set of variables or by freezing the value of the BO V, at some limiting value at negative displacements. The first choice has been developed for applications to the Mg+ HF reaction for which a variable falling off at both short and long distances has been used [ 41. For the Li+HF reaction the second choice was adopted. Such a choice takes into account the fact that at all 8 values the minimum energy path is located at internuclear distances larger than LiF and HF equilibrium distances and that, as a consequence, at distances less than the equilibrium the potential starts being repulsive. A freezing of V, at is equilibrium value is equivalent (apart from a constant shift) to the imposition of a diatomic-like repulsive behaviour to the short-range part of the fitted PES. In this way, in addition to ensuring a physically realistic behaviour of the surface at all distances, the fit of the potential energy can be optimized to reproduce the features of the low-energy region of the PES. Plots of the resulting BO PES (including the shortrange correction to V3) are reported in figs. 1 and 2. In the first figure minimum energy paths calculated at fixed values of B are plotted as a function of the angle a! formed by the straight line used for cutting the reactive channel for the local search of the minimum of the potential energy. As is the case of other reactive systems the MEP of the BO PES well reproduces the features of calculated potential energy values and of the corrections introduced for having a more accurate surface without introducing spurious structures. In particular, the lower left-hand-side (lhs) panel of fig. 1 shows that the collinear MEP goes through a small well before climbing the high barrier (of about 20 kcal/mole) to the products. This barrier lowers and the early well deepens when moving to perpendicular geometries (second panel from the bottom at the lhs). The minimum barrier height is reached at 8 = 74” (top lhs panel) while at this 6 value the well has already become unimportant. The right-hand-side (rhs) panels show MEPs related to collision for more bent geometries. Going from the lower panel up the MEPs for ti9=45,35 and 30” are reported. All show a high barrier (in increasing or177

Volume 143, number 2

CHEMICAL PHYSICS LETTERS

8 January 1988

35-

_5;,.,.,,,,,,,,,,,. 0.0

-511 0.0

20.0

,

,,I 40.0

,,',,,!/,.!,, 20.0 40.0

60.0

60.0

55~;7w.,:.771~.

80.0

,.,! 80.0 ,"J !

/ 1 r 1

i 1

c

i

1 15p P

j

t : 3=-A-;

1

-5 0.0

Q

I

20.0

40.0

60.0

80.0

clegree

Fig. 1. Minimum energy paths of the BO potential energy surface for B = 180, 90 and 74” (from bottom up Ihs panels) and t9=45, 35 and 30” (from bottom up rhs panels) reported as a function of o.

der) which makes strongly bent attacks more classically forbidden. A more pictorial description of the potential energy surface is given in fig. 2, where isoenergetic contour maps of the PES have been drawn as a function of the reactant and product diatom internuclear distances. From these contours the overall smoothness of the surface can be better appreciated. In this way, more quantitative information can be obtained about the steepness of the walls of the reaction channel and other critical features of the PES. As an example, from this figure both the gradual lowering of the bar178

rier with 19(in the interval 180”> r9> 90” going from the lower to the upper lhs panel) and the decrease (followed by an increase) of the early well in the same interval can be clearly seen. From the variation with 19the location and steepness of the barrier can also be estimated. In the rhs panels of the figure, going from the lower to the upper panel, contours are reported for 8 =45, 30 and O”, respectively, and it is apparent that the introduced correction to Vs is able to impose a realistic behaviour to the short-range part of the BO PES so that HF and LiF channels separate smoothly when r9 decreases.

8 January I988

CHEMICAL PHYSICS LETTERS

Volume 143, number 2

0.5 1.0

1.5

2.0

2.5

3.0

S.5

4.0

4.5

5.0

1.0

R(LiF) /

1.5

2.0

2.5

5.0

3.5

4.0

4.5

5.0

A

Fig. 2. Isoenergetic contour maps of the BO potential energy surface for r3= 180,90 and 74” (rhs panels from bottom up) and B = 4530 and 0’ (lhs panels from bottom). Energy cuts have been drawn at 4.5,11.5,34.5 and 85 kcal/mole. The energy zero was set at tbe entrance channel asymptote. Potential energy values higher than 85 kcal/mole have been shaded.

4. Summary

In this paper we have presented a new tit to the LiFH potential energy surface using a bond-order variable expansion. Such an expansion ensures that

and long-range regions of the fitted PES behave smoothly. To make the fitted surface suitable for 3D quantum-mechanical calculations it has been necessary to smooth the three-body term of the BO expansion at short distances. Values to be fitted were

intermediate-

179

Volume 143, number 2

CHEMICAL PHYSICS LETTERS

derived from ab initio points by forcing the reproduction of asymptotic properties and by lowering the transition state barrier to a more realistic value. The correction of the ab initio points in the transition state region has been guided by theoretical considerations based on the limitations of the basis set and the CI space as well as by approximate dynamical (classical trajectories and 10s quantum) calculations. The graphical analysis has confirmed that his PES well reproduces (modified) ab initio values and behaves in a physically reasonable way suitable for three-dimensional calculations.

Acknowledgement One of us (EG) thanks the Basque Government for a fellowship.

References [ 1 ] E. Garcia and A. Laganii, Mol. Phys. 52 (1984) 1115. [ 21 J.M. Alvarifio, M.L. Hemandez, E. Garcia and A. Lagana, J. Chem. Phys. 84 (1986) 3059. [ 31 E. Garcia and A. Lagana, Mol. Phys. 56 (1985) 629. [4] A. Lagana, M. Paniagua and J.M. Alvariflo, in preparation. [ 51M.M.L. Chen and H.F. Schaefer III, J. Chem. Phys. 72 (1980) 4376. [ 61 C.H. Becker, P. Casavecchia, P.W. Tiedemann, J.J. Valentini and Y.T. Lee, J. Chem. Phys. 73 (1980) 2833. [ 71 M. Baer and D.J. Kouri, in: The theory of chemical reaction dynamics, ed. DC. Clary (Reidel, Dordrecht, 1986) p. 167.

180

8 January 1988

[ 81 J.M. Bowman and A.F. Wagner, in: The theory of chemical reaction dynamics, ed. D.C. Clary (Reidel, Dordrecht, 1986) p. 47. [ 91 D.C. Clary and J.P. Henshaw, in: The theory of chemical reaction dynamics, ed. D.C. Clary (Reidel, Dordrecht, 1986) p. 359. IO] R.B. Walker and A.F. Hayes, in: The theory of chemical reaction dynamics, ed. D.C. Clary (Reidel, Dordrecht, 1986) p, 105. 111 D.L. Miller and R.E. Wyatt, J. Chem. Phys. 86 (1987) 5557. 121 B.M.D.D. Jansen op de Haar and G.G. Balint-Kurti, J. Chem. Phys. 85 (1987) 2614. 131 K.S. Sorbic and J.N. Murrell, Mol. Phys. 29 (1975) 1387. 141 J.N. Murrell, S. Carter, P. Huxley, SC. Farantos and A.J.C. Varandas, Molecular potential energy functions (Wiley, New York, 1984). 15 ] A. Lagan&, R. T Pack and G.A. Parker, in preparation. 161 A. Kuppermann and P.G. Hypes, J. Chem. Phys. 84 (1986) 5962. [ 171 G.A. Parker, R.T. Parker, B.J. Archer and R.B. Walker, J. Chem. Phys., to be published. [ 181 F. Webster and J.C. Light, J. Chem. Phys. 85 (1984) 4744. [ 191 K. Haug, D.W. Schwenke, Y. Shima, D.G. Truhlar, J. Zhang andD.J. Kouri, J. Phys. Chem. 90 (1986) 6757. [20] A. Kuppermann, G.C. Schatz and J.P. Dwyer, Chem. Phys. Letters45 (1977) 71. [2 1 ] J. Lindenberg and B. Vessal, Intern. J. Quantum Chem. 3 1 (1987) 65. [22] L. Wolniewicz and J. Hinze, J. Chem. Phys. 85 (1986) 2012. [231 P. Palmieri, E. Garcia and A. Lagani, J. Chem. Phys., to be published. [ 241 J.J. Margarido, M.L. Hemandez, J.M. Alvariflo and A. Lagana, in preparation. [ 251 S. Carter and J.N. Murrell, Mol. Phys. 41 (1980) 567. [26] A. Lagana and E. Garcia, Chem. Phys. Letters 139 (1987) 140. [27] D.C. Clary, Mol. Phys. 44 (1981) 1083; Chem. Phys. 81 (1983) 379. [28] E. Garcia and A. Laganh, Mol. Phys. 56 (1985) 621.