Approximation Methods: The Variational Theorem
Copyright (c) RDMCHEM LLC 2019
Overview
Interactive Concepts
Applications in Chemistry
A tenant of the Born interpretation of quantum mechanics is that contained within the wavefunction, ψ, is all the information one can possibly know about a quantum system, where for time-invariant systems, ψ is the solution to the Schrodinger equation. In this chapter, we consider only time-invariant Hamiltonians and look at solutions to the time-independent Schrodinger equation (TISE), H∧ψ = Eψ . For all but a handful of systems, which include the particle in a box, the harmonic oscillator, and the hydrogen atom, exact analytical solutions to the TISE are not known, and one must resort to approximate solutions. In this chapter, we explore an extremely important theorem in quantum mechanics known as the Variational Theorem, which provides a means for finding an upper bound to the true energy of a system.
Consider that we can rewrite the TISE as
E = ∫ψ∗H∧ ψ dx∫ψ∗ ψ dx
where for a normalized wavefunction, the denominator is equal to 1. Since we do not in general know ψ, let us choose a trial wavefunction, ϕ, as a guess to the solution, where ϕ must obey the same boundary conditions as ψ. The variational theorem says:
Variational Theorem: For a trial wavefunction, ϕ, that obeys the boundary conditions of the TISE, then the expectation value of the energy for the trial wavefunction, hereafter referred to as the variational energy, Ef , given by
Ef = ∫φ∗H∧ φ dx∫φ∗ φ dx
is greater than or equal to the true energy, E. That is, Ef ≥ E, where the equality holds only for φ = ψ.
That is, the energy calculated using the above equation is always greater than the true energy. If we choose better and better guesses for the trial wavefunction φ, then Ef will continue to decrease until it converges (within some acceptable tolerance) to the true energy, E!
In the Interactive Concepts section below, we focus on the particle in a one-dimensional box and consider increasingly more flexible trial wavefunctions to show the power of the variational principle. In the Applications to Chemistry section, we introduce the Morse oscillator, a more realistic potential for describing vibrational energies of diatomics and apply the variational theorem to approximate the allowed energies. We also see how electronic structure methods based on the variational method can be used to calculate the energies of atoms and molecules.
To show the variational theorem in action, we want to consider a quantum system for which we know the correct solutions with which to compare. For this we choose the beloved 'particle in a box'! (For more on this topic, see the Particle in a Box exercise). The particle in a box entails a particle with mass m moving in a one-dimensional potential V(x)
V(x) = 0−L2<x<L2∞else,
where L is the length of the box. The solution to the TISE yields normalized wavefunctions and energies given by
ψ(x) = 2LsinnπxL+nπ2
and
En = n2h28 mL2,
where h = Planck's constant and n = 1, 2, 3, .... For simplicity, we will work in atomic units, in which me = Z = e = 4πe0 = 1. The unit for length is the Bohr radius, a0, the unit for mass is the mass of the electron, me, and the unit for energy is the Hartree, Eh. Therefore, for an electron in a box of length L = 2 a0, we have
En = n2π28
Now, what should we use as an approximate wavefunction, ϕ? In what follows we consider as an initial guess the function ϕ(x) =1 - x2 and then explore different methods for improving upon this initial guess.
Trial Wavefunction #1: φx=1−x2
Consider the trial wavefunction φ(x) =1 - x2. Does 4 obey the boundary conditions of the problem?
Answer
Yes!! By inspection, the trial wavefunction goes to zero at x = ± L/2 = ±1. It is also continuous with a continuous first derivative, making it an acceptable trial wavefunction for this problem!
Using the Maple input below, calculate the energy, Ef, of this trial wavefunction and compare it to the true ground state energy. You may also use the plot function to compare ϕ to the true wavefunction, ψ.
restart : withplots:
φx≔1−x2:
Etrial≔intφx⋅diffdiff−12⋅ φx,x,x,x=−1..1intφx⋅φx,x=−1..1
Etrial≔54
Etrue≔π28
evalfEtrial−Etrue
0.016299450
PercentError≔evalfEtrial−EtrueEtrue⋅100
PercentError≔1.321183653
A≔plotsinPi⋅x2+Pi2,x=−1..1,color=Red: B≔plotphix,x=−1..1,color=Blue: displayA,B
In accord with the variational theorem, the energy of the trial wavefunction is higher than the true ground state energy! From the plot of the true wavefunction (red) and the trial wavefunction (blue), we can see that there is room for improvement! Can you think of another trial wavefunction that could be tried? Enter your trial wavefunction for ϕ(x) in the Maple input above to see if your trial wavefunction is better than 1-x2. Remember, your trial wavefunction must obey the original boundary conditions!
Trial Wavefunction #2: φx=1−x2+ βx2−x4
One way to add flexibility to a trial wavefunction is to add a 'variational parameter'. Consider the trial wavefunction given by
φ(x ; β) =1 - x2+ β (x2 - x4)
where (x ; β) implies that ϕ is a function of x and parameterized by the variational parameter β. That is, we want to find the value of β that gives us the best trial wavefunction of the form given above. Given the infinite possibilities for β, how are we to proceed? A conceptual way is to use Maple's interactive plot feature to plot both ψ(x) (the true wavefunction, red) and ϕ(x ; β) (blue), vary β until the two plots match, and then use the variational method to calculate the energy!
φx,β≔1−x2+β⋅x2−x4:
ψx ≔ sinπ⋅x2+π2:
#interactivepsix,phix,beta:
Exploreplotsmultipleplot,cosπx2,x=−1..1,color=Red,1−x2+β−x4+x2,x=−1..1,color=Blue,labels=x,,parameters=β=−10.0..10.0,initialvalues=β=0.
β
Enter the value of β that allowed for the best match between the trial and true wavefunctions below to see the new trial energy:
beta≔−0.256 :
Etrial≔intφx,β⋅diffdiff−12⋅ φx,β,x,x,x=−1..1intφx,β⋅φx,β,x=−1..1;
Etrial≔1.234165113
Etrue≔π28;
evalfEtrial−Etrue;
0.000464563
PercentError≔Etrial−EtrueEtrue⋅100;
PercentError≔0.03765605842
Were you able to improve upon the original trial energy?
Yes! With the optimized beta from the plot above, we find the percent error is only 0.04%, as opposed to the 1.32% from Trial Wavefunction#1.
While varying β and comparing differences in the two wavefunctions is a simple, intuitive approach, it is not generally useful since we do not in practice know the true wavefunction! A more mathematically rigorous approach is to find the minimum of trial with respect to β by solving
dEφdβ=0
Use the Maple input below to solve for b that minimizes dEφ/db and calculate the resulting variational energy to compare with the value that you were able to find from the plot above.
restart :
Etrial≔intφx,β⋅diffdiff−12⋅ φx,β,x,x,x=−1..1intφx,β⋅φx,β,x=−1..1:
dE≔diffEtrial,β:
g≔solvedE=0,β:
evalfg
−0.220749972,−7.317711566
Here we see two solutions, one very close to what was seen from the plot approach above! Plug in each value to the expression for the variational energy to see which is the 'global' minimum!
β≔−0.22075:
Etrial≔1.233718703
Etrue≔Pi28;
0.000018153
Error≔Etrial−EtrueEtrue⋅100;
Error≔0.001471426757
Which value of b yielded the minimum energy?
We find that β = -0.22075 is the absolute best value for β, and that for the functional form φ(x ; β) =1 - x2+ β (x2 - x4) for the trial wavefunction, 0.00147% error is the best that we can do!
Trial Wavefunction #3: φx=α1−x2+ βx2−x4
In order to improve further upon our trial wavefunction, we can increase the number of variational parameters. This method will be used in the Applications in Chemistry section below to estimate the ground and excited state energies for the Morse Oscillator. The equations generated for this multivariational approach are provided in the subsection below. The Maple input relies on the expressions introduced therein.
Linear Combinations of Basis Functions for the Trial Wavefunction, 4
This subsection explores the equations encountered when the trial wavefunction is a linear combination of basis functions. It is the background information for Interactive Concepts: Trial Wavefunction #3 and Chemical Applications: The Morse Oscillator.
Consider as a trial wavefunction a linear combination of N basis functions g1, g2, ..., gN, each of which satisfying the boundary conditions of the system:
φ(x) = c1 g1(x) + c2 g2 (x) + ⋅⋅⋅ + cN gN (x).
The set of expansion coefficients c1, c2, ..., cN, become variational parameters! It can be shown that the above expression for φ(x) leads to a generalized matrix equation:
Hc = Ef Sc
where Ef is the variational energy, c is the column vector with elements c1, c2, ..., cN, H and S are N × N matrices with matrix elements given by
Hij = ∫gi H∧gj ⅆx
Sij = ∫gi gj ⅆx.
The S matrix is known as the overlap matrix and reduces to the identity matrix for an orthonormal basis set. The matrix equation above leads to a set of linear equations that have a non-trivial solution if and only if
det [ H - Ef S] = 0
For a linear combination of N basis functions, one gets an N-order polynomial expression for the energy Eφ, and it can be shown that the ordered set of Eφ is an approximation to the ordered set of Etrue ground and excited states!
Consider the trial wavefunction
φ(x) = c1( 1 - x2) + c2 ( x2 - x4),
We know we must solve the set of equations given by
det H11 − Eφ S11 H12 − Eφ S12H21 − Eφ S21H22 − Eφ S22= 0
The solution to the above equation gives rise to a quadratic equation for Ef . The smaller of the two solutions is an upper-bound to the ground state, and the larger of the two solutions is an upper bound to the first even excited state! Why is it not an approximation to the very first excited state?
Since (1 - x2) and (x2 - x4) are both even functions, and any linear combination of two even functions is still even. The first excited state is an odd function, which cannot be approximated by an even function!
Use the Maple input below to calculate the ground and first even excited states of the particle in a box and compare them to the true values.
restart:withLinearAlgebra:
g1≔1−x2:
g2≔ x2− x4:
H11≔intg1⋅−12⋅diffdiffg1,x,x,x=−1..1:
H12≔intg1⋅−12⋅diffdiffg2,x,x,x=−1..1:
H21≔intg2⋅−12⋅diffdiffg1,x,x,x=−1..1:
H22≔intg2⋅−12⋅diffdiffg2,x,x,x=−1..1:
S11≔intg1⋅g1,x=−1..1:
S22≔intg2⋅g2,x=−1..1:
S12≔intg1⋅g2,x=−1..1:
M≔MatrixH11−Etrial⋅S11,H12−Etrial⋅S12,H21−Etrial⋅S12,H22−Etrial⋅S22:
solveDeterminantM=0.,Etrial
12.76628130,1.233718703
Etrue≔Pi28.
Etrue≔1.233700550
Error≔1.233718703−EtrueEtrue⋅100
Compare this answer to the optimal value of Trial Wavefunction #2. Is it better? How good is the 12.766 value to the first even (n = 3) excited state?
Etrue≔32⋅Pi28.:
Error≔12.76628130 − EtrueEtrue⋅100
Error≔14.97730953
We find that this is a very poor approximation to the n = 3 state (or the first even excited state)! What could be done to improve this result??
Use more basis functions in the linear combination!
Morse Oscillator and the Vibrational Energies of HCl
In the section above, we used the variational theorem to approximate the energies of the particle in a box, and we saw that the use of variational parameters can be used to get us quite close to the true energy of the system. In this section, we consider a slightly less trivial system: the Morse oscillator, which can be used to model a vibrating diatomic molecule. We consider the use of trial wavefunctions that are linear combinations of underlying basis functions, as seen in Interactive Concepts: Trial Wavefunction #3. While the math may be more involved, the underlying ideas of the variational theorem hold: as long as a trial wavefunction ϕ satisfies the boundary conditions of the Hamiltonian, then the expectation value of the energy for ϕ will be higher than the true energy, and as more basis functions and variational parameters are used in the approximation, the better the energy will be.
The Morse potential is an analytical function for the potential energy that includes anharmonicity for diatomic vibrations and is given by
V(r) = De { 1 - exp(-β (r - r0))}2
where De is the well depth, r - r0 is the displacement from the equilibrium bond length r0, and β is a molecular parameter related to the force constant of the molecule. Table 1 contains a list of these parameters for HCl, which are determined by fitting the potential to spectroscopic data. Also contained in Table 1 are the first two actual energies of the Morse oscillator and the corresponding fundamental frequency in wavenumbers. This can be compared with the observed fundamental frequency.
Table 1: Morse Parameters for HCl
r0
2.409 a0
0.9890 a0-1
De
0.1697 Eh
μ
1785.7 me
E0
0.006749 Eh
E1
0.01984 Eh
υ~⋅Morse
2827.2916 cm-1
υ~⋅obs
2885.9 cm-1
Use the Maple input below and the information from Table 1 to generate the Morse potential for HCl.
withplots: plot0.1697⋅1 − exp−0.9890⋅r − 2.409 2,r=1..10,y=0..0.5
The Hamiltonian for the Morse oscillator in atomic units is given by
H∧= −12 μd2dx2+De 1 − exp−β x2
where x is the displacement from equilibrium distance, and m is the reduced mass of HCl (1H35). While an analytic series solution does exist for the Morse oscillator, we will use the variational method to calculate approximate solutions. We will choose as a trial wavefunction a linear combination of N basis functions:
φ(x) = c1g1(x) + c2g2(x) + ... + cNgN(x)
As shown in Interactive Concepts: Trial Wavefunction #3, a linear combination of N basis functions leads to N approximate solutions. What should we use as basis functions? Given that we are interested in the ground and low lying vibrational states for which the potential is still rather harmonic, let's try harmonic oscillator wavefunctions as our basis!
Use the Maple input below to calculate the lowest N states of the Morse oscillator using a basis of N harmonic oscillator functions. (Recall that the vibrational quantum numbers for the first N states of the harmonic oscillator are v = 0, 1 , 2, 3, ..., N - 1.) Begin with N = 2 to get an approximation to the ground and first excited state of the Morse oscillator and compare them with the analytical values provided in Table 1. Once the energies have been calculated, use them to calculate the fundamental IR vibrational frequency of the HCl molecule in wavenumbers (cm-1) and compare this value with the analytical value from Table 1, 2827.2916 cm-1.
restart:withLinearAlgebra:withorthopoly:withplots:
N≔10:
r0≔2.409: β≔0.9890: De≔0.1697: μ≔1785.7 :k≔2⋅De⋅β2:α≔sqrtk⋅μ:ω≔sqrtkμ:
Nvv≔2v⋅v!−12⋅απ14:
gx,v≔Nvv⋅exp−α⋅x22⋅Hv,sqrtα⋅x:
Vx≔De⋅1−exp−β⋅x2:
for i from 0 to N−1 by 1 do gpx,i≔−12⋅μ⋅diffdiffgx,i,x,x+Vx⋅gx,i end do:
M≔MatrixN,N,shape=symmetric :
for i from 0 to N−1 by 1 do for j from i to N−1 by 1 do Mi+1,j+1≔intgx,i⋅gpx,j,x=−infinity..infinity end do end do
E≔EigenvaluesMatrixM,shape=symmetric:for i from 1 to N by 1 do Ei end do
0.00674891921590597
0.0198360575752205
0.0323804258588314
0.0444821633770340
0.0566871381495674
0.0704423866035275
0.0864225594292201
0.103813903970456
0.134191601102211
0.213316627270026
PercentError0≔abs0.006749−E10.006749⋅100
PercentError0≔0.00119697872338527
PercentError1≔abs0.01984−E20.01984⋅100
PercentError1≔0.0198710926400233
Fundamental≔E2−E1⋅219476.329478
Fundamental≔2872.31709047309
PercentError≔abs2827.2916− Fundamental2827.2916⋅100
PercentError≔1.59253083317739
With only 2 basis functions, we actually get a pretty good approximation to both the ground and 1st excited state, with 0.5% and 6.9%, respectively! Why do you think the agreement is so good?
For low lying states, the Morse potential is fairly harmonic, so the harmonic oscillator basis functions are very good approximations to the ground and first excited states, though we do see that as we increase in energy, the agreement diminishes since more anharmonicity is added to the Morse potential. Note the approximation to the fundamental frequency is approximately 12%, reflecting the larger error in the first excited state.
Now what happens as you increase N? Do you get better approximations to the ground and first excited state, consistent with the Variational Theorem? To answer this, N = 5 in the Maple input and recalculate the approximate energies.
As N increases, we not only get more approximations to excited states, but the approximations get better and better for lower lying states, in accord with the variational theorem. The new errors are only 0.0040% and 0.22% for the ground and first excited state, respectively. We also note that we get a better approximation to the fundamental frequency, with now only 2% error!
Now let N = 15 in the Maple input above and recalculate the energies. Note that this might take some time! Once the energies are calculated, use the Maple input below to explore some of the key quantitative and qualitative differences between the harmonic oscillator and the Morse oscillator.
In the plot below, we show the first 5 states of the harmonic oscillator (red) and the Morse oscillator (blue).
EHO≔Vector5: for i from 0 to 4 by 1 do EHOi+1≔ω⋅i+0.5 end do:
X≔Vector5,1:Y≔Vector5,2:B≔plotX,EHO, x=0..3,y=0..0.07, tickmarks=0,0,symbolsize=15, style=point,symbol=circle,color=Red:C≔plotY,E, symbolsize=15, style=point,symbol=circle,color=Blue: plotsdisplayB,C
X≔VectorN,r0:A≔plotDe,x=0..8 ,y=0..De, tickmarks=0,0, color=Black:B≔plotX,E, symbolsize=15, style=point,symbol=circle,color=Blue:C≔plotDe⋅1−exp−β⋅x−r02,x=0..8,y=0..De: plotsdisplayA,B,C
1) Is the spacing between adjacent energy levels evenly spaced, as seen in the harmonic oscillator? What happens to this spacing as energy increases?
We see that the energy levels are not evenly space as seen in the harmonic oscillator. We note that at the bottom of the well, where the approximate energies are more accurate, the spacing decreases for increasing energies. If we were to use an infinite basis set, we would notice this trend to continue for all states in the well. The anharmonicity in the potential decreases the energy.
2) Are there an infinite number of 'bound' vibrational energy levels, as seen in the harmonic oscillator? If not, how many bounds states does HCl have?
While there is an infinite number of allowed energy levels in the Morse oscillator, only 12 of these are 'bound'. One of the key qualitative features of the Morse potential is that it describes bond breaking. For vibrational states with an energy below De (below the black line in the plot), the wavefunction must decay to zero at large values of x, which implies that the vibration is 'bound'. That is, the atoms do not have sufficient kinetic energy to overcome the potential attraction, and the bond remains intact. However, for energies above De, where the kinetic energy is greater than the potential energy at large x, the wavefunction may continue to oscillate at large x (similar to the free particle), and the bond breaks. In contrast, the harmonic oscillator bond cannot break since its potential energy diverges for increasing x. This gives rise to an infinite number of bound states for the harmonic oscillator.
Variational Calculations of HCl Potential Energy at Equilibrium Bond Distance
In the previous section, we used an analytical fit to the PES for HCl and used the variational principle to calculate the vibrational manifold. In this section, we actually calculate the potential energy for HCl at its equilibrium bond distance, 2.409 a0. We will use the Hartree-Fock method, which approximates the true electronic wavefunction as being a single Slater determinant constructed from molecular (spin) orbitals. The MOs are optimized by assuming each electron interacts only through an average electric potential due to all other electrons in the system. Each MO can be represented in a given basis set, similar to how the vibrational energy levels of HCl were represented in a basis of harmonic oscillator functions. As the number of basis functions increases, we should see the electronic energy decrease variationally.
We will explore two series of related basis sets, each of increasing size: (sto-3g, 6-31g, and 6-311g) and (cc-pvdz, cc-pvtz, cc-pvqz, and cc-pv5z).
withQuantumChemistry:Digits≔15:
HCl≔H,0,0,0,Cl,0,0, 2.409:
AObasis≔sto-3g:
E≔EnergyHCl,method=HartreeFock,basis=AObasis;
E≔−454.90257928
AObasis≔6-31g:
E≔−459.87664400
AObasis≔6-311g:
E≔−459.90417048
AObasis≔cc-pvdz:
E≔−459.90375490
AObasis≔cc-pvtz:
E≔−459.91888469
AObasis≔cc-pvqz:
E≔−459.92307201
AObasis≔cc-pv5z:
E≔−459.92379791
For each series, does the energy decrease variationally as the size of the basis increases? Which series of basis sets gave rise to better energies? What did notice regarding the time required to carry out each calculation as the size of the basis set increased?
Yes, for the STO-3G, 6-31G, and 6-311G, we notice the energy decreases. The energy is converged to approximately 0.03 Hartrees, or 17 kcal/mol.
For the cc-pVDZ, cc-pVTZ, cc-pVQZ, and cc-pV5Z series, we notice that the cc-pVDZ is consistent with the 6-311G calculation but continues to decrease for increasing basis sizes. The cc-pV5Z calculation is converged to approximately 0.005 Hartrees, or only 3 kcal/mol and is 12 kcal/mol smaller than the best 6-311G calculation.
The price one pays for accuracy, however, is time. The calculations for the larger basis sets were considerably longer. While the time was not intractable for the HF method, for other more accurate electronic structure methods, this increase can be intractable!
Download Help Document
