  We propose an algorithm for accurate, systematic and scalable computation of
interatomic forces within the auxiliary-field Quantum Monte Carlo (AFQMC)
method. The algorithm relies on the Hellman-Fenyman theorem, and incorporates
Pulay corrections in the presence of atomic orbital basis sets. We benchmark
the method for small molecules by comparing the computed forces with the
derivatives of the AFQMC potential energy surface, and by direct comparison
with other quantum chemistry methods. We then perform geometry optimizations
using the steepest descent algorithm in larger molecules. With realistic basis
sets, we obtain equilibrium geometries in agreement, within statistical error
bars, with experimental values. The increase in computational cost for
computing forces in this approach is only a small prefactor over that of
calculating the total energy. This paves the way for a general and efficient
approach for geometry optimization and molecular dynamics within AFQMC.
