Are Pulay forces expensive to compute?

Question

In a thesis I am reading, it is said that one of the reasons for using plane-wave basis sets for first-principles molecular dynamics (aka ab initio MD) is that the Pulay forces^[1,2] that arise from an MD using atomic basis sets are computationally expensive to calculate.

While I understand that having additional terms means more code to write, i.e. they make writing the software harder to write, is it true that they are CPU-intensive to compute? The criterion I would use to quantify this subjective statement is:

Given that you have already calculated the energy and forces at that point, you have already computed a large number of integrals required for this task and involving basis functions: overlap integrals, terms of the form $\left\langle\phi_\alpha\left|\hat{A}\right|\phi_\beta\right\rangle$ where operator $\hat{A}$ is either the hamltonian, its gradient, or any other operator necessary for the calculation of energy or forces. Would the calculation of Pulay forces require any more integrals to be evaluated, or can it be straightforwardly computed from those previously-calculated integrals alone?

---

[1] P. Pulay, Molec. Phys. 19, 197 (1969)
[2] See also slide 6 of this

Answer 1

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<\psi\left|H\right|\psi\right>$. Applying the chain rule to this:

$$-F = \nabla E = \left<\psi \left\vert \nabla H \right\vert \psi\right> + 2\left< \nabla\psi \left\vert H \right\vert \psi\right>$$

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<\chi_i^x\vert\chi_j\right>$ 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.