Add `Topology.solvate`

[mattwthompson]

Is your feature request related to a problem? Please describe.

People often want to solvate their systems (small molecule, protein, combination thereof), often with water. It would be nice to have a toolkit-level API to make this easy

Describe the solution you'd like

Add Topology.solvate

Describe alternatives you've considered

Work through to OpenMM and possibly convert back

Additional context

@Yoshanuikabundi has lots of thoughts on this, @j-wags does too

[Yoshanuikabundi]

TL;DR

My proposal is

class Topology:     
    def solvate(
        self,
        *,
        box_vectors: Quantity | None = None,
        box_shape: np.ndarray | None = None,
        box_padding: Quantity | None = None,
        box_width: Quantity | None = None,
        box_n_waters: int | None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        """
        Solvate the topology with a pre-equilibrated box of water.
        
        Parameters
        ==========
        box_vectors
            Box vectors to solvate. Mutually exclusive with all other box arguments.
        box_shape
            Shape of the generated box. Mutually exclusive with ``box_vectors``.
        box_padding
            Padding distance between atoms in ``self`` and the edge of the box.
            Mutually exclusive with all other box arguments except ``box_shape``.
        box_width
            Width of the box. Mutually exclusive with all other box arguments
            except ``box_shape``.
        box_n_waters
            Number of waters in the final solvated topology, including any in 
            ``self``. The box will be sized to fit this many waters. Mutually 
            exclusive with all other box arguments except ``box_shape``.
        salt_conc
            The concentration of the specified salt in the final topology. 
        positive_ion
            A monovalent, monoatomic cation to neutralise and salt the topology.
        negative_ion
            A monovalent, monoatomic anion to neutralise and salt the topology.
        neutralize
            ``True`` to add or remove ions via SLTCAP to ensure the solvated
            topology is electronically neutral, ``False`` to guarantee the salt
            concentration is as accurate as possible.
        """

        ...

This allows a fairly thin wrapper around Modeller, but not quite as thin as we discussed (see final section). What do you think?

Background info

I've previously designed two solvation APIs for our PackMol wrapper - one for pack_box, which computes the box from the target density, number of molecules, and a box shape, and the second for solvate_topology, which is more oriented towards solvation of an existing topology rather than simply generating a box and computes the box from a padding distance around the existing topology. solvate_topology still requires a target density, but this is used to compute the number of solvent molecules rather than the size of the box.

Both of these are a bit inappropriate here, as the density is really defined by the pre-equilibrated box and isn't a free parameter.

Desired features for initial implementation

At our solvation power hour we decided that the following features would form a good basis that we all agreed were both valuable to and unlikely to be misused by users:

• Limit solvents to water with a few ions
• Solvate a topology by tiling a pre-equilibrated box and removing clashing waters/ions
• Implement by using OpenMM Modeller's addSolvent method with Sage via OpenMMForcefields' SMIRNOFFTemplateGenerator
• Takes a topology as self and either modifies that topology in place or returns a new one
• Possibly extend to other solvents, user-provided solvent boxes, other methods in the future

We did not specifically discuss box generation, but I think the shared assumption was that this would be a part of this method.

Proposed API designs

I can see two basic paths here:

Design 1: Require the box vectors to already be set on the Topology, and fill that box with solvent.

Pros: Simple, clean separation of concerns, can provide multiple simple methods for producing boxes, simple linters can easily tell user if they've done something wrong (because they will be calling the wrong box generation method with the wrong arguments), very obvious to users how to extend box generation code (just set the box!), autocompletion can tell users what box generation methods are available (....generate_box_with<tab>). What to do with existing box vectors on topology is always obvious. Code is very self-documenting.

Cons: More code for us to write in box generation methods, less one-stop-shop feel in solvate method, slightly more user code, no obvious way to exactly specify number of waters in final box and maintain density

Example user code:

# 0.1 M saline in a dodecahedral box whose least atomic image distance is <solute length> + 11 angstroms
topology.generate_box_with_padding(shape=RHOMBIC_DODECAHEDRON, distance=Quantity(11, "angstrom"))
topology.solvate(ionic_strength=Quantity(0.1, "molar"))

# Pure water in a cubic box with width 4 nm
topology.generate_box_from_width(shape=CUBE, width=Quantity(4, "nanometer"))
topology.solvate()

# Solvate a box with a user-defined shape and size with 1 M cesium iodide
topology.box_vectors = Quantity([[...], [...], [...]], "nanometer")
topology.solvate(
    ionic_strength=Quantity(1, "molar"), 
    positive_ion=Molecule.from_smiles("[Cs+]"),
    negative_ion=Molecule.from_smiles("[I-]"),
)

# Rhombic dodecahedral box just big enough to fit 750 waters around the solute
scipy.optimizer...

# Raise an exception
topology.box_vectors = None
topology.solvate()

Questions: What box generation methods do we provide, if any? Would it suffice to just provide constants RHOMBIC_DODECAHEDRON and CUBE that specify 3x3 box vectors of each shape with unit least atomic image distance, and a Topology.compute_bounding_sphere_diameter() method, and allow the user to compute the box themselves? eg, the generate_box_with_padding call above could be topology.box_vectors = RHOMBIC_DODECAHEDRON * (Quantity(11, "angstrom") + topology.compute_bounding_sphere_diameter())

Design 2: Have a solvate method with a complex signature that provides several mutually exclusive sets of arguments for specifying box generation

Pros: Simple for the user. Pythonic. With careful argument design and appropriate @overloads, type checkers would likely be able to check user input. Small API surface area. If we adopt the (very reasonable) box specification algorithm of the Modeller.addSolvent() method, initial implementation becomes extremely simple.

Cons: Complex function signature with mutually exclusive arguments. Difficult to add new box generation methods in the future. Simple linters likely cannot check validity of user code. Users have to read the docstring to understand what's happening.

Example user code:

# 0.1 M Saline in a dodecahedral box whose least atomic image distance is <solute length> + 11 angstroms
topology.solvate(
    box_padding=Quantity(11, "angstrom"),
    box_shape=RHOMBIC_DODECAHEDRON,
    salt_conc=Quantity(0.1, "molar"),
)

# Pure water in a cubic box with width 4 nm
topology.solvate(
    box_shape=CUBE,
    box_width=Quantity(4, "nanometer"),
)

# Solvate a box with a user-defined shape and size with 1 M cesium iodide
topology.solvate(
    box_vectors=Quantity([[...], [...], [...]], "nanometer"),
    salt_conc=Quantity(1, "molar"), 
    positive_ion=Molecule.from_smiles("[Cs+]"),
    negative_ion=Molecule.from_smiles("[I-]"),
)
# OR
topology.box_vectors = Quantity([[...], [...], [...]], "nanometer")
topology.solvate(
    salt_conc=Quantity(1, "molar"), 
    positive_ion=Molecule.from_smiles("[Cs+]"),
    negative_ion=Molecule.from_smiles("[I-]"),
)

# Rhombic dodecahedral box just big enough to fit 750 waters around the solute
topology.box_vectors = None
topology.solvate(
    box_shape=RHOMBIC_DODECAHEDRON,
    box_n_waters=750,
)

# Raise an exception
topology.box_vectors = None
topology.solvate()

Questions: What do we do when the topology has box vectors that conflict with the requested box parameters? Do we have a default box generation method when no parameters are provided? What box parameters should be include initially and how do we validate and apply them?

My proposal

I'm really drawn to design 1, with the significant caveat that I don't know how I'd implement n_waters (numAdded in the Modeller method). If anyone can think of a way to do that cleanly, I'd love to see it, but otherwise I'd propose design 2. To implement SLTCAP, we'd need to use modeller to produce a box of pure water, and then replace the appropriate number of waters with ions ourselves. I'd propose we use a signature like this:

class Topology:
    def solvate(
        self,
        *,
        box_vectors: Quantity | None = None,
        box_shape: np.ndarray | None = None,
        box_padding: Quantity | None = None,
        box_width: Quantity | None = None,
        box_n_waters: int | None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ):
        ...

With overloads and comments describing the mutually exclusive arguments, this looks like:

class Topology:
    # If you specify box_vectors, none of the other box arguments may be
    # specified
    @overload
    def solvate(
        self,
        *,
        box_vectors: Quantity,
        box_shape: None = None,
        box_padding: None = None,
        box_width: None = None,
        box_n_waters: None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...
        
    # If you specify box_padding, you must also specify box_shape, and none of
    # the other box arguments may be specified
    @overload
    def solvate(
        self,
        *,
        box_shape: Quantity,
        box_padding: Quantity,
        box_vectors: None = None,
        box_width: None = None,
        box_n_waters: None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...
        
    # If you specify box_width, you must also specify box_shape, and none of
    # the other box arguments may be specified
    @overload
    def solvate(
        self,
        *,
        box_shape: Quantity,
        box_width: Quantity,
        box_vectors: None = None,
        box_padding: None = None,
        box_n_waters: None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...
        
    # If you specify n_waters, you must also specify box_shape, and none of
    # the other box arguments may be specified. box_n_waters is the final total
    # number of waters; we need to subtract the number of waters in the initial
    # topology and the number of ions we want before we pass it to numAdded
    @overload
    def solvate(
        self,
        *,
        box_shape: Quantity,
        box_n_waters: Quantity,
        box_vectors: None = None,
        box_padding: None = None,
        box_width: None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...
        
    # If you specify none of the box arguments, the topology's box vectors are 
    # used
    @overload
    def solvate(
        self,
        *,
        box_vectors: None = None,
        box_shape: None = None,
        box_padding: None = None,
        box_width: None = None,
        box_n_waters: None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...
        
    def solvate(
        self,
        *,
        box_vectors: Quantity | None = None,
        box_shape: np.ndarray | None = None,
        box_padding: Quantity | None = None,
        box_width: Quantity | None = None,
        box_n_waters: int | None = None,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
        neutralize: bool = True,
    ) -> None:
        ...

Note that this design implies that we do a little more than just a thin wrapper around Modeller.addSolvent():

• We use the SLTCAP method, so we must compute the numbers of each ion and we must add them ourselves by replacing water molecules from the pure water box returned by OpenMM
• If the user specifies a total number of waters, we can't just pass this on to Modeller. First, we subtract the number of waters in the topology so that the final count includes these. Then, we add the number of ions we will add so that these additional waters can be replaced with ions.
• This design supports arbitrary molecules as ions; at least initially, we should restrict this to monatomic monovalent ions by erroring if the provided molecule is not that.

I have considered a few alternatives and details:

• An alternative to accepting molecules as ions would be for the arguments to be positive_ion: Literal['Cs+', 'K+', 'Li+', 'Na+', 'Rb+'] = 'Na+' and negative_ion: Literal['Cs+', 'K+', 'Li+', 'Na+', 'Rb+'] = "Cl-", or for there to be no such arguments and to only support NaCl. Starting with Molecule allows us to extend this to small polyatomics like NH4+ and multivalent ions like Ca2+ without a breaking change.
• We could consider making RHOMBIC_DODECAHEDRON the default value for box_shape. If we do this, we should do it by interpreting box_shape=None as RHOMBIC_DODECAHEDRON when exactly one of the box methods is set, and documenting this clearly (this also changes several of the overloads above). This is so we can clearly interpret the case where the user specifies box_shape=RHOMBIC_DODECAHEDRON but no other box arguments as an error (rather than using the topology's box vectors).
• We could have an initial design that does not use SLTCAP and only accepts the ions supported by Modeller - this would allow Topology.solvate() to be a very thin wrapper indeed. However, Chapin says we should use SLTCAP (at least in protein benchmarks), which appears definitive.
• It would be reasonable to NEVER use the topology's box vectors (and require the user to specify box_vectors=topology.box_vectors if they want them). I'm not sure if this would be more surprising or less, but I can imagine it saving a user if they forget to specify a box argument and have inappropriate box vectors from, say, a PDB file.
• I have no strong preference between modifying the topology in place, or returning a copy that has been solvated. The latter is generally my preference, but the former is more in line with the Toolkit's existing API.
• box_n_waters could be box_n_waters_added, but I think conceptually you generally want to control the total number of waters, not the number added. If I imagine a solute topology with waters that I want to preserve (say, crystallographic waters embedded in a protein), I'm doing it because I expect those water positions to be important in the equilibrium distribution, and so I'm accelerating equilibration by putting them in there to start with. I don't really care, at equilibrium, where those waters came from or if those waters are exchanging with the solvent. I can imagine having multiple crystal structures of the same protein in different configurations with different crystallographic waters that I want to be able to interconvert with REMD or something much more easily than any use case where I only care about the number of waters I added at some obscure step in system prep.
• Along similar line of logic, I think SLTCAP should take ions of the specified species already present in the topology into account when calculating how many ions to add - ie, find all the specific ions already present, subtract their charge from the topology's charge when computing the number of ions to add, then add that many fewer ions than you otherwise would. Not sure what we should do if the existing topology already has more ions of one species or the other than SLTCAP would add; I guess this would trigger a fallback where initial ions are treated as part of the solute.

[mattwthompson]

completely independent of the changes here, I think the toolkit needs better support (documentation and testing) for non-cubic boxes. I'm pretty sure it can be done, but I'm not sure if it's thoroughly tested and I also don't think we have quality documentation to point users to

---

Thanks for writing all of this up! Overall I'm a big fan of these ideas and the value they can provide to users. I think we're already at a point of picking between viable options and details, so (modulo the implementation) I'm happy we're in this state.

I favor option 1

• The overall amount of user code is non-trivial, and will require some handholding. But I like that determining the box vectors and solvating the topology are separate steps.
  • Some topologies will come pre-loaded with box vectors and that step can be skipped (or overriden)
  • This gels with my mental model of steps in a system prep pipeline
• Nice behavior: just calling the method without arguments just works (in the sense that it gives a sensible result with zero cognitive effort of customization)
• I think people come to the toolkit wanting their ligand/protein/etc. to be solvated. I don't think users care deeply if the number of waters is 4095 or 4096. Not sure if/how to deal with translating this to OpenMM's numAdded and a little lost as to whether this is an implementation detail or an important user-facing design question
• Prefer keeping the Topology.make_me_a_new_fancy_box methods as simple as possible starting out, and ideally as few of them as we need for 80% of users. (We can add more methods in later, and getting the last 0.1% of use cases covered should not be a blocker for the first implementation.) If we add nothing here (not necessarily the best idea) users are still able to specify box vectors via the setter, it's just not the easiest.

Thoughts on option 2

• As alluded to above, there's a lot going on inside of the parenthesis. For the simplest use cases it's more than I'd want and for more complex use cases it's ... well, it's a lot to wrangle if I don't want to be a power-user.

More scattered thoughts:

• There are innumerable ways in which the input topology could be funky and misbehave with a well-written first implementation of this feature. I think we should be prudent and error out in a few ways but not worry about covering absolutely every possible way in which the inputs could be too wonky to handle
• If it's called Topology.solvate and only works with water/aqueous solutions ... that probably makes total sense to a biophysicist but a chemical engineer would think that's incorrect.
• Chapin has a SLTCAP implementation in OpenFF code Implement SLTCAP in solvate_topology openff-interchange#1273 - why can't we port this over? Is it because that wraps Packmol and this code wouldn't
• Starting with Literal[...] or otherwise restricting what ions could be passed through in the first implementation is a good idea to me
• As much as I hate the consequences on type-checking, it might be useful to have the ion arguments to be str | Molecule | list[str]. If I was a user given a fatal including of you passed 'Molecule with name '' and SMILES '[Na+]'' but expected '[Na+]' ... I would not be so happy.
  • I assume here we'd only take a single cation and single anion. Lists of each are a niche use case to my knowledge
• I think we could make things significantly easier with modest reduction of feature scope by taking in box information or number of waters.
• Calling Topology.sovlate on a box-less topology could link to a nice user guide in the error message
• Prefer modifying topology in-place for the same reasons described
• What to do when OpenMM returns an error? Pass the error through, abort and don't modify the object, or something else altogether?
• Intuitively, I think that the value of Topology.box_vectors matters and we should treat it as intentionally there. The realistic scenario of a PDB file having different periodicity than the user wants is a data cleanliness issue and I'd love to see that dealt with in another way (input handling or something like that). My hope could be unrealistic

[Yoshanuikabundi]

I would love to go with option 1. I have personally been using packmol recently specifically because I needed to exactly define the number of waters - this was for initializing a replica exchange simulation with different conformations of a solute molecule at different rungs (I didn't realize Modeller provided this). So I'll gently push back that there are situations where it matters whether there are 4095 or 4096 waters. It might not be a use case we need to implement in a first pass, but I'd be limited in where I could dogfood the solvation method until we implemented it.

Quickly responding to some of your other thoughts:

• re: non-aqueous solvents - with option 1 a very natural extension to the API would be to add a tile: Topology argument that takes a pre-equilibrated box as a topology and tiles it instead of the default water box. We'd need to re-implement most/all of Modeller's method at that point though
• I think we can port Chapin's SLTCAP implementation over - but AFAIK all SLTCAP does is choose the number of ions to add, so we still need to do the adding ourselves. Modeller doesn't have a n_cation argument or similar.
• Taking in box info or number of waters is a cool idea - via option 1, solvate() could have an optional n_solvent argument, but it only works when the box is unspecified on the topology (allowing it to choose the box)
• Would have to think more about what to do when OpenMM raises an error, but initial thought would be wrap it in an OpenFF error to leave the option open of doing our own error handling in the future without a breaking change (for example, if we need to re-implement openmm stuff to handle non-aqueous solvents or custom tiles)

I broadly agree with all your other points!

[j-wags]

The key thing is that this should be really simple and easy to maintain. Every bell and whistle that we add invites people with very particular use cases to become reliant on specific behaviors, which is the opposite of what we want (see AMBER's addIons and addIons2). We want a very simple "I have a protein ligand system, give me a one-stop shop to solvate it in an unobjectionable way so I can get on with running a simulation".

• I'm in favor of a really minimal Option 1 -
  • I was leaning toward 2 because the box-setting methods are even more API points to add, but Matt made a great point that we already have a field for box vectors so we should refer to that
  • We could start with just Topology.generate_box_with_padding as @Yoshanuikabundi described and see if requests come in for other things or consider upstreaming popular implementations.
• I'm against having a place for users to specify n_waters.
  • This would be the only non-volume-based addition method and require its own branches of logic (might this lead to weird edge cases with SLTCAP?)
  • If we had more personnel we could entertain more custom use cases. But this should be able to form the basis of a lot of workflows. For example, folks who want a specific number of waters in the final system can use our functionality to add too many waters, and then remove some
• I'm weakly in favor of writing our own ion addition logic instead of using Modeller's, since that would let us use SLTCAP. Per @Yoshanuikabundi's suggestion, we should place the ions by replacing random waters (Modeller's docs says this is what they do too)
  • But I'd like to explicitly declare out of scope any sort of "smart" ion addition where it eg places Mg2+ ions to counter nucleic acid backbone charge. This kind of thing is well and good for specific circumstances, but it's opinionated and has tradeoffs, and it will make our implied behavior surface huge. The only thing we should offer is random ion placement.

...accepting molecules as ions...

One advantage of accepting (and perhaps ONLY accepting) Molecules in the ions fields is that the users could set their atom/residue info as desired.

Would have to think more about what to do when OpenMM raises an error, but initial thought would be wrap it in an OpenFF error to leave the option open of doing our own error handling in the future without a breaking change

👍 Agree

re: non-aqueous solvents - with option 1 a very natural extension to the API would be to add a tile: Topology argument that takes a pre-equilibrated box as a topology and tiles it instead of the default water box. We'd need to re-implement most/all of Modeller's method at that point though

👍 Agree, this would be a great option to consider adding in the future.

[Yoshanuikabundi]

Okay sounds like we've settled on Option 1! As an initial implementation, something like:

RHOMBIC_DODECAHEDRON: Array = ... # Like in openff-packmol, this is the box vectors with unit width
...

class Topology:
    def solvate(
        self,
        salt_conc: Quantity = Quantity(0.0, "mole/liter"),
        positive_ion: Molecule = Molecule.from_smiles("[Na+]"),
        negative_ion: Molecule = Molecule.from_smiles("[Cl-]"),
    ) -> Topology: ...

    def with_box_from_padding(
        self,
        distance: Quantity, # Should this have a default?
        shape: Array = RHOMBIC_DODECAHEDRON,
    ) -> Topology: ...

In that case, I'd question if we want to depend on openmmforcefields. It adds an enormous layer of possible bugs that we only have some control over, just to get a few LJ parameters that we already have into a format that Modeller understands. While playing with an implementation for instance, I found that I had to strip residue information out of my topology and then add it back in after solvation, as otherwise the force field would try to assign parameters to it residue-by-residue - which obviously failed as I didn't provide a protein force field. The simplest local implementation also requires a round trip through OpenMM, though we could work around this. If instead we depended on OpenMM, we could still use the water boxes they provide, but do the tiling, trimming, and clashing water molecule deleting ourselves. This is more code for us in the short term, but it's very simple code, and much of it could just be vendored/copied from Modeller. It would also give us an independent extant code path for testing.

We should also decide if solvate() is an in-place transformation or if it returns a modified copy of the topology. I have a weak preference for return a modified copy because it facilitates chaining - eg, top = molecule.to_topology().generate_box_with_padding(...).solvate(...).

I've dug into SLTCAP a bit more, I found a paper that has a really nice break down, an alternative form that can be computed from the number of waters instead of the volume, and a simplifying approximation that becomes exact when the number of waters is much greater than the solute charge (ie, almost always). I'd advocate we use this equivalent alternative form (not the approximation), as it's quite common in my experience for box sizes to change during constant-pressure equilibration, especially when the tile was equilibrated under a water model with a different density or at a different ionic strength (eg, 0). It seems to be at least related to how Chapin implemented it in the packmol wrapper too. Unfortunately neither this nor the original SLTCAP paper really addresses the fact that waters are generally replaced by ions rather than ions being added, but I think this is kind of fine since the actual parameter that SLTCAP is derived from is the volume of solvent, which is roughly constant under transformation from water to ion - it's at least better than computing the volume from the number of waters and then adding ions on top.

Regarding n_waters - I don't think this would complicate SLTCAP but I guess at the end of the day Jeff's right, my use case would be solved by just putting all conformations in boxes of the same size and then deleting waters until they all have the same number. So I guess I'm happy to leave it out for now.

One advantage of accepting (and perhaps ONLY accepting) Molecules in the ions fields is that the users could set their atom/residue info as desired.

This raises a point about the water - should we add atom/residue information to it? I'd suggest we at the very least set atom.name and atom.metadata["residue_name"], as these are important for visualization. And if we're doing that for water, there's no reason we couldn't do it for salt as well.

---

Did anyone like the compute_bounding_sphere_diameter idea? It would replace with_box_from_padding at definition with:

def bounding_sphere_diameter(self) -> Quantity: ... #  Returns float with units of length

And at use with:

topology.box_vectors = RHOMBIC_DODECAHEDRON * (topology.bounding_sphere_diameter() + 2 * 1.4 * units.nm)

Would be a bit less code for us, avoids setting the with_box precedent, and is more flexible for the user. It also avoids the "should padding be my cutoff distance, or twice my cutoff distance" ambiguity.

[j-wags]

I found that I had to strip residue information out of my topology and then add it back in after solvation, as otherwise the force field would try to assign parameters to it residue-by-residue

Stripping residue info will be unnecessary in the next release of OpenMMForceFields. This addition of multi-residue support to the OMMFFs stack is what made me open to the idea of this wrapper in our API, and I'm comfortable relying on the external implementation.

We should also decide if solvate() is an in-place transformation or if it returns a modified copy of the topology

I don't have a strong preference for this particular method. I'm thinking about how to be consistent with our existing API:

• Molecule.generate_conformers and Molecule.assign_partial_charges modify in place
• Molecule.enumerate_stereoisomers/tautomers/protomers return lists of molecules (because those operations unavoidably create distinct molecules)
• Topology.add_molecule(s) modifies in place

So I guess I'm in favor of in-place modification since it would be most consistent, and that's the biggest factor for me.

This raises a point about the water - should we add atom/residue information to it?

Modeller will have some default behavior for water residue assignment and I recommend we just use what it gives us. We'll still need to make up some scheme for assigning residue number+chain to our ions since those would be added as a different step. No strong preference on how we do that, though "assign ion residue+chain info similar to how Modeller would have done it" would be a good approach (with user overrides if they provide custom residue or atom names on the input ions)

Did anyone like the compute_bounding_sphere_diameter idea?

I'm slightly in favor of your compute_bounding_sphere_diameter proposal, but would be happy with either solution, as long as we pick just one to keep the API surface minimal.