DL_FFLUX is a force field based on quantum chemical topology that can perform molecular dynamics for flexible molecules endowed with polarisable atomic multipole moments (up to hexadecapole). Using the machine learning method kriging (aka Gaussian Process Regression) DL_FFLUX has access to atomic properties (energy, charge, dipole moment, …) with quantum mechanical accuracy. Newly optimised and parallelised using domain-decomposition MPI, DL_FFLUX is now able to deliver this rigorous methodology at scale while still in reasonable time frames. DL_FFLUX is delivered as an add-on to the widely distributed molecular dynamics code DL_POLY 4.08. For the systems studied here (10 3 -10 5 atoms), DL_FFLUX is shown to add minimal computational cost to the standard DL_POLY package. In fact, the optimisation of the electrostatics in DL_FFLUX means that, when high rank multipole moments are enabled, DL_FFLUX is up to 1.25x faster than standard DL_POLY. The parallel DL_FFLUX preserves the quality of the scaling of the MPI implementation in standard DL_POLY. For the first time it is feasible to use the full capability of DL_FFLUX to study systems that are large enough to be of real world interest. For example, a fully flexible, high-rank polarised (up to and including quadrupole moments) 1 ns simulation of a system of 10,125 atoms (3,375 water molecules) takes 30 hours (wall time) on 18 cores.