Volume 153, number
L
CHEMICAL
PHYSICS
LETTERS
2 December
1988
THE COLLOCATION METHOD FOR BOUND SOLUTIONS OF THE SCHRODINGER EQUATION Weitao YANG and Andrew C. PEET Department of Chemisrry, University of California, and Mater& and Chemical Sciences Division, Lawrence Rerkeiey Lahoralory, Berkeley. CA 94 720, USA Received 5 August 1988; in final form 20 September
The collocation method for obtaimng require the evaluation of integrals and used recently, but has the advantage of on two Morse oscillator problems and conjunction with a distributed Gaussian
1988
the bound solutions of the Schrb;dinger equation is investigated. The technique dots not is very simple to implement. It is closely connected with other pointwise representations requiring less effort to construct the algebraic eigenvalue equations. The method is tested found to give results which are as accurate as the conventional variational approach. In basis the collocation method is shown to he capable of describing highly excited states.
1. Introduction The conventional way to obtain the approximate eigcnsolutions of the Schriidinger equation is to expand the wavefunction as a linear combination of basis functions, then use the Rayleigh-Ritz variational principle to reduce the original differential equation to a set of algebraic eigenvalue equations. Recently there has been a surge of interest in alternative techniques for setting up algebraic eigenvalue equations. The idea is to replace the potential by a pointwise representation but to retain a basis set for defining the kinetic energy term. Such methods have been known for many years but have only recently been shown to provide a useful framework for solving real quantum mechanical problems [ l-161. In particular, Light and co-workers have had much success in developing the discrete variable representation (DVR) which they have applied to both scattering and bound state problems [l-4]. In addition, Friesner has constructed a pseudospectral method for the solution of the Hartree-Fock equations [ 5-7 ] and managed to achieve significant savings in computational effort. In this Letter we consider perhaps the simplest of the pointwise approximations, the collocation method. This technique is a classic one for solving differential equations and is widely used in engi98
neering problems [ 12,13 1. However, to our knowledge it has only once been applied to the calculation of molecular bound states. Unfortunately, this previous study produced results inferior to those of the variational basis set method [ 111. Here we investigate the properties of the collocation method more completely and find it capable of producing highly accurate results. We also suggest reasons for the lack of accuracy in the earlier work.
2. The collocation method In essence the collocation method [ 12,131 is as follows. A trial wavefunction Y is expanded as a linear combination of Nbasis functions, ,Xc;6,. Y is then taken to satisfy the Schriidinger equation at Npoints, rl, .. .rN, in space. This yields N algebraic equations which have the structure of a general eigenvalue problem and give the eigenvalues and coefficients. 2.1. The relation to conventional variational
methods To give some insight into the relation of the collocation method to conventional variational techniques, we present the different approaches in parallel using Dirac notation. The general aim of all the
Volume 153,number 1
CHEMICALPHYSICSLETTERS
methods is to select coefficients which optimise the trial wavefunction. The accuracy of the chosen function may be measured in terms of the residue function IR) = (H-E) 1Y}, where His the Hamiltonian and E the eigenvalue. The most popular technique is the Rayleigh-Ritz variational principle. The residue function IR) is converted into the functional ( !PlH-E I Y} and this is minimised with respect to variation of the c,. This approach yields the well known algebraic eigenvalue equations
,c,
((@,IH-EI$j,,-))Cj=O,
i=l,.-.N.
(1)
Another method, due to Galerkin [ 12,131, is to force IR) to be orthogonal to the basis functions used, (@iIR)=(q$lH-EI!P)=O,
i=l,...N.
(2)
Here, this condition gives exactly the same algebraic equations as the Rayleigh-Ritz method, eq. ( 1). However, the Galerkin derivation is more illuminating due to its close connection with the collocation method. In the collocation method 1R) is forced to be orthogonal to the N position vectors I r,) (rj]R)=(rillrl-EIY)=O, yielding the collocation f
((rjIH_E($j))Cj=O,
i=l,...N,
(3)
equations i=l,*..N.
(4)
j=l
We may now see the relation between the collocation method and the Rayleigh-Ritz principle. Multiplying the ith equation of eq. (4) by (@j Iri> w, and summing over i gives
2 December 1988
tion points {I,} should be chosen to make the N-point quadrature (5) as accurate as possible. For a basis set of global orthogonal functions, such as harmonic oscillator functions, the associated Gaussian quadrature scheme may be used to provide a reliable set of collocation points. For localised functions the choice is not as clear and no standard procedure exists for determining an optimum quadrature scheme. With distributed Gaussians we have found that the centres of the functions provide a satisfactory set of collocation points. The collocation method may be thought of as a fitting procedure in which the coefficients of a trial function I R) are selected to give an exact tit at the collocation points. Consequently, only the value of 1R) at these points is required. This is a great simplification over the variational methods which require 1R) to be evaluated at a sufficient number of points to allow integrals of the type (#, I R > to be converged. Another way of viewing this property is to note that the Hamiltonian enters into the collocation equations (4) as a matrix in the mixed representation of {r,} and {&}. Although the collocation method provides a particularly simple way of setting up the algebraic eigenvalue equations, it does lack some preferable features of the Rayleigh-Ritz-Galerkin formulation. The technique is not variational and so the approximate eigenvalues are no longer upper bounds of the exact ones. Furthermore, the Hamiltonian is represented by a non-symmetric matrix and so complex eigenvalues and non-orthogonal eigenvectors may occur. However, this should only be a formal difftculty as the collocation equations are equivalent to the Rayleigh-Ritz-Galerkin equations with an Npoint quadrature approximation to the integrals, a formulation which should be nearly symmetric. 2.2. The relation to other pointwise representations
which is an N-point quadrature approximation to eq. (2). Hence, we see that the collocation method is equivalent to the Rayleigh-Ritz principle with all the integrals evaluated by the same N-point quadrature (see also page 137 of ref. [ 121). It is interesting to note that the energy levels and wavefunctions given by eq. ( 5 ) are totally independent of the weights {wZ} used. From the above analysis, we see that the colloca-
To make the relation of the collocation method to other pointwise methods as transparent as possible we consider the case of a particle moving in one dimension. The collocation equations (4) now become ,$, (-
$
i=l,...N,
@Ytri) + V(ri>
4jCri)
-E$j(ri))c,=o,
(6) 99
Volume 153, number
CHEMICAL
1
where 4; (r, ) is the second derivative of @jevaluated at the collocation point ri. The quantity fi is the mass of the particle and V(r) is the potential function. Using the matrices R,=#Jr,),
Gil=-
$$‘(r,),
v, = Vr,) &,,
(7) we may write the collocation equations in the form (G+VR-ER)c=O.
(8)
The eigenvalues are then given by setting the determinant of the matrix in parentheses to zero. The connection with Friesner’s pseudospectral representation may now easily be obtained by using a property of determinants, IGR-‘+V-El]
]RI =O.
(9)
Assuming that the determinant of R is not zero then yields the pseudospectral prescription for the eigenvalues. From this analysis it is apparent that the collocation and pseudospectral methods will give identical results for the same basis set and points. We note that the pseudospectral approach of Orszag is the collocation method [ 15,161. Whereas the pseudospectral representation of Friesner may be considered as a rearrangement of the collocation equations, the discrete variable representation (DVR) is based upon a pointwise representation of the Rayleigb-Ritz-Gale&in equations ( 1) with an N-point quadrature approximation for the potential matrix. Due to the reliance of the method upon quadrature, the matrix corresponding to R in eq. (7) now contains weights and is given the form, T!,= ,,& @,(rI;). The Rayleigh-Ritz-Galerkin equations may then be written (approximately) as (K+T=VT-ES)c=O,
(10)
where K is the kinetic energy matrix and S is the overlap matrix, both evaluated exactly in the basis set representation, In all realistic studies presently reported, the DVR has been used in conjunction with a Gaussian quadrature scheme and the corresponding orthogonal functions. In this case the matrix T is orthogonal and S becomes the identity. Using these properties and 100
2 December1988
PHYSICSLETTERS
proceeding as before we obtain the condition for the eigenvalues ITT1 ITKT’SV-El1
ITI =O,
(11)
which, if the determinant of T is not zero, gives the DVR equations. Comparing eqs. (9) and (11) it is clear that the Gaussian version of the DVR differs from the pseudospectral representation only by the manner in which the kinetic energy matrix is set up. In the DVR, the kinetic energy matrix is evaluated exactly, whereas in the pseudospectral representation an Npoint quadrature approximation is implied. In general, the errors in the methods will be dominated by the approximation to the potential matrix and so the DVR is not likely to produce results which are significantly more accurate than pseudospectral technique. We again note that the pseudospectral and collocation approaches are of equivalent accuracy, see eqs. (8) and (9). A generalised DVR for treating orthogonal basis functions with non-Gaussian quadratures has also been proposed [ 31. The technique relies upon a quadrature scheme which is not of the type used in eq. ( 5) (the weights now depend upon two indices) and is not simpiy related to the collocation method. To our knowledge a form of the DVR which can deal with a non-orthogonal basis has not yet been proposed, We note that the collocation method may be applied to such basis sets (as long as a suitable set of collocation points may be chosen). Although the collocation, pseudospectral and DVR techniques may yield results of equivalent accuracy, the equations which must be solved are not the same. Different amounts of effort may, therefore, be required to set up and solve the respective equations. Both the DVR and pseudospectral approaches have the disadvantage that a matrix inverse must be computed in order to determine the kinetic energy matrix in the pointwise representation (although in the DVR this is merely obtained as a matrix transpose if the method is restricted to a set of orthogonal functions and a corresponding Gaussian quadrature). The collocation method does not require such an opcration and, in general, will provide the simplest way of setting up the algebraic eigenvalue equations. As these equations are invariably solved by using a “canned” routine, a very attractive feature of the
Volume 153, number 1
CHEMICAL PHYSICS LETTERS
collocation method is the simplicity with which it may be programmed up. The collocation equations are in the form of a general eigenvalue problem. One way to solve eq. (8) is first to use R-’ to convert it into Friesner’s pseudospectral equation ( 9 ) . However, more efficient methods exist for solving general eigenvalue problems. We used the routine F02BJF from the NAG library. A disadvantage of the collocation method is that the algebraic eigenvalue problem is inherently unsymmetric with the consequent disadvantages discussed in section 2.1. Similar lack of symmetry does not occur in the DVR and only occurs in the pseudospectral representation due to inaccuracies in the computation of the kinetic energy matrix, a problem which may be overcome by averaging the off-diagonal elements.
3. Test cases - the Morse oscillator 3.1. Harmonic
oscillator basis functions
First, we consider vibrational states which are described well by simple harmonic oscillator functions. This gives the natural choice of Gauss-Hermite points for the collocation. Such functions and points have been adopted in other pointwise approaches, allowing a close comparison of the methods [ 31. Also, the use of the collocation method with an orthogonal basis and the corresponding Gaussian points is known to provide an accurate framework for solv-
ing certain
2 December 1988
types of ordinary
differential
equations
1121. We use the Morse oscillator model of the (J=O) HF molecule considered by Light and co-workers [ 3,171. The Hamiltonian for this problem is H=
---
h2 d2 2,~ dr2
+D{exp[-2a(r-r,)]-2exp[-cu(r-r,)]}, (12) with r,= 1.75 a,, fl= 19/20 amu, 0=5.716 eV and a= 1.2121 a,‘. For this model problem, there are 23 bound states. The basis set of harmonic oscillator eigenfunctions are defined by the parameters r,-,=2.18 a0 and p= (pk/fi’) ‘14=4.6 a,y’, which were chosen to have their (approximately ) optimum values [ 3 1. We first consider the accuracy of the collocation method with a relatively small basis set (NE 12). The results for the first 8 levels of the Morse oscillator are given in column C of table 1. In columns A and B we give the values obtained from the full variational calculation [ 31 and the DVR of Light et al. [ 31, respectively. Comparison of the three sets of results clearly indicates that (for this problem) the variational, DVR and collocation methods are all of similar accuracy. We also note that the non-symmetric eigenvalue problem, which is inherent to the collocation method, presents no practical difficulty. Complex eigenvalues, if any, are found to occur only well above the set of states converged by the basis set. In column D of table I we give the values ob-
Table 1 Comparison of the errors in the eigenvalues (d,=E,(calc) -E,(exact) ) for Morse oscillator A using a simple harmonic oscillator basis. Columns A-C N= 12, A: variational, B: DVR with Gauss-Hem&e (GH) points and weights (ref. [ 31)) C: collocation with GH points; D, E N=25, D: collocation with GH points, E: collocation with equal spaced points between Rmin= I .O and J&,=3.3. Exact & (cV) - 5.4620 -4.9714 -4.5038 -4.0594 -3.6380 - 3.2397 -2.8645 -2.5123
A,,, N= 12
A,,, N=25
A
B
C
D
0.0000 0.0003 0.0014 0.0033 0.0095 0.0312 0.0689 0.1131
0.0000 0.0004 -0.0015 -0.0125 -0.0126 0.0294 0.0914 -0.0057
0.0000 0.0003 -0.0015 -0.0115 -0.0104 0.0303 0.0869 - 0.0202
0.13 0.51 0.16 -0.23 -0.17 -0.28 0.89 0.50
E (-9) (-8) (-7) (-6) (-5) (-5) (-5) (-4)
0.88 0.23 0.36 0.26 0.67 -0.37 -0.44 -0.44
(-11) (-9)
(-8) (-7) (-7) (-6) (-6)
(-5)
101
Volume 153, number 1
2 December I988
CHEMICAL PHYSICS LETTERS
tained by the collocation method with a larger basis set (N= 25). We also carried out a calculation with 25 collocation points placed equally between values of r which were chosen to be in the classically forbidden regions. The results, shown in column E, are seen to be as accurate as those obtained using the Gauss-Hermite points. However, care should be taken in this respect if N is very small, as the method is known to be unreliable in this case for non-Gaussian points (see page 97 of ref. [ 121).
collocation and variational calculations. For simplicity we distribute the functions evenly and choose r ,,,),,= - 3.0 and r,,,,,=6.0, values. The widths of the Gaussians are determined by the parameters Ai, where $,(~)=(2A,/z)‘/~exp[-A,(r-rj)2],
(13)
and these are selected by the Hamilton and Light [ 17 1,
prescription
of
A, =c2/(~~-r,)*,
3.2. Distributed Gaussian functions - convergence with basis set
Aj=4C2/(ri+1
The basis set which has had perhaps the most success in describing the whole spectrum of energy levels is the distributed Gaussian basis proposed by Hamilton and Light [ 171. Here we combine these local basis functions with the collocation method and compare the procedure with the variational approach. We look at a Morse oscillator (B) which has 24 bound states and is defined by ,u= 6, D= 12, (x=0.5 and r,=O.O [17,18]. In order to make the comparison as close as possible we choose the positions of the Gaussians to be the same in both the
The quantity c is left as a parameter which we choose separately for the two methods. In this study the positions of the collocation points are simply taken to be the centres of the Gaussians. In figs. la and lb we present the errors in the energy of the ground, 9th excited and 14th excited states for the variational and collocation methods, respectively. The parameters c were chosen in each case to yield optimum sets of results while avoiding the problems associated with using very broad Gaussians [ 171; their choice is discussed more fully in
.
-r,_l)*,
AN=C2/(rN-rf+-l)2.
(14)
’
. d
1
.
A
L .
I
.
*
. .
.
. .
.
. l
. .
A
.
.
.
1
.
. .
.
.
.
.
.
.
.
.
b
a 25
35
$
25
N Fig. 1. Plot of log,, absolute error in the ground (squares), 9th excited (circles) and 14th excited for (a) the variational method, ~~0.7 and (b) collocation method, c= 0.5.
102
35
N (triangles)
45
i
states, against basis set size
Volume 153, number
1
CHEMICAL
PHYSICS LETTERS
section 3.3. From figs. la and lb it is seen that the collocation and variational methods have almost identical convergence behaviour. It is also interesting to consider the accuracy of the eigenfunctions. As collocation is a fitting technique, the solutions are correct at the chosen points. A measure of the error in the eigenfunctions is thus given by the value of the residue (rjR) between these points. In view of this we calculated the sum of squares of ( r IR) evaluated at the mid-distances between the collocation points. The results indicate that the eigenfunctions corresponding to converged bound states are of good accuracy.
2 December
1988
bidden regions. Again, the basis functions are given the same positions for both collocation and variational approaches, and the collocation points are the centres of the Gaussian functions. In table 2 we give the errors in the level energies for the lowest and highest five states of the Morse oscillator, obtained using the variational and collocation methods. The results demonstrate the similar accuracy of the collocation and variational methods. The difference in detail between the variational results of column A and those of ref. [ 171 we attribute to the use of slightly different positions for the Gaussian functions. It is interesting to note that the results of the collocation method with ~~0.35 are somewhat better than those obtained with c= 0.5. This preference for broader Gaussians may be rationalised by remembering the equivalence of the collocation equations to the Rayleigh-Rilz-Galerkin equations with an Npoint quadrature approximation to the integrals (see eq. (6) ). Increasing the width of the Gaussians may then be considered to improve the collocation results by increasing the accuracy of the integrals. The use of broad Gaussians is essential to the successful application of the collocation technique and we have invariably found the results to improve as the value of c is decreased. This observation also provides an explanation for the adverse comments concerning the collocation method made by Shore [ Ill. In his work, Shore used spline functions which were much nar-
3.3. Distributed Gaussian functions - highly excited states
In the final comparison we examine the ability of the collocation method combined with distributed Gaussian functions to describe highly cxcitcd states. In order to prevent the basis set becoming unnecessarily large we now use (approximately) optimum positions for the Gaussians. As shown previously [ 17,lS 1, a basis set of 67 functions is then sufXcient to give all the bound states for Morse oscillator B. We place 61 of the functions semiclassically at the energy of the 24th bound state using the method of Hamilton and Light [ 17 1. The remaining 6 Gaussians are then placed evenly in the classically for-
Table 2 Compartson of errors in the eigenvalues (d,) for Morse oscillator variational c=O.S, C: collocation c=O.S, D: collocation c=O.35 Exact E”
n
18 19 20 21 22 23
Gaussian
basis set. A: variational
c=O.7,
B:
A,, N=67 A
0 1 2 3 4 5
B using a distributed
B
D
C
- 11.5052 - 10.5469 - 9.6302 -8.7752 -7.9219 -7.1302
4.7 3.4 1.7 2.4 5.2 2.4
(-8) (-8) (-7) (-7) (-7) (-7)
7.8 -4.7 5.4 8.6 -3.5 -1.6
( - 10) (-9) (-9) (-9) (-9) (-8)
-0.6302 -0.4219 -0.2552 -0.1302 -0.0469 -0.0052
7.9 9.4 1.2 2.7 4.2 3.5
(-8) ( -8) (-7) (-7) (-6) (-5)
-3.6 (-8) -2.4 (-8) -2.2 (-8) -1.2(-8) 3.7 (-7) 4.7 (-6)
2.3 -4.3 5.0 -2.0 -4.7 -1.4
(-8) (-8) (-8) (-7) (-7) (-6)
5.4 -5.0 -6.4 3.8 4.8 -5.6
(-10) (-9) (-9) (-8) (-8) (-8)
-1.1 -9.3 -1.0 -1.1 7.1 -1.9
(-5) (-6) (-5) (-4) (-4) (-3)
2.1 -3.6 -4.9 7.2 -7.7 -4.5
(-7) (-7) (-7) (-7) (-7) (-4)
103
Volume 153, number
1
CHEMICAL
PHYSICS LETTERS
rower than the Gaussians of the present studies. In relation to the width of the local functions, we also note that caution should be exercised when very broad Gaussians are used, since problems may be encountered with the basis set. This may be illustrated by reference to the variational results in table 2. The eigenvalues obtained using c=O.7 with the variational method are all strictly upper bounds to the correct values. However, this is not the case for ~~0.5 and for values of c less than this the variational method totally breaks down. The phenomenon is related to the loss of linear independence of the basis set in the numerical calculations and has been discussed previously by Hamilton and Light [ 171. The existence of such problems was the reason for selecting c=O.7 in the calculations presented in fig. la. The problems due to loss of linear independence do not appear to be as severe for the collocation method and generally the algebraic eigenvalue equations may still be solved for very small values of c to yield good results. However, with c < 0.4 spurious eigenvalues have been found to appear very occasionally amongst the physically correct ones. Such eigenvalues are easily discarded as their values are highly sensitive to the basis set used. It should be said that spurious eigenvalues do not occur very often and, in fact, we did not encounter any in the calculations of table 2. Moreover, for values of c greater than 0.4 no spurious eigenvalues have been found in any calculations we have yet carried out. For this reason a value of c= 0.5 was used in the convergence study of section 3.2.
4. Conclusions We have investigated the collocation method for calculating bound states of molecular systems. The technique is simple to implement and totally avoids the need to evaluate integrals over the basis set. The method was tested on the Morse oscillator problem and shown to give accurate results using both simple harmonic oscillator basis functions, which are global, and distributed Gaussian functions, which are local. In particular, the combination of the collocation method with a distributed Gaussian basis was found
104
2 December
1988
to be capable of describing even the highest excited levels of the oscillator. Most importantly, in no case investigated were the variational results significantly more accurate than those of the collocation method, for a corresponding basis. The comparative simplicity of constructing the collocation equations thus makes the method very attractive. This will be particularly true for large systems where multidimensional integrals over the potential must be evaluated if the variational method is to be used.
Acknowledgement It is a pleasure to acknowledge the encouragement of Bill Miller throughout this work and we thank him for many valuable discussions. WY was funded by the National Science Foundation, grant CHE 8416345. ACP gratefully acknowledges the support of a NATG/SERC fellowship.
References [l] J.V. Lill, G.A. Parker and J.C. Light, Chem. Phys. Letters 89 (1982) 483. [2] R.W. Heatherand JCLight, I. Chem. Phys. 79 (1983) 147. 131 J.C. Light, 1-P. Hamilton and J.V. Lill, J. Chem. Phys. 82 (1985) 1400. [4] Z. Bacic and J.C. Light, J. Chem. Phys. 85 (1986) 4594; 86 ( 1987) 3065. [ 51R. Friesner, Chem. Phys. Letters I I6 (1985) 39.
[ 61R. Friesner, J. Chem. Phys. 85 (1986) 1462. [ 71 R. Fricsner, J. Chcm. Phys. 86 ( 1987) 3522. [8] D.O. Harris, G.G. Engerholm and W.D. Gwinn, J. Chem. Phys.43 (1965) 1515. [9] P.F. Endres, J. Chem. Phys. 47 (1967) 798. [lo] A.S. Dickins0nandP.R. Certain, J. Chem. Phys. 49 ( 1968) 4209. [ II] B.W. Shore, J. Chem. Phys. 58 (I 973) 3855. [ 121 B.A. Finlayson, The method ofweighted residuals and variational principles (Academic Press, New York, 1972). [ 131 L. Collatz, The numerical treatment of diffcrcntial cquations (Springer, Berlin, 1960). [ 141 B.W. Shore, J. Chem. Phys. 59 (1973) 6450. [ 151 S.A. OrsLag, Studies Appl. Math. 51 (1972) 253. [ 161 D. Gottlicb and S.A. Orszag. Numerical analysis ofspectral methods: theory and applications (SIAM, Philadelphia, 1977). [ 171 I.P. Hamilton and J.C. Light, I. Chem. Phys. 84 (1986) 306. [ 181M.J. Davisand E.J. Heller, J. Chem. Phys. 71 (1979) 33X3.