Add truncated octahedral solvation option for Modeller.addSolvent()

[jchodera]

It turns out that it's quite difficult to automate solvation within a truncated octahedron.

Could we add an option to Modeller.addSolvent() that would allow the user to select between standard box shapes?

cc : @glass-w, who can provide more information on use cases.

[peastman]

In what way is it difficult? You just specify appropriate box vectors. Here's the code from OpenMM-Setup that does it.

https://github.com/openmm/openmm-setup/blob/931bddaf9fa185e44c79af2c638c22eed3cbe7d4/openmmsetup/openmmsetup.py#L249-L256

[jchodera]

@glass-w can provide more information about the difficulties here, but I think the bigger picture question is that you're asking "why is it hard for a new user to just do

vectors = mm.Vec3(1,0,0), mm.Vec3(1/3,2*sqrt(2)/3,0), mm.Vec3(-1/3,sqrt(2)/3,sqrt(6)/3)                
boxVectors = [(maxSize+geompadding)*v for v in vectors]

?" then the point is already made. :)

If we can easily encapsulate these lines into Modeller.addSolvent() under the control of an optional box_geometry kwarg, we can make it easier for anyone to just specify a supported geometry and not have to worry about this.

Related to #3113, are there any assumptions in Modeller.addSolvent() that assume a particular box geometry? There are some lines that do not appear to be agnostic.

[peastman]

I was thinking more along the lines of, "Why is it hard for a new user to just select the 'truncated octahedron' option in OpenMM-Setup."

[jchodera]

Because we need to automate it, and OpenMM-Setup is intended to be used as a web app and not a user-accessible library like openmm.app.Modeller is intended.

[zhang-ivy]

I've just tried solvating (with a truncated octahedron) my protein complex with the following code (mostly copied from https://github.com/openmm/openmm-setup/blob/931bddaf9fa185e44c79af2c638c22eed3cbe7d4/openmmsetup/openmmsetup.py#L249-L256):

# Create a modeller
modeller = app.Modeller(topology, positions)

# Solvate     
geompadding = 0.9 * unit.nanometers
maxSize = max(max((pos[i] for pos in positions))-min((pos[i] for pos in positions)) for i in range(3))
vectors = openmm.Vec3(1,0,0), openmm.Vec3(1/3,2*np.sqrt(2)/3,0), openmm.Vec3(-1/3,np.sqrt(2)/3,np.sqrt(6)/3)
boxVectors = [(maxSize+geompadding)*v for v in vectors]
modeller.addSolvent(self.system_generator.forcefield, model=water_model, boxVectors=boxVectors, ionicStrength=ionic_strength)

# Write to a PDB file
solvated_topology = modeller.getTopology()
solvated_positions = modeller.getPositions()
app.PDBFile.writeFile(solvated_topology, solvated_positions, open("trunc_oct_box_vecs_complex.pdb", "w"), keepIds=True)

And when I visualize my system in pymol, the solvent box looks rectangular:

When I use supercell to visualize the unit cells, it looks like the unit cell isn't actually rectangular, but has rhombic nature:

Is this the desired behavior? Or should the solvent box/unit cells actually look like a truncated octahedron?

[peastman]

That's correct. See the last paragraph under http://docs.openmm.org/latest/userguide/theory.html#periodic-boundary-conditions.

[zhang-ivy]

Got it, thanks!

[peastman]

I'm looking into how to implement this. Here is how we currently implement the padding option:

openmm/wrappers/python/openmm/app/modeller.py

Lines 393 to 394 in 4142e19

	3. You can give a padding distance. The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
	box of size (largest dimension)+2*padding is used.

An obvious implementation would be to add another option for box shape. We would build a box of the specified shape with the specified padding distance. However, the way we currently compute box size (largest dimension along x, y, or z) doesn't really make sense for a non-rectangular box. Would it be better to instead find the minimal bounding sphere containing all particles, then set the box width to (sphere diameter)+padding? This could produce boxes either larger or smaller than the current calculation. It would guarantee what you really want: no matter how the molecule rotates, it will never come closer than the specified padding to another copy. The current calculation doesn't guarantee that. And it makes sense for any box shape, not just rectangular ones.

[zhang-ivy]

@peastman : Sounds good to me, tagging @jchodera to see what he thinks.

[jchodera]

This sounds like a great approach. Ideally, we would exclude solvent from the bounding sphere computation (eg crystal waters and salts would not contribute), but this is a minor detail.