Solvent padding description

Sorry for the long explanation - the actual issue here is incredibly pointed and I want to make sure I'm highlighting only it.

OpenFE's docstring of the padding option in solvation settings is:

https://github.com/OpenFreeEnergy/openfe/blob/695d03208bc429fc282e4332ce668b6911111889/openfe/protocols/openmm_utils/omm_settings.py#L53

This setting appears to be passed 1:1 to OpenMM without modification

OpenMM's documentation of a similar argument is a little more verbose

You can give a padding distance. A bounding sphere containing the solute is determined, and the box size is set to (sphere diameter)+(padding). This guarantees no atom in the solute will come closer than the padding distance to any atom of another periodic copy. If the sphere diameter is less than the padding distance, the box size is set to 2*(padding) to ensure no atom is closer than the padding distance to two periodic copies of any other atom.

For some small molecules, however, this produces results as if the bounding sphere representing the solute is size zero, which doesn't intuitively hold up to me. OpenFF has a similar implementation (and documented argument) that goes through Packmol after determining the expected system size, and for similar systems it frequently produces boxes a few Angstroms larger. I'm running small tests with a benzene-in-water system, and OpenFE's code path produces box vectors of 2.4 nm (simply twice the padding value of 1.2) whereas the OpenFF code path produces numbers closer to 2.8 nm (twice the padding value plus a little bit of fuzziness representing the size of the molecule). Maybe for a single atom it makes sense, but for a 3D molecule with some size to it, I don't think the box vectors should be twice the padding value in all three dimensions. Unfortunately this seems to be what OpenMM does (see below).

Both OpenFF and OpenFE code paths produce sensible results (no crashes) as far as I've observed, so I don't think it makes sense to ask for a behavior change. The documentation doesn't communicate that there should be different behavior, so I think that could be updated in some way. I have no reason to believe results would be different, but accidentally packing into a larger-than-needed-box is an easy way to slow down simulations with no benefit ((2.8 / 2.4) ** 3 ~= 1.6, or 60% slower compute!).

Here's a rough demonstration of how the OpenMM code path produces different system sizes with molecules of different size

import openmm
import openmm.app
import openmmforcefields
import openmm.unit

from openff.toolkit import Molecule

from openmmforcefields.generators import GAFFTemplateGenerator


def smiles_to_box_vectors(smiles: str):
    molecule = Molecule.from_smiles(smiles)

    molecule.generate_conformers(n_conformers=1)

    gaff = GAFFTemplateGenerator(molecules=molecule)

    forcefield = openmm.app.ForceField("tip3p.xml")
    forcefield.registerTemplateGenerator(gaff.generator)

    modeller = openmm.app.Modeller(
        molecule.to_topology().to_openmm(),
        molecule.conformers[0].to_openmm(),
    )

    modeller.addSolvent(
        forcefield,
        padding=1.0 * openmm.unit.nanometers,
    )

    a, b, c = modeller.getTopology().getPeriodicBoxVectors()

    return {round(value._value, 3) for value in [a[0], b[1], c[2]]}


for smiles in [
    "N",
    "CCO",
    "c1ccccc1",
    "C12C3C4C1C5C2C3C45",
    "CCCC[N+](CCCC)(CCCC)CCCC",
    20 * "C",
]:
    print(smiles, smiles_to_box_vectors(smiles))

# N {2.0}
# CCO {2.0}
# c1ccccc1 {2.0}
# C12C3C4C1C5C2C3C45 {2.0}
# CCCC[N+](CCCC)(CCCC)CCCC {2.226}
# CCCCCCCCCCCCCCCCCCCC {2.767}

If it would be useful I can

• Run a MWE of a system setup using OpenFE's API (I have this elsewhere, not stripped down, and the behavior is consistent with
• Show what I referred to earlier with OpenFF's Packmol code and why the results are a little different
• Push this up to OpenMM, though updating the documentation there might take longer with slower release cycles

Start with the padding docstring in openfe/protocols/openmm_utils/omm_settings.py at the linked line, then compare its behavior with the pass-through in openfe/protocols/openmm_utils/system_creation.py. Update the description to explain how OpenMM determines the box size for small and larger solutes, and verify that the documented behavior matches the cited OpenMM documentation.