5

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.

Sankaran
  • 215
  • 1
  • 5
  • 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). – tpg2114 Jan 16 '13 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? – Pedro Jan 16 '13 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. – Sankaran Jan 16 '13 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. – Sankaran Jan 16 '13 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 !! – user35952 Oct 08 '15 at 05:55

2 Answers2

4

Did you check the centre-of-mass velocity of your particles, i.e. did you set it to zero at the beginning of your simulation? The centre-of-mass velocity should be preserved throughout the simulation, and should be zero to get correct temperatures (see Flying Ice Cube).

Pedro
  • 9,573
  • 1
  • 36
  • 45
  • Interesting point! I did set my center of mass velocity to be zero (machine precision 1e-17) initially (was not aware of flying ice cube though, thanks!). It seems to jump to 1e-9 around the same place, although not quite at the same time step. Although non-zero and a several orders of mag jump, 1e-9 is still pretty small for velocity to jump from a 3 to 3e4! – Sankaran Jan 16 '13 at 21:30
  • Update: In many runs I could reproduce the same problem with no change in center of mass velocity (1e-17 to 1e-16). If I filter out the rouge particles (about 2-3 usually and always less than 10 in 100) the remaining distribution in Gaussian looking. Just that these sit way higher in velocity. – Sankaran Jan 16 '13 at 21:55
  • Then it is very possible that your time step prediction is wrong. What happens when you use a sufficiently small, constant time step? – Pedro Jan 16 '13 at 22:21
  • 1
    Yes you are right. I was so focussed on getting the energy correct that I was coarsening the time step in the parts when energy was being conserved with no problem. My guess is that when the situation was not strongly attracting or repelling, the coarse time step allowed some particles to step over potential attraction and therefore interaction. Now of course I don't get the problem, but of course not able to proceed to much larger times. I think these particles enter the highly repulsive regions and then the time step required becomes very small! Thanks again. – Sankaran Jan 17 '13 at 19:23
0

I think you should capture these few particles behind a door in a separate compartment. There would doubtlessly be some energy that can be gained from them.

I'm sure Maxwell would be proud of such a demonic scheme if you succeeded!

Wolfgang Bangerth
  • 55,373
  • 59
  • 119