Abstract In this work, we study the damage in crystalline molybdenum material samples due to neutron bombardment in a primary knock-on atom (PKA) range of 0.5–10 keV at room temperature. We perform classical molecular dynamics (MD) simulations using a previously derived machine learning (ML) interatomic potential based on the Gaussian approximation potential (GAP) framework. We utilize a recently developed software workflow for fingerprinting and visualizing defects in damaged crystal structures to analyze the Mo samples with respect to the formation of point defects during and after a collision cascade. As a benchmark, we report results for the total number of Frenkel pairs (a self-interstitial atom and a single vacancy) formed and atom displacements as a function of the PKA energy. A comparison to results obtained using an embedded atom method (EAM) potential is presented to discuss the advantages and limits of the MD simulations utilizing ML-based potentials. The formation of Frenkel pairs follows a sublinear scaling law as ξ b where b is a fitting parameter and ξ = E PKA/E 0 with E 0 as a scaling factor. We found that the b = 0.54 for the GAP MD results and b = 0.667 for the EAM simulations. Although the average number of total defects is similar for both methods, the MD results show different atomic geometries for complex point defects, where the formation of crowdions by the GAP potential is closer to the DFT-based expectation. Finally, ion beam mixing results for GAP MD simulations are in a good agreement with experimental mixing efficiency data. This indicates that the modeling of atom relocation in cascades by machine learned potentials is suited to interpret the corresponding experimental findings.