I haven’t followed the details of this thread.
But I think intrinsically that fix heat and fix shake
are not fully compatible. When you just add random
energy to atoms in a water molecule that is being constrained
by SHAKE, some of the velocity change you induce
will be directed along the bonds, which means the
next iteration of fix shake will take that motion
out via the constraint forces. So I don’t think
the total energy of the system will be conserved, i.e.
fix shake will remove energy. Which isn’t something
it normally does if you just run NVE with conservative
forces between molecules.
Paul may want to comment.
Steve