Laplacian of the Hartree Fock second-order reduced density matrix

I am attempting to replicate and understand the Colle-Salvetti method for approximating the electron correlation energy.[1] I have successfully replicated the method, but there is one part of the derivation that I cannot correctly derive, the Laplacian of the Hartree–Fock, second-order reduced density matrix. 


With some trial and error I was able to modify my derived form for the Laplacian operator to then exactly match the correlation energies calculated by Colle and Salvetti for the helium atom. My derivation of the following term is as follows
$$
\nabla^2\_\mathbf{r} P\_{2 \text{HF}}\left(\mathbf{R} - \frac{\mathbf{r}}{2}, \mathbf{R} + \frac{\mathbf{r}}{2}\right)\_{r=0},
\tag1
$$
where $\mathbf{R} = \frac{1}{2}(\mathbf{r}\_1 + \mathbf{r}\_2)$ and $\mathbf{r}=|\mathbf{r}\_1 - \mathbf{r}\_2|$. I found a useful article written by Colle and Salvetti in 1983,[2] which states a form for $\nabla^2\_\mathbf{r}$ in terms of the inter-particle coordinates $\mathbf{r\_1}$ and $\mathbf{r\_2}$
$$
\nabla\_\mathbf{r}^2 = \frac{1}{4}\left(\nabla\_{\mathbf{r\_1}}^2 + \nabla\_{\mathbf{\mathbf{r\_2}}}^2 - 2 \nabla\_{\mathbf{r\_1}} \cdot \nabla\_{\mathbf{\mathbf{r\_2}}}\right). 
\tag2\label{eq:two}
$$


This form is useful as the Hartree Fock wavefunction for helium I am using is in terms of $\mathbf{r\_1}$ and $\mathbf{r\_2}$. I apply the chain rule to the derivatives in terms of $\mathbf{r\_1}$ and $\mathbf{r\_2}$
\begin{aligned}
 \frac{\partial}{\partial \mathbf{r\_1}} 
 &= \frac{\partial \mathbf{r}}{\partial \mathbf{r\_1}} 
 \frac{\partial }{\partial \mathbf{r}} 
 + \frac{\partial \mathbf{R}}{\partial \mathbf{r\_1}} 
 \frac{\partial }{\partial \mathbf{R}} \\ 
 &= \frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}},
\end{aligned}
and
\begin{aligned}
 \frac{\partial}{\partial \mathbf{r\_2}} 
 &= \frac{\partial \mathbf{r}}{\partial \mathbf{r\_2}} 
 \frac{\partial }{\partial \mathbf{r}} 
 + \frac{\partial \mathbf{R}}{\partial \mathbf{r\_2}} 
 \frac{\partial }{\partial \mathbf{R}} \\ 
 &= -\frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}.
\end{aligned}


The second order derivatives take the form
\begin{aligned}
 \frac{\partial^2}{\partial \mathbf{r\_1}^2} 
 &= \left(\frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}\right)
 \left(\frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}\right) \\
 &= \frac{\partial^2}{\partial \mathbf{r}^2} 
 + \frac{\partial^2}{\partial \mathbf{r} \partial \mathbf{R}} 
 + \frac{1}{4}\frac{\partial^2}{\partial \mathbf{R}^2},
\end{aligned}
and
\begin{aligned}
 \frac{\partial^2}{\partial \mathbf{r\_2}^2} 
 &= \left(-\frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}\right)
 \left(-\frac{\partial}{\partial \mathbf{r}} 
 + \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}\right) \\
 &= \frac{\partial^2}{\partial \mathbf{r}^2} 
 - \frac{\partial^2}{\partial \mathbf{r} \partial \mathbf{R}} 
 + \frac{1}{4}\frac{\partial^2}{\partial \mathbf{R}^2}.
\end{aligned}


The Laplacian is evaluated when the inter-electronic distance, $\mathbf{r}$ is $0$, thus all these derivatives simplify
$$
\frac{\partial}{\partial \mathbf{r\_1}} 
 = \frac{\partial}{\partial \mathbf{r\_2}} 
 = \frac{1}{2}\frac{\partial}{\partial \mathbf{R}},
$$
and
$$
\frac{\partial^2}{\partial \mathbf{r\_1}^2} 
 = \frac{\partial^2}{\partial \mathbf{r\_2}^2} 
 = \frac{1}{4}\frac{\partial^2}{\partial \mathbf{R}^2},
$$


Substituting these into equation \eqref{eq:two} gives
\begin{aligned}
 \nabla\_\mathbf{r}^2 
 &= \frac{1}{4}\left[
 \frac{1}{4}\frac{\partial^2}{\partial \mathbf{R}^2} 
 + \frac{1}{4}\frac{\partial^2}{\partial \mathbf{R}^2} 
 - 2\left(\frac{1}{2}\frac{\partial}{\partial \mathbf{R}}\right) 
 \cdot \left(
 \frac{1}{2}\frac{\partial}{\partial \mathbf{R}}
 \right)
 \right] \\
 &= \frac{1}{4}\left[
 \frac{1}{2}\frac{\partial^2}{\partial \mathbf{R}^2} 
 - \frac{1}{2}\left(
 \frac{\partial}{\partial \mathbf{R}}
 \right)^2
 \right]
\end{aligned}


Now act this operator on the two-electron reduced density matrix
$$
\nabla^2\_\mathbf{r} P\_{2 \text{HF}}\left(\mathbf{R} 
 - \frac{\mathbf{r}}{2}, \mathbf{R} + \frac{\mathbf{r}}{2}\right)\_{r=0} 
= \frac{1}{4}\left[
 \frac{1}{2}\frac{\partial^2}{\partial \mathbf{R}^2} 
 - \frac{1}{2}\left(
 \frac{\partial}{\partial \mathbf{R}}
 \right)^2
 \right] 
 \rho\_2^\text{HF} (\mathbf{r\_1}, \mathbf{r\_2})\rvert\_{
 \mathbf{r\_1} = \mathbf{r\_2} = \mathbf{R}
 },
$$
where I assume the two-electron density is a product of the two, one-electron densities
$$
\rho\_2^\text{HF} (\mathbf{r\_1}, \mathbf{r\_2})\rvert\_{\mathbf{r\_1} 
= \mathbf{r\_2} 
= \mathbf{R}} 
= \frac{1}{2}\rho(\mathbf{r\_1})\rho(\mathbf{r\_2})\rvert\_{\mathbf{r\_1} 
= \mathbf{r\_2} 
= \mathbf{R}} 
= \frac{1}{2}\rho(\mathbf{R})\rho(\mathbf{R})
$$


We now have
$$
\nabla^2\_\mathbf{r} P\_{2 \text{HF}}\left(\mathbf{R}, \mathbf{R}\right) 
= \frac{1}{8}\left[\frac{1}{2}
 \frac{\partial^2 \rho(\mathbf{R})\rho(\mathbf{R})}{\partial \mathbf{R}^2} 
 - \frac{1}{2}\left(\frac{\partial \rho(\mathbf{R})\rho(\mathbf{R})}{\partial \mathbf{R}}\right)^2
 \right]
$$


Apply the product rule and factor out $\rho(\mathbf{R})$
$$
\nabla^2\_\mathbf{r} P\_{2 \text{HF}}\left(\mathbf{R}, \mathbf{R}\right) 
= \frac{1}{8}\rho(\mathbf{R})\left[
 \frac{\partial^2 \rho(\mathbf{R})}{\partial \mathbf{R}^2} 
 + \frac{1}{\rho(\mathbf{R})}\left(
 \frac{\partial \rho(\mathbf{R})}{\partial \mathbf{R}}
 \right)^2 
 - 2\rho(\mathbf{R})\left( 
 \frac{\partial \rho(\mathbf{R})}{\partial \mathbf{R}}
 \right)^2
 \right]
$$


This gives the incorrect answer. However by dividing the last term, $2\rho(\mathbf{R})\left( \frac{\partial \rho(\mathbf{R})}{\partial \mathbf{R}}\right)^2$ by $\rho(\mathbf{R})^2$ allows for Colle and Salvetti's correlation energy to be replicated exactly; and you can combine the $\left( \frac{\partial \rho(\mathbf{R})}{\partial \mathbf{R}}\right)^2$ terms. As this example is for the helium atom there is a spherical symmetry to exploit so the vector notation is not strictly required. 


I cannot see what I have done wrong in this derivation. It is very close to being correct with a seemingly rogue $\rho(\mathbf{R})^2$ in the final term. Any assistance is greatly appreciated in finding what I have done wrong.


### References:


1. Colle, R.; Salvetti, O. Approximate calculation of the correlation energy for the closed shells. *Theor. Chim. Acta* **1975,** *37* (4), 329–334. [DOI: 10.1007/BF01028401](https://doi.org/10.1007/BF01028401).
2. Colle, R.; Salvetti, O. A general method for approximating the electronic correlation energy in molecules and solids. *J. Chem. Phys.* **1983,** *79* (3), 1404–1407. [DOI: 10.1063/1.445899](https://doi.org/10.1063/1.445899).


