Additional FixShake::end_of_step() routine

Dear LAMMPS Developers,

I am interested in carrying out a NVE simulation of SPC/E water with a custom end_of_step()-Fix. That fix makes a tiny modification to velocities and coordinates and I ran into a problem when trying to extend FixShake in such a way that it can take the changes into account. Let’s assume that the fix has the following form:

void FixSomefix::end_of_step() {
// modify x and v a tiny bit
domain->pbc();
}

In FixShake::post_force() the forces are modified such that the coordinates will satisfy the constraints after the next integration step exactly (w.r.t. the specified tolerance). However, if FixSomefix::end_of_step() modifies velocities and/or coordinates after that, the forces need to be adjusted again to take the modifications into account. If user-defined fixes modify the forces before FixShake::post_force() is invoked, everything should be fine. However, in my case I want to do that at the end of the timestep to avoid other complications.

I was thinking of the following procedure:

Comments below.

Steve

Steve, thank you for your comments and your help! Actually, I got it to work right now after I removed domain->pbc() and another bug (not in the methods of the previous posts, they were OK).

With regard to the motivation of doing this, I completely agree with you and I realized that even without calling shake twice, the geometry of the molecules was still preserved up to a tolerance of 1.e-5 (max. error). I am more than happy with that, but I needed to know that I can control the error, if it was necessary. It was more like a proof of concept. Now everything is fine!

Please allow me to ask you two more questions:

  1. When you said “better do that check in fix_shake”, were you referring to testing that the constraints are satisfied or to checking whether there are any other end_of_step()-fixes after fix_shake?

  2. How can I make sure that fix_shake::end_of_step() is the last end_of_step()-fix without messing around in modify, i.e. fix[list_end_of_step[i]]->end_of_step()? Should I just change the order in the input file or is there another way?

Best wishes,

Peter

Comments below.

Steve