Fix deform screw restarts and what to do with it?

I encountered to what seems to me two bugs in restarting with fix deform command.

To understand them a quick recapitulation of what lammps is doing during a timestep for what is important to us in this bug report:

  1. Update velocities: v(t+dt/2) = v(t) + f(t)*dt/2 <= This assumes we know the forces from the previous timestep!
    Update positions: x(t+dt) = x(t) + v(t+dt/2)*dt

  2. Compute forces f(t+dt) using the velocities v(t+dt/2) and positions x(t+dt) computed in 1)

  3. Update velocities: v(t+dt) = v(t+dt/2) + f(t+dt)*dt/2

  4. Change the box size according to fix deform rules

  5. Write restart data to file

Lammps do not write forces to the restart files (I am actually not sure why?). Therefore, if we restart a simulation, before being able to do 1), lammps need to recompute the forces. When lammps do that it will obtain forces that are different from the one it computed at the previous timestep in 2) for two reasons. First, the velocities will be different by half a timestep. This is relevant only to the granular pair styles (where forces can depend on velocities), is not depended on fix deform and is correctly pointed out in the manual. So no point of talking about it here. However, there is a second problem, the box size (or shape) has been changed. The magnitude of the force error, which will arise from this, will obviously depend on the type and speed of deformation as well as the pair style. I assume it is very small in most cases. However, in my case of simulating a compression of loose granular particles, the smaller box create many new “spurious” contacts that arise between particles at the periodic border. The force can be sometime overestimated by several orders of magnitude!

I suppose that this bug was at the origin of this reported abnormality problems about "restart+fix deform" and the one I reported myself a long time ago Perfect restart for granular pair style . Though certainly I did not understand at that time what the origin of it was.

This bug seems not to be directly related to the one reported here restarts + fix deform and here assumed bug in fix deform erate upon restart and which has already been addressed.

The dirty hack solution in the script file that I found is to change the box size one step back, compute the forces, change box size once again and then continue the simulations. For my case of a true rate compression in 2D with rate gamma and timestep dt, an example script will read as follows:

read_restart restart_file
fix iso_comp all deform 1 x trate {gamma} y trate {gamma} units box
run 0 #To be able to change box size after reading restart
variable back_scale equal exp(-{gamma}*{dt})
change_box all x scale {back_scale} y scale {back_scale} units box
run 0 post no #Compute the correct pressure
variable forward_scale equal exp({gamma}*{dt})
change_box all x scale {forward_scale} y scale {forward_scale} units box
run 1000 pre no #Continue simulation without recomputing the forces

I suppose that at minimum this problem should be mentioned in the manual both in the fix deform section and read_restart section. In addition, the previous script “solution” could be mentioned in the manual. A better approach, of course, will be somehow to modify the code to resolve this issue. One approach will be for fix deform to detect the same fix, and scale the simulation box size back for the first force computation. However, this will not work if I wanted to change the type of deformation. For example, I actually want to do first some compression and then shearing. In that case the fix deform commands will be different and the previous approach will not work. However, once again it could be combined with a script hack, by first specifying the old deformation, computing the forces with a “run 0” the changing the fix deform and continuing the simulations.

There is another bug if one makes a restart while using fix deform, which only apply when one use lammps for granular simulations. The latter can have a force, which depends on the relative difference in velocity between two particles, which means that “remap velocity” options should always be used with periodic borders, else huge spurious forces will arise between contacting particles at the border. Lammps correctly remap the velocity for forces calculations during a run. However, when the force is computed at the beginning of a restart, the velocities are not remaped. And depending on the box size, this can add huge incorrect forces.

Contrary to the previous case, I do not see how this can be remediated by some Band-Aid in the script file. The only script solution that I see is to store the forces and reapply them at the beginning of restart. I assume that it can be done (have not yet tried it) by using fix property/atom to store the forces in a restart file. Then at restart I could use fix setforce to write these forces to the atoms with a “run 0”, then unfix setforce and continue the simulation without recomputing the forces. However, this will only work for future simulations where I would have stored the forces, no way I could use that for my previously saved ones. If somebody will come with a better solution for this one (especially one that will allow me to use the old saved configurations), I will be happy to hear it.


Sorry, I don’t get what you are after outside of the fact, that LAMMPS was not designed to handle restarts for granular simulations. Perhaps @sjplimp has some more specific comments and suggestions.

Hi Alex,

The first problem affects any type of pair correlations, not just granular simulations.

The box size (shape) is different for the initial force computation after a restart than the one used to compute the forces if no restart have been performed. This means that the distance between neighboring particles interacting through the periodic border will be different. This will change the force computed by any distance dependent potential.


It is difficult to discuss this without having a proper example demonstrating the issue.
Ideally a very simple input with just a few atoms and a lj/cut pair style where one then can dump all positions, velocities and forces and compare. Also the proper place to report this is the github issue tracker

Ok, I will prepare an example script and create an issue on github.

Thanks for the details explanations Anton. You mention 3 problems.

  1. inexact restart with velocity-dependent potentials, e.g. granular pair styles.
    As you say this is a known issue mentioned in the doc pages (read_restart and granular pair styles).
    I think usually this is a very minor inexactness but otherwise statistically identical results.
    E.g. akin to restarting with a different random number seed for DPD or a Langevin thermostat.
  2. change in box size upon restart when using fix deform. This one I do not understand.
    Fix deform makes its changes to box size/shape at the end of the timestep. So the new
    box size/shape should be saved in the restart file and you should thus restart with the
    same size/shape as if the simulation had continued for more steps. Is this due to some
    discontinuity in the formula for deformation that fix deform is using as a function of time
    in the old run versus the new run?
  3. issue with velocity remap in fix deform when restarting. This one I also don’t understand
    clearly. Are you saying the new fix deform in the new restarted simulation does not
    remap velocities when it creates ghost atoms?

Simple example scripts would be helpful.


Hi Steve,

We agree on point 1.

For the point 2, the new box shape will be saved to the file, and this is the box shape that will be used to compute the initial forces after the restart. However, the box shape that was used to compute the forces of the previous timestep was different. Thus, the forces that will be used at the beginning of the run after a restart are different from the ones that would have been used, if the run would not had been stopped.

For the point 3, the new fix deform will remap the velocity of the ghost particles, except for the initial forces computation.

I have prepared scripts to illustrate the problem, but I am not authorized to upload them here. I will open an issue on github and upload them there.