Jump in atom position when it is over 21900

Hi Lammps community

I use the 10Feb15 version of Lammps to simulate diffusion of gas molecules with a simulation time of 50 ns and compute msd.

The problem i have is, msd numbers get crazy (example: 220 to 6*10^7 in 1 ps) at certain points. When i checked the trajectory, i have seen that this certain points are when a position of atom reaches 21900, because the next step is -21900 instead of positive.

What does this indicate? Is this a limitation?

Marcel

Hi Lammps community

I use the 10Feb15 version of Lammps to simulate diffusion of gas molecules
with a simulation time of 50 ns and compute msd.

The problem i have is, msd numbers get crazy (example: 220 to 6*10^7 in 1
ps) at certain points. When i checked the trajectory, i have seen that this
certain points are when a position of atom reaches 21900, because the next
step is -21900 instead of positive.

What does this indicate? Is this a limitation?

most likely a misunderstanding on your side, but it is impossible to
say anything with confidence from such little information.
can you provide the header section of your dump file and any dump
related commands in your input?

axel.

The dump command i use is;

group CO2 molecule > 81

compute msdc CO2 msd com yes

dump 2 CO2 custom 10000 co2.lammpstrj id type xu yu zu vx vy vz

and i am attaching part of my dump file which is at the moment when the switches begin.

co2.lammpstrj (3.29 MB)

The dump command i use is;

group CO2 molecule > 81
compute msdc CO2 msd com yes
dump 2 CO2 custom 10000 co2.lammpstrj id type xu yu zu vx vy vz

and i am attaching part of my dump file which is at the moment when the
switches begin.

ok, what is happening is that your atom/molecules have passed through
the box so many times that the image flags wrap around.
image flags are normally stored together in single 32-bit integer,
leaving space for a signed 10-bit integer which can store numbers from
-512 to 511.
dividing your x dimension by the size of the box shows that your
system has reached that point.

what is suspicious is that your coordinates are not quite as far
scattered as the image flags are needed, which hints at the fact that
you may have a significant drift of your center of mass and thus may
be subject to the so-called flying ice cube syndrome, i.e. that your
system is actually behaving as if it was much colder than it is, since
the temperature compute does not consider the bias resulting from your
center of mass drift.

when i analyze the velocity information of your trajectory segment, i
find an average drift of (0.0111,-0.064,-0.014) which is consistent
with the x dimension wrapping around first.

to get a meaningful output, you will have to find the origin of this
drift and remove it.

axel.

If you compile LAMMPS with -DLAMMPS_BIGBIG

the image flags are 20 bits, instead of 10, So you

can wrap-around up to 1M times. See Section 2,2

of the manual.

Steve

Axel, thank you for your quick help.

I have three polymers, all having this center of mass drift which we could not find a solution and as our COM substracted MSD results showed accuracy with literature, we ignored it for now.

MD and MC results of all polymers showed good accuracy with experimental results (Tg, sorption coefficients of several gases, Permeabilities) and as we already used Nose-Hoover instead of Berendsen, we could not find any solution and went on.

Steve, thank you for your help also. Is that an option i need to use while making Lammps on command line or is it something i need to input into makefile?

Marcel

If you compile LAMMPS with -DLAMMPS_BIGBIG
the image flags are 20 bits, instead of 10, So you
can wrap-around up to 1M times. See Section 2,2
of the manual.

yeah, but that doesn't solve the drift problem, only postpones it.

...and unless marcel uses a dissipative thermostat (unlikely or
otherwise the drift would dissipate as well) the drift will only grow
over time and thus make the error in the computed temperature larger
and the MSD worse and worse.

axel.

Axel, thank you for your quick help.

I have three polymers, all having this center of mass drift which we could
not find a solution and as our COM substracted MSD results showed accuracy
with literature, we ignored it for now.

ignorance is rarely a good idea. :wink:

many times the COM drift originates from the initial temperature
assignment or is an artifact of equilibration. so using momentum
removal option of the velocity command often helps or at least reduces
the issue. best would be to do it twice: right at the beginning and
then after equilibration and then just rescale the velocities to
regain the equilibrium temperature.

other potential reasons are too aggressive a time step, too short a
cutoff, too sloppy kspace convergence and similar or a systematic
inhomogeniety of your system. the impact of most of those can be
reduced by careful testing of parameters. the last is pathological.

if there is still a significant drift forming, after all other options
have been properly exploited, you can use fix momentum. normally, the
this would just remove kinetic energy from your system, and thus is
undesirable, if the amount of energy is too large. i have a modified
version of fix momentum in my LAMMPS-ICMS branch that has an added
option to redistribute the COM drift kinetic energy to all other DOFs.
with fix momentum, the trick is to have it run infrequently, but not
so infrequent that the energy to be redistributed becomes significant.

HTH,

axel.

p.s.: steves suggestion requires a full recompiled of LAMMPS. the flag
he mentions is a preprocessor define.

Axel, thank you for your quick help.

I have three polymers, all having this center of mass drift which we could
not find a solution and as our COM substracted MSD results showed accuracy
with literature, we ignored it for now.

MD and MC results of all polymers showed good accuracy with experimental
results (Tg, sorption coefficients of several gases, Permeabilities) and as
we already used Nose-Hoover instead of Berendsen, we could not find any
solution and went on.

please note that not all properties are equally affected by such an
error, and that also data in the published literature can be off.
especially more recently with more and more of the simulation work
being done by less well trained and less experienced people. subtle
issues like a moderate COM drift cannot be spotted when reviewing
simulation data in a publication. if this is a systemic problem, then
perhaps the authors of the publication would have to deal with the
same problem.

axel.

We actually tried some of the things you recommended and couldn’t make any difference.

Before we begin our 50 ps NVT run, we tried to randomly distribute atom velocities to remove effects from prior runs.
We tried to use “fix 1 all momentum 1 linear 1 1 1” again in the beginning of the simulation.

Both did not help us. I will try everything you proposed, again i thank you for your quick and detailed help. I will update when i solve this problem.

Marcel

We actually tried some of the things you recommended and couldn't make any
difference.

Before we begin our 50 ps NVT run, we tried to randomly distribute atom
velocities to remove effects from prior runs.

that doesn't make sense, just remove the COM velocity. by random
reseeding velocities you increase the likeliness to have a COM drift
right from the beginning.
also, you should check what kind of value your image flags have before
you start your production and reset them.

We tried to use "fix 1 all momentum 1 linear 1 1 1" again in the beginning
of the simulation.

applying fix momentum for every step is a bad idea. it will confuse
the NH-thermostat and disallow necessary fluctuations. if you have a
systemic COM drift issue, you definitely want to use my modified fix
momentum to put less strain on the NH-thermostat to keep the
temperature where you want it.

axel.

Yes, we are resetting the image flags at the beginning. I will try removing COM velocity and look if the problem persists.

Then i will try the LAMMPS-ICMS from your website. With this much help, i will need to cite you more than one if this work somehow gets to publishment stage :slight_smile:

Again, thanks a lot. Will let you know.

Marcel

Yes, we are resetting the image flags at the beginning. I will try removing
COM velocity and look if the problem persists.

Then i will try the LAMMPS-ICMS from your website. With this much help, i

the change to fix momentum is independent, so you should be good by
just downloading the fix_momentum.h and fix_momentum.cpp and drop it
into your existing LAMMPS copy. LAMMPS-ICMS has currently several
other changes that are not in the mainline code, so if you are running
production, you probably don't want to risk importing a difference (i
do test my changes, but mistakes can happen to everybody...).

axel.

One last question is, are there any problem in using this new fix momentum with GPU package?

One last question is, are there any problem in using this new fix momentum
with GPU package?

not at all. the GPU package is moving updated data back and forth
during each step and thus has no practically compatibility issue with
fixes or computes (except for those requiring access to neighbor
lists).

since you mention the GPU package, are you running in single
precision? that may have a significant impact.

which reminds me of another thing: i noticed that your origin is in
the corner of the box. try moving your system so that the origin is in
the center. due to the nature of floating point math, you are
importing a bias from the different resolution of numbers. centering
the box around the origin will distribute the error more evenly. this
becomes particularly visible when using single precision math.

axel.

If you guys allow me I will scoop in with a small observation that is not so much aim at you but at the future newbies reading this post:

The ultimate cause that allows formation of such drifts (leading to the many times reported flying ice-cube problem & Co) is simply the fact that the unit cell representation chosen by lammps eliminates the box rotation eigenvalues but not the translational ones, i.e, atoms in reduced coordinates can be translated within the simulation box while keeping the system invariant. This is something one has to deal with when transition state searches are performed in the enthalpy landscape at finite pressures. The most general solution would be to remove these modes from the dynamics but of course this is something that would only make sense as an optional feature given that a lot of people use lammps in simulations where the translational modes are indeed wanted.

Carlos

I have been using a rented cluster for my runs before and GPU was not enabled then. I am using my Tesla 40K GPU for smaller runs, I use it with mixed precision.

As an update, removing linear momentum in the beginning of simulation seems to solve the problem for now. No significant COM shift is seen in 50 ns.

Thanks for the help.