Address change in padding behavior for addSolvent

[jchodera]

#3480 significantly changed the behavior of addSolvent, which will currently require all users to update their code to avoid issues like choderalab/perses#949

While it may be possible to address this with an extensive plan to update the user guide recommendations, be clear about how this change may break existing code in the release nodes, and add more useful warnings on what likely went wrong when MonteCarloBarostat raises an exception, there may be a simpler solution:

We can retain old behavior for boxShape=None by default, but allow the new behavior to be used when 'cube', 'dodecahedron', and 'octahedron' are specified.

This would ensure we don't break any existing code but use the newer method when users specify a particular shape.

We should also update the user guide with information on these options.

[jchodera]

@peastman : What do you think of this proposed minor change that would enable us to make sure existing code doesn't break but add these new capabilities for people who want to use them?

[peastman]

The new behavior is better in every way. The old behavior used an ad hoc heuristic that often produced unnecessarily large boxes yet also didn't guarantee your box was large enough. Most users don't need to change their code. Now it just does what they probably assumed it did all along.

The problems with the perses test cases are because they're very atypical situations. They build tiny water boxes around molecules that are smaller than the cutoff distance in every dimension.

[jchodera]

@peastman : We can change the recommendations, but I'm worried that this sudden change will literally will break a lot of existing code. It certainly breaks ours.

We can issue a warning about the default option to suggest folks change, and that the old default may be changed in future releases, but not giving folks a heads-up is going to cause users/developers a lot of pain.

[jchodera]

I'll quote from @peastman just a few days ago:

If code that used to work starts throwing exceptions, users are going to start screaming at me.

I do think we should adopt a consistent philosophy of graceful deprecation when possible:

• Introduce the new options, but make sure valid code still runs
• Introduce a warning and update user guide recommendations
• Later, after sufficient warning, deprecate or update the old default behavior

[peastman]

I'll quote from @peastman #3501 (comment):

You left out the sentence that came immediately after that. :) In that case, you were the one arguing for making a change that would cause existing code to start throwing exceptions. I was pointing out the contradiction! But in that case, it would have happened in a very common situation affecting lots of users, rather than an obscure situation affecting very few people. Also, in that situation it would have created a new default behavior that was worse than the previous default, unlike this situation where the new default is better than the previous one.

I do think we should adopt a consistent philosophy of graceful deprecation when possible:

I thought about all of these issues very carefully before changing the behavior, including soliciting feedback about the idea. (You replied that it sounded like "a great approach".) There are complicated tradeoffs in any change. Don't break things without a good reason, but also don't be afraid to make changes when it's clearly for the better. In this case, it's unlikely many users actually knew what the exact box size calculation was, and the new behavior is probably closer to what they assumed it was doing all along.

[jchodera]

If we're committing to breaking code in the wild, what's our plan for engaging with users to help them with transitioning their code to this new approach so they aren't surprised? There's no easy upgrade path where they can prepare for the new code because the optional arguments are not available in earlier OpenMM versions. If they just increase their box padding, things will be much slower (perhaps unusably so) with existing OpenMM versions. Can we at least design this in a way to give them an easy onramp?

[peastman]

They don't need to change their code. That's the point. Aside from some edge cases like the perses unit tests, which mostly aren't representative of real world use, existing code will continue to work as it is. To the extent that the behavior changes, it's strictly an improvement. If the amount of water added increases, that means it really wasn't adding enough water before and they weren't getting as much padding as they asked for.

[jchodera]

@peastman : The perses test failure is simply the canary in the coal mine: The test adds 9A of padding around small molecule solutes, which has been standard practice for production calculations for over a decade. This test is representative of real world use.

Previously, this resulted in a system that was large enough to ensure that barostat changes wouldn't result in a system that shrunk slightly below the cutoff. This is no longer going to be the case for a lot of production code, and we need to provide developers a smoother path to adopting the new modes provided without just breaking their code without explanation.

[peastman]

The problem in those tests was not that 9A of padding is insufficient in general. The problem is that they're applying it to tiny molecules, getting a water box that is too small to simulate with the cutoff they specified (which must be less than half the box width). For a protein, or pretty much anything larger than the tiny molecules used in the tests, that wouldn't be a problem.

This is actually a good illustration of one of the problems with the old method: it didn't do what you thought it did. You asked for 9A of padding, but actually got more and didn't realize it. The test code was broken, but passed anyway for the wrong reason. In other cases you would get less than you asked for, and again wouldn't realize it. The new method guarantees you get the actual amount of padding you asked for, not more and not less.

[peastman]

Let's look at a concrete example. I think that will make it much clearer than just words. Let's create a very simple system, a straight rod that's 3 nm long. (In the example I represent it as two separate atoms at the endpoints, but that's just to keep the code simple.) We'll tell Modeller to create a water box around it with 0.9 nm of padding. And we'll do it twice for two initial orientations of the rod, either along the x axis or along the x-y-z diagonal.

from openmm import *
from openmm.app import *
from openmm.unit import *
import math

ff = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
top = Topology()
chain = top.addChain()
res = top.addResidue('Na', chain)
top.addAtom('Na', element.sodium, res)
res = top.addResidue('Cl', chain)
top.addAtom('Cl', element.chlorine, res)

pos = [Vec3(0, 0, 0), Vec3(3, 0, 0)]
modeller = Modeller(top, pos)
modeller.addSolvent(ff, padding=0.9*nanometers)
print(modeller.topology.getUnitCellDimensions())

pos = [Vec3(0, 0, 0), 3*Vec3(1, 1, 1)/math.sqrt(3.0)]
modeller = Modeller(top, pos)
modeller.addSolvent(ff, padding=0.9*nanometers)
print(modeller.topology.getUnitCellDimensions())

In this case it's obvious what the right answer is: the box should be 3.9 nm on all sides. If it's smaller than that, the rod can rotate so that the endpoints are less then 0.9 nm apart. That would be bad and could lead to simulation artifacts. If the box is larger than that, the user is getting a bigger box than they asked for, which makes their simulation run slower than it needs to. And the box size shouldn't depend on the initial orientation of the molecule. It can freely rotate during the simulation.

Using the new method, we get exactly what we expected:

Vec3(x=3.9, y=3.9, z=3.9) nm
Vec3(x=3.9000000000000004, y=3.9000000000000004, z=3.9000000000000004) nm

But with the old version we get something quite different:

Vec3(x=4.8, y=4.8, z=4.8) nm
Vec3(x=3.5320508075688775, y=3.5320508075688775, z=3.5320508075688775) nm

When we build the box for the axis aligned orientation, it's much bigger than required. And when we build it for the diagonal orientation, it's smaller than required. The endpoints can come within 0.53 nm of each other!

It gets worse as the molecule gets longer. For a 5 nm rod, it's actually possible to get a "padded" water box that is less than 5 nm on a side, allowing copies of the rod to directly contact each other!

[jchodera]

@peastman: I'm 100% onboard with the idea that the new code is superior. My concerns are instead: 1. How do we help developers whose code will break do the right thing to fix "broken code" to help them use the improved methods. Their code (like ours) may have worked before, but will break. What can we do to help them address this, rather than just pin to old OpenMM versions (like many of our users do when we break their code) 2. How do we help them adjust their code now in advance in a way that will work with both 7.7 and the next release? The approach of pushing the next release and forcing everybody's code to break is not appealing---there should be a better path. I don't think any of your responses have explicitly addressed either of these aspects. Let's assume we all agree that the new method is great (we do!). How do we address 1 and 2 above?

[peastman]

The approach of pushing the next release and forcing everybody's code to break is not appealing

Why do you say everybody's code will break? Most code will not break. It will work better than before, with no modification. In fact, I'd say that any code that doesn't work with the new method is already broken, and they just don't realize it. They will now get the amount of padding they asked for instead of a different amount than they asked for.

But in terms of how we communicate it, we will of course

• Discuss it in the release notes
• Point it out in the release announcement
• Describe the new behavior in the API documentation

What else do you suggest?

[peastman]

Based on our offline conversation earlier today, I have a hypothesis about the continuing problems with perses tests. @jchodera said that the padding distance in the tests has been increased from 0.9 nm to 2.0 nm. This caused it to become many times slower, but it still fails with the same error.

Here's a guess: is there a single variable that sets both the padding distance and the cutoff distance? While increasing the padding distance, did the cutoff distance also inadvertently get increased to 2.0 nm?

[jchodera]

Based on our offline conversation earlier today, I have a hypothesis about the continuing problems with perses tests. @jchodera said that the padding distance in the tests has been increased from 0.9 nm to 2.0 nm. This caused it to become many times slower, but it still fails with the same error.
Here's a guess: is there a single variable that sets both the padding distance and the cutoff distance? While increasing the padding distance, did the cutoff distance also inadvertently get increased to 2.0 nm?

@mikemhenry is investigating this in choderalab/perses#953, but we suspect this arises from using a different strategy to call Modeller.addSolvent in some tests that were missed in the first-pass attempt to enlarge padding. @mikemhenry can comment further.

@mikemhenry : See #3537, which ensures the box size is never less than 2*padding, which should make tests and simulations with small box sizes much more robust to failure. This may obviate the need for choderalab/perses#953 once merged, though we will likely want to do some systematic studies of how padding impacts computed binding free energies once the new OpenMM is released.

[mikemhenry]

though we will likely want to do some systematic studies of how padding impacts computed binding free energies once the new OpenMM is released.

This is where the work that @ijpulidos has put in will really pay off :) We can run our benchmarks with different padding and check how things are impacted energy and performance wise.

I will take another pass on choderalab/perses#953, I know missed a few spots in the tests + some of the slowdown was from a regression that has since been fixed (the bug was introduced right before I started working on the fix and was misdiagnosed, so bad luck)

[ijpulidos]

I run the simulations using a padding of 20 Angstroms instead of 9. And of course this has a significant impact in terms of the performance. While the free energy calculations remain the same, the performance is half for the complex phase. Edges that took 1-1.5 hours are now taking 2-3 hours using the same GPUs.

[peastman]

20 A padding is way more than you need. And with the latest changes, you may not need to increase it at all.

[mikemhenry]

I was able to fix our tests by going from 9 A to 11 A padding in almost all cases and the tests don't take much longer (some of the noise in CI spawning is larger than test execution)

[peastman]

If it's any slower at all, that means you're building a bigger water box than you were before. You can decrease the padding further.

[jchodera]

The fixed behavior is almost there, but we're encountering issues with it when solvating simple dipeptides.

Here's an example of a simple example that rapidly fails after just a few steps:

"""                                                                                                                                                                                                                                
Example illustrating failure of new padding algorithm for Modeller.addSolvent()                                                                                                                                                    
"""

# Create alanine dipeptide in vacuum                                                                                                                                                                                               
from openmmtools import testsystems
t = testsystems.AlanineDipeptideVacuum()
topology = t.topology
positions = t.positions

# Create a forcefield                                                                                                                                                                                                              
from openmm import app
forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3p.xml')

# Solvate the peptide                                                                                                                                                                                                              
from openmm import unit
padding = 9.0 * unit.angstroms
modeller = app.Modeller(topology, positions)
modeller.addSolvent(forcefield, padding=padding)
topology = modeller.topology
positions = modeller.positions

# Create System                                                                                                                                                                                                                    
system = forcefield.createSystem(topology, nonbondedMethod=app.PME)

# Add a barostat                                                                                                                                                                                                                   
import openmm
pressure = 1.0 * unit.atmospheres
temperature = 300 * unit.kelvin
barostat = openmm.MonteCarloBarostat(pressure, temperature)
system.addForce(barostat)

# Simulate                                                                                                                                                                                                                         
n_iterations = 1000
n_steps_per_iteration = 250
timestep = 4.0 * unit.femtoseconds
collision_rate = 1.0 / unit.picoseconds
integrator = openmm.LangevinMiddleIntegrator(temperature, collision_rate, timestep)
context = openmm.Context(system, integrator)
context.setPositions(positions)
for iteration in range(n_iterations):
    print(iteration)
    integrator.step(n_steps_per_iteration)

Here's the error:

$ python padding-failure.py 
Traceback (most recent call last):
  File "/lila/home/chodera/buffer-example/padding-failure.py", line 39, in <module>
    context = openmm.Context(system, integrator)
  File "/lila/home/chodera/miniconda/envs/openmm-dev/lib/python3.9/site-packages/openmm/openmm.py", line 16230, in __init__
    _openmm.Context_swiginit(self, _openmm.new_Context(*args))
openmm.OpenMMException: NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.

[peastman]

What change do you suggest?

It appears to me that addSolvent() is working as intended. It creates a box that's big enough. The problem comes when you use a barostat to apply pressure to the system. That causes the box to shrink so it's no longer big enough.

You could increase the padding. But for solvating very small molecules like this, the boxSize argument is probably more appropriate than padding.

Or we could allow you to specify both boxSize and padding, and it would use whichever produced a larger box?

[jchodera]

What change do you suggest?

The philosophy behind our compromise was that we would alter the padding calculation behavior for small solutes so that things that used to "just work" with the current released OpenMM behavior would continue to work with the new behavior, but for large solutes, the more efficient, better approach to computing the padding would be used.

I think we might just need to increase the minimum here.

The default cutoff is 10A, and so the minimum box size for padding=9A is 18A now, so it's easy to still shrink just below twice the cutoff.

What if we made our minimum size 2.5 times the padding? That should take care of these issues.

[jchodera]

@peastman : Can we just try increasing the minimum until this test robustly passes?

[peastman]

I really think that's the wrong approach. The test is making invalid assumptions. The solution is to fix the test, not to make the code conform to the invalid assumptions.

[jchodera]

@peastman: You're missing the point here. The test was making assumptions that were perfectly valid for how OpenMM defined padding for the previous decade, but are now invalid with the new definition of padding for the unreleased OpenMM. Remember that we agreed that we could address both the redefinition and provide users with a much better definition of padding for larger systems as a compromise for accepting this redefinition of padding without an upgrade path.

[peastman]

We had an offline meeting about this and decided on the following plan.

1. I'll create an FAQ entry discussing this error message and the things that can cause it.
2. We'll modify the message for this exception (and probably some others as well) to include a link to the FAQ.

[jchodera]

Thanks!