Does anyone know of a paper, textbook, or other resource which outlines the expression for the analytic gradient of the HF energy (restricted or otherwise) in terms of Dunning's contracted Gaussian-type orbitals $\chi^{CGTO}_\mu(r)$? I'm talking about an expression that someone could simply sit down and program as a free-standing code using only a converged HF calculation as input. I.e. has someone written this down in a form where the only arguments of the expression would be the primitive GTO exponents, the contraction coefficients, the atom centers, the nuclear charges, and the converged eigenvalues and eigenvectors of the Hartree Fock SCF?
Plenty of textbooks (Szabo and Ostlund, Yamaguchi et al, etc) describe the theory of analytic gradients using equations independent of basis set type, i.e. in an abstract way. Likewise, plenty of standard codes (open source or otherwise) perform these calculations with Dunning CGTOs and/or other basis set types. But because of how huge and modular (i.e. calling across many routines to read) these codes are, interpreting them takes a great deal of time.
The form of the derivative of primitive GTOs is given by Szabo and Ostlund page 442 between equations C.12 and C.13. Since the contraction coefficients are not spatially dependent, this form should be (hopefully) simple to extend to a Dunning basis. From here, the form of the various derivative integrals (such as the matrix elements $\frac{\partial V_{N N}}{\partial X_A}_{\mu \nu}$ and $\frac{\partial V_{N e}}{\partial X_A}_{\mu \nu}$ given in the middle of page 442) to obtain the HF gradient should be doable analytically in Mathematica etc. Thus, the explicit form of all these integrals and the overall gradient expression should be possible to write down in a way someone (maybe even someone without electronic structure knowledge) could just sit down and code. The expression would be long, but easily interpreted.
I'm just wondering if this full expression exists anywhere, as I haven't been able to find it anywhere online and would like to before "recreating the wheel." In particular, I will be using output from a density-fitted Hartree Fock code for my analytic gradient implementation, but any expression will do, as extending to density-fit Coulomb and exchange is probably the least tedious part of this process. Thank you!