  We present a comprehensive methodology to enable addition of van der Waals
(vdW) corrections to machine learning (ML) atomistic force fields. Using a
Gaussian approximation potential (GAP) [Bart\'ok et al., Phys. Rev. Lett. 104,
136403 (2010)] as baseline, we accurately machine learn a local model of atomic
polarizabilities based on Hirshfeld volume partitioning of the charge density
[Tkatchenko and Scheffler, Phys. Rev. Lett. 102, 073005 (2009)]. These
environment-dependent polarizabilities are then used to parametrize a screened
London-dispersion approximation to the vdW interactions. Our ML vdW model only
needs to learn the charge density partitioning implicitly, by learning the
reference Hirshfeld volumes from density functional theory (DFT). In practice,
we can predict accurate Hirshfeld volumes from the knowledge of the local
atomic environment (atomic positions) alone, making the model highly
computationally efficient. For additional efficiency, our ML model of atomic
polarizabilities reuses the same many-body atomic descriptors used for the
underlying GAP learning of bonded interatomic interactions. We also show how
the method enables straightforward computation of gradients of the observables,
even when these remain challenging for the reference method (e.g., calculating
gradients of the Hirshfeld volumes in DFT). Finally, we demonstrate the
approach by studying the phase diagram of C$_{60}$, where vdW effects are
important. The need for a highly accurate vdW-inclusive reactive force field is
highlighted by modeling the decomposition of the C$_{60}$ molecules taking
place at high pressures and temperatures.
