@brief Two-electron derivative-Fock contraction tr(M . F^x[P]) for the CPHF nuclear right-hand side, built natively on top of the validated 2e gradient driver (grd2_driver) with no changes to the Rys internals and no libint.
The 2e gradient driver contracts the derivative ERIs d(uv|ls)/dx with a four-index density product supplied by a grd2_compute_data_t extension. The standard (energy-gradient) extension forms D (x) D. Here we instead form a MIXED product M (x) P, so the same driver returns, for each nuclear coordinate x, g_x = sum_{uvls} d(uv|ls)/dx * [ 4 c M_uv P_ls - x_hf ( M_ul P_vs + M_us P_vl ) ] which is exactly sum_uv M_uv F^x_uv[P] for the closed-shell response Fock F^x[P] = J^x[P] - 1/2 K^x[P] (Coulomb scaled by c, exchange by x_hf=HFscale), summed over the two equivalent index orderings that the driver already exploits. M is the "probe" matrix; for a CPHF RHS element B^x_{ia} the probe is the symmetric AO matrix C_{.,i} C_{.,a}^T + C_{.,a} C_{.,i}^T.
This is the F^x building block of the native CPHF chain. It is validated by the trace identity tr(P . F^x[P]) = (2e part of dE/dx), i.e. against the already-validated grd2_driver energy gradient (exact, non-iterative).