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.

  • 1
    $\begingroup$ Technically your Eq. 1 is the just the second-quantized version of my Eq. 1 here: mattermodeling.stackexchange.com/a/6162/5. You could actually just do the integrals on my Eq. 2 there instead of on Eq. 1 (i.e. work with the 2nd quantized version of my Eq. 2). What is your goal with this? You want to do geometry optimization, or DBOC? $\endgroup$ Jun 15, 2021 at 20:14
  • 1
    $\begingroup$ If you're interested in doing DBOC, I believe the first paper to do it at HF level with analytic derivatives was this one by Handy et al.. In 2006 Gauss et al. went beyond HF while still using analytic derivatives. $\endgroup$ Jun 15, 2021 at 20:21
  • $\begingroup$ The goal is to compute the Hessian of the ground state energy function $E(R)$, for which one needs the second derivative of the Hamiltonian given above. $\endgroup$ Jun 16, 2021 at 0:42
  • 2
    $\begingroup$ The geometrical Hessian is available in most electronic structure software packages, especially because of it's use in geometry optimization. OpenMOLCAS can do it. There's a PyMOLCAS that comes with it, and we had extensive discussions last year about transitioning to a more Python-friendly software, but the consensus was not to do it (for various reasons, and you should probably also consider being open to more languages). I've never explicitly asked a program to print me the Hessian, but it's used all the time for geometry opt. $\endgroup$ Jun 16, 2021 at 16:03
  • $\begingroup$ Are you able to elaborate on why you want to calculate the Hessian, rather than just doing the geometry optimization (or whatever you need the Hessian for)? $\endgroup$ Jun 21, 2021 at 23:12

1 Answer 1


PySCF and Psi4 both implement analytical geometric Hessians for Hartree-Fock calculations. The computation of the geometric Hessian indeed requires solving the coupled perturbed Hartree-Fock equations: while one can evaluate the force acting on the nuclei just with the SCF wave function and the integral derivatives, the second derivative also contains a response term for the relaxation of the orbitals upon the perturbation.

You also stated interest in getting the derivative integrals. These are naturally also available in both PySCF and Psi4; the former computes integrals with the libcint library by Qiming Sun (original author of PySCF), the latter uses the libint library by Ed Valeev.

  • $\begingroup$ +1. Pretty much every major electronic structure software will have the Hessian for geometry optimization (OpenMOLCAS and GAMESS were mentioned in the comments), but I suppose the OP would have to go through the code and add print statements to output the elements of the Hessian matrix, since this is not typically printed? $\endgroup$ Jun 28, 2021 at 22:55
  • 1
    $\begingroup$ @NikeDattani e.g. Psi4 returns the Hessian as a matrix, so no code modifications are needed psicode.org/psi4manual/master/freq.html $\endgroup$ Jun 29, 2021 at 23:21
  • $\begingroup$ Excellent, I suppose the OP has the answer to their question then! Does PySCF not print out the Hessian without the user having to add more print statements? $\endgroup$ Jun 29, 2021 at 23:23
  • 1
    $\begingroup$ PySCF also returns the Hessian as a matrix github.com/pyscf/pyscf.github.io/blob/master/examples/hessian/… $\endgroup$ Jun 30, 2021 at 17:04

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .