Added option for hydrogen mass repartitioning

[peastman]

There's now an option to enable HMR. Based on openmm/openmm#2864, I have it enabled by default with a hydrogen mass of 1.5 amu. @jchodera does that seem reasonable?

[jchodera]

@peastman: Published methods use 3*amu:

• AMBER: https://pubs.acs.org/doi/abs/10.1021/ct5010406
• Membrane simulations: https://pubs.acs.org/doi/abs/10.1021/acs.jctc.9b00160

In addition, our in-house data suggests masses near 3 to be optimal:

I know of no precedent for a hydrogen mass of 1.5*amu, so I would suggest 3*amu by default (to match the literature reports) unless we want to make a research project out of this.

[peastman]

Using that large a mass significantly slows down sampling. We found that when selecting parameters before, remember? openmm/openmm#2520 (comment) My goal here is to use the very smallest possible mass that avoids instability.

Also, I based this on the data you posted before, which showed a very different curve: openmm/openmm#2520 (comment). On that graph the stability peak is at 1.5.

[jchodera]

I'm not sure that data is compelling enough to break from literature precedent, nor is the difference between 1.5 and 3 amu that large in the other hydrogen mass scan.

Let's look at it this way: if we go with 3amu, we just need to cite a couple of papers. If we go with 1.5amu, one of us is going to have to write some documentation explaining the experiment and findings that led to this choice. The first path is easier, but the second path is also practical but more work (though it would be better if we could use thresholds to examine the stability more systematically across a variety of systems).

[peastman]

It isn't an original result. It's well known that HMR slows down conformational search. See, for example, https://doi.org/10.1002/(SICI)1096-987X(199906)20:8%3C786::AID-JCC5%3E3.0.CO;2-B which reached a similar conclusion.

nor is the difference between 1.5 and 3 amu that large in the other hydrogen mass scan.

Remember that "doing nothing" means a mass of 1. We're talking about the difference between increasing it by 50% and tripling it. That's a big difference. Also consider the effect on the heavy atom it's bonded to. Consider a carbon with three hydrogens. In one case, we decrease its mass from 12 to 10.5. In the other, we decrease its mass from 12 to 6. That's a really big difference.

Given the known problems with HMR, I think we need to follow the principle of, "first do no harm". That means our options are either disable HMR by default, or increase the hydrogen mass by only a small amount by default. Increasing it from 1 to 1.5 showed a big benefit in your lysosyme data, and also was enough to fix the problems in the ApoA1 benchmark, so I think we could justify that as a default setting. Increasing it any more than that would require strong justification.

[jchodera]

It isn't an original result. It's well known that HMR slows down conformational search. See, for example, https://doi.org/10.1002/(SICI)1096-987X(199906)20:8%3C786::AID-JCC5%3E3.0.CO;2-B which reached a similar conclusion.

Hm, but that paper doesn't say anything about 1.5 amu---it only evaluates 1 amu and 4 amu.

Remember that "doing nothing" means a mass of 1. We're talking about the difference between increasing it by 50% and tripling it. That's a big difference. Also consider the effect on the heavy atom it's bonded to. Consider a carbon with three hydrogens. In one case, we decrease its mass from 12 to 10.5. In the other, we decrease its mass from 12 to 6. That's a really big difference.

Given the known problems with HMR, I think we need to follow the principle of, "first do no harm". That means our options are either disable HMR by default, or increase the hydrogen mass by only a small amount by default. Increasing it from 1 to 1.5 showed a big benefit in your lysosyme data, and also was enough to fix the problems in the ApoA1 benchmark, so I think we could justify that as a default setting. Increasing it any more than that would require strong justification.

I'm with you on the concept, but I think this requires we do a real systematic study that is thoroughly documented in the manual (and is reproducible by others) if we are going to go this route, since there is no literature justification for this specific choice of 1.5*amu yet. It may be the right thing to do, but it's certainly going to be more work.

[peastman]

Clearly we don't have the time to do that in the next few days before releasing 7.5. So you're arguing HMR should be off by default? I could agree to that, but only if we also reduce the default time step to 3.5 fs, since we have concrete cases where 4 fs is unstable without HMR.

I don't understand your argument though. We know that HMR with 1.5 amu is more stable than no HMR. We have lots of evidence on that, and I don't think you're questioning it. We also know that 1.5 amu will have less impact on dynamics than pushing it all the way to 3 amu. I don't think you're questioning that either. So what do you see as the problem with using that value? Certainly, we don't know whether it's the optimal value. That would take more study to determine. But it clearly seems to strike a reasonable balance between the two goals of not making simulations blow up and not impacting dynamics and sampling any more than necessary.

[jchodera]

I'm thinking about the poor students who are going to have to write the Methods section of a paper. What will they say to justify why the chose 1.5 amu for HMR if there is no literature precedent? We just need to put something in the manual that will satisfy reviewers, otherwise we're doing them a disservice (even if this is the best choice for the science).

We can still possibly come up with something in the next few days---openmmtools has a variety of test systems, or we could even just use the benchmark systems in the OpenMM benchmarks tool. Could we just run the thresholds tool on those systems and toss the data in an appendix in the manual to show that 1.5 amu helps improve the maximum stable timestep for BAOAB over 1 amu, to give them something to cite?

[peastman]

Sure. Can you run it on the DHFR and ApoA1 systems from the benchmark script? That's 5dfr_solv-cube_equil.pdb and apoa1.pdb in the examples directory.

[peastman]

Are you still planning to do these tests? We need to get this merged soon.

[peastman]

Merging this, since we need it for the release.

[jchodera]

I asked around, but nobody in the lab had time to run this, and I'm seriously backlogged.
We should make sure to do this and document soon. Perhaps you could figure out how to run it?