Research: Science & Education
Applications of the Variation Method Stephen K. Knudson Department of Chemistry, College of William and Mary, PO Box 8975, Williamsburg, VA 23187-8795 The variation principle plays a significant role in a number of chemical applications. With this stimulus, two examples emphasizing the fundamentals of the method are presented. One example uses the linear variation function version of the method; the other, the nonlinear variation approach. Both use basis sets selected to provide simple integral calculations. Several accounts and applications of the variation principle have appeared in this Journal. Some of these have stressed the chemical significance of the application; other accounts have stressed aspects of the variation method itself. Parson presents a geometrical interpretation of the linear variation method (1). While computational methods are not the focus of the work, examples are presented using Mathematica, the first report of a symbolic math tool in this context. Bendazolli employs a limit form of exponentials for a one-term nonlinear treatment (2). Keeports applies the linear variation method to the particle-in-the-box (PIB), generating wave functions from Fortran calculations (3); later work by the same author considers potentials of the form V(x) = Axn and provides a package of computer code (4). Reed and Murphy point out that the variational principle, being based on an energy criterion, does not necessarily lead to the best possible values for properties, particularly those properties dependent on the amplitude in regions with relatively small influence on the energy (5). Their calculations are all analytic. Fucaloro uses approximate product wave functions for PIB excited states, emphasizing properties in addition to the energy (6). Sims and Ewing apply the linear variation method to a number of model chemical systems; their examples employ Fortran codes from the QCPE (7). Finally, Robiette presents a nonlinear variational treatment of the one-electron hydrogen molecule ion system in the Born–Oppenheimer approximation, determining the energy as a function of internuclear distance using a Fortran subroutine (8). In this note we first illustrate the application of the linear variation method to PIB by use of a basis selected to provide particularly simple integrals, as an aid to pedagogical use of the variation method. The symmetry of the model is exploited to reduce the size of the secular equation. The matrix elements may be computed by hand, by a spreadsheet, by a symbolic math program, or of course by a higher level language. Convergence using the basis is rapid, and the accuracy of the wave function correspondingly high. PIB is considered because it is a popular choice for examples in introductory quantum chemistry texts (9), because there is a readily accessible analytic solution and PIB serves as a model system for physical as well as pedagogical purposes. The second example, for the nonlinear variation method, describes the harmonic oscillator wave function in terms of an exponential function. This example (i) involves simple variational wave functions and integrals; (ii) demonstrates the importance of choosing variational functions that satisfy the constraints of acceptable wave functions; and (iii) speaks to issues in which Gaussian wave functions are used to approximate exponentials.
930
Example 1. The Linear Variational Calculation As commonly defined, the PIB potential is zero in the range 0 ≤ x ≤ L, and infinite outside that range. As a result, the wave function is nonzero only within the box, the range to which we henceforth confine ourselves. The PIB Hamiltonian is then just the kinetic energy operator:
H ={
2
2
h d 2M dx 2
(1)
where M is the mass of the particle. The analytic results for the eigenvalues En and the eigenfunctions ψn(x) are well known: 2
En =
n 2h ; 2 8 ML
ψn(x) = 2 L
1/2
sin nπx L
(2)
where n, which takes on positive integer values, is the quantum number. To emphasize the symmetry of the Hamiltonian, we transform the independent variable from x to t:
t = 2x – 1 L
(3)
In this variable the potential vanishes in the range {1 ≤ t ≤ 1 and is infinite elsewhere, the Hamiltonian operator is given by
H ={
2
2
2h d 2 2 ML dt
(4)
and the analytic solution is
ψn(t) = 2 L
cos nπt , 2 × nπt , sin 2
1/2
n odd
(even parity)
n even
(odd parity)
(5)
The function χ(x) = x(L – x) is frequently used to illustrate the variation method (3). However, the function contains no variational parameters, either linear or nonlinear. To apply the linear variation method, consider the expansion of the variational approximation to the wave function ψ(t) in the set of symmetry adapted basis functions {φj(t)}:
ψn(t) =
Σ c j φ2j+p(t); N
j =0
φk (t) = N k 1 – t 2 t k
(6)
where the {ck} are the variational constants to be determined, N is the number of functions in the basis set expansion, p is a parity index that is either 0 (even parity) or 1 (odd parity), and Nk is the normalization constant, taken to be real. Normalized basis functions are used to facilitate interpretation of the linear coefficients {ck }. The factor of (1 – t2) = (1 – t)*(1 + t) ensures satisfaction of the boundary conditions ψ(t = ±1) = 0. The exact solution van-
Journal of Chemical Education • Vol. 74 No. 8 August 1997
Research: Science & Education ishes linearly at each boundary, since, from eq 2,
Table 1: Overlap Matrix S
sin nπx ~ nπ x, x → 0; L L (7)
sin nπx ~ {1 L
n+1 nπ
L
L – x , x →L
Each basis function φj(t) also vanishes linearly at t = ±1. The basis function φ0 = N0 (1 – t2) is identical (except for normalization) to the single term approximation χ(x) mentioned above. Associated Legendre polynomials of order 2 provide a set of orthogonal polynomials appropriate for an expansion basis, but here the advantages of an orthogonal basis do not compensate for the more complicated integrals it generates. The orthogonal basis would be preferred for computational reasons if large values of N were required. In the notation of this paper, Keeports (3) uses basis functions of the form φk = (1 – t 2k+2 )tp, where p is a parity index as above; this basis leads to slightly more complicated integral matrix elements than the basis presented here. Fucaloro’s single-term approximations have the form φk ∝ ∏k (t – 1 – 2k/n) (4), which again lead to more complicated integrals. The variational solution is given from the secular equation
j\k 1 2 4 6 8
0
2
1
(8)
6
8
0.4411
0.3218
0.2478
0.6547
1
0.9188
0.7940
0.6815
0.4411
0.9188
1
0.9630
0.8913
0.3218
0.7940
0.9630
1
0.9787
0.2478
0.6815
0.8913
0.9787
1
Table 2. Hamiltonian Matrix H
j\k
0
2
4
6
8
0
5.0000
4.5826
3.9698
3.5395
3.2220
2
4.5826
33.0000
42.4476
45.7115
46.5663
4
3.9698
42.4476
72.4286
91.4803
103.3080
6
3.5395
45.7115
91.4803
128.2727
155.9178
8
3.2220
46.5663
103.3080
155.9178
200.2000
2
H – WS c = 0; H jk = φj H φj , S jk = φj |φk
4
0.6547
H jk = –
2h 2 ML
1 {1
2
1–t t
J+K–2
2
K K– 1 – K+ 2 K+ 1 t dt
(12) where W denotes the variational approximation to the energy, c denotes the vector of linear coefficients, and all states are restricted to the same parity manifold (have a given value of the parity index p). The solution requires the determination of the overlap and Hamiltonian matrix elements Sjk and Hjk. All of these matrix elements can be expressed in terms of the simple integral (10) 1
2 , n even n +1 z dz = 0 , n odd n
In =
(9)
{1
The overlap matrix element for the nonorthogonal functions is proportional to the primitive integral 1
s jk =
1–t {1
2
2
t
J+K
dt =
16 J+K+ 1 J+K+ 3 J+K+ 5
(10)
where J = 2j + p and K = 2k + p in eqs 10–12. The diagonal elements of s are the reciprocal squares of the normalization constants Nj, so the nonvanishing overlap integral in the normalized basis is given by
(2J+ 1)(2J+ 3)(2J+ 5)(2K+ 1)(2K+ 3)(2K+ 5) S jk =
(J+ K+ 1)(J+ K+ 3)(J+ K+ 5)
and the Hamiltonian matrix element is found as
1 2
(11)
=
16h
2
ML
2
N jN k
2 JK+ J+ K– 1 J+ K– 1 J+ K+ 1 J+ K+ 3
Matrix elements for the even parity basis set with N = 5 are listed in Tables 1 (S) and 2 (H). Calculation of matrix elements and of individual eigenvalues can be done with a spreadsheet; EISPAK Fortran routines were used for the computations reported. Results for PIB For convenience, define the relative variation energy as the variational energy W divided by the exact energy of the ground state. Convergence of the variational calculation with increasing N, the size of the basis set, is rapid, as can be seen from Figure 1, which plots the relative variation energies for the four even parity states of lowest energy as a function of N. Numerical values are listed in Table 3 for even parity states and in Table 4 for odd parity states. These energies and the wave functions agree with those of Keeports (3), since for a given N the basis functions are equivalent. The one term approximate wave function x(L – x) for the ground state of the PIB, shown by the entry for the N = 1 even parity expansion, provides an energy in error by only 1.3%. The energy of the lowest state of odd parity is in error by 6%. For linear variation calculations, the variational bound property is extended: each variational energy is an upper bound to the energy of successive excited states. That fact and the slower convergence of more highly excited states are displayed in the table. Associated with each energy obtained from the solution to the secular equation is a set of coefficients {cj}. The coefficients for the even parity basis set with N = 8 are listed in Table 5, as are the variational and exact energies. The rapid convergence of the ground state can be seen in the rapid decrease in the values of the coefficients. Rapid convergence persists for several states, but then is lost for states with
Vol. 74 No. 8 August 1997 • Journal of Chemical Education
931
Research: Science & Education 100
.
90
n =0 n =2 n =4 n =6
1 .0
80 70
Wr
.
50 40
.
.
.
ψ
0 .5
60
0 .0
.
30
.
.
.
.
.
-0 .5
-1 .0
20
.
. .
. .
. .
. .
. .
. .
. .
N=1
N=2
N=3
N=4
N=5
N=6
N=7
N=8
10 0
-1 .0
-0 .8
-0 .6
-0 .4
-0 .2
0 .0
0 .2
0 .4
0 .6
0 .8
1 .0
t Figure 1. Relative variational energies, Wr, for even parity states. N designates the number of basis functions in the expansion. The relative energy Wr is the variational energy divided by the energy of the ground state, h2/8 mL2. The exact results for Wr for the four states shown are 1, 9, 25, and 49.
larger quantum number. The Taylor series expansion of the sine function requires more terms as the quantum number increases, and the coefficients show that behavior. But convergence is so rapid that even a three-state basis could be profitably used, leading to a 3 × 3 secular determinant that could be solved without a computer. Figure 2 plots the variational approximations to the four even-parity wave functions of lowest energy, using the eight-state basis listed in Table 3. The exact solutions are indistinguishable from the variational approximations on the scale of the figure. Example 2. Nonlinear Variational Calculation The harmonic oscillator also provides a useful model system since it normally has been considered before the variation principle is introduced. The Hamiltonian is
Figure 2. The four lowest energy, even parity variational wave functions, using a basis set with eight functions. The independent variable t is defined in eq 3.
ψv x = N vH v a x exp {α x 2/ 2
(14)
where the parameters for the oscillator are
ω=
K ; M
α = Mω h
(15)
and Nv is a normalization constant. The function selected to represent the wave function in a variational treatment is an exponential multiplied by a power:
φn x = A nx nexp { γ x
(16)
where x is the extension of the oscillator from its equilibrium position, with the range of x taken to be -∞ ≤ x ≤ ∞; K is the force constant, and M is the mass of the oscillator. The eigenfunctions of this Hamiltonian are products of Hermite polynomials and a “Gaussian”:
where n is a positive integer, An is the normalization constant, and γ is regarded as the nonlinear variational parameter. The variation function is symmetric (n even) or antisymmetric (n odd) in x, to take advantage of the symmetry of the Hamiltonian. The variational function φ vanishes at both large positive and large negative values of x, as is required in order that the wave function can be normalized. As is well known (9), the simple exponential function with n = 0 does not properly describe the ground state; the function with n = 1 does properly describe the lowest energy state of odd parity, v = 1. While the variational method is
Table 3. Relative Variational Energies, Wr , of Even Parity States
Table 4. Relative Variational Energies, Wr , of Odd Parity States
2
H =–
Basis Size, N
E1
2
2 h d + 1 Kx 2M dx 2 2
E3
E5
1
1.0132
2
1.0000
10.3480
3
1.0000
9.0352
35.5594
4
1.0000
9.0003
25.7779
5
1.0000
9.0000
6
1.0000
7 8
932
(13)
E7
Basis Size, N
E2
1
4.2555
E4
E6
E8
2
4.0023
20.3147
3
4.0000
16.2106
57.8067
89.0495
4
4.0000
16.0044
38.1448
25.0307
53.8696
5
4.0000
16.0000
36.1370
73.6350
9.0000
25.0005
49.4508
6
4.0000
16.0000
36.0039
65.1919
1.0000
9.0000
25.0000
49.0204
7
4.0000
16.0000
36.0000
64.0801
1.0000
9.0000
25.0000
49.0004
8
4.0000
16.0000
36.0000
64.0027
Journal of Chemical Education • Vol. 74 No. 8 August 1997
131.5240
Research: Science & Education normally presented as providing an upper bound to the ground state, in fact it provides an upper bound to the lowest energy state of each applicable symmetry. The formulas below are developed for any nonnegative integer value of n. The integrals required to determine the variational energy
φn H φn
E var =
(17)
φn|φn
∞ 0
x k e {ax dx = k! a k+1
2n + 2 2n + 1 2γ
2
An analytic expression for the optimum value of the parameter γ, γm, is found from the variation condition:
dW γ |γ = γm = 0; γ 2m = 1 dγ 2
2n – 1 2n + 1 2n + 2 α (23)
The variational energy, E var, is given analytically by
can be expressed in terms of the familiar integral (10)
Ik a =
V nn = K 2
(18)
E var = 1 2
2n + 1 2n + 2 2n – 1
hω
(24)
The overlap integral is given by 0
φn|φn =
{∞ ∞
=2
x 2ne +2γ x dx +
∞
x 2ne {2γ x dx
0
(19) 2n {2γ x
x e 0
dx = 2I 2n 2γ
γ 2m = 3α ;
We suppress the argument 2γ henceforth in the notation. The Hamiltonian matrix element is found to be
H nn = 2
2 h 1 γ2 + K 2n + 2 2n + 1 I 2n 2 2M 2n – 1 2 2γ
(20)
by use of the recursion formula for the integral,
I k+1 2γ = k + 1 I k 2γ 2γ
(21)
Combining these gives the expression for the variational energy
W=
H nn = T nn + V nn S nn 2
T nn =
In practice, a spreadsheet program determines the value numerically by minimization of the energy expression; spreadsheet programs normally include such optimization tools. For n = 1, both the analytic and the numerical approaches give
h 1 γ2 2M 2n – 1
(22)
φ1 x = x exp {
E var = 3hω
(25)
3α x ; ψ1 x = N 1 x exp {0.5αx 2 (26)
The exact energy is 1.5hω , so the exponential variational function leads to a 15% error in the energy for the odd-parity ground state. This approximate wave function and the exact wave function are shown in Figure 3. The agreement between the variational and exact energies is not particularly good. In confirmation, only very close to the origin do the wave functions closely resemble each other. Curiously, the virial theorem remains satisfied for the exponential variation function: Tvar = Vvar for all n. The simple exponential function with n = 0 fails to provide a description of the ground state. The difficulty is not that the description is poor: the formula for the energy above applied with n = 0 reports the variational energy to be imaginary. The reason for the failure is that for n = 0 the function does not have a continuous derivative at x = 0; instead, it has the famous cusp at the origin, familiar to chemists in association with s-states. Such a function does not satisfy the quantum mechanical postulate that the derivative of the wave function be continuous. A more physical
Table 5. Eigenvalues, E , and Eigenvectors for N = 8 1
2
3
4.93480
44.4132
123.370
1.0328
1.0328
2:
{0.0527
{2.2770
3:
0.0022
1.1645
{14.3456
E cij 1:
4 241.808
{1.03278
1.03047
6.72538
{13.32281 60.36186
5 400.956 {0.96503 19.94911 {144.8569
6 644.3710 0.67557 {18.65494 {178.53
7
8
1294.45
4845.70
{0.40501 12.83853 142.1366
{0.30300 9.93234 {115.2166
4:
-6.215e -05
{0.3311
13.6809
{126.9997
481.3502
{766.66
5:
1.201e -06
0.0615
{7.7770
147.5515
{843.7623
{1690.80
6:
{1.692e -08
-0.0081
2.9187
{102.1427
816.9829
{1992.07
2451.728
2442.418
7:
1.792e -10
{0.7119
40.3953
{417.5012
{1193.95
{1675.093
{1820.452
88.4659
{286.83
452.474
8: {1.109e -12 E (exact) 4.93480
7.731e -04 {4.67e -05 44.41320
0.0886 123.370
{7.10019 241.805
399.719
597.111
707.9575 {1807.763
833.982
608.4604 {1664.857
540.591 1110.33
Vol. 74 No. 8 August 1997 • Journal of Chemical Education
933
Research: Science & Education 0.8
ψ v =1 φ n =1
0 .9
ψ v=0 φ n =0
0.6 0.4
0 .8 0 .7 0 .6
0.2
0 .5
-3
-2
-1
0.0
0
1
2
0 .4
3
0 .3
-0.2
0 .2
-0.4
0 .1
-0.6
-3
-2
-1
0
1
2
3
x/α1/2
-0.8 Figure 3. The exponential variational approximate wave function φn =1 for the lowest energy odd parity state of the harmonic oscillator (solid line). The exact wave function, ψv=1, is shown for comparison (dotted line). The dimensionless independent variable is x/α1/2.
Figure 4. The exponential variational approximate wave function φn =0 for the lowest energy even parity state of the harmonic oscillator (solid line). The exact wavefunction, ψv=0 , is shown for comparison (dotted line).
way to express this is that the average kinetic energy, Tnn in eq 22 , is negative for n = 0. A modified exponential function for the ground state of even parity defined as
Literature Cited
φ0 x = exp {γ x – 1 exp {βγ x β
1. 2. 3. 4.
(27)
(where β is any positive constant) removes the cusp (generates a positive average kinetic energy). The parameter β may also be considered to be an additional variational parameter. Analytic formulas are available, but we have chosen to evaluate β via nonlinear optimization with a spreadsheet, finding
5. 6. 7. 8. 9.
γ = 0.664893α 1/2, β = 19.00943, Evar = 0.511842hω A comparison with the exact result for the energy, E = 0.5ω, shows that the exponential approximation to the Gaussian leads to a 22% error in the ground state energy. The exact and variational approximate wave functions for the ground state are shown in Figure 4, and the comments pertinent for the first excited state also hold for the ground state.
934
10.
Journal of Chemical Education • Vol. 74 No. 8 August 1997
Parson, R J. Chem. Educ. 1993, 70, 115–119. Bendazzoli, G. L. J. Chem. Educ. 1993, 70, 912–913. Keeports, D. J. Chem. Educ. 1989, 66, 314–318. Keeports, D. Am. J. Phys. 1990, 58, 230–234; Keeports, D.; Furth, J. Physics Academic Software; American Institute of Physics: New York, 1993. Reed, L. H.; Murphy, A. R. J. Chem. Educ. 1986, 63, 757– 759. Fucaloro, A. F. J. Chem. Educ. 1986, 63, 579–581. Sims, J. S.; Ewing, G. E. J. Chem. Educ. 1979, 56, 546–550. Robiette, A. G. J. Chem. Educ. 1975, 52, 95–96. For example, Levine, I.; Quantum Chemistry, 4th. ed.; Prentice Hall: New York, 1991; Atkins, P. W.; Molecular Quantum Mechanics, 2nd ed.; Oxford Student Edition, 1983. Dwight, H. B. Tables of Integrals and other Mathematical Data, 4th ed., Macmillan: New York, 1964.
