Potential change for logical improvements: move control of constraints purely to the .mdp
Here's a suggestion that could control flow control within the program logic, though it might not be worth doing because of screwing up old build status.
Right now, whether constraints are applied are controlled in two places: the .mdp, and the topology. If [ settles ] are defined in the topology, then they are used, independently of what constraints is set as in the mdp.
Logically, it would make the most sense to have decisions about how to run the simulation be under the .mdp. So in this case, [ settles ] would only be applied if the proper keywords on constraints was selected.
Exactly how to handle defaults, warning messages, deprecations, etc. can be discussed if people think it's actually worth doing. Making control only through the .mdp would make some of the integrator tasks as little simpler, but it's not required.
This is related to an email that Berk brought up before:
Specifically Berk said:
"For future Gromacs versions I would like to have the [ settle ] sections
in the water topology replaced with [ constraints ] and let mdrun
decide if to use settle.
From a user point of view, this will not change much though,
since settle is only useful for water and a user (nearly) never
writes a topology for water."
#1 Updated by Mark Abraham over 3 years ago
Logically, it would make the most sense to have decisions about how to run the simulation be under the .mdp
It's sensible to have a topology (perhaps with programmatic modifications), and to apply algorithms that calculate forces and propagate coordinates with those forces. I don't yet see an intrinsic advantage to localizing the decision to treat bonds with constraints to either the topology or the algorithms.
From the point of view of the current implementation of grompp, it doesn't matter a great deal, because the
constraints mdp field specifies what bonded interactions should be converted to constraints. That could just as readily be preprocessing, like all the standard water .itp files have for flexible vs settle waters. But there's a distinction between
constraints and, say,
pcoupletype inasmuch as the former doesn't get expressed directly in the tpr.
It's conceptually nicer if mdrun (or grompp) can decide how to implement constraints, based on recognizing patterns, so that we don't have the question of implementing and documenting special topology directives for particular kinds of rigid bodies.
#3 Updated by Mark Abraham almost 2 years ago
There's several related issues here.
Berk wants the choice of settle vs lincs/shake to be up to mdrun, given that something has expressed that water molecules use constraints. I support that, but we haven't done it.
The user can sensibly want to specify particular constrained interactions in the topology, so we must support that.
Most water models we use are designed as rigid, and so be default they should be.
Currently the .mdp has a keyword to convert interactions to constraints (which perhaps Michael has misunderstood/forgotten). That makes sense to couple to things like the timestep choice also in the .mdp file. Setting the FLEXIBLE preprocessor symbol to de-constrain the waters is also a reasonable approach.
I suggest we reject most aspects of this change, as we are doing things that make sense.
#4 Updated by Michael Shirts almost 2 years ago
Here's my current thoughts. It makes sense if mdrun makes decisions, AS LONG AS the .log file clearly states, perhaps with * text *, or something, that clearly states the decision is being made, and it's possible to override. So if mdrun decides to use settle vs lincs, that's fine, as long as a message is printed to the log file, and there exists some MDP option that can force that change one way or the other. One can imagine that if there are two valid choices, someone will find a circumstance where they would rather have one or the other.
So in this case - [ constraints ] are part of the topology, and should be states in the topology. The algorithm to run can be decided by mdrun, but there should be some way to override it (bPreferSettle vs. bPreferLincs).
#5 Updated by Mark Abraham almost 2 years ago
The action of both FLEXIBLE and the constraints keyword are preprocessing stages. We could formalize that as a dedicated transformation of the topology (e.g.
gmx constrain -f topol.top -constrain selectionstring) so that there is an intuitive way to use a flexible topology during system preparation and to transition to a constrained one as we approach production stages. This fits a data model where the user chooses to have holonomic constraints on certain degrees of freedom, or not. Later, we could implement that as a transformation of the hierarchical key-value data structure that represents the system topology. Since the relevant thing is the mass of the particle, a selection string like "mass < 5" or "all" can be interpreted as choosing a set of atoms that trigger the conversion of all associated bonds to constraints. By default, water models are constrained, but one could use @gmx constrain -f -unconstrain selectionstring for selections like "resname SOL", or similar. (Also,
gmx convert-vsites could work on similar principles, rather than being bundled in pdb2gmx.)
The .mdp should retain the ability for the user to choose parameters for any constraint algorithms that (by default) mdrun chooses to implement the required constraints. We already have constraint-algorithm to choose shake or lincs, which is recorded in the .log file.
We should add something like "use-analytic-constraints = selectionstring" which permits the use of settle/rattle (according to the integrator flavour) for the constraints expressed in the topology that match the selectionstring. Add a line in the log file that confirms the resulting decision. This approach seems like it would generalize to any future special case implementations for molecules or integrators.