Analytical gradient of the Hartree-Fock nuclear-electronic repulsion term
There are a few points to discuss:


Since there are $3N$ possible $\{X_{A}\}$, each term where an $X_{A}$ appears will result in $3N$ matrix elements. In your first equation, that will be $S$, $T$, $V$, and the 2-electron contribution, which should really be rewritten into a Fock matrix term $F$ where the ket is precontracted with the density.
There is nothing wrong with this solution, though it is a little strange because it's in the MO basis, and we normally prefer to work in the AO basis to avoid any unnecessary transformations. I'm going to drop the $Ne$ from here on out since $V$ is understood in this context to only be the one-electron nuclear attraction term.
I think you're asking what happens to the density matrix. The derivative of $V$ expands into multiple parts, as we will see below. I will reference a bunch of equations that use $h$, where $\left( h = H = H^{\text{core}} \right) \equiv T_{e} + V_{Ne}$, and differentiating into a sum should hopefully be clear. What follows should work for any real-valued one-electron operator. I will use the book by Yamaguchi with liberties taken. In many places I use $p,q$ rather than $i,j$, as these indices run over all MOs, not just occupied ones.




Expanding $\frac{\partial V_{ij}}{\partial X_{A}}$ using the product rule gives

$$
\begin{align*}
\frac{\partial V_{ij}}{\partial X_{A}} &amp;= \frac{\partial}{\partial X_{A}} \left( \sum_{\mu\nu}^{\text{AO}} C_{\mu i} C_{\nu j} V_{\mu\nu} \right) \tag{3.80} \\
&amp;= \sum_{\mu\nu}^{\text{AO}} \left( \frac{\partial C_{\mu i}}{\partial X_{A}} C_{\nu j} V_{\mu\nu} + C_{\mu i} \frac{\partial C_{\nu j}}{\partial X_{A}} V_{\mu\nu} + C_{\mu i} C_{\nu j} \frac{\partial V_{\mu\nu}}{\partial X_{A}} \right), \label{3.81}\tag{3.81}
\end{align*}
$$

where the third/last term is the true AO integral derivative, and the first two terms, the MO coefficient derivatives, come from differentiating the density matrix.



The simplest term to write out, but the most difficult to implement, is the AO integral derivative. I will use $\mu,\nu$ rather than $\chi_{\mu},\chi_{\nu}$, so they refer to both AO basis functions and their matrix indices.

$$
\begin{align*}
\frac{\partial V_{\mu\nu}}{\partial X_{A}} &amp;= \frac{\partial}{\partial X_{A}} \left&lt; \mu | \hat{V} | \nu \right&gt; \tag{3.24} \\
&amp;= \left&lt; \frac{\partial \mu}{\partial X_{A}} | \hat{V} | \nu \right&gt; + \left&lt; \mu | \frac{\partial \hat{V}}{\partial X_{A}} | \nu \right&gt; + \left&lt; \mu | \hat{V} | \frac{\partial \nu}{\partial X_{A}} \right&gt; \label{3.25}\tag{3.25}
\end{align*}
$$

What you have presented,

$$
\frac{\partial \hat{V}}{\partial X_{A}} = -Z_{A} \sum_{i} \frac{X_{i} - X_{A}}{r_{iA}^{3}},
$$

is only the derivative of the operator. The basis function derivatives are also required, as the AOs in this case depend on the nuclear positions. Contrast this with an electric field perturbation, where only the derivative of the operator is required, or magnetic field derivatives, where basis functions may be perturbation-dependent if gauge-including atomic orbitals (GIAOs) are used. Also important is that the index $i$ here refers to an electron, not an occupied MO. This is an electric field integral.



Now for the derivative of the MO coefficients/density matrix. Clearly they don't appear in the final HF gradient expression. They disappear through the magic of Wigner's $2n + 1$ rule. From page 25:


  When the wavefunction is determined up to the $n$th order, the expectation value (electronic energy) of the the system is resolved, according to the results of perturbation theory, up to the $(2n+1)$st order. This principle is called Wigner's $2n+1$ theorem [29-31].


More explicitly, we have the zeroth-order wavefunction, so we must be able to calculate the first-order correction to the energy. Worded differently, any first derivative of the energy must be easily calculated without differentiating MO coefficients, which is only required for second derivatives, such as the molecular Hessian or the dipole polarizability.

First, rewrite the MO coefficient derivatives

$$
\frac{\partial C_{\mu i}}{\partial X_{A}} = \sum_{m}^{\text{MO}} U_{mi}^{X_{A}} C_{\mu m}, \tag{3.7}
$$

where the index $m$ runs over all occupied and unoccupied/virtual MOs. The same goes for $p,q$. I've replaced $a$ from the text with $X_{A}$, but this holds for any general perturbation (such as $\lambda$, pick your favorite unused index). The key insight is that we can write the effect of a perturbation on the MO coefficients as the contraction of the unmodified MO coeffcients with a unitary matrix describing single-particle excitations from occupied to virtual MOs, as well as deexcitations from virtual to occupied MOs. In matrix form, this is

$$
\mathbf{C}^{(X_{A})} = \mathbf{C}^{(0)} \left( \mathbf{U}^{(X_{A})} \right)^{T},
$$

where it is usually easiest to have the dimension of $\mathbf{U}$ be $[N_{\text{orb}}, N_{\text{orb}}]$, with only the occ-virt and virt-occ blocks being non-zero. The problem now becomes solving for $\mathbf{U}$ and eliminating it from the final gradient expression. I will walk through parts of section 4.3 to show how this is done.

Given the energy expression for an RHF wavefunction,

$$
E = 2 \sum_{i}^{\text{d.o.}} h_{ii} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj) - (ij|ij) \right], \tag{4.1}
$$

the first derivative with respect to $X_{A}$ is

$$
\frac{\partial E_{\text{elec}}}{\partial X_{A}} = 2 \sum_{i}^{\text{d.o.}} h_{ii}^{X_{A}} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj)^{X_{A}} - (ij|ij)^{X_{A}} \right] + 4 \sum_{m}^{\text{all}} \sum_{i}^{\text{d.o.}} U_{mi}^{X_{A}} F_{im}, \label{4.16}\tag{4.16}
$$

where the Fock matrix is defined as

$$
\begin{align*}
F_{pq} &amp;= h_{pq} + \sum_{k}^{\text{d.o.}} \left[ 2(pq|kk) - (pk|qk) \right] \tag{4.6} \\
&amp;= h_{pq} + 2J_{pq} - K_{pq},
\end{align*}
$$

and I've introduced the Coulomb and exchange matrices $\mathbf{J}$ and $\mathbf{K}$ as well. I skipped all the steps in expanding the first derivative, and $\eqref{4.16}$ is what results upon collecting terms with $\mathbf{U}$.

Using the RHF variational conditions, the Fock matrix from a converged calculation is diagonal in the MO basis, corresponding to the MO energies

$$
F_{pq} = \delta_{pq} \epsilon_{pq}, \tag{4.7}
$$

so $\eqref{4.16}$ simplifies to

$$
\frac{\partial E_{\text{elec}}}{\partial X_{A}} = 2 \sum_{i}^{\text{d.o.}} h_{ii}^{X_{A}} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj)^{X_{A}} - (ij|ij)^{X_{A}} \right] + 4 \sum_{m}^{\text{all}} \sum_{i}^{\text{d.o.}} U_{mi}^{X_{A}} \epsilon_{im}, \tag{4.17}
$$

which can be further simplified as

$$
\frac{\partial E_{\text{elec}}}{\partial X_{A}} = 2 \sum_{i}^{\text{d.o.}} h_{ii}^{X_{A}} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj)^{X_{A}} - (ij|ij)^{X_{A}} \right] + 4 \sum_{i}^{\text{d.o.}} U_{ii}^{X_{A}} \epsilon_{ii}. \tag{actual 4.17}
$$

Now we use one of the most important tricks in quantum chemistry. Given the orthonormality of the MOs,

$$
S_{pq} = \delta_{pq}, \tag{3.44}
$$

we must have

$$
\frac{\partial S_{pq}}{\partial X_{A}} \overset{!}{=} 0. \tag{3.45}
$$

This is where not using general notation in the above is a bit confusing because it seems to conflict with your original expression, but remember that similar to $\eqref{3.81}/\eqref{3.25}$, this is in fact multiple terms: two for the basis functions, and two for the MO coefficients, giving

$$
\begin{align*}
\frac{\partial S_{pq}}{\partial X_{A}} &amp;= \sum_{\mu\nu}^{\text{AO}} C_{\mu p} C_{\mu q} \frac{\partial S_{\mu\nu}}{\partial X_{A}} + \sum_{m}^{\text{all}} \left( U_{mp}^{X_{A}} S_{mq} + U_{mq}^{X_{A}} S_{pm} \right) \tag{3.40 + 3.43} \\
                                       &amp;= S_{pq}^{X_{A}} + \sum_{m}^{\text{all}} \left( U_{mp}^{X_{A}} S_{mq} + U_{mq}^{X_{A}} S_{pm} \right). \tag{3.43}
\end{align*}
$$

The sum over all MOs can be eliminated by reusing the orthonormality condition, so in the first term $m \overset{!}{=} q$ and for the second term $m \overset{!}{=} p$, and the overlap matrix in the MO basis is unity for those terms, giving

$$
\frac{\partial S_{pq}}{\partial X_{A}} = S_{pq}^{X_{A}} + U_{qp}^{X_{A}} + U_{pq}^{X_{A}} \overset{!}{=} 0. \tag{3.46}
$$

Recognizing that we only need diagonal terms, this can be rewritten as

$$
U_{pp}^{X_{A}} = -\frac{1}{2} S_{pp}^{X_{A}}, \tag{4.20}
$$

which is then plugged back into the first derivative expression to give

$$
\begin{align*}
\frac{\partial E_{\text{elec}}}{\partial X_{A}} &amp;= 2 \sum_{i}^{\text{d.o.}} h_{ii}^{X_{A}} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj)^{X_{A}} - (ij|ij)^{X_{A}} \right] + 4 \sum_{i}^{\text{d.o.}} \left( -\frac{1}{2} S_{ii}^{X_{A}} \right) \epsilon_{ii} \\
                                                &amp;= 2 \sum_{i}^{\text{d.o.}} h_{ii}^{X_{A}} + \sum_{ij}^{\text{d.o.}} \left[ 2(ii|jj)^{X_{A}} - (ij|ij)^{X_{A}} \right] - 2 \sum_{i}^{\text{d.o.}} S_{ii}^{X_{A}} \epsilon_{ii}. \tag{4.21}
\end{align*}
$$

Rewrite the last term

$$
\begin{align*}
\sum_{i}^{\text{d.o.}} S_{ii}^{X_{A}} \epsilon_{ii} &amp;= \sum_{i}^{\text{d.o.}} \sum_{\mu\nu}^{\text{AO}} C_{\mu i} C_{\mu i} \frac{\partial S_{\mu\nu}}{\partial X_{A}} \epsilon_{ii} \\
                                                    &amp;= \sum_{i}^{\text{d.o.}} \sum_{\mu\nu}^{\text{AO}} C_{\mu i} C_{\mu i} \epsilon_{ii} \frac{\partial S_{\mu\nu}}{\partial X_{A}} \\
                                                    &amp;= \sum_{\mu\nu}^{\text{AO}} W_{\mu\nu} \frac{\partial S_{\mu\nu}}{\partial X_{A}} \tag{4.24}
\end{align*}
$$

to use the energy-weighted density matrix $\mathbf{W}$, which in your expression is called $\mathbf{Q}$.

Again, the elimination of the $\mathbf{U}$ matrix is one of the most important results in quantum chemistry, as it means the coupled-perturbed SCF equations do not need to be solved for first derivatives of SCF wavefunctions. This is why you do not see any density or MO coefficient derivatives in the gradient expression.

References


Yamaguchi, Yukio; Goddard, John D.; Osamura, Yoshihiro; Schaefer III, Henry F. A New Dimension to Quantum Chemistry: Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory (International Series of Monographs on Chemistry); Oxford University Press: 1994.
Epstein, Saul T. General Remainder Theorem. J. Chem. Phys. 1968, 48, 4725.
Epstein, Saul T. Constraints and the $V^{2n+1}$ theorem. Chem. Phys. Lett. 1980, 70, 311.


If you don't have access to the Yamaguchi book, the seminal article on HF derivatives contains a condensed version of this derivation, but be aware of typos.


Pople, John A.; Raghavachari, Krishnan; Schlegel, H. Bernhard (Berny); Binkley, J. Stephen. Derivative Studies in Hartree-Fock and Møller-Plesset Theories. Int. J. Quantum Chem. 1979, 16 (num. S13), 225-241.