The molecular Hamiltonian operator (in the molecular frame) consists of the kinetic energy operators of nuclei and electrons, the Coulomb electron-electron and nuclear-nuclear repulsion, and electron-nuclear attraction operators,

$$ H=T_n+T_e+V_\text{ee} + V_{nn} + V_{en} = \\ = \frac{1}{2}\sum_a m_a\dot{\mathbf{r}}_a^2+ \frac{1}{2}m_e\sum_i \dot{\mathbf{r}}_i^2

• \frac{e^2}{4\pi\epsilon_0}\sum_{i\neq j}\frac{1}{r_{ij}}
• \frac{e^2}{4\pi\epsilon_0}\sum_{a\neq b}\frac{Z_aZ_b}{r_{ab}}

• \frac{e^2}{4\pi\epsilon_0}\sum_{i,a}\frac{Z_a}{r_{ai}},

$$

where indices $a$ and $i$ refer to nuclei and electrons, respectively.

The quantum mechanical expression is obtained by replacing the classical Cartesian momentum operators $m\dot{\mathbf r}$ by their quantum-mechanical equivalents $-i\hbar\partial/\partial\mathbf{r}=-i\hbar\nabla$. Using atomic units, the expression for the Hamiltonian reads:

$$ H= -\frac{1}{2}\sum_a \frac{\nabla_a^2}{m_a} -\frac{1}{2}\sum_i \nabla_i^2

• \sum_{i\neq j}\frac{1}{r_{ij}}
• \sum_{a\neq b}\frac{Z_aZ_b}{r_{ab}}

• \sum_{i,a}\frac{Z_a}{r_{ai}}. $$

By taking the kinetic energy of nuclei out of the equation, we get the electronic Hamiltonian $H_e$, for which the energies $E_n(r_a)$ and wave functions $\phi_n(r_i;r_a)$ may be calculated for several electronic states $n=1,2,3...$ using the electronic structure techniques developed by quantum chemistry,

$$ H_e= -\frac{1}{2}\sum_i \nabla_i^2

• \sum_{i\neq j}\frac{1}{r_{ij}}
• \sum_{a\neq b}\frac{Z_aZ_b}{r_{ab}}

• \sum_{i,a}\frac{Z_a}{r_{ai}}, \\ H_e\phi_n(r_i;r_a)=E_n(r_a)\phi_n(r_i;r_a). $$

The electronic wave functions are functions of electronic coordinates, but they also have a parametric dependency on nuclei coordinates due to the Coulomb potential terms' dependence on nuclear positions $r_{ia}$.

Now, being able to calculate the electronic solutions $\phi_n$, we seek the solutions for the total Hamiltonian $H=-\frac{1}{2}\sum_a \frac{\nabla_a^2}{m_a} + H_e$ in the following form (also called as Born-Huang ansatz):

$$ \Psi_p(r_i,r_a) = \sum_n c_n^{(p)}(r_a)\phi_n(r_i;r_a), $$

i.e., we expand the total wave function into a linear combination of different electronic states with linear coefficients depending on nuclear coordinates. We plugin this representation of wave function into the Schrödinger equation for the total Hamiltonian. First, let's have a look at the effect of the nuclear kinetic energy operator when applied to the wave function, i.e.

$$ -\frac{1}{2}\sum_a \frac{\nabla_a^2}{m_a}\Psi_p(r_i,r_a) = -\frac{1}{2}\sum_a\frac{1}{m_a}\sum_n\nabla_a\cdot\nabla_a \left(c_n^{(p)}(r_a)\phi_n(r_i;r_a)\right) = \\ -\frac{1}{2}\sum_a\frac{1}{m_a}\sum_n\nabla_a\cdot\left[\phi_n(r_i;r_a)\nabla_a c_n^{(p)}(r_a) + c_n^{(p)}(r_a)\nabla_a\phi_n(r_i;r_a)\right] = \\ -\frac{1}{2}\sum_a\frac{1}{m_a}\sum_n\left[\phi_n(r_i;r_a)\nabla_a^2 c_n^{(p)}(r_a) + c_n^{(p)}(r_a)\nabla_a^2\phi_n(r_i;r_a)+ 2\nabla_a\phi_n(r_i;r_a)\cdot\nabla_a c_n^{(p)}(r_a) \right]

$$

Combining it with the electronic Hamiltonian, the Schrödinger equation can be written as

$$ \left(-\frac{1}{2}\sum_a \frac{\nabla_a^2}{m_a} +H_e\right)\Psi_p(r_i,r_a) = E_p\Psi_p(r_i,r_a) \\ \sum_n \left[ c_n^{(p)}(r_a) E_n(r_a)\phi_n(r_i;r_a) -\frac{1}{2}\sum_a\frac{1}{m_a}\phi_n(r_i;r_a)\nabla_a^2 c_n^{(p)}(r_a) \right.- \\ \left. -\frac{1}{2}\sum_a\frac{1}{m_a}c_n^{(p)}(r_a)\nabla_a^2\phi_n(r_i;r_a) -\sum_a\frac{1}{m_a}\nabla_a\phi_n(r_i;r_a)\cdot\nabla_a c_n^{(p)}(r_a) \right] = \\ =E_p\sum_n c_n^{(p)}(r_a)\phi_n(r_i;r_a).

$$

By taking the inner product of the left and right hand sides with $\phi_l(r_i;r_a)$ the equation becomes

$$ c_l^{(p)}(r_a) E_l(r_a) -\frac{1}{2}\sum_a\frac{1}{m_a}\nabla_a^2 c_l^{(p)}(r_a)- \\ -\frac{1}{2}\sum_n\sum_a\frac{1}{m_a}c_n^{(p)}(r_a)\langle \phi_l|\nabla_a^2|\phi_n\rangle -\sum_n\sum_a\frac{1}{m_a}\langle \phi_l|\nabla_a|\phi_n\rangle\cdot\nabla_a c_n^{(p)}(r_a) = \\ =E_pc_l^{(p)}(r_a)

$$

The Born-Oppenheimer approximation may now be presented; it consists in ignoring terms containing derivatives of electronic wave functions with respect to nuclear coordinates. These terms are grouped together on the second line of the above equation and are referred to as non-adiabatic derivative coupling terms. In the Born-Oppenheimer approximation, the Schrödinger equation for coupled electronic states reduces to a set of independent equations describing motions of nuclei in each electronic state,

$$ c_l^{(p)}(r_a) E_l(r_a) -\frac{1}{2}\sum_a\frac{1}{m_a}\nabla_a^2c_l^{(p)}(r_a) =E_pc_l^{(p)}(r_a). $$

It consists of two operators: electronic energy of the $l$-th electronic state, which is a function of nuclear coordinates, and nuclei kinetic energy. The electronic energy is referred to as the molecular Potential Energy Surface, and its correct computation is the primary job of quantum chemistry theories.