PHY.K02UF Molecular and Solid State Physics

Molecular orbitals of the molecular ion H2+

The molecular ion H2+ consists of one electron and two protons. The molecular orbital Hamiltonian in this case is,

\[ \begin{equation} H_{\text{mo}}^{H_2^+}= - \frac{\hbar^2}{2m_e}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_A|}- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_B|} , \end{equation} \]

where $\vec{r}_A$ and $\vec{r}_B$ are the positions of the two protons. Consider a linear combination of the two 1s orbitals, $\psi_{\text{mo}}=c_1\phi_{\text{1s}}(\vec{r}-\vec{r}_A) +c_2\phi_{\text{1s}}(\vec{r}-\vec{r}_B)$. The time independent Schrödinger equation is,

\begin{equation} H_{\text{mo}}\psi_{\text{mo}}=E\psi_{\text{mo}} . \end{equation}

Multiply the Schrödinger equation from the left by each of the atomic orbitals and integrate over all space. This results in a set of 2 algebraic equations called the Roothaan equations.

\[ \begin{equation} \begin{matrix} \langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi_{\text{1s}}(\vec{r}-\vec{r}_A)|\psi_{\text{mo}} \rangle \\ \langle \phi_{\text{1s}}(\vec{r}-\vec{r}_B) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi_{\text{1s}}(\vec{r}-\vec{r}_B)|\psi_{\text{mo}} \rangle\\ \end{matrix} \end{equation} \]

By substituting in the form for $\psi_{\text{mo}}$ from above, the Roothaan equations can be written in matrix form,

\[ \begin{equation} \left[ \begin{matrix} H_{11} & H_{12} \\ H_{12} & H_{11} \end{matrix} \right] \left[ \begin{matrix} c_1 \\ c_2 \end{matrix} \right] = E \left[ \begin{matrix} S_{11} & S_{12} \\ S_{12} & S_{11} \end{matrix} \right]\left[ \begin{matrix} c_1 \\ c_2 \end{matrix} \right], \end{equation} \]

where

\[ \begin{equation} \begin{matrix} H_{11}=\langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A)|H_{\text{mo}}|\phi_{\text{1s}}(\vec{r}-\vec{r}_A)\rangle, \\ H_{12}=\langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A)|H_{\text{mo}}|\phi_{\text{1s}}(\vec{r}-\vec{r}_B)\rangle, \\ S_{11}=\langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A)|\phi_{\text{1s}}(\vec{r}-\vec{r}_A)\rangle=1, \\ S_{12}=\langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A)|\phi_{\text{1s}}(\vec{r}-\vec{r}_B)\rangle. \end{matrix} \end{equation} \]

From experience with $2\times 2$ matrices of this form, we know that the eigenvectors of both the Hamiltonian matrix and the overlap matrix are,

\[ \begin{equation} \left[ \begin{matrix} c_1 \\ c_2 \end{matrix} \right] = \left[ \begin{matrix} 1 \\ 1 \end{matrix} \right],\left[ \begin{matrix} 1 \\ -1 \end{matrix} \right] \end{equation} \]

It is easy to check that these are the correct eigenvectors but letting $\textbf{H}$ and $\textbf{S}$ operate on them. The energies are,

\[ \begin{equation} E_{+}= \frac{H_{11}+ H_{12}}{1+ S_{12}}, \qquad E_{-}= \frac{H_{11}- H_{12}}{1- S_{12}}. \end{equation} \]

The normalized molecular orbitals are,

\[ \begin{equation} \psi_{\pm}= \frac{1}{\sqrt{2}}\left(\phi_{\text{1s}}(\vec{r}-\vec{r}_A) \pm \phi_{\text{1s}}(\vec{r}-\vec{r}_B)\right) \end{equation} \]

Both $H_{11}$ and $H_{12}$ will turn out to be negative so the bonding orbital $\psi_+$, has a lower energy than the antibonding orbital $\psi_-$.

Determining the matrix elements $H_{11}$, $H_{12}$, $S_{11}$, and $S_{12}$

The 1s atomic orbital has the form,

\begin{equation} \phi_{1s}= \frac{1}{\sqrt{\pi a^3_0}}\exp\left(-\frac{r}{a_0}\right). \end{equation}

$a_0^{3/2}\phi_{1s}$

$\frac{x}{a_0}$

For normalized wavefunctions, the diagonal elements of the overlap matrix are equal to 1, $S_{11}=\langle \phi_{\text{1s}}(\vec{r}-\vec{r}_A)|\phi_{\text{1s}}(\vec{r}-\vec{r}_A)\rangle=1$.

The overlap integral $S_{12}$ of two 1s orbitals located at positions $\vec{r}_1$ and $\vec{r}_2$ is,

\begin{equation} S_{12}=\langle \phi_{1s}(\vec{r}-\vec{r}_1)|\phi_{1s}(\vec{r}-\vec{r}_2)\rangle . \end{equation}

For a hydrogen molecule, $\vec{r}_1=-0.38\,\hat{x}$ Å and $\vec{r}_2=0.38\,\hat{x}$ Å. Below $\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$ is plotted along the $x$-axis and along the $y$-axis.

$a_0^3\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$

$\frac{x}{a_0}\,,\frac{y}{a_0}$

The code below uses a Monte-Carlo method to integrate $\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$ and calculate $S_{12}$.

Press the 'Execute' button a few times and notice that the answer keeps changing. By doing this you can estimate the error in the calculation. The error should decrease like $1/\sqrt{N}$ where $N$ is the number of random numbers chosen. By setting $x_1=x_2=0$ in the code and pressing 'Execute', you calculate $S_{11} = \langle \phi_{1s}(\vec{r})|\phi_{1s}(\vec{r})\rangle$ which should equal 1 if the wave function is properly normalized.

The diagonal Hamiltonian matrix element with two 1s orbitals located at position $\vec{r}_1$ is,

\begin{equation} H_{11}=\Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{11}=\Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}\right|\phi_{1s}(\vec{r}-\vec{r}_1) \Big \rangle + \Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

The wave function $\phi_{1s}(\vec{r}-\vec{r}_1)$ is an eigenfunction of the atomic orbital Hamiltonian in the first term $H\phi_{1s}(\vec{r}-\vec{r}_1) = E_1 \phi_{1s}(\vec{r}-\vec{r}_1)$, so the first term is easily evaluated,

\begin{equation} H_{11}=E_1 - \Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left| \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

The second term has a singularity at $\vec{r}_2$ which makes it difficult to evaluate numerically. We break the second term into an integral over a spherical volume of radius $\delta$ centered around $\vec{r}_2$ and a second integral outside that volume.

\begin{equation} H_{11}=E_1 - \int\limits_{|\vec{r}-\vec{r}_2| < \delta} \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}(\vec{r}-\vec{r}_1) d^3r - \int\limits_{|\vec{r}-\vec{r}_2| > \delta} \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}(\vec{r}-\vec{r}_1) d^3r. \end{equation}

Close to $\vec{r}_2$, $\exp\left(\frac{ -|\vec{r}-\vec{r}_1|}{a_0}\right)\approx \exp\left(\frac{ -|\vec{r}_2-\vec{r}_1|}{a_0}\right)$. Using this approximation, the first integral which includes the singularity can be performed analytically for small $\delta$.

\begin{equation} H_{11}=E_1 - \frac{e^2\delta^2}{2\pi\epsilon_0a_0^3} \exp(-2|\vec{r}_1-\vec{r}_2|/a_0) - \int H(|\vec{r}-\vec{r}_2|-\delta ) \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}(\vec{r}-\vec{r}_1) d^3r, \end{equation}

The second integral integrates over all space but a Heaviside step function has been introduced. $H(|\vec{r}-\vec{r}_2|-\delta ) = 0$ for $|\vec{r}-\vec{r}_2| < \delta$ and is 1 otherwise. The second integral contains no singularity and can be evaluated numerically.

$H(|\vec{r}-\vec{r}_2|-\delta )\frac{a_0^3e\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_1)}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|}$

$\frac{x}{a_0}$

The integrand of the matrix element plotted along the $x$-axis for $\delta = a_0/10$.

The code below uses a Monte-Carlo method to and calculate $H_{11}$.

Similarly, the off-diagonal Hamiltonian matrix element with two 1s orbitals located at positions $\vec{r}_1$ and $\vec{r}_2$ is,

\begin{equation} H_{12}=\Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{12}=\Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|}\right|\phi_{1s}(\vec{r}-\vec{r}_2) \Big \rangle + \Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left|- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \right|\phi_{1s}(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

The wave function $\phi_{1s}(\vec{r}-\vec{r}_2)$ is an eigenfunction of the atomic orbital Hamiltonian in the first term $H\phi_{1s}(\vec{r}-\vec{r}_2) = E_1 \phi_{1s}(\vec{r}-\vec{r}_2)$, so the first term is easily evaluated,

\begin{equation} H_{12}=E_1S_{12} - \Big \langle \phi_{1s}(\vec{r}-\vec{r}_1)\left| \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \right|\phi_{1s}(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

The second term has a singularity at $\vec{r}_1$ which makes it difficult to evaluate numerically. We break the second term into an integral over a spherical volume of radius $\delta$ centered around $\vec{r}_1$ and a second integral outside that volume.

\begin{equation} H_{12}=E_1S_{12} - \int\limits_{|\vec{r}-\vec{r}_1| < \delta} \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}(\vec{r}-\vec{r}_2) d^3r - \int\limits_{|\vec{r}-\vec{r}_1| > \delta} \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}(\vec{r}-\vec{r}_2) d^3r. \end{equation}

Close to $\vec{r}_1$, $\exp\left(\frac{ -Z|\vec{r}-\vec{r}_2|}{a_0}\right)\approx \exp\left(\frac{ -Z|\vec{r}_1-\vec{r}_2|}{a_0}\right)$. Using this approximation, the first integral which includes the singularity can be performed analytically for small $\delta$.

\begin{equation} H_{12}=E_1S_{12} - \frac{Z^4e^2}{\pi a_0^2\epsilon_0} (a_0-\exp (-\delta Z/a_0) (a_0+\delta Z)) \exp(-Z|\vec{r}_1-\vec{r}_2|/a_0) - \int H(|\vec{r}-\vec{r}_1|-\delta ) \phi_{1s}(\vec{r}-\vec{r}_1) \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}(\vec{r}-\vec{r}_2) d^3r, \end{equation}

The second integral integrates over all space but a Heaviside step function has been introduced. $H(|\vec{r}-\vec{r}_1|-\delta ) = 0$ for $|\vec{r}-\vec{r}_1| < \delta$ and is 1 otherwise. The second integral contains no singularity and can be evaluated numerically.

$H(|\vec{r}-\vec{r}_1|-\delta )\frac{a_0^3e\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}$

$\frac{x}{a_0}$

The integrand of the matrix element plotted along the $x$-axis for $\delta = a_0/10$.

The code below uses a Monte-Carlo method to and calculate $H_{12}$.

A Matlab script to calculate the overlap matrix and the Hamiltonian matrix is calcHamiltonian.m. This script uses the atomic orbitals defined in atorbit.m. Matlab files to plot the molecular orbitals and the bond potential can be found here.

To calculate the bond potential, the energy of the ground-state wavefunction is evaluated with the Hamiltonian that includes the proton-proton interaction energy.

\[ \begin{equation}\label{eq:helec} \begin{matrix} H^{H_2^+}=- \frac{\hbar^2}{2m_e}\nabla - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_A|}- \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_B|} +\frac{e^2}{4\pi\epsilon_0 |\vec{r}_A-\vec{r}_B|}. \end{matrix} \end{equation} \]

For a given distance between the protons, $|\vec{r}_A-\vec{r}_B|$, the energy is evaluated by calculating,

\[ \begin{equation} E= \frac{\langle \psi_+ |H^{H_2^+}|\psi_+\rangle}{\langle \psi_+|\psi_+\rangle}. \end{equation} \]

A plot of the energy as a function of interatomic distance is shown below.

$E\text{ [eV]}$

$r$ [Å]

The experimental value for the bond length is $r_e = 1.06$ Å and the dissociation energy is 2.79 eV. To get a more accurate calculation, a larger basis set is needed. The plot below shows the results of this calculation for larger basis sets. The calculation labeled 'Large Basis Set' uses functions beyond just the atomic orbitals which are needed for an accurate calculation. Accurate calculations are computationally intensive.

$E\text{ [eV]}$

$r$ [Å]