J . Phys. Chem. 1988, 92, 3091-3096 number. All three of these predictions are in agreement with what has been seen experimentally. The more general conclusions of this analysis include the following: (1) highly non-Franck-Condon vibrational propensities due to the appearance of d/dR and R2factors in the radial and angular nonadiabatic coupling matrix elements (see 1l a and 11b); (2) strong J dependence for one component of each A-doublet and weak J dependence for the other; (3) propensity for small changes in the rotational qu_antum number, with the minimum charge determined by the 1 value of the ejected electron; (4) weak J dependence for the vibration-induced rates and strong J dependence for rotation-induced rates.
3091
Acknowledgment. We acknowledge the financial support fo the National Science Foundation (Grants No. CHE-8206845 and CHE-8511307) and the donors of the Petroleum Research Fund, administered by the American Chemical Society. We also acknowledge the Harris Corp. for their generous computer system grant and the National Science Foundation for its San Diego Supercomputer time award. We also thank the San Diego Supercomputer Center for their computer time grant. G.C. thanks the Polish Academy of Sciences within the program CPBPOl.12 for their support. R.A.K. thanks the University of Utah for their partial support through the Graduate Research Fellowship program.
An Automatic Grid Generation Scheme for Pseudospectral Self-Consistent Field Calculations on Polyatomic Molecules Richard A. Friesner Department of Chemistry, The University of Texas at Austin, Austin, Texas 78712 (Received: August 17. 1987; In Final Form: November 13, 1987)
We describe a stable and efficient approach to generating a reliable polyatomic grid for self-consistent electronic structure calculations using the pseudospectralmethod. Systematicconvergence of spectroscopic constants to the Roothaan HartreeFock values as the grid size is increased is demonstrated. The approach offers an estimated 1-2 order of magnitude reduction over alternative numerical self-consistent field algorithms in the number of grid points per atom required to produce accurate results.
I. Introduction In a series of previous papers,'-3 a new algorithm for solution of self-consistent electronic structure equations has been proposed and implemented. The algorithm involves a synthesis of conventional quantum chemical techniques with the pseudospectral m e t h ~ d an , ~ efficient numerical method originally applied to large-scale hydrodynamic simulations5 and recently employed in various contexts in chemical physics problems.6 The essential feature of the approach is use of both a basis set and a numerical grid representation of the solution, allowing accurate evaluation of derivatives and integrals but also permitting multiplications to be carried out efficiently on the physical space grid. The net results for the Hartree-Fock equations are elimination of twoelectron integrals and an scaling of both integral evaluation and Fock matrix assembly, thus potentially yielding 2-3 orders of magnitude reduction in computation time for large molecules. In ref 3, quantitatively accurate results (as compared to equivalent calculations done via the Roothaan procedure') for the Friesner, R. Chem. Phys. Lett. 1985, 116, 39. Friesner, R. J. Chem. Phys. 1986, 85, 1462. Friesner, R. J . Chem. Phys. 1987, 86, 3522. Orszag, S. A. Stud. Appl. Math. 1972, 51, 253. Orszag, S. A. Stud. Appl. Math. 1971, 50, 293. Orszag, S.A. Phys. Reu. Lett. 1971, 26, 1100. Gottlieb, D.; Orszag, S.Numerical Analysis of Spectral Methods; Theory and Application; SIAM: Philadelphia, 1977. (5) Brachet, M. E.; Meiron, D. E.; Orszag, S. A,; Nickel, B. G.; Morf, R. H.; Frisch, V. J . Fluid Mech. 1983,130,411. Patera, A. T. J. Comput. Phys. 1984, 54, 468. Marcus, P. S.J. Fluid Mech. 1984, 446, 45. Orszag, S.; Patera, A. T. Phys. Rev. Lett. 1980,45, 989. Yakhot, V.; Orszag, S.A,; Pelz, R. B. In Proceedings of the 9th International Conference on Numerical Merhods in Fluid Dynamics; Springer-Verlag: New York, 1985. Tuckerman, L. S.; Marcus, P. S.Ibid. Springer-Verlag: New York, 1985. (6) Kosloff, D.; Kosloff, R. J. Comput. Phys. 1983, 52, 35. Gerber, R. B.; Yinnon, A. T.; Kosloff, R. Chem. Phys. Lett. 1984, 105, 523. Light, J. C.; Hamilton, I. P.; Lill, J. V. J. Chem. Phys. 1985, 82, 1400. Selloni, A,; Carnevali, P.; Car, R.; Parrinello, M. Phys. Reu. Lett., in press. (7) Roothaan, C. C. J. Rev. Mod. Phys. 1951, 23, 69. (1) (2) (3) (4)
0022-3654/88/2092-3091$01 S O / O
water molecule were obtained for the dissociation energy, force constants, and equilibrium geometry. This demonstrated that the pseudospectral method was capable of achieving the high accuracy reequired for quantum chemical calculations on polyatomic molecules. However, a cumbersome nonlinear optimization procedure was required to achieve these results. Furthermore, this procedure is not easily generalized to large polyatomics. Consequently, a more automated and universal implementation is required if the pseudospectral method is to become a useful quantum chemical approach which is easily applied to a userspecified molecular system. In this paper, the most difficult step in that direction is taken. The key problem identified in ref 3 is that of automatic grid generation for a polyatomic system. A practical solution to that problem is described below. There are still details of the scheme which can be improved; we expect that such improvement will allow reduction of the number of grid points required. However, the method as depicted should be directly applicable to an arbitrary polyatomic system. The explicit demonstration of this (along with any improvements developed along the way) will appear in subsequent publications. Here, we report results for the water molecule which are analogous in quality to those in ref 3. In contrast to the work reported in that paper, the grid design implemented here is a systematic procedure which should not require reoptimization for each molecular species. No new timing results are presented here; the principal concern is description of the grid generation algorithm and establishment of accuracy for that algorithm. However, the conclusions regarding the number of grid points per atom, which is a critical parameter in determining the efficiency of any numerical selfconsistent field (SCF) method, are in accord with those in ref 3. The requisite number of points is 1 or more orders of magnitude less than that used in previously described numerical schemes, whether for the Hartree-Fock (HF) equations or in the less demanding context of local density theory. The reason such large 0 1988 American Chemical Society
3092 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 reductions are possible will be discussed in detail in the conclusion. 11. Pseudospectral Method for the Hartree-Fock Equations We first present a brief review of the basic principles of our pseudospectral S C F approach. The description given below is identical with that found in section I11 of ref 3; it is repeated here for the convenience of the reader, who may be unfamiliar with the mathematical structure of the method. Further details can be found in ref 1-3. The basic idea of the pseudospectral method is to use both a physical space (Le., grid) representation +(r) and a spectral (basis set) representation c of the occupied molecular oribtals and transform between them. Producing function values on a grid given the basis set and coefficients is straightforward. The inverse transform is accomplished via a least-squares procedure in which additional functions (dealiasing functions) are included and then discarded; this step is the only approximation (other than basis set incompleteness) made in the method. The number of grid points must be greater than or equal to the number of basis functions. The use of dealiasing functions has proved to be crucial in obtaining accurate results with the pseudospectral method. The pseudospectral Fock operator F, (see below) acts on a coefficient vector c and produces a function on the physical space grid. This function can be well-represented by the basis set, but there are some components produced (the alias) which are outside of the basis set. If left alone, these components are misinterpreted by the basis set, leading to fluctuations as a function of geometry which degrade the accuracy of the relative energy. By including additional functions, the great majority of the alias is filtered from F,c, rendering the pseudospectral calculation nearly equivalent to a conventional one (in which the alias is exactly removed by analytic integration). More discussion of this can be found in ref 2. The basis set coefficients satisfy the usual Fock equation Fc = S c E
(1)
where S is the spectral overlap matrix, and F is given by
F = Ho
+ 25-
K
(2)
The operators H,,(the standard one-electron Hamiltonian) and S are constructed by conventional analytic evaluation of integrals over the basis set. The two-electron operators J and K are formed via the equations J = Q*Jps
(3)
K = Q*Kps
(4)
where J, and K, produce on the physical space grid the result of acting with the two-electron Coulomb and exchange operators on a coefficient vector c. The operator F,, is defined as F, = 25, -
Kps.
The operators J, and Kp are constructed from the three-center, one-electron integrals
where xn and x m are the usual atomic basis functions. Formulas for J,, and K, are given in ref 1 and 2. The inverse operator Q is obtained from the normal leastsquares equations, which give the best fit of the basis and dealiasing functions to the result of the operation Fp&, i.e. Q = S[R+wR]-'R+w
(6)
Here, R is the matrix of basis and dealiasing functions on the grid, R+ the transpose, w a diagonal matrix of least-squares weights, and S the usual analytic overlap matrix. The dealiasing functions are initially orthogonalized to the basis set by using analytic overlap integrals and hence are simply ignored in subsequent calculation. Note that the S operator is needed in eq 6 to project the fitted expansion onto the individual basis functions. This operator was
Friesner inadvertently omitted in ref 3. (It was employed in the code used . . to generate the results of that paper.) An important modification of the above scheme is to construct all one-center, two-electron integrals analytically and use these to assemble the purely atomic parts of the Fock operator. This guarantees separation to analytic HF atoms and removes the largest errors near the atomic nuclei. Once the Fock operator is constructed, the usual iteration of the pseudoeigenvalue equations is used to obtain a solution. Canonical orthogonalization is employed to deal with the overlap matrix. Because the pseudospectral Fock operator is not necessarily Hermitian, a projection operator technique is used to enforce orbital orthogonality. The latter methodology is unimportant; test calculations show that a straightfoward symmetrization of the Fock matrix (averaging the off-diagonal elements) also yields acceptable results. 111. Automatic Grid Generation Procedure
A. General Principles. In this section, a systematic procedure is described for constructing a reliable and efficient pseudospectral grid which should be applicable to an arbitrary polyatomic molecule. Refinements of the procedure will allow reduction in the number of grid points; in general, we have utilized a sufficient number of points to secure convergence without complicated optimization procedures. The results nevertheless provide an idea of the number of points which will be required in a general polyatomic calculation. The procedure has two distinct parts. First, an atomic grid is generated for each atom. In future work, optimized atomic grids associated with each element of the periodic table and corresponding atomic basis set will be developed for molecular applications. This approach is analogous to the design of atomic basis sets in conventional H F methods. The principles described below appear to be quite general and efficient and present few conceptual or numerical difficulties. The second stage consists of merging the atomic grids to form a molecular grid. This is a much more complex task, particularly because one must worry about accuracy on several different levels. The reader should recall that our goal is to reproduce the Roothaan (variational) results; this means that we would like our numerical integration scheme to be as close as possible to analytic integration. On the other hand, it is important to avoid introducing spurious variations in the total energy as a function of geometry when calculating force constants. Balancing these two imperatives (in general, accurate integration schemes require moving the grid as one moves the nuclei) is a subtle design problem. The solution presented below appears to be both general and reasonably accurate; however, tests in which the same algorithm is applied to a series of molecules will be required before the method can be used with confidence. B. Atomic Grid. As in ref 3, an atomic grid is constructed as a series of radial shells. The optimization of the radial shell positions is an essentially quantitative problem which can easily be addressed by carrying out and optimizing atomic H F calculations, designing the radial grid to obtain the best agreement with the Roothaan results for the basis set of interest. In the calculations which follow, we utilize a hyperbolic tangent spacing of radial points8 in each of three segments, the latter corresponding to (roughly) the core, bonding, and extended regions of the atom. On the other hand, we have discovered an angular point distribution which is a qualitative improvement over that used previously. The idea is actually rather simple; the quadrature scheme for each radial shell is designed to perform exact integration on spherical harmonics up to and including order L. The appropriate pseudospectral grid for spherical harmonics has been worked out by G l a t ~ m a e i r . ~There are ( L + 1)* spherical harmonics to be handled; this means that ( L + 1)2grid points on the (8) Thompson, J. F.; Warsi, Z.; Mastin, C. Numerical Grid Generation; North-Holland: New York, 1985. (9) Glatzmaier, G . A. J . Comput. Phys. 1984, 55, 461; personal communication.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3093
Pseudospectral Calculations on Polyatomic Molecules sphere are required. The optimal distribution is to have L + 1 rings of constant 0, each containing L 1 equally spaced points in the I#J direction. The phases of the rings (Le., the initial location of the first I$ point) are staggered so as to produce maximum coverage of the sphere. The location of the 0 rings is given by the zeros of the Legendre polynomial of order L 1. The quadrature weight of each point is obtained by solving the collocation equations on the sphere and determining .the coefficient of the L = 0 spherical harmonic (the only one which survives integration on the sphere) for a function which is unity at the grid point in question and zero elsewhere. Formally, if the points on the sphere are labeled ( P j ) y z ! )then 2,
+
+
(7)
+
+
Here Y is a ( L 1) X ( L 1) matrix of the ( L + 1)2 normalized spherical harmonics evaluated on the ( L + 1)' grid points positioned as described above, with Yoo occupying the first column of Y. [Y-lIlj is the magnitude of the Yoo component contained in a function which is unity at Pj and zero elsewhere on the sphere. The constant K is then given by
K =
lTl'" Yoo(O,q) sin 8 d%d p =
(2i~)l/~
(8)
An atomic grid is completely specified by the number of radial shells, the radial position of each shell, the value of L assigned to each shell, and the overall phase of each shell (i.e., its orientation in three-dimensional space). The appropriate set of L values can again be determined empirically by experimenting with a series of molecular calculations. In general, one would expect that small L values would be sufficient near the nucleus, while the largest values would be required in the bonding region. The shell orientations should be designed to reflect molecular symmetry. In this paper, we make only a minimal attempt to carry out optimization of the atomic grid. Our main purpose is to show that the method is stable and accurate and can be made to converge reliably as the grid size is increased. To this end, a few calculations are shown in which the values of L and number of radial shells are systematically increased. C. Merging of Atomic Grids. Two fundamental strategies exist for merging atomic grids. Most simply, one can just combine the grids as they are without further adjustment, moving each atomic grid as its position is varied. This approach was adopted in ref 3. It has the advantage of smooth dissociation to the separated atom limit. On the other hand, it is difficult to compute the appropriate renormalization of quadrature weights required to achieve accurate numerical integration with such a merged grid. This can lead not only to deviations of the total energy from the variational results but also to instabilities in the least-squares procedure. In ref 3, these instabilities were minimized by nonlinear optimization; however, this procedure is not easy to automate efficiently. An alternative approach is to truncate each atomic grid, e.g., by discarding points which are closer to another atom than to the "parent" atom. This is equivalent to demanding that each atomic grid perform integration inside its own Voronoi polyhedron. As will be described below, a simple and reasonably accurate algorithm to take care of boundary effects can easily be produced. The principal difficulty with this method is that the grid is constantly changing as a function of geometry, both in the number of points and in their weights. Probably, it is possible to design a procedure in which this has no adverse consequences; however, this is a difficult task in practice. The effect on relative energies is rather small, affecting only the weakest force constants. Nevertheless, some way of dealing with this is essential if accurate force fields are to be produced with the pseudospectral approach. The general solution adopted here is to make the grid as consistent as possible in the region of the potential minimum for the purpose of calculating force constants. (There is only a negligible effect on the equilibrium geometry.) The perturbations introduced by the normal grid alteration procedure are important only in this region, as their absolute size (on the order of au) is ordinarily swamped by the correct change in energy as a function of geom-
etry. It is only near the minimum that these errors are at all significant (and then, only at the level of the second derivatives). Furthermore, it is only angular variation that presents difficulties; radial perturbation of the grid is easily accomplished in an unbiased fashion. The reasons for this will be discussed in detail below. A decisive point favoring the above scheme is that it leads naturally to efficient implementation of analytic gradient techniques. That is, we show that it is appropriate to neglect variation of the grid weights with geometry when computing derivatives. In essence, this approach maximizes cancellation of error in the computation of relative energies near the minimum. Such an approach greatly reduces the amount of work that will be necessary to develop gradient algorithms, because one does not then have to consider differentiation of the least-squares fitting matrix with respect to the nuclear coordinates. The point is that the errors induced by this neglect are locally quite small-indeed, smaller than would be produced by most algorithms designed to generate the appropriate corrections. Implementation of the distance criterion for neglecting points is straightforward. For atom A and a grid point at a given angular orientation (%,I$), the critical radius R , beyond which one enters the Voronoi region of a neighboring atom B is R, =
RAB
2 cos E
(9)
where RAB is the distance between atoms A and B and is the angle between a unit vector connecting A and B and a unit vector along the direction (e,$). In order to produce a smooth interpolation scheme in which points can be created and destroyed as a function of geometry, we retain a grid point if it is one radial shell past the critical radius computed from its angular direction. The weight of this point is determined by its contribution to the integration inside the Voronai boundary, which is calculated in a trapezoid rule type of integration scheme via Lagrange interpolation.8 That is, the radial weight of a point just past R , is given by
r Rc R - R,,
Here Rj-l is the radial value just below R,. The total weight of a point is a product of the radial weight multiplied by the appropriate angular weight obtained from the spherical harmonic analysis described above (Le., eq 7). As eq 10 is appropriate for one-dimensional integration, we must also multiply by R,' to obtain a proper three-dimensional volume element, i.e.
The integral in eq 10 is easily evaluated to yield an analytic formula for wR. Similarly, the weight of a point at a value immediately below R , must be corrected; a similar Lagrange interpolation scheme yields
Thus, a point Rj enters the grid with zero weight (when R , = R j - ] )and the weight smoothly grows to the full value it would attain in the untruncated atom. There is a discontinuity in the derivative at the initial entrance into the grid; this is a very small numerical effect, having consequences only for force constant computation. A more serious objection is the rapid variation in R, for small values of cos E. Because there is a maximal value of R for each atom, the singularity in eq 9 at cos E = 0 is not relevant. However, the change in the grid as a function of angular rotation causes unwanted biasing of the angular force constants. (Note that, in contrast, the variation of R, with RAB is rather slow and unproblematic.) As mentioned above, the solution is to eliminate the dependence of R , on cos E in the region of the minimum. In the
3094
Friesner
The Journal of Physical Chemistry, Vol. 92, No. 11, I988
TABLE I: Dealiasing Functions for All Pseudospectral Calculations"
atom
L
zmin
zmsx
0 0 0 0 H H H
S
0.1 0. I 0.1 0.1 0.05 0.1 0.1
40.0 13.5 2.4 2.7 13.5 27.0 3.0
P d
f S
P d
"Exponents have ratios of roughly a factor of 3, beginning with Zmi, and ending with Z,,,. Exponents already in the basis set (e.g., for oxygen, there is a d function which has Z = 0.8) are not repeated. present code, one chooses an initial geometry and simply leaves all the 4 values fixed when varying the geometry around the minimum. The procedure is so stable that the results for the force constants are virtually independent of the starting geometry, within the range of values utilized. When the switch is made to analytic gradient methods, this methodology will be automatically built into force constant computation, so that the grid generation algorithm can proceed in the normal fashion for all geometry points. The generality of the procedure, of course, can only be demonstrated when computations are carried out for many molecules. Nevertheless, the present results are quite encouraging, particularly when contrasted to other solutions to this problem which we have investigated. To implement this scheme in the present calculations, we must choose an initial starting geometry to fix the cos 4 values. We select the geometry R = 1.78, 0 = 106.0' (very close to the equilibrium geometry) for this purpose. Numerical experiments (not shown below) confirm that the choice of starting point is unimportant (as it should be, if the arguments made above justifying the approach are correct). The entire grid merging procedure described above is easily automated in a few hundred lines of computer code for an arbitrary polyatomic. The computation time to generate the grid and appropriate least-squares weights is negligible.
IV. Results A . Description of the Calculations. As in ref 3, we study the water molecule in symmetric geometries. However, in the present code, symmetry adapted functions are not employed; instead, symmetry is built in through the grid design. For the oxygen grid, for example, each radial shell is constructed with C,, symmetry. A consequence of this is that only odd values of L can be employed for oxygen. The basis set is again the 6-31G** basis set of Pople and co-workers.1° We should mention that slightly erroneous exponents were used in ref 3; this discrepancy has now been corrected, and all one-electron matrix elements computed in our program agree with those in GAUSSIAN 8 2 to many significant figures. In ref 3, the set of dealiasing functions had to be carefully optimized in order to avoid instabilities. The substantial increase in accuracy of the new integration grid means that this is no longer a crucial issue. To demonstrate this, we employ here complete (with regard to angular components) sets of dealiasing functions whose only characteristic features are that the values of the exponents are spaced by a factor of 3. The minimum and maximum values of the exponents are chosen to cover the range of interest determined by the basis set; the precise values of these choices are not crucial to obtaining accurate results. Angular sets are included up to one value of L beyond the largest value in the atomic basis set, e.g., f functions are incorporated on oxygen and d functions on hydrogen. Obviously, exponents which already appear in the basis set are not duplicated. Table I lists the minimum and maximum exponents for each atom and orbital shell (s, p, d, f). Future work will undoubtedly involve optimization of these functions via atomic calculations; however, we thought it important to show in this paper that such optimization is no longer necessary. (10) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J . Chem. Phys. 1972, 56. 2257. Harihan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213.
TABLE 11: Comoarison of Pseudosoectral Results with GAUSSIAN 82
parameter
pseudospectral
GAUSSIAN 82
Emin
-7 6.0235 8 1.784 106.3 0.632 0.030 0.178
-76.02362 1.783 106.0 0.635 0.030 0.178
Re 0, ~ R R
fR8 fat
Parameters for the pseudospectral grid are those of Table V. All quantities are reported in atomic units, with angles measured in radians. Force constants and equilibrium geometries were obtained by fitting the energy evaluated at nine water geometries to a quadratic polynomial in the bending and symmetric stretch coordinates. To ensure that there were no artifacts in the fitting procedure, results from GAUSSIAN 82 were also fit at these same geometries. The results reported for GAUSSIAN 82 in Table I1 were evaluated via this procedure. Agreement with published values for the equilibrium geometry" is quantitative, while the force constants agree quite well with near H F limit calculations.'2 B. General Properties of the Results. In the following section, we report specific results for selected grid configurations utilizing the algorithm described above. Here, we want to make some general comments on the behavior of the method, so as to give the reader a feeling for which aspects of the results are most sensitive to details of the implementation, where potential problems in future applications lie, and possible alternate directions which it would be fruitful to pursue. First, the total energy, orbital eigenvalues, and wave function coefficients are invariably extremely close to the Roothaan H F results. In hundreds of test calculations with a variety of, e.g., different approaches to treating the problem of interfacing the atomic grids at their boundaries, the total energy was always within 0.001 au of the Roothaan result at any geometry, and typically somewhat closer. Orbital eigenvalues displayed a similar agreement, accurately reflecting variation as a function of geometry as well. The normalized basis set coefficients agreed to within 0.001 (again, typically more accurately). These results imply that the pseudospectral wave functions are in every sense an isomorphic replacement for conventional Roothaan results and that all of the existing folklore concerning the applicability of H F models at various basis set levels can be taken over. With regard to the energy, bond energies and barrier heights can clearly be evaluated to yield results entirely equivalent to the conventional ones, as the discrepancy between the two methods at any geometry is smaller than the uncertainty induced by the HF approximation even when it is most reliable (1-2 kcal/mol). The equilibrium geometry is also computed with high reliability without great sensitivity to the methodology of grid interfacing. The equilibrium bond length is almost always within 0.003 au of the Roothaan value, while the bond angle agrees to fl'. An increase in the number of grid points (see below) procedures results that are even closer to the conventional ones. Similar stability and accuracy are observed for the radial force constantf,, which is typically within 2% or 3% of the Roothaan figure. (Note that this translates into even better agreement for vibrational frequencies, which depend on the square root of the force constant.) The force constants which depend upon the angular variable 8,oH VRsand fss) are more sensitive to the grid design. Convergence of these quantities to the appropriate values requires in general more grid points (in both the angular and radial degrees of freedom) and use of a suitable grid merging procedure (i.e., that described above). Nevertheless, careful attention to these details does permit systematic convergence of these force constants to the Roothaan values. Other grid merging procedures typically led to uncontrollable errors of 0.01-0.02 au (angular units expressed in radians). ( 1 1 ) Szabo, A,; Ostlund, N . Modern Quantum Chemistry; MacMillan: New York, 1982. (12) Ermler, W. C.; Kern, C. W. J . Chem. Phys. 1971, 55, 4851.
Pseudospectral Calculations on Polyatomic Molecules
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3095
axis is used to define 0.) 4 is the angle made with the z axis; Le., 4 = 0 lies along the positive z axis. The angular definitions apply to both the local atomic and global coordinate system, with the positions z = 0 and y = 0 moved to the appropriate locations for the former. Transforming from the local to the global coordinate system is straightfoward. To define the angular point distribution, we specify three parameters designated R L , , RL2, and RL,. From the first radial shell up to RL,, we set L = 2 (Le., nine points per shell). From RL, to RL2 we employ L = 4 ( 2 5 points per shell) while from RL, to RL, L = 6 (49 points per shell) is utilized. For shells greater than R L , we revert back to L = 4. This is a straightforward scheme which builds in the essential intuition about how L should vary as a function of R and permits one to systematically increase the angular resolution of the grid by, e.g., increasing RL3 or decreasing RL2. The orientation of the oxygen shells is determined by the symmetry constraints. Only two orientations are possible; one grid point on each odd-numbered ring is either pointing directly between the two hydrogen atoms or is 180' from this orientation (as explained above, the even rings are 180' out of phase). Because the bonding region is by far the most important one (as well as the most difficult to converge numerically), we choose the former of the two above-mentioned orientations. (There are more odd than even rings, and the odd ring in the y z plane is the most important one.) The orientation of the hydrogen grids is constrained only by reflection symmetry through the y z plane at x = 0, and by the requirement that each hydrogen be the mirror image of the other across the xz plane a t y = 0. Again, we locate points in the bonding region by orienting the odd rings so that one grid point is along the 0-H bond axis. Note that the entire issue of spherical shell orientation does have to be dealt with systematically for the general polyatomic case; the success of the straightforward approach taken here suggests that this should easily be accomplished. One further point is that the same grid is used for both hydrogen and oxygen. This is clearly not the optimal arrangement; the way to improve matters is to make the cell filled by the oxygen grid larger, Le., to give the oxygen a larger finite element volume. One could then design a specially optimized grid for hydrogen. In a version of the program intended for practical chemical applications, this sort of technical optimization will be implemented extensively. Here we avoid considering the problem, as timing results are not our main concern. However, in estimating the scaling of the method with increasing basis set size, this point should be taken into account. Table I1 compares converged results obtained from a large grid (after merging, the total number of grid points is roughly 2000, or 650 points per atom) for spectroscopic constants, equilibrium position, and the total energy with those obtained from GAUSSIAN 82 (using the same 6-31G** basis set) and with near HF limit calculations. Agreement is obviously quantitative. The grid generation parameters used to produce this grid are presented in Table V; the actual values of the radial shell positions and L values for each shell are displayed in Table VI. Table I11 displays results obtained by fixing the radial grid and varying R L , systematically. Convergence to the correct values is apparent. Table IV presents similar results obtained by holding the R L parameters at values sufficient to secure convergence and varying the radial grid parameter NR2. The values of the Rj and the 6, are held constant throughout at values (those specified in Table V) found to produce a reasonable rate of convergence. (These are the parameters which will require detailed optimization via atomic calculations in future work.) The above results demonstrate that the pseudospectral S C F method can be used in a rational, unbiased fashion to produce
TABLE 111: Pseudospectral Results as a Function of Radial Resolution' parameter NR, = 8 NR, = 10 NR, = 12 Emin
-76.02433 1.782 106.4 0.629 0.034 0.185
Re 9, fRR
fR0
he
-76.02379 1.783 106.3 0.639 0.030
0.182
-76.02358 1.784 106.3 0.632 0.030 0.178
"All grid parameters other than NR2 are as in Table V.
TABLE I V Pseudospectral Results as a Function of Resolution" parameter RL? = 15 RL3 = 20 E,," -76.02382 -76.02354 Re 1.784 1.784
4
104.5 0.645 0.033 0.213
fRR fR0
he
106.4 0.63 1 0.030 0.176
Angular RL3 = 25 -76.02358 1.784 106.3 0.632 0.030 0.178
"All grid parameters other than RL, are as in Table V. In what follows, the calculated spectroscopic constants of the water molecules are explicitly displayed as a function of one radial and one angular grid generation parameter, selected because the results exhibit the greatest dependence on them. We have made no serious attempt to optimize the grid; we simply wish to demonstrate convergence as a function of increasing grid size. Thus, the number of grid points used here can almost certainly be reduced substantially. Furthermore, preliminary results indicate that extremely small grids (as little as 75 points per atom), while inadequate to compute force constants, can generate reasonably accurate total energies. This means that such grids can be used while locating the potential minimum and securing most of the iterative convergence, whereupon a switch to a more accurate grid is made. A multigrid technique designed along these lines could provide an extremely efficient method of determining molecular structures at the a b initio level with a large basis set.
V. Numerical Results All computations were performed on the University of Texas Cray X-MP supercomputer. Energies were converged to lo-' au. We do not explicitly report eigenvectors or eigenvalues below; as stated previously, these quantities agree quite closely with the Roothaan results for all cases discussed here. The following parameters are used to characterize the grid, which is then generated from these parameters as described in section 111. There are three contiguous radial regions for each atom, the boundaries of-which are delineated by four radial values R , , R2,R,, R,. To construct the hyperbolic tangent grid in each region, we specify two values of a spacing parameter 6, one for each extremum of the region. To simplify the number of parameters, the low end of region J is constrained to have the same spacing parameter as the high end of region j - 1. There are thus also four spacing parameters, labelled in a manner analogous to the R values. Additionally, each radial region is assigned a specified number N R j of radial shells. Before proceeding to a detailed discussion of the angular grid, it will be useful to specify the atom-centered and global coordinate systems more precisely. We follow the usual convention for the water molecule in which the molecule lies in the y z plane at x = 0, with the oxygen atom at (O,O,O)and the water molecules at equal positive z displacements and equal in magnitude but opposite in sign y displacements. The angle 0 is defined as the angle made by a coordinate vector with the x axis. (Note that this differs from the usual spherical coordinate convention, in which the z TABLE V: Grid Parameters for Converged Calculation" R1 R2 R3 R4 81 62 0.01
0.5
3.5
8.0
0.01
0.2
63
84
NR,
NR2
NR3
RL,
RL2
RLj
0.4
2.0
10
12
6
7
10
25
"All parameters are defined in the text; results for this grid are shown in Table 11.
3096
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
TABLE VI: Radial Shell Positions and L Values for Atomic Grid Specified in Table V shell no. radius L shell no. radius L 1
2 3 4 5 6 7 8 9 10 II I2 13
0.010 0.017 0 028 0.047 0.078 0.126 0.197 0.289 0.395 0.500 0.698 0.907 1.129
-7 7
2 2 2
-
1
2 4 4 4 6 6 6
14 15 16 17 18 19 20 21 22 23 25 24 26
1.363 1.612 1.877 2.159 2.460 2.782 3.128 3.5 3.952 4.533 5.307 6.386 8.0
6 6 6 6 6 6 6 6 6 6 6 6 6
quantitatively accurate results which are equivalent to those generated by conventional ab initio H F programs. The remaining issue is the magnitude of the timing advantage; we hope to have results for this (beyond the conclusions of ref 3) in a few months. VI. Conclusion The results presented here, even without optimization of the radial grid, represent a remarkable improvement in the efficiency of numerical quadrature in molecular H F or density functional (DF) calculations. For example, the DVM of Ellis and cow o r k e r ~requires ~~ over an order of magnitude more grid points per atom to achieve converged results for molecular force constants using a D F approach. Similarly, least-squares numerical H F calculations of TholeI4 were unable even to achieve reasonable convergence for eigenvalues (this is accomplished effortlessly here) for a diatomic using thousands of grid point per atom restricted to two dimensions. As a final example, the extremely complicated numerical DF method of Feibleman15 utilized roughly 30000 points to study a diatomic (again in two dimensions). Furthermore, it is important to realize that the number of grid points per atom is a crucial determinant in the efficiency of a computational S C F method. To assert that one has achieved N3 scaling is meaningless if the “scaling” prefactor is so large that the value of N required to realize an advantage in an actual calculation is untreatable in any case. If, for example, DF advocates (who have adopted numerical grid methods perforce, not by choice) could really achieve accurate results at the H F level in orders of magnitude less computation time than their competitors (the implication of claims regarding N3 scaling), the DF approach would have long since replaced H F technology. In fact, the large prefactor has up to now made all grid-based methods inferior for computations which demand high numerical precision. (It is a different story if one is interested in an electronvolt as opposed to wavenumber level of accuracy; then, D F approaches have been able to get away with coarse grids, although still not as sparse as the ones used here.) The approach described above, in contrast, is specifically making claims regarding overall computation times. Of course, these claims will have to be rigorously demonstrated for a series of (1 3) Case, D.A. Annu. Reo. Phys. Chem. 1983,33, 1 5 1. Baerends, E.J.; Ellis, D.E.; Ros, P.Chem. Phys. 1973, 2, 41. (14) Thole, B. T. Inr. J . Quantum. Chem. 1985, 28, 535. (15) Feibelman, P.J. J Chem. Phys. 1984, 81, 5864.
Friesner molecules in head-to-head competition with GAUSSIAN 82. Our preliminary indications are that a substantial advantage will be obtained even for small molecules once every section of the code is brought up to a reasonable level of efficiency. We believe that the power and flexibility of the modified pseudospectral algorithm and the variety of available augmentational strategies (e.g., the multigrid approach mentioned above) allow one to surpass both the older grid-based methods and the Roothaan algorithm. In effect, the approach advocated here returns electronic structure theory to the mainstream of modern large-scale numerical computation, where the type of efficiency issues raised above and the type of solution strategies adopted are part of the normal methodology. It is worth reviewing once more the combination of techniques which yield such high accuracy with such a sparse three-dimensional mesh. Firstly, the use of analytic integration whenever it is inexpensive is crucial, and avoids, e.g., the necessity of using pseudopotentials in the core regions. (Pseudopotentials can, of course, still be employed, but it is probably not a good idea to do so for, eg., first-row atoms.) The fact that all purely atomic terms are evaluated analytically is also appealing from a conceptual point of view; the pseudospectral calculations determine only the twoelectron, two-or-more-center terms. Secondly, the benefits of using least squares and dealiasing instead of straight numerical quadrature are substantial. One knows perfectly well that high angular quantum numbers and high-frequency radial terms are unimportant, and failure to exploit this information is very expensive. (Ordinary quadrature is not “smart” enough to utilize the information fully.) Furthermore, our numerical methods do extensively use overlap integrals (which are also relatively inexpensive), thus introducing some analytically based elements into the quadrature scheme. Thirdly, the carefully designed grid generation scheme introduced here is specially tailored for molecular problems. The electronic structure problem is so hard, the level of accuracy required so demanding, that generic numerical approaches are inevitably doomed to failure. This is why numerous trial-and-error computations with a variety of schemes had to be made (the number of unpublished calculations performed by this author between ref 3 and this paper is quite substantial), and why the method is in far from what will be its final form. The next step will follow the basic approach taken by Pople and co-workers, extensive and systematic optimization for each element of the periodic table. To make the sort of investment of human and computation time required to carry out these optimizations, one must be convinced that the underlying algorithm is sufficiently superior to its already established competitors to be worth the effort. Unlike the grid generation procedure described in ref 3 (which was presented with a disclaimer regarding its general utility), we believe that the procedures described above satisfy this criterion. Acknowledgment. We thank Laurette Tuckerman for many useful discussions and for pointing out to us the work of Glatzmaier. We also thank Youngdo Won for carrying out computations with the GAUSSIAN 82 system of programs. This work was supported by a grant (DMB8416842) from the National Science Foundation. R.A.F. is an Alfred P. Sloan Foundation fellow and a Camille and Henry Dreyfus Teacher-Scholar. Use of the Cray X-MP supercomputer at the University of Texas Center for High Performance Computing is gratefully acknowledged.
