I'm currently writing a program which evaluates the first derivatives of molecular integrals (over gaussian basis sets) with respect to the cartesian coordinates of the nucleus. I'm somewhat stuck with the derivative of the nuclear-electron attraction integral, i.e.: $$ \frac{\partial V_{\mu\nu}}{\partial X_A}=\frac{\partial }{\partial X_A} \left\langle \mu| \hat V |\nu\right\rangle $$ where $\mu$ and $\nu$ are Gaussian-Type Orbitals (GTOs) and $\hat V $ is the operator: $$ \hat V = -\sum_A \frac{Z_A}{|r-R_A|} $$ As was previously discussed on this forum, the upper derivative can be expressed as a sum of derivatives applying the chain rule: $$ \frac{\partial }{\partial X_A} \left\langle \mu| \hat V |\nu\right\rangle=\left\langle\frac{\partial \mu}{\partial X_A}|\hat V|\nu\right\rangle+\left\langle \mu|\frac{\partial \hat V}{\partial X_A}|\nu\right\rangle+\left\langle \mu|\hat V|\frac{\partial \nu}{\partial X_A}\right\rangle$$ The first and latter terms are relatively easy to evaluate, since the derivative of a Gaussian-Type Orbital with respect to one of the coordinates of the point where it's centered is a sum of GTOs with higher and lower angular momentum. The evaluation of the middle term doesn't seem so obvious. According to Szabo & OustlundOstlund (pg. 442), the derivative of the $\hat V$ operator is: $$ \frac{\partial \hat V}{\partial X_A}=-Z_A\sum_i \frac{X_i-X_A}{r_{iA}^3} $$ where $\mu$ and $\nu$ are Gaussian-Type Orbitals (GTOs) So we are dealing with a new type of molecular integral: $$ \left\langle \mu\left|-Z_A\sum_i \frac{X_i-X_A}{r_{iA}^3}\right|\nu\right\rangle $$ After searching for some time, I couldn't find any reference which discuss this topic. Could you provide some insight on the evaluation of this particular integral? Thanks in advance.