From 6a0195915462885dd5a1a2f89211226d1a72bff4 Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Wed, 19 Apr 2023 11:33:33 +0800 Subject: [PATCH 01/12] Modifying parameter in the same manner as in Espaloma --- .../Modifying Parameters in System.ipynb | 322 ++++++++++++++++++ ...uerying Charges and Other Parameters.ipynb | 12 +- 2 files changed, 330 insertions(+), 4 deletions(-) create mode 100644 notebooks/cookbook/Modifying Parameters in System.ipynb diff --git a/notebooks/cookbook/Modifying Parameters in System.ipynb b/notebooks/cookbook/Modifying Parameters in System.ipynb new file mode 100644 index 0000000..1be8b13 --- /dev/null +++ b/notebooks/cookbook/Modifying Parameters in System.ipynb @@ -0,0 +1,322 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5cc09957", + "metadata": {}, + "source": [ + "## Modifying Parameters in System\n", + "\n", + "Some machine learning force field construction schemes, such as [Espaloma](https://pubs.rsc.org/en/content/articlehtml/2022/sc/d2sc02739a), would output parameters which can be picked up by the parameter-holding `System` object of OpenMM, which could be further used in molecular dynamics (MD) simulations just like a vanilla `System`.\n", + "In this tutorial, we show how you can modify the parameters of `System` object so you can write your own machine learning force field API.\n", + "\n", + "This example is inspired by the [Espaloma deployment](https://github.com/choderalab/espaloma/blob/master/espaloma/graphs/deploy.py) implementation." + ] + }, + { + "cell_type": "markdown", + "id": "5003093c", + "metadata": {}, + "source": [ + "First, let's load a PDB file and use Amber14 force field as placeholder parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "9c163ab2", + "metadata": {}, + "outputs": [], + "source": [ + "from openmm.app import *\n", + "from openmm import *\n", + "\n", + "pdb = PDBFile('ala_ala_ala.pdb')\n", + "forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')\n", + "system = forcefield.createSystem(pdb.topology)" + ] + }, + { + "cell_type": "markdown", + "id": "0c5c2c35", + "metadata": {}, + "source": [ + "Now we loop through the forces to get bonded components of the energy." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f92baca2", + "metadata": {}, + "outputs": [], + "source": [ + "forces = list(system.getForces())" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e8b4bef3", + "metadata": {}, + "outputs": [], + "source": [ + "harmonic_bond_force = [force for force in forces if \"HarmonicBondForce\" in force.__class__.__name__].pop()\n", + "harmonic_angle_force = [force for force in forces if \"HarmonicAngleForce\" in force.__class__.__name__].pop()\n", + "periodic_torsion_force = [force for force in forces if \"PeriodicTorsionForce\" in force.__class__.__name__].pop()" + ] + }, + { + "cell_type": "markdown", + "id": "0d8cde7e", + "metadata": {}, + "source": [ + "We can inspect each component to get a list of indices and parameters." + ] + }, + { + "cell_type": "markdown", + "id": "425b070e", + "metadata": {}, + "source": [ + "- Bonds:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "b09324e6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Bond between 10 and 4 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", + "Bond between 10 and 11 has equilibrium length 0.12290000000000001 nm and force constant 476975.9999999999 kJ/(nm**2 mol).\n", + "Bond between 4 and 6 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", + "Bond between 4 and 5 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 4 and 0 has equilibrium length 0.1471 nm and force constant 307105.5999999999 kJ/(nm**2 mol).\n", + "Bond between 6 and 7 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 6 and 8 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 6 and 9 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 1 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", + "Bond between 2 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", + "Bond between 3 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", + "Bond between 10 and 12 has equilibrium length 0.1335 nm and force constant 410031.99999999994 kJ/(nm**2 mol).\n", + "Bond between 20 and 14 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", + "Bond between 20 and 21 has equilibrium length 0.12290000000000001 nm and force constant 476975.9999999999 kJ/(nm**2 mol).\n", + "Bond between 14 and 16 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", + "Bond between 14 and 15 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 14 and 12 has equilibrium length 0.1449 nm and force constant 282001.5999999999 kJ/(nm**2 mol).\n", + "Bond between 16 and 17 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 16 and 18 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 16 and 19 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 13 and 12 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", + "Bond between 20 and 22 has equilibrium length 0.1335 nm and force constant 410031.99999999994 kJ/(nm**2 mol).\n", + "Bond between 30 and 24 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", + "Bond between 30 and 31 has equilibrium length 0.125 nm and force constant 548940.7999999999 kJ/(nm**2 mol).\n", + "Bond between 30 and 32 has equilibrium length 0.125 nm and force constant 548940.7999999999 kJ/(nm**2 mol).\n", + "Bond between 24 and 26 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", + "Bond between 24 and 25 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 24 and 22 has equilibrium length 0.1449 nm and force constant 282001.5999999999 kJ/(nm**2 mol).\n", + "Bond between 26 and 27 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 26 and 28 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 26 and 29 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", + "Bond between 23 and 22 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n" + ] + } + ], + "source": [ + "for idx in range(harmonic_bond_force.getNumBonds()):\n", + " idx0, idx1, r0, k = harmonic_bond_force.getBondParameters(idx)\n", + " print(f\"Bond between {idx0} and {idx1} has equilibrium length {r0} and force constant {k}.\")" + ] + }, + { + "cell_type": "markdown", + "id": "9d7dbc26", + "metadata": {}, + "source": [ + "- angles:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "fa82f562", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Angle among 0, 4, and 5 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 0, 4, and 6 has equilibrium angle 1.9408061282176945 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 0, 4, and 10 has equilibrium angle 1.9408061282176945 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 1, 0, and 2 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 1, 0, and 3 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 1, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 2, 0, and 3 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 2, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 3, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 4, 6, and 7 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 4, 6, and 8 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 4, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 4, 10, and 11 has equilibrium angle 2.101376419401173 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 4, 10, and 12 has equilibrium angle 2.035053907825388 rad and force constant 585.76 kJ/(mol rad**2).\n", + "Angle among 5, 4, and 6 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 5, 4, and 10 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 6, 4, and 10 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", + "Angle among 7, 6, and 8 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 7, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 8, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 10, 12, and 13 has equilibrium angle 2.0943951023931953 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 10, 12, and 14 has equilibrium angle 2.1275563581810877 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 11, 10, and 12 has equilibrium angle 2.1450096507010312 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 12, 14, and 15 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 12, 14, and 16 has equilibrium angle 1.9146261894377796 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 12, 14, and 20 has equilibrium angle 1.9216075064457567 rad and force constant 527.184 kJ/(mol rad**2).\n", + "Angle among 13, 12, and 14 has equilibrium angle 2.0601866490541068 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 14, 16, and 17 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 14, 16, and 18 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 14, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 14, 20, and 21 has equilibrium angle 2.101376419401173 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 14, 20, and 22 has equilibrium angle 2.035053907825388 rad and force constant 585.76 kJ/(mol rad**2).\n", + "Angle among 15, 14, and 16 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 15, 14, and 20 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 16, 14, and 20 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", + "Angle among 17, 16, and 18 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 17, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 18, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 20, 22, and 23 has equilibrium angle 2.0943951023931953 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 20, 22, and 24 has equilibrium angle 2.1275563581810877 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 21, 20, and 22 has equilibrium angle 2.1450096507010312 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 22, 24, and 25 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 22, 24, and 26 has equilibrium angle 1.9146261894377796 rad and force constant 669.44 kJ/(mol rad**2).\n", + "Angle among 22, 24, and 30 has equilibrium angle 1.9216075064457567 rad and force constant 527.184 kJ/(mol rad**2).\n", + "Angle among 23, 22, and 24 has equilibrium angle 2.0601866490541068 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 24, 26, and 27 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 24, 26, and 28 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 24, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 24, 30, and 31 has equilibrium angle 2.0420352248333655 rad and force constant 585.76 kJ/(mol rad**2).\n", + "Angle among 24, 30, and 32 has equilibrium angle 2.0420352248333655 rad and force constant 585.76 kJ/(mol rad**2).\n", + "Angle among 25, 24, and 26 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 25, 24, and 30 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", + "Angle among 26, 24, and 30 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", + "Angle among 27, 26, and 28 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 27, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 28, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", + "Angle among 31, 30, and 32 has equilibrium angle 2.199114857512855 rad and force constant 669.44 kJ/(mol rad**2).\n" + ] + } + ], + "source": [ + "for idx in range(harmonic_angle_force.getNumAngles()):\n", + " idx0, idx1, idx2, r0, k = harmonic_angle_force.getAngleParameters(idx)\n", + " print(f\"Angle among {idx0}, {idx1}, and {idx2} has equilibrium angle {r0} and force constant {k}.\")" + ] + }, + { + "cell_type": "markdown", + "id": "01e07231", + "metadata": {}, + "source": [ + "Interestingly, one can also modify such parameters in-place.\n", + "In a machine learning framework, usually such modification is done gloabally.\n", + "We here only show how to modify a few bonds and angles, whereas global re-parametrization can be easily achieved in a similar manner with a larger `PyTree` holding all parameters." + ] + }, + { + "cell_type": "markdown", + "id": "ca4e40f9", + "metadata": {}, + "source": [ + "Say we would want to modify the bond parameters of the bond between atom 0 and 1." + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "cbdc43ec", + "metadata": {}, + "outputs": [], + "source": [ + "for idx in range(harmonic_bond_force.getNumBonds()):\n", + " idx0, idx1, r0, k = harmonic_bond_force.getBondParameters(idx)\n", + " if idx0 == 1 and idx1 == 0:\n", + " new_r0 = unit.Quantity(2666, unit.nanometer)\n", + " new_k = unit.Quantity(1984, unit.kilojoule / (unit.nanometer ** 2 * unit.mole))\n", + " harmonic_bond_force.setBondParameters(idx, idx0, idx1, new_r0, new_k)" + ] + }, + { + "cell_type": "markdown", + "id": "bb0f5768", + "metadata": {}, + "source": [ + "Now if you query again you will see that the bond parameter has been changed." + ] + }, + { + "cell_type": "markdown", + "id": "723167ec", + "metadata": {}, + "source": [ + "Similarly, for angles:" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "e2417613", + "metadata": {}, + "outputs": [], + "source": [ + "for idx in range(harmonic_angle_force.getNumAngles()):\n", + " idx0, idx1, idx2, r0, k = harmonic_angle_force.getAngleParameters(idx)\n", + " if idx0 == 1 and idx1 == 0 and idx2 == 2: \n", + " new_r0 = unit.Quantity(2666, unit.radian)\n", + " new_k = unit.Quantity(1984, unit.kilojoule / (unit.radian ** 2 * unit.mole))\n", + " harmonic_angle_force.setAngleParameters(idx, idx0, idx1, idx2, new_r0, new_k)" + ] + }, + { + "cell_type": "markdown", + "id": "3c2f89e5", + "metadata": {}, + "source": [ + "Note that we can use the `unit` module to control the units of the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "222d4f49", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/cookbook/Querying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying Charges and Other Parameters.ipynb index 6ad40e3..2a89907 100644 --- a/notebooks/cookbook/Querying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying Charges and Other Parameters.ipynb @@ -101,9 +101,8 @@ } ], "metadata": { - "tags": ["force field", "inspection", "forces"], "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -117,8 +116,13 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.1" - } + "version": "3.11.3" + }, + "tags": [ + "force field", + "inspection", + "forces" + ] }, "nbformat": 4, "nbformat_minor": 5 From e7b035c7764b7fa19eb51f5920380e0c9f3abf5f Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Wed, 19 Apr 2023 14:24:05 +0800 Subject: [PATCH 02/12] title change --- ...ystem.ipynb => Modifying Parameters in an OpenMM System.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename notebooks/cookbook/{Modifying Parameters in System.ipynb => Modifying Parameters in an OpenMM System.ipynb} (100%) diff --git a/notebooks/cookbook/Modifying Parameters in System.ipynb b/notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb similarity index 100% rename from notebooks/cookbook/Modifying Parameters in System.ipynb rename to notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb From 4266d19c22f54cfcf4a03f991dc3aef88928dfac Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Wed, 19 Apr 2023 19:06:57 +0800 Subject: [PATCH 03/12] start implementing own functional forms --- .DS_Store | Bin 0 -> 6148 bytes notebooks/cookbook/=0.2.0 | 0 notebooks/cookbook/=2.5 | 0 notebooks/cookbook/=4.10 | 56 ++++++++++ ...ur MM Functional Form Implementation.ipynb | 102 ++++++++++++++++++ 5 files changed, 158 insertions(+) create mode 100644 .DS_Store create mode 100644 notebooks/cookbook/=0.2.0 create mode 100644 notebooks/cookbook/=2.5 create mode 100644 notebooks/cookbook/=4.10 create mode 100644 notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..3044d3ec7b4713506bc4f7bef7b7297c7cfadc54 GIT binary patch literal 6148 zcmeHKIc~#13?vg54&118xnIZ+7J~c&eZUB8xJV!axoTB@m!Ib00mFivCP12CD3G&D zaaPb2q9`KTe(jz`8WCB-4drBEZg$^%W)GQBARK4B$Pw48e0=w_>iY@fKFI*qr|fU} z@59@1IK;{KsmxLVDnJFO02QDDzfiz>FKoFAWTXOAfC^j{u 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mopenff\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtoolkit\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtopology\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Molecule\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'openff'" + ] + } + ], + "source": [ + "import openmm\n", + "from openff.toolkit.topology import Molecule" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "3fd0e951", + "metadata": {}, + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'openff'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[2], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mopenff\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtoolkit\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtopology\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Molecule\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'openff'" + ] + } + ], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31d99d79", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 56819415ccbd795bc7c77f4da54d287b3583f038 Mon Sep 17 00:00:00 2001 From: yuanqing-wang Date: Wed, 19 Apr 2023 07:35:27 -0400 Subject: [PATCH 04/12] get energy components --- ...ur MM Functional Form Implementation.ipynb | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb b/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb index c94a57a..f90aff5 100644 --- a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb +++ b/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb @@ -28,7 +28,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "6347e55b", "metadata": {}, "outputs": [ @@ -49,6 +49,25 @@ "from openff.toolkit.topology import Molecule" ] }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d0deedd0", + "metadata": {}, + "source": [ + "### Get Energy Components with OpenMM" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "951dc8ab", + "metadata": {}, + "source": [ + "Querying the energy contributions from OpenMM is somewhat non-intuitive.\n", + "We would need to set different force groups first." + ] + }, { "cell_type": "code", "execution_count": 2, @@ -67,7 +86,10 @@ ] } ], - "source": [] + "source": [ + "for idx, force in enumerate(system.getForces()):\n", + " force.setForceGroup(idx)" + ] }, { "cell_type": "code", From a052c9752eba6815b2348d6f4e6d5921f40d1efd Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Wed, 19 Apr 2023 19:48:13 +0800 Subject: [PATCH 05/12] trash clean --- .DS_Store | Bin 6148 -> 6148 bytes notebooks/.DS_Store | Bin 0 -> 8196 bytes notebooks/cookbook/=0.2.0 | 0 notebooks/cookbook/=2.5 | 0 notebooks/cookbook/=4.10 | 56 -------------------------------------- 5 files changed, 56 deletions(-) create mode 100644 notebooks/.DS_Store delete mode 100644 notebooks/cookbook/=0.2.0 delete mode 100644 notebooks/cookbook/=2.5 delete mode 100644 notebooks/cookbook/=4.10 diff --git a/.DS_Store b/.DS_Store index 3044d3ec7b4713506bc4f7bef7b7297c7cfadc54..43d58cd6e91ad0af1d447c5805cba49adcb61cfa 100644 GIT binary patch delta 311 zcmZoMXfc=|#>B)qu~2NHo}wrd0|Nsi1A_nqLlQ$i5N0zJCzWs9xSVmafe1@ICqo`k ztb`#IS+=+!DJMS(sBcG7K~83IiGjg2MkZz!RyKAHb`EZi*x-!(^5BxhlG0+Q#G+^r zFC;%dCke(*ObW|PEsqxvan8>xNzBYkEdp!EOi2YQi3!ilOUW;H$}i1JDUOanlHuUw z;EWfLsIE4$G}TcsF)*vuQK+^wG61rT&1!2oIYgE9t%KsTb8_?YyMT@Z0!E<6zI0!J<G**W+*fB~^_<9FuC{35y{AX$)(1_%w-yg5W<1M|cN761$yO$7h| delta 88 zcmZoMXfc=|#>CJzu~2NHo}wrt0|NsP3otMwG2{bbHbZe)@WeuOM#jz0m=swzD=FPWS_>3pvRlh5SnReGq#(F#%S&28wxth>bss}Jb~;mcX16?2 zQ)7IaXnZC9B0hgMYBXXDzJ8hbil!lvL}JuL{NpeFU;=vX+*v|jKQzXG&fVO5=G=47 znK@^^bMKyA#uyU$Y&~ODj4_!yr&<{eGc+#eeOi-(nNp%4ea6zvq0eA?#2uctI-Ce2 z5Jn)3Kp25A0$~KMg$U4@&6_;Sxi6&Q7)BtB!2dD=e1C{i=QJMBF+u&Ug9g6@Aj(Ss zztJ_-0lrQ&pz(l?3F@noO;J4{P(@H;K&X>E%9j(32XstOq0S)G8G<__s8HbFo%Evq za)!8|;TT3BjKK5=@br;cM%3+B*Y9b|Nt5hUde|~uFBUrksibty+4T;*(x+86JGPfHZ6oI^m|26C^mdwd))`4Sc{}U+x=A6B z6A2 z+wIfJGbt;Rw^E+5-Lm*%np@~*a*k!~cU;pm9lN{BGmd&iLz6^ZJ+71Untfw^nwEEs z0|wbL3q}iBBvTdqZu^G08TX)(bBWf(W5P{ke8G~)iZ$yRZo56cV^7=N6AQKS3bj&Q z+*vSPGi@0iBWBhc=uWvA!_Jy^e-B02_RJ%ukyFZ1tF+sxA!DMVYU#3@mh1Y|)kP~^ z`Lua|-t>-XHL@)3QLWb~+PiXZmL9Ur$#RQXU!qs*t7Z8RpKs0?G>U0;tzN5WDY@OH z00|RyYQ3VlO8X$eK&UpV8)Yr4bXcjhVTq|Gb&IU+RXPj8b1_+|-XSYp8E0sSoN>im zLT{0^&*XdfjayT$w{O_6i;tAhcgu2*JCMp5`$hl8$ED1(RS;??Rwou8FS~ucRivHnhN(|%ZY)Prj;|VNJ@YU5fs}XGurH)`%ttNIZrHx>3 zsf|Q5iLzC&vHFOrl~67!zV3!vVi7P+n^jdSrK}XPZ91_mU=qz_hUb|4Pow=Cc9#9X z&a(^bBIcqD8X{POThU19-ih6W@K&@D!n@FeL+FKpewetA@b002qd11+cnA;UQ9Opn z3GL701-yupc$pCY8ex7MZ{jVyjSp}dAK_zs;$wdT-{VL8GzHCVlhKR|GoiVV%-NP> z_mixc1(VgTE3QO5Z56Nocg+6#|8?b4_^@FF!U+7e2%xkj*^(fg+4eMEYe%WyOPx2q wZcI>Lg$BNa&->*#(JOx#(seYD#REDfC`qXO>puki-LKr?`5&JD|Dp5$PaG Date: Wed, 19 Apr 2023 23:59:13 +0800 Subject: [PATCH 06/12] Update Use OpenMM to Verify Your MM Functional Form Implementation.ipynb --- ...ur MM Functional Form Implementation.ipynb | 500 +++++++++++++++++- 1 file changed, 473 insertions(+), 27 deletions(-) diff --git a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb b/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb index f90aff5..079bb4c 100644 --- a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb +++ b/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb @@ -13,7 +13,7 @@ "id": "024e7580", "metadata": {}, "source": [ - "Recently, there have been efforts where MM functional forms are implemented within machine learning packages using tensor-accelerating framework to assist the differentiable parametrization of force fields, such as [Espaloma](https://github.com/choderalab/espaloma) and [JAX-MD](https://github.com/jax-md/jax-md).\n", + "Recently, there have been efforts where MM functional forms are implemented within machine learning packages using tensor-accelerating framework to assist the differentiable parametrization of force fields, such as [Espaloma](https://github.com/choderalab/espaloma), [JAX-MD](https://github.com/jax-md/jax-md), and [TimeMachine](https://github.com/proteneer/timemachine/tree/master).\n", "When you implement your own MM functional forms, you might want to make sure that they match what are used in, say, OpenMM, exactly.\n", "Since the nonbonded part of the energy requires a lot more love to implement exactly, in this tutorial we first show you how to check your bonded energy implementation." ] @@ -23,34 +23,153 @@ "id": "7808e632", "metadata": {}, "source": [ - "For simplicity, we only work with small molecules for now." + "For simplicity, we only work with small molecules for now, which means that we would need to import `openff.toolkit` for chemoinformatics modeling." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 238, "id": "6347e55b", "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import openmm\n", + "from openff.toolkit.topology import Molecule" + ] + }, + { + "cell_type": "markdown", + "id": "23dbd437", + "metadata": {}, + "source": [ + "### Get a Toy Water Molecule System" + ] + }, + { + "cell_type": "markdown", + "id": "d9e3fa86", + "metadata": {}, + "source": [ + "We use a water molecule to focus first on bonds and angles." + ] + }, + { + "cell_type": "code", + "execution_count": 150, + "id": "079df807", + "metadata": {}, + "outputs": [], + "source": [ + "molecule = Molecule.from_smiles(\"[H]O[H]\")" + ] + }, + { + "cell_type": "code", + "execution_count": 151, + "id": "960825b2", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 151, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "molecule.visualize()" + ] + }, + { + "cell_type": "markdown", + "id": "f915510d", + "metadata": {}, + "source": [ + "It's also worth noting the index-element correspondance in this molecule, where the 0th and 2nd atoms are hydrogens whereas the 1st is oxygen." + ] + }, + { + "cell_type": "code", + "execution_count": 152, + "id": "38f4b37d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Atom(name=, atomic number=1),\n", + " Atom(name=, atomic number=8),\n", + " Atom(name=, atomic number=1)]" + ] + }, + "execution_count": 152, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "molecule.atoms" + ] + }, + { + "cell_type": "markdown", + "id": "a5f83c48", + "metadata": {}, + "source": [ + "To create OpenMM systems, we use an OpenFF force field." + ] + }, + { + "cell_type": "code", + "execution_count": 153, + "id": "6c5c757b", + "metadata": {}, + "outputs": [], + "source": [ + "from openff.toolkit.typing.engines.smirnoff import ForceField\n", + "forcefield = ForceField(\"openff_unconstrained-1.2.0.offxml\")" + ] + }, + { + "cell_type": "code", + "execution_count": 154, + "id": "e168c674", + "metadata": {}, "outputs": [ { - "ename": "ModuleNotFoundError", - "evalue": "No module named 'openff'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mopenmm\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mopenff\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtoolkit\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtopology\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Molecule\n", - "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'openff'" + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/wangy1/miniconda3/envs/malt/lib/python3.10/site-packages/openff/interchange/components/interchange.py:339: UserWarning: Automatically up-converting BondHandler from version 0.3 to 0.4. Consider manually upgrading this BondHandler (or section in an OFFXML file) to 0.4 or newer. For more details, see https://openforcefield.github.io/standards/standards/smirnoff/#bonds.\n", + " warnings.warn(\n" ] } ], "source": [ - "import openmm\n", - "from openff.toolkit.topology import Molecule" + "molecule.assign_partial_charges(\"mmff94\")\n", + "system = forcefield.create_openmm_system(molecule.to_topology(), charge_from_molecules=[molecule])\n", + "topology = molecule.to_topology().to_openmm()" ] }, { - "attachments": {}, "cell_type": "markdown", "id": "d0deedd0", "metadata": {}, @@ -59,7 +178,6 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "id": "951dc8ab", "metadata": {}, @@ -70,31 +188,359 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 155, "id": "3fd0e951", "metadata": {}, + "outputs": [], + "source": [ + "for idx, force in enumerate(system.getForces()):\n", + " force.setForceGroup(idx)" + ] + }, + { + "cell_type": "markdown", + "id": "af83459f", + "metadata": {}, + "source": [ + "Next, we create a simulation so that we can inspect each energy component." + ] + }, + { + "cell_type": "code", + "execution_count": 156, + "id": "31d99d79", + "metadata": {}, + "outputs": [], + "source": [ + "integrator = openmm.VerletIntegrator(0.0)\n", + "simulation = openmm.app.Simulation(topology,system, integrator)" + ] + }, + { + "cell_type": "markdown", + "id": "bab7fa2a", + "metadata": {}, + "source": [ + "Let's also set the initial positions to one of the conformations of the water molecule." + ] + }, + { + "cell_type": "code", + "execution_count": 157, + "id": "b320c078", + "metadata": {}, + "outputs": [], + "source": [ + "molecule.generate_conformers() # generate molecule conformers\n", + "\n", + "# NOTE:\n", + "# OpenMM and OpenFF use two sets of unit systems,\n", + "# hence the awkward conversion\n", + "simulation.context.setPositions(\n", + " openmm.unit.nanometer * molecule._conformers[0].magnitude,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 162, + "id": "65898ea4", + "metadata": {}, "outputs": [ { - "ename": "ModuleNotFoundError", - "evalue": "No module named 'openff'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[2], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mopenff\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtoolkit\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtopology\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Molecule\n", - "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'openff'" + "name": "stdout", + "output_type": "stream", + "text": [ + "Bond Energy 368275.4375 kJ/mol\n", + "Angle Energy 1.3650074005126953 kJ/mol\n" ] } ], "source": [ "for idx, force in enumerate(system.getForces()):\n", - " force.setForceGroup(idx)" + " state = simulation.context.getState(getEnergy=True, groups={idx})\n", + " if \"HarmonicBondForce\" in force.getName():\n", + " bond_energy = state.getPotentialEnergy()\n", + " if \"HarmonicAngleForce\" in force.getName():\n", + " angle_energy = state.getPotentialEnergy()\n", + "\n", + "print(\"Bond Energy\", bond_energy)\n", + "print(\"Angle Energy\", angle_energy)" + ] + }, + { + "cell_type": "markdown", + "id": "a9555d75", + "metadata": {}, + "source": [ + "### Implement our own distance and angle functions\n", + "Here we show a very simple example.\n", + "This is usually done in an elaborate tensor-accelerating framework." + ] + }, + { + "cell_type": "markdown", + "id": "6d5d0f5e", + "metadata": {}, + "source": [ + "The bond length can be computed as the root of the squared difference between two bonded atoms." + ] + }, + { + "cell_type": "code", + "execution_count": 166, + "id": "b466f554", + "metadata": {}, + "outputs": [], + "source": [ + "def get_distance_vector(x0, x1):\n", + " return x1 - x0" + ] + }, + { + "cell_type": "code", + "execution_count": 194, + "id": "94e76b8c", + "metadata": {}, + "outputs": [], + "source": [ + "def get_norm(x01):\n", + " return np.linalg.norm(x01, ord=2, axis=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": 195, + "id": "abe40b29", + "metadata": {}, + "outputs": [], + "source": [ + "def get_distance(x0, x1):\n", + " return get_norm(get_distance_vector(x0, x1))" + ] + }, + { + "cell_type": "code", + "execution_count": 208, + "id": "cc038688", + "metadata": {}, + "outputs": [], + "source": [ + "def get_angle(x0, x1, x2):\n", + " x10 = get_distance_vector(x1, x0)\n", + " x12 = get_distance_vector(x1, x2)\n", + " c012 = np.arctan2(\n", + " get_norm(np.cross(x10, x12)),\n", + " (x10 * x12).sum(-1)\n", + " )\n", + " return c012" + ] + }, + { + "cell_type": "markdown", + "id": "638fa9f4", + "metadata": {}, + "source": [ + "### Implement our own MM forces" + ] + }, + { + "cell_type": "markdown", + "id": "f42fbbd0", + "metadata": {}, + "source": [ + "Let's implement harmonic bond and angle forces from scratch:" + ] + }, + { + "cell_type": "code", + "execution_count": 209, + "id": "25c71478", + "metadata": {}, + "outputs": [], + "source": [ + "def harmonic(x, k, r0):\n", + " return 0.5 * k * (x - r0) ** 2" + ] + }, + { + "cell_type": "markdown", + "id": "2f422194", + "metadata": {}, + "source": [ + "### Compute angle and bond energy" + ] + }, + { + "cell_type": "markdown", + "id": "2137744d", + "metadata": {}, + "source": [ + "One of the nice things about OpenMM is that it works with units, whereas most tensor-accelerating frameworks don't.\n", + "Here our crappy implementation doesn't take care of units either.\n", + "So we work with numeric values only---everything has the default of OpenMM units: `nanometer` for distances, `radian` for angles, and `kJ/mol` for energies." + ] + }, + { + "cell_type": "markdown", + "id": "ad357717", + "metadata": {}, + "source": [ + "We first compute the bonds and angles using the functions defined above." + ] + }, + { + "cell_type": "code", + "execution_count": 210, + "id": "e23a808d", + "metadata": {}, + "outputs": [], + "source": [ + "conformer = molecule._conformers[0].magnitude\n", + "x0, x1, x2 = conformer" + ] + }, + { + "cell_type": "code", + "execution_count": 227, + "id": "e6abf994", + "metadata": {}, + "outputs": [], + "source": [ + "r01 = get_distance(x0, x1)\n", + "r12 = get_distance(x1, x2)" + ] + }, + { + "cell_type": "code", + "execution_count": 214, + "id": "b3057adc", + "metadata": {}, + "outputs": [], + "source": [ + "c012 = get_angle(x0, x1, x2)" + ] + }, + { + "cell_type": "markdown", + "id": "7c6d91c2", + "metadata": {}, + "source": [ + "Then we grab the parameters: bond equilibrium length `r0`, bond force constant `k_r`, equilibrium angle `theta_0`, angle force constant `k_theta`, from the OpenMM `System`.\n", + "Note that here we know that there is only one kind of bond and one kind of angle, saving us the need for bookkeeping.\n", + "In real life one would need to carefully manage bond and angle types." + ] + }, + { + "cell_type": "code", + "execution_count": 230, + "id": "44ac59b5", + "metadata": {}, + "outputs": [], + "source": [ + "for force in system.getForces():\n", + " if \"HarmonicBond\" in force.getName():\n", + " _, __, r0, k_r = force.getBondParameters(0)\n", + " if \"HarmonicAngle\" in force.getName():\n", + " _, __, ___, theta0, k_theta = force.getAngleParameters(0)" + ] + }, + { + "cell_type": "markdown", + "id": "31117ca9", + "metadata": {}, + "source": [ + "We have to erase the units, unfortunately." + ] + }, + { + "cell_type": "code", + "execution_count": 231, + "id": "2769d77d", + "metadata": {}, + "outputs": [], + "source": [ + "r0, k_r, theta0, k_theta = r0._value, k_r._value, theta0._value, k_theta._value" + ] + }, + { + "cell_type": "code", + "execution_count": 239, + "id": "aadb3640", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "368275.4577602855" + ] + }, + "execution_count": 239, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bond_energy_computed = harmonic(r01, k_r, r0) + harmonic(r12, k_r, r0)\n", + "bond_energy_computed" + ] + }, + { + "cell_type": "code", + "execution_count": 240, + "id": "c753c122", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.3650163207594144" + ] + }, + "execution_count": 240, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "angle_energy_computed = harmonic(c012, k_theta, theta0)\n", + "angle_energy_computed" + ] + }, + { + "cell_type": "markdown", + "id": "1e4a33be", + "metadata": {}, + "source": [ + "### Check consistency" + ] + }, + { + "cell_type": "code", + "execution_count": 246, + "id": "e8080be0", + "metadata": {}, + "outputs": [], + "source": [ + "assert np.allclose(bond_energy_computed, bond_energy._value)\n", + "assert np.allclose(angle_energy_computed, angle_energy._value)" + ] + }, + { + "cell_type": "markdown", + "id": "839ea7ab", + "metadata": {}, + "source": [ + "### Postlude\n", + "Here we have only checked the _bond_ and _angle_ energy consistency for _one_ molecule with _one_ type of bond and angle each.\n", + "But with the help of careful index matching and bookkeeping with `PyTree` containers, we are able to scale things up dramatically to check the consistency for other terms among considerable chemical and conformational spaces.\n", + "Check out unit tests [here](https://github.com/choderalab/espaloma/blob/master/espaloma/mm/tests/test_openmm_consistency.py) to see how this is done." ] }, { "cell_type": "code", "execution_count": null, - "id": "31d99d79", + "id": "08590304", "metadata": {}, "outputs": [], "source": [] @@ -116,7 +562,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.10.4" } }, "nbformat": 4, From 08af90856560c8734c4cb8839b0e57d7891054d0 Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Thu, 20 Apr 2023 00:01:22 +0800 Subject: [PATCH 07/12] trash cleaning --- .DS_Store | Bin 6148 -> 0 bytes notebooks/.DS_Store | Bin 8196 -> 0 bytes 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 .DS_Store delete mode 100644 notebooks/.DS_Store diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index 43d58cd6e91ad0af1d447c5805cba49adcb61cfa..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHK%}T>S5Z-O8O(;SR3OxqA7K~L3#Y>3w1&ruHr6we3Xv~%-wTDv3SzpK}@p+ut z-H4?XJc-zuF#FBUPnP{Q>}DBb+?z&+j9H8^0~E1lLbFCNj=CZ>?Lp-79U&9&j94Dd zR3@qnP>#{K;@QyJJrh|GBw0L>1xXY`%I$TMMPlZOX%^+O)>i?;G|ZvZp3l4e zla77t^cEd^-tRhKAN3Xs)7aZTIK3E-(x*(kC^|WOc~W*XR`3eOnwdTOlPne414OIJ zs)CRhAO?tmO<_Rqe@1IlCQS1r28e+lGl2Vp4T|V$EEUSF0~-8&#CQV{1#G-aAle#T zjio|}fN)g`s7krMVsKRsep~0c8cT($oN+layho4B^$mr~(ZO%abjDqU)Di>4z%m1M z)wS^aKmY#zzg$E;Vt^RX10P`|O_=}y diff --git a/notebooks/.DS_Store b/notebooks/.DS_Store deleted file mode 100644 index 733a4423f801be68384fe45b4044adeaa046ebf2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeHMYitx%6u#fIz>FPWS_>3pvRlh5SnReGq#(F#%S&28wxth>bss}Jb~;mcX16?2 zQ)7IaXnZC9B0hgMYBXXDzJ8hbil!lvL}JuL{NpeFU;=vX+*v|jKQzXG&fVO5=G=47 znK@^^bMKyA#uyU$Y&~ODj4_!yr&<{eGc+#eeOi-(nNp%4ea6zvq0eA?#2uctI-Ce2 z5Jn)3Kp25A0$~KMg$U4@&6_;Sxi6&Q7)BtB!2dD=e1C{i=QJMBF+u&Ug9g6@Aj(Ss zztJ_-0lrQ&pz(l?3F@noO;J4{P(@H;K&X>E%9j(32XstOq0S)G8G<__s8HbFo%Evq za)!8|;TT3BjKK5=@br;cM%3+B*Y9b|Nt5hUde|~uFBUrksibty+4T;*(x+86JGPfHZ6oI^m|26C^mdwd))`4Sc{}U+x=A6B z6A2 z+wIfJGbt;Rw^E+5-Lm*%np@~*a*k!~cU;pm9lN{BGmd&iLz6^ZJ+71Untfw^nwEEs z0|wbL3q}iBBvTdqZu^G08TX)(bBWf(W5P{ke8G~)iZ$yRZo56cV^7=N6AQKS3bj&Q z+*vSPGi@0iBWBhc=uWvA!_Jy^e-B02_RJ%ukyFZ1tF+sxA!DMVYU#3@mh1Y|)kP~^ z`Lua|-t>-XHL@)3QLWb~+PiXZmL9Ur$#RQXU!qs*t7Z8RpKs0?G>U0;tzN5WDY@OH z00|RyYQ3VlO8X$eK&UpV8)Yr4bXcjhVTq|Gb&IU+RXPj8b1_+|-XSYp8E0sSoN>im zLT{0^&*XdfjayT$w{O_6i;tAhcgu2*JCMp5`$hl8$ED1(RS;??Rwou8FS~ucRivHnhN(|%ZY)Prj;|VNJ@YU5fs}XGurH)`%ttNIZrHx>3 zsf|Q5iLzC&vHFOrl~67!zV3!vVi7P+n^jdSrK}XPZ91_mU=qz_hUb|4Pow=Cc9#9X z&a(^bBIcqD8X{POThU19-ih6W@K&@D!n@FeL+FKpewetA@b002qd11+cnA;UQ9Opn z3GL701-yupc$pCY8ex7MZ{jVyjSp}dAK_zs;$wdT-{VL8GzHCVlhKR|GoiVV%-NP> z_mixc1(VgTE3QO5Z56Nocg+6#|8?b4_^@FF!U+7e2%xkj*^(fg+4eMEYe%WyOPx2q wZcI>Lg$BNa&->*#(JOx#(seYD#REDfC`qXO>puki-LKr?`5&JD|Dp5$PaG Date: Thu, 20 Apr 2023 10:16:03 +0800 Subject: [PATCH 08/12] merged modifying parameters into querying parameters --- ...fying Parameters in an OpenMM System.ipynb | 322 ------------------ ...ifying Charges and Other Parameters.ipynb} | 50 ++- 2 files changed, 49 insertions(+), 323 deletions(-) delete mode 100644 notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb rename notebooks/cookbook/{Querying Charges and Other Parameters.ipynb => Querying and Modifying Charges and Other Parameters.ipynb} (75%) diff --git a/notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb b/notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb deleted file mode 100644 index 1be8b13..0000000 --- a/notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb +++ /dev/null @@ -1,322 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "5cc09957", - "metadata": {}, - "source": [ - "## Modifying Parameters in System\n", - "\n", - "Some machine learning force field construction schemes, such as [Espaloma](https://pubs.rsc.org/en/content/articlehtml/2022/sc/d2sc02739a), would output parameters which can be picked up by the parameter-holding `System` object of OpenMM, which could be further used in molecular dynamics (MD) simulations just like a vanilla `System`.\n", - "In this tutorial, we show how you can modify the parameters of `System` object so you can write your own machine learning force field API.\n", - "\n", - "This example is inspired by the [Espaloma deployment](https://github.com/choderalab/espaloma/blob/master/espaloma/graphs/deploy.py) implementation." - ] - }, - { - "cell_type": "markdown", - "id": "5003093c", - "metadata": {}, - "source": [ - "First, let's load a PDB file and use Amber14 force field as placeholder parameters." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "9c163ab2", - "metadata": {}, - "outputs": [], - "source": [ - "from openmm.app import *\n", - "from openmm import *\n", - "\n", - "pdb = PDBFile('ala_ala_ala.pdb')\n", - "forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')\n", - "system = forcefield.createSystem(pdb.topology)" - ] - }, - { - "cell_type": "markdown", - "id": "0c5c2c35", - "metadata": {}, - "source": [ - "Now we loop through the forces to get bonded components of the energy." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "f92baca2", - "metadata": {}, - "outputs": [], - "source": [ - "forces = list(system.getForces())" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "e8b4bef3", - "metadata": {}, - "outputs": [], - "source": [ - "harmonic_bond_force = [force for force in forces if \"HarmonicBondForce\" in force.__class__.__name__].pop()\n", - "harmonic_angle_force = [force for force in forces if \"HarmonicAngleForce\" in force.__class__.__name__].pop()\n", - "periodic_torsion_force = [force for force in forces if \"PeriodicTorsionForce\" in force.__class__.__name__].pop()" - ] - }, - { - "cell_type": "markdown", - "id": "0d8cde7e", - "metadata": {}, - "source": [ - "We can inspect each component to get a list of indices and parameters." - ] - }, - { - "cell_type": "markdown", - "id": "425b070e", - "metadata": {}, - "source": [ - "- Bonds:" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "b09324e6", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Bond between 10 and 4 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", - "Bond between 10 and 11 has equilibrium length 0.12290000000000001 nm and force constant 476975.9999999999 kJ/(nm**2 mol).\n", - "Bond between 4 and 6 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", - "Bond between 4 and 5 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 4 and 0 has equilibrium length 0.1471 nm and force constant 307105.5999999999 kJ/(nm**2 mol).\n", - "Bond between 6 and 7 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 6 and 8 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 6 and 9 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 1 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", - "Bond between 2 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", - "Bond between 3 and 0 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", - "Bond between 10 and 12 has equilibrium length 0.1335 nm and force constant 410031.99999999994 kJ/(nm**2 mol).\n", - "Bond between 20 and 14 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", - "Bond between 20 and 21 has equilibrium length 0.12290000000000001 nm and force constant 476975.9999999999 kJ/(nm**2 mol).\n", - "Bond between 14 and 16 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", - "Bond between 14 and 15 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 14 and 12 has equilibrium length 0.1449 nm and force constant 282001.5999999999 kJ/(nm**2 mol).\n", - "Bond between 16 and 17 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 16 and 18 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 16 and 19 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 13 and 12 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n", - "Bond between 20 and 22 has equilibrium length 0.1335 nm and force constant 410031.99999999994 kJ/(nm**2 mol).\n", - "Bond between 30 and 24 has equilibrium length 0.1522 nm and force constant 265265.6 kJ/(nm**2 mol).\n", - "Bond between 30 and 31 has equilibrium length 0.125 nm and force constant 548940.7999999999 kJ/(nm**2 mol).\n", - "Bond between 30 and 32 has equilibrium length 0.125 nm and force constant 548940.7999999999 kJ/(nm**2 mol).\n", - "Bond between 24 and 26 has equilibrium length 0.1526 nm and force constant 259407.99999999994 kJ/(nm**2 mol).\n", - "Bond between 24 and 25 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 24 and 22 has equilibrium length 0.1449 nm and force constant 282001.5999999999 kJ/(nm**2 mol).\n", - "Bond between 26 and 27 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 26 and 28 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 26 and 29 has equilibrium length 0.10900000000000001 nm and force constant 284511.99999999994 kJ/(nm**2 mol).\n", - "Bond between 23 and 22 has equilibrium length 0.101 nm and force constant 363171.19999999995 kJ/(nm**2 mol).\n" - ] - } - ], - "source": [ - "for idx in range(harmonic_bond_force.getNumBonds()):\n", - " idx0, idx1, r0, k = harmonic_bond_force.getBondParameters(idx)\n", - " print(f\"Bond between {idx0} and {idx1} has equilibrium length {r0} and force constant {k}.\")" - ] - }, - { - "cell_type": "markdown", - "id": "9d7dbc26", - "metadata": {}, - "source": [ - "- angles:" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "id": "fa82f562", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Angle among 0, 4, and 5 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 0, 4, and 6 has equilibrium angle 1.9408061282176945 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 0, 4, and 10 has equilibrium angle 1.9408061282176945 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 1, 0, and 2 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 1, 0, and 3 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 1, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 2, 0, and 3 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 2, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 3, 0, and 4 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 4, 6, and 7 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 4, 6, and 8 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 4, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 4, 10, and 11 has equilibrium angle 2.101376419401173 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 4, 10, and 12 has equilibrium angle 2.035053907825388 rad and force constant 585.76 kJ/(mol rad**2).\n", - "Angle among 5, 4, and 6 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 5, 4, and 10 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 6, 4, and 10 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", - "Angle among 7, 6, and 8 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 7, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 8, 6, and 9 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 10, 12, and 13 has equilibrium angle 2.0943951023931953 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 10, 12, and 14 has equilibrium angle 2.1275563581810877 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 11, 10, and 12 has equilibrium angle 2.1450096507010312 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 12, 14, and 15 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 12, 14, and 16 has equilibrium angle 1.9146261894377796 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 12, 14, and 20 has equilibrium angle 1.9216075064457567 rad and force constant 527.184 kJ/(mol rad**2).\n", - "Angle among 13, 12, and 14 has equilibrium angle 2.0601866490541068 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 14, 16, and 17 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 14, 16, and 18 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 14, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 14, 20, and 21 has equilibrium angle 2.101376419401173 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 14, 20, and 22 has equilibrium angle 2.035053907825388 rad and force constant 585.76 kJ/(mol rad**2).\n", - "Angle among 15, 14, and 16 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 15, 14, and 20 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 16, 14, and 20 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", - "Angle among 17, 16, and 18 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 17, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 18, 16, and 19 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 20, 22, and 23 has equilibrium angle 2.0943951023931953 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 20, 22, and 24 has equilibrium angle 2.1275563581810877 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 21, 20, and 22 has equilibrium angle 2.1450096507010312 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 22, 24, and 25 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 22, 24, and 26 has equilibrium angle 1.9146261894377796 rad and force constant 669.44 kJ/(mol rad**2).\n", - "Angle among 22, 24, and 30 has equilibrium angle 1.9216075064457567 rad and force constant 527.184 kJ/(mol rad**2).\n", - "Angle among 23, 22, and 24 has equilibrium angle 2.0601866490541068 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 24, 26, and 27 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 24, 26, and 28 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 24, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 24, 30, and 31 has equilibrium angle 2.0420352248333655 rad and force constant 585.76 kJ/(mol rad**2).\n", - "Angle among 24, 30, and 32 has equilibrium angle 2.0420352248333655 rad and force constant 585.76 kJ/(mol rad**2).\n", - "Angle among 25, 24, and 26 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 25, 24, and 30 has equilibrium angle 1.911135530933791 rad and force constant 418.40000000000003 kJ/(mol rad**2).\n", - "Angle among 26, 24, and 30 has equilibrium angle 1.9390607989657 rad and force constant 527.184 kJ/(mol rad**2).\n", - "Angle among 27, 26, and 28 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 27, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 28, 26, and 29 has equilibrium angle 1.911135530933791 rad and force constant 292.88 kJ/(mol rad**2).\n", - "Angle among 31, 30, and 32 has equilibrium angle 2.199114857512855 rad and force constant 669.44 kJ/(mol rad**2).\n" - ] - } - ], - "source": [ - "for idx in range(harmonic_angle_force.getNumAngles()):\n", - " idx0, idx1, idx2, r0, k = harmonic_angle_force.getAngleParameters(idx)\n", - " print(f\"Angle among {idx0}, {idx1}, and {idx2} has equilibrium angle {r0} and force constant {k}.\")" - ] - }, - { - "cell_type": "markdown", - "id": "01e07231", - "metadata": {}, - "source": [ - "Interestingly, one can also modify such parameters in-place.\n", - "In a machine learning framework, usually such modification is done gloabally.\n", - "We here only show how to modify a few bonds and angles, whereas global re-parametrization can be easily achieved in a similar manner with a larger `PyTree` holding all parameters." - ] - }, - { - "cell_type": "markdown", - "id": "ca4e40f9", - "metadata": {}, - "source": [ - "Say we would want to modify the bond parameters of the bond between atom 0 and 1." - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "id": "cbdc43ec", - "metadata": {}, - "outputs": [], - "source": [ - "for idx in range(harmonic_bond_force.getNumBonds()):\n", - " idx0, idx1, r0, k = harmonic_bond_force.getBondParameters(idx)\n", - " if idx0 == 1 and idx1 == 0:\n", - " new_r0 = unit.Quantity(2666, unit.nanometer)\n", - " new_k = unit.Quantity(1984, unit.kilojoule / (unit.nanometer ** 2 * unit.mole))\n", - " harmonic_bond_force.setBondParameters(idx, idx0, idx1, new_r0, new_k)" - ] - }, - { - "cell_type": "markdown", - "id": "bb0f5768", - "metadata": {}, - "source": [ - "Now if you query again you will see that the bond parameter has been changed." - ] - }, - { - "cell_type": "markdown", - "id": "723167ec", - "metadata": {}, - "source": [ - "Similarly, for angles:" - ] - }, - { - "cell_type": "code", - "execution_count": 34, - "id": "e2417613", - "metadata": {}, - "outputs": [], - "source": [ - "for idx in range(harmonic_angle_force.getNumAngles()):\n", - " idx0, idx1, idx2, r0, k = harmonic_angle_force.getAngleParameters(idx)\n", - " if idx0 == 1 and idx1 == 0 and idx2 == 2: \n", - " new_r0 = unit.Quantity(2666, unit.radian)\n", - " new_k = unit.Quantity(1984, unit.kilojoule / (unit.radian ** 2 * unit.mole))\n", - " harmonic_angle_force.setAngleParameters(idx, idx0, idx1, idx2, new_r0, new_k)" - ] - }, - { - "cell_type": "markdown", - "id": "3c2f89e5", - "metadata": {}, - "source": [ - "Note that we can use the `unit` module to control the units of the parameters." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "222d4f49", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.3" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/notebooks/cookbook/Querying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb similarity index 75% rename from notebooks/cookbook/Querying Charges and Other Parameters.ipynb rename to notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb index 2a89907..e74e031 100644 --- a/notebooks/cookbook/Querying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb @@ -5,7 +5,7 @@ "id": "3a4fdeed", "metadata": {}, "source": [ - "## Querying Charges and Other Parameters\n", + "## Querying and Modifyinig Charges and Other Parameters\n", "\n", "Sometimes you want to inspect the charges or other parameters of the particles or bonds in a System. Force field parameters are stored in the Force objects added to a System. As an example, let's load a PDB file and model it using the Amber14 force field." ] @@ -98,6 +98,54 @@ " if particle1 == 0 or particle2 == 0:\n", " print(f'Particles ({particle1}, {particle2}), length = {length}, k = {k}')" ] + }, + { + "cell_type": "markdown", + "id": "50f5157c", + "metadata": {}, + "source": [ + "In some cases you may want to modify those parameters to make a bond behave differently from what the `ForceField` assigned.\n", + "We show you here that to modify these parameters are also easy." + ] + }, + { + "cell_type": "markdown", + "id": "cd5396c5", + "metadata": {}, + "source": [ + "Let's say we would want to modify the force constant of the bond connecting particle 1 and 0 to a new number `_k`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "0dc0db66", + "metadata": {}, + "outputs": [], + "source": [ + "k_ = 2666 * unit.kilojoule / ( unit.nanometer ** 2 * unit.mole )" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "086c7629", + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(bonded.getNumBonds()):\n", + " particle1, particle2, length, k = bonded.getBondParameters(i)\n", + " if particle1 == 1 and particle2 == 0:\n", + " bonded.setBondParameters(i, particle1, particle2, length, k_)" + ] + }, + { + "cell_type": "markdown", + "id": "8f9728ad", + "metadata": {}, + "source": [ + "Now if you query again you can see new parameters." + ] } ], "metadata": { From 9c47ba769fdc6b9fc0adfbd2b1560ce21fb7efac Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Tue, 25 Apr 2023 12:21:34 +0800 Subject: [PATCH 09/12] moved benchmarking script into another PR; fixed typo; added querying again --- ...difying Charges and Other Parameters.ipynb | 24 +- ...ur MM Functional Form Implementation.ipynb | 570 ------------------ 2 files changed, 23 insertions(+), 571 deletions(-) delete mode 100644 notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb diff --git a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb index e74e031..319cebc 100644 --- a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb @@ -5,7 +5,7 @@ "id": "3a4fdeed", "metadata": {}, "source": [ - "## Querying and Modifyinig Charges and Other Parameters\n", + "## Querying and Modifying Charges and Other Parameters\n", "\n", "Sometimes you want to inspect the charges or other parameters of the particles or bonds in a System. Force field parameters are stored in the Force objects added to a System. As an example, let's load a PDB file and model it using the Amber14 force field." ] @@ -146,6 +146,28 @@ "source": [ "Now if you query again you can see new parameters." ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "68498ff5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Particles (1, 0), length = 0.101 nm, k = 2666.0 kJ/(nm**2 mol)\n" + ] + } + ], + "source": [ + "bonded = [f for f in system.getForces() if isinstance(f, HarmonicBondForce)][0]\n", + "for i in range(bonded.getNumBonds()):\n", + " particle1, particle2, length, k = bonded.getBondParameters(i)\n", + " if particle1 == 1 and particle2 == 0:\n", + " print(f'Particles ({particle1}, {particle2}), length = {length}, k = {k}')" + ] } ], "metadata": { diff --git a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb b/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb deleted file mode 100644 index 079bb4c..0000000 --- a/notebooks/cookbook/Use OpenMM to Verify Your MM Functional Form Implementation.ipynb +++ /dev/null @@ -1,570 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "62ee8371", - "metadata": {}, - "source": [ - "## Use OpenMM to Verify Your MM Functional Form Implementation" - ] - }, - { - "cell_type": "markdown", - "id": "024e7580", - "metadata": {}, - "source": [ - "Recently, there have been efforts where MM functional forms are implemented within machine learning packages using tensor-accelerating framework to assist the differentiable parametrization of force fields, such as [Espaloma](https://github.com/choderalab/espaloma), [JAX-MD](https://github.com/jax-md/jax-md), and [TimeMachine](https://github.com/proteneer/timemachine/tree/master).\n", - "When you implement your own MM functional forms, you might want to make sure that they match what are used in, say, OpenMM, exactly.\n", - "Since the nonbonded part of the energy requires a lot more love to implement exactly, in this tutorial we first show you how to check your bonded energy implementation." - ] - }, - { - "cell_type": "markdown", - "id": "7808e632", - "metadata": {}, - "source": [ - "For simplicity, we only work with small molecules for now, which means that we would need to import `openff.toolkit` for chemoinformatics modeling." - ] - }, - { - "cell_type": "code", - "execution_count": 238, - "id": "6347e55b", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import openmm\n", - "from openff.toolkit.topology import Molecule" - ] - }, - { - "cell_type": "markdown", - "id": "23dbd437", - "metadata": {}, - "source": [ - "### Get a Toy Water Molecule System" - ] - }, - { - "cell_type": "markdown", - "id": "d9e3fa86", - "metadata": {}, - "source": [ - "We use a water molecule to focus first on bonds and angles." - ] - }, - { - "cell_type": "code", - "execution_count": 150, - "id": "079df807", - "metadata": {}, - "outputs": [], - "source": [ - "molecule = Molecule.from_smiles(\"[H]O[H]\")" - ] - }, - { - "cell_type": "code", - "execution_count": 151, - "id": "960825b2", - "metadata": {}, - "outputs": [ - { - "data": { - "image/svg+xml": [ - "\n", - "\n", - " \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 151, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "molecule.visualize()" - ] - }, - { - "cell_type": "markdown", - "id": "f915510d", - "metadata": {}, - "source": [ - "It's also worth noting the index-element correspondance in this molecule, where the 0th and 2nd atoms are hydrogens whereas the 1st is oxygen." - ] - }, - { - "cell_type": "code", - "execution_count": 152, - "id": "38f4b37d", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[Atom(name=, atomic number=1),\n", - " Atom(name=, atomic number=8),\n", - " Atom(name=, atomic number=1)]" - ] - }, - "execution_count": 152, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "molecule.atoms" - ] - }, - { - "cell_type": "markdown", - "id": "a5f83c48", - "metadata": {}, - "source": [ - "To create OpenMM systems, we use an OpenFF force field." - ] - }, - { - "cell_type": "code", - "execution_count": 153, - "id": "6c5c757b", - "metadata": {}, - "outputs": [], - "source": [ - "from openff.toolkit.typing.engines.smirnoff import ForceField\n", - "forcefield = ForceField(\"openff_unconstrained-1.2.0.offxml\")" - ] - }, - { - "cell_type": "code", - "execution_count": 154, - "id": "e168c674", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/wangy1/miniconda3/envs/malt/lib/python3.10/site-packages/openff/interchange/components/interchange.py:339: UserWarning: Automatically up-converting BondHandler from version 0.3 to 0.4. Consider manually upgrading this BondHandler (or section in an OFFXML file) to 0.4 or newer. For more details, see https://openforcefield.github.io/standards/standards/smirnoff/#bonds.\n", - " warnings.warn(\n" - ] - } - ], - "source": [ - "molecule.assign_partial_charges(\"mmff94\")\n", - "system = forcefield.create_openmm_system(molecule.to_topology(), charge_from_molecules=[molecule])\n", - "topology = molecule.to_topology().to_openmm()" - ] - }, - { - "cell_type": "markdown", - "id": "d0deedd0", - "metadata": {}, - "source": [ - "### Get Energy Components with OpenMM" - ] - }, - { - "cell_type": "markdown", - "id": "951dc8ab", - "metadata": {}, - "source": [ - "Querying the energy contributions from OpenMM is somewhat non-intuitive.\n", - "We would need to set different force groups first." - ] - }, - { - "cell_type": "code", - "execution_count": 155, - "id": "3fd0e951", - "metadata": {}, - "outputs": [], - "source": [ - "for idx, force in enumerate(system.getForces()):\n", - " force.setForceGroup(idx)" - ] - }, - { - "cell_type": "markdown", - "id": "af83459f", - "metadata": {}, - "source": [ - "Next, we create a simulation so that we can inspect each energy component." - ] - }, - { - "cell_type": "code", - "execution_count": 156, - "id": "31d99d79", - "metadata": {}, - "outputs": [], - "source": [ - "integrator = openmm.VerletIntegrator(0.0)\n", - "simulation = openmm.app.Simulation(topology,system, integrator)" - ] - }, - { - "cell_type": "markdown", - "id": "bab7fa2a", - "metadata": {}, - "source": [ - "Let's also set the initial positions to one of the conformations of the water molecule." - ] - }, - { - "cell_type": "code", - "execution_count": 157, - "id": "b320c078", - "metadata": {}, - "outputs": [], - "source": [ - "molecule.generate_conformers() # generate molecule conformers\n", - "\n", - "# NOTE:\n", - "# OpenMM and OpenFF use two sets of unit systems,\n", - "# hence the awkward conversion\n", - "simulation.context.setPositions(\n", - " openmm.unit.nanometer * molecule._conformers[0].magnitude,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 162, - "id": "65898ea4", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Bond Energy 368275.4375 kJ/mol\n", - "Angle Energy 1.3650074005126953 kJ/mol\n" - ] - } - ], - "source": [ - "for idx, force in enumerate(system.getForces()):\n", - " state = simulation.context.getState(getEnergy=True, groups={idx})\n", - " if \"HarmonicBondForce\" in force.getName():\n", - " bond_energy = state.getPotentialEnergy()\n", - " if \"HarmonicAngleForce\" in force.getName():\n", - " angle_energy = state.getPotentialEnergy()\n", - "\n", - "print(\"Bond Energy\", bond_energy)\n", - "print(\"Angle Energy\", angle_energy)" - ] - }, - { - "cell_type": "markdown", - "id": "a9555d75", - "metadata": {}, - "source": [ - "### Implement our own distance and angle functions\n", - "Here we show a very simple example.\n", - "This is usually done in an elaborate tensor-accelerating framework." - ] - }, - { - "cell_type": "markdown", - "id": "6d5d0f5e", - "metadata": {}, - "source": [ - "The bond length can be computed as the root of the squared difference between two bonded atoms." - ] - }, - { - "cell_type": "code", - "execution_count": 166, - "id": "b466f554", - "metadata": {}, - "outputs": [], - "source": [ - "def get_distance_vector(x0, x1):\n", - " return x1 - x0" - ] - }, - { - "cell_type": "code", - "execution_count": 194, - "id": "94e76b8c", - "metadata": {}, - "outputs": [], - "source": [ - "def get_norm(x01):\n", - " return np.linalg.norm(x01, ord=2, axis=-1)" - ] - }, - { - "cell_type": "code", - "execution_count": 195, - "id": "abe40b29", - "metadata": {}, - "outputs": [], - "source": [ - "def get_distance(x0, x1):\n", - " return get_norm(get_distance_vector(x0, x1))" - ] - }, - { - "cell_type": "code", - "execution_count": 208, - "id": "cc038688", - "metadata": {}, - "outputs": [], - "source": [ - "def get_angle(x0, x1, x2):\n", - " x10 = get_distance_vector(x1, x0)\n", - " x12 = get_distance_vector(x1, x2)\n", - " c012 = np.arctan2(\n", - " get_norm(np.cross(x10, x12)),\n", - " (x10 * x12).sum(-1)\n", - " )\n", - " return c012" - ] - }, - { - "cell_type": "markdown", - "id": "638fa9f4", - "metadata": {}, - "source": [ - "### Implement our own MM forces" - ] - }, - { - "cell_type": "markdown", - "id": "f42fbbd0", - "metadata": {}, - "source": [ - "Let's implement harmonic bond and angle forces from scratch:" - ] - }, - { - "cell_type": "code", - "execution_count": 209, - "id": "25c71478", - "metadata": {}, - "outputs": [], - "source": [ - "def harmonic(x, k, r0):\n", - " return 0.5 * k * (x - r0) ** 2" - ] - }, - { - "cell_type": "markdown", - "id": "2f422194", - "metadata": {}, - "source": [ - "### Compute angle and bond energy" - ] - }, - { - "cell_type": "markdown", - "id": "2137744d", - "metadata": {}, - "source": [ - "One of the nice things about OpenMM is that it works with units, whereas most tensor-accelerating frameworks don't.\n", - "Here our crappy implementation doesn't take care of units either.\n", - "So we work with numeric values only---everything has the default of OpenMM units: `nanometer` for distances, `radian` for angles, and `kJ/mol` for energies." - ] - }, - { - "cell_type": "markdown", - "id": "ad357717", - "metadata": {}, - "source": [ - "We first compute the bonds and angles using the functions defined above." - ] - }, - { - "cell_type": "code", - "execution_count": 210, - "id": "e23a808d", - "metadata": {}, - "outputs": [], - "source": [ - "conformer = molecule._conformers[0].magnitude\n", - "x0, x1, x2 = conformer" - ] - }, - { - "cell_type": "code", - "execution_count": 227, - "id": "e6abf994", - "metadata": {}, - "outputs": [], - "source": [ - "r01 = get_distance(x0, x1)\n", - "r12 = get_distance(x1, x2)" - ] - }, - { - "cell_type": "code", - "execution_count": 214, - "id": "b3057adc", - "metadata": {}, - "outputs": [], - "source": [ - "c012 = get_angle(x0, x1, x2)" - ] - }, - { - "cell_type": "markdown", - "id": "7c6d91c2", - "metadata": {}, - "source": [ - "Then we grab the parameters: bond equilibrium length `r0`, bond force constant `k_r`, equilibrium angle `theta_0`, angle force constant `k_theta`, from the OpenMM `System`.\n", - "Note that here we know that there is only one kind of bond and one kind of angle, saving us the need for bookkeeping.\n", - "In real life one would need to carefully manage bond and angle types." - ] - }, - { - "cell_type": "code", - "execution_count": 230, - "id": "44ac59b5", - "metadata": {}, - "outputs": [], - "source": [ - "for force in system.getForces():\n", - " if \"HarmonicBond\" in force.getName():\n", - " _, __, r0, k_r = force.getBondParameters(0)\n", - " if \"HarmonicAngle\" in force.getName():\n", - " _, __, ___, theta0, k_theta = force.getAngleParameters(0)" - ] - }, - { - "cell_type": "markdown", - "id": "31117ca9", - "metadata": {}, - "source": [ - "We have to erase the units, unfortunately." - ] - }, - { - "cell_type": "code", - "execution_count": 231, - "id": "2769d77d", - "metadata": {}, - "outputs": [], - "source": [ - "r0, k_r, theta0, k_theta = r0._value, k_r._value, theta0._value, k_theta._value" - ] - }, - { - "cell_type": "code", - "execution_count": 239, - "id": "aadb3640", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "368275.4577602855" - ] - }, - "execution_count": 239, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "bond_energy_computed = harmonic(r01, k_r, r0) + harmonic(r12, k_r, r0)\n", - "bond_energy_computed" - ] - }, - { - "cell_type": "code", - "execution_count": 240, - "id": "c753c122", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "1.3650163207594144" - ] - }, - "execution_count": 240, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "angle_energy_computed = harmonic(c012, k_theta, theta0)\n", - "angle_energy_computed" - ] - }, - { - "cell_type": "markdown", - "id": "1e4a33be", - "metadata": {}, - "source": [ - "### Check consistency" - ] - }, - { - "cell_type": "code", - "execution_count": 246, - "id": "e8080be0", - "metadata": {}, - "outputs": [], - "source": [ - "assert np.allclose(bond_energy_computed, bond_energy._value)\n", - "assert np.allclose(angle_energy_computed, angle_energy._value)" - ] - }, - { - "cell_type": "markdown", - "id": "839ea7ab", - "metadata": {}, - "source": [ - "### Postlude\n", - "Here we have only checked the _bond_ and _angle_ energy consistency for _one_ molecule with _one_ type of bond and angle each.\n", - "But with the help of careful index matching and bookkeeping with `PyTree` containers, we are able to scale things up dramatically to check the consistency for other terms among considerable chemical and conformational spaces.\n", - "Check out unit tests [here](https://github.com/choderalab/espaloma/blob/master/espaloma/mm/tests/test_openmm_consistency.py) to see how this is done." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "08590304", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.4" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From f82102948131b518625f5ef859f1e460b1756cda Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Wed, 26 Apr 2023 13:32:03 +0800 Subject: [PATCH 10/12] minor fixes; update parameters in `Context` --- ...difying Charges and Other Parameters.ipynb | 75 ++++++++++++++----- 1 file changed, 57 insertions(+), 18 deletions(-) diff --git a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb index 319cebc..af2e4ec 100644 --- a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "id": "9323aaae", "metadata": {}, "outputs": [], @@ -25,6 +25,28 @@ "system = forcefield.createSystem(pdb.topology)" ] }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "6d26c481", + "metadata": {}, + "source": [ + "Let's also create a `Simulation` with this system." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "552481f3", + "metadata": {}, + "outputs": [], + "source": [ + "integrator = LangevinMiddleIntegrator(\n", + " 300 * unit.kelvin, 1 / unit.picosecond, 0.004 * unit.picoseconds)\n", + "simulation = Simulation(pdb.topology, system, integrator)\n", + "simulation.context.setPositions(pdb.positions)" + ] + }, { "cell_type": "markdown", "id": "3c6c73f4", @@ -35,7 +57,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "id": "1734190f", "metadata": {}, "outputs": [], @@ -53,7 +75,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "id": "bf38fc9a", "metadata": {}, "outputs": [], @@ -76,7 +98,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 6, "id": "933928f9", "metadata": {}, "outputs": [ @@ -109,21 +131,12 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "id": "cd5396c5", "metadata": {}, "source": [ - "Let's say we would want to modify the force constant of the bond connecting particle 1 and 0 to a new number `_k`." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "0dc0db66", - "metadata": {}, - "outputs": [], - "source": [ - "k_ = 2666 * unit.kilojoule / ( unit.nanometer ** 2 * unit.mole )" + "Let's say we would want to modify the force constant of the bond connecting particle 1 and 0 to a new number." ] }, { @@ -136,7 +149,10 @@ "for i in range(bonded.getNumBonds()):\n", " particle1, particle2, length, k = bonded.getBondParameters(i)\n", " if particle1 == 1 and particle2 == 0:\n", - " bonded.setBondParameters(i, particle1, particle2, length, k_)" + " bonded.setBondParameters(\n", + " i, particle1, particle2, length, \n", + " 2666 * unit.kilojoules_per_mole/unit.nanometer**2\n", + " )" ] }, { @@ -157,7 +173,10 @@ "name": "stdout", "output_type": "stream", "text": [ - "Particles (1, 0), length = 0.101 nm, k = 2666.0 kJ/(nm**2 mol)\n" + "Particles (4, 0), length = 0.1471 nm, k = 307105.5999999999 kJ/(nm**2 mol)\n", + "Particles (1, 0), length = 0.101 nm, k = 2666.0 kJ/(nm**2 mol)\n", + "Particles (2, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)\n", + "Particles (3, 0), length = 0.101 nm, k = 363171.19999999995 kJ/(nm**2 mol)\n" ] } ], @@ -165,9 +184,29 @@ "bonded = [f for f in system.getForces() if isinstance(f, HarmonicBondForce)][0]\n", "for i in range(bonded.getNumBonds()):\n", " particle1, particle2, length, k = bonded.getBondParameters(i)\n", - " if particle1 == 1 and particle2 == 0:\n", + " if particle1 == 0 or particle2 == 0:\n", " print(f'Particles ({particle1}, {particle2}), length = {length}, k = {k}')" ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "47729e2e", + "metadata": {}, + "source": [ + "Modifying `Force` objects will affect any new `Simulations` or `Contexts` you create, but it will have no effect on ones that already exist. \n", + "If you want your modifications to apply to an existing `Simulation`, you can copy them over by calling" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "77075eed", + "metadata": {}, + "outputs": [], + "source": [ + "bonded.updateParametersInContext(simulation.context)" + ] } ], "metadata": { From eafe1dfab4c63a8fbbaa5dee81e61b7f7161a6af Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Fri, 28 Apr 2023 11:41:40 +0800 Subject: [PATCH 11/12] remove line char limit; get rid of `Context` --- ...difying Charges and Other Parameters.ipynb | 39 +------------------ 1 file changed, 2 insertions(+), 37 deletions(-) diff --git a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb index af2e4ec..a1cc1c8 100644 --- a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb @@ -25,28 +25,6 @@ "system = forcefield.createSystem(pdb.topology)" ] }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "6d26c481", - "metadata": {}, - "source": [ - "Let's also create a `Simulation` with this system." - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "552481f3", - "metadata": {}, - "outputs": [], - "source": [ - "integrator = LangevinMiddleIntegrator(\n", - " 300 * unit.kelvin, 1 / unit.picosecond, 0.004 * unit.picoseconds)\n", - "simulation = Simulation(pdb.topology, system, integrator)\n", - "simulation.context.setPositions(pdb.positions)" - ] - }, { "cell_type": "markdown", "id": "3c6c73f4", @@ -149,10 +127,7 @@ "for i in range(bonded.getNumBonds()):\n", " particle1, particle2, length, k = bonded.getBondParameters(i)\n", " if particle1 == 1 and particle2 == 0:\n", - " bonded.setBondParameters(\n", - " i, particle1, particle2, length, \n", - " 2666 * unit.kilojoules_per_mole/unit.nanometer**2\n", - " )" + " bonded.setBondParameters(i, particle1, particle2, length, 2666 * unit.kilojoules_per_mole/unit.nanometer**2)" ] }, { @@ -195,17 +170,7 @@ "metadata": {}, "source": [ "Modifying `Force` objects will affect any new `Simulations` or `Contexts` you create, but it will have no effect on ones that already exist. \n", - "If you want your modifications to apply to an existing `Simulation`, you can copy them over by calling" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "77075eed", - "metadata": {}, - "outputs": [], - "source": [ - "bonded.updateParametersInContext(simulation.context)" + "If you want your modifications to apply to an existing `Simulation`, you can copy them over by calling `bonded.updateParametersInContext(simulation.context)`." ] } ], From 4c9438a19af1c3eddd9417927e62eaf821a163c3 Mon Sep 17 00:00:00 2001 From: Yuanqing Wang Date: Mon, 1 May 2023 23:30:05 +0200 Subject: [PATCH 12/12] removed redundant line --- .../Querying and Modifying Charges and Other Parameters.ipynb | 1 - 1 file changed, 1 deletion(-) diff --git a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb index a1cc1c8..430c61f 100644 --- a/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb +++ b/notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb @@ -156,7 +156,6 @@ } ], "source": [ - "bonded = [f for f in system.getForces() if isinstance(f, HarmonicBondForce)][0]\n", "for i in range(bonded.getNumBonds()):\n", " particle1, particle2, length, k = bonded.getBondParameters(i)\n", " if particle1 == 0 or particle2 == 0:\n",