Update Packmol defauls

[mattwthompson]

For historical reasons, the Packmol wrappers by default trim the box and then expand it by 10%, which is a bit funky. Before making these functions part of the public API (#1170), I would like to propose reverting both details of this behavior such that the box_vectors argument is used for the box to pack molecules into and the requested density target_density is not modified.

• For relatively simple use cases, the packing result will be strictly different (different box and positions) but functionally similar (should equilibrate to the same similar thermodynamic state)
• These functions are currently in the private API and have no stability guarantees
• Currently, there is no guarantee that the returned box is at the requested density
• The existing behavior could be maintained by changing argument values, passing a different-sized box and/or a lower target density
• Users are expected to be responsible for equilibrating their systems before production runs

Josh goes into more detail here; this issue is mostly paraphrasing their suggestion

Pinging some stakeholders who would be affected - @IAlibay @atravitz @lilyminium @j-wags @hannaomi @orionarcher (or another Atomate2 contact, maybe @JaGeo?) @JMorado @JHoeflich1 @jaketanderson

I am strongly in favor of this change since prior to locking in a 1.0 API, but am open to reasons why this behavior would be worse and not worth the change.

[jaclark5]

I think that's a reasonable proposition provided that the periodic option is enabled and the version of packmol is restricted to be equal or greater than the version where periodic boundaries are introduced.

[IAlibay]

Background: OpenFE doesn't use packmol in it's first class Protocols. This may change post-SystemGenerator replacement, but I suspect that would be a long term change (post Interchange v1.0). So my opinions here are solely considering the more immediate term needs of pontibus.

functionally similar (should equilibrate to the same similar thermodynamic state)

The should here is doing most of the heavily lifting. From my perspective, the proposal is fine if things are shown to still work in practice. I.e. for a suitable selection within the range of solvents we are interested in, things still successfully pack (there's part of a test for this with SFEs somewhere in pontibus), and simulations don't crash / run to an acceptable outcome. Note: running just packing wouldn't be sufficient here, you'd have to actually run MD to make sure things don't crash out during NPT (we've seen cases in the past where things would only get caught that late).

This is probably a medium-low amount of work, but something that would squarely fall within the remit of pontibus maintenance, i.e. someone from OpenFF or OpenFE would have to do it.

provided that the periodic option is enabled

I may be wrong here but I was under the impression that triclinic boxes was not supported by the packmol PBC option?

Pontibus simulations that use water as a solvent are encouraged to use dodecahedron boxes. This is particularly important for the upcoming RBFE Protocol since we'd otherwise be leaving a ton of performance behind. Ideally we'd need to know that these cases would still work.

[lilyminium]

I think echoing Irfan here -- this change seems fine conceptually, and my main concern would be a) does this functionally reproduce previous behaviour (or not decrease from previous behaviour) and b) if this would be sufficiently robust to generalize. Skimming quickly over the history it looks like the scaling up the box for general stability was inherited or has been in use since before 2014 in openmoltools (choderalab/openmoltools@3621ccb); if packmol has updated since then for PBC and other stability to no longer be a concern then this seems reasonable.

The existing behavior could be maintained by changing argument values, passing a different-sized box and/or a lower target density

Would it be possible to expose the box_scaleup_factor argument in pack_box for backwards compatibility?

[hannaomi]

These changes seem fine to me as long as the PBC are still implemented by pack_box . We don't currently use the target_density argument

[mattwthompson]

Thanks all for the feedback -

I'm not planning to change anything re: Packmol's current PBC support¹.

I'm apprehensive about Interchange being responsible for the reliability and stability of MD equilibration routines; I never envisioned this to be in scope of the project. I'm iffy of the value of running a large suite of equilibration workflows (which will not land in this source code) and have always seen this as the practitioner's responsibility. My hope is that users who want the behavior of a 10% scaled-up box could just reduce the target density by the same value, a change which would be documented.

I will make an effort to be open to being convinced otherwise, but I fear this sort of project has the potential to balloon into a large amount of work without advancing the key objectives of this project. (Applying force fields to chemical topologies is already a bit of a headache, ensuring behavior of downstream simulation protocols is more.)

Footnotes

1. I forgot, or never knew, that it used to not support periodicity, if that was the case ↩

[Yoshanuikabundi]

Skimming quickly over the history it looks like the scaling up the box for general stability was inherited or has been in use since before 2014 in openmoltools (choderalab/openmoltools@3621ccb)

I think this is a coincidence; when I first adapted this implementation to support triclinic boxes it didn't have a scaleup factor.

I may be wrong here but I was under the impression that triclinic boxes was not supported by the packmol PBC option?

For as long as it has been in Interchange, pack_box has worked with triclinic boxes by computing the rectangular brick representation of the triclinic box, and then passing that as a non-periodic box to packmol. Packmol is completely unaware of the PBCs and just fills that box. There's no proposal to change this behavior.

Historical info

When I initially wrote the code, I was worried about clashes along the edges of the box. I imagined this could happen because packmol is not aware of our PBCs, so the box could have, say, two water oxygen atoms that form no clash in packmol but do clash once the PBCs are added. I "fixed" this by adding a small buffer region between the box filled by packmol and the brick representation of the triclinic box. No-one else seems to be worried about this issue and its never come up (energy minimization is probably able to deal with it fine), but people don't seem to like the vacuum created by the buffer region.

If the buffer region is not desired, we should just remove the code that adds it rather than try to compensate by growing the box. But if we grow the box for other reasons, there might be reason to keep it.

My API design concerns

I think if people want a lower density packing that can then be compressed during NVP equilibration, I think it's better API design to actually fulfill the requested density but have the user request a lower density. This preserves functionality for users who want to use NVT equilibration. I think it would be better to just reduce the default density than actually provide a density 10% lower than what is requested. Another possible approach might be to grow the solute by ten percent, pack at a lower density, and then shrink everything back down