Update Packmol defaults (#1332)

* API: Update default packing behavior

* FIX: Fix box scale-up test

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* FIX: Fix tests

* Lint

* DOC: Update packed box example

* API: Avoid changing private function default

* ENH: Use PBCs in Packmol wrappers

* DOC: Document behavior changes

* FIX: Try to be backwards-compatible with older Packmol

* FIX: Remove temporary code

* DOC: Update release history

* MAINT: Try a build without AmberTools?

* Revert "MAINT: Try a build without AmberTools?"

This reverts commit f16de82.

* MAINT: Debug

* Fix?

* FIX: Fix parsing output if no error

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* LINT

* REF: Only use Packmol PBCs with rectangular boxes

* Apply suggestion from @Yoshanuikabundi

Co-authored-by: Josh A. Mitchell <[email protected]>

* Apply suggestion from @Yoshanuikabundi

Co-authored-by: Josh A. Mitchell <[email protected]>

* Apply suggestion from @Yoshanuikabundi

Co-authored-by: Josh A. Mitchell <[email protected]>

* FIX: Fix double-uninstall

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* DOC: Update docstrings

* FIX: Simpler version parsing

* Apply suggestions from code review

Co-authored-by: Josh A. Mitchell <[email protected]>

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Josh A. Mitchell <[email protected]>

Author: 3 people

docs/releasehistory.md

@@ -33,6 +33,16 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 
 ## 0.4.8 - 2025-10-09
 
+### Behavior changes
+
+* #1332 Introduces the following behavior changes to the private Packmol wrappers:
+  * Packmol version 20.15.0 or newer is recommended.
+  * Periodic boundary conditions are accounted for when placing molecules if the box is orthorhombic and Packmol version 20.15.0 or newer is installed; this is the minimum version supporting this feature.
+    * This is functionally the same as the previous behavior, so no workaround is needed to recover it.
+  * The `tolerance` parameter is subtracted from computed box lengths when placing molecules if a Packmol version older than 20.15.0 is installed; this is the previous behavior.
+  * The `target_density` is used in box size calculations with modification; previously, box volumes were scaled up by a factor of 1.1.
+    * The previous behavior can be restored by passing scaling the `target_density` argument down by a factor of 1.1.
+
 ### Bug fixes
 
 * #1340 Fixes JSON deserialization issues when charges come from `NAGLCharges` or `ChargeIncrementModel`.

examples/packed_box/packed_box.ipynb

@@ -107,14 +107,6 @@
    "cell_type": "markdown",
    "id": "9",
    "metadata": {},
-   "source": [
-    "We can get the positions as an array from the PDB file object:"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "10",
-   "metadata": {},
    "source": [
     "## Parametrize with Interchange\n",
     "\n",
@@ -124,7 +116,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "11",
+   "id": "10",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -135,7 +127,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "12",
+   "id": "11",
    "metadata": {},
    "source": [
     "We can visualize it (though, since we can't see the stored physics parameters, it'll look the same):"
@@ -144,7 +136,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "13",
+   "id": "12",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -153,16 +145,34 @@
   },
   {
    "cell_type": "markdown",
+   "id": "13",
+   "metadata": {},
+   "source": [
+    "Packmol places molecules at random positions in the box, satisfying geometric constraints but no thermodynamic considerations. The resulting configurations need to be subject to energy minimization and equilibration routines before being used in production simulations. A full equilibration routine is outside the scope of this example, but a brief energy minimization immediately after packing lessens the likelihood of unfavorable energetics or crashes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
    "id": "14",
    "metadata": {},
+   "outputs": [],
+   "source": [
+    "interchange.minimize(max_iterations=1000)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15",
+   "metadata": {},
    "source": [
     "And we can calculate and compare energies with different MD engines! (The LAMMPS exporter isn't optimized yet for large systems, so we're only looking at OpenMM, GROMACS, and Amber.)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "15",
+   "id": "16",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -187,7 +197,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.13.7"
   }
  },
  "nbformat": 4,

openff/interchange/_tests/unit_tests/components/test_packmol.py

@@ -51,15 +51,24 @@ class TestPackmolWrapper:
             Quantity(234, "angstrom**3"),
         ],
     )
-    def test_scale_box(self, box, volume):
+    @pytest.mark.parametrize(
+        "factor",
+        [
+            0.9,
+            1.0,
+            1.1,
+            2.0,
+        ],
+    )
+    def test_scale_box(self, box, volume, factor):
         """Test that _scale_box() produces a box with the desired volume."""
-        scaled_box = _scale_box(box, volume)
+        scaled_box = _scale_box(box, volume, factor)
         a, b, c = scaled_box
 
         # | (a x b) . c | is the volume of the box
         # _scale_box uses numpy.linalg.det instead
         # linear dimensions are scaled by 1.1, so volumes are scaled by 1.1 ** 3
-        assert numpy.isclose(numpy.abs(numpy.dot(numpy.cross(a, b), c)), volume * 1.1**3)
+        assert numpy.isclose(numpy.abs(numpy.dot(numpy.cross(a, b), c)), volume * factor**3)
 
         assert str(scaled_box.u) == "angstrom"
 
@@ -487,6 +496,7 @@ def test_solvate_topology_highly_charged_low_nacl_conc(self, n_monomers):
 
         topology = solvate_topology(
             solute.to_topology(),
+            target_density=Quantity(0.8, "gram / milliliter"),
             nacl_conc=Quantity(1e-3, "molar"),
             padding=Quantity(2, "nanometer"),
         )

openff/interchange/components/_packmol.py

@@ -14,6 +14,7 @@
 from numpy.typing import ArrayLike, NDArray
 from openff.toolkit import Molecule, Quantity, RDKitToolkitWrapper, Topology
 from openff.utilities.utilities import requires_package, temporary_cd
+from packaging.version import Version
 
 from openff.interchange.exceptions import PACKMOLRuntimeError, PACKMOLValueError
 
@@ -345,10 +346,14 @@ def _box_from_density(
     total_mass = sum(sum([atom.mass for atom in molecule.atoms]) * n for molecule, n in zip(molecules, n_copies))
     volume = total_mass / target_density
 
-    return _scale_box(box_shape, volume)
+    return _scale_box(
+        box=box_shape,
+        volume=volume,
+        box_scaleup_factor=1.0,
+    )
 
 
-def _scale_box(box: NDArray, volume: Quantity, box_scaleup_factor=1.1) -> Quantity:
+def _scale_box(box: NDArray, volume: Quantity, box_scaleup_factor: float = 1.1) -> Quantity:
     """
     Scale the parallelepiped spanned by ``box`` to the given volume.
 
@@ -442,12 +447,33 @@ def _create_molecule_pdbs(molecules: list[Molecule]) -> list[str]:
     return pdb_file_names
 
 
+def _get_packmol_version() -> Version:
+    """Return the version of packmol installed."""
+    with temporary_cd(tempfile.mkdtemp()):
+        try:
+            result = subprocess.run(
+                _find_packmol(),
+                input="",
+                stdout=subprocess.PIPE,
+                stderr=subprocess.STDOUT,
+                text=True,
+            )
+
+            found_version = Version(result.stdout.split("Version ")[1].split(" ")[0])
+
+        except Exception as error:
+            raise PACKMOLRuntimeError(f"Unexpected error ({type(error)}) while running Packmol.") from error
+
+    return found_version
+
+
 def _build_input_file(
     molecule_file_names: list[str],
     molecule_counts: list[int],
     structure_to_solvate: str | None,
     box_size: Quantity,
     tolerance: Quantity,
+    rectangular: bool = False,
 ) -> tuple[str, str]:
     """
     Construct the packmol input file.
@@ -462,10 +488,11 @@ def _build_input_file(
         The path to the structure to solvate.
     box_size
         The lengths of each side of the box we want to fill. This is the box
-        size of the rectangular brick representation of the simulation box; the
-        packmol box will be shrunk by the tolerance.
+        size of the rectangular brick representation of the simulation box
     tolerance
         The packmol convergence tolerance.
+    rectangular
+        Whether the box is rectangular (True) or triclinic (False).
 
     Returns
     -------
@@ -475,7 +502,12 @@ def _build_input_file(
         The path to the output file.
 
     """
-    box_size = (box_size - tolerance).m_as("angstrom")
+    if rectangular:
+        PACKMOL_USE_PBC = _get_packmol_version() >= Version("20.15.0")
+    else:
+        PACKMOL_USE_PBC = False
+
+    box_size = box_size.m_as("angstrom") if PACKMOL_USE_PBC else (box_size - tolerance).m_as("angstrom")
     tolerance = tolerance.m_as("angstrom")
 
     # Add the global header options.
@@ -487,6 +519,9 @@ def _build_input_file(
         "",
     ]
 
+    if PACKMOL_USE_PBC:
+        input_lines.append(f"pbc 0. 0. 0. {box_size[0]} {box_size[1]} {box_size[2]}")
+
     # Add the section of the molecule to solvate if provided.
     if structure_to_solvate is not None:
         input_lines.extend(
@@ -508,7 +543,7 @@ def _build_input_file(
             [
                 f"structure {file_name}",
                 f"  number {count}",
-                f"  inside box 0. 0. 0. {box_size[0]} {box_size[1]} {box_size[2]}",
+                "" if PACKMOL_USE_PBC else f"  inside box 0. 0. 0. {box_size[0]} {box_size[1]} {box_size[2]}",
                 "end structure",
                 "",
             ],
@@ -633,14 +668,20 @@ def pack_box(
 
     Notes
     -----
-    Returned topologies may have smaller or larger box vectors than what would be defined by the
-    target density if the box vectors are determined by `target_density`. When calling Packmol, each
-    linear dimension of the box is scaled up by 10%.  However, Packmol by default adds a small
-    buffer (defined by the `tolerance` argument which otherwise defines the minimum distance,
-    default 2 Angstrom) at the end of the packed box, which causes small voids when tiling copies of
-    each periodic image. This void is removed in hopes of faster equilibration times but the box
-    density is slightly increased as a result. These changes may cancel each other out or result in
-    larger or smaller densities than the target density, depending on argument values.
+    Packmol places molecules at random positions in the box, satisfying
+    geometric constraints but no thermodynamic considerations. The resulting
+    configurations need to be subjected to energy minimization and equilibration
+    routines before being used in production simulations.
+
+    Orthorhombic periodic boundary conditions are accounted for during packing
+    if using Packmol version 20.15.0 or newer.
+
+    Non-orthorhombic triclinic boxes, and all boxes in Packmol versions older
+    than 20.15.0, are accounted for by packing the
+    `"brick" representation <https://manual.gromacs.org/current/reference-manual/algorithms/periodic-boundary-conditions.html#fig-pbc>`__
+    after shrinking each dimension by the ``tolerance``. This introduces a small
+    void at the edges of the box to avoid clashes. This void will quickly be
+    filled in during equilibration.
 
     """
     # Make sure packmol can be found.
@@ -663,6 +704,8 @@ def pack_box(
         target_density,
     )
 
+    is_rectangular = numpy.all(box_shape == numpy.diag(numpy.diagonal(box_shape)))
+
     # Estimate the box_vectors from mass density if one is not provided.
     if target_density is not None:
         box_vectors = _box_from_density(
@@ -683,8 +726,7 @@ def pack_box(
             + "05_other_features.html#periodic-boundary-conditions",
         )
 
-    # Compute the dimensions of the equivalent brick - this is what packmol will
-    # fill
+    # Compute the dimensions of the equivalent brick - this is what packmol will fill
     brick_size = _compute_brick_from_box_vectors(box_vectors)
 
     # Center the solute
@@ -721,6 +763,7 @@ def pack_box(
             solute_pdb_filename,
             brick_size,
             tolerance,
+            rectangular=is_rectangular,
         )
 
         with open(input_file_path) as file_handle:
@@ -831,7 +874,7 @@ def solvate_topology(
     """
     Add water and ions to neutralise and solvate a topology.
 
-    The [SLTCAP](10.1021/acs.jctc.7b01254)  method is used to determine the number of each ion to add.
+    The [SLTCAP](10.1021/acs.jctc.7b01254) method is used to determine the number of each ion to add.
 
     Parameters
     ----------
@@ -876,8 +919,20 @@ def solvate_topology(
 
     Notes
     -----
-    Returned topologies may have larger box vectors than what would be defined
-    by the target density.
+    Packmol places solvent molecules at random positions in the box, satisfying
+    geometric constraints but no thermodynamic considerations. The resulting
+    configurations need to be subjected to energy minimization and equilibration
+    routines before being used in production simulations.
+
+    Orthorhombic periodic boundary conditions are accounted for during packing
+    if using Packmol version 20.15.0 or newer.
+
+    Non-orthorhombic triclinic boxes, and all boxes in Packmol versions older
+    than 20.15.0, are accounted for by packing the
+    `"brick" representation <https://manual.gromacs.org/current/reference-manual/algorithms/periodic-boundary-conditions.html#fig-pbc>`__
+    after shrinking each dimension by the ``tolerance``. This introduces a small
+    void at the edges of the box to avoid clashes. This void will quickly be
+    filled in during equilibration.
 
     Returned topologies may have ion concentrations higher than the value of the
     the `nacl_conc` argument.
@@ -1040,8 +1095,20 @@ def solvate_topology_nonwater(
 
     Notes
     -----
-    Returned topologies may have larger box vectors than what would be defined
-    by the target density.
+    Packmol places solvent molecules at random positions in the box, satisfying
+    geometric constraints but no thermodynamic considerations. The resulting
+    configurations need to be subjected to energy minimization and equilibration
+    routines before being used in production simulations.
+
+    Orthorhombic periodic boundary conditions are accounted for during packing
+    if using Packmol version 20.15.0 or newer.
+
+    Non-orthorhombic triclinic boxes, and all boxes in Packmol versions older
+    than 20.15.0, are accounted for by packing the
+    `"brick" representation <https://manual.gromacs.org/current/reference-manual/algorithms/periodic-boundary-conditions.html#fig-pbc>`__
+    after shrinking each dimension by the ``tolerance``. This introduces a small
+    void at the edges of the box to avoid clashes. This void will quickly be
+    filled in during equilibration.
     """
     _check_box_shape_shape(box_shape)