I am a newbie in the field, thus please excuse me if my problem might appear trivial.
I am trying to simulate (NPT) a set of rigid bodies (11 spheres arranged in a crescent-like manner) via WCA potential. Below I attach the input script. The thing that bothers me is that around the 4x10^6 timestep the whole system collectively starts to perform the translation move across the simulation box (see the video). Frankly, I have no clue what might be the cause. Any hints/advices are highly appreciated. initial_script.in (2.3 KB)
This is not as unusual as you might think.
The phenomenon is generally referred to as “flying ice cube syndrome”. This is because the thermostat cannot remove the center of mass drift, so to maintain the total temperature more and more kinetic energy will be transferred from vibrations to translation until the system freezes and at the same time accelerate a lot as a whole.
Its origin lies in the use of floating-point math and the fact that it summing floating-point numbers is inexact and that their resolution depends on their magnitude. Thus something that (theoretically) should sum up to (exactly) zero does not. Most of the time, this is not a problem. There will be some fluctuations around the origin of the cell, but it would mostly cancel out. But sometimes a center of mass drift will appear and then problems start, since there usually is nothing in a typical simulation that can dissipate it.
There are multiple contributing factors.
the way how the initial temperature is set (but you used the mom yes option, that will remove it)
the initial system geometry (you have an “asymmetric” geometry in your constituent molecules), an atomic system is less prone
the equilibration protocol: it is often a good idea to let a system equilibrate until the worst potential energy “hot spots” are well dissipated and then use “velocity all scale mom yes” or similar
using a variable cell time integration (because it converts all positions to fractional coordinates and conventional coordinates multiple times during a timestep and that adds numerical noise). this is a mechanism that becomes more dominant the larger your center of mass drift is
using a rigid body fix: there are additional issues with how time integration for rigid bodies work, specifically there are different errors and timestep requirements for translation and rotation. using a shorter timestep can help a lot
using a non-dissipative thermostat (worst are temp/rescale or temp/berendsen).
using GPU acceleration with mixed or single precision
using (too) aggressive compiler optimization
bugs in the code: try with the latest LAMMPS patch version (27 May 2021 as of this writing)
There are some cures, but since you are using a rigid fix, your options are limited, since you cannot use most of them (the rigid fix stores per rigid object velocities and positions). You can break down your simulation into multiple chunks and issue a velocity command in between to remove any residual center of mass drift. After that it is crucial that you repeat the fix rigid command (you do not have to unfix it, if you use the same fix) with the same settings. That will trigger initialization of the rigid body velocities (you will get a message in the screen and logfile output).
Unfortunately, the problem persists (even for twofold shorter timestep than initially).
Onward, I have included in the script:
fix mm1 all momentum 1000000 linear 1 1 1 rescale
and run the simulation for 500 rigid bodies (x11, ergo 5500 “atoms”). Below I attach the results for the components of system’s COM speed (v_x, v_y, v_z), which were monitored throughout the simulation via:
As one can see up to 1.2*10^7 Timestep there is no drift, afterwards “the ice cube syndrom starts to kick-in”. “Horizontal lines” are results of the fix mm1, however it is visible that after each “reset” the system builds up to the state from before the “reset” and “step-by-step” the whole system gradually drift away. Frankly, I am in a blind spot.
Fix momentum cannot have a lasting effect as the per-rigid-body velocities are cached in the rigid body fix. Only if you break your run into sections, reissue the fix rigid command, and remove the center of mass momentum with the velocity command, should you be able to properly remove it.
So “technically” speaking, I should follow the scheme:
(INITIALL STUFF)
fix rigid all rigid/npt/small molecule temp 1.0 1.0 1.0 iso 1.0 1.0 1.0
run 10000000
unfix rigid
velocity all zero linear
fix rigid all rigid/npt/small molecule temp 1.0 1.0 1.0 iso 1.0 1.0 1.0
run 10000000
etc.
since you are overwriting the same fix, there is no need to unfix it, just reinitialize it
I would remove the momentum more often so there is no significant drop in kinetic energy
with index variables their value can be overwritten from the command line with the -var flag
depending on how you run, you may need to replace SELF with the file name
variable ntotal index 70000000
variable nchunk index 1000000
variable nloop index $(round(v_ntotal/v_nchunk))
variable i loop ${nloop}
label loop
run ${nchunk} post no
velocity all zero linear
fix rigid all rigid/npt/small molecule temp 1.0 1.0 1.0 iso 1.0 1.0 1.0
next i
jump SELF loop
