Granular LAMMPS bug?

Dear users of the granular package

I have found a problem when I run the file I supply below (using the 15Mar11 version with the granular package enabled). I have added comments in to make it easier to inderstand, but I basically create a packing of overlapping spheres bounded by granular walls. I then apply viscous damping to remove any unwanted forces, and I am monitoring the grain z coordinate. Then I remove damping completely.

What happens is that I register different behaviour if I do run 10000 instead of run 5000, run 5000. This behaviour might be difficult to spot unless you run with zero damping, but is probably a bug.

So what I think is happening is that the run command deletes all the forces and recalculates them, and I think that something is not done properly in the way the shear force gets calculated for each interaction pair (ie in the way the last value of shear[] gets passed on to the pairs at the initial step after the command run Nsteps is found). I looked into the code a bit and indeed the function Verlet::setup seems to be running, which as I understand it deletes all forces and recalculates them, but also builds the neighbour list again and maybe (not sure) remaps the atoms to processors etc..

I am emailing because I have spent quite a lot of time trying to figure out how to do the force recalculation properly and I got lost- maybe recalculation needs to be disabled? I am not sure.. Help is gratefully appreciated.



# Created on 17/3/2011 for LAMMPS 15Mar version

I'll look into this - sorry I didn't respond to your earlier
email ...


I just posted a patch for this. As you suspected there was some
shear history information being lost when continuing a run with
a second run command. Also, a 2nd run computes the forces
again, which was updating the shear history as though additional
time had elapsed, so I put in a flag to prevent that. Now it
looks like your test gives the same answers with 1 run or 2.

I think in general these are subtle effects, just causing a slight
perturbation at the start of a 2nd run.


Dear Steve,

Thank you for your help. My input file now runs fine using pair_gran_hooke_history, and as you said the difference in behaviour induced by the bug is very small.

I had a look at the patch you posted for pair_gran_hertz_history.cpp. Lines 56-60 repeat the same if else statement. I wonder if what you had in mind of inserting there was the following, as you did in pair_gran_hooke_history.cpp

+ int shearupdate = 0;
+ if (update->ntimestep > laststep) shearupdate = 1;

yes - I made the change - it will be in the next patch.