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.
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.
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?
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
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?