Molecular dynamics simulations are usually O(N log(N)), because that's the time complexity of the Fast Fourier Transform, which allows computation of the N^2 charge interactions in less than N^2 time (e.g. Ewald Summation https://en.wikipedia.org/wiki/Ewald_summation).

For short-range interactions like Van der Waals forces, MD is O(N) like you said, because you only have to evaluate a finite radius around each atom. But charge interactions, decaying as 1/r, don't converge in any finite radius.