# Are Pulay forces expensive to compute?

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)

• can you please clarify a little what you mean by Pulay forces? Are you referring to Pulay stress or something different? – Richard Terrett May 7 '12 at 3:36
• @RichardTerrett Well, Pulay forces are similar in nature to Pulay stress… even at constant cell volume, the Hellmann-Feynman theorem doesn't hold if the basis set is not fixed, and the extra term appearing is called the Pulay forces. I added a few links, to the original paper and to the first result of a Google search for “Pulay forces”. – F'x May 7 '12 at 6:46
• Stress is used in materials science (haven't really encountered it elsewhere, so correct me if this is an isolated or universal thing) to indicate a force normalized to the cross-section area, so the two are comparable I'd think, right? – Emmie MC May 11 '12 at 21:33

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.

• It turns out that the MO coefficient derivatives are usually by far the more expensive piece. Solving for these analytically is usually done with coupled-perturbed self-consistent field (CPSCF) equations: for KS DFT one has CPKS, and for HF, CPHF, etc. – Jiahao Chen May 11 '12 at 5:04
• In response to AcidFlask: MO coefficient derivatives, whilst non-zero, aren't required for the first derivative of the energy as a result of the energy being minimised with respect to them. See Szabo and Ostlund, p. 440: $$\frac{d E}{d x}=\frac{\partial E}{\partial x}+\sum_{i a}\frac{\partial E}{\partial C_{i a}}\frac{\partial C_{i a}}{\partial x},$$ since $\partial E/\partial C_{i a}=0$, you don't need to compute $\partial C_{ia}/\partial x$. – user346 Jun 18 '12 at 13:15