SHAKE | Does FixShake::post_force() calculate the exact correction forces when combined with Velocity-Verlet?

Dear LAMMPS developers,

I have a question concerning the computation of the force correction in the FixShake::post_force() routine. In order to explain the problem, let us consider a NVE simulation of SPC/E water:

Say we are at the beginning of the update step n->n+1 (input: v^n, x^n; output v^n+1, x^n+1). In FixNVE::initial_integrate() the velocity is integrated by half a timestep, from v^n to v^n+1/2. The coordinates are then integrated by a full timestep using the velocity v^n+1/2, i.e. from x^n to x^n+1.

Afterwards, in FixShake::post_force() the idea is to calculate the Lagrange parameters and force corrections, such that the SHAKE constraints are satisfied after the next integration step, from x^n+1 to x^n+2. (x^n+2 should satisfy the constraints.) In FixShake::post_force() an estimate of x^n+2 is calculated ( = xshake) and on that basis the Lagrange parameters and force corrections are determined. Am I right?

After FixShake::post_force() the velocity is integrated twice for half a timestep again (from v^n+12 to v^n+1 and from v^n+1 to v^n+3/2), before the coordinates are eventually integrated again, i.e. from x^n+1 to x^n+2.

The thing that confuses me is that the correction forces in FixShake::post_force() are calculated on the basis of the velocity v^n+12. However the coordinates are eventually updated by v^n+3/2 and therefore I cannot see, why this procedure should satisfy the SHAKE constraints exactly for x^n+2. With other words, xshake only approximates x^n+2 and therefore I would assume that the coordinates x^n+2 satisfy the constraints only approximately.

Is this intended? Is there anything else going on, which I did not notice?
A while ago I implemented SHAKE in combination with normal Verlet and there the coordinates always satisfied the constraints exactly.

Any comments are much appreciated!


Have a look at FixShake::unconstrained_update()

The update is to adjust forces so the positions will be correct after the next timestep’s position
update (at the beginning of the step), which is a full step of velocity
from now.

You can also turn on fix shake output and see how close the atom coords

are to the exact result. Since SHAKE is iterative it will never be “exact”
but it gets pretty close.


Hi Steve,

thanks a lot for your quick reply. Actually, I had a look at "FixShake::unconstraint_update() which motivated the question:

for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
dtfmsq = dtfsq / mass[type[i]];
xshake[i][0] = x[i][0] + dtvv[i][0] + dtfmsqf[i][0];
xshake[i][1] = x[i][1] + dtvv[i][1] + dtfmsqf[i][1];
xshake[i][2] = x[i][2] + dtvv[i][2] + dtfmsqf[i][2];
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;

Here xshake is calculated and it should approximate the coordinates after the next (unconstrained) integration step. Based on this approximation the Lagrange parameters are then calculated and the forces adjusted accordingly such that the updated coordinates (r^n+2) should obey the constraints after the update step.

If the velocity would not change in-between FixShake::unconstrained_update() and FixNVE::initial_integrate(), the constraints should be satisfied exactly. However, xshake is not the same as the unconstrained coordinates you would get in the next update step without correction forces because the velocity changes and therefore I still don’t see how this should work.

Could you please suggest any references according to which the algorithm was implemented?


PS: I mean exact w.r.t. the specified tolerance.

The NVE update at the end of step is

v’ = v + dt/2m * f

The NVE update at the beginning of the next step is

x’ = x + dt v’

That’s exactly the equation for the unconstrained update
in your email, is it not? What am I missing?

As the header of the method states, the update assumes NVE.


Hi Steve,

did you possibly forget about the additional velocity update in FixNVE::initial_integrate() right before the coordinate update?

1.) post_force: fShake

sorry, I was sloppy.

See these lines of code in setup()

// half timestep constraint on pre-step, full timestep thereafter
dtv = update->dt;
dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
dtfsq = update->dt * update->dt * force->ftm2v;

dtfsq is a full step, not half-step, i.e. there is no factor of 1/2

used in unconstrained_update(), except on the setup constraint before
the 1st step.


Hi Steve,

thanks for pointing out that dtfsq changes, I didn’t realize that.

However, I cannot see that this actually resolves the issue. Before post_force(vflag), we have dtfsq = 0.5update->dtupdate->dt*force->ftm2v, and therefore

xshake = x + dt * v + dt*dt/2m * f

as assumed previously.


Sorry, Steve. Did you mean that after the first half-timestep (setup), reset_dt is invoked at some point and throughout the simulation we have dtfsq = update->dt * update->dt * force->ftm2v?

Best wishes,


Dear Steve,

you were right and the problem is therefore resolved.
Thanks for your patience!

Best wishes,