- Email: [email protected]
A spline-fitted potential surface for bent triatomic molecules
Volume 48, number 1
CHEMICAL PHYSICS LETTERS
15 Lfay 1977
A SPLINE-FITTED POTENTIAL SURFACE FOR BENT TRIATOMIC MOLECULES Stephen
K. GRAY, James S. WRIGHT $
Department of Chemistry, Carleton University, Ottawa, Canada KIS iB6
and Xavier CHAPLJJSAT Laboratoire de Chimie Thebrique $$, Universitt de Paris-Sud, Centre d’Orsay. 91405 Orsay, France Received 24 January 1977
The rotated-Morse curve cubic spliie method developed previously is extended to bent triatomic molecules by usinS 25 cubic splmes. The Yates-Lester potential for bent Hs is shown to be accurately fitted over the entire raqe of the three internal coordinates, with a standard deviation of less than 1 kcal mol-’ . The spline method compares favorably in computa-
tional speed with the analytical potentiaL
1. Introduction A general procedure to interpolate ab initio potential energy surface data is needed for studies in molecular dynamics [l] . Wall and Porter [2] developed a rotated Morse curve @MC) representation for producing collinear triatomic potential surfaces. Connor et al. 131 and Bowman and Kuppermann [4] expressed the RMC parameters as cubic spline functions, and by so doing obtained smooth potential surfaces with desired reaction path and saddle point characteristics. Gray and Wright [5] showed that the PMC-cubic spline method could give an accurate fit to an ab initio potential surface for collinear Hs, and that classical trajectories on the spline-fitted surface were in good agreement with those obtained by using the analytical SSMK function [6]. Gray and Wright [5] also found the trajectory computation time to be comparable for splme-fitted and analytical surfaces, thus demonstrating the attractiveness of the method for trajectory calculations involving only two internal coordinates RA, and
ordinates were deveioped by. lMcLaugh!in and Thompson [7] and extended by Sathyamurthy et ai. [S,9] _ En each case purely numerical grids along the internal coordinates were used. The 15 X 15 X 15 grid chosen by Sathyamurthy and Raff [8] to describe the systems (He + Hz) and (D -I- HCI) by 3D cubic spline interpotation proved to be relatively inefficient, both with respect to accuracy and computational speed of the classical trajectories- Some purely numerical studies using cardinal splines have also been reported [IO,1 I]. In this paper we outline an extension of the co& linear RMC-cubic spline method given previously [3--51 and apply it to the case of three internal coordinates [12] potentiaI RAB, RBC, and R,c- The Yates-Lester surface for bent H, is used as an ihustrative exampfe. A previous RMC potential surface for the Hz system was used by Csizmadia et al. [ 131 and included nonlinear geometries; however, no details were given and a large standard deviation was reported. The use of generahzed Morse functions [ 14,l S] to interpolate potentiaf surface data is also relevant to the present work.
RBC-
Sphne-fitted
potential
surfaces for three internal
co2. Method of calculation
$ To whom correspondence should be addressed. $* Associated with the CNRS (ERA no. 549).
For the three-particle
system
A... B...C let y be the
CHEMICALPHYSICSLETTERS
Volume48, number 1
ABC angle. The three internal coordinates are then y, RAB, andRBC, and the allowed range of y is from 180” (Linear) to 60°. For angles much smaller than 60° any RMC surface must become pathological [ 151; this problem can be avoided by redefining the internal coordinates so that 7 > 60” aIways. This is always possible since one of the interior angles of triangle ABC is necessarily greater than or equal to 60’. For our purposes here it is sufficient to consider just geometries contained within the range 60° G y < 180°. A direct extension of the collinear RMCcubic spline potential [3--51 is outhned beIow. At each specific 7, an (RAB, RBC) plot is drawn using orthogonal RAB and RBc axes, exactly as in the co&near case. Now 6 is the swing angle for the rotated Morse curve, defined as the clockwise angle of rotation about the swing point (RAB = RBC = 5 A) from the reference line RBc = 5 A. The potentiaI energy is: F(R,B,RBc,R_&= =D@,
F(e, Y, Z)
y)Y)Icl - exp[P@,
yXZ - Z,(e,
face were given by Porter and Karplus [ 171. These surfaces are smooth at y = 180° and 120° but show discontinuities in the derivatives (kinks in the potential) at +y= 60”. As before [S] ,3 I vahres of 0 were used ranging from 0 to 90’. The spacing of nodes was every 5’ from O-35” and every lo from 38-45O (and similarly for 46-900); this choice of nodes gave reasonable incremeni in the potential energy. It was found that only 5 values of y were needed to accurately fit the surface: 60”, ZZO”, 180”, 240”, and 300”. The Iatter two nodes are identical by symmetry to the first two and are used to give the proper curvature to the spline function at y = 180°; a similar argument was used previously [5 ] at e = 45O. The resulting 5 X 31 arrays for o(e, y), p(f7, y) and Z,(0,7) were fitted by 2D natural cubic spline functions ([7] , and references therein)_
3. Results and discussion y))l I2 - 11 ,
where tan 8 = (R& - RBC)/(RiB - RAB), 0 G 0 < 90”, (RiBI R&) = origin of (8, I) coordinate system (swing point), and C= [(RLB - RAB)2 f (REc - R,C)2] ‘j2. The Morse parameters D, P and Zeq, which in the collinear case were functions only of 6, now become func-
tions of 8 and -y. For various discrete values of y we obtain the specific values for D, ~3and Zes in an identical manner to that used for the cohinear geometry [S] . Briefly, we choose a discrete set of f3 and scan down rays of constant 0 to find Zes and D; the force constant k is obtained numerically at Zes and, from this, p = (k/2D)‘12. When this process is completed we have the required two-diiensional arraysD(B,7),Z,_&,?) and P@,y). We now fit these arrays to two-dimensional naturaI cubic splines using the
method given by McLaughIin and Thompson (71. To fuhy determine the sphne functions it is necessary to specify either the first derivativespr the second derivatives along the grid boundaries in e and 7. For simplicity, we set the second derivatives equal to zero along these grid boundaries. AS a test for the method we chose the YatesLester [12] analytical surface for H3, which is based on the collinear SCFCI calculation of Liu [ 16 ] and uses a Porter-Karplus [17] formulation to ahow for non-linear geometries. Contour maps for a similar sur156
15 May 1977
2000 points were generated randomly in the range 0.5
A comparison of the Yates-Lester (YL) potential and spline potential (VSPLINE) at these points showed a standard deviation of 0.70 kcal mol-I. To investigate the %hemicaIIy interesting” region of collision energies below about 2 eV, we discarded points leading to a YL energy above -55 kcal mol-r. The remaining 685 points showed a standard deviation of 0.38 kcai mol-l. This can be considered to be very good agreement, particularly since only 5 nodes in y were used. For comparison, the standard deviation of Csizmadia et al. [ 131 on the Hz surface was reported to be 6.65 kcal mol-l. Sathyamurthy and Raff [8], using 3375 nodal points to fit the He + H$ surface (we use only 155 nodes for the H + H2 surface) obtained a standard deviation of 0.57 kcal mol-1. This good agreement of the RMC-spline method with the YE function suggests that the spline fit may be accurate even when the cusp in the potential is present at 7 = 60”. To demonstrate that this is indeed the case, contour maps of the YL potential and VSPLINE are shown in fig. 1, for contours -105, -90, -75, -60 and -45 kcal mol- 1. It can be seen from fig_ 1 that the cusp at 8 = 45O is accurately reproduced and that the contour maps are essentially superimposable. Thus we find that a spline-fitted RMC function is capable
Volume 48, number 1
CHEMICAL PHYSICS LETTERS
15 May 1977
represented by a rotated Morse curve. This will not be the case if a potential well is present aIong a reaction coordinate containing a barrier; such systems are, however, rare. When they occur, they may be fitted by adding an additiona gaussian function to simulate the well, as was done by Kwei 1181 for the collinear K-Cl-Na system. Other difficult cases may be imagined, but the methdd outIined herein shouId be applicable to the vast majority of ground and excited state triatomic potential surfaces. The extension of the method to include a direct tit of ABC ab initio potential surface data, as well as dynamics on the fitted surface, is in progress.
References [I]
!2] [ 31 [4] Fig. 1. (a) Contour map for Yates-Lester potential surface at 7 = 60”, showing the cusp at 0 = 4S0 _Contour units: kcal mol-t . (b) Contour map for RMC-cubic spline potential sxface at y = 60’. The cusp is reproduced and the contours are very cIose to those in (a).
of fitting the cusp, contrary to the assumption of Bowman and Kuppermann [4] _ The RMC spline method given here is computationally attractive due to the small number of nodes (155 in this case) needed to evaluate the potential. (Due to the symmetry of the H, surface, only 16 X 3, or 48 of these nodes are unique. This fact can be used to reduce the interval search time.) Relative CPU times for VSPLINE : YL potential function evaluations are approximately 1.3 : 1 .O. For comparison, Sathyamurthy and Raff 183 report trajectories using the 3375 node spline grid to require a factor of 5 more CPU time than trajectories using the corresponding semi-empirical potential function (the potential and its three derivatives are evaluated at each point of the trajectory). The RMC spline method should be capable of fitting any triatomic system with good accuracy, provided that at any discrete -y the system can indeed be
[S]
[6] [7] [8] [9] [IO] [ 11) [ 121 [ 131 [ 141 [ 15 ] [ 161 [ 171 [ 181
hf. Karplus, in: The world ofquantum chemistry, eds.
R. Daudel and B. Pullman (Reidel, Dordrecht, 1974) p. 102. F.T. Wall and R.N. Porter, J. Chem. Phys. 36 (1962) 3256. J.N.L. Connor, W. Jakubetz and J. Xfanz, XfoL Phyr 29 (1975) 347. J.M. Bowman and A. Kuppermann, Chem. Phys. Letters 34 (1975) 523. S.K. Gray acd J.S. Wright, J. Chem. Phyr, to be published. I. Shavitt, R&I. Stevens, F.L. hIinn and hf. Karplus, 5. Chem. Phys. 49 (1968) 2700; I. Shavitt, J. Chem. Phys. 49 (1968) 4048. D.R. McLaughlin and D.L. Thompson, I. Chem. Phys: 59 (1973) 4393. N. Sathyamurthy and L-M. Raff, J. Chem. Phyr 63 (1975) 464. N. Sathyanmrthy, G.E. Kellerhals and L.&I. Raff, J. Chem. Phys. 64 (1976) 2259. X. Chapuisat and Y. Jean, Topics Current Chem. 68 (1977) 1. C. Leforestier, Thi%e de 3gme cycle, Facultk des Sciences, Orsay, France (1976). A.C. Yates and W. Lester Jr. , Cbem. Phys. Letters 24 (1974) 305. LG. Csizmadia, J.C. Polanyi, AC. Roach and W.H. Wong, Can. J. Chem. +7 (1969) 4097_ G.D. Carney and R.N. Porter, J. Chem. Phys. 60 (L9i4) 4251. G.D. Camey and R.N. Porter, J. Chem. Phys 65 (1976) 3547. B. Liu, J. Chem. Phys 58 (1973) 1925. R.N. Porter and hl. Karplus, J. Chem. Phys 40 (1964) 1105. G. Kwei. J. Chem. Phys. 58 (1973) 1722.
CHEMICAL PHYSICS LETTERS
15 Lfay 1977
A SPLINE-FITTED POTENTIAL SURFACE FOR BENT TRIATOMIC MOLECULES Stephen
K. GRAY, James S. WRIGHT $
Department of Chemistry, Carleton University, Ottawa, Canada KIS iB6
and Xavier CHAPLJJSAT Laboratoire de Chimie Thebrique $$, Universitt de Paris-Sud, Centre d’Orsay. 91405 Orsay, France Received 24 January 1977
The rotated-Morse curve cubic spliie method developed previously is extended to bent triatomic molecules by usinS 25 cubic splmes. The Yates-Lester potential for bent Hs is shown to be accurately fitted over the entire raqe of the three internal coordinates, with a standard deviation of less than 1 kcal mol-’ . The spline method compares favorably in computa-
tional speed with the analytical potentiaL
1. Introduction A general procedure to interpolate ab initio potential energy surface data is needed for studies in molecular dynamics [l] . Wall and Porter [2] developed a rotated Morse curve @MC) representation for producing collinear triatomic potential surfaces. Connor et al. 131 and Bowman and Kuppermann [4] expressed the RMC parameters as cubic spline functions, and by so doing obtained smooth potential surfaces with desired reaction path and saddle point characteristics. Gray and Wright [5] showed that the PMC-cubic spline method could give an accurate fit to an ab initio potential surface for collinear Hs, and that classical trajectories on the spline-fitted surface were in good agreement with those obtained by using the analytical SSMK function [6]. Gray and Wright [5] also found the trajectory computation time to be comparable for splme-fitted and analytical surfaces, thus demonstrating the attractiveness of the method for trajectory calculations involving only two internal coordinates RA, and
ordinates were deveioped by. lMcLaugh!in and Thompson [7] and extended by Sathyamurthy et ai. [S,9] _ En each case purely numerical grids along the internal coordinates were used. The 15 X 15 X 15 grid chosen by Sathyamurthy and Raff [8] to describe the systems (He + Hz) and (D -I- HCI) by 3D cubic spline interpotation proved to be relatively inefficient, both with respect to accuracy and computational speed of the classical trajectories- Some purely numerical studies using cardinal splines have also been reported [IO,1 I]. In this paper we outline an extension of the co& linear RMC-cubic spline method given previously [3--51 and apply it to the case of three internal coordinates [12] potentiaI RAB, RBC, and R,c- The Yates-Lester surface for bent H, is used as an ihustrative exampfe. A previous RMC potential surface for the Hz system was used by Csizmadia et al. [ 131 and included nonlinear geometries; however, no details were given and a large standard deviation was reported. The use of generahzed Morse functions [ 14,l S] to interpolate potentiaf surface data is also relevant to the present work.
RBC-
Sphne-fitted
potential
surfaces for three internal
co2. Method of calculation
$ To whom correspondence should be addressed. $* Associated with the CNRS (ERA no. 549).
For the three-particle
system
A... B...C let y be the
CHEMICALPHYSICSLETTERS
Volume48, number 1
ABC angle. The three internal coordinates are then y, RAB, andRBC, and the allowed range of y is from 180” (Linear) to 60°. For angles much smaller than 60° any RMC surface must become pathological [ 151; this problem can be avoided by redefining the internal coordinates so that 7 > 60” aIways. This is always possible since one of the interior angles of triangle ABC is necessarily greater than or equal to 60’. For our purposes here it is sufficient to consider just geometries contained within the range 60° G y < 180°. A direct extension of the collinear RMCcubic spline potential [3--51 is outhned beIow. At each specific 7, an (RAB, RBC) plot is drawn using orthogonal RAB and RBc axes, exactly as in the co&near case. Now 6 is the swing angle for the rotated Morse curve, defined as the clockwise angle of rotation about the swing point (RAB = RBC = 5 A) from the reference line RBc = 5 A. The potentiaI energy is: F(R,B,RBc,R_&= =D@,
F(e, Y, Z)
y)Y)Icl - exp[P@,
yXZ - Z,(e,
face were given by Porter and Karplus [ 171. These surfaces are smooth at y = 180° and 120° but show discontinuities in the derivatives (kinks in the potential) at +y= 60”. As before [S] ,3 I vahres of 0 were used ranging from 0 to 90’. The spacing of nodes was every 5’ from O-35” and every lo from 38-45O (and similarly for 46-900); this choice of nodes gave reasonable incremeni in the potential energy. It was found that only 5 values of y were needed to accurately fit the surface: 60”, ZZO”, 180”, 240”, and 300”. The Iatter two nodes are identical by symmetry to the first two and are used to give the proper curvature to the spline function at y = 180°; a similar argument was used previously [5 ] at e = 45O. The resulting 5 X 31 arrays for o(e, y), p(f7, y) and Z,(0,7) were fitted by 2D natural cubic spline functions ([7] , and references therein)_
3. Results and discussion y))l I2 - 11 ,
where tan 8 = (R& - RBC)/(RiB - RAB), 0 G 0 < 90”, (RiBI R&) = origin of (8, I) coordinate system (swing point), and C= [(RLB - RAB)2 f (REc - R,C)2] ‘j2. The Morse parameters D, P and Zeq, which in the collinear case were functions only of 6, now become func-
tions of 8 and -y. For various discrete values of y we obtain the specific values for D, ~3and Zes in an identical manner to that used for the cohinear geometry [S] . Briefly, we choose a discrete set of f3 and scan down rays of constant 0 to find Zes and D; the force constant k is obtained numerically at Zes and, from this, p = (k/2D)‘12. When this process is completed we have the required two-diiensional arraysD(B,7),Z,_&,?) and P@,y). We now fit these arrays to two-dimensional naturaI cubic splines using the
method given by McLaughIin and Thompson (71. To fuhy determine the sphne functions it is necessary to specify either the first derivativespr the second derivatives along the grid boundaries in e and 7. For simplicity, we set the second derivatives equal to zero along these grid boundaries. AS a test for the method we chose the YatesLester [12] analytical surface for H3, which is based on the collinear SCFCI calculation of Liu [ 16 ] and uses a Porter-Karplus [17] formulation to ahow for non-linear geometries. Contour maps for a similar sur156
15 May 1977
2000 points were generated randomly in the range 0.5
A comparison of the Yates-Lester (YL) potential and spline potential (VSPLINE) at these points showed a standard deviation of 0.70 kcal mol-I. To investigate the %hemicaIIy interesting” region of collision energies below about 2 eV, we discarded points leading to a YL energy above -55 kcal mol-r. The remaining 685 points showed a standard deviation of 0.38 kcai mol-l. This can be considered to be very good agreement, particularly since only 5 nodes in y were used. For comparison, the standard deviation of Csizmadia et al. [ 131 on the Hz surface was reported to be 6.65 kcal mol-l. Sathyamurthy and Raff [8], using 3375 nodal points to fit the He + H$ surface (we use only 155 nodes for the H + H2 surface) obtained a standard deviation of 0.57 kcal mol-1. This good agreement of the RMC-spline method with the YE function suggests that the spline fit may be accurate even when the cusp in the potential is present at 7 = 60”. To demonstrate that this is indeed the case, contour maps of the YL potential and VSPLINE are shown in fig. 1, for contours -105, -90, -75, -60 and -45 kcal mol- 1. It can be seen from fig_ 1 that the cusp at 8 = 45O is accurately reproduced and that the contour maps are essentially superimposable. Thus we find that a spline-fitted RMC function is capable
Volume 48, number 1
CHEMICAL PHYSICS LETTERS
15 May 1977
represented by a rotated Morse curve. This will not be the case if a potential well is present aIong a reaction coordinate containing a barrier; such systems are, however, rare. When they occur, they may be fitted by adding an additiona gaussian function to simulate the well, as was done by Kwei 1181 for the collinear K-Cl-Na system. Other difficult cases may be imagined, but the methdd outIined herein shouId be applicable to the vast majority of ground and excited state triatomic potential surfaces. The extension of the method to include a direct tit of ABC ab initio potential surface data, as well as dynamics on the fitted surface, is in progress.
References [I]
!2] [ 31 [4] Fig. 1. (a) Contour map for Yates-Lester potential surface at 7 = 60”, showing the cusp at 0 = 4S0 _Contour units: kcal mol-t . (b) Contour map for RMC-cubic spline potential sxface at y = 60’. The cusp is reproduced and the contours are very cIose to those in (a).
of fitting the cusp, contrary to the assumption of Bowman and Kuppermann [4] _ The RMC spline method given here is computationally attractive due to the small number of nodes (155 in this case) needed to evaluate the potential. (Due to the symmetry of the H, surface, only 16 X 3, or 48 of these nodes are unique. This fact can be used to reduce the interval search time.) Relative CPU times for VSPLINE : YL potential function evaluations are approximately 1.3 : 1 .O. For comparison, Sathyamurthy and Raff 183 report trajectories using the 3375 node spline grid to require a factor of 5 more CPU time than trajectories using the corresponding semi-empirical potential function (the potential and its three derivatives are evaluated at each point of the trajectory). The RMC spline method should be capable of fitting any triatomic system with good accuracy, provided that at any discrete -y the system can indeed be
[S]
[6] [7] [8] [9] [IO] [ 11) [ 121 [ 131 [ 141 [ 15 ] [ 161 [ 171 [ 181
hf. Karplus, in: The world ofquantum chemistry, eds.
R. Daudel and B. Pullman (Reidel, Dordrecht, 1974) p. 102. F.T. Wall and R.N. Porter, J. Chem. Phys. 36 (1962) 3256. J.N.L. Connor, W. Jakubetz and J. Xfanz, XfoL Phyr 29 (1975) 347. J.M. Bowman and A. Kuppermann, Chem. Phys. Letters 34 (1975) 523. S.K. Gray acd J.S. Wright, J. Chem. Phyr, to be published. I. Shavitt, R&I. Stevens, F.L. hIinn and hf. Karplus, 5. Chem. Phys. 49 (1968) 2700; I. Shavitt, J. Chem. Phys. 49 (1968) 4048. D.R. McLaughlin and D.L. Thompson, I. Chem. Phys: 59 (1973) 4393. N. Sathyamurthy and L-M. Raff, J. Chem. Phyr 63 (1975) 464. N. Sathyanmrthy, G.E. Kellerhals and L.&I. Raff, J. Chem. Phys. 64 (1976) 2259. X. Chapuisat and Y. Jean, Topics Current Chem. 68 (1977) 1. C. Leforestier, Thi%e de 3gme cycle, Facultk des Sciences, Orsay, France (1976). A.C. Yates and W. Lester Jr. , Cbem. Phys. Letters 24 (1974) 305. LG. Csizmadia, J.C. Polanyi, AC. Roach and W.H. Wong, Can. J. Chem. +7 (1969) 4097_ G.D. Carney and R.N. Porter, J. Chem. Phys. 60 (L9i4) 4251. G.D. Camey and R.N. Porter, J. Chem. Phys 65 (1976) 3547. B. Liu, J. Chem. Phys 58 (1973) 1925. R.N. Porter and hl. Karplus, J. Chem. Phys 40 (1964) 1105. G. Kwei. J. Chem. Phys. 58 (1973) 1722.