fix_shake: a hybrid method between SHAKE and RATTLE?

I don’t see any basic problems with what you are doing. Since FixRattle derives from
FixSHAKE, you don’t need to prefix anything with FixSHAKE:: like
FixShake::dtfsq = dtfsq_full/2.;

unless you are calling the same parent method from within a child method.

Since you are defining an initial_integrate() that adjusts forces, how are you
guaranteeing it will be called before fix nve’s initial_integrate (or fix nvt, etc)?

I also don’t really understand why another update is done at initial_integrate().
Could it be done at the previous step’s end_of_step() or even post_force(), or

does that turn it back into SHAKE?

Another thing to look at, is how FixSHAKE does something
different during setup (before the 1st step) versus during timestepping.
I.e. this chunk of code in setup()

// half timestep constraint on pre-step, full timestep thereafter

if (strstr(update->integrate_style,“verlet”)) {
dtv = update->dt;
dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
post_force(vflag);
dtfsq = update->dt * update->dt * force->ftm2v;

I don’t know if there is a similar need for RATTLE.

Steve

Dear Steve,

thanks for looking at the code and for your suggestions. This time I will
reply in-line.

I don't see any basic problems with what you are doing. Since FixRattle

derives from
FixSHAKE, you don't need to prefix anything with FixSHAKE:: like
FixShake::dtfsq = dtfsq_full/2.;
unless you are calling the same parent method from within a child method.

Sure, I just wanted to make it easier for others to read the code. I will
remove the prefix at the end.

Since you are defining an initial_integrate() that adjusts forces, how are
you
guaranteeing it will be called before fix nve's initial_integrate (or fix
nvt, etc)?

This is something I will have to check explicitly at the beginning of the
simulation. I don't really like it and would prefer a
*pre_initial_integrate()* method, but as far as I am aware there is none.

I also don't really understand why another update is done at
initial_integrate().
Could it be done at the previous step's end_of_step() or even
post_force(), or
does that turn it back into SHAKE?

Yes, the coordinate constraining forces could also be calculated in the
previous step's *end_of_step()* function. The reason why I decided to go
for *initial_integrate()* was that I thought it would be conceptually
simpler if everything that is relevant for the current timestep happens
within the current timestep. However, this is just my personal taste.

I don't think that it is possible to do everything within the *post_force()*
method alone. In fact, I tried to do that, but I encountered the following
problem: Consider the sequence of steps in *FixRattle::post_force()* for
the integration step *t^n -> t^n+1*:

1.) The coordinates *r^n+1 *already satisfy the constraints at this stage.
We can call *FixShake::post_force()* to calculate constraining forces,*
g_r^n+1*, for the next coordinate update (*r^n+2*) and add them to the
force, i.e. f^n+1 := f^n+1 + g^n+1

2.) We can then carry out an unconstrained velocity update, *v^n+1/2 -> vp*,
using the updated forces. That allows us to calculate the constraining
forces *g_v^n+1*. Now we have the problem that we cannot modify the forces
any more, because then the coordinates at time *n+2* wouldn't be integrated
correctly. The only thing we can do is to correct *vp* with *g_v^n+1* and
integrate half a timestep backwards, i.e. redefine *v^n+1/2 := vp - dt/2m
(f*^n+1)*.

3.) With this scheme we have achieved that *(r^n+1,v^n+1)* with satisfy the
constraints at the end of the first timestep. However, since we modified
the velocities in step 2.), our predicted unconstrained coordinate update
(relevant for *FixShake::post_force()*) is not correct any more and
therefore the coordinates will not satisfy the constraints after the next
*initial_integration()* is called, i.e. *r^n+2* would be incorrect.

Another possible scenario would be to change the order of velocity update
and coordinate update, i.e. swap steps 1.) and 2.):

1.) Carry out an unconstrained velocity update, *v^n+1/2 -> vp*, using the
forces *f^n+1*. Calculate the constraining forces* g_v^n+1* and update the
forces, i.e. *f^n+1 := f^n+1 + g_v^n+1*. *FixNVE::final_integrate()* would
then yield velocities that satisfy the constraints.

2.) If we now call *FixShake::post_force()*, the forces will be modified
such that the coordinates *r^n+2* will be fine. But if we do that, the
velocities will not satisfy the constraints any more at time n+1.

Therefore, my conclusion was that the corrections have to happen at
different stages, unless we introduce some additional logic. If anybody had
a concrete suggestion for how you could do it in one go with the available
components, I would be happy to try it.

Another thing to look at, is how FixSHAKE does something
different during setup (before the 1st step) versus during timestepping.
I.e. this chunk of code in setup()

  // half timestep constraint on pre-step, full timestep thereafter

  if (strstr(update->integrate_style,"verlet")) {
    dtv = update->dt;
    dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
    post_force(vflag);
    dtfsq = update->dt * update->dt * force->ftm2v;

I don't know if there is a similar need for RATTLE.

I took that into account. It's now *if (!rattle) post_force(vflag);* There
is no need to do that with *RATTLE*, if we stick to the implementation I
proposed.

Thanks,

Peter

Dear Steve,

Here is one more thing I realised when testing RATTLE:

If you use SHAKE with NVT or NPT integration, the coordinate constraints are actually not satisfied up to the target precision (unless you set the tolerance to a sufficiently low value, like 10^-4 for NPT). Is this the expected behaviour?

I will probably not look into that, but wanted to point it out in case this should not affect SHAKE.
For the same reason you can only expect RATTLE to satisfy the constraints correctly for NVE integration, as it uses SHAKE.

Best regards,
Peter

yes, the reason for that is that FixShake has no simple way to

predict the future positions based on what fix nve or nvt or npt will

do (since they modify velocities or the box size themselves).
So it assumes fix nve, and hopefully the difference versus nvt/npt

is small. I don’t see an easy workaround for that, since FixShake

does not even know which atoms will be integrated by which fixes.
E.g. the user could have some atoms integrated by fix nve, and others
by fix nvt.

Steve

Since you are defining an initial_integrate() that adjusts forces, how are you
guaranteeing it will be called before fix nve’s initial_integrate (or fix nvt, etc)?

This is something I will have to check explicitly at the beginning of the simulation. I don’t really like it and would prefer a >pre_initial_integrate() method, but as far as I am aware there is none.

You can look at init() in fix_recenter.cpp to see how it insures it is the last
fix to call initial_integrate(), and do something similar to insure fix rattle

is the first.

Therefore, my conclusion was that the corrections have to happen at different stages, unless we introduce some >additional logic. If anybody had a concrete suggestion for how you could do it in one go with the available components, I >would be happy to try it.

As I understand it, the motivation for fix rattle is to insure two things.
First that the velocities at the end of the timestep (after final integrate) are
correct. Second that the positions after initial integrate (next step) are correct,

These corrections are both affected by the newly computed forces, which are
available as early as post_force.

So it seems like you have 3 choices. (a) is what you are doing now (I think?).

(b) and © might allow for both corrections to be done together
and possibly allow for one communication operation instead of two

(a) in post_force, change the forces so that the final_integrate will
produce correct velocities, then just before initial_integrate, change
the forces so the NVE initial_integrate will produce correct positions

(and also the correct half-step velocities?)

(b) in post_force(), change both the forces and velocities, so that

the NVE final_integrate() and initial_integrate() will produce correct

velocity (at end of step) and positions (next step). It seems like you are trying
to control two outcomes (velocity and position) and you have 2

variables to do it with (force and current velocity).

© similar to (b) except do the operation in end_of_step(). The
velocities have already been updated in final_integrate (wrongly), but you can
correct them, and then change the force so that the

intitial_integrate on the next step will give correct positions

If (b) or © is workable, the choice might be affected by ease of

implementation or whether the virial contribution can be easily

computed. All things being equal, (b) is probably better, b/c

other fixes can be invoked at end_of_step, like diagnostic
calculations, and it would be better to assume that velocity
and force are already “correct” at that point in the timestep.

Steve

Dear Steve,

thanks for your suggestions.

Dear Steve,

I just wanted to add that my approach is indeed equivalent to your suggestion (a) and I get identical results if I modify solely the forces. I quite like that, thanks!
So to summarize up: At the moment I am using (a), but I would prefer to use (b), if it worked. Any more thoughts on that would be highly welcome.

Best regards,

Peter

Peter - not following all of your argument as to why (b)
will not work to enforce the constraints for both x,v in post-force.

I was thinking (b) could work, b/c of this logic:

a) at post-force, you have a correct x, current v, and new f
b) you want to insure a correct vnew after final-integrate (fi)
c) you want to insure a correct xnew after initial-integrate (ii)
d) at post-force, you can calculate the desired vnew and xnew
and you are free to adjust both v and f to insure (b) and ©

Assuming NVE updates, we know vnew,xnew will be:

(fi) vnew = v + dt/2m f
(ii) xnew = x + dt (v + dt/m f)

where (v + dt/m f) is the result of 2 half-step updates of v.

Eqs (fi) and (ii) are simply 2 eqs with 2 unknowns (v,f). So if you
adjust the current v,f to those solved-for values the NVE operations
will produce the desired vnew,xnew.

Basically, SHAKE was only trying to insure xnew is correct, so it only
adjusts the force in post-force. If you want RATTLE to insure xnew
and vnew are correct, you can adjust both the force and velocity in
post-force.

Am I missing something - or making a wrong assumption?

If so, why don’t we take the subsequent emails off-line until

we resolve things …

Steve