Faster computation of the Pulay error vector for the direct inversion in the iterative subspace acceleration scheme is achieved by expanding the density matrix in terms of the occupied molecular orbitals and by performing matrix multiplications in an optimal order. The predicted speed-up over the usual implementation is N/O, where N is the size of the atomic orbital basis set, and 0 is number of occupied molecular orbitals, which is confirmed by test calculations. (c) 2005 Elsevier B.V. All rights reserved.