The Born-Oppenheimer Approximation

C. David Sherrill
School of Chemistry and Biochemistry
Georgia Institute of Technology
March 2005

We may write the nonrelativistic Hamiltonian for a molecule as a sum of five terms:

$${\hat H} = - \frac{\hbar^2}{2m} \sum_{i} \nabla^2_i - \s... ..._0 R_{AB}} + \sum_{i>j} \frac{e^2}{4 \pi \epsilon_0 r_{ij}},$$ (1)

where $i, j$ refer to electrons and $A, B$ refer to nuclei. In atomic units, this is just

$${\hat H} = - \frac{1}{2} \sum_{i} \nabla^2_i - \sum_{A} ... ..._{A>B} \frac{Z_A Z_B}{R_{AB}} + \sum_{i>j} \frac{1}{r_{ij}}.$$ (2)

The Schrödinger equation may be written more compactly as

$${\hat H} = {\hat T}_N(\mathbf{R}) + {\hat T}_e(\mathbf{r}) ... ...R}) + {\hat V}_{NN}(\mathbf{R}) + {\hat V}_{ee}(\mathbf{r}),$$ (3)

where ${\mathbf R}$ is the set of nuclear coordinates and ${\mathbf r}$ is the set of electronic coordinates. If spin-orbit effects are important, they can be added through a spin-orbit operator ${\hat H}_{so}$.

Unfortunately, the ${\hat V}_{eN}(\mathbf{r, R})$ term prevents us from separating ${\hat H}$ into nuclear and electronic parts, which would allow us to write the molecular wavefunction as a product of nuclear and electronic terms, $\Psi(\mathbf{r,R}) = \Psi(\mathbf{r}) \chi(\mathbf{R})$. We thus introduce the Born-Oppenheimer approximation, by which we conclude that this nuclear and electronic separation is approximately correct. The term ${\hat V}_{eN}(\mathbf{r, R})$ is large and cannot be neglected; however, we can make the ${\mathbf R}$ dependence parametric, so that the total wavefunction is given as $\Psi(\mathbf{r;R}) \chi(\mathbf{R})$. The Born-Oppenheimer approximation rests on the fact that the nuclei are much more massive than the electrons, which allows us to say that the nuclei are nearly fixed with respect to electron motion. We can fix $\mathbf{R}$, the nuclear configuration, at some value $\mathbf{R_a}$, and solve for the electronic wavefunction $\Psi(\mathbf{r; R_a})$, which depends only parametrically on $R$. If we do this for a range of $\mathbf{R}$, we obtain the potential energy curve along which the nuclei move.

We now show the mathematical details. Initially, ${\hat T}_N(\mathbf{R})$ can be neglected since ${\hat T}_N$ is smaller than ${\hat T}_e$ by a factor of $M_A / \mu_e$, where $\mu_e$ is the reduced mass of an electron. Thus for a fixed nuclear configuration, we have

$${\hat H}_{el} = {\hat T}_e(\mathbf{r}) + {\hat V}_{eN}(\... ...}) + {\hat V}_{NN}(\mathbf{R}) + {\hat V}_{ee}(\mathbf{r})$$ (4)

such that

$${\hat H}_{el} \Psi(\mathbf{r; R}) = E_{el} \Psi(\mathbf{r; R})$$ (5)

This is the ``clamped-nuclei'' Schrödinger equation. Quite frequently ${\hat V}_{NN}(\mathbf{R})$ is neglected in the above equation, which is justified since in this case $\mathbf{R}$ is just a parameter so that ${\hat V}_{NN}(\mathbf{R})$ is just a constant and shifts the eigenvalues only by some constant amount. Leaving ${\hat V}_{NN}(\mathbf{R})$ out of the electronic Schrödinger equation leads to a similar equation,

$${\hat H}_e = {\hat T}_e(\mathbf{r}) + {\hat V}_{eN}(\mathbf{r;R}) + {\hat V}_{ee}(\mathbf{r})$$ (6)

$${\hat H}_e \Psi(\mathbf{r;R}) = E_e \Psi_e(\mathbf{r;R})$$ (7)

For the purposes of these notes, we will assume that ${\hat V}_{NN}(\mathbf{R})$ is included in the electronic Hamiltonian. Additionally, if spin-orbit effects are important, then these can be included at each nuclear configuration according to

$${\hat H}_0 = {\hat H}_{el} + {\hat H}_{so}$$ (8)

$${\hat H}_0 \Psi(\mathbf{r;R}) = E_0 \Psi(\mathbf{r;R})$$ (9)

Consider again the original Hamiltonian (1). An exact solution can be obtained by using an (infinite) expansion of the form

$$\Psi(\mathbf{r, R}) = \sum_k \Psi_k(\mathbf{r;R}) \chi_k(\mathbf{R}),$$ (10)

although, to the extent that the Born-Oppenheimer approximation is valid, very accurate solutions can be obtained using only one or a few terms. Alternatively, the total wavefunction can be expanded in terms of the electronic wavefunctions and a set of pre-selected nuclear wavefunctions; this requires the introduction of expansion coefficients:

$$\Psi^i(\mathbf{r, R}) = \sum_{kl} c^i_{kl} \Psi_k(\mathbf{r;R}) \chi_{kl}(\mathbf{R})$$ (11)

where the superscript $i$ has been added as a reminder that there are multiple solutions to the Schrödinger equation.

Expressions for the nuclear wavefunctions $\chi_k(\mathbf{R})$ can be obtained by inserting the expansion (10) into the total Schrödinger equation yields

$${\hat H} \sum_{k'} \Psi_{k'}(\mathbf{r;R}) \chi_{k'}(\mathbf... ... E \sum_{k''} \Psi_{k''}(\mathbf{r;R}) \chi_{k''}(\mathbf{R})$$ (12)

or

$$\int d\mathbf{r} \Psi_k^*(\mathbf{r;R}) {\hat H} \sum_{k'}... ...(\mathbf{r;R}) \chi_{k'}(\mathbf{R}) = E \chi_{k}(\mathbf{R})$$ (13)

if the electronic functions are orthonormal. Simplifying further,

$$E \chi_{k}(\mathbf{R}) = \sum_{k'} \int d\mathbf{r} \Psi_{k... ...H}_{so} \right] \Psi_{k'}(\mathbf{r;R}) \chi_{k'}(\mathbf{R})$$ (14)

$$= \sum_{k'} \left\{ \langle \Psi_{k}(\mathbf{r;R}) \vert {\hat H}_{... ...t H}_{so} \vert \Psi_{k'}(\mathbf{r;R}) \rangle \right\} \chi_{k'}({\mathbf R})$$

$$+ \sum_{k'} \int d\mathbf{r} \Psi_{k}^*(\mathbf{r;R}) {\hat T}_N \Psi_{k'}(\mathbf{r;R}) \chi_{k'}(\mathbf{R})$$

The last term can be expanded using the chain rule to yield

$$\sum_A \frac{-1}{2M_A} \left[ \nabla^2_A \chi_{k}(\mathbf... ... \vert \nabla^2_A \vert \Psi_k(\mathbf{r;R}) \rangle \right]$$ (15)

$$- \sum_A \frac{1}{2M_A} \sum_{k' \neq k} \left[ 2 \langle \Psi_k(\m... ...k(\mathbf{r;R}) \vert \nabla^2_A \vert \Psi_{k'}(\mathbf{r;R}) \rangle \right].$$

At this point, a more compact notation is very helpful. Following Tully [1], we introduce the following quantities:

$$U_{kk'}(\mathbf{R}) = U^{(el)}_{kk'}(\mathbf{R}) + U^{(so)}_{kk'}(\mathbf{R})$$ (16)

$$U^{(el)}_{kk'}(\mathbf{R}) = \langle \Psi_k(\mathbf{r;R}) \vert {\hat H}_{el} \vert \Psi_{k'}(\mathbf{r;R}) \rangle$$ (17)

$$U^{(so)}_{kk'}(\mathbf{R}) = \langle \Psi_k(\mathbf{r;R}) \vert {\hat H}_{so} \vert \Psi_{k'}(\mathbf{r;R}) \rangle$$ (18)

$${\hat T}_{kk'}'(\mathbf{R}) = \sum_A \frac{-1}{M_A} \mathbf{d}_{kk'}^{(A)}(\mathbf{R}) \cdot \nabla _A$$ (19)

$$T''_{kk'}(\mathbf{R}) = \sum_A \frac{-1}{2M_A} D_{kk'}^{(A)}(\mathbf{R})$$ (20)

$${\mathbf{d}}_{kk'}^{(A)}(\mathbf{R}) = \langle \Psi_k({\math... ...R}}) \vert \nabla _A \vert \Psi_{k'}({\mathbf{r;R}}) \rangle$$ (21)

$$D_{kk'}^{(A)}(\mathbf{R}) = \langle \Psi_k({\mathbf{r;R}}) \vert \nabla^2_A \vert \Psi_{k'}({\mathbf{r;R}}) \rangle$$ (22)

Note that equation (18) of reference [1] should not contain a factor of 1/2 as it does. Now we can rewrite equations (14) and (15) as

$$\left[ {\hat T}_N + \sum_A \left( \frac{-1}{2M_A} \right) ... ...k}^{(A)} \right\} + U_{kk} - E \right] \chi_k({\mathbf{R}})$$ (23)

$$= - \sum_{k \neq k'} \left[ U_{kk'} + \sum_{A} \left( \frac{-1}{2M_... ...}_{kk'}^{(A)} \nabla _A + D_{kk'}^{(A)} \right\} \right] \chi_k'({\mathbf{R}}),$$

or

$$\left[ {\hat T}_N + {\hat T}'_{kk} + T''_{kk} + U_{kk} - E ... ...+ {\hat T}'_{kk'} + T''_{kk'} \right] \chi_{k'}({\mathbf{R}})$$ (24)

This is equation (14) of Tully's article [1]. Tully simplifies this equation by one more step, removing the term ${\hat T}'_{kk}$. By taking the derivative of the overlap of $\langle \Psi_k({\mathbf{r;R}}) \vert \Psi_k({\mathbf{r;R}}) \rangle$ it is easy to show that this term must be zero when the electronic wavefunction can be made real. If we use electronic wavefunctions which diagonalize the electronic Hamiltonian, then the electronic basis is called adiabatic, and the coupling terms $U_{kk'}$ vanish.¹ This is the general procedure. However, the equation above is formally exact even if other electronic functions are used. In some contexts it is preferable to minimize other coupling terms, such as ${\hat T}'_{kk'}$; this results in a diabatic electronic basis. Note that the first-derivative nonadiabatic coupling matrix elements ${\hat T}'_{kk'}$ are usually considered more important than the second-derivative ones, $T''_{kk'}$.

In most cases, the couplings on the right-hand side of the preceeding equation are small. If they can be safely neglected, and assuming that the wavefunction is real, we obtain the following equation for the motion of the nuclei on a given Born-Oppenheimer potential energy surface:

$$\left[{\hat T}_N + T''_{kk} + U_{kk} \right] \chi_k({\mathbf{R}}) = E \chi_k({\mathbf{R}})$$ (25)

This equation clearly shows that, when the off-diagonal couplings can be ignored, the nuclei move in a potential field set up by the electrons. The potential energy at each point is given primarily by $U_{kk}$ (the expectation value of the electronic energy), with a small correction factor $T''_{kk}$. Following Steinfeld [2], we can estimate the magnitude of the term $T''_{kk}$ as follows: a typical contribution has the form $1/(2M_A) \nabla^2_A \Psi({\mathbf{r;R}})$, but $\nabla _A \Psi({\mathbf{r;R}})$ is of the same order as $\nabla _i \Psi({\mathbf{r;R}})$ since the derivatives operate over approximately the same dimensions. The latter is $\Psi({\mathbf{r;R}})p_e$, with $p_e$ the momentum of an electron. Therefore $1/(2M_A) \nabla^2_A \Psi({\mathbf{r;R}}) \approx p_e^2/(2M_A) = (m/M_A)E_e$. Since $m/M_A \sim 1/10000$, this term is expected to be small, and it is usually dropped. However, to achieve very high accuracy, such as in spectroscopic applications, this term must be retained.

The Born-Oppenheimer Diagonal Correction In a perturbation theory analysis of the Born-Oppenheimer approximation, the first-order correction to the Born-Oppenheimer electronic energy due to the nuclear motion is the Born-Oppenheimer diagonal correction (BODC),

$$E_{BODC} = \langle \Psi({\mathbf r}; {\mathbf R}) \vert {\hat T}_n \vert \Psi({\mathbf r}; {\mathbf R}) \rangle ,$$ (26)

which can be applied to any electronic state $\Psi({\mathbf r}; {\mathbf R})$. The BODC is also referred to as the ``adiabatic correction.'' One of the first systematic investigations of this effect was a study by Handy, Yamaguchi, and Schaefer in 1986 [3]. In this work, the authors evaluated the BODC using Hartree-Fock self-consistent-field methods (and, where relevant, two-configuration self-consistent-field) for a series of small molecules. One interesting finding was that the BODC changes the singlet-triplet splitting in methylene by 40 cm$^{-1}$, which is small on a ``chemical'' energy scale but very relevant for a spectroscopic energy scale. Inclusion of the BODC is required to accurately model the very dense rovibrational spectrum of the water molecule observed at high energies, and these models were a critical component in proving the existence of water on the sun [4,5].

For many years, it was only possible to compute the BODC for Hartree-Fock or multiconfigurational self-consistent-field wavefunctions. However, in 2003 the evaluation of the BODC using general configuration interaction wavefunctions was implemented [6] and its convergence toward the ab initio limit was investigated for H$_2$, BH, and H$_2$O. This study found that the absolute value of the BODC is difficult to converge, but errors in estimates of the BODC largely cancel each other so that even BODC's computed using Hartree-Fock theory capture most of the effect of the adiabatic correction on relative energies or geometries. Table 1 displays the effect of the BODC on the barrier to linearity in the water molecule and the convergence of this quantity with respect to basis set and level of electron correlation. Although the absolute values of the BODC's of bent and linear water change significantly with respect to basis set and level of electron correlation, their difference does not change much as long as a basis of at least cc-pVDZ quality is used. For the cc-pVDZ basis, electron correlation changes the differential BODC correction by about 1 cm$^{-1}$. Table 2 displays the effect of the BODC on the equilibrium bond lengths and harmonic vibrational frequencies of the BH, CH$^+$, and NH molecules [7] and demonstrates somewhat larger changes to the spectroscopic constants than one might have expected, particularly for BH.

Table 1: Adiabatic correction to the barrier to linearity of water in the ground state (in cm$^{-1}$)$^{\rm a}$

Basis	Method	$C_{2v}$	$D_{\infty h}$	$\Delta E_e^{\rm b}$
DZ	RHF	613.66	587.69	-25.97
DZ	CISD	622.40	596.43	-25.97
DZ	CISDT	623.62	597.56	-26.06
DZ	CISDTQ	624.56	598.28	-26.28
DZ	CISDTQP	624.61	598.32	-26.29
DZP	RHF	597.88	581.32	-16.56
cc-pVDZ	RHF	600.28	585.20	-15.08
cc-pVDZ	CISD	615.03	599.15	-15.88
cc-pVDZ	CISDT	616.82	600.62	-16.20
cc-pVTZ	RHF	596.53	581.43	-15.10
cc-pVTZ	CISD	611.89	596.73	-15.16
cc-pVQZ	RHF	595.57	580.72	-14.85
$^{\rm a}$Data from Valeev and Sherrill [6].
$^{\rm b}$The difference between the adiabatic correction
for the $C_{2v}$ and $D_{\infty h}$ structures.

Table 2: Adiabatic corrections to bond length and harmonic frequencies of BH, CH$^{+}$, and NH$^{\rm a}$

	BH	CH$^+$	NH
$\Delta r_e$	0.00066	0.00063	0.00027
$\Delta \omega_e$	-2.25	-2.81	-1.38
$^{\rm a}$Data from Temelso, Valeev, and
Sherrill [7].