run 0 alters trajectory of granular triclinic simulation

Dear LAMMPS users,

I am running granular simulations and using run 0 to output the state of my system at a given point. I realized doing so altered my results. I could reduce it down to the following MWE (ran on version 07Aug2019):

units si
dimension 3
boundary p p p
atom_style sphere
comm_modify vel yes
newton off
neigh_modify delay 0 every 1 check yes # THERMO OUTPUT ALSO DIFFERS WHEN USING check no
neighbor 0.2 bin
region box prism -3 3 -3 3 -3 3 0 0 0 # NO DIFFERENCE IN THERMO OUTPUT WHEN USING: region box block -3 3 -3 3 -3 3
create_box 1 box
lattice sc 1
create_atoms 1 region box
variable drand atom random(0.25,1,1234) # NO DIFFERENCE IN THERMO OUTPUT WHEN USING: variable drand atom random(0.95,0.95,1234)
set group all diameter v_drand density 1.0
pair_style granular
pair_coeff * * hooke 1e3 0 tangential linear_nohistory 0 0
fix 1 all nve
fix 2 all deform 1 x trate -0.05 y trate -0.05 z trate -0.05 remap x
thermo_style custom step pxx pyy pzz pxy pxz pyz lx ly lz
thermo 10000
timestep 0.00001
run 1000000
unfix 2
run 0 # THERMO OUTPUT DIFFERS WHEN THIS LINE IS COMMENTED OUT
run 100000

That simulation isotropically compresses a 3D periodic cubic lattice during a first run. A second run letting the system evolve with no additional deformation follows. Between the runs, I do/don’t call ‘run 0’. The thermo output of the first run is identical in both cases. The thermo outputs of the second run are initially identical and after some time, there are some differences in the pressures at some of the timesteps whether or not I call ‘run 0’ e.g. on step 1100000 :

  • with run 0 : 1100000 73.190118 69.024216 74.024811 0.45726143 0.4994176 0.22467665 3.639184 3.639184 3.639184
  • without run 0 : 1100000 71.491376 68.62511 75.650158 0.27721563 0.03264701 1.7231778 3.639184 3.639184 3.639184

I have performed many variations of the previous MWE and noticed that:

  • using an orthogonal box, outputs are identical in both cases but different from the triclinic case (issue related to triclinic box ?)

  • using monodisperse particles of diameter 1, outputs are identical in both cases (issue related to polydispersity / cutoff ?)

  • using ‘neigh_modify check yes’, outputs differ in both cases (not a neighbor list problem ?)

Looking into the log files (attached), the simulations that differ during the second run have one extra ‘Neighbor list builds’ during that run. Comparing the cases with triclinic and orthogonal boxes, the ‘Neighbor list builds=’ for the successive runs are:

  • triclinic with run 0: 212, 0, 121

  • triclinic without run 0: 212, 122

  • orthogonal with run 0: 211, 0, 117

  • orthogonal without run 0: 211, 117

However, differences in thermo output are also observed when using ’ neigh_modify check yes’ although the neighbor list is built every timestep.

To my understanding, ‘run 0’ command or neighbor list building with no time integration should not alter the trajectory/precision of the computations. How to explain the fact that pressure differs also during the first run when using orthogonal or triclinic boxes although no tilt occurs? Is there something very obvious that I am missing or could there be a glitch with triclinic simulations of granular systems ?

I very much appreciate your input and help.
Jibril B. Coulibaly

Postdoctoral Fellow, Northwestern University

log_ortho_with.lammps (19.2 KB)

log_tri_with.lammps (19.2 KB)

log_tri_without.lammps (18 KB)

log_ortho_without.lammps (18 KB)

log_tri_checkno_without.lammps (18.1 KB)

log_tri_checkno_with.lammps (19.3 KB)

log_tri_mono_without.lammps (18.2 KB)

log_tri_mono_with.lammps (19.3 KB)