limit on absolute value of atom coordinate? (fix npt, boundary p p p)

Dear LAMMPS users,

I was running a simulation of a polymer gel, and I made a silly mistake, forgetting to zero the velocity of the center of mass of the system during a preliminary run.

The result is that during my production run the whole system was moving at quite a high speed. This wouldn’t normally be a problem since the center of mass velocity can be always subtracted when analyzing the coordinates.

However, I noticed something unexpected (to me): every time the coordinate of an atom becomes larger than approximately 500 times the box size, its sign gets changed, i.e., it’s like some ulterior PBC are applied to the system (see attached plot, showing the coordinates of an atom, rescaled by dividing for the box length, vs the number of steps. Notice that I’m running an NPT simulation using “fix nvt”, therefore the box length is not constant).

I report below my input script. If some of the commands I am using implicitly enforces this limitation on the atom coordinates, I am not aware of it.

Could someone help me understand what is causing this?

I am using the version LAMMPS 16 Mar 2018.

Thank you very much,

Valerio

INPUT SCRIPT:

variable sr_two equal (2^(1/6.)) variable mytemp equal 1.0 variable mypress equal 0.0 variable tnpt equal 5e8 variable thermodump equal 1000 variable trajdump equal (v_tnpt/200)
variable trestart equal (v_tnpt/5) variable tstep equal 0.003 variable Tdamp equal 100 variable pdamp equal 1000 variable fname string data variable thermofile string thermo.dat variable simname string npt_nose_T{mytemp}P${mypress}
variable log_times_file file log_times
${tnpt}_200_5.times
variable logdump equal “next(log_times_file)”

units lj
boundary p p p
atom_style bond

read_restart ${fname}
reset_timestep 0

group mono type 1
group crosslink type 2

mass 1 1.0
mass 2 1.0

pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 ${sr_two}
pair_modify shift yes

special_bonds fene
bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

neighbor 0.4 bin
neigh_modify every 10 delay 0 check yes

variable step equal step
variable temp equal temp
variable etotal equal etotal
variable ke equal ke
variable pe equal pe
variable ebond equal ebond
variable vol equal vol
variable press equal press
variable pxx equal pxx
variable pxy equal pxy
variable pxz equal pxz
variable pyy equal pyy
variable pyz equal pyz
variable pzz equal pzz

thermo_style custom step ke pe temp press vol
thermo ${thermodump}

timestep ${tstep}

fix 0 all print {thermodump} "{step} {ke} {pe} {ebond} {temp} {press} {vol} {pxx} {pxy} {pxz} {pyy} {pyz} {pzz}" file ${thermofile} screen no title “#1:step 2:ke 3:pe 4:ebond 5:temp 6:press 7:vol 8:pxx 9:pxy 10:pxz 11:pyy 12:pyz 13:pzz”

fix 1 all npt temp {mytemp} {mytemp} (v_tstep*v_Tdamp) iso {mypress} {mypress} (v_tstep*v_pdamp) fixedpoint 0.0 0.0 0.0

dump 1 all custom {trajdump} dump/dump.{simname}.*.lammpstrj id type xu yu zu
dump_modify 1 sort id
dump_modify 1 format line “%d d .8f .8f .8f”
dump_modify 1 every v_logdump

dump 2 all custom {trajdump} dumplin/dump.{simname}.*.lammpstrj id type xu yu zu
dump_modify 2 sort id
dump_modify 2 pad 10
dump_modify 2 format line “%d d .8f .8f .8f”

restart {trestart} restart/{simname}.*.restart

run ${tnpt}

unfix 0
unfix 1

write_data data.{simname} write_restart restart/restart.{simname}

moving_com.pdf (32.7 KB)

Dear LAMMPS users,

I was running a simulation of a polymer gel, and I made a silly mistake, forgetting to zero the velocity of the center of mass of the system during a preliminary run.
The result is that during my production run the whole system was moving at quite a high speed. This wouldn't normally be a problem since the center of mass velocity can be always subtracted when analyzing the coordinates.

actually, it is a problem if you are using a thermostat, since the
actual temperature will be biased by the high center of mass velocity,
and thus you are seeing what is often referred to as the "flying
icecube syndrome". so the effective temperature of your system will be
reduced by the kinetic energy in the of the center of mass motion.

However, I noticed something unexpected (to me): every time the coordinate of an atom becomes larger than approximately 500 times the box size, its sign gets changed, i.e., it's like some ulterior PBC are applied to the system (see attached plot, showing the coordinates of an atom, rescaled by dividing for the box length, vs the number of steps. Notice that I'm running an NPT simulation using "fix nvt", therefore the box length is not constant).

I report below my input script. If some of the commands I am using implicitly enforces this limitation on the atom coordinates, I am not aware of it.

Could someone help me understand what is causing this?

that is because LAMMPS stores coordinates as coordinates inside the
box plus image flags. every time an atom is crosses through a periodic
boundary, its image flag for that direction is incremented or
decremented accordingly. this is useful, as otherwise the resolution
for representing coordinates would be reduced every time the absolute
value of a coordinate doubles.

you are seeing the sign change, because those image flags are under
normal circumstances 10 bit signed integers. so after 512 passes
through the box, the sign flips.
if you compile LAMMPS with -DLAMMPS_BIGBIG, those image flags are
expanded to 21 bit signed integers and you can have instead 1048576
passes through the periodic boundaries before it wraps around.

axel.