I'm an undergraduate Chemistry student and I'm currently writing my Final Dissertation.
I've written a form of the Born-Oppenheimer approximation with the only initial assumption of parametric dependence on nuclear coordinates. Let me know what do you think about it. Thanks in advance for any help you are able to provide.

For a given molecular system, within the non-relativistic approximation, the molecular Schrodinger equation can be written as $$\hat{\mathcal{H}}_{\mathrm{tot}}\Psi\left(\mathbf{R},\,\mathbf{r}\right)=\mathcal{E}_{\mathrm{tot}}\Psi\left(\mathbf{R},\,\mathbf{r}\right)\quad\left(1\right)$$ where $$\hat{\mathcal{H}}_{\mathrm{tot}}=\hat{\mathcal{T}}_{\mathrm{e}}\left(\mathbf{r}\right)+\hat{\mathcal{T}}_{\mathrm{N}}\left(\mathbf{R}\right)+\hat{\mathcal{V}}_{\mathrm{ee}}\left(\mathbf{r}\right)+\hat{\mathcal{V}}_{\mathrm{NN}}\left(\mathbf{R}\right)+\hat{\mathcal{V}}_{\mathrm{eN}}\left(\mathbf{R},\mathbf{r}\right)=\hat{\mathcal{T}}_{\mathrm{N}}+\hat{\mathcal{H}}_{\mathrm{e}}$$

1. Without loss of generality, if there exists a complete basis $\left\{ \Psi_{k}\left(\mathbf{R},\,\mathbf{r}\right)\right\}$ , the molecular wave function $\Psi\left(\mathbf{R},\,\mathbf{r}\right)$ can be expanded as $$\Psi\left(\mathbf{R},\,\mathbf{r}\right)={\displaystyle \sum_{k}}c_{k}\Psi_{k}\left(\mathbf{R},\,\mathbf{r}\right)$$

2. The nuclei in a molecule are heavy enough that, to a good approximation, they behave as classical particles. Thus, it is completely meaningful to talk about nuclear configurations, or in other words to assume a parametric dependence from nuclear coordinates. $$\hat{\mathcal{H}}_{\mathrm{tot}}\Psi\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}_{\mathrm{tot}}\left(\mathrm{R}\right)\Psi\left(\mathrm{R},\,\mathbf{r}\right)$$

3. For a given complete basis set $\left\{\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\}$, the molecular wavefunction is determined only by the coefficients $c_{k}$. Thus, the coefficients $c_{k}$ depends on nuclear configuration: $$\Psi\left(\mathrm{R},\,\mathbf{r}\right)={\displaystyle \sum_{k}}c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\quad\left(2\right)$$ Now by substituting eq.(1) into eq.(2)

$$\begin{array}{c} {\displaystyle \sum_{k}}\left(\hat{\mathcal{T}}_{\mathrm{N}}+\hat{\mathcal{H}}_{\mathrm{e}}\right)c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}_{\mathrm{tot}}{\displaystyle \sum_{k}}c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\\ {\displaystyle \sum_{k}}\hat{\mathcal{T}}_{\mathrm{N}}c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)+{\displaystyle \sum_{k}}c_{k}\left(\mathrm{R}\right)\hat{\mathcal{H}}_{\mathrm{e}}\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}_{\mathrm{tot}}{\displaystyle \sum_{k}}c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right) \end{array}$$

1. First term of eq.(3) $$\begin{array}{c} {\displaystyle \sum_{k}}\hat{\mathcal{T}}_{\mathrm{N}}c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)=-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ \nabla_{\mathrm{N}_{a}}\cdot\nabla_{\mathrm{N}_{a}}\left[c_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right]\right\} =\\ =-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ \nabla_{\mathrm{N}_{a}}\cdot\left[c_{k}\left(\mathrm{R}\right)\nabla_{\mathrm{N}_{a}}\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)+\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\nabla_{\mathrm{N}_{a}}c_{k}\left(\mathrm{R}\right)\right]\right\} =\\ =-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ c_{k}\left(\mathrm{R}\right)\nabla_{\mathrm{N}_{a}}^{2}\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)+\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\nabla_{\mathrm{N}_{a}}^{2}c_{k}\left(\mathrm{R}\right)+2\left[\nabla_{\mathrm{N}_{a}}c_{k}\left(\mathrm{R}\right)\nabla_{\mathrm{N}_{a}}\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right]\right\} \end{array}$$ Multiplying from the left by $\Psi_{j}^{\star}\left(\mathrm{R},\,\mathbf{r}\right)$ and integrating over electron coordinates, remembering that $\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle =\delta_{jk}$ $$\begin{array}{c} -{\displaystyle \sum_{a}}\dfrac{\hbar^{2}}{2M_{a}}\nabla_{\mathrm{N}_{a}}^{2}c_{k}\left(\mathrm{R}\right)-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ c_{k}\left(\mathrm{R}\right)\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}^{2}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle +2\nabla_{\mathrm{N}_{a}}c_{k}\left(\mathrm{R}\right)\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \right\} =\\ =\hat{\mathcal{T}}_{\mathrm{N}}c_{k}\left(\mathrm{R}\right)-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ c_{k}\left(\mathrm{R}\right)\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}^{2}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle +2\nabla_{\mathrm{N}_{a}}c_{k}\left(\mathrm{R}\right)\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \right\} \end{array}$$

2. Second term of eq.(3)
Multiplying by $\Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)$ and integrating over electron coordinates $${\displaystyle \sum_{k}}c_{k}\left(\mathrm{R}\right)\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\hat{\mathcal{H}}_{\mathrm{e}}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle =\mathcal{E}_{k}\left(\mathrm{R}\right)c_{k}\left(\mathrm{R}\right)$$

3. Final form of the molecular electronic function $$\left(\hat{\mathcal{T}}_{\mathrm{N}}+\mathcal{E}_{k}\left(\mathrm{R}\right)\right)c_{k}\left(\mathrm{R}\right)-{\displaystyle \sum_{a,k}}\dfrac{\hbar^{2}}{2M_{a}}\left\{ C_{2a,jk}c_{k}\left(\mathrm{R}\right)+2C_{1a,jk}\nabla_{\mathrm{N}_{a}}c_{k}\left(\mathrm{R}\right)\right\} =\mathcal{E}_{\mathrm{tot}}c_{k}\left(\mathrm{R}\right)$$
where $$\begin{array}{c} C_{1a,jk}=\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}^{2}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \\ C_{2a,jk}=\left\langle \Psi_{j}\left(\mathrm{R},\,\mathbf{r}\right)\left|\nabla_{\mathrm{N}_{a}}\right|\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)\right\rangle \end{array}$$
By neglecting the coupling terms $$\left(\hat{\mathcal{T}}_{\mathrm{N}}+\mathcal{E}_{k}\left(\mathrm{R}\right)\right)c_{k}\left(\mathrm{R}\right)=\mathcal{E}_{\mathrm{tot}}c_{k}\left(\mathrm{R}\right)$$

Final conclusion

Within this picture, the $c_{k}\left(\mathrm{R}\right)$ behave, in a natural way, as a nuclear eigenfunction. $$\begin{array}{c} \left(\hat{\mathcal{T}}_{\mathrm{N}}+\mathcal{E}_{k}\left(\mathrm{R}\right)\right)\Phi_{k}\left(\mathrm{R}\right)=\mathcal{E}_{\mathrm{tot}}\Phi_{k}\left(\mathrm{R}\right)\\ \hat{\mathcal{H}}_{\mathrm{e}}\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)=\mathcal{E}_{k}\left(\mathrm{R}\right)\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right) \end{array}$$

Questions

1. While computing the Second term of eq.(3), I assumed $\Psi_{k}\left(\mathrm{R},\,\mathbf{r}\right)$ as an eigenfunction of the electronic Hamiltonian. How can this assumption be demonstrated?

2. Do you think that the Assumption 3 is correct?

• Hi Im just wondering how you got your final two equations? I understand your proof up to the point where you neglect the coupling terms but then how do you generate two final results? Didn't you already use the electronic Hamiltonian when you took its eigenvalue $\epsilon _k$ into the nuclear expression you quote? Commented Nov 14, 2014 at 13:31
• @RedPen: The second equation follow from the fact that, as pointed out by Wildcat, we choose to expand, without any approximation, the molecular wave function as $${\displaystyle \sum_{ij}}c_{ij}\Psi_{\mathrm{N},i}\left(\mathbf{R}\right)\Psi_{\mathrm{el},j}\left(\mathbf{r}\right)$$ where the electronic wavefunction is by definition an eigenfunction of the electronic Hamiltonian Commented Nov 14, 2014 at 19:12

With respect to your question, we do not assume that $\Psi_k(R,r)$ is an an eigenfunction of the electronic Hamiltonian, we know that for sure. If you look at equation (2) then essentially what is going on here is the following. Any well-behaved function of two independent variables $\Psi(R,r)$ can always be expanded over the complete set of functions $\{ g_{i}(R) \}$ and $\{ h_{j}(r) \}$ in the following way $$\Psi(R,r) = \sum\limits_{i} \sum\limits_{j} c_{ij} g_{i}(R) h_{j}(r) \, ,$$ or by defining $$c_{j}(R) = \sum\limits_{i} c_{ij} g_{i}(R) \, ,$$ as $$\Psi(R,r) = \sum\limits_{j} c_{j}(R) h_{j}(r) \, .$$ Now as a complete set $\{ h_{j}(r) \}$ of functions of electronic coordinates we can choose the set of electronic wave functions $\Psi_k(R,r)$ which are by definition eigenfunctions of electronic Hamiltonian. The fact that these function carry parametric dependence on nuclear coordinates does not matter here.
• @DavideLaVardera The parametric dependence of electronic wave function on nuclear coordinates is not an assumption but rather a consequence of an assumption that we could treat the motion of a relatively fast moving subsystem (electrons) as if it depends only on the static configuration of a relatively slow moving one (neclei). Once we have made this assumption we obtain the electronic Schrödinger equation and we know that the electronic Hamiltonian is different for different nuclear configurations due to $V_{nn}$ and $V_{en}$ terms. Commented Nov 13, 2014 at 13:55