issue with pressure after write_restart

LAMMPS LAMMPS Mailing List Mirror
1

Hi-
When I write a restart file, the pressure afterwards is different, even though nothing has changed. Attached is a somewhat simple script demonstrating this. From what I can tell, it only happens when I use SHAKE. I added a dump command before and after and find that some of the atoms have moved very slightly. The attached diff.lammpstrj shows this. If anyone can help me figure out what is going on, that would be great.

Thanks, Steve

diff.lammpstrj (52.1 KB)

in.lammps (1.16 KB)

log.lammps (7.29 KB)

spce.mol (475 Bytes)

2

Hi-
When I write a restart file, the pressure afterwards is different, even
though nothing has changed. Attached is a somewhat simple script
demonstrating this. From what I can tell, it only happens when I use SHAKE.
I added a dump command before and after and find that some of the atoms
have moved very slightly. The attached diff.lammpstrj shows this. If anyone
can help me figure out what is going on, that would be great.

​the difference is quite small and, as you already speculated, due to fix
shake.

there are some general problems with fix shake when doing velocity verlet
time integration, as shake was derived for a regular verlet integrator,
where velocities don't show up in the time integration. thus fix shake does
a regular-verlet-like unconstrained ​position update before applying the
constraints resulting in sometimes in somewhat inconsistent velocities.
this is usually rectified after the first full MD step and the error is
small from then on. it is the largest, when you have initial coordinates
that are not fully representing the constrained geometry. the more correct
solution would be to also apply constraints to the velocities, which is
done with RATTLE (see fix rattle), but considering that you are only
solving the constraint equations with finite accuracy, it will probably be
not perfect either.

​axel.​

3

I see, thanks for the help.
-Steve

4

Peter Wirnsberger can also comment on this thread.

The following used to be the case, but I’m not 100%
clear that with Peter’s reformulation of fix shake in LAMMPS
it is still the case.

What is different

about restarting versus continuing a run,

is that a restart applies SHAKE for a half-step

prediction of the needed force constraints
before the run starts for use with velocity-Verlet.

Whereas in a continuing run, a full-step

constraint is applied once per step.

So due to the iterative nature of SHAKE, the

result will not be exactly the same for pressure
(SHAKE force contribution to virial), and the resulting

atom positions after the 1st step will be
slightly different, etc. But it should all be

withing the accurarcy tolerance specified

by fix shake.

Independent of restarts, I think this all may apply
for doing two “run 1000” commands versus

a single “run 2000” command as well. Peter can

verify.

This is noted on the read_restart doc page,

where fix shake is mentioned as a fix

that does not restart exactly.

Steve

5

Hi Steven, Steve & Axel!

I can confirm that the difference in pressure goes away as you increase
the SHAKE tolerance. To test this, I reran the input script with 'kspace
ewald 1e-12' (just to be sure) and two different SHAKE precisions
yielding the following results:

shake 10^-4: dp = 5.8 * 10^-1
shake 10^-12: dp = -4.2* 10^-9

As you can see, increasing the SHAKE precision by 8 orders of magnitude
decreases the difference in pressure, dp, by 8 orders of magnitude.

I therefore agree with Steve's explanation and think that you get
slightly different values for the Lagrange multipliers for the
half-timestep estimate at the beginning of the new run compared to the
full-timestep estimate at the end of the old run.

The constraining forces will therefore be slightly different and, as
they enter the calculation of the virial, the pressure will be different
as well. However, at timestep 0 the atom positions are not affected at all.

Since the entire fixed point iteration is a function of the timestep, I
don't see why you could expect to get exactly the same answer in both
cases for a finite number of iterations. It looks all reasonable to me.

Another cross-check would be to look at dimers with a single bond: In
this case you really shouldn't see any difference because
FixShake::shake() calculates the constraining forces analytically.

Best wishes,

Peter

6

Cool, thanks everyone.
-Steve