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