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ók 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 Londondispersion 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.