I was wondering, why does fix shake have to come before fix NPT? Normally, the shake command has to come last in order to ensure that the constraint is applied appropriately as stated in the documentation, and I have been unable to find the reason for this special case. The description of the error message states that, “else the SHAKE fix contribution to the pressure virial is incorrect.” Why is this so?
"If you define fixes (e.g. fix efield) that add additional force to
the atoms after fix shake operates, then this fix will not take them
into account and the time integration will typically not satisfy the
SHAKE constraints. The solution for this is to make sure that fix
shake is defined in your input script after any other fixes which add
or change forces (to atoms that fix shake operates on)."
I don't really think of fix npt this way.
In any case, fix shake does make modifications to the forces acting on
the particles. (Consequently, it also updates the virial. See
fix_shake.cpp inear lines 1394, 1559, 1802, 2035.) I think the
calculation for the virial must be completed before any barostat (like
fix npt, nph, or press/berendsen) can be applied.
Anyone who knows more should correct me.
Different fixes operate at different points during the timestep. For ones
that operate at the same point, they operate in the order they are specified
in the input script.
This is why SHAKE needs to be the last fix, among those that apply
force constraints, e.g. fix efield. This is because SHAKE takes the current
force (including due to other fixes that changed it) and applies new forces
to satisfy the SHAKE constraints. If you don't satisfy this, by adding other
forces after SHAKE has, then the SHAKE constraints (e.g. bond lengths)
are not guaranteed to be satisfied.
All of this has nothing to do with fix NPT. It operates at different points
of the timestep, so there is no ordering issue between the 2.
I can't remember now why the order wrt fix NPT matters. It's possible
this is an out-dated error, since the virial is now computed somewhat
differently than in the past. The way to test would be to comment out
the error message and run both ways and see if the pressure changes
in the first few timesteps (in parallel or serial). If so, there must
be a simple reason that
I am forgetting ...
I just did the test I mentioned with bench/in.rhodo and saw no
difference in pressure or dynamics when swapping the order of the