  The computational bottleneck of molecular dynamics is the pairwise additive
long-range interactions between particles. The random batch Ewald (RBE) method
provides a highly efficient and superscalable solver for long-range
interactions, but the stochastic nature of this algorithm leads to unphysical
self-heating effect during the simulation. We propose an energy stable scheme
(ESS) for particle systems by employing a Berendsen-type energy bath. The
scheme removes the notorious energy drift which exists due to the force error
even when a symplectic integrator is employed. Combining the RBE and the ESS,
the new method provides a perfect solution of the computational bottleneck of
molecular dynamics at the microcanonical ensemble. Numerical results for
primitive electrolyte and all-atom pure water systems demonstrate the
attractive performance of the algorithm including its dramatically high
accuracy, linear complexity and overcoming the energy drift for long-time
simulations.
