Update conformer energies example (#2101)

* Update conformer energies example

* Unify some prose

* Avoid units diversion

Author: mattwthompson

examples/conformer_energies/conformer_energies.ipynb

@@ -21,24 +21,9 @@
    "metadata": {
     "tags": []
    },
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "cfce2d5082774a069a379ab1dd3237eb",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": []
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "import openmm\n",
-    "from openff.interchange import Interchange\n",
-    "from openff.units.openmm import from_openmm\n",
+    "from openff.interchange.drivers.openmm import get_openmm_energies\n",
     "from rdkit.Chem import rdMolAlign\n",
     "\n",
     "from openff.toolkit import ForceField, Molecule\n",
@@ -67,7 +52,27 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "089a1b3c120542d79f178ee6179c8cf4",
+       "model_id": "3dd98e55d74b4334adb0f4e207bd0874",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/mattthompson/.local/share/mamba/envs/openff-toolkit-test-rdkit/lib/python3.12/site-packages/nglview/__init__.py:12: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.\n",
+      "  import pkg_resources\n"
+     ]
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "f792d6accb2846c1b4dcee7c9d80a849",
        "version_major": 2,
        "version_minor": 0
       },
@@ -129,18 +134,17 @@
    },
    "outputs": [],
    "source": [
-    "# Load the openff-2.1.0 force field appropriate for vacuum calculations (without constraints)\n",
-    "forcefield = ForceField(\"openff_unconstrained-2.2.0.offxml\")"
+    "# Load the openff-2.2.1 force field appropriate for vacuum calculations (without constraints)\n",
+    "forcefield = ForceField(\"openff_unconstrained-2.2.1.offxml\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We'll now apply the Sage force field to the molecule to produce an [`Interchange`], and then an OpenMM [`Simulation`] object. Interchange can produce inputs for a whole host of other MD engines though!\n",
+    "We'll now apply the Sage force field to the molecule to produce an [`Interchange`], which stores the result of parametrizing with this force field\n",
     "\n",
-    "[`Interchange`]: https://docs.openforcefield.org/projects/interchange/en/stable/_autosummary/openff.interchange.Interchange.html#openff.interchange.Interchange\n",
-    "[`Simulation`]: http://docs.openmm.org/latest/api-python/generated/openmm.app.simulation.Simulation.html#openmm.app.simulation.Simulation"
+    "[`Interchange`]: https://docs.openforcefield.org/projects/interchange/en/stable/_autosummary/openff.interchange.Interchange.html#openff.interchange.Interchange"
    ]
   },
   {
@@ -161,11 +165,8 @@
    ],
    "source": [
     "print(f\"Parametrizing {molecule.name} (may take a moment to calculate charges)...\")\n",
-    "interchange = Interchange.from_smirnoff(forcefield, [molecule])\n",
-    "print(\"Done.\")\n",
-    "\n",
-    "integrator = openmm.VerletIntegrator(1 * openmm.unit.femtoseconds)\n",
-    "simulation = interchange.to_openmm_simulation(integrator)"
+    "interchange = forcefield.create_interchange(molecule.to_topology())\n",
+    "print(\"Done.\")"
    ]
   },
   {
@@ -191,21 +192,20 @@
     "minimized_molecule = Molecule(molecule)\n",
     "minimized_molecule.conformers.clear()\n",
     "\n",
+    "\n",
     "for conformer in molecule.conformers:\n",
-    "    # Tell the OpenMM Simulation the positions of this conformer\n",
-    "    simulation.context.setPositions(conformer.to_openmm())\n",
+    "    # Use this conformer to update the positions of the Interchange object\n",
+    "    interchange.positions = conformer\n",
     "\n",
-    "    # Keep a record of the initial energy\n",
-    "    initial_energies.append(simulation.context.getState(getEnergy=True).getPotentialEnergy())\n",
+    "    # Get the (total) initial energy from this conformer and store it\n",
+    "    initial_energies.append(get_openmm_energies(interchange).total_energy)\n",
     "\n",
-    "    # Perform the minimization\n",
-    "    simulation.minimizeEnergy()\n",
+    "    # Minimize using Interchange.minimize, which wraps OpenMM\n",
+    "    interchange.minimize(engine=\"openmm\")\n",
     "\n",
     "    # Record minimized energy and positions\n",
-    "    min_state = simulation.context.getState(getEnergy=True, getPositions=True)\n",
-    "\n",
-    "    minimized_energies.append(min_state.getPotentialEnergy())\n",
-    "    minimized_molecule.add_conformer(from_openmm(min_state.getPositions()))"
+    "    minimized_energies.append(get_openmm_energies(interchange).total_energy)\n",
+    "    minimized_molecule.add_conformer(interchange.positions.to(\"angstrom\"))"
    ]
   },
   {
@@ -225,7 +225,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "a4750d0da648475bb17cf557ef89b691",
+       "model_id": "c7b7dec6bc7b4c498bab35c31cedb4c3",
        "version_major": 2,
        "version_minor": 0
       },
@@ -261,16 +261,16 @@
      "text": [
       "ruxolitinib: 10 conformers\n",
       "Conformer         Initial PE        Minimized PE        RMSD\n",
-      "    1 /    10 :   -88.682 kcal/mol -124.600 kcal/mol    0.456 Å\n",
-      "    2 /    10 :   -83.492 kcal/mol -121.554 kcal/mol    0.466 Å\n",
-      "    3 /    10 :   -43.748 kcal/mol -121.709 kcal/mol    0.951 Å\n",
-      "    4 /    10 :   -81.885 kcal/mol -123.365 kcal/mol    0.695 Å\n",
-      "    5 /    10 :   -83.609 kcal/mol -120.864 kcal/mol    0.674 Å\n",
-      "    6 /    10 :   -85.390 kcal/mol -123.760 kcal/mol    0.604 Å\n",
-      "    7 /    10 :   -82.167 kcal/mol -121.679 kcal/mol    0.973 Å\n",
-      "    8 /    10 :   -91.164 kcal/mol -123.829 kcal/mol    0.413 Å\n",
-      "    9 /    10 :   -89.107 kcal/mol -124.535 kcal/mol    0.808 Å\n",
-      "   10 /    10 :   -86.407 kcal/mol -124.416 kcal/mol    0.610 Å\n"
+      "    1 /    10 :   -88.258 kcal/mol -123.916 kcal/mol    0.452 Å\n",
+      "    2 /    10 :   -82.884 kcal/mol -120.706 kcal/mol    0.433 Å\n",
+      "    3 /    10 :   -43.067 kcal/mol -120.913 kcal/mol    0.949 Å\n",
+      "    4 /    10 :   -81.524 kcal/mol -122.732 kcal/mol    0.693 Å\n",
+      "    5 /    10 :   -82.897 kcal/mol -120.303 kcal/mol    0.789 Å\n",
+      "    6 /    10 :   -84.764 kcal/mol -122.897 kcal/mol    0.588 Å\n",
+      "    7 /    10 :   -81.634 kcal/mol -120.972 kcal/mol    1.004 Å\n",
+      "    8 /    10 :   -90.702 kcal/mol -123.087 kcal/mol    0.374 Å\n",
+      "    9 /    10 :   -88.588 kcal/mol -123.813 kcal/mol    0.757 Å\n",
+      "   10 /    10 :   -85.671 kcal/mol -123.590 kcal/mol    0.634 Å\n"
      ]
     }
    ],
@@ -322,8 +322,8 @@
     "    output.append(\n",
     "        [\n",
     "            i + 1,\n",
-    "            init_energy.value_in_unit(openmm.unit.kilocalories_per_mole),\n",
-    "            min_energy.value_in_unit(openmm.unit.kilocalories_per_mole),\n",
+    "            init_energy.m_as(\"kilocalories_per_mole\"),\n",
+    "            min_energy.m_as(\"kilocalories_per_mole\"),\n",
     "            minimization_rms,\n",
     "        ]\n",
     "    )\n",
@@ -361,7 +361,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.12.11"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {

examples/conformer_energies/conformer_energies.py

@@ -1,30 +1,26 @@
 import argparse
 
-import openmm
-from openff.interchange import Interchange
-from openff.units.openmm import from_openmm
+from openff.interchange.drivers.openmm import get_openmm_energies
 from rdkit.Chem import rdMolAlign
 
-from openff.toolkit import ForceField, Molecule, RDKitToolkitWrapper
+from openff.toolkit import ForceField, Molecule
+from openff.toolkit.utils import get_data_file_path
 
 
 def compute_conformer_energies_from_file(filename):
-    # Load in the molecule and its conformers.
-    # Note that all conformers of the same molecule are loaded as separate Molecule objects
-    # If using a OFF Toolkit version before 0.7.0, loading SDFs through RDKit and OpenEye may provide
-    # different behavior in some cases. So, here we force loading through RDKit to ensure the correct behavior
-    rdktkw = RDKitToolkitWrapper()
-    loaded_molecules = Molecule.from_file(filename, toolkit_registry=rdktkw)
-    # The logic below only works for lists of molecules, so if a
-    # single molecule was loaded, cast it to list
+    # First, load conformers from an SDF file.
+    loaded_molecules = Molecule.from_file(
+        get_data_file_path("molecules/ruxolitinib_conformers.sdf"),
+    )
+
+    # Normalize to list
     try:
         loaded_molecules = [*loaded_molecules]
     except TypeError:
         loaded_molecules = [loaded_molecules]
 
-    # Collatate all conformers of the same molecule
-    # NOTE: This isn't necessary if you have already loaded or created multi-conformer molecules;
-    # it is just needed because our SDF reader does not automatically collapse conformers.
+    # from_file loads each entry in the SDF into its own molecule,
+    # so collapse conformers into the same molecule
     molecule = loaded_molecules.pop(0)
     for next_molecule in loaded_molecules:
         if next_molecule == molecule:
@@ -45,13 +41,12 @@ def compute_conformer_energies_from_file(filename):
         + f" ({molecule.name})"
     )
 
-    # Load the openff-2.1.0 force field appropriate for vacuum calculations (without constraints)
-    forcefield = ForceField("openff_unconstrained-2.1.0.offxml")
+    # Load the openff-2.2.1 force field appropriate for vacuum calculations (without constraints)
+    forcefield = ForceField("openff_unconstrained-2.2.1.offxml")
+    # Create an Interchange object, which stores the result of parametrizing with this force field
     print(f"Parametrizing {molecule.name} (may take a moment to calculate charges)...")
-    interchange = Interchange.from_smirnoff(forcefield, [molecule])
+    interchange = forcefield.create_interchange(molecule.to_topology())
     print("Done.")
-    integrator = openmm.VerletIntegrator(1 * openmm.unit.femtoseconds)
-    simulation = interchange.to_openmm_simulation(integrator)
 
     # We'll store energies in two lists
     initial_energies = []
@@ -62,20 +57,18 @@ def compute_conformer_energies_from_file(filename):
     minimized_molecule.conformers.clear()
 
     for conformer in molecule.conformers:
-        # Tell the OpenMM Simulation the positions of this conformer
-        simulation.context.setPositions(conformer.to_openmm())
+        # Use this conformer to update the positions of the Interchange object
+        interchange.positions = conformer
 
-        # Keep a record of the initial energy
-        initial_energies.append(simulation.context.getState(getEnergy=True).getPotentialEnergy())
+        # Get the (total) initial energy from this conformer and store it
+        initial_energies.append(get_openmm_energies(interchange).total_energy)
 
-        # Perform the minimization
-        simulation.minimizeEnergy()
+        # Minimize using Interchange.minimize, which wraps OpenMM
+        interchange.minimize(engine="openmm")
 
         # Record minimized energy and positions
-        min_state = simulation.context.getState(getEnergy=True, getPositions=True)
-
-        minimized_energies.append(min_state.getPotentialEnergy())
-        minimized_molecule.add_conformer(from_openmm(min_state.getPositions()))
+        minimized_energies.append(get_openmm_energies(interchange).total_energy)
+        minimized_molecule.add_conformer(interchange.positions.to("angstrom"))
 
     n_confs = molecule.n_conformers
     print(f"{molecule.name}: {n_confs} conformers")
@@ -122,8 +115,8 @@ def compute_conformer_energies_from_file(filename):
         output.append(
             [
                 i + 1,
-                init_energy.value_in_unit(openmm.unit.kilocalories_per_mole),
-                min_energy.value_in_unit(openmm.unit.kilocalories_per_mole),
+                init_energy.m_as("kilocalories_per_mole"),
+                min_energy.m_as("kilocalories_per_mole"),
                 minimization_rms,
             ]
         )