Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Not only querying, but also modifying system parameters #18

Merged
merged 14 commits into from
May 1, 2023
Merged
Show file tree
Hide file tree
Changes from 8 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
322 changes: 322 additions & 0 deletions notebooks/cookbook/Modifying Parameters in an OpenMM System.ipynb
Original file line number Diff line number Diff line change
@@ -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."
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This description is overly specific. Nothing in this notebook is really specific to machine learning. You're just showing how to examine and change the parameters in a System. How about something like this?

Sometimes it is useful to examine the parameters in a System. For example, you might want to check what parameters a ForceField has assigned to a particular bond or angle. In other cases you may want to modify those parameters to make a bond behave differently from what the ForceField assigned.

]
},
{
"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",
yuanqing-wang marked this conversation as resolved.
Show resolved Hide resolved
"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."
yuanqing-wang marked this conversation as resolved.
Show resolved Hide resolved
]
},
{
"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",
yuanqing-wang marked this conversation as resolved.
Show resolved Hide resolved
" 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
}
12 changes: 8 additions & 4 deletions notebooks/cookbook/Querying Charges and Other Parameters.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,8 @@
}
],
"metadata": {
"tags": ["force field", "inspection", "forces"],
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -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
Expand Down
Loading