Are Pulay forces expensive to compute?
Yes, you need additional quantities beyond the minimum necessary for calculating energies and the Hellmann-Feynman piece of the force when the wavefunction is not variational.

Here is a very rough sketch of why:

Puláy forces arise from applying the chain rule for calculating forces and was first discussed in the context of applying the Hellman-Feynman theorem.

Recall that the force is the change of energy with changes in (nuclear) coordinate and the energy is the expectation value of the (electronic) Hamiltonian, $E = \left&lt;\psi\left|H\right|\psi\right&gt;$. Applying the chain rule to this:

$$-F = \nabla E = \left&lt;\psi \left\vert \nabla H \right\vert \psi\right&gt; +  2\left&lt; \nabla\psi \left\vert H \right\vert \psi\right&gt;$$

The first term is what you get from Hellman-Feynman and is the expectation of a one-electron operator of the first you have listed. The second term goes away in Hellmann-Feynman only because it assumes the wavefunction is variational.

If you expand out the orbitals in the Slater determinant in some AO basis $\chi$:

$$ \psi(r_1, r_2, ..., r_N) = \phi(r_1) \wedge \phi(r_2) \wedge \dots \wedge \phi(r_N) $$
$$ \phi(r_1) = \sum_i c_i \chi_i(r_i) $$

then it is clear that $\nabla\psi$ generates terms which are MO coefficient derivatives (usually denoted $c_i^x$) and also AO derivatives ($\chi_i^x$). If the wavefunction were variational, then by definition all these derivatives are zero. If you work through the algebra, you will find that there are some new quantities that need to be calculated in order to work out the non-Hellmann-Feynman theorem, most notably the one-sided AO overlap derivative matrices $\left&lt;\chi_i^x\vert\chi_j\right&gt;$ which do not normally show up otherwise. These terms turn out to be critical for getting correct ab initio molecular dynamics since the wavefunction is rarely (if ever) variational in time.

There are tricks one can do to reduce the cost of calculating MO coefficient derivatives , but they cannot be obtained for free in general. And there's really no way around calculating the overlap derivative terms.