I'm currently working on a quantum chemistry project, I am dealing with a second-quantized electronic Hamiltonian (corresponding to some molecule), written in terms of fermionic creation and annihilation operators of the form:
$$H(R) = \displaystyle\sum_{p, q} h_{pq}(R) c_p^{\dagger} c_q + \frac{1}{2} \displaystyle\sum_{p, q, r, s} h_{pqrs}(R) c^{\dagger}_p c^{\dagger}_q c_r c_s\tag{1}$$
where $R$ is the nuclear coordinates of the molecule, and:
$$h_{pq} = \displaystyle\int dx \phi_p(x)^{*} \left( -\frac{\nabla_r^2}{2} - \displaystyle\sum_{I} \frac{Z_I}{|r - R_I|} \right) \phi_q(x)\tag{2}$$
$$h_{pqrs}(R) = \displaystyle\int dx_1 dx_2 \frac{\phi_p(x_1)^{*} \phi_q(x_2)^{*} \phi_r(x_2) \phi_s(x_1)}{|r_1 - r_2|}\tag{3}$$
where each $\phi$ is a molecular orbital, obtained from performing Hartree-Fock for the molecular Hamiltonian of interest (which depends implicitly on $R$). Computing these orbitals, and therefore the integrals, is straightforward through use of most electronic structure program suites.
I am interested in calculating the first and second derivatives of $H(R)$, with respect to $R$, so in effect, I need to compute the first and second derivatives of the integral outlined above. Is it possible to use PySCF or Psi4 for this task as well? I've looked around, and it seems like this problem of computing derivatives has to do with solving the coupled-perturbed Hartree-Fock equations, although I am unsure.
Addendum: I did some more looking around, and it seems like GAMESS might be able to find the derivatives I'm looking for, as it can do analytical calculations of the Hessian, but I'm still unsure how to access these derivatives in the output file.