Explicit form of Hartree Fock analytic gradient with contracted GTOs? With or without density fitting

Question

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!

Answer 1

Molecular Electronic-Structure Theory, by Helgaker et al., is an extremely comprehensive textbook reference for electronic structure on finite systems. I take the following from Section 10.8.3 (pp. 481–482).

This section provides an expression for the gradient in terms of the generalized Fock matrix

$$ F_{mn} = \sum_\sigma \langle \text{CSF} \vert a_{m\sigma}^\dagger [ a_{n\sigma}, H] \vert \text{CSF} \rangle; $$

for real orbitals, at least, the gradient can be written as $$ E_{mn}^{(1)} = 2(F_{mn} - F_{nm}). $$

The following equations, (10.8.19–10.8.23), provide expressions for the generalized Fock matrix in terms of one- and two-electron integrals. For the Hamiltonian $$ H = \sum_{pq} h_{pq} E_{pq} + \frac12 \sum_{pqrs} g_{pqrs} e_{pqrs} + h_{nuc}, $$ recall that the Fock matrix is defined in terms of $a_{m\sigma}^\dagger[a_{n\sigma}, H]$, so that $$ F_{mn} = \langle \text{CSF} \vert \sum_\sigma a_{m\sigma}^\dagger[a_{n\sigma}, H] \vert \text{CSF} \rangle = \langle \text{CSF} \vert \sum_q h_{nq} E_{mq} + \sum_{qrs} g_{nqrs} e_{mqrs} \vert \text{CSF} \rangle, $$ then $$ F_{mn} = \sum_q D_{mq} h_{nq} + \sum_{qrs} d_{mqrs} g_{nqrs}, $$ where by definition $$ D_{pq} = \langle \text{CSF} \vert E_{pq} \vert \text{CSF} \rangle, \\ d_{pqrs} = \langle \text{CSF} e_{pqrs} \vert \text{CSF} \rangle.. $$

If you prefer, you can also get the gradient in terms of Roothan–Hall-style equations (Eq. 10.7.80).

Note (Section 2.5 of 1, pp. 51–53) that $\lvert \text{CSF} \rangle$ refers to a configuration state function, which is a spin-adapted modification of the Slater determinant. Construction of these is detailed in Section 2.6.

Answer 2

If you want to use Mathematica, simply write down the expression for the total energy and take the derivative with respect to the nuclear coordinates.

The explicit form of the expressions can be found in standard textbooks and the literature; see e.g. the seminal work of Pulay. All you need is density matrix elements, and the molecular (gradient) integrals. Automated differentiation is another option, and you can also find these works easily in the literature.