BAOAB integrator

[peastman]

Following up on https://github.com/openmm/openmm/issues/2386#issuecomment-531438507, I've been looking into creating a new BAOAB integrator. As a starting point, I wanted to get a sense for how well LangevinIntegrator (which uses the Langevin Leapfrog algorithm of https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4308582/) compares, and to consider different ways of formulating it. So I tried to reproduce the tests described in https://royalsocietypublishing.org/doi/full/10.1098/rspa.2016.0138. I simulated alanine dipeptide in a box of water, using rigidWater=True and constraints=HBonds. I tried various timesteps and three different integrators:

• LangevinIntegrator
• BAOABIntegrator from openmmtools
• GeodesicBAOABIntegrator from openmmtools with K_r=2.

For each one I computed the average potential energy and temperature over a set of simulations.

The first thing I found was that there was no significant difference in stability between them. All the integrators were stable at 4 fs. None was stable at 6 fs.

Here is the average potential energy as a function of step. (Ignore the "simplified" line. I'll get to that shortly.)

The BAOAB integrators do much better than LangevinIntegrator as you increase the time step. An interesting difference is that the average energy increase with LangevinIntegrator (the system is too hot) but decreases with BAOAB (the system is too cold). The difference between the two BAOAB integrators is very small, but GeodesicBAOABIntegrator seems to be slightly worse than BAOABIntegrator.

Now here is the average temperature.

The situation is reversed. LangevinIntegrator has no trouble holding the temperature right at 300K. The BAOAB integrators show significant deviation even at 2 fs, and by 4 fs the temperature is about 9 degrees too low.

How should I interpret this? The paper on BAOAB said it was more accurate for configurational quantities than kinetic ones, but I wasn't expecting it to be quite that inaccurate. Should BAOAB be avoided for simulations where you want to compute kinetic quantities?

Next I wanted to consider the best way of implementing it. Contrary to what the names would make you think, BAOABIntegrator and GeodesicBAOABIntegrator are both geodesic integrators. The only difference is that BAOABIntegrator hardcodes K_r=1. It applies constraints seven times in each step, compared to only once in LangevinIntegrator, which makes it much slower in some situations. I've said before that it doesn't actually need to do it that many times. So I removed four of the five calls to addConstrainVelocities(). That is the "simplified" integrator shown in the graphs above. As you can see, the results are totally indistinguishable from the openmmtools integrator.

Of the four constraint applications I removed, two are totally extraneous and have no effect at all on the algorithm. The other two cause the position updates to not move quite as far away from the constraint manifold, which in principle could affect the following step that projects it back onto the manifold. But in practice that effect seems to be negligible.

[maxentile]

Thanks for looking into this @peastman ! Excited at the prospect of seeing an optimized implementation of this integrator in OpenMM.

The first thing I found was that there was no significant difference in stability between them. All the integrators were stable at 4 fs. None was stable at 6 fs.

I've not detected a major difference in stability between these either. The stability at longer timesteps is mainly claimed / observed when including multiple timestepping with solvent solute splitting.

Here is the average potential energy as a function of step.

Nice! This is consistent with our observations.

Now here is the average temperature. The situation is reversed. LangevinIntegrator has no trouble holding the temperature right at 300K.

There are different ways to measure temperature.

I think you'd find that although the average kinetic energy appears accurate, the configurational temperature would be corrupted under the default Langevin integrator (note that in your previous plot, the default Langevin integrator samples higher-energy configurations on average at large timesteps).

@c-matthews may be able to comment more on this.

How should I interpret this? The paper on BAOAB said it was more accurate for configurational quantities than kinetic ones, but I wasn't expecting it to be quite that inaccurate. Should BAOAB be avoided for simulations where you want to compute kinetic quantities?

This is consistent with our observations.

Maybe @c-matthews could comment on how to interpret this -- I'll initially note that getting the marginal velocity distribution wrong doesn't necessarily mean that long-timescale kinetic properties will be wrong.

However, kinetic properties from a Langevin simulation under any discretization should probably be interpreted with caution, and will depend sensitively on the collision rate.

Next I wanted to consider the best way of implementing it. Contrary to what the names would make you think, BAOABIntegrator and GeodesicBAOABIntegrator are both geodesic integrators. The only difference is that BAOABIntegrator hardcodes K_r=1.

Sorry, that was confusing naming on my part!

The BAOABIntegrator in OpenMMTools is indeed an implementation of g-BAOAB whenever the system contains constraints. If the system does not contain constraints, it is just BAOAB.

It applies constraints seven times in each step, compared to only once in LangevinIntegrator, which makes it much slower in some situations. I've said before that it doesn't actually need to do it that many times. So I removed four of the five calls to addConstrainVelocities(). That is the "simplified" integrator shown in the graphs above. As you can see, the results are totally indistinguishable from the openmmtools integrator. Of the four constraint applications I removed, two are totally extraneous and have no effect at all on the algorithm.

Thanks! We should also adopt these simplifications. We had erred on the side of including extraneous constraint calls in the initial implementation to be extra conservative with regard to supporting different user-specified Langevin splitting strings, but should have removed these.

The other two cause the position updates to not move quite as far away from the constraint manifold, which in principle could affect the following step that projects it back onto the manifold. But in practice that effect seems to be negligible.

Probably another good point for @c-matthews to potentially chime in on!

[jchodera]

@peastman: For a deeper explanation of what is going on when some integrators preserve configurational properties even at large timesteps (such as BAOAB) while others seem preserve the kinetic energy distribution (like OBABO), see Josh's excellent paper that examines the ideas of @c-matthews in biomolecular systems.

The key idea is illustrated here in plotting the errors in configuration space and phase space this quartic potential: BAOAB (called VRORV in our paper) lines up the errors to cancel when you marginalize over velocities. That's fine, since nearly everyone is interested in configuration space properties---the velocity distribution is analytically tractable, after all! It may corrupt some kinetics, but Langevin will do that anyway as you increase the collision rate.

There's lots of disagreement over exactly which integrator to use for more complex kinetic properties like time-correlation functions, but it's very clear that we should be recommending people use BAOAB whenever they're interested in configurational properties, free energies, etc. (which is much of the time).

[peastman]

I'm already visualizing the flood of posts I'll have to deal with saying, "My temperature is coming out too low, what's wrong with my simulation?" (Or more likely, "What's wrong with your code?") I'm not sure how we can deal with that, but we definitely need to be ready for it!

[peastman]

We should also adopt these simplifications.

All you need to do is remove all the calls to addConstrainVelocities() except the last one. You can think of this as closely analogous to RATTLE: instead of requiring the velocity constraints to be satisfied at every intermediate point within a step, you only require them to be satisfied at the end of the step.

My tests were only for the basic BAOABIntegrator. With the geodesic integrator with more substeps, it might make more of a difference. That will need testing. But at the very least whenever you have a series of steps of the form

v = ... constrain velocities v = ... constrain velocities

you can remove the first constraint call. Projecting velocities onto the constraint manifold is a linear operation, so the result will be identical up to numeric precision.

One other observation. BAOAB is not just a terrible name to try to figure out how to pronounce. It's also a very awkward one to type. I wonder if we can think of a different name...

[jchodera]

I'm already visualizing the flood of posts I'll have to deal with saying, "My temperature is coming out too low, what's wrong with my simulation?" (Or more likely, "What's wrong with your code?") I'm not sure how we can deal with that, but we definitely need to be ready for it!

We should educate people on how integrators work and what they should care about. We can do that by pointing them to clear, easy to understand material and literature. I'm happy to write a section for the manual!

[peastman]

We need to pick a name for this. As far as I can tell, the only MD package that currently supports this method is Amber. They refer to it as the "middle" method, following the naming from https://arxiv.org/abs/1707.04685. Also, they use a leapfrog version of the method, so the actual name they use in the documentation is LFMiddle.

I also checked the documentation for CHARMM, NAMD, Gromacs, LAMMPS, and Tinker. None of them appears to support it.

Here are some possible names:

• BAOABIntegrator (+) Matches the name in openmmtools (+) Accurately describes what it does (-) Most users won't know what it means (-) Awkward to type (-) No one knows how to pronounce it

• BAOABLangevinIntegrator (+) Slightly more informative (-) Still awkward to type and unpronounceable

• LangevinIntegrator2 (+) Makes it clear this is roughly equivalent to LangevinIntegrator (-) Gives no information about how it's different

• LangevinMiddleIntegrator (+) Same advantages as LangevinIntegrator2, but slightly more informative (+) Matches terminology used in Amber

• LMLangevinIntegrator (LM=Leimkuhler-Matthews) (+) Another slightly more informative alternative to LangevinIntegrator2

• Something else?

[peastman]

Any thoughts on this?

[jchodera]

I think BAOABIntegrator is the cleanest choice for now. It clearly and unambiguously links the integrator to the literature.

I'd honestly like to see this replace LangevinIntegrator, but this would create problems with reproducibility. I think it would be fine to do some renaming at a major release (which allows API breaks in standard semantic versioning schemes), but we argue about that later.

Have we updated the user guide as well? Do you want me to provide some text for this with appropriate citations?

[peastman]

It clearly and unambiguously links the integrator to the literature.

In the literature this method variously gets referred to a BAOAB, VRORV, and "middle", so it's no more unambiguous than some of the other choices.

I think it would be fine to do some renaming at a major release (which allows API breaks in standard semantic versioning schemes), but we argue about that later.

We're not just talking about breaking API compatibility (which is something I try hard not to do). We're talking about causing existing scripts to silently start using a different algorithm that produces significantly different results, and that may violate assumptions made in those scripts. I would never do that under pretty much any circumstances.

Have we updated the user guide as well?

I haven't even started implementing it yet.

[jchodera]

In the literature this method variously gets referred to a BAOAB, VRORV, and "middle", so it's no more unambiguous than some of the other choices.

From the authors (Ben Leimkuhler and Charles Matthews), "BAOAB" is the canonical name. VRORV was a name we used for consistency with an earlier paper, but we concur that BAOAB should be the canonical name. "middle" is not uniquely descriptive.

We're not just talking about breaking API compatibility (which is something I try hard not to do). We're talking about causing existing scripts to silently start using a different algorithm that produces significantly different results, and that may violate assumptions made in those scripts. I would never do that under pretty much any circumstances.

We can discuss this in a separate thread. Fundamentally, Langevin dynamics is simply a numerical method for generating samples of positions. (Velocities are technically an intermediate variable.) These different splittings are simply different numerical integration approximation schemes. We make changes to the approximation behavior of other numerical schemes all the time.

Regarding implementation: One of the more useful aspects of the openmmtools version is that it supports alternative splittings (perhaps not so useful other than for research) as well as multiple force groups (which could be very useful in speeding things up). I'm tagging @maxentile, who has been experimenting with this (and solute-solvent splitting) with Charles Matthews, to share thoughts on what might be useful for supporting in an OpenMM native version.

[peastman]

How do you feel about BAOABLangevinIntegrator? It's equally explicit while being more informative.

Fundamentally, Langevin dynamics is simply a numerical method for generating samples of positions.

Yes, but what does that have to do with breaking API compatibility and causing scripts to behave differently?

Another important question is whether we should use the leapfrog or velocity-Verlet version of the algorithm? The former matches the behavior of LangevinIntegrator, and also is the version implemented in Amber. The latter matches the openmmtools version.

[jchodera]

How do you feel about BAOABLangevinIntegrator? It's equally explicit while being more informative.

Sounds great.

Yes, but what does that have to do with breaking API compatibility and causing scripts to behave differently?

The thing below.

Another important question is whether we should use the leapfrog or velocity-Verlet version of the algorithm? The former matches the behavior of LangevinIntegrator, and also is the version implemented in Amber. The latter matches the openmmtools version.

I think it is significantly easier for users to use synchronized velocities and positions, since the half-timestep-lagged version trips people up all the time. (This would be the major API change if we changed the behavior of LangevinIntegrator from lagged to synchronized velocities.)

[maxentile]

I also checked the documentation for CHARMM, NAMD, Gromacs, LAMMPS, and Tinker. None of them appears to support it.

Huh, although it's not mentioned in the documentation(!), it looks like it's in the NAMD source as of ~2012 (see the boolean simulation parameter langevin_useBAOAB).

Regarding implementation: One of the more useful aspects of the openmmtools version is that it supports alternative splittings (perhaps not so useful other than for research) as well as multiple force groups (which could be very useful in speeding things up). I'm tagging @maxentile, who has been experimenting with this (and solute-solvent splitting) with Charles Matthews, to share thoughts on what might be useful for supporting in an OpenMM native version.

I think the way we support multiple splittings and multiple force groups in openmmtools is useful for research, but is probably not mature enough to be used in openmm proper. In ongoing work, we are interested in quantifying the sampling errors introduced by different kinds of force-group splittings.

I think the BAOAB integrator without force-group splitting is a very solid method though, and I am excited to see an optimized implementation included with OpenMM.

How do you feel about BAOABLangevinIntegrator? It's equally explicit while being more informative.

That sounds good to me.

Another important question is whether we should use the leapfrog or velocity-Verlet version of the algorithm? The former matches the behavior of LangevinIntegrator, and also is the version implemented in Amber. The latter matches the openmmtools version.

I'm not sure I follow here. Do you mean changing how the position and velocity updates are computed, or do you mean holding the sequence of computations fixed, but changing what is returned by .getVelocities()?

[peastman]

Do you mean changing how the position and velocity updates are computed, or do you mean holding the sequence of computations fixed, but changing what is returned by .getVelocities()?

Velocity Verlet and Leapfrog both involve exactly the same sequence of calculations. The only difference is where you draw the lines between time steps. If you take the velocity half step at the end of one time step and merge it with the one at the beginning of the next time step, you get a leapfrog algorithm.

[peastman]

Yes, but what does that have to do with breaking API compatibility and causing scripts to behave differently?

The thing below.

Exactly my point! You can't just swap one for the other and expect people's code to work, even if they're both "Langevin dynamics".

[jchodera]

Velocity Verlet and Leapfrog both involve exactly the same sequence of calculations. The only difference is where you draw the lines between time steps. If you take the velocity half step at the end of one time step and merge it with the one at the beginning of the next time step, you get a leapfrog algorithm.

It's so much easier for all users (new and old) to deal with synchronized (x,v) pairs without having to worry about correcting for velocities that are halfway through the integration cycle.

What about adding a VelocityVerletIntegrator at the same time as the BAOABLangevinIntegrator? Both could use synchronized (x,v).

[peastman]

What about adding a VelocityVerletIntegrator at the same time as the BAOABLangevinIntegrator? Both could use synchronized (x,v).

That's certainly an option. The one disadvantage is that it's a bit slower. Leapfrog integrators only need to apply constraints once per step, but velocity Verlet has to do it twice.

[jchodera]

If you're requesting the integrator run N steps of velocity Verlet, could you skip the constraint at the end of steps 1...(N-1) under the hood? That would have the same effect.

That would just be slower for integrator.step(1) calls.

[peastman]

Yes, that would be possible, but it would mean that repeated calls to step(1) would produce a different result from a single call to take many steps. SHAKE and RATTLE are both accurate to the same order, but they aren't the same.

[jchodera]

Does the API explicitly guarantee that these will produce numerically identical results?

[peastman]

I think I didn't explicitly state it only because it seemed so obvious, it didn't occur to me it would ever need to be stated. That would certainly surprise a lot of users if it weren't true.

[peastman]

Implemented by #2436.

[peastman]

I ran some tests to see how much impact the extra constraint evaluations have on performance. Here are the numbers for the DHFR benchmark.

CUDA

• LangevinIntegrator: 423
• BAOABLangevinIntegrator: 387

CPU

• LangevinIntegrator: 15.1
• BAOABLangevinIntegrator: 14.6

So it's roughly 10% slower.

[jchodera]

Did you want to explore the idea of allowing the user to skip the extra constraint step when calling integrator.step(N)? Or do you think it's not worth it if this is no longer numerical identical to N calls to step(1)?

[peastman]

I think it's better not to.