  Accurate treatment of the long-range electron correlation energy, including
van der Waals (vdW) or dispersion interactions, is essential for describing the
structure, dynamics, and function of a wide variety of systems. Among the most
accurate models for including dispersion into density functional theory (DFT)
is the range-separated many-body dispersion (MBD) method [A. Ambrossetti et
al., J. Chem. Phys. 140, 18A508 (2014)], in which the correlation energy is
modeled at short-range by a semi-local density functional and at long-range by
a model system of coupled quantum harmonic oscillators. In this work, we
develop analytical gradients of the MBD energy with respect to nuclear
coordinates, including all implicit coordinate dependencies arising from the
partitioning of the charge density into Hirshfeld effective volumes. To
demonstrate the efficiency and accuracy of these MBD gradients for geometry
optimizations of systems with intermolecular and intramolecular interactions,
we optimized conformers of the benzene dimer and isolated small peptides with
aromatic side-chains. We find excellent agreement with the wavefunction theory
reference geometries of these systems (at a fraction of the computational cost)
and find that MBD consistently outperforms the popular TS and D3(BJ) dispersion
corrections. To demonstrate the performance of the MBD model on a larger system
with supramolecular interactions, we optimized the C$_{60}$@C$_{60}$H$_{28}$
buckyball catcher host-guest complex. Finally, we find that neglecting the
implicit nuclear coordinate dependence arising from the charge density
partitioning, as has been done in prior numerical treatments, leads to an
unacceptable error in the MBD forces, with relative errors of ~20% (on average)
that can extend well beyond 100%.
