Once solutions of the electronic Schrödinger equation {ψe,i(r⃗ e)}∞i=1 and {Ee,i}∞i=1 are found for a number of different fixed nuclear configurations r⃗ n, we took the second step: the task is to solve an equation of motion for nuclei which is obtained as follows. Note first that since H^e is self-adjoint, its eigenfunctions {ψe,i(r⃗ e,[r⃗ n])}∞i=1 form a complete orthonormal set. Now we use the assertion that a function of two independent variables ψ(r⃗ e,r⃗ n) can be expanded over a complete set of functions of one variable {ψe,i(r⃗ e,[r⃗ n])}∞i=1 in the following way ψ(r⃗ e,r⃗ n)=∑i=1∞ψe,i(r⃗ e,[r⃗ n])ψn,i(r⃗ n).(2) Substitution of (2) in the Schrödinger equation (1), multiplication of both sides by ψ¯¯¯e,j, and integration over electronic coordinates r⃗ e only, leads to

(T^n+Ee,j+V^nn−E)ψn,j−∑i=1∞Λ^jiψn,i=0for j=1,2,…,(3)

where operator Λ^ji is defined as follows

Λ^ji=∑α=1ν12mα(2⟨ψe,j|∇α|ψe,i⟩∇α+⟨ψe,j|∇2α|ψe,i⟩).

and arguments r⃗ e and r⃗ n were dropped for brevity. The last equation describes a system of differential equations (one equation for each value of j), and once the system is solved and a set of functions ψn,i is obtained, one could write down the exact solution of the Schrödinger equation (1) in the form we introduced above (2).

Operators Λ^ji may be regarded as the elements of a square matrix Λ and of course the fact that we equate these elements with zero is an approximation: in the first case, when we equate with zero all off-diagonal elements Λ^ji, we introduce the adiabatic approximation, while in the second case, when we equate with zero all elements Λ^ji, we introduce the Born–Oppenheimer approximation.