Water ($H_2 O$) consists of one oxygen atom and two hydrogen atoms. The molecular orbital Hamiltonian in this case is,
\begin{equation} H_{\text{mo}}= - \frac{\hbar^2}{2m_e}\nabla^2 - \frac{8e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_O|} - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_{H_1}|} - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_{H_2}|} , \end{equation}where $\vec{r}_O$, $\vec{r}_{H_1}$, $\vec{r}_{H_2}$ are the positions of the oxygen and the hydrogen atoms. We consider a trial wavefunction that is a linear combination of 7 atomic orbitals. These are the 1s, 2s and 2p orbitals of the oxygen atom and the 1s orbitals of the hydrogen atoms.
The following formulas were implemented in python in order to calculate the molecular orbitals.
Atomic orbital wavefunction :
\begin{equation} \psi_{nlm} (r,\theta,\phi) = R_{nl}(r) \, Y_{lm}(\theta, \phi) \end{equation}Radial part of the wavefunction:
\begin{equation} R_{nl}(r) = \sqrt{\left(\frac{2 Z}{n a_0}\right)^3 \frac{(n - l -1)!}{2 n (n+l)!}} \, e^{-\rho / 2} \rho^l \, L_{n - l -1}^{2l + 1}(\rho) \\ \rho = \frac{2 Z r}{n a_0} \end{equation}Elements of the hamiltonian matrix:
\begin{equation} \langle \phi_i^c (\vec{r} - \vec{r}_c) \, \vert \, \hat{H} \, \vert \, \phi_j^b(\vec{r}) \rangle = \langle \phi_i^c (\vec{r} - \vec{r}_i) \, \vert \, \frac{\hat{P}^2}{2 m_e} + \sum_a \hat{V}_a \, \vert \, \phi_j^b(\vec{r}) \rangle \\ \\ = \langle \phi_i^c (\vec{r} - \vec{r}_c) \, \vert \, \frac{\hat{P}^2}{2 m_e} + \hat{V}_b + \sum_{a \neq b} \hat{V}_a \, \vert \, \phi_j^b(\vec{r}) \rangle\\ \\ = \langle \phi_i^c (\vec{r} - \vec{r}_c) \, \vert \, \frac{\hat{P}^2}{2 m_e} + \hat{V}_b \, \vert \, \phi_j^b(\vec{r}) \rangle + \langle \phi_i^c (\vec{r} - \vec{r}_c) \, \vert \, \sum_{a \neq b} \hat{V}_a \, \vert \, \phi_j^b(\vec{r}) \rangle\\ \\ = S_{ij}^{ab} E_j^b + \sigma_{ij}^{cb} \end{equation}Here $\phi_i^a$ is the $i$-th atomic orbital of atom $a$, $\vec{r}_c$ is the position of atom $c$, $\hat{V}_a$ is the potential operator caused by atom $a$, $\hat{P}$ is the momentum operator, $S_{ij}^{ab}$ is the overlap integral of orbital $i$ around atom $a$ and orbital $j$ around atom $b$, $E_j^b$ is the energy of orbital $j$ around atom $b$ and $\sigma_{ij}^b$ is the integral of $\psi_i^c$ and $\psi_j^b$ over all potential operators except the one belonging to atom $b$. The coordinate system is centered around atom $b$ in this case. This way of calculating the hamiltonian matrix elements is described in more detail here .
The functions are saved in the file mo_functions.py .
The file mo.py is used to calculate and save the solutions to the roothaan equations as well as structural information about the molecule as .txt files.
orbitalplot2d.py and orbitalplot3d.py are dedicated to visualizing the molecular orbitals.
The program was written in python 2.7. The python distribution was anaconda with spyder 3.2.8 as the programming environment.
The following packages are needed in order to run the code:
numpy 1.11.3 (numerical calculations, vectorization)
scipy 1.1.0 (numerical integration and differentiation)
matplotlib 1.5.1 (2d plots)
mayavi 4.5.0 (3d isosurface plots)
The file mo_env.txt can be used to recreate the exact environment the program was written in.
The first step ist to open mo.py and to define the molecule.
This is done by first creating the atoms of the molecule with the atom class, for example a hydrogen atom:
H = atom( [ [$n_1$, $l_1$, $m_1$], [$n_2$, $l_2$, $m_2$], ... , [$n_k$, $l_k$, $m_k$] ], $Z$ , [$x_0$, $y_0$, $z_0$], symbol = "H" )
Here $n_i$, $l_i$ and $m_i$ are the quantumnumbers of the atomic orbitals, $Z$ is the atomic number of the atom
and [$x_0$, $y_0$, $z_0$] is the position of the atom in carthesian coordinates. The optional parameter symbol is used for labeling in the 2d plots.
The molecule class is used to create a molecule object:
H2O = molecule([H1, H2, O], "H2O")
The atom objects H1, H2 and O aswell as the name of the molecule are passed to the molecule class.
The name is needed for saving the results of the calculations.
In order to start the calculation of the eigenfunctions and eigenvalues the solve() function of the H2O molecule object is used:
H2O.solve()
This process takes several minutes. After a matrix element is calculated successfully
the element (Sij ord Hij) and the time it took to calculate it is displayed. The total time needed is shown after the function is done.
Now the results can be saved using:
H2O.save_molecule_data()
This function creates two .txt files with the following layout:
"name_construction.txt"
atom_1 (name of atom 1)
$x_{1:0}$ $y_{1:0}$ $z_{1:0}$ (position of atom 1)
$Z_1$ (atomic number of atom 1)
$n_1$ $l_1$ $m_1$ $n_2$ $l_2$ $m_2$ ... $n_{k_1}$ $l_{k_1}$ $m_{k_1}$ (quantum numbers of the $k_1$ orbitals of atom 1)
atom_2
$x_{2:0}$ $y_{2:0}$ $z_{2:0}$
$Z_1$
$n_1$ $l_1$ $m_1$ $n_2$ $l_2$ $m_2$ ... $n_{k_2}$ $l_{k_2}$ $m_{k_2}$
. .
. .
. .
atom_j
$x_{j:0}$ $y_{j:0}$ $z_{j:0}$
$Z_j$
$n_1$ $l_1$ $m_1$ $n_2$ $l_2$ $m_2$ ... $n_{k_j}$ $l_{k_j}$ $m_{k_j}$
"name_solution.txt"
$E_1$ (lowest energy)
$v_1^1$ $v_1^2$ ... $v_1^k$ (eigenvector with lowest energy, $k$ is the total number of orbitals
$E_1$
$v_2^1$ $v_2^2$ ... $v_2^k$
. .
. .
. .
$E_k$ (highest energy)
$v_k^1$ $v_k^2$ ... $v_k^k$ (eigenvector with highest energy)
In order to plot the result of the calculation open either orbitalplot2d.py or orbitalplot3d.py .
In orbitalplot2d.py and orbitalplot3d.py the variable name has to be change to the name of the molecule defined in the first step.
If orbitalplot3d.py is used the index of the molecular orbital ef_idx has to be specified as well.
Now the scripts can be run to plot the molecular orbitals.
The calculation of the matrix elements for the hamiltonian and the overlap matrix takes very long. The results of these integrals also strongly depend on the upper integration limit of r. For $H_2 O$ the best results were obtained using $r_{max} = 50 a_0$. For lower values of $r_{max}$ some of the resulting molecular orbitals were asymmetric which is not possible given the symmetry of the $H_2 O$ molecule.