Numerical Generation of Optimized Discrete Variable Representations

Numerical Generation of Optimized Discrete Variable Representations

J. J. Soares Neto and L. S. Costa

Universidade de Brasília, Departamento de Física

70.910-900 Brasí lia- DF, Brazil

Received January 17, 1997. Revised November 3, 1997.

We develop a procedure for calculating an optimized Discrete Variable Representation (DVR) optimized for a given potential. The method leads to an efficient and compact way to obtain numerical solutions of quantum mechanical problems. The procedure is applied to several physical problems. To illustrate the strength of the algorithm in dealing with multidimensional calculations, we obtain accurate levels up to 19,000 cm^-1 for the vibrational energies of the water molecule.

The efficiency of quantum mechanical numerical calculations rely in great part on the right choice of basis functions. The solution of multidimensional problems requires the manipulation of large matrices and a bad choice of basis functions could easily turn the problem unfeasible. We should also be concerned with the matrix structure and the computer time for the computation of the matrix elements. The calculation of the elements of a matrix makes use of Gaussian quadratures and an excessive number of points may lead to prohibitive time consuming codes. This paper shows a procedure for optimizing all these aspects. The strength of the approach is easily seen when we deal with three-dimensional problems.

There are several ways of expanding the quantum mechanical wave function. The more traditional procedure is the expansion in terms of global basis functions [1-3]. Generally, we choose the functions that diagonalize the kinetic energy operator and integrate numerically the potential using Gaussian quadratures. Another family of procedures are the finite element methods [4-12]. The basis functions are local and we have to discretize the space in many elements to implement the calculation [6]. Finally, we may choose a set of basis functions that have the feature of diagonalizing the potential. The great advantage in this case is that the kinetic energy operator may be integrated analytically, leading to a great efficiency during the evaluation of the matrix elements. The last procedure is called the "Discrete Variable Representation" (DVR) [13-23].

The Discrete Variable Representation has been used extensively by Light and collaborators [13-16] who performed accurate calculations of the ro-vibrational energies of several triatomic molecules. Choi and Light [24] also used this procedure to study Van der Waals molecules. Le Queré and Leforestier [25] used a time dependent approach and a DVR implementation for studying the photodissociation of Ozone. More recently, the method for expanding the multidimensional quantum mechanical wave function [26-32] has appeared in other calculations.

This paper develops a procedure for obtaining a Numerically Generated Discrete Variable Representation (NG-DVR). All publications mentioned before, except the work of Echave and Clary [26], do not use a numerically optimized Discrete Variable Representation for the potential of the system under consideration. In this way, the number of basis functions is not optimized for the particular problem to be solved. The procedure proposed by Echave and Clary [26] does not calculate explicitly the weight of the optimized Gaussian quadrature while the method developed here obtain both points and weights. This feature gives more flexibility for using the NG-DVR. To show its features we apply the NG-DVR for some standard one-dimensional problems. However, the strength of this methodology appears when we treat multidimensional calculations. A three-dimensional code is developed and applied to the calculation of the vibrational modes of H₂O.

The paper is organized as follows. In Section II we develop the theory of the Discrete Variable Representation. In Section III we discuss the numerical procedure for optimizing the DVR for a given potential, and use this methodology to solve some one dimensional problems. Section IV is concerned with the application of the methodology for a multidimensional problem and the calculation of the vibrational energies of the H₂O molecule. Finally, in the Conclusion we present several remarks.

I The discrete variable representation

Most methodologies for dealing with multidimensional problems are based on direct products of one-dimensional sets of basis functions. In view of this, we will treat all basic concepts of DVR using a one-dimensional terminology. Section IV shows an application for a multidimensional problem. We stress that the use of N basis functions in each direction of the space will lead to a matrix with dimension Nd where d is the dimension of the space. It is essential to spend as few as possible basis functions for each direction of the space in order to keep the matrix problem treatable.

The Hamiltonian operator,

is composed of the kinetic energy and potential operators. Let us expand the wave function Y in a set of basis functions f_i(x),i = 1, .., N belonging to a L²-space,

where a_i; i = 1, ..., N are the coefficients of the expansion. We will require that

where x_l; l = 1, .., N are points of a Gaussian quadrature rule. Later, we will develop an algorithm to obtain an optimized Gaussian quadrature.

We call a set of functions obeying Eq. (2) a Discrete Variable Representation (DVR). In order to generate a numerically optimized DVR we start with another set of L² functions g_j(x),j = 1, .., N. The wave function Y(x) may also be expanded in terms of this set,

The finite set of functions g_j(x), j = 1, .., N obey the relation

We will call g_j(x), j = 1, ..., N the primitive functions, and the set f_i(x), i = 1,..,N is called the Numerically Generated Discrete Variable Representation (NG-DVR). The g_j are given functions and the NG-DVR will be obtained, for any potential V(x), by the proposed algorithm.

We use Eq.(4) to obtain

The integral ág_j|f_iñ may be evaluated applying a Gaussian quadrature rule of points and weights x_l, w_l; l = 1, .., N. This quadrature will also be numerically optimized for the potential of the system in question. At this stage, we assume that the optimized quadrature rule has been done and we obtain

Finally, using Eq. 2, we have

We can see that Eq. 2 is automatically verified for the functions f_i(x), i = 1, .., N generated by Eq. 7. It is straightforward to obtain an optimized quadrature rule using Eq. 2 and Eq. 7. For a given quadrature point x_l, we have

From the theory of Gaussian quadrature rules, we know that the points x_l are the eigenvalues of the matrix with elements

where x is the position operator. As it will be shown in the applications in the next two sections, an efficient NG-DVR is obtained when we use as primitive functions, g_j(x), j = 1, .., N, the eigenfunctions of the one-dimensional Hamiltonian . For many real systems, there is no analytical solution for the Schrödinger equation HY = EY, so our procedure is all done numerically.

II Simple examples

Simple examples will be treated in this Section to emphasize the main features of the algorithm. As we mentioned earlier, the strength of the method will become clear for multidimensional problems, and the next Section will deal with a fully three-dimensional calculation. The simplest problem one could treat is the particle in a box. We may solve Eqs. 7- 9 analytically and obtain the quadrature points and weights as well as the DVR's. Muckerman [23] and Colbert and Miller [22] have done that and generated the Discrete Variable Representation associated with the particle in a box. We show in Fig. 1-A one of these functions and in Fig. 1-B five of them. The quadrature points associated with this DVR are equally spaced. From now on we will use this set of DVR to expand the primitive g_j(x), j = 1, .., N defined in Eq. 3.

Figure 1A.
Particle in a box DVR. The system is defined for a box ranging from - 10 to 10 bohr.

Figure 1B.
A set of 5 particle in a box DVR's. The quadrature points are equaly spaced in this case.

In order to have a clear idea of the accuracy of the method we treat now the harmonic oscillator. As the Gaussian quadrature points can be analytically calculated for this problem, we can see how well we are able to reproduce them. The numerical procedure is as follows. First, we use 100 particles in a box DVR's to generate 5 eigenfunctions of the Hamiltonian . These eigenfunctions constitute the set g_j(x), j = 1, .., 5. Eq. 9 is applied to obtain the quadrature points. Table 1 compares the exact 5 Gauss-Hermite quadrature points with the numerically calculated ones. We also provide a 9 point Gauss-Hermite quadrature rule in Table 1. As we notice, the methodology is able to offer highly accurate numbers, and they are accurate up to the tenth figure. Finally, Eq. 7 gives us the NG-DVR optimized for the harmonic oscillator. Fig. 2A shows one NG-DVR associated with the harmonic oscillator and Fig. 2B pictures the full set of five. One may generate as many quadrature points and NG-DVR as necessary to deal with the problem in hand. We also point out that the property expressed in Eq. 2 is accurately satisfied.

Figure 2A.
A Numerically Generated DVR (NG-DVR) for the harmonic oscillator.

Figure 2B.
A set of 5 NG-DVR for the harmonic oscillator.

Finally, we treat the H₂ molecule and also use 100 particle in a box DVR's to calculate 5 eigenfunctions of H₂. These eigenfunctions are our primitive set g_j(x), j = 1, .., 5. Applying the same procedure used for the harmonic oscillator, we obtain the quadrature points optimized for this potential and the NG-DVR. Fig. 3 shows the 5 NG-DVR's as well as the potential for H₂ and the quadrature points.

Figure 3.
Five NG-DVR for the H₂ molecule. The Figure also shows the H₂ potential energy surface and the five quadrature points used to generate the DVR's.

III The H₂O molecule. A three-dimensional example

Now the concepts developed in the previous Sections will be applied to a multidimensional problem. The generalization of the algorithm is straightforward and easy to perform. The ability of the algorithm to generate high accurate quadratures points and NG-DVR will allow us to use the minimum number of basis functions and integration points for solving the multidimensional problem.

The potential energy surface of a triatomic system depends on the three internuclear distances and ro-vibrational calculations have to deal with the expansion of the wave function in these directions. The choice of coordinates and basis functions are crucial for the success of any algorithm. The most commonly used are the hyperspherical coordinates, Jacobi coordinates and mass weighted Jacobi coordinates. We are using a particular set of hyperspherical coordinates that describe the full space by means of a hyper-radius and five hyper-angles. The hyper-radius range goes from zero to infinity while the hyper-angles have finite ranges. Other details about the coordinates we are using are given in previous publications and we refer the reader to them [5,31-35].

Accurate calculations of ro-vibrational levels of triatomic molecules are now feasible due to the speed and memory of the modern computers. Efforts have been done by several groups to obtain the accurate spectrum of tetra-atomic molecules. For this purpose, we should run a full six-dimensional calculation and pay attention in each step of the calculation in order to keep the dimension of the matrix and the number of operations within the limits of the computers. Even for three-atomic systems, the optimization of basis functions is essential in order to calculate high vibrational modes with high accuracy. The wave functions associated with excited states may sample a large region of the three-dimensional configuration space and the basis functions being used should be able to accurately represent such states. Moreover, the number of basis functions must be small in order to keep the matrix dimension manageable.

The optimization procedure we use is as follows. First, we optimize the functions for the hyper-radial direction and, following, hyper-angular basis functions are obtained. We obtain the basis functions for the hyper-radial direction choosing hyper-angles such that the hyper-radial line passes through the potential bottom. The Hamiltonian constrained to the hyper-radial direction defines an one-dimensional problem. The procedure described in Sections II and III are them applied to generate a NG-DVR for the hyper-radial direction. We point out that this set of basis functions brings information about the potential and they are efficient for expanding the total wave function.

The numerically optimized basis functions for the hyper-angular space is obtained applying a two-dimensional finite element method approach. The two-dimensional hyper-angular space is discretized by means of triangles and the Hamiltonian matrix is diagonalized. The eigenvectors are utilized as two- dimensional basis functions for the hyper-angular space. This contraction method has been used before and we refer the reader to the Ref. 5. We keep a number of 110 hyper-angular eigenvectors in the present calculation.

We used the same potential of Miller and Tennyson, Wei and Carrington and Choi and Light. This is an empirical potential energy surface obtained by Jensen [36] who collected a large amount of experimental data to fit it. Table 2 shows states of the water molecule up to 19,000cm^-1. We also compare the calculated numbers with experimental values and calculations by Fernley, Miller and Tennyson (FMT) [27], Wei and Carrington (WC) [30], Choi and Light (CL) [29] and J. J. Soares Neto and F. V. Prudente (NP) [5] in Table 2 and conclude that our calculated numbers agree well with other calculations. We also notice that this calculation improves considerably our previous attempt [5], in special, for highly excited states, where the present basis set describes the wave function more accurately. There are three key parameters that are responsible for the convergence of the calculation and we use the following values for them: 22 hyper-radial basis functions, 3095 triangular finite elements describing the hyper-angular space and the maximum number of 150 hyper-angular contraction basis functions. The dimension of the matrices associated with these parameters is the maximum we can manage with our present computer resources. The final contracted matrix for this calculations is of dimension 1500 ×1500 and the complete run takes 5 hours of a CPU of a Alpha 3000 (Digital) workstation.

IV Concluding remarks

We developed a numerical procedure to optimize a DVR for a given system. The methodology has been applied to the calculation of the ro-vibrational states of the water molecule. We obtained high vibrational levels of this molecule with good accuracy compared to other three-dimensional calculations. The results of this publication are considerably better than those that we have published [5] using another approach for making the expansions along the hyper-radial direction. The application to other three-atomic systems is straightforward and does not require any change in the program since the optimization procedure is automatic. The algorithm deals with both the generation of basis functions and the calculation of quadratures points for the problem.

Öhrn and Linderberg [37] developed a set of hyperspherical coordinates for tetra-atomic systems. Presently, we are making efforts to develop an algorithm to calculate bound states of tetra-atomic molecules using these coordinates and the Numerically Generated DVR described in this paper. We are also developing a contraction scheme similar to the one used for the water molecule.

Acknowledgments

This work was supported by "Conselho Nacional de Desenvolvimento Científico e Tecnológico" (CNPq) through grants to J. J. Soares Neto and L. S. Costa. The computer facilities were provided by the "Laboratório de Cálculo Cientí fico" (LCCFIS) at Department of Physics, University of Brasília, Brazil.

References

1. P. L. Goodfriend, Chem. Phys. Lett. 186, 509 (1991).

2. T. D. Crawford and W. Yang, Chem. Phys. Lett. 192, 45 (1992).

3. A. Preiskorn, D. Frye and E. Clementi, J. Chem. Phys. 94, 7204 (1991).

4. J.J. Soares Neto and J. Linderberg, Comp. Phys. Comm. 66, 55 (1991).

5. J. J. Soares Neto and F. V. Prudente, Theor. Chim. Acta 89, 415 (1994).

6. P. G. Ciarlet, The Finite Element Method for Elliptic Problem, North- Holland, Amsterdam, 1978.

7. D. H. Norrie and G. de Vries, The Finite Element Method, Academic Press, NewYork, 1973.

8. K. W. Morton, Comp. Phys. Rep. 6, 1 (1987).

9. R. Jaquet, Theor. Chim. Acta 71, 425 (1987).

10. R. Jaquet, Chem. Phys. 118, 17 (1987).

11. N. Sato and S. Iwata, J. Comp. Chem. 9, 222 (1988).

12. L. Yang, D. Heinemann and D. Kolb, Chem. Phys. Lett. 192, 499 (1992).

13. Z. Bacic and J. C. Light, J. Chem. Phys., J. Chem. Phys. 85, 4594 (1986).

14. Z. Bacic, R. M. Whitnell, D. Brown and J. C. Light, Comp. Phys. Comm. 51, 35 (1988).

15. Z. Bacic and J. C. Light, J. Chem. Phys. 86, 3065 (1987).

16. J. C. Light and Z. Bacic, J. Chem. Phys. 87, 4008 (1987).

17. M. Mladenovic and Z. Bacic, J. Chem. Phys.93, 3039 (1990).

18. Z. Bacic and J. Z. Zhang, Chem. Phys. Lett. 184, 513 (1991).

19. Z. Bacic and J. Z. Zhang, J. Chem. Phys. 96, 3707 (1992).

20. V. Szalay, J. Chem. Phys. 99, 1978 (1993).

21. C. Leforestier, 4th Topsfe Summer School on Time Dependent Methods in Quantum Mechanics, Copenhagen, 1991.

22. D. T. Colbert anq W. H. Miller, J. Chem. Phys. 96, 1982 (1992).

23. J. T. Muckerman, Chem. Phys. Lett. 173, 200 (1990).

24. E. Choi Seung and J. C. Light, J. Chem. Phys. 92, 2129 (1990).

25. F. Le Quéré and C. Leforestier, J. Chem. Phys. 94, 1118 (1991).

26. D. Wang and J. M. Bowman, J. Chem. Phys. 100, 1021 (1994).

27. J. A. Fernley, S. Miller and J. Tennyson, J. Mol. Spectr. 150, 597 (1991).

28. J. Tennyson, J. Chem. Phys. 98, 9658 (1993).

29. S. E. Choi and J. C. Light, J. Chem. Phys. 97, 7031 (1992).

30. H. Wei and T. Carrington, Jr., J. Chem. Phys. 97, 3029 (1992).

31. J. J. Soares Neto and J. Linderberg, J. Comp. Chem. 12, 1237 (1991).

32. J. J. Soares, S. B. Padkjaer and J. Linderberg, Int. J. Quantum Chem. S23, 127 (1989).

33. J. Linderberg and B. Vessal, Int. J. Quantum Chem. 31, 65 (1987).

34. M. Mishra, J. Linderberg and Y. Ohrn, Chem. Phys. Lett. 111, 439 (1984).

35. J. J. Soares Neto and J. Linderberg, J. Chem. Phys. 95, 8022 (1991).

36. P. Jensen, J. Mol. Spectrosc. 133, 438 (1989).

37. Y. Ohrn and J. Linderberg, Mol. Phys. 49, 53 (1983).

Publication Dates

• Publication in this collection
  28 Apr 1999
• Date of issue
  Mar 1998

History

• Reviewed
  03 Nov 1997
• Received
  17 Jan 1997