NVE MD simulation of inert gas: Problem maintaining equilibrium

I have been trying to simulate a simple problem of taking about 100-1000 Ar molecules in a NVE (fixed vol, energy) system with equal speed but randomized velocity, and evolving them to obtain a Maxwellian speed distribution. I am choosing a Lennard Jones pair potential with no cutoff (despite being computationally expensive), periodic boundary conditions, and a varying integration time step to keep numerical error on energy to below 2%. I observe that my system reaches Maxwellian-like speed distribution (visually looks like a good Gaussian) from the initial rather uniform speed distribution. But very soon after some (about 10 in 100) particles accelerate to very large speeds even if my total energy is being held constant (to within 2% of initial energy) skewing the distribution to one side.

I have tried to vary the no. of particles, increased temperature to allow re-equilibration of what appears like a growing fluctuation but can't make this effect go away. I am puzzled because I am managing to keep total energy with a 2% bound. Has anyone else has observed such long-time behavior in simple MD simulations, and/or is it well-known behavior? I am using the velocity Verlet scheme. I am using a well-scaled system (i.e., not trying to code actual atomic dimensions) and using a reasonable well distributed initial condition.

• Just for record keeping, I posted this comment on the post in Phys.SE before it was closed: If your total energy remains the same but some particles go really fast it means others must be going to zero somewhere (otherwise your total energy computation is wrong). So I'd start looking at where those occur and see if something in the neighborhood is causing numerical stability issues (ie. check what your time step restriction needs to be and see if you are satisfying it at all times and locations). Jan 16, 2013 at 19:18
• What happens if you run your simulation with half the time step? Does the same thing happen at more or less the same (simulation) time? Jan 16, 2013 at 19:36
• I vary my time step at each iteration to keep total energy change to less than 2%, i.e., the code reduces time step by an order until the energy converges to within 2% of initial value (I have also tried 1%). There is noting noteworthy when this happens. The time step does not plummet or peak, the total energy fluctuation is smaller than other times, there are no unphysical velocities (e.g., negative), and oddly there are many plummeting time-steps at later stages but that does not seem to affect. Jan 16, 2013 at 19:49
• contd... Also, when I see the particle trajectories when this happens, I can see a few particles speeding up and then not ever interacting with others until the end. What surprises me is that this occurs repeatably but without apparent predictability. Well if it is just a bug I am surprised the distribution goes to reasonably Maxwellian to begin with. Jan 16, 2013 at 19:52
• @Sankaran : Can you please what is the potential you have used in this problem, I mean the LJ parameters. I have been trying to solve the same problem I have not been able to get it right !! Oct 8, 2015 at 5:55