Merge remote-tracking branch 'upstream/develop' into preset-charges

.github/workflows/ci.yaml

@@ -98,10 +98,10 @@ jobs:
  python -m pytest $COV openff/interchange/ -r fExs -n logical --durations=10
 
  - name: Run small molecule regression tests
- if: ${{ matrix.python-version == '3.10' && matrix.openeye == true && matrix.openmm == true }}
+ if: ${{ matrix.openeye == true && matrix.openmm == true }}
  run: |
- micromamba install deepdiff rich click -c conda-forge
- python -m pip install git+https:/openforcefield/interchange-regression-testing.git
+ micromamba install "deepdiff =5" rich click -c conda-forge
+ python -m pip install git+https:/openforcefield/interchange-regression-testing.git@fa07b62824c39af591142b402bd68103f4043f8b
 
  create_openmm_systems \
  --input "regression_tests/small-molecule/input-topologies.json" \

docs/releasehistory.md

@@ -22,6 +22,10 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 * Several classes and methods which were deprecated in the 0.3 line of releases are now removed.
 * Previously-deprecated examples are removed.
 * `ProperTorsionKey` no longer accepts an empty tuple as atom indices.
+* Default densities for Packmol wrapper functions have been lowered to 0.9 g/cc and 0.6 g/cc for water and non-water solvent, respectively.
+* Argument `mass_density` to some packing functions has been renamed to `target_density` for consistency and to better reflect its usage.
+* Topologies returned by packing functions have boxes scaled up by 10% in linear dimensions compared to the size implied by the target density.
+* An error is now raised when HMR would result in an OpenMM particle (aside from virtual sites) having negative (or zero) mass.
 * Fixes a regression in which some `ElectrostaticsCollection.charges` properties did not return cached values.
 * The `charge_from_molecules` argument must include only molecules that contain partial charges and are non-isomorphic with each other.
 * The `charge_from_molecules` argument as used by the OpenFF Toolkit is handled internally as `molecules_with_preset_charges`.

docs/using/edges.md

@@ -30,6 +30,17 @@ The charges specified by the force field can be overridden by providing molecule
 
 Using preset charges with virtual sites is discouraged as it can provide surprising results.
 
+## Quirks of core OpenFF objects
+
+Future refactors may remove the side effects of these quirks, but currently there are some
+surprising inconsistencies in some API points between different OpenFF tools.
+
+### Contents of `Interchange.topology` and `Interchange` may not always be in sync
+
+Currently, the `Interchange.topology` attribute is defined by the OpenFF Toolkit's `Topology` object, which is feature-rich in cheminformatics functionality but not designed for MM interoperability. Most importantly, that representation does not include virtual sites (because molecules do not have virtual sites/dummy atoms). As a result, functionality involving virtual sites must go through `Interchange` API points instead of `Interchange.topology`.
+
+For example, `Interchange.topology.get_positions()` will never include positions of virtual sites. To get the positions of a system with virtual sites included, use `Interchange.get_positions()`. The default behavior of `Interchange.positions` is also to return positions without virtual sites, but this may change in the future.
+
 ## Quirks of `from_openmm`
 
 ### Modified masses are ignored

examples/lammps/lammps.ipynb

@@ -34,7 +34,6 @@
  "from openmm.app import PDBFile\n",
  "\n",
  "from openff.interchange import Interchange\n",
- "from openff.interchange.components.mdconfig import MDConfig\n",
  "\n",
  "# Read a structure from the Toolkit's test suite into a Topology\n",
  "pdbfile = PDBFile(get_data_file_path(\"systems/packmol_boxes/propane_methane_butanol_0.2_0.3_0.5.pdb\"))\n",
@@ -87,14 +86,14 @@
  "metadata": {},
  "outputs": [],
  "source": [
- "interchange.to_lammps(\"interchange.lmp\")"
+ "interchange.to_lammps_datafile(\"interchange.lmp\")"
  ]
  },
  {
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Now we need to write an input file for LAMMPS. Parts of these input files depend on force field parameters, so we should use a sample input file written for our interchange as a starting point. We can generate such a sample file from `MDConfig`:"
+ "Now we need to write an input file for LAMMPS. Parts of these input files depend on force field parameters, so we should use a sample input file written for our interchange as a starting point. We can generate such a sample file with the `to_lammps_input()` method:"
  ]
  },
  {
@@ -103,17 +102,32 @@
  "metadata": {},
  "outputs": [],
  "source": [
- "mdconfig = MDConfig.from_interchange(interchange)\n",
- "mdconfig.write_lammps_input(input_file=\"auto_generated.in\", interchange=interchange)\n",
- "with open(\"auto_generated.in\") as f:\n",
+ "interchange.to_lammps_input(\"interchange_pointenergy.in\", data_file=\"interchange.lmp\")\n",
+ "with open(\"interchange_pointenergy.in\") as f:\n",
  " print(f.read())"
  ]
  },
  {
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "That sample file will only perform a single point energy calculation; here's a more complete file that includes the above parameters but will run an actual MD simulation. Note that `out.lmp` is changed to `interchange.lmp` to reflect the filename we used earlier:\n",
+ "Note that the `read_data` line must match the data file name - in this case, `interchange.lmp`. That's why `to_lammps_input` needs the data file name. We could alternatively use the `to_lammps` method to write both files in a consistent way:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "interchange.to_lammps(\"interchange\") # writes interchange.lmp and interchange_pointenergy.in"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "That sample file will only perform a single point energy calculation; here's a more complete file that includes the above parameters but will run an actual MD simulation. \n",
  "\n",
  "<div class=\"alert alert-warning\" style=\"max-width: 700px; margin-left: auto; margin-right: auto;\">\n",
  " <b>⚠️ Don't use example input files in production</b><br />\n",
@@ -208,6 +222,13 @@
  "traj.make_molecules_whole()\n",
  "nglview.show_mdtraj(traj)"
  ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
  }
  ],
  "metadata": {
@@ -227,7 +248,7 @@
  "name": "python",
  "nbconvert_exporter": "python",
  "pygments_lexer": "ipython3",
- "version": "3.10.14"
+ "version": "3.11.0"
  }
  },
  "nbformat": 4,

examples/packed_box/packed_box.ipynb

@@ -79,7 +79,7 @@
  "topology = pack_box(\n",
  " molecules=molecules,\n",
  " number_of_copies=4 * [200],\n",
- " mass_density=850 * unit.kilogram / unit.meter**3,\n",
+ " target_density=850 * unit.kilogram / unit.meter**3,\n",
  " box_shape=UNIT_CUBE,\n",
  ")"
  ]

openff/interchange/_annotations.py

@@ -3,6 +3,7 @@
 from typing import Annotated, Any
 
 import numpy
+from annotated_types import Gt
 from openff.toolkit import Quantity
 from pydantic import (
  AfterValidator,
@@ -13,6 +14,8 @@
  WrapValidator,
 )
 
+PositiveFloat = Annotated[float, Gt(0)]
+
 
 def _has_compatible_dimensionality(
  quantity: Quantity,

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

@@ -2,6 +2,8 @@
 Units tests for openff.interchange.components._packmol
 """
 
+import pathlib
+
 import numpy
 import pytest
 from openff.toolkit.topology import Molecule
@@ -53,9 +55,12 @@ def test_scale_box(self, box, volume):
  """Test that _scale_box() produces a box with the desired volume."""
  scaled_box = _scale_box(box, volume)
  a, b, c = scaled_box
+
  # | (a x b) . c | is the volume of the box
  # _scale_box uses numpy.linalg.det instead
- assert numpy.isclose(numpy.abs(numpy.dot(numpy.cross(a, b), c)), volume)
+ # 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 scaled_box.u == unit.angstrom
 
  @pytest.mark.parametrize(
@@ -149,6 +154,7 @@ def test_packmol_bad_box_shape(self, molecules):
  molecules[0].to_topology(),
  solvent=Molecule.from_smiles("CCCCCCO"),
  box_shape=20 * numpy.identity(4) * unit.angstrom,
+ target_density=1.0 * unit.grams / unit.milliliter,
  )
 
  def test_packmol_underspecified(self, molecules):
@@ -167,7 +173,7 @@ def test_packmol_overspecified(self, molecules):
  pack_box(
  molecules,
  number_of_copies=[1],
- mass_density=1.0 * unit.grams / unit.milliliter,
+ target_density=1.0 * unit.grams / unit.milliliter,
  box_vectors=20 * numpy.identity(3) * unit.angstrom,
  )
 
@@ -192,7 +198,7 @@ def test_packmol_water(self, molecules):
  topology = pack_box(
  molecules,
  [10],
- mass_density=1.0 * unit.grams / unit.milliliter,
+ target_density=1.0 * unit.grams / unit.milliliter,
  )
 
  assert topology is not None
@@ -397,7 +403,7 @@ def test_noninteger_serial_error(self):
  number_of_copies=[11112],
  box_shape=UNIT_CUBE,
  tolerance=1.0 * unit.angstrom,
- mass_density=0.1 * unit.grams / unit.milliliters,
+ target_density=0.1 * unit.grams / unit.milliliters,
  )
 
  @pytest.mark.slow
@@ -409,7 +415,7 @@ def test_noninteger_serial_fallback(self):
  number_of_copies=[11112],
  box_shape=UNIT_CUBE,
  tolerance=1.0 * unit.angstrom,
- mass_density=0.1 * unit.grams / unit.milliliters,
+ target_density=0.1 * unit.grams / unit.milliliters,
  )
 
  @pytest.mark.slow
@@ -420,5 +426,24 @@ def test_load_100_000_atoms(self):
  number_of_copies=[11112],
  box_shape=UNIT_CUBE,
  tolerance=1.0 * unit.angstrom,
- mass_density=0.1 * unit.grams / unit.milliliters,
+ target_density=0.1 * unit.grams / unit.milliliters,
  )
+
+ @pytest.mark.parametrize("use_local_path", [False, True])
+ def test_save_error_on_convergence_failure(self, use_local_path):
+ with pytest.raises(
+ PACKMOLRuntimeError,
+ match="failed with error code 173",
+ ):
+ pack_box(
+ molecules=[Molecule.from_smiles("CCO")],
+ number_of_copies=[100],
+ box_shape=UNIT_CUBE,
+ target_density=1000 * unit.grams / unit.milliliters,
+ working_directory="." if use_local_path else None,
+ )
+
+ if use_local_path:
+ assert "STOP 173" in open("packmol_error.log").read()
+ else:
+ assert not pathlib.Path("packmol_error.log").is_file()

openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py

@@ -8,7 +8,6 @@
 
 from openff.interchange import Interchange
 from openff.interchange._tests import MoleculeWithConformer, needs_lmp
-from openff.interchange.components.mdconfig import MDConfig
 from openff.interchange.drivers import get_lammps_energies, get_openmm_energies
 
 rng = numpy.random.default_rng(821)
@@ -105,24 +104,17 @@ def test_unique_lammps_mol_ids(
  topology.box_vectors = Quantity([4, 4, 4], "nanometer")
 
  with temporary_cd():
- lammps_input_path = Path.cwd() / "temp.in"
- lammps_data_path = Path.cwd() / "out.lmp"
+ lammps_prefix = Path.cwd() / "lammps_test"
 
  interchange = sage_unconstrained.create_interchange(topology)
- interchange.to_lammps(lammps_data_path)
-
- mdconfig = MDConfig.from_interchange(interchange)
- mdconfig.write_lammps_input(
- interchange=interchange,
- input_file=lammps_input_path,
- )
+ interchange.to_lammps(lammps_prefix)
 
  # Extract molecule IDs from data file
  with lammps.lammps(
  cmdargs=["-screen", "none", "-log", "none"],
  ) as lmp:
  lmp.file(
- "temp.in",
+ "lammps_test_pointenergy.in",
  )
  written_mol_ids = {
  mol_id

openff/interchange/_tests/unit_tests/interop/openmm/test_hmr.py

@@ -1,9 +1,9 @@
 import random
 
 import pytest
-from openff.toolkit import unit
+from openff.toolkit import Molecule, unit
 
-from openff.interchange.exceptions import UnsupportedExportError
+from openff.interchange.exceptions import NegativeMassError, UnsupportedExportError
 
 
 @pytest.mark.parametrize("reversed", [False, True])
@@ -70,6 +70,16 @@ def test_hmr_not_applied_to_water(sage, water):
  assert system.getParticleMass(particle_index) == element.hydrogen.mass
 
 
+def test_mass_is_positive(sage):
+ pytest.importorskip("openmm")
+
+ with pytest.raises(
+ NegativeMassError,
+ match="Particle with index 0 would have a negative mass after.*5.0",
+ ):
+ sage.create_interchange(Molecule.from_smiles("C").to_topology()).to_openmm(hydrogen_mass=5.0)
+
+
 def test_virtual_sites_unsupported(tip4p, water):
  with pytest.raises(UnsupportedExportError, match="with virtual sites"):
  tip4p.create_interchange(water.to_topology()).to_openmm(hydrogen_mass=2.0)