A water molecule consists of one oxygen atom and two hydrogen atoms. The O-H bond length is 0.097 nm and the H-O-H bond angle is 104.5°.

The molecular orbital Hamiltonian for water is,

\begin{equation} H_{\text{mo}}= - \frac{\hbar^2}{2m_e}\nabla^2 - \frac{8e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_{\text{O}}|} - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_{\text{H1}}|} - \frac{e^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_{\text{H2}}|} . \end{equation}

A trial wave function consisting of the 1s, 2s and 2p orbitals of the oxygen atom and the 1s orbitals of the hydrogen atoms is used.

$$\psi_{\text{mo}} = c_1\phi^{\text{O}}_{\text{1s}}+c_2\phi^{\text{O}}_{\text{2s}}+c_3\phi^{\text{O}}_{\text{2px}}+c_4\phi^{\text{O}}_{\text{2py}}+c_5\phi^{\text{O}}_{\text{2pz}}+c_6\phi^{\text{H1}}_{\text{1s}}+c_7\phi^{\text{H2}}_{\text{1s}}.$$

This trial wave function is substituted into the time independent Schrödinger equation,

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

Multiplying the Schrödinger equation from the left by each of the atomic orbitals and integrating over all space results in the Roothaan equations.

\[ \begin{equation} \langle \phi^{\text{O}}_{\text{1s}}(\vec{r}-\vec{r}_\text{O}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{O}}_{\text{1s}}(\vec{r}-\vec{r}_\text{O})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{O}}_{\text{2s}}(\vec{r}-\vec{r}_\text{O}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{O}}_{\text{2s}}(\vec{r}-\vec{r}_\text{O})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{O}}_{\text{2px}}(\vec{r}-\vec{r}_\text{O}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{O}}_{\text{2px}}(\vec{r}-\vec{r}_\text{O})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{O}}_{\text{2py}}(\vec{r}-\vec{r}_\text{O}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{O}}_{\text{2py}}(\vec{r}-\vec{r}_\text{O})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{O}}_{\text{2pz}}(\vec{r}-\vec{r}_\text{O}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{O}}_{\text{2pz}}(\vec{r}-\vec{r}_\text{O})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{H1}}_{\text{1s}}(\vec{r}-\vec{r}_\text{H1}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{H1}}_{\text{1s}}(\vec{r}-\vec{r}_\text{H1})|\psi_{\text{mo}} \rangle \\ \langle \phi^{\text{H2}}_{\text{1s}}(\vec{r}-\vec{r}_\text{H2}) | H_{\text{mo}} |\psi_{\text{mo}}\rangle =E\langle\phi^{\text{H2}}_{\text{1s}}(\vec{r}-\vec{r}_\text{H2})|\psi_{\text{mo}} \rangle \\ \end{equation} \]

The Roothaan equations can be written in matrix form,

\[ \begin{equation} \left[ \begin{matrix} H_{11} & H_{12} & H_{13} & H_{14} & H_{15} & H_{16} & H_{17}\\ H_{21} & H_{22} & H_{23} & H_{24} & H_{25} & H_{26} & H_{27}\\ H_{31} & H_{32} & H_{33} & H_{34} & H_{35} & H_{36} & H_{37}\\ H_{41} & H_{42} & H_{43} & H_{44} & H_{45} & H_{46} & H_{47}\\ H_{51} & H_{52} & H_{53} & H_{54} & H_{55} & H_{56} & H_{57}\\ H_{61} & H_{62} & H_{63} & H_{64} & H_{65} & H_{66} & H_{67}\\ H_{71} & H_{72} & H_{73} & H_{74} & H_{75} & H_{76} & H_{77}\\ \end{matrix} \right] \left[ \begin{matrix} c_1 \\ c_2 \\ c_3 \\ c_4 \\ c_5 \\ c_6 \\ c_7 \end{matrix} \right] = E \left[ \begin{matrix} S_{11} & S_{12} & S_{13} & S_{14} & S_{15} & S_{16} & S_{17}\\ S_{21} & S_{22} & S_{23} & S_{24} & S_{25} & S_{26} & S_{27}\\ S_{31} & S_{32} & S_{33} & S_{34} & S_{35} & S_{36} & S_{37}\\ S_{41} & S_{42} & S_{43} & S_{44} & S_{45} & S_{46} & S_{47}\\ S_{51} & S_{52} & S_{53} & S_{54} & S_{55} & S_{56} & S_{57}\\ S_{61} & S_{62} & S_{63} & S_{64} & S_{65} & S_{66} & S_{67}\\ S_{71} & S_{72} & S_{73} & S_{74} & S_{75} & S_{76} & S_{77}\\ \end{matrix} \right] \left[ \begin{matrix} c_1 \\ c_2 \\ c_3 \\ c_4 \\ c_5 \\ c_6 \\ c_7 \end{matrix} \right], \end{equation} \]

where the matrix elements are,

\[ \begin{equation} H_{ij}=\langle \phi_i|H_{\text{mo}}|\phi_j\rangle, \\ S_{ij}=\langle \phi_i|\phi_j\rangle. \end{equation} \]

The Hamiltonian matrix elements $H_{ij}$ and the overlap matrix elements $S_{ij}$ were determined numerically using Python and two versions in matlab: matlab1, matlab2 (including 3d orbitals). The calculations agree on the overlap matrix but produce different Hamiltonian matrices.

The Hamiltonian matrix and the overlap matrix calculated using the Python code are,

\[ \begin{equation} H= \left[ \begin{matrix} -832.396 & -4.928 & 0.0 & 0.0 & 0.0 & -68.355 & -68.355 \\ -4.928 & -99.915 & 0.0 & 0.0 & 0.0 & 79.004 & 79.004 \\ 0.0 & 0.0 & -101.112 & 0.841 & 0.0 & 40.116 & -40.116 \\ 0.0 & 0.0 & 0.841 & -101.112 & 0.0 & -40.116 & 40.116 \\ 0.0 & 0.0 & 0.0 & 0.0 & -97.618 & 0.0 & 0.0 \\ -49.921 & 56.966 & 20.687 & -20.687 & 0.0 & -129.73 & -57.09 \\ -49.922 & 56.966 & -21.195 & 21.195 & 0.0 & -57.09 & -129.73 \\ \end{matrix} \right], \qquad S = \left[ \begin{matrix} 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.06 & 0.06 \\ 0.0 & 1.0 & 0.0 & 0.0 & 0.0 & -0.542 & -0.542 \\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & -0.192 & 0.192 \\ 0.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.192 & -0.192 \\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0 \\ 0.06 & -0.542 & -0.192 & 0.192 & 0.0 & 1.0 & 0.371 \\ 0.06 & -0.542 & 0.192 & -0.192 & 0.0 & 0.371 & 1.0 \\ \end{matrix} \right]. \end{equation} \]

The molecular orbitals were then found by solving the eigen value problem,

$$\textbf{S}^{-1}\textbf{H}\,\vec{c}=E\,\vec{c}.$$

The molecular orbitals calculated by the Python code are:

$$E_1=-832.65\,\text{eV},\quad\left[ \begin{matrix} -0.9999 \\ -0.0121 \\ 0 \\ 0 \\ 0 \\ -0.0053 \\-0.0053 \end{matrix} \right] ;\qquad E_2=-123.42\,\text{eV}\quad\left[ \begin{matrix} 0.0927 \\ 0.5976 \\ 0 \\ 0 \\ 0 \\ -0.5631 \\ -0.5631 \end{matrix} \right] ;\qquad E_3=-110.7\,\text{eV}\quad\left[ \begin{matrix} -0.0029 \\ -0.0707 \\ -0.6868 \\ 0.6868 \\ 0 \\ 0.1793 \\ -0.1398 \end{matrix} \right] ;\qquad E_4=-100.27\,\text{eV}\quad\left[ \begin{matrix} 0 \\ 0 \\ 0.7071 \\ 0.7071 \\ 0 \\ 0 \\ 0 \end{matrix} \right] ;\qquad \\ \text{HOMO}=E_5=-97.62\,\text{eV}\quad\left[ \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ -1.0 \\ 0 \\ 0 \end{matrix} \right] ;\qquad \text{LUMO}=E_6=-91.34\,\text{eV}\quad\left[ \begin{matrix} 0.0303 \\ -0.9797 \\ 0 \\ 0 \\ 0 \\ -0.1400 \\ -0.1400 \end{matrix} \right] ;\qquad E_7=-76.0\,\text{eV}\quad\left[ \begin{matrix} -0.0026 \\ 0.0433 \\ -0.6295 \\ 0.6295 \\ 0 \\ -0.3068 \\ 0.3339 \end{matrix} \right] .$$