Many-body electronic responses such as dispersion and polarization (at and beyond dipole order) present fundamental challenges in the simulation of materials at the molecular scale. To address these, an emerging strategy employing embedded quantum oscillators as a coarse-grained representation of such responses has been effective in predicting material properties with remarkable accuracy. However, applications have so far been limited to relatively small system sizes. Here we introduce strategies enabling efficient implementation of this framework on high-performance, heterogeneous CPU–GPU (with multiple Graphic Processing Units) computing platforms thereby increasing significantly the scale of accessible problems. Physical properties are reported for a benchmark system of 104 water molecules — to our knowledge, the largest system yet simulated using molecular dynamics with electronically-derived forces.