Key breakthrough: Schrödinger's equation
\begin{equation} \hat{H} \Psi = E \Psi. \end{equation}Examples of systems that can be solved exactly:
All are quantized.
The Schrödinger equation splits into two parts: the kinetic energy piece $\hat{T}$ and the potential energy piece $\hat{V}$, giving
\begin{equation} \left( \hat{T} + \hat{V} \right) \Psi = E \Psi. \end{equation}All are operators, taking one function (the wavefunction here) to another function.
The wavefunction is a function of the $3N$ coordinates, where there are $N$ particles. This will be written as
\begin{equation} \Psi(x_1, y_1, z_1, x_2, y_2, z_2, \dots, x_N, y_N, z_N) = \Psi({\bf r}_1, {\bf r}_2, \dots, {\bf r}_N). \end{equation}For a hydrogen molecule (two hydrogen nuclei connected by a bond), we consider the nuclei as one particle each, with a single coordinate. For the purposes of counting for now, we note that we have two H nuclei, plus two electrons, for four particles, hence 12 coordinates.
For a benzene molecule (six carbon atoms connected by various bonds, each connected to a hydrogen atom) we have six carbon nuclei and six hydrogen nuclei, and each carbon atom has six electrons (for 36 electrons), and each hydrogen atom has six electrons. This gives 54 particles for 162 coordinates.
Observables or expectation values correspond to experimental measurements. According to quantum theory we compute the outcome of any measurement by "averaging" the appropriate operator as
\begin{equation} < x > = \frac{\int_{-\infty}^{\infty} \Psi^*(x) \hat{x} \Psi(x) \, \text{d}x}{\int_{-\infty}^{\infty} \Psi^*(x) \Psi(x) \, \text{d}x}. \end{equation}This particular example gives the expectation value, or average value, of the position along the $x$-axis, where the operator for position $\hat{x}$ is the multiplicative operator $\hat{x} \Psi(x) = x \Psi(x)$.
Write
\begin{equation} \int \Psi^*({\bf x}) \, \hat{C} \, \Phi({\bf x}) \, \text{d}{\bf x} \longrightarrow \langle \Psi | \hat{C} | \Phi \rangle. \end{equation}The "bra" is
\begin{equation} \langle \Psi | \equiv \int \text{d}{\bf x} \, \Psi^*({\bf x}). \end{equation}The "ket" is
\begin{equation} | \Phi \rangle \equiv \Phi({\bf x}). \end{equation}The operator doesn't appear explicitly; it's the $\hat{C}$ in "bra"$\hat{C}$"ket".
As an example, the Schrödinger equation is written $\hat{H} | \Psi \rangle = E | \Psi \rangle$.
The key result of the examples is that
\begin{equation} E = \frac{\langle \Psi | \hat{H} | \Psi \rangle}{\langle \Psi | \Psi \rangle}. \end{equation}Most operators are obvious: take the classical expression directly (eg potential $V(x)$ has operator $V(x)$). However, momentum is different:
\begin{equation} m v_x \longrightarrow - i \hbar \frac{\text{d}}{\text{d}x} \end{equation}Now, the Hamiltonian is the sum of the kinetic and potential energy operators. Kinetic energy is
\begin{align} T &= \frac{1}{2} m v_x^2 \\ &= \frac{(m v_x)^2}{2 m} \\ &= \frac{1}{2m} p_x^2. \end{align}Hence the previous result gives
\begin{align} \hat{T} &= \frac{1}{2m} \hat{p}_x^2 \\ &= \frac{1}{2m} \hat{p}_x \hat{p}_x \\ &= \frac{1}{2m} \left( -i \hbar \frac{\text{d}}{\text{d}x} \right) \left( -i \hbar \frac{\text{d}}{\text{d}x} \right) \\ & = -\frac{\hbar^2}{2m} \frac{\text{d}^2}{\text{d}x^2}. \end{align}Hence the Hamiltonian operator is
\begin{equation} \hat{H} = -\frac{\hbar^2}{2m} \frac{\text{d}^2}{\text{d}x^2} + V(x). \end{equation}Simplest case is a hydrogen atom. The location of the nucleus is given by ${\bf R}$. The location of the electron is given by ${\bf r}$. The vector connecting the nucleus to the electron is therefore ${\bf r} - {\bf R}$.
We remember that the force between two charged particles separated by distance $r$ is given by Coulomb's Law $F = q_1 q_2 / r^2$ and has potential $v({\bf r}) = 1 / (4 \pi \epsilon_0) q_1 q_2 / |{\bf r}|$.
The Hamiltonian operator for this single-atom system is
\begin{equation} \hat{H} = -\frac{\hbar^2}{2M} \nabla^2_{\bf R} -\frac{\hbar^2}{2m} \nabla^2_{\bf r} - \frac{1}{4 \pi \epsilon_0} \frac{e^2}{|{\bf r} - {\bf R}|}, \end{equation}where $m$ is the mass of the electron and $M$ the mass of the nucleus.
Choose units so that $\hbar^2 / m_e = 1 = e^2 / (4 \pi \epsilon_0)$. This simplifies the Hamiltonian to
\begin{equation} \hat{H} = -\frac{1}{2M} \nabla^2_{\bf R} -\frac{1}{2} \nabla^2_{\bf r} - \frac{1}{|{\bf r} - {\bf R}|}, \end{equation}where $M$ is now written in terms of electron masses, and energies are now in Hartrees.
The final term is the repulsion between the two electrons. The factors of two follow as the nucleus now contains two protons, so has charge 2.
Notation is such that ${\bf R}_1 = {\bf R}_O$, ${\bf R}_2 = {\bf R}_{H_1}$, and so on. Note how the double sums only count each term once.
Nuclei are much heavier than electrons ($m_p \approx 1800 m_e$), so it's a good approximation to assume that the electrons move much faster than nuclei: the electrons instantaneously equilibriate to the locations of the nuclei.
This is the core of the Born-Oppenheimer approximation: Solve Schrödinger's equation only in the electronic coordinates, with the nuclear coordinates held fixed. We then split the full Hamiltonian into the nuclear piece and the electronic piece $\hat{H}_{\text{elec}}$, and only solve for
\begin{align} \hat{H}_{\text{elec}} \Phi_{\text{elec}} &= E_{\text{elec}} \Phi_{\text{elec}} \\ \Phi_{\text{elec}} &= \Phi_{\text{elec}} \left( \left\{ {\bf r}_i \right\}; \left\{ {\bf R}_I \right\} \right) \\ E_{\text{elec}} &= E_{\text{elec}} \left( \left\{ {\bf R}_I \right\} \right) \end{align}That is, the solutions that we get depend parametrically on the coordinates of the nuclei.
For the water molecule, the Born-Oppenheimer approximation gives
\begin{equation} \hat{H}_{\text{elec}} = - \sum_{i=1}^{10} \frac{1}{2} \nabla^2_{{\bf r}_i} - \sum_{i=1}^{10} \frac{8}{|{\bf r}_i - {\bf R}_O|} - \sum_{i=1}^{10} \frac{1}{|{\bf r}_i - {\bf R}_{H_1}|} - \sum_{i=1}^{10} \frac{1}{|{\bf r}_i - {\bf R}_{H_2}|} + \sum_{i=1}^{10} \sum_{j=i+1}^{10} \frac{1}{|{\bf r}_i - {\bf r}_j|} + \sum_{I=1}^{3} \sum_{J=I+1}^{3} \frac{Z_I Z_J}{|{\bf R}_I - {\bf R}_J|}. \end{equation}That is, terms given by varying the location of the nuclei vanish. In addition, the final nuclear interaction term becomes constant.
Now the approximations get more serious.
A wavefunction for a single electron is called a molecular orbital (MO). MOs with spatial and spin coordinates are called spin orbitals and are products of a spatial orbital and a spin function. There are only two kinds of spin function: "up" ($\alpha$) and "down" ($\beta$). We write
\begin{align} \chi^{\uparrow}({\bf x}) &= \psi({\bf r}) \alpha(\omega), \\ \chi^{\downarrow}({\bf x}) &= \psi({\bf r}) \beta(\omega) \end{align}In order to ensure that our (approximate) full wavefunction has the right symmetries, but is made out of single electron MOs, we use the Slater determinant, eg
\begin{equation} \Psi({\bf x}_1, {\bf x}_2) = \begin{vmatrix} \chi_1({\bf x}_1) & \chi_2({\bf x}_1) \\ \chi_1({\bf x}_2) & \chi_2({\bf x}_2) \end{vmatrix} \end{equation}This can be extended to an arbitrary number of MOs by increasing the number of terms in the determinant. However, it should include a $N!^{-1/2}$ normalization so that $\langle \Psi_{\text{SD}} | \Psi_{\text{SD}} \rangle = 1$.
You will often see the representation
\begin{equation} \Psi_{\text{SD}}({\bf x}_1, {\bf x}_2, \dots {\bf x}_N) = | \chi_1({\bf x}_1) \chi_2({\bf x}_2) \dots \chi_N({\bf x}_N) \rangle \end{equation}The spin orbitals are always orthonormal $\langle \chi_i | \chi_j \rangle = \delta_{ij}$.
A determinant is zero if two or more of its rows or columns are identical: therefore a spin orbital can be included only once in a Slater determinant (otherwise the wavefunction is zero everywhere, which isn't physically acceptable). This is the Pauli exclusion principle successfully encoded.
The electronic energy of a single SD is
\begin{align} E_{\text{SD}} &= \frac{ \langle \Psi_{\text{SD}} | \hat{H}_{\text{elec}} | \Psi_{\text{SD}} \rangle }{ \langle \Psi_{\text{SD}} | \Psi_{\text{SD}} \rangle } \\ & = \langle \Psi_{\text{SD}} | \hat{H}_{\text{elec}} | \Psi_{\text{SD}} \rangle \end{align}with
\begin{equation} \hat{H}_{\text{elec}} = - \sum_{i=1}^{N_{\text{elec}}} \frac{1}{2} \nabla_i^2 - \sum_{I=1}^{N_{\text{nuclei}}} \sum_{i=1}^{N_{\text{elec}}} \frac{Z_I}{|{\bf r}_i - {\bf R}_I|} + \sum_{i=1}^{N_{\text{elec}}} \sum_{j=i+1}^{N_{\text{elec}}} \frac{1}{|{\bf r}_i - {\bf r}_j|}. \end{equation}After manipulations we find
\begin{equation} E_{\text{kin}} = \sum_i \langle \chi_i | - \nabla_i^2 | \chi_i \rangle \end{equation}and the nuclear energies give
\begin{equation} E_{e-n} = \sum_i \langle \chi_i | \sum_I \frac{-Z_I}{|{\bf r}_i - {\bf R}_I|} | \chi_i \rangle \end{equation}The nuclear attraction potential for each electron is
\begin{equation} \hat{v}_{\text{ext}} = \end{equation}The difficult term is the electron-electron repulsion energy, or the two-electron integrals
\begin{equation} E_{e-e} = \frac{1}{2} \sum_i \sum_j \left[ \int\int \dots - \int\int \dots \right] \end{equation}where $r_{12} = |{\bf r}_1 - {\bf r}_2|$.
Remember that here $1, 2$ are dummy indices for the integration, not the first and second electrons.
The computational cost of computing these integrals goes as the fourth power of the number of electrons.
The part of the energy that doesn't include the two integral terms is written $\sum_i \langle \chi_i | \hat{h} | \chi_i \rangle$, where $\hat{h} = - \frac{1}{2} \nabla^2 + \sum_I \frac{-Z_I}{|{\bf r}_i - {\bf R}_I|}$ is the core Hamiltonian, which neglects the electron-electron coupling.
The spin terms can be integrated out, so that the sums over each term can become $2 \sum_{i=1}^{N_{\text{elec}/2}}$ (integer division). The energy of the Slater determinant then breaks down into the core Hamiltonian $h_i$, the Coulomb integral $J_{ij}$ which is symmetric in the coordinates ${\bf r}_{1,2}$, and the Exchange integral $K_{ij}$, which is not (and doesn't have a physical interpretation).
In summary,
\begin{equation} E_{\text{SD}} = 2 \sum_{i=1}^{N_{\text{elec}}/2} h_i + \sum_{i=1}^{N_{\text{elec}}/2} \sum_{j=1}^{N_{\text{elec}}/2} ( 2 J_{ij} - K_{ij}). \end{equation}Assume we have two Slater determinants for a molecule, $\Psi^{A,B}_{\text{SD}}$, with energies $E^{A,B}_{\text{SD}} = \langle \Psi^{A,B}_{\text{SD}} | \hat{H} | \Psi^{A,B}_{\text{SD}} \rangle$. Then if $E^A_{\text{SD}} > E^B_{\text{SD}}$ then, according to teh variation principle, we take $\Psi^B$ to be a better approximation to the wavefunction.
By performing functional variation we can show that the minimal value occurs when
\begin{equation} \hat{f} \psi_i({\bf r}) = \epsilon_i \psi_i({\bf r}) \end{equation}This is an eigenfunction equation involving the Fock operator $\hat{f}$ which is made of three pieces,
\begin{equation} \hat{f} = \hat{h} + \sum_i ( 2 \hat{J}_i - \hat{K}_i ) \end{equation}where $h, J, K$ correspond to the core Hamiltonian, Coulomb and Exchange operators.
The eigenfunctions of the Fock operator are the optimal MOs, and the eigenvalues their associated energies.
An example: the Hartree-Fock eigenvalue equation for a molecule with 6 electrons is
\begin{equation} \left[ \hat{h} + \sum_{i=1}^3 \left( 2 \hat{J}_i - \hat{K}_i \right) \right] \Psi_j({\bf r}) = \epsilon_j \Psi_j({\bf r}) \end{equation}The index $j$ can be anything. The Fock operator is made of the first three orbitals $\Psi_{1,2,3}$.
So we solve iteratively using the Self-Consistent Field procedure. First guess the MOs $\psi_i$. Construct the Fock operator and solve the eigenfunction problem to find a new set of MOs. If the new MOs are close enough to the old MOs the procedure has converged; otherwise, construct the Fock operator and repeat.