Hong Nguyen and I are seeing a noticeable center-of-mass drift in long NPT simulations of bead-spring polymer melts, and wondering how to suppress it.
We’re using the 9Sep2014 version of LAMMPS.
The model we’re simulating is similar to the Kremer-Grest model: LJ pair interactions, harmonic bonds, and cosine angles. Our systems are smallish: 500 chains of 25 monomers each --> total number of beads = 12500.
The relevant commands (those controlling time-integration) we’re using are:
fix 1 all npt temp 0.62 0.62 100 iso 0 0 1000
timestep .005
run 300000000
Over the course of the 300 million timestep (1.5 million Lennard-Jones time unit) run, the system picks up a net momentum that can be clearly seen in VMD movies. This effect is reproducible - we also see it in other long simulations with different initial states and different T.
The question I have is: is there a way to suppress this drift velocity using e.g. flags in fix_npt? Looking at the online docs, there doesn’t seem to be one. I’d thought there used to be a way to do this in LAMMPS, but I may be misremembering. (In the past I’ve mostly used Langevin thermostats where one simply adds a flag to zero the net momentum.)
It is a know artifact of long time integration with thermal control, called the “flying ice cube” artifact.
I found that the use of the “fix recenter” only hide the problem. Check the RDF before and after the artifact happens. If the peaks of the RDF become suddenly sharper without any reason, it is the “flying ice cube” artifact.
In this case, the only solution i found is to reset velocity periodicaly. Every 3000000 steps for example.
It is a know artifact of long time integration with thermal control,
called the “flying ice cube” artifact.
I found that the use of the “fix recenter” only hide the problem. Check
the RDF before and after the artifact happens. If the peaks of the RDF
this is a very important point. fix recenter doesn't remove the center of
mass momentum. fix momentum does, and with the recently added "rescale"
flag, you can also do it in a fashion that conserves the energy.
I agree, fix_momentum (with rescale) does sound like the better choice (since it operates on momenta rather than positions). We’ll definitely check into this…
Hi again, guys. We’re still seeing rather finicky behavior in our simulations, and so have some followup questions. I understand why one would want to use the rescale flag. However…
If one zeros out the system’s center-of-mass momentum (COMm) at the beginning of the simulation, and there are no body forces on the system, than any subsequent changes in COM are “errors” of one type or another. These can be numerical ones arising from the finite integration timestep or physical ones arising from the design of the thermostatting algorithm.
Since fix_momentum WITHOUT rescaling just enforces Newton’s laws, there shouldn’t be any problem with applying the fix arbitrarily often. Momentum is strictly conserved only if the fix is applied every single timestep, so from a pure physics POV (ignoring the communication-time issues arising from the MPI_Allreduces and so on), it makes sense to apply it every timestep.
Fix_momentum WITH rescaling is a different story because it fixes kinetic energy. This affects the system’s dynamics and thermodynamics. Applying the fix every timestep would now be equivalent (afaict) to applying an isokinetic thermostat. We don’t want that.
The code in fix_momentum.cpp rescales system kinetic energy to its value in before the COMm is removed. It’s “time-local” in the sense fix_momentum.cpp is just using comparing ekin_old and ekin_new, which are evaluated at different points during the same integration timestep (specifically, before and after COMm is removed).
This is where I get confused. Points (3-4) together seem to imply that one shouldn’t apply the with-rescaling version of fix_momentum every timestep UNLESS one wants isokinetic systems. This makes sense by itself, but then what about point (2) and the points raised by Xavier and Axel below? Is it really a case that eliminating the flying ice cube effect only makes physical sense if you apply energy conservation infrequently? And if so, wouldn’t it make more sense for the fix to “remember” the kinetic energy of the system during the LAST timestep it was applied, instead of or in addition to the kinetic energy earlier in the same timestep? In other words. if the fix is applied every N steps, wouldn’t one want the “ekin_old” variable in fix_momentum.cpp to be the kinetic energy at time (t-Ndt) rather than the “present” time t? It seems the (t-Ndt) option would do a better job of maintaing kinetic energy over long times.
And one can’t resolve this issue by applying fix_momentum infrequently, because then COMm drift will occur.
Consequently we’re having a hard time deciding what’s a reasonable application frequency for fix_momentum with the rescaling option. Xavier mentions every applying it 3 million steps, but of course this is a system-dependent number. We’re seeing significant N-dependence in measures of dynamical relaxation that are expected to be sensitive to COMm drift and kinetic energy fluctuations. Of course such issues aren’t unusual in MD, but in this case I don’t have a good sense of the right/accepted answer.
The NH thermostat and barostat implementation in LAMMPS have some nice
mathematical properties that I believe include conserving COM velocity
(at least for the special case of zero initial COM velocity). But
LAMMPS is not mathematics, it is code executed on a finite-precision
computer, so after 300 million timesteps, I am not surprised that you
are seeing some drift. The only question is, what to do about it. I
can think of three choices:
1. Suppress it every timestep (fix momentum command).
2. Suppress it every large number of timesteps (using set velocity
zero linear command).
3. Allow drift, but use temp/com with npt fix to avoid flying ice cube issues.
Option 1 is probably the best, especially if the effect is very small,
as I suspect. I don't think it should matter whether you use the
rescale option, assuming the momentum adjustment per timestep is much
smaller than the adjustments due to the thermostat and barostat.
Note: fix momentum with rescale will *not* result in an isokinetic
trajectory. The kinetic energy is conserved across the fix momentum
operation at the end of the timestep, not across the entire timestep.
No matter what route you choose, it is a good idea to construct a
particularly badly drifting system to test whether everything is being
done as you expect in your script. You don't want to have to run 300
million timesteps on a large system to find out that something was not
right.
Thanks, Aidan! I see the source of my confusion now. For some reason I’d assumed the rescaling produced a isokinetic trajectory. But of course it only means that FIX_MOMENTUM doesn’t change the kinetic energy — other fixes still change KE in the same way they would in fix_momentum’s absence. This makes sense, of course. I’ve played with LAMMPS enough over the years that I should have realized this beforehand. Oops!
We’ll try applying fix_momentum/rescale every step (+ various rapid frequencies) to a badly drifting system and see what the differences are. If we find anything interesting, I’ll report back.
option 4. use a weak tether on the total center of mass (of the mobile
part) with fix spring (will fail, if the number of particles in the
fix group varies)
option 5. use fix momentum frequently, but not in every step. i would
expect that calling fix momentum in every step will make matters
worse, since you enforce truncation to the floating point resolution,
and thus you create noise. if you apply fix momentum every 10, 20, 50,
100 steps, you allow for some cancellation of the errors due to the
coarseness of floating point numbers and non-associative summation,
but don't allow to build a significant momentum.
Ah. Thanks, Axel, that makes sense. Option 5 seems like a good way to go.
BTW, I replied to Aidan’s message using my USF email instead of this one, and got an automated reply from lammps-users-bounces. So I’m not sure my reply to Aidan was posted to lammps-users. If not, someone please let me know and I’ll resend it :).
NOTE: This implementation has been shown to conserve linear momentum up to machine precision under NVT dynamics. Under NPT dynamics, for a system with zero initial total linear momentum, the total momentum fluctuates close to zero. It may occasionally undergo brief excursions to non-negligible values, before returning close to zero. Over long simulations, this has the effect of causing the center-of-mass to undergo a slow random walk. This can be mitigated by resetting the momentum at infrequent intervals using the fix momentum command.
