Reading few sources on the Born–Oppenheimer approximation I don't understand one particular thing.
If you look for example here (PDF, 70 KB) and focus attention on equations 14 and 15 than it is clear that
$$ \nabla_{A}^2 \left( \psi_k(\mathbf{r}; \mathbf{R}) \chi_k(\mathbf{R}) \right) = \psi_k(\mathbf{r}; \mathbf{R}) \nabla_{A}^2 \chi_k(\mathbf{r}; \mathbf{R}) + 2 \nabla_{A} \psi_k(\mathbf{r}; \mathbf{R}) \nabla_{A} \chi_k(\mathbf{r}; \mathbf{R}) + \chi_k(\mathbf{R}) \nabla_{A}^2 \psi_k(\mathbf{r}; \mathbf{R}) $$
where
$$ \nabla_A^2 = \frac{\partial^2}{\partial X_A^2} + \frac{\partial^2}{\partial Y_A^2} + \frac{\partial^2}{\partial Z_A^2} $$
and $$ \mathbf{R} = \{ \mathbf{R_i} \}_{i=1}^N = \{ (X_i, Y_i,Z_i) \}_{i=1}^N $$ is a set of all nuclear coordinates.
But honestly I don't undestand why is that so. The fact is that $\psi_k$ depends implicitly only on $\mathbf{r}$ and parametrically on $\mathbf{R}$ (that's why I think they are delimited by $;$ and not just $,$). As far as I know this parametric dependence means that for each set of nuclear coordinates $\mathbf{R}$ there is a complete set of electronic wave functions $\{ \psi_k(\mathbf{r}) \}_{k}$ which are functions of electronic coordinates only. And then of course when you differentiate $\psi_k(\mathbf{r}) \chi_k(\mathbf{R})$ twice with respect to $\mathbf{R_A}$ you get just $\psi_k(\mathbf{r}) \nabla_{A}^2 \chi_k(\mathbf{R})$ because $\psi_k(\mathbf{r})$ is constant with respect to $\mathbf{R}$.
And one more thing - the linked resource (and many others) claimed that the chain rule is used going from 14 to 15. I don't see any usage of chain rule but I see a usage of product rule.
Seems like I don't understand what's going on here but this is critical step because non-adiabatic coupling terms come from this expansion.