H bond constraint

[kexul]

Hi, recently I read a paper Best Practices for Alchemical Free Energy Calculations, which said that:

One issue which requires particular care is the use of constraints. Commonly, bonds involving hydrogen are constrained to a fixed length using algorithms such as SHAKE or LINCS, allowing the use of longer timesteps . However, in single topology relative free energy calculations, the atoms involved might be mutated to other atom types for example, in a mutation of methane to methanol, one hydrogen might become an oxygen atom. Typical molecular dynamics engines are not set up to recognize this change, or at least not to correctly include contributions to the free energy from changing constraints/constraint length, so results for a transformation would usually be erroneous. At present the most general solution to this problem is simply to avoid the use of constraints (and thus use a smaller timestep if necessary, usually of around 1 fs) in any relative free energy calculation involving a transformation of a constrained bond, as done by GROMACS.

I think decreasing the time step is easy to do in GROMACS, but how about avoid the use of constraints? Do I need to modify the topology file? Or should I set the constraints=none in mdp file? Have these things been handled when preparing the system by ProtoCaller?

Thanks!

[msuruzhon]

constraints=none in the mdp file will be enough to achieve this. If you are using ProtoCaller to run your simulations, then you can pass constraint='no' to the Protocol initialiser, as listed here or alternatively, directly into Simulation.runSimulation(), which automatically instantiates a Protocol before the run. If you do this, you should probably also decrease your timestep to 1 fs, since 2 fs are likely to crash without hydrogen constraints.