3.2. Linear variational calculations

In this section, we describe two quantum mechanical problems that can be analyzed numerically with a linear variational calculation. In both cases, a generalised matrix eigenvalue problem (3.13) must be solved, which can easily be done using a program like MATLAB.

3.2.1 The infinitely deep potential well

The potential well with infinite barriers is given by $$V (x) = \left\{ \begin{array}{ll} \infty \phantom{xx}& {\rm for} \phantom{xx} |x|>|a|, \\ 0 & {\rm for} \phantom{xx} |x|\leq |a|. \end{array} \right.$$ It forces the wave function to vanish at the boundaries of the well ($x=\pm a$). The exact solution for this problem is known and treated in every textbook on quantum mechanics (see for example Griffiths). Here we discuss a linear variational approach to be compared with the exact solution. We take $a=1$ and use natural units such that $\hslash^2/2m = 1$.

As basis functions we take simple polynomials that vanish on the boundaries of the well: $$\psi_n (x) = x^n (x-1)(x+1),\; n=0, 1, 2, \ldots$$ The reason for choosing this particular form of basis functions is that the relevant matrix elements can easily be calculated analytically. We start with the matrix elements of the overlap matrix, defined by $$S_{mn} = \left\langle\psi_n | \psi_m \right\rangle= \int_{-1}^{1} \psi_n(x) \psi_m(x) dx.$$ There is no complex conjugate with the $\psi_n$ in the integral because we use real basis functions. Working out the integral gives $$S_{mn} = \frac{2}{n+m+5}-\frac{4}{n+m+3}+\frac{2}{n+m+1}$$ for $n+m$ even; otherwise $S_{mn} = 0$.

We can also calculate the Hamilton matrix elements -- you can check that they are given by: $$\begin{aligned} H_{mn} & = \left\langle\psi_n |p^2 |\psi_m \right\rangle= \int_{-1}^{1} \psi_n(x) \left(-\frac{d^2}{dx^2} \right) \psi_m(x) dx \\ & = -8 \left[ \frac{1-m-n-2mn}{(m+n+3)(m+n+1)(m+n-1)} \right] \end{aligned}$$ for $m+n$ even, else $H_{mn} = 0$.

The exact solutions are given by $$\tag{3.19}\label{Eq:EigDeep} \psi_n(x) = \left\{ \begin{array}{l} \cos(k_n x) \phantom{xxxx} n {\text{ odd}}\\ \sin(k_n x) \phantom{xxxx} n {\text{ even and positive}} \end{array} \right.$$ with $k_n = n \pi/2$, $n=1,2, \ldots$, with corresponding energies $$E_n = k_n^2 = \frac{n^2 \pi^2}{4}.$$ For each eigenvector $\bf C$, the function $\sum_{p=1}^N C_p \chi_p(x)$ should approximate an eigenfunction \eqref{Eq:EigDeep}. The variational levels are shown in table below, together with the analytical results.

$$\begin{array}{rrrrr} \hline \hline N=5 & N=8 & N=12 & N=16 & \text{Exact} \\ \hline 2.4674 & 2.4674 & 2.4674 & 2.4674 & 2.4674 \\ 9.8754 & 9.8696 & 9.8696 & 9.8696 & 9.8696 \\ 22.2934 & 22.2074 & 22.2066 & 22.2066 & 22.2066 \\ 50.1246 & 39.4892 & 39.4784 & 39.4784 & 39.4784 \\ 87.7392 & 63.6045 & 61.6862 & 61.6850 & 61.6850 \\ \hline \hline \end{array}\\ $$

The table above lists energy levels of the infinitely deep potential well. The first four columns show the variational energy levels for various numbers of basis states N. The last column shows the exact values. The exact levels are approached from above as in Figure (3.1).

3.2.2 Variational calculation for the hydrogen atom

As we shall see in further on in this course, one of the main problems of electronic structure calculations is the treatment of electron--electron interactions. Here we develop a scheme for solving the Schrödinger equation for an electron in a hydrogen atom for which the many-electron problem does not arise, so that a direct variational treatment of the problem is possible which can be compared to the analytic solution.

The electronic Schrödinger equation for the hydrogen atom reads: $$\left[ -\frac{\hslash^2}{2m} \nabla^2 -\frac{e^2}{4\pi \epsilon_0} \frac{1}{r}\right] \psi({\bf r}) = E \psi({\bf r})$$ where the second term in the square brackets is the Coulomb attraction potential of the nucleus. The mass $m$ is the reduced mass of the proton--electron system which is approximately equal to the electron mass. The ground state is found at energy $$E=-\frac{m}{\hslash^2} \left( \frac{e^2}{4\pi \epsilon_0} \right)^2 \approx -13.6058~ e{\rm V}$$ and the wave function is given by $$\label{Eq:ExHydr} \psi({\bf r}) = \frac{1}{a_0^{3/2}\sqrt{\pi}} e^{-r/a_0},$$ in which $a_0$ is the Bohr radius, $$a_0 =\frac{4\pi \epsilon_0 \hslash^2}{me^2} \approx 0.529\,18 \unicode{x212B}.$$

Units

When performing a calculation for such an equation, it is convenient to use units such that equations take on a simple form, involving only coefficients of order 1. Standard units in electronic structure physics are so-called atomic units: the unit of distance is the Bohr radius $a_0$, masses are expressed in the electron mass $m_e$ and the charge is measured in unit charges ($e$). The energy is finally given in 'Hartrees' ($E_{\text{H}}$), given by $m_e c^2 \alpha^2$ ($\alpha$ is the fine-structure constant and $m_e$ is the electron mass) which is roughly equal to 27.212 $e{\rm V}$. In these units, the Schrödinger equation for the hydrogen atom assumes the following simple form: $$\left[ -\frac{1}{2} {\bf \nabla}^2 - \frac{1}{r} \right] \psi({\bf r}) = E \psi ({\bf r}). \label{simpschrod}$$

We try to approximate the ground state energy and wave function of the hydrogen atom in a linear variational procedure. We use Gaussian basis functions. For the ground state, we only need angular momentum $l=0$ functions (s-functions) -- they have the form: $$\chi_p (r) = e^{-\alpha_p r^2}$$ centred on the nucleus (which is thus placed at the origin). We have to specify the values of the exponents $\alpha_p$. A large $\alpha_p$ defines a basis function which is concentrated near the nucleus, whereas small $\alpha_p$ characterises a function with a long tail. Optimal values for the exponents $\alpha_p$ can be found by solving the non-linear variational problem including the linear coefficients $C_p$ and the exponents $\alpha_p$. Several numerical methods for solving such non-linear optimisation problem exist and the solutions can be found in textbooks or documentation of quantum chemical software packages.

We shall use known, fixed values of the exponents: $$\begin{aligned} \alpha_1 &= 13.007\,73 \nonumber\\ \alpha_2 &= 1.962\,079 \nonumber\\ \alpha_3 &= 0.444\,529 \\ \alpha_4 &= 0.121\,949\,2, \nonumber \end{aligned}$$ but relax the values of the coefficients $C_p$. The wave function therefore has the form $$\left\langle{\bf r}| \psi \right\rangle= \sum_{p=1}^4 C_p e^{-\alpha_p r^2}$$ with the $\alpha_p$ listed above. We now discuss how to find the best values of the linear coefficients $C_p$. To this end, we need the overlap and Hamiltonian matrix. The advantage of using Gaussian basis functions is that analytic expressions for these matrices can be found. In particular, it is not so difficult to show that the elements of the overlap matrix $\bf S$, the kinetic energy matrix $\bf T$ and the Coulomb matrix $\bf A$ are given by: $$\label{Eq:HydrMat} \begin{aligned} S_{pq} & = \int d^3r e^{-\alpha_p r^2} e^{-\alpha_q r^2} = \left(\frac{\pi}{\alpha_p+\alpha_q}\right)^{3/2}, \\ T_{pq} & = - \frac{1}{2} \int d^3r e^{-\alpha_p r^2} \nabla^2 e^{-\alpha_q r^2} = 3 \frac{\alpha_p \alpha_q \pi^{3/2}}{(\alpha_p +\alpha_q)^{5/2}}, \\ A_{pq} & = - \int d^3r e^{-\alpha_p r^2} \frac{1}{r} e^{-\alpha_q r^2} = -\frac{2 \pi}{\alpha_p +\alpha_q}. \end{aligned}$$ These expressions can be put into a computer program which solves the generalised eigenvalue problem. The resulting ground state energy is $-0.499278$ Hartree, which is amazingly close to the exact value of $-1/2$ Hartree, which is $-13.6058$ eV. We conclude that four Gaussian functions can be linearly combined into a form which is surprisingly close to the exact ground state wave function which is known to have the so-called Slater-type form $\exp(-r/a)$ rather than a Gaussian! This is shown in the figure below (3.2)

Figure (3.2): The best fit of a Gaussian basis set consisting of 1, 2, 3 and 4 basis functions, to the exact Slater solution $e^{-r/a}$ of the hydrogen atom. The mnemonic STO-nG denotes a fit using $n$ Gaussian functions to a Slater Type Orbital (STO).

3.2.3. Exploiting symmetry

We have seen that the solution of a stationary quantum problem using linear variational calculus, in the end, boils down to solving a (generalised) matrix eigenvalue problem. Finding the eigenvalues (and eigenvectors) of an $N\times N$ matrix, or solving a generalised eigenvalue problem, requires a number of floating point operations in the computer proportional to $N^3$. This means that if we double the size of the basis set used in the variational analysis, the computer time goes up by a factor of 8. As it turns out, we are often interested in problems having some symmetry. We shall now briefly sketch how this can be used to significantly reduce the computer time for variational calculations.

In subsection 3.2.1, we considered a problem having a very simple symmetry: replacing $x$ by $-x$ does not change the potential, and therefore the Hamiltonian is insensitive to this transformation. Let us denote the operation $x \rightarrow -x$ by $\mathcal R$. Because flipping the sign of $x$ twice leaves the space invariant, we have $${\mathcal R}^2 = \unicode{120793},$$ where $\unicode{120793}$ is as usual the identity operator. Let us consider the eigenvalues $\lambda$ of this operator. From ${\mathcal R}^2 = \unicode{120793}$ we have that $\lambda^2 = 1$. Therefore, $\lambda=\pm 1$. Furthermore, ${\mathcal R}$ commutes with the Hamiltonian: $${\mathcal R} H - H {\mathcal R} = 0,$$ since the Hamiltonian is not affected by ${\mathcal R}$. We know (or should know!) that if an operator commutes with $H$ we can always find eigenvalues which are eigenvalues of $H$ and of that operator. This means that we can divide the eigenfunctions of $H$ into two classes: one of symmetric eigenfunctions (symmetric meaning having eigenvalue $\lambda = +1$ when acting on it with $\mathcal R$) and one of antisymmetric eigenfunctions ($\lambda = -1$).

Now suppose we construct our variational basis set such that it can be divided into two classes, that of symmetric and that of anti-symmetric basis functions. Let us calculate the inner product of a symmetric and an anti-symmetric eigenfunction. Using antisymmetry of the product of the two may immediately convince you that this vanishes. To illustrate the more general procedure, we consider two eigenfunctions, $\left| \phi_1\right\rangle$ and $\left| \phi_2\right\rangle$ with different eigenvalues $\lambda_1$ and $\lambda_2$ for the symmetry operation ${\mathcal R}$. Then we can write $$\left\langle\phi_1 | {\mathcal R} |\phi_2 \right\rangle= \lambda_1 \left\langle\phi_1 |\phi_2 \right\rangle= \lambda_2 \left\langle\phi_1 |\phi_2 \right\rangle,$$ where we first let ${\mathcal R}$ act on the left, and then on the right function. The fact that $\lambda_1 \neq \lambda_2$ leads to the well-known theorem saying that two eigenvectors of a Hermitian operator with different eigenvalues, are orthogonal (if the eigenvalues are the same, the wave functions are either identical or they can be chosen orthogonal).

The key result is that a similar conclusion can be drawn for the expectation value of the Hamiltonian. $$\left\langle\phi_1 \left| {\mathcal R} H \right| \phi_2 \right\rangle= \lambda_1 \left\langle\phi_1 \left| H \right| \phi_2 \right\rangle= \left\langle\phi_1 \left| H {\mathcal R} \right| \phi_2 \right\rangle = \lambda_2 \left\langle\phi_1 \left| H \right| \phi_2 \right\rangle,$$ which, as $\lambda_1 \neq \lambda_2$ directly gives $$\left\langle\phi_1 | H | \phi_2 \right\rangle= 0.$$

For an orthonormal basis set, we see that, if we would order the basis functions in our set with respect to their eigenvalue of ${\mathcal R}$, the Hamiltonian becomes block-diagonal. For our simple reflection-symmetric example, denoting the symmetric basis functions by $\left| \phi_{ps} \right\rangle$ and the antisymmetric ones by $\left| \phi_{pa}\right\rangle$ where $p$ runs from 1 to $M=N/2$, we have $$H = \left( \begin{array}{ccc|ccc} H_{1s, 1s} & \ldots & H_{1s, Ms} & 0 \ldots & 0 \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ H_{Ms, 1s} & \ldots & H_{Ms, Ms} & 0 \ldots & 0 \\ \hline 0 & \ldots & 0 & H_{1a, 1a} & \ldots & H_{1a, Ma} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & \ldots & 0 & H_{Ma, 1a} & \ldots & H_{Ma, Ma} \end{array}\right)$$ We can diagonalise the two blocks on the diagonal independently. This takes $2 (N/2)^3$ steps (up to a multiplicative constant), which is 4 times less than $N^3$ (up to the same constant) required to diagonalise the full Hamiltonian! If there are additional symmetries, they can be used to reduce the work required even further.

If the basis is non-orthogonal, it still holds that basis functions having different eigenvalues under the symmetry-operator ${\mathcal R}$ are orthogonal and that the matrix elements of the Hamiltonian between them vanishes. This means that the Hamiltonian matrix and the overlap matrix have the same block-diagonal structure. Therefore, the respective generalised eigenvalues for the blocks can be dealt with independently of each other and we achieve the same speed-up.

What we have touched upon is an example of the application of group theory in physics, which is an important topic on its own.