Atom drift (perhaps due to the flying ice cube effect)

Efrem_Braun1 · June 29, 2016, 9:31pm

Hello,

I’m running some simulations of benzene in a flexible porous framework. I have two groups of atoms, those of the adsorbate and those of the framework. I use separate time integration fixes for these two groups (since the benzene molecules are rigid), and then I set velocities (zeroing out linear momentum) and run.

The atoms then begin to drift continuously in a particular direction. I’m trying to figure out why this drift occurs, what I need to do to fix it, and whether I need to remove the velocity bias before thermostatting to the temperature I want.

The thermostat is what’s causing the drift. When I ran the simulation in the NVE ensemble, no drift occurred. (Note that energy in the NVE ensemble was conserved, so the timestep is not too large, and drift occurs in the NVT ensemble even when I lowered the timestep from 0.5 to 0.1 fs, in response to the suggestion of http://lammps.sandia.gov/threads/msg34724.html.)

I suspect that the problem is due to the flying ice cube effect, though I’m not completely sure. My working through the solutions proposed in http://lammps.sandia.gov/threads/msg53543.html:

Timestep: tested (see above), not the issue.
Cutoff: The cutoff of 12.0 that I’m using is almost as large as it can be with a 25.8 Angstrom box, and it’s also more than 3X bigger than the greatest LJ sigma, so I can’t imagine this is the issue.
Kspace convergence: I turned off long-range electrostatics, and the drift still occurs.

Note that when I have just the benzene atoms in the system, or when I have just framework atoms in the system, no drift occurs. It’s only when both are there that I see drift. Why would the flying ice cube effect only be present if the full system is in place, but not when each half is there separately?

Finally, if it is the flying ice cube effect, I’m not quite sure how to fix it since my simulation has rigid bodies. Per http://lammps.sandia.gov/threads/msg56336.html, “fix momentum” would be the thing to use, but per http://lammps.sandia.gov/threads/msg34731.html, “fix momentum” doesn’t work when rigid bodies are present. I’m guessing that “set velocity zero linear” has the same problem. Any way to zero out the linear momentum with rigid molecules?

A minimal working example is below, with the data file attached. Any pointers would be greatly appreciated.

Efrem Braun

#settings
units real
atom_style full
neigh_modify delay 0 every 1
pair_style lj/cut/coul/long 12.0
pair_modify mix arithmetic shift yes
special_bonds lj/coul 0.0 0.0 1.0
kspace_style pppm 1.0e-6
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style cvff

#initialize atomic positions and FF
read_data “data_init.data”

#settings
group Benzene type 1 5
group MOF type 2 3 4 6 7 8 9
timestep 0.5
thermo 200
dump 1 all custom 200 traj.dump id mol type x y z ix iy iz xu yu zu

#run
fix 1 Benzene rigid/nvt molecule temp 300 300 100.0
fix 2 MOF nvt temp 300 300 100.0
run 0
velocity all create 300 38564
run 400000

data_init.data (348 KB)

akohlmey · June 29, 2016, 9:41pm

Hello,

I'm running some simulations of benzene in a flexible porous framework. I
have two groups of atoms, those of the adsorbate and those of the framework.
I use separate time integration fixes for these two groups (since the
benzene molecules are rigid), and then I set velocities (zeroing out linear
momentum) and run.

The atoms then begin to drift continuously in a particular direction. I'm
trying to figure out why this drift occurs, what I need to do to fix it, and
whether I need to remove the velocity bias before thermostatting to the
temperature I want.

The thermostat is what's causing the drift. When I ran the simulation in the
NVE ensemble, no drift occurred. (Note that energy in the NVE ensemble was
conserved, so the timestep is not too large, and drift occurs in the NVT
ensemble even when I lowered the timestep from 0.5 to 0.1 fs, in response to
the suggestion of http://lammps.sandia.gov/threads/msg34724.html.)

I suspect that the problem is due to the flying ice cube effect, though I'm
not completely sure. My working through the solutions proposed in
http://lammps.sandia.gov/threads/msg53543.html:
1. Timestep: tested (see above), not the issue.
2. Cutoff: The cutoff of 12.0 that I'm using is almost as large as it can be
with a 25.8 Angstrom box, and it's also more than 3X bigger than the
greatest LJ sigma, so I can't imagine this is the issue.
3. Kspace convergence: I turned off long-range electrostatics, and the drift
still occurs.

Note that when I have just the benzene atoms in the system, or when I have
just framework atoms in the system, no drift occurs. It's only when both are
there that I see drift. Why would the flying ice cube effect only be present
if the full system is in place, but not when each half is there separately?

Finally, if it is the flying ice cube effect, I'm not quite sure how to fix
it since my simulation has rigid bodies. Per
http://lammps.sandia.gov/threads/msg56336.html, "fix momentum" would be the
thing to use, but per http://lammps.sandia.gov/threads/msg34731.html, "fix
momentum" doesn't work when rigid bodies are present. I'm guessing that "set
velocity zero linear" has the same problem. Any way to zero out the linear
momentum with rigid molecules?

A minimal working example is below, with the data file attached. Any
pointers would be greatly appreciated.

a simple solution to suppress a drift of your framework would be to
attach a weak position restraint to it using fix spring.

axel.

Efrem_Braun1 · June 29, 2016, 11:11pm

Thanks Axel, that should solve the problem.

That said, I’d still be interested in knowing more about the problem to make sure there’s not some pathological error that this patch hides but doesn’t fix.

One reason I’m worried about this is because I ran the simulation 10 times with 10 different random number seeds for the “velocity create” lines, and each time the atoms drifted in the same direction. See attached rendering of the final snapshot; the turquoise atoms are the benzenes in the initial unit cell (they were wrapped by LAMMPS due to them being rigid molecules, but in truth they drifted with the framework atoms), and you can see that all frameworks drifted in the negative x and y directions (and the same z direction too).

If the drifting is due to the flying ice cube effect, they should be going in random directions, no? The linear momentum would build up in the direction of a small random deviation from 0 linear momentum that happened to occur at the right time.

vmdscene.tga (768 KB)

Efrem_Braun1 · June 30, 2016, 7:14am

Thanks Axel, that should solve the problem.

That said, I’d still be interested in knowing more about the problem to make sure there’s not some pathological error that this patch hides but doesn’t fix.

One reason I’m worried about this is because I ran the simulation 10 times with 10 different random number seeds for the “velocity create” lines, and each time the atoms drifted in the same direction. See attached rendering of the final snapshot; the turquoise atoms are the benzenes in the initial unit cell (they were wrapped by LAMMPS due to them being rigid molecules, but in truth they drifted with the framework atoms), and you can see that all frameworks drifted in the negative x and y directions (and the same z direction too).

If the drifting is due to the flying ice cube effect, they should be going in random directions, no? The linear momentum would build up in the direction of a small random deviation from 0 linear momentum that happened to occur at the right time.

akohlmey · June 30, 2016, 10:51am

Thanks Axel, that should solve the problem.

That said, I'd still be interested in knowing more about the problem to make
sure there's not some pathological error that this patch hides but doesn't
fix.

One reason I'm worried about this is because I ran the simulation 10 times
with 10 different random number seeds for the "velocity create" lines, and
each time the atoms drifted in the same direction. See attached rendering of
the final snapshot; the turquoise atoms are the benzenes in the initial unit
cell (they were wrapped by LAMMPS due to them being rigid molecules, but in
truth they drifted with the framework atoms), and you can see that all
frameworks drifted in the negative x and y directions (and the same z
direction too).

If the drifting is due to the flying ice cube effect, they should be going
in random directions, no? The linear momentum would build up in the
direction of a small random deviation from 0 linear momentum that happened
to occur at the right time.

due to having a part interconnected (across PBC), part free flowing
system, you may also have momentum buildup for each subsystem due to
their initial geometry. the total system momentum will then be
preserved. that is not what is usually referred to by the flying
icecube syndrom, where your whole system will pick up momentum and
start drifting, where the thermostat will then lead to the relative
temperature of the whole system going down until it suddenly freezes
and all kinetic energy will be in global translation.

what you are observing, has sometimes been referred to as a "system
creep". it can be minimized by using the mom yes flag to the velocity
command, a dissipative thermostat - especially during equilibration -
and then removing the center of mass velocity selectively for both the
stationary and fluid parts during equlibration.

at any rate, using fix spring on the stationary framework is the right
approach here, and that will make any use of fix momentum or other
tricks during production run (using fix nvt) unnecessary. the force
constant required should be quite small, so the impact minimal. if the
fluid part will "migrate" in a particular direction, that should not
be a big concern, as the net speed of this process will be extremely
slow and thus should not affect results negatively.

apropos speed. in case you are running in parallel, you may want to
consider using fix rigid/small/nvt for better parallel scaling.

axel.

Efrem_Braun1 · June 30, 2016, 1:32pm

Axel, thanks again for the response. I’m going to have to acknowledge you in my dissertation.

Unfortunately, the total linear momentum of the system is not being conserved. Both the adsorbates and the framework have their centers of mass drifting in the same direction (this after I zero out the system’s linear momentum in the beginning of the run with ‘velocity all create 300 38564’ which uses the default keyword ‘mom yes’). If they were drifting in opposite directions and the COM of the overall system stayed in place, then I agree that I wouldn’t need to be worried about a more pathological problem. This indicates that the problem is not “system creep”, but something else, no?

On a side note, although I’ve been reading in the mail list archives that the flying ice cube effect is blamed for drift of the system COM (e.g. http://lammps.sandia.gov/threads/msg53543.html), I don’t think that the flying ice cube effect occurs when the Nose-Hoover thermostat is used since it doesn’t use a velocity rescaling algorithm. I tested this by recreating the model diatomic system given in the initial publication of the flying ice cube effect (doi.org/10.1002/(SICI)1096-987X(199805)19:7<726::AID-JCC4>3.0.CO;2-S). I found that the bond energy decays to 0 while the COM kinetic energy increases when using a Berendsen thermostat, but everything looks fine using the default Nose-Hoover thermostat of ‘fix nvt.’ See attached for my input files and plots of the results. Am I misunderstanding something?

system.data (335 Bytes)

system.in.old (1.18 KB)

athomps · June 30, 2016, 3:55pm

Efrem,

I suspect that both the Berendsen and NVT simulations are invalid, for two reasons:

Average COM velocity is not slowly growing but quickly rises to a large value (of order 100 m/s ~ atomic thermal speed)
The PE distributions are very different than NVE.
For NVT, the fluctuations in COM velocity are correlated with fluctuations in PE

None of these things should happen in simulations that are even halfways decent.

I am guessing that in the case of NVT, you are getting some kind of interaction between the two thermostats:

fix 1 Benzene rigid/nvt molecule temp 300 300 100.0
fix 2 MOF nvt temp 300 300 100.0

One thing to try would be after some short equilibration, zeroing the COM of both groups separately. I suspect the culprit is rigid/nvt, since this is less widely tested. If you replaced that with straight nvt, the problem might go away.

Aidan