I have identified a bug with LIGGGHTS, which is causing problems every time a run command is encountered in the input script. What I have seen 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 (ie in the way the last value of shear gets passed on to the pairs at the initial step after run Nsteps is found).
The file I attach reproduces this behaviour. It creates a text file which outputs the coordinate of one grain. If at the end of the input file you merge the run 5000 and run 5000 to a single run 10000 command you will see that the output in the text file is different, with a jump in particle position at the beginning of the second run 5000 command.
I have included comments in my input file, but basically what I am doing is just creating a packing, choosing walls and setting grain properties. I chose this configuration as I am not sure how clear the difference would be using different packings. I basically set the parameters to strange values so that I make sure I get no damping and I am not sure if the timestep is 100% correct, but this should not be an issue.
In order to solve this problem we need to check what happens with the fixes and neighbour list when a simulation is restarted with a run command.
A bit of looking into the code has pointed me to the function Verlet::setup . This gets run and 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..
Ideas/help on how to get rid of this bug are gratefully appreciated. I have already spent some time on it but did not get very far.