fix_shake: a hybrid method between SHAKE and RATTLE?

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users and developers,

I realise that in the past questions about the exact LAMMPS implementation of SHAKE were raised for several times, (1-3), but unfortunately my question remains open.

(1) http://lammps.sandia.gov/threads/msg10867.html

(2) http://lammps.sandia.gov/threads/msg24771.html

(3) http://lammps.sandia.gov/threads/msg44269.html

Judging by the fix_shake implementation in LAMMPS (and by Steve’s answers in (1) and (3)), it seems that the Lagrange parameters are calculated such that the coordinates satisfy the constraints at time n+1, as originally suggested in Ref. [1] for the Verlet algorithm. However, as far as I understand, the advantage of RATTLE (Ref. [2]) is that the time derivatives of the constraints are satisfied as well for which they introduce additional Lagrange multipliers (see attached screenshot “RATTLE.jpg” of the relevant section in [2]). This is important for the velocity integration.

To me it seems that the algorithm in LAMMPS is a hybrid between SHAKE and RATTLE, for velocity Verlet is used, but the time derivatives of the constraints are not taken into account. If this was indeed the case, the velocity integration would be slightly inconsistent. I was wondering, are there any studies on the effects of the additional velocity correction in RATTLE that would justify ignoring the velocity correction? It might be a small effect anyway, but I was hoping that somebody could comment on this question.

Best wishes,

Peter

References:

RATTLE.jpg

2

I was hoping someone else would answer.
I don’t recall the differences between SHAKE and RATTLE.
But I don’t see what the issue would be with velocities.
SHAKE, as implemented in LAMMPS, simply
corrects forces, so that the next
velocity/position update will satisfy the distance
constraints. I think it is a correct implementation
for velocity Verlet.

Steve

3

I was hoping someone else would answer.
I don’t recall the differences between SHAKE and RATTLE.

As far as i recall: shake applies only to regular Verlet integration, velocity Verlet requires rattle in addition to maintain the constraint.

4

Dear Steve,

thanks for your reply. The questions really is, does LAMMPS explicitly correct the velocities as suggested in the Andersen paper?

I would appreciate, if you could have a quick look at the relevant page of that paper:

I have the impression that fix_shake only uses g_{RR} as suggested in (2.7). That guarantees that the coordinates r_i(t+h) satisfy the imposed constraints. For Verlet that is enough because the velocities follow from central differences along the trajectories. However, according to the paper another constraining force, g_{RV}, is necessary for Velocity Verlet in order to make sure that the integrated velocities satisfy the time derivatives of the constrained coordinates. Either way the coordinates will always be fine, because g_{RR} does its job. However, that does not mean that the velocities are integrated correctly without g_{RV}.

Maybe I am missing something, but I think fix_shake only satisfies (2.7) but not (2.8), in which case the velocities would not satisfy the correct constraints. In this case the velocities would be slightly off and just by monitoring the imposed bond lengths one could not tell that anything is wrong.

My apologies if I am misunderstanding how fix_shake works. I just think that if I was right, this might be a problem.

Peter

5

My apologies for posting this twice, but I forgot to attach the paper in the previous post.

Peter

rattle.pdf (59.8 KB)

6

Dear Axel,

thank you for your reply. That is also what I understood and therefore I was hoping that LAMMPS would use RATTLE even though the fix is called shake. However, the implementation actually suggests that it uses SHAKE, if I interpreted the code correctly.

Peter

7

I looked at the page of the Anderson paper you sent,
but don’t get the full context as to what he is indicating
is a problem. I don’t have access to the full paper on travel.

Fix shake as implemented in LAMMPS is doing the
following. It adds a constraint force in the middle
of the timestep (after interatomic forces and other constraint
forces have been computed, such as wall forces). That
constraint force is “exact” in that the next time positions
are updated (at the beginning of the next time step), the
SHAKE distance contraints will be satisfied (to high precision).
B/c that force is added to the existing forces it is applied
at the end of the current time step to change the velocities
(by 1/2 step), and again at the beginning of the next time step
to change the velocities (by another 1/2 step), before those
velocities are used to update the positions.

So the SHAKE constraint forces change both the positions
and velocities.

Steve

8

Dear Steve/Axel,

thank you for having a look. Sorry that I didn’t attach the full paper (see andersen_1983.pdf).

Steve, the updated forces do affect the velocities as well, but not in the correct way. Since the molecules are rigid, the velocity component along the bond vector should vanish. This is the condition (A2) in the attached paper and the RATTLE algorithm makes sure that this is indeed the case.

I carried out several short simulations with LAMMPS (1Feb14) to show that there is indeed a non-negligible velocity component along the bonds in a simulation of 108 SPC/E water molecules with a SHAKE tolerance of 1.e-10. The error is quadratic in the timestep, which suggests that it is due to the missing velocity correction in the RATTLE algorithm. Please have a look at shake_problem.pdf.

Best regards,

Peter

shake_problem.pdf (33.2 KB)

andersen_1983.pdf (613 KB)

9

ok - makes sense. Is it correct to say that LAMMPS
implements SHAKE (as advertised) within velocity Verlet,
but not RATTLE?

If someone (you?) wants to implement a fix rattle, perhaps
as a derived class of fix shake, that adds the velocity
correction in Q, go for it.

Steve

10

Dear Steve,

thanks for your reply. I don’t want to cause any confusion, where there should not be one:
Of course, LAMMPS does exactly what is advertised in the documentation. It combines SHAKE with velocity Verlet.

I don’t really know, how the velocity correction in RATTLE affects the dynamics of the system. I am sure that many people have carefully validated their results with literature data and if there is indeed an effect, then it must be very small. For example, the radial distribution function and the VACF of SPC/E water do not seem to be affected at all, as the comparison with my in-house code, which uses Verlet & SHAKE, shows (please see gr.eps, vacf.eps).

I think it should not be too much work to implement RATTLE based on the current SHAKE implementation. I can think about it and give it a go. However, I don’t know the LAMMPS code well. In particular, I don’t know exactly at which point the velocity correction should be carried out, so that it does not interfere with other fixes. I might need some assistance for that. If anyone of the more experienced developers would be willing to prepare the scaffold for the new fix to speed things up a little, I would be happy to think about the implementation.

It would definitely be a good idea to include the RATTLE algorithm. Maybe this adds a bit too much to the wish list now, but it might also be a good time for including the SETTLE algorithm, which is the analytical solution for 3 constraints. It is faster and more accurate than RATTLE and SHAKE.

Thanks again for your help,

Peter

gr.eps (36.1 KB)

vacf.eps (25.4 KB)

11

sure, we can help with an implementation of a fix rattle.

If RATTLE is effectively SHAKE + a velocity correction,
then you could do this:

a) create a fix rattle that derives from fix shake
b) fix rattle will only need one method
c) figure out where in the time step the velocity
correction should occur - probably at the same
time the SHAKE forces are added, which is
treated as a force correction by LAMMPS, which
occurs right after the interatomic forces (and other
constraint forces, like wall forces) are calculated
d) this means fix rattle implements simply a post_force()
method - if could call the parent (SHAKE) post_force()
method, which will compute the SHAKE forces, then
calculate the velocity correction and add it to the velocities

that might be all that is needed

Steve

12

Thanks, Steve!

I will give it a go, but it might take a few weeks until I can find time for that.
Will get back to you as soon as I run into difficulties.

Peter

13

Dear Steve, Axel and all other LAMMPS developers,

It has been a while since we talked about the implementation of RATTLE in LAMMPS. Finally, I found the time to look into that. The velocity constraints are now satisfied exactly (analytically) at the end of the timestep and the coordinate constraints up to the desired tolerance. SHAKE still works normally.

I decided to follow a slightly different approach to the one you suggested, Steve, because it turned out that you cannot correct coordinates and velocities in a single post_force method. Here is what I have done so far:

a) I derived FixRattle from FixShake and changed the isolation level of the member variables from private to protected.
b) I implemented FixRattle::initial_integrate(), which simply calls FixShake::post_force() and takes care of the constraining forces for the coordinates.
c) I implemented FixRattle::post_force(), which applies the velocity corrections.
d) For testing purposes I also implemented FixRattle::end_of_step() to check whether the constraints are actually satisfied.

At this stage RATTLE only works in serial and the tests I have carried out so far only include the constraints and the pair correlation function. The code is written in parallel, of course, with FixShake as a blueprint, but I am getting an error with atom->map. To be more precise, if I run on more than one core the indices I get from

int i0 = atom->map(shake_atom[m][0]);
int i1 = atom->map(shake_atom[m][1]);
int i2 = atom->map(shake_atom[m][2]);

sometimes exceed the number of atoms consequently leading to a segmentation fault (please find more details in the README file in the attached archive). Presumably, I am missing a synchronisation or mapping step somewhere. This should be fairly obvious for someone who knows LAMMPS well, but I am having a hard time finding the bug. I attached an archive containing a README file, documented source files and input scripts for a SPC/E water test simulation. For my tests I used lammps-9Dec14. In case I was not entirely sure about certain lines of the code, I left a comment “// QUESTION”.

Could anyone please help me look into this?

Best regards and thanks,

Peter

PS: Needless to say that this is just a preliminary version of RATTLE which should not be used for real simulations until tested extensively.

rattle_v1.tar (800 KB)

14

Dear Steve, Axel and all other LAMMPS developers,

[...]

At this stage RATTLE only works in serial and the tests I have carried out
so far only include the constraints and the pair correlation function. The
code is written in parallel, of course, with FixShake as a blueprint, but I
am getting an error with atom->map. To be more precise, if I run on more
than one core the indices I get from

  int i0 = atom->map(shake_atom[m][0]);
  int i1 = atom->map(shake_atom[m][1]);
  int i2 = atom->map(shake_atom[m][2]);

sometimes exceed the number of atoms consequently leading to a segmentation
fault (please find more details in the README file in the attached archive).
Presumably, I am missing a synchronisation or mapping step somewhere. This
should be fairly obvious for someone who knows LAMMPS well, but I am having
a hard time finding the bug. I attached an archive containing a README file,

it looks like you don't have any Fix::grow_arrays(),
Fix::copy_arrays(), Fix::update_arrays(), Fix::set_arrays() methods in
your class.

your vp[][] array needs to be grown and updated as atoms migrate
between domains (just as xshake[][] in fix shake).

axel

15

yup. grow_arrays is what you are missing. please try out the attached
modified code.

axel.

fix_rattle_mods.tar.gz (7.24 KB)

16

Dear Axel,

Thanks for your quick and helpful reply! Yes, you are right. I did not overwrite the methods grow_arrays, copy_arrays,… etc., but only created vp once in the constructor of FixRattle:

memory->create(vp, atom->nmax,3,“rattle:vp”);

What I still don’t understand when I run your updated version is the value of nmax when FixRattle::grow_arrays() is called. Here is the call stack (running on 1 core only):

(1) FixShake::FixShake()
(2) FixShake::grow_arrays() -> nmax = 768
(3) FixRattle::FixRattle()
…
(x) FixRattle::grow_arrays() is called from somewhere -> nmax = 16384

I thought nmax should be smaller or equal to the number of atoms in the system. Even with SHAKE alone, at step (x) FixShake::grow_arrays() is called with nmax = 16384 from somewhere outside the class.

Does nmax = 16384 make sense for a system consisting of 256 water molecules (768 atoms)?

Thanks,
Peter

PS: It’s surprising that I didn’t get a segmentation fault all the time…

17

Dear Axel,

Thanks for your quick and helpful reply! Yes, you are right. I did not
overwrite the methods grow_arrays, copy_arrays,... etc., but only created vp
once in the constructor of FixRattle:

memory->create(vp, atom->nmax,3,"rattle:vp");

What I still don't understand when I run your updated version is the value
of nmax when FixRattle::grow_arrays() is called. Here is the call stack
(running on 1 core only):

(1) FixShake::FixShake()
(2) FixShake::grow_arrays()
-> nmax = 768
(3) FixRattle::FixRattle()
........
(x) FixRattle::grow_arrays() is called from somewhere -> nmax =
16384

I thought nmax should be smaller or equal to the number of atoms in the
system. Even with SHAKE alone, at step (x) FixShake::grow_arrays() is called

nope. nmax has to be larger or equal than the largest number of
*local* atoms across all domains. it is grown in increments of 16384.
please see AtomVec::grow_nmax()

with nmax = 16384 from somewhere outside the class.

Does nmax = 16384 make sense for a system consisting of 256 water molecules
(768 atoms)?

yes. since the avec->grow() is triggered somehow for parallel
calculations. all fixes that maintain similarly distributed per atom
data (like xshake or vp) have to be notified of that, hence you have
atom->add_callback(0); in the constructor. this registers the fix with
the atom vec class and then has the fix->grow_arrays() method being
called at part of avec->grow(). check out the block:

  if (atom->nextra_grow)
    for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
      modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);

at the end of AtomVecFull::grow().

so it all makes sense, when you unfold this step-by-step. this
complexity is the price to be paid for the flexibility of the overall
architecture in LAMMPS that makes it so easy to extend its
functionality with code that (for the most part) only has to worry
about the feature at hand.

HTH,
     axel.

18

Dear Axel,

Thanks for the detailed explanation! I will have a close look at these methods.

Best wishes,
Peter

19

nmax has to be larger or equal than the largest number of
local atoms across all domains. it is grown in increments of 16384.
please see AtomVec::grow_nmax()

If by “across domains” you mean processors, then nmax

is defined by individual processors. It is just the length
of the allocated per-atom arrays on each proc, and can
be different on each proc. It grows in chunks of 16384 just

so that there are not continual reallocations. I.e. bump
it up once in a while be a reasonable amount when needed.
It never shrinks.

Steve

20

Dear Steve and Axel,

Thanks for your comments and suggestions! I finally got my head around the communication of arrays in LAMMPS using forward_pack and forward_unpack. At this stage, RATTLE works now in serial and parallel and gives identical results.

However, my RATTLE implementation requires twice as much communication as SHAKE. To be precise, during each timestep the arrays xshake (unconstrained coordinate update) and vp (unconstrained velocity update) are communicated. Before implementing the missing methods (vrattle2, vrattle3, vrattle4), I wanted to ask you, if either of you could please have a quick look at the code and tell me, whether you see any mistakes or potential for optimisation. I also added some in-line documentation with additional questions (just grep for QUESTION to get to the relevant sections) and made some minor changes in fix_shake.h/cpp.

Thanks,

Peter

fix_rattle.cpp (17.7 KB)

fix_shake.cpp (82.1 KB)

fix_rattle.h (2.2 KB)

fix_shake.h (8.41 KB)