Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

DEMO: Solving H^{+}_2 via Variational Method

Populating the HH and SS Matrices

To determine the variationally optimal energy and wave function, we need to determine the variationally optimal linear coefficients c1c_1 and c2c_2 of the trial wave function. The procedure for this is to:

  1. Posit the trial wave function (done)

  2. Populate the HH and SS matrices

  3. Diagonalize S−1HS^{-1}H.

Here will will populate the 2×22\times 2 matrices.

The Overlap Matrix SS

We start with SS:

S=[S11S12S12S22]\begin{align} S = \begin{bmatrix} S_{11} & S_{12} \\ S_{12} & S_22 \end{bmatrix} \end{align}
(1)

where

S11=⟨1sA∣1sA⟩S12=⟨1sA∣1sB⟩S22=⟨1sB∣1sB⟩\begin{align} S_{11} &= \left\langle 1s_A | 1s_A \right\rangle \\ S_{12} &= \left\langle 1s_A | 1s_B \right\rangle \\ S_{22} &= \left\langle 1s_B | 1s_B \right\rangle \end{align}
(2)

From normalization conditions, we know that

S11=1S22=1\begin{align} S_{11} &= 1 \\ S_{22} &= 1 \end{align}
(3)

We now set about determining S12S_{12}:

S12=⟨1sA∣1sB⟩=∫1sA(r⃗)1sB(r⃗)dr⃗\begin{align} S_{12} &= \left\langle 1s_A | 1s_B \right\rangle \\ &= \int 1s_A(\vec{r})1s_B(\vec{r})d\vec{r} \end{align}
(4)

This integral is somewhat involved but can be solved. Rather than doing so here I just provide the solution

S12=e−R(1+R+R23)S_{12} = e^{-R}\left(1+R+\frac{R^2}{3}\right)
(5)

We note the S12S_{12} is a function of nuclear separation RR. This makes sense as the degree of overlap between the two 1s function will depend on then separation distance between the nuclei. We should expect no overlap when the nuclei are infinitely separated and maximal overlap when the nuclei are sitting on top of each other. Below is a plot of S12(R)S_{12}(R).

From this information we can populate the SS matrix as a function of RR. For example, if we choose R=2.5R=2.5 a0a_0 we get and SS matrix of

S(2.5)=[10.4580.4581]S(2.5) = \begin{bmatrix} 1 & 0.458 \\0.458 & 1 \end{bmatrix}
(6)

Below is a code snippet to populate this matrix.

The Hamiltonian Matrix HH:

The Hamiltonia matrix HH is more challenging. We start with the definition

H=[H11H12H12H22]\begin{align} H = \begin{bmatrix} H_{11} & H_{12} \\ H_{12} & H_{22} \end{bmatrix} \end{align}
(7)

where

H11=⟨1sA∣H^∣1sA⟩H12=⟨1sA∣H^∣1sB⟩H22=⟨1sB∣H^∣1sB⟩\begin{align} H_{11} &= \left\langle 1s_A | \hat{H} | 1s_A \right\rangle \\ H_{12} &= \left\langle 1s_A | \hat{H} | 1s_B \right\rangle \\ H_{22} &= \left\langle 1s_B | \hat{H} | 1s_B \right\rangle \end{align}
(8)

Each of these will need to be determined but we will be able to determine H22H_{22} by analogy to H11H_{11} since all we are doing in swapping the identity of the electron.

H11=⟨1sA∣H^∣1sA⟩=∫1sA(r⃗)(−12∇12−1r1A−1r1B+1R)1sA(r⃗)dr⃗=∫1sA(r⃗)(−12∇12−1r1A)1sA(r⃗)dr⃗+∫1sA(r⃗)(−1r1B+1R)1sA(r⃗)dr⃗=∫1sA(r⃗)(−12∇12−1r1A)1sA(r⃗)dr⃗+∫1sA(r⃗)(−1r1B)1sA(r⃗)dr⃗+1R∫1sA(r⃗)1sA(r⃗)dr⃗=E1s+J+1R\begin{align} H_{11} &= \left\langle 1s_A | \hat{H} | 1s_A \right\rangle \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} - \frac{1}{r_{1B}} + \frac{1}{R} \right) 1s_A(\vec{r})d\vec{r} \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} \right) 1s_A(\vec{r})d\vec{r} + \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}} + \frac{1}{R} \right) 1s_A(\vec{r})d\vec{r} \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} \right) 1s_A(\vec{r})d\vec{r} + \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}} \right) 1s_A(\vec{r})d\vec{r} + \frac{1}{R}\int 1s_A(\vec{r})1s_A(\vec{r})d\vec{r}\\ &= E_{1s} + J + \frac{1}{R} \end{align}
(9)

where we recognize that the first integral is simply the 1s energy of a hydrogen atom and we define the second integral as JJ, what is also called the Coulomb integral. That is

J=∫1sA(r⃗)(−1r1B)1sA(r⃗)dr⃗=−∫1sA(r⃗)1r1B1sA(r⃗)dr⃗\begin{align} J &= \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}} \right) 1s_A(\vec{r})d\vec{r} \\ &= -\int 1s_A(\vec{r})\frac{1}{r_{1B}}1s_A(\vec{r})d\vec{r} \end{align}
(10)

where that last integral becomes 1R\frac{1}{R} because RR, the nuclear separation, is constant w.r.t r⃗\vec{r} and the 1sA1s_A function is normalized.

For the H2+H_2^+ case the remaining integral in JJ can be solved analytically. Doing so yields

J=e−2R(1+1R)−1R\begin{align} J = e^{-2R}\left(1+\frac{1}{R}\right) - \frac{1}{R} \end{align}
(11)

Thus, we get that

H11=E1s+J+1R=−12+e−2R(1+1R)H22=E1s+J+1R=−12+e−2R(1+1R)\begin{align} H_{11} & = E_{1s} + J + \frac{1}{R}= -\frac{1}{2} + e^{-2R}\left(1+\frac{1}{R}\right) \\ H_{22} & = E_{1s} + J + \frac{1}{R}= -\frac{1}{2} + e^{-2R}\left(1+\frac{1}{R}\right) \end{align}
(12)

For H12H_{12}:

H12=⟨1sA∣H^∣1sB⟩=∫1sA(r⃗)(−12∇12−1r1A−1r1B+1R)1sB(r⃗)dr⃗=∫1sA(r⃗)(−12∇12−1r1A)1sB(r⃗)dr⃗+∫1sA(r⃗)(−1r1B+1R)1sB(r⃗)dr⃗=∫1sA(r⃗)(−12∇12−1r1A)1sB(r⃗)dr⃗+∫1sA(r⃗)(−1r1B)1sB(r⃗)dr⃗+1R∫1sA(r⃗)1sB(r⃗)dr⃗=E1sS12+K+S12R\begin{align} H_{12} &= \left\langle 1s_A | \hat{H} | 1s_B \right\rangle \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} - \frac{1}{r_{1B}} + \frac{1}{R} \right) 1s_B(\vec{r})d\vec{r} \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} \right) 1s_B(\vec{r})d\vec{r} + \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}} + \frac{1}{R} \right) 1s_B(\vec{r})d\vec{r} \\ &= \int 1s_A(\vec{r}) \left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} \right) 1s_B(\vec{r})d\vec{r} + \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}}\right)1s_B(\vec{r})d\vec{r} + \frac{1}{R}\int 1s_A(\vec{r}) 1s_B(\vec{r})d\vec{r} \\ &= E_{1s}S_{12} + K + \frac{S_{12}}{R} \end{align}
(13)

where we recognize that (−12∇12−1r1A)1sB=E1s1sB\left(-\frac{1}{2}\nabla_1^2 - \frac{1}{r_{1A}} \right) 1s_B = E_{1s}1s_B, S12=⟨1sA∣1sB⟩S_{12} = \langle 1s_A | 1s_B\rangle, and have defined the exhange integral, KK, to be

K=∫1sA(r⃗)(−1r1B)1sB(r⃗)dr⃗=−∫1sA(r⃗)1r1B1sB(r⃗)dr⃗=−e−R(1+R)\begin{align} K &= \int 1s_A(\vec{r}) \left(- \frac{1}{r_{1B}} \right) 1s_B(\vec{r})d\vec{r} \\ &= -\int 1s_A(\vec{r})\frac{1}{r_{1B}}1s_B(\vec{r})d\vec{r} \\ &= -e^{-R}\left(1+R\right) \end{align}
(14)

where the last step I have used the solution for the simple case of H2+_2^+ but have not shown it to be true.

We can now populate the HH matrix as

H=[E1s+J+1RE1sS12+K+S12RE1sS12+K+S12RE1s+J+1R]=(E1s+1R)[1S12S121]+[JKKJ]\begin{align} H &= \begin{bmatrix} E_{1s} + J + \frac{1}{R} & E_{1s}S_{12} + K + \frac{S_{12}}{R} \\ E_{1s}S_{12} + K + \frac{S_{12}}{R} & E_{1s} + J + \frac{1}{R}\end{bmatrix} \\ &= \left(E_{1s} + \frac{1}{R}\right)\begin{bmatrix} 1 & S_{12} \\ S_{12} & 1\end{bmatrix} + \begin{bmatrix} J & K \\ K & J \end{bmatrix} \end{align}
(15)

For a nuclear separation value of 2.5 Bohr, for example, this yields

H(R=2.5)=[−0.491−0.333−0.333−0.491]H(R=2.5) = \begin{bmatrix} -0.491 & -0.333 \\ -0.333 & -0.491\end{bmatrix}
(16)

Below is the code to compute this.

Combined Solution

Now that we have equations for HH and SS in terms of RR (the nuclear separation; S12S_{12}, JJ, and KK are all functions of RR), we can compute the eigenvalues and eigenvectors of S−1HS^{-1}H for a given nuclear separation distance.

Let’s choose 2.5 Bohr, for example, in which case we have

S(2.5)=[10.4580.4581]H(2.5)=[−0.491−0.333−0.333−0.491]S(2.5)−1H(2.5)=[−0.428−0.138−0.138−0.428]\begin{align} S(2.5) &= \begin{bmatrix} 1 & 0.458 \\0.458 & 1 \end{bmatrix} \\ H(2.5) &= \begin{bmatrix} -0.491 & -0.333 \\ -0.333 & -0.491\end{bmatrix} \\ S(2.5)^{-1}H(2.5) &= \begin{bmatrix} -0.428 & -0.138 \\ -0.138 & -0.428\end{bmatrix} \end{align}
(17)

Computing the eigenvalues of S(2.5)−1H(2.5)S(2.5)^{-1}H(2.5) yields:

E(2.5)=−0.565 and −0.291E(2.5) = -0.565\text{ and } -0.291
(18)

The ground state energy is the first (minimum) of these. The corresponding eigenvectors (or vectors of variational coefficients) are

v=[−0.7070.707−0.707−0.707]v = \begin{bmatrix} -0.707 & 0.707 \\ -0.707 & -0.707\end{bmatrix}
(19)

where these are listed as column vectors meaning the that coefficients for the ground state energy are

c1=−0.707c2=−0.707\begin{align} c_1 &= -0.707 \\ c_2 &= -0.707 \end{align}
(20)

Due to the symmetry of this problem and normalization requirements, it turns out that the coefficients must always be

c1=±12(1−S12)c2=±12(1−S12)\begin{align} c_1 &= \pm \frac{1}{\sqrt{2(1-S_{12})}} \\ c_2 &= \pm \frac{1}{\sqrt{2(1-S_{12})}} \end{align}
(21)

That is the must always be the same magnitude (symmetry) and their values are dictated by normalization (which depends on the overlap of the two basis functions).

The minimum energy state is the linear combination in which the coefficients both have the same sign. The higher energy (antibonding orbital) has the coefficients with opposing sign. These can be described as

ϕ+=12(1−S12)(1sA+1sB)ϕ−=12(1−S12)(1sA−1sB)\begin{align} \phi_+ = \frac{1}{\sqrt{2(1-S_{12})}}\left(1s_A + 1s_B\right) \\ \phi_- = \frac{1}{\sqrt{2(1-S_{12})}}\left(1s_A - 1s_B\right) \end{align}
(22)

Solution as a Function of RR

In order to investigate things such as the predicted bond length and bonding energy of this model of the H2+_2^+ molecule, we must perform the above calculation (population of the HH and SS matrices and diagonalization of the S−1HS^{-1}H matrix) for various RR values. The bond distance will be the value of RR at the minimum energy and the bonding energy will be the difference between the energy at that position and the energy at infinite separation of the nuclei (-0.5 Hartree in this case).

Below I make a plot of energies as a function of RR for this model of the H2+_2^+ molecule.

We see from the plot and code output above that the predicted bond distance for this model of H2+_2^+ is 2.493 Bohr with a bond energy of -0.06483 Hartree (-40.7 kcal/mol).

These are not actually in great agreement with the “exact” values of 1.999 Bohr and -0.1 Hartree. But the simple approach outline here does give qualitative agreement and is a useful model for generalizing to other molecules.

We could get better agreement in a number of ways. These include:

  1. Increasing the number of basis functions

  2. Improving the basis functions (either by modifying them or including additional variational parameters)

Orbital/Electron Densities

Appendix C: Python and Interactive Demos
Demo: Hartree-Fock with PySCF
Appendix C: Python and Interactive Demos
DEMO: HF calculations on small organic molecules by PySCF