This post is the first part in an upcoming series explaining how to perform electronic structure calculations. The focus of this series is how to set-up an electronic structure calculation from scratch. We are not going to use any of the available quantum chemical codes, but write our very own code. This way, we will obtain a better understanding of how quantum chemical calculations are being performed. The way I am going to present the topic is by first providing the underlying theory and then further explaining the theory by giving a set of examples.

The outline of this series is as follows: We start of with explaining some basics such as the (linear) variational principle, the Hartree product and the Slater determinant. From there on, we derive the Hartree-Fock equations and show how to set-up a simple calculation for the H_{2} molecule. Next, we explain about basis sets, Gaussian integrals and the Hartree-Fock self-consistent field procedure. Finally, we talk about atomic and molecular orbitals, the construction of MO diagrams and the visualization of orbitals.

For this series, I am going to assume that you have a basic understanding of quantum mechanics. You know what a wave function is and its statistical interpretation. Furthermore, I assume that you have passed introductory courses in calculus and linear algebra.

Given a normalized wave function $\lvert\phi\rangle$ that satisfies the appropriate boundary conditions, then the expectation value of the Hamiltonian is an upper bound to the exact ground state energy. In other words, the lowest energy we can find, gives us the best wave function. Thus,

if $ \langle \phi \rvert \phi \rangle=1 $, then $ \langle \phi \lvert \hat{H} \rvert \phi \rangle \geq E_{0} $.

The above integral can be divided into two seperate integrals of the following forms $$ I_1 = - \left( \frac{4 \pi N^{2}}{2} \right) \int_0^\infty r^2 \exp \left( -\alpha r^2 \right) \left[ r^{-2} \frac{\textrm{d}}{\textrm{d} r} \left( r^{2} \frac{\textrm{d}}{\textrm{d} r} \right) \exp \left( -\alpha r^2 \right) \right] \textrm{d}r $$ and $$ I_2 = - \left( \frac{4 \pi N^{2}}{2} \right) \int_0^\infty r^2 \exp \left( -\alpha r^2 \right) \frac{1}{r} \exp \left( -\alpha r^2 \right) \textrm{d}r .$$ Since the wavefunction has no angular dependencies, the azimuthal ($\theta$) and polar ($\varphi$) parts have already been integrated in the above equations. These give the factor $4\pi$. To solve the first integral ($I_{1}$), we need to apply the Laplacian. This gives $$ I_1 = - \left( \frac{4 \pi N^{2}}{2} \right) \int_0^\infty \exp \left( -\alpha r^2 \right) \left[ -6 \alpha^{2} + 4 \alpha^{2}r^{4} \right] \exp \left( -\alpha r^2 \right) \textrm{d}r .$$ Slightly rearranging this equation gives us $$ I_1 = - 2 \pi N^{2} \int_0^\infty \left( 4 \alpha^{2}r^{4} -6 \alpha^{2} \right) \exp \left( -2 \alpha r^2 \right) \textrm{d}r .$$ This integral can be solved by performing integration by parts, but here we will use just use the following general formulas: $$\int_{0}^{\infty} \textrm{d}r \, r^{2m} \exp \left(-\alpha r^{2} \right) = \frac{(2m)!\pi^{1/2}}{2^{2m+1}m!\alpha^{m+1/2}} $$ and $$\int_{0}^{\infty} \textrm{d}r \, r^{2m+1} \exp \left(-\alpha r^{2} \right) = \frac{m!}{2\alpha^{m+1}} .$$ Applying these formulas for $I_1$ yields $$I_{1} = -2 \pi N^{2} \left[ 4 \alpha^{2} \frac{4! \pi^{1/2}}{2^5 2! (2 \alpha)^{5/2}} - 6 \alpha \frac{2!\pi^{1/2}}{2^3 1! (2\alpha)^{3/2}} \right] = N^{2} \frac{3 \pi^{3/2}} {4 \sqrt{2\alpha}} $$ and for $I_{2}$: $$I_{2} = -4 \pi N^{2} \int_{0}^{\infty} r \exp \left( -2 \alpha r^{2} \right) \textrm{d}r = -4 \pi N^{2} \frac{1}{4 \alpha} = \frac{-\pi}{\alpha} N^{2} .$$ Summing these two parts gives $$\langle \phi \lvert H \rvert \phi \rangle = I_1 + I_2 = N^2 \left( \frac{3 \pi ^{3/2}}{4 \sqrt{2 \alpha}} - \frac{\pi}{\alpha} \right) .$$ Using the same procedure, we can also evaluate the integral $ \langle \phi \rvert \phi \rangle $. $$ \langle \phi \rvert \phi \rangle = 4 \pi N^{2} \int_{0}^{\infty} r^2 \exp \left( -2 \alpha r^2 \right) \textrm{d}r = N^{2} \frac{\pi^{3/2}}{2 \sqrt{2}\alpha^{3/2}} $$ In order to evaluate the expectation value for the energy $E$ for the trial wave function $\lvert \varphi \rangle$, we need to evaluate $$ E = \frac{\langle \phi \lvert \hat{H} \rvert \phi \rangle }{\langle \phi \rvert \phi \rangle} = \frac{3 \alpha}{2} - \frac{2 \sqrt{2} \alpha ^{1/2}}{\pi ^{1/2}} .$$ The reason we need to divide by $\langle \phi \rvert \phi \rangle$ is because $\lvert \phi \rangle$ is not normalized yet. If $\lvert \phi \rangle$ would be normalized, then the integral $\langle \phi \rvert \phi \rangle$ would be equal to one and thus drop out. The best value for $\alpha$ can then be found by differentiating $E$ with respect to $\alpha$ and finding that value for $\alpha$ where the differential is zero. $$ \frac{\textrm{d}E}{\textrm{d}\alpha} = 0 .$$ Thus $$ \frac{\textrm{d}E}{\textrm{d}\alpha} = \frac{3}{2} - \frac{\sqrt{2}}{\pi ^{1/2} \alpha^{1/2}} $$ and $$ \alpha = \frac{8}{9 \pi}.$$ This gives us the following upper bound for the energy $E$: $$E = \frac{3}{2} \frac{8}{9 \pi} - \frac{2 \sqrt{2} \left( \frac{8}{9 \pi} \right)^{1/2}}{\pi ^{1/2}} = -\frac{4}{3 \pi} \approx -0.4244 .$$

The exact energy for the 1s orbital in hydrogen is -0.5 HT with the corresponding wave function $\lvert \phi \rangle = \frac{1}{\sqrt{\pi}} \exp \left( - r \right)$.

Although the exact answer to the energy of the 1s orbital in hydrogen is known, let us assume that we do *not* know it and would like to obtain it. One way is to pick various trial wave functions with one independent variable and optimize that variable. The lowest energy we are going to find gives us then the best approximation of the ground state wave function. This approach would be quite tedious (and how would you know what trial wave function to choose in the first place anyway?).

It would be nice if we could have just a single trial wave function that has a lot of 'variational freedom'. Here, I propose to use a trial wave function that is a linear set of fixed wave function like so:

$$ \lvert \phi \rangle = \sum_{i=1}^N c_{i} \lvert \psi_{i} \rangle $$

If we can assume that the wave functions $\lvert \psi_i \rangle$ are orthonormal in the sense that $\langle \psi_i \lvert \psi_j \rangle = \delta_{ij}$, then the energy is

$$E = \sum_{ij} \langle \psi_i \lvert \hat{H} \rvert \psi_j \rangle = \sum_{ij} c_{i}c_{j} H_{ij}$$

and we are tasked with finding the best set of coefficients $c_{i}$ that minimizes $E$. It are exactly these coefficients (and accompanying wave functions) that provide us with the necessary variational freedom. (Note that $H_{ij} = \langle \psi_i \lvert \hat{H} \rvert \psi_j \rangle$)

The procedure to find this set is by using Langrange's method of undetermined multiplies. Herein

$$L \left( c_1, c_2, \cdots, E \right) = \langle \phi \lvert \hat{H} \rvert \phi \rangle - E \left(\langle \phi \lvert \phi \rangle - 1 \right) $$ and we need to find that particular set of $c_k$ where

$$ \frac{\partial L}{\partial c_{k}} = 0 $$

for each $k$.

For each $k$, this gives the following equation

$$ \sum_{j} c_j H_{kj} + \sum_i c_i H_{ik} - 2 E c_{k} = 0 $$.

Since $\langle \psi_a \lvert \hat{H} \rvert \psi_b \rangle = \langle \psi_b \lvert \hat{H} \rvert \psi_a \rangle = H_{ab} = H_{ba}$, we can simplify the above equation to

$$ \sum_{j} H_{ij} c_j - E c_{i} = 0. $$

The above expression is just the matrix expression

$$ H\vec{c} = E\vec{c} .$$

This last equation looks probably very familiar to you and relates to the matrix eigenvalue problem. It is a very typical problem in linear algebra. To solve the matrix eigenvalue problem, one needs to diagonalize the $H$ matrix. This gives a set of eigenvalues and corresponding eigenvalues. The eigenvector corresponding to the lowest eigenvalue is the best set of $c_k$ we are interested in.

In order to find the best energy and approximation to the ground state wave function, we need to diagonalize the H-matrix. This can be easily done by hand, or using the very handy Wolfram Alpha website. The result we obtain is $$D = \begin{bmatrix} 1.09554 & 0 \\ 0 & -0.378614 \end{bmatrix}$$ for the diagonal (eigenvalue) matrix with corresponding $$X = \begin{bmatrix} -0.0126156 & -0.99992 \\ 0.99992 & -0.0126156 \end{bmatrix}$$ for the eigenvector matrix. The first column of $X$ contains the coefficients corresponding to the first eigenvalue (energy) of the $D$-matrix and the second column of $X$ for the second eigenvalue. Since the second eigenvalue is lower ($-0.378614 < 1.09554$), the second column of $X$ contains the coefficients we are looking for. Thus, the best approximation to the ground state wave function is: $$\lvert \phi \rangle = c_{1} \lvert \psi_{1} \rangle + c_{2} \lvert \psi_{2} \rangle = -0.99992 \lvert \psi_{1} \rangle - 0.0126156 \lvert \psi_{2} \rangle.$$ The corresponding energy for this wave function is -0.378614 HT. This result is actually worse than what we obtained using the procedure of example 1 (which was -0.4244 HT). The reason for this is that the variational freedom using a linear combination of only two wave functions is actually less than one where we can vary a coefficient inside an exponent. To get better results, we would need to use a larger set of wave functions in the linear combination.

In Example 2 it was shown how we can obtain the best energy given a trial wave function which is a linear combination of a set of wave functions. Let us call these wave functions in the linear combination the basis functions. In order for the procedure in example 2 to work, the basis functions need to be orthonormal to each other. Finding a set of orthonormal functions by hand is somewhat tedious (try it for yourself if you do not believe me!). It would therefore be very handy if we could apply the above methodology (albeit somewhat changed) using a set of non-orthonormal basis functions.

I already showed you how you can use Lagrange's method of undetermined multipliers in order to derive the matrix equation $Hc = Ec$. When the basis functions are no longer orthonormal, this matrix equation changes slightly and becomes

$$ H\vec{c} = ES\vec{c}.$$

Herein, $S$ is termed the overlap matrix where each term in the matrix is

$$ S_{ij} = \langle \psi_{i} \lvert \psi_{j} \rangle.$$

The introduction of the overlap matrix $S$ makes the matrix diagonalization a tiny bit more complicated. There are several procedures for obtaining the eigenvalues and -vectors. The procedure we are going to use follows the canonical orthogonalization scheme. It works as follows

- Construct the overlap matrix $S$.
- Diagonalize the overlap matrix $S$ to obtain the eigenvalues of the overlap matrix and its corresponding eigenvector matrix $U$.
- Construct a transformation matrix $X = U s^{-1/2}$.
- Transform the Hamiltonian matrix to the canonically orthonormalized basis set $H' = X^{\dagger} H X$.
- Diagonalize the transformed Hamiltonian matrix to obtain $D$ and $C'$ (where C' are the transformed coefficients).
- Obtain the coefficients of the original basis by calculating $C = XC'$.

The above probably looks a bit complicated, but really, it is not. Let me elaborate by showing you another example.

(

(

This does not necessarily imply that the underlying contributions of the basis functions are the same. If you would examine the coefficients and the corresponding basis functions, you would notice a particular difference. The contribution of each basis functions in each set is visualized in the graph below. We can clearly see that the contributions differ.

If we now sum up the contributions of the two basis functions for each linear combination, we note that the two wave functions in example 2 and this example are in principle the same.

Thus, we obtain the same result for each wave function.

Although we have now the tools to find the best approximation of the ground state energy and its corresponding wave function using any set of basis functions, the result so far was not very satisfactory. Using only a single trial wave function as for instance given in example 1 gave better results. In order to improve on the previous results, we need to expand our set of basis functions. A larger basis set should in principle give a better result.

Further improvement by expanding the basis set to 10 basis functions gives a result of -0.498481 HT.

In this blogpost I have shown you that in order to find the best approximation to the ground state energy, the variational principle can be employed. By introducing some kind of variety in a trial wave function and then minimizing its corresponding energy, gives the best approximation to the ground state energy can be found.

In the linear variational method, the variational freedom is provided by the specific set of basis functions used. In this method, the best coefficients in the linear combination of said basis functions that minimize the energy can be found. Here, I have shown that a sufficiently large set of basis functions composed of a Gaussian part (the exponent) multiplied with a polynomial of increasing order is able to approach the true ground state energy of the hydrogen atom.

Finally, I would like to acknowledge Bart Zijlstra for thoroughly proofreading this post and pointing out the mistakes.

If you have questions or comments, feel free to drop a line!
Like what you read? Share this page with your friends and colleagues.