Closed ppxasjsm closed 4 years ago
I haven't investigated this in depth personally. I think it is important to distinguish fully flexible setups vs setups where only the bonds involving hydrogen atoms are constrained. The rule of thumb seems to be that the timestep should be <x10 the period of the fastest vibrational motions.
The rule of thumb seems to be that the timestep should be <x10 the period of the fastest vibrational motions.
This is all I know as well.
a timestep of 1 fs is more an upper limit for MD simulations involving non-constrained bonds, I would rather suggest a value of 0.5 fs.
I don't know that this has been investigated totally thoroughly, since most people at least put constraints on H-bonds and run at 2 fs (or even repartition the hydrogen mass and run at 4 fs/step). Maybe a recommendation that 0.5 fs is safer, but it's not clear there are definite problems with 1.0 fs in practice (I've seen some analyses that suggest that 0.8 fs may be good enough for all intents and purposes - I think Bill Swope suggested this, but I may be misremembering).
My general understanding of stiff integrators is that if they are numerically stable, the answers are generally close to the right answer.
I agree with the points above -- basically 2 fs with constraints on hbonds, 4 fs with HMR, and 0.5-1 fs without constraints of hbonds, plus the rule of thumb.
Do we need to add this explicitly somewhere?
I am happy to add this recommendation but without a citation?
I had a chat with Ben Leimkuhler about this who has been looking at integrator stability alot. 1 fs is a good recommendation in terms of stability threshold, but accuracy can already be bad at 1 fs. It really depends on the type of integrator used and the implementation of it in the software package from my understanding of our discussion. So this may be more of a 'this is something you might want to look out for' setting. Maybe the recommendation is: Look up what the user manual for the MD code you are using recommends in terms of a time step? (This is based on the naive assumption that this information could be found in a user manual).
but accuracy can already be bad at 1 fs
And I think there is a matter of criteria here. Detectable variation from the "true" answer may be well below the limits of the statistical uncertainty of a free energy calculation. Maybe ask Ben what he means by "bad" in this case. If it's causing errors of 0.2 kcal/mol in a binding free energy calculation, that's a problem. If it's causing an error of 0.01 kcal/mol in both ends of a free energy calculation that end up cancelling, it's not.
I am not sure anyone has looked at this for alchemical type calculations in terms of how accuracy is affected. Ben certainly has not looked at this, but maybe this might something that warrants further investigation. The deviations beyond 1 fs for a water box they looked at are quite large, but then we are asking about small deviations around 1 fs timesteps. Here is the reference: https://royalsocietypublishing.org/doi/full/10.1098/rspa.2016.0138#d3e3484
Hmm. It looks at 2 fs like maybe 3-4 kcal/mol for 216 water molecules, with a total value of ~2700 kcal/mol, so 0.1% relative error. Exactly how that propagates into a free energy calculation is unknown, but I expect, that it would be more like 0.1% error for a free energy calculation; which if it was 10 kcal/mol, that would be 0.01 kcal/mol, which would be acceptable in most cases.
But that last study was with constraints, I think.
I wrote a response in the "response" file, but realized I didn't know where this should go in the main document. There's not an obvious place for "general simulation settings that could affect free energy calculations". Should there be?
I'll do something about this. Thank you for the current input!
I added a small whole new section now and have hopefully addressed this sufficiently.
Thoughts and opinions on this?
Constraints and relative free energy calculations: for me, a timestep of 1 fs is more an upper limit for MD simulations involving non-constrained bonds, I would rather suggest a value of 0.5 fs.