- Email: [email protected]
Method for semiclassical calculation of grid of eigenvalues
Volume
126, number
2
CHEMICAL
PHYSICS
LETTERS
2 May 1986
METHOD FOR SEMICLASSICAL CALCULATION OF GRID OF EIGENVALUES * Bobby G. SUMPTER’ Department
of Chemrstty,
Oklahoma
State University, Sttllwater,
OK 74078, USA
and D.W. NOID ’ Chemrstty Division, Oak Rtdge Natronal Luboratory, Oak Rrdge, TN 37831, USA and Department of Chemrstty, Unrversrty of Tennessee, KnoxutIle, TN 37916-1600, Received
17 December
1985; in final form 1 February
USA
1986
A method utilizing integration along curves on the Poincare surfaces of section is described for the semiclassical calculation of a grid of eigenvalues. The method is extremely fast and the eigenvalues calculated for anharmonically coupled oscillators are in excellent agreement with the exact quantum eigenvalues. In the example studied only 28 trajectories were used to obtain 42 eigenvalues.
1. Introduction Semiclassical eigenvalues of non-separable systems have been calculated by a variety of methods, which were recently reviewed [l] . Many of these methods utilize classical trajectories, while others use perturbative-iterative techniques. The path integral method has its origins in Einstein’s old quantum [2] theory and later in WKB-type semiclassical theory [ 3,4] . For this purpose, one computes, in the case of the quasiperiodic region, N-dependent canonical invariants (fi = 1)
0) Ci
where the 3 applies for a system of non-degenerate os-
Research sponsored by the U.S. Department of Energy under contract DE-AC05-84OR21400 with Martin Marietta Energy Systems, Inc. Research supported by Army Research Office. Present address: Institute for Defense Analysis, Alexandria, VA 22311. USA.
cillators, the Qk ,Pk are coordinates and their conjugate momenta, the ni are integers, and ci are topologitally independent closed paths. The paths Ci need not be along actual trajectories. The earliest path integral method was developed by Eastes and Marcus for a two-dimensional system [5 ] . The method involved locating the caustics of a trajectory. Semiclassically, the caustics separated the classically allowed from the classically non-allowed regions. Integrations along each of two sides of the four-sided caustic yielded the two independent path integrals. Using a procedure to interpolate data from three trajectories, the eigenvalues were determined by requiring both n 1 and n2 in eq. (1) to be integers. This method was successful for trajectories which gave rise to an obvious caustic but did not appear appropriate for complicated caustics. Noid et al. [6-l l] have developed a more general method in which the independent paths were taken as the Poincare surfaces of section. The method was applied successfully to a system of two dimensions, and has recently been applied to three-dimensional systems [9] and systems with higherorder resonances [ 1 l] . This method involves the numerical integration of a trajectory, and the behavior of the trajectory in a plane of 181
Volume 126, number 2
CHEMICAL PHYSICS LETTERS
phase space is examined. A section of phase space is defined by a plane, for example, Q1 = 0, and each time the trajectory crosses the plane in a given direction (PI > 0), the point corresponding to the value of Q2 and P2 was marked, thus generating a closed “curve” of P2 versus Q2. The same procedure is followed to generate the P, versus Q1 closed “curve”. The two “surfaces of section” are then used as paths to evaluate the phase integrals. In this paper, the method is modified so that a grid of eigenvalues is computed. This method is very easy to program and the computational time is relatively short. The calculated eigenvalues are in excellent agreement with the exact quantum eigenvalues. Software exists for fitting scattered data in three dimensions [ 121. Thus the method presented in this paper for calculating a grid of eigenvalues can be easily extended to three-dimensional systems. In section 2, the Hamiltonian and the method of the calculation are discussed. The results and a discussion are given in section 3.
2 May 1986
A grid of surfaces of section was computed by evaluating N trajectories as described above. The initial conditions were taken from (fi = 1) E, = wl(n; ++),I E2 = w2(n; t;),
(4)
where values of rry varied from 0 to 6 and n; varied from 0 to 5. The trajectories used to calculate the eigenvalues were those corresponding to the initial conditions (n ‘2= 0 to 6 and ny = 0,2,4,5). The quantum conditions are, using topologically independent paths cl and c2 in units of fi = 1:
W3
+$=f
(PIdQl +P2dQ2),
Cl
2W2
+$
=j-
(PI
dQl +P, de,).
(5)
CZ Then1 and n2 are integers when the energy is an eigenvalue . The curves on the surfaces of section Q1 = 0 and Q2 = 0 are topologically equivalent to the paths c2 and cl. Thus eq. (5) yields in units of R = 1:
2. Method A two-dimensional oscillator with no resonances is studied here by way of example. The model Hamiltonian is
where the terms in the square brackets describe the uncoupled normal coordinate Hamiltonian for two normal modes. The parameters used for the calculations are ~0: = 0.49, wi = 1.69, h = -0.10, and g = 0.10. The Hamiltonian yields trajectories which are “box like”. Hamilton’s equations for the Hamiltonian in eq. (2) were integrated using the program DEROOT [ 131. This routine integrates the trajectory and automatically determines points on the trajectory which are roots of an equation provided by the user. The root-finding part of DEROOT was set to return points when the equation
QlQ2=o
(3)
was satisfied. These points were saved for later use in evaluating the surfaces of section. 182
2n(n2 ++)=$P2
dQ,,
Q, = 0.
(6)
The points are first transformed to action-angle variables (Jr?, Wjo) given below with standard convention on the relation between the phase Wj and the sign of Pj [14] : 27rlVF= tanV1(tijQj/Pj),
J: = rr(Pf t o$Q~)/~j,(7)
where Wj is the zerothorder angular frequency of the ith coordinate. The points are next sorted in increasing order of WF. To the first point (IV:), 1 is added and the path integral in eq. (6) becomes (FJ= 1) 1+wy s
Ji”
dk’f” = 27r(nj t 4).
(8)
WY The phase integral in eq. (8) is evaluated using a quadrature routine [ 151. Twenty points are generally sufficient for accuracy to six digits. This reduces the total time needed to generate a value for a nj in eq. (8) to about two cpu seconds. Each eigenvalue is then found
CHEMICAL PHYSICS LETTERS
Volume 126, number 2
Table 1 Comparison of semiclassical eigenvalues for the Hamiltonian given in eq. (2)
0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 3 3 3 3 3 3 4 4 4 4 4 4 4 5 5 5 5 5 5 5
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6
0.9955 1.6870 2.3750 3.0595 3.7404 4.4176 5.0910 2.2782 2.9584 3.6348 4.3071 4.9153 5.6393 6.2989 3.5480 4.2164 4.8803 5.5397 6.1944 6.8443 7.4891 4.8045 5.4601 6.1108 6.7564 1.3961 8.0315 8.6607 6.0468 6.6889 1.3255 7.9564 8.5814 9.2001 9.8125 1.2142 7.9019 8.5235 9.1387 9.7472 10.3489 10.9434
0 0 0 0.003 0 0 0.002 0.004 0 0.003 0.005 0.008 0.014 0.22 0.003 0.005 0.008 0.013 0.021 0.035 0.055 0.004 0.007 0.015 0.027 0.043 0.070 0.109 0.008 0.016 0.029 0.050 0.083 0.130 0.191 0.017 0.029 0.053 0.090 0.145 0.23 0.35
0.9955 1.6870 2.3750 3.0595 3.7404 4.4175 5.0908 2.2781 2.9584 3.6341 4.3069 4.9749 5.6385 6.2974 3.5480 4.2162 4.8800 5.5391 6.1931 6.8418 7.4847 4.8043 5.4598 6.1100 6.7545 7.3935 8.0258 8.6511 6.0463 6.6881 1.3225 7.9506 8.5742 9.1880 9.7929 7.2731 7.8999 8.5 194 9.1313 9.7319 10.3250 10.9085
0.9955 1.6870 2.3750 3.0596 3.7404 4.4176 5.0909 2.2781 2.9584 3.6341 4.3069 4.9149 5.6385 6.2975 3.5479 4.2162 4.8199 5.5390 6.1931 6.8419 7.4850 4.8043 5.4591 6.1099 6.7546 7.3935 8.0259 8.6513 6.0463 6.6878 7.3234 7.9524 8.5113 9.1882 9.7932 7.2730 7.8996 8.5190 9.1305 9.7331 10.3253 10.9053
0 0 0 0.003 0 0.002 0.002 0 0 0 0 0 0 0.0016 0.003 0 0.002 0.002 0 0.002 0.004 0 0.002 0.002 0.001 0 0.001 0.002 0 0.004 0.01 0.02 0.001 0.002 0.003 0.001 0.004 0.005 0.009 0.01 0.003 0.03
2 May 1986
by using a recent code for the Sheppard [ 161 method of constructing a surface from scattered data. Twentyeight trajectories were needed to construct a surface in which the forty-two interpolated eigenvalues in table 1 were accurate to six digits. This allows the simultaneous calculation of a grid of eigenvalues in a very short time, and clearly required less than one trajectory per eigenvalue. (Almost two eigenvalues per trajectory were obtained.)
3. Results and discussion A grid of eigenvalues is calculated for the Hamiltonian in eq. (2). The results are given in table 1 and are compared to calculations using the BirkoffGustavson normal form [ 171 and the exact quantum eigenvalues. As can be seen, the eigenvalues are in excellent agreement with the exact quantum eigenvalues. The largest deviation (0.03) occurs for the last eigenvalue which was at the edge of the grid. The present semiclassical method provides more accurate eigenvalues for the Hamiltonian in eq. (2) than the eigenvalues calculated using Birkoff-Gustavson normal form. The biggest advantage of this technique is that a grid of eigenvalues can be simultaneously determined in a relatively short computational time. The eigenvalues given in table 1 were calculated from twenty-eight trajectories. The total number of trajectories required approximately 55 s on an IBM 3033 and the interpolation took 4 s. If less accuracy is desired, the time can be reduced by using fewer trajectories to calculate the eigenvalues. For example, a calculation using twenty-one trajectories required approximately 40 s and the eigenvalues were accurate to five digits. The only restriction of the method is that the motion must be quasiperiodic. The chaotic trajectories of an N-dimensional system do not have N isolating integrals of the motion. Good action-angle variables no longer exist, and accordingly, these trajectories cannot be quantized with the present method.
a) Calculated by Swimm and Delos [ 171.
b) AE% = (IE- EQMI/EQM) X 100, where E is either EBGNF or ESC. c) Values calculated using the method described in the text. d) Exact quantum mechanical eigenvalues.
Acknowledgement
One of us (BGS) is pleased to acknowledge that he was a Visiting Scientist at Oak Ridge National Laboratory while working on this paper. 183
Volume 126, number 2
CHEMICAL PHYSICS LETTERS
References 111D.W. Noid, M.L. Koszykowski and R.A. Marcus. Ann. Rev. Phys..Chem. 32 (1981) 267. (21 A. Einstein, Verh. Deutsch. PhysIk. Ges. 19 (1917) 82. (31 J.B. Keller, Ann. Phys. 4 (1958) 1801. [4]R.A. Marcus, Discussions Faraday Sot. 55 (1973) 34. [S] W.A. Eastes and R.A. Marcus, J. Chem. Phys. 61 (1974) 4301. [6] D.W. Noid and R.A. Marcus, J. Chem. Phys. 62 (1975) 2119. (71D.W. Noid and R.A. Marcus, J. Chem. Phys. 67 (1977) 559. [8] D.W. Noid, M.L. Koszykowskiand R.A. Marcus, J. Chem. Phys. 71 (1979) 2864. [ 91 D.W. Noid, M.L. Koszykowski and R.A. Marcus, J. Chem. Phys. 73 (1980) 3911.
184
2 May 1986
[lo] D.W. Noid and M.L. Koszykowski, Chem. Phys. Letters 73 (1980) 114. [ 111 D.W. Noid and R.A. Marcus, J. Chem. Phys., _ submitted for publication. fl21 R. Franke, private communication. 1131L.F. Shampine and M.K. Gordon, Computer solutions of ordinary diiferential equations: the initial value problem (Freeman, San Francisco, 1975). [141 M. Born, Mechanics of the atom (Unger, New York, 1960). iI51VA. Dixon, in: Solutions for numerical mathematics, ed. DJ. Evans (Academic Press, New York, 1974) pp. 105-113; P.E. GiII and G.F. Miller, Comput. J. 15 (1972) 80. (161 R. Franke, Math. Comp. 38 (1982) 181. [171 R.T. Swimm and J.B. Delos, J. Chem. Phys. 71 (1979) 1706.
126, number
2
CHEMICAL
PHYSICS
LETTERS
2 May 1986
METHOD FOR SEMICLASSICAL CALCULATION OF GRID OF EIGENVALUES * Bobby G. SUMPTER’ Department
of Chemrstty,
Oklahoma
State University, Sttllwater,
OK 74078, USA
and D.W. NOID ’ Chemrstty Division, Oak Rtdge Natronal Luboratory, Oak Rrdge, TN 37831, USA and Department of Chemrstty, Unrversrty of Tennessee, KnoxutIle, TN 37916-1600, Received
17 December
1985; in final form 1 February
USA
1986
A method utilizing integration along curves on the Poincare surfaces of section is described for the semiclassical calculation of a grid of eigenvalues. The method is extremely fast and the eigenvalues calculated for anharmonically coupled oscillators are in excellent agreement with the exact quantum eigenvalues. In the example studied only 28 trajectories were used to obtain 42 eigenvalues.
1. Introduction Semiclassical eigenvalues of non-separable systems have been calculated by a variety of methods, which were recently reviewed [l] . Many of these methods utilize classical trajectories, while others use perturbative-iterative techniques. The path integral method has its origins in Einstein’s old quantum [2] theory and later in WKB-type semiclassical theory [ 3,4] . For this purpose, one computes, in the case of the quasiperiodic region, N-dependent canonical invariants (fi = 1)
0) Ci
where the 3 applies for a system of non-degenerate os-
Research sponsored by the U.S. Department of Energy under contract DE-AC05-84OR21400 with Martin Marietta Energy Systems, Inc. Research supported by Army Research Office. Present address: Institute for Defense Analysis, Alexandria, VA 22311. USA.
cillators, the Qk ,Pk are coordinates and their conjugate momenta, the ni are integers, and ci are topologitally independent closed paths. The paths Ci need not be along actual trajectories. The earliest path integral method was developed by Eastes and Marcus for a two-dimensional system [5 ] . The method involved locating the caustics of a trajectory. Semiclassically, the caustics separated the classically allowed from the classically non-allowed regions. Integrations along each of two sides of the four-sided caustic yielded the two independent path integrals. Using a procedure to interpolate data from three trajectories, the eigenvalues were determined by requiring both n 1 and n2 in eq. (1) to be integers. This method was successful for trajectories which gave rise to an obvious caustic but did not appear appropriate for complicated caustics. Noid et al. [6-l l] have developed a more general method in which the independent paths were taken as the Poincare surfaces of section. The method was applied successfully to a system of two dimensions, and has recently been applied to three-dimensional systems [9] and systems with higherorder resonances [ 1 l] . This method involves the numerical integration of a trajectory, and the behavior of the trajectory in a plane of 181
Volume 126, number 2
CHEMICAL PHYSICS LETTERS
phase space is examined. A section of phase space is defined by a plane, for example, Q1 = 0, and each time the trajectory crosses the plane in a given direction (PI > 0), the point corresponding to the value of Q2 and P2 was marked, thus generating a closed “curve” of P2 versus Q2. The same procedure is followed to generate the P, versus Q1 closed “curve”. The two “surfaces of section” are then used as paths to evaluate the phase integrals. In this paper, the method is modified so that a grid of eigenvalues is computed. This method is very easy to program and the computational time is relatively short. The calculated eigenvalues are in excellent agreement with the exact quantum eigenvalues. Software exists for fitting scattered data in three dimensions [ 121. Thus the method presented in this paper for calculating a grid of eigenvalues can be easily extended to three-dimensional systems. In section 2, the Hamiltonian and the method of the calculation are discussed. The results and a discussion are given in section 3.
2 May 1986
A grid of surfaces of section was computed by evaluating N trajectories as described above. The initial conditions were taken from (fi = 1) E, = wl(n; ++),I E2 = w2(n; t;),
(4)
where values of rry varied from 0 to 6 and n; varied from 0 to 5. The trajectories used to calculate the eigenvalues were those corresponding to the initial conditions (n ‘2= 0 to 6 and ny = 0,2,4,5). The quantum conditions are, using topologically independent paths cl and c2 in units of fi = 1:
W3
+$=f
(PIdQl +P2dQ2),
Cl
2W2
+$
=j-
(PI
dQl +P, de,).
(5)
CZ Then1 and n2 are integers when the energy is an eigenvalue . The curves on the surfaces of section Q1 = 0 and Q2 = 0 are topologically equivalent to the paths c2 and cl. Thus eq. (5) yields in units of R = 1:
2. Method A two-dimensional oscillator with no resonances is studied here by way of example. The model Hamiltonian is
where the terms in the square brackets describe the uncoupled normal coordinate Hamiltonian for two normal modes. The parameters used for the calculations are ~0: = 0.49, wi = 1.69, h = -0.10, and g = 0.10. The Hamiltonian yields trajectories which are “box like”. Hamilton’s equations for the Hamiltonian in eq. (2) were integrated using the program DEROOT [ 131. This routine integrates the trajectory and automatically determines points on the trajectory which are roots of an equation provided by the user. The root-finding part of DEROOT was set to return points when the equation
QlQ2=o
(3)
was satisfied. These points were saved for later use in evaluating the surfaces of section. 182
2n(n2 ++)=$P2
dQ,,
Q, = 0.
(6)
The points are first transformed to action-angle variables (Jr?, Wjo) given below with standard convention on the relation between the phase Wj and the sign of Pj [14] : 27rlVF= tanV1(tijQj/Pj),
J: = rr(Pf t o$Q~)/~j,(7)
where Wj is the zerothorder angular frequency of the ith coordinate. The points are next sorted in increasing order of WF. To the first point (IV:), 1 is added and the path integral in eq. (6) becomes (FJ= 1) 1+wy s
Ji”
dk’f” = 27r(nj t 4).
(8)
WY The phase integral in eq. (8) is evaluated using a quadrature routine [ 151. Twenty points are generally sufficient for accuracy to six digits. This reduces the total time needed to generate a value for a nj in eq. (8) to about two cpu seconds. Each eigenvalue is then found
CHEMICAL PHYSICS LETTERS
Volume 126, number 2
Table 1 Comparison of semiclassical eigenvalues for the Hamiltonian given in eq. (2)
0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 3 3 3 3 3 3 4 4 4 4 4 4 4 5 5 5 5 5 5 5
0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6
0.9955 1.6870 2.3750 3.0595 3.7404 4.4176 5.0910 2.2782 2.9584 3.6348 4.3071 4.9153 5.6393 6.2989 3.5480 4.2164 4.8803 5.5397 6.1944 6.8443 7.4891 4.8045 5.4601 6.1108 6.7564 1.3961 8.0315 8.6607 6.0468 6.6889 1.3255 7.9564 8.5814 9.2001 9.8125 1.2142 7.9019 8.5235 9.1387 9.7472 10.3489 10.9434
0 0 0 0.003 0 0 0.002 0.004 0 0.003 0.005 0.008 0.014 0.22 0.003 0.005 0.008 0.013 0.021 0.035 0.055 0.004 0.007 0.015 0.027 0.043 0.070 0.109 0.008 0.016 0.029 0.050 0.083 0.130 0.191 0.017 0.029 0.053 0.090 0.145 0.23 0.35
0.9955 1.6870 2.3750 3.0595 3.7404 4.4175 5.0908 2.2781 2.9584 3.6341 4.3069 4.9749 5.6385 6.2974 3.5480 4.2162 4.8800 5.5391 6.1931 6.8418 7.4847 4.8043 5.4598 6.1100 6.7545 7.3935 8.0258 8.6511 6.0463 6.6881 1.3225 7.9506 8.5742 9.1880 9.7929 7.2731 7.8999 8.5 194 9.1313 9.7319 10.3250 10.9085
0.9955 1.6870 2.3750 3.0596 3.7404 4.4176 5.0909 2.2781 2.9584 3.6341 4.3069 4.9149 5.6385 6.2975 3.5479 4.2162 4.8199 5.5390 6.1931 6.8419 7.4850 4.8043 5.4591 6.1099 6.7546 7.3935 8.0259 8.6513 6.0463 6.6878 7.3234 7.9524 8.5113 9.1882 9.7932 7.2730 7.8996 8.5190 9.1305 9.7331 10.3253 10.9053
0 0 0 0.003 0 0.002 0.002 0 0 0 0 0 0 0.0016 0.003 0 0.002 0.002 0 0.002 0.004 0 0.002 0.002 0.001 0 0.001 0.002 0 0.004 0.01 0.02 0.001 0.002 0.003 0.001 0.004 0.005 0.009 0.01 0.003 0.03
2 May 1986
by using a recent code for the Sheppard [ 161 method of constructing a surface from scattered data. Twentyeight trajectories were needed to construct a surface in which the forty-two interpolated eigenvalues in table 1 were accurate to six digits. This allows the simultaneous calculation of a grid of eigenvalues in a very short time, and clearly required less than one trajectory per eigenvalue. (Almost two eigenvalues per trajectory were obtained.)
3. Results and discussion A grid of eigenvalues is calculated for the Hamiltonian in eq. (2). The results are given in table 1 and are compared to calculations using the BirkoffGustavson normal form [ 171 and the exact quantum eigenvalues. As can be seen, the eigenvalues are in excellent agreement with the exact quantum eigenvalues. The largest deviation (0.03) occurs for the last eigenvalue which was at the edge of the grid. The present semiclassical method provides more accurate eigenvalues for the Hamiltonian in eq. (2) than the eigenvalues calculated using Birkoff-Gustavson normal form. The biggest advantage of this technique is that a grid of eigenvalues can be simultaneously determined in a relatively short computational time. The eigenvalues given in table 1 were calculated from twenty-eight trajectories. The total number of trajectories required approximately 55 s on an IBM 3033 and the interpolation took 4 s. If less accuracy is desired, the time can be reduced by using fewer trajectories to calculate the eigenvalues. For example, a calculation using twenty-one trajectories required approximately 40 s and the eigenvalues were accurate to five digits. The only restriction of the method is that the motion must be quasiperiodic. The chaotic trajectories of an N-dimensional system do not have N isolating integrals of the motion. Good action-angle variables no longer exist, and accordingly, these trajectories cannot be quantized with the present method.
a) Calculated by Swimm and Delos [ 171.
b) AE% = (IE- EQMI/EQM) X 100, where E is either EBGNF or ESC. c) Values calculated using the method described in the text. d) Exact quantum mechanical eigenvalues.
Acknowledgement
One of us (BGS) is pleased to acknowledge that he was a Visiting Scientist at Oak Ridge National Laboratory while working on this paper. 183
Volume 126, number 2
CHEMICAL PHYSICS LETTERS
References 111D.W. Noid, M.L. Koszykowski and R.A. Marcus. Ann. Rev. Phys..Chem. 32 (1981) 267. (21 A. Einstein, Verh. Deutsch. PhysIk. Ges. 19 (1917) 82. (31 J.B. Keller, Ann. Phys. 4 (1958) 1801. [4]R.A. Marcus, Discussions Faraday Sot. 55 (1973) 34. [S] W.A. Eastes and R.A. Marcus, J. Chem. Phys. 61 (1974) 4301. [6] D.W. Noid and R.A. Marcus, J. Chem. Phys. 62 (1975) 2119. (71D.W. Noid and R.A. Marcus, J. Chem. Phys. 67 (1977) 559. [8] D.W. Noid, M.L. Koszykowskiand R.A. Marcus, J. Chem. Phys. 71 (1979) 2864. [ 91 D.W. Noid, M.L. Koszykowski and R.A. Marcus, J. Chem. Phys. 73 (1980) 3911.
184
2 May 1986
[lo] D.W. Noid and M.L. Koszykowski, Chem. Phys. Letters 73 (1980) 114. [ 111 D.W. Noid and R.A. Marcus, J. Chem. Phys., _ submitted for publication. fl21 R. Franke, private communication. 1131L.F. Shampine and M.K. Gordon, Computer solutions of ordinary diiferential equations: the initial value problem (Freeman, San Francisco, 1975). [141 M. Born, Mechanics of the atom (Unger, New York, 1960). iI51VA. Dixon, in: Solutions for numerical mathematics, ed. DJ. Evans (Academic Press, New York, 1974) pp. 105-113; P.E. GiII and G.F. Miller, Comput. J. 15 (1972) 80. (161 R. Franke, Math. Comp. 38 (1982) 181. [171 R.T. Swimm and J.B. Delos, J. Chem. Phys. 71 (1979) 1706.