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