This is another opportunity to quote George Box:
It is inappropriate to be concerned about mice when there are tigers abroad.
I am sure there is a “correct” order in which thermostatting and rigid constraints should be applied. But both this order and the order’s effect will vary tremendously from system to system, especially on its phase and response characteristics (which are often the very properties to be measured, so what then?). This paper lists forty-two different orderings for the various manipulations in an NVT-SHAKE scheme.
So I would argue that it is impractical for LAMMPS to try and police these relative orderings, other than simply offering a common convention as described in the fix shake documentation. We could simply make additions (in bold) to the current text:
If you define fixes (e.g. fix efield, and including thermostats like fix nvt and fix langevin) 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.
and then just count on scientists to do their best. The impact of force field approximations, floating point truncation, PPPM inaccuracy, and others is likely to far outweigh inaccuracies from thermostat-constraint non-commutativity. I also suspect that most choices of timestep are either conservative enough that the effects will be limited, or large enough that the entire simulation destabilizes, and thus one would have to be unlucky (and careless) to land a timestep that is both small enough to return a largely stable simulation and large enough that flipping thermostat and constraints would give different answers.
Finally – these days I am seeing so many peer-reviewed papers sans their input scripts, and it drives me to despair. If I can’t see the script that was used to generate data, it’s really immaterial whether fix langevin came before or after fix shake, because neither choice is testable or reproducible relative to the paper data. It’s a little like running OSHA checks on the sturdiness of the deck chairs on the Titanic. The greatest need now is for strong community awareness and support for the most basic of standards – no peer review acceptance without at least a pretense of putting a script online – especially when it’s not long before someone at Google starts drafting their next Nature paper on “autonomous high-throughput non-equilibrium molecular dynamics to discover novel rheological regimes”, where they proudly tour a menagerie of LAMMPS simulations that meticulously sample complicated chemicals modelled in thoroughly unphysical ensembles. I’m personally too upset by this sort of thing to care so much any more about the finer points of O(t^2) errors from non-commutativity of integrator substeps.