Hi!

I’m a PhD candidate at the Technical University of Madrid (Spain). We are running LAMMPS simulations on granular systems. Our protocol applies linear compression to 2D systems of viscous and Hookean particles. Then we analyze the statistical distribution of some stress components (we obtain them by using “compute stress/atom virial” command).

We have realized that if we use the same inputs and data file but we only change the timestep (*) we obtain solutions that have very different statistical distributions of the stress components (all the chosen timesteps are below the critical value).

We think that this result is an immediate consequence of the velocity-verlet integration scheme, the non-linear nature of these systems and the finite precision of computers: when timesteps are very small, after each time integration some overlappings are also very small and cannot be considered and, in consequence, the resulting sharing out of forces and the dynamic evolution of the system are different. We have prepared a paper dealing with this issue, and one of the referees has suggested that it could be a bug in the LAMMPS code, but we do not think so. Indeed, we have implemented the velocity-verlet algorithm with other code to model another (simpler) non-linear system and we have checked how its dynamic evolution does depend on the timestep. The thing is that in this particular case the stress distribution could be strongly determined by the timestep.

What do you think about possible bugs in the LAMMPS code?

Thank you very much!

Ignacio

(*): We have changed the number of integration steps too, so that the represented real time, force rate and total force are always the same.

e.g.:

xxxxxxxxxxxxxxxxxxxxx SIMULATION A:

units lj

dimension 2

[…]

pair_style gran/hooke 202200 NULL 550.0 NULL 0. 0

read_data inital.grn

pair_coeff * *

neigh_modify exclude type 2 3

neigh_modfy… […]

group balls type 1

group … […]

timestep 0.001

variable rampforce equal -500.0*(step/100000)

communicate single vel yes

fix 1 all enforce2d

fix 2 piston aveforce 0.0 v_rampforce 0.0

fix 3 rightwall aveforce v_rampforce 0.0 0.0

fix 4 wall freeze

fix 5 all nve

compute 3 all stress/atom virial

dump output all custom 10000 packings.out x y id vx vy c_3[1] c_3[2] c_3[4]

run 100000

xxxxxxxxxxxxxxxxxxxx SIMULATION B:

["]

timestep 0.0001

variable rampforce equal -500.0*(step/1000000)

["]

run 1000000

xxxxxxxxxxxxxxxxxxxxxxxx

```
```