A hybrid neural network (NN) and Lennard-Jones (LJ) model driver for multilayer graphene systems. The NN term models short-range intralayer and orbital overlap interactions and the theoretically-motivated LJ term models long-range dispersion. The potential is trained against a large dataset of monolayer graphene, bilayer graphene, and graphite configurations obtained from ab initio total-energy calculations based on density functional theory (DFT). The potential provides accurate energy and forces for both intralayer and interlayer interactions, correctly reproducing DFT results for structural, energetic, and elastic properties such as the equilibrium layer spacing, interlayer binding energy, elastic moduli, and phonon dispersions to which it was not fit.