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:
-
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 -
Compute forces f(t+dt) using the velocities v(t+dt/2) and positions x(t+dt) computed in 1)
-
Update velocities: v(t+dt) = v(t+dt/2) + f(t+dt)*dt/2
-
Change the box size according to fix deform rules
-
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.
Anton