Update for 2026

MM-GBSA.ipynb

@@ -12,9 +12,9 @@
     "\n",
     "For a good review, see: [Genheden and Ryde](https://www.tandfonline.com/doi/full/10.1517/17460441.2015.1032936#d1e147).\n",
     "\n",
-    "In summary, the MM-GBSA method can be used to post=process an MD simulation of a receptor-ligand complex and extract a predicted binding free energy. The steps are:\n",
+    "In summary, the MM-GBSA method can be used to post-process an MD simulation of a receptor-ligand complex and extract a predicted binding free energy. The steps are:\n",
     "\n",
-    "1. Three further simulation systems are parameterized, as if for implicit solvent simulations:\n",
+    "1. Starting with the original **explicitly-solvated** system for the receptor-ligand complex, three further simulation systems are **parameterized**, as if for **implicit solvent** simulations:\n",
     "    1. The complex (receptor + ligand, no explicit ions or solvent)\n",
     "    2. The receptor alone.\n",
     "    3. The ligand alone.\n",
@@ -29,10 +29,10 @@
     "\n",
     "![Figure 1](Figure1.png)\n",
     "\n",
-    "## MM-GBSA calculations using *Ambertools* and *OpenMM*\n",
+    "## MM-GBSA calculations using *ParmEd* and *OpenMM*\n",
     "\n",
     "This worshop demonstrates how an MM-GBSA analysis can be run in a Python notebook using a combination of [OpenMM](https://openmm.org) and [parmed](https://parmed.github.io/ParmEd/html/index.html). \n",
-    "The data comes from a short simulation of the complex between the protein Mcl-1 and ligand 5X3 (PDB code [5FDR](https://www.rcsb.org/structure/5FDR)), run with *OpenMM*. \n",
+    "The data comes from a short simulation of the complex between the protein Abl and the imatinib-analogue ligand that you generated earlier in the workshop. \n",
     "\n",
     "### 1. Load the required Python packages\n",
     "\n",
@@ -62,7 +62,7 @@
    "metadata": {},
    "source": [
     "### 2. Load the trajectory data\n",
-    "In the folder you will find `5fdr_A_solvated.prmtop`, the Amber-format *prmtop* file for the system, and `5fdr_A_solvated.xtc`, the trajectory file in the compressed *Gromacs* format. We load this using *MDTraj*, then generate two lists: one of the indices of the atoms that are part of the receptor, and one for atoms that are part of the ligand. From these two we can generate a list of the atoms that form the complex:"
+    "In the /data folder you will find `abl_ligand.prmtop`, the Amber-format *prmtop* file for the system, and `abl_ligand_prod.nc`, the trajectory file generated previously. We load this using *MDTraj*, then generate two lists: one of the indices of the atoms that are part of the receptor, and one for atoms that are part of the ligand. From these two we can generate a list of the atoms that form the complex:"
    ]
   },
   {
@@ -72,11 +72,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "t = mdt.load('5fdr_A_solvated.xtc', top='5fdr_A_solvated.prmtop')\n",
+    "t = mdt.load('abl_ligand_prod.nc', top='abl_ligand.prmtop')\n",
     "receptor_atoms = t.topology.select('protein')\n",
-    "ligand_atoms = t.topology.select('resname \"5X3\"')\n",
+    "ligand_atoms = t.topology.select('resname LIG') # or maybe 'resname UNL', depending on how you set up your own simulation\n",
     "\n",
-    "complex_atoms = np.sort(np.concatenate((receptor_atoms, ligand_atoms)))"
+    "complex_atoms = np.sort(np.concatenate((receptor_atoms, ligand_atoms)))\n",
+    "\n",
+    "print(f'There are {len(receptor_atoms)} atoms in the receptor and {len(ligand_atoms)} in the ligand.')"
    ]
   },
   {
@@ -86,7 +88,7 @@
    "source": [
     "### 3. Generate the implicit solvent simulation systems\n",
     "\n",
-    "We can use *parmed* to manipulate the original parameter/topology file, and create simulation systems for complex, receptor, and ligand. Each is modelled with an implicit solvent model (here the Hawkins, Cramer, Truhlar model, but others are possible), and as is usual for implicit solvent simulations, no nonbonded cutoff is used:"
+    "We can use *parmed* to manipulate the original parameter/topology file, and create simulation systems for complex, receptor, and ligand. Each is modelled with an implicit solvent model (here \"HCT\" - the Hawkins, Cramer, Truhlar model - but others are possible), and as is usual for implicit solvent simulations, no nonbonded cutoff is used:"
    ]
   },
   {
@@ -96,7 +98,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "prmtop = AmberParm('5fdr_A_solvated.prmtop')\n",
+    "prmtop = AmberParm('abl_ligand.prmtop')\n",
     "\n",
     "# Slice the prmtop file as required:\n",
     "c_prmtop = prmtop[complex_atoms]\n",
@@ -187,7 +189,7 @@
    "id": "05de0592-be10-4c58-a12a-52765ed749bb",
    "metadata": {},
    "source": [
-    "From these numbers, calculating $Ei$ (given by $Ec$ - ($Er$ + $El$)) should give a value of about -130 kJ/mol. If you have made an error defining which atoms belong to the receptor and/or ligand at the start of the process, you will most likely end up with nonsensical numbers instead.\n",
+    "From these numbers, calculating $Ei$ (given by $Ec$ - ($Er$ + $El$)) should give a value of about -250 kJ/mol. If you have made an error defining which atoms belong to the receptor and/or ligand at the start of the process, you will most likely end up with nonsensical numbers instead.\n",
     "\n",
     "### 6. Gather the MM-GBSA data for each snapshot\n",
     "\n",
@@ -251,7 +253,7 @@
    "source": [
     "### 8. Explore the data in more detail\n",
     "\n",
-    "To dig a bit deeper into the numbers, a graph can help. Let's plot the individual *Ei* values for each snapshot. Not that *Matplotlib* doesn't understand *OpenMM*'s concept of *Quantities* so we have to strip these off before plotting:"
+    "To dig a bit deeper into the numbers, a graph can help. Let's plot the individual *Ei* values for each snapshot. Note that *Matplotlib* doesn't understand *OpenMM*'s concept of *Quantities* so we have to divide the units away before plotting:"
    ]
   },
   {
@@ -273,7 +275,7 @@
    "id": "e1b60e96-478a-4d82-ba0c-92302e052d02",
    "metadata": {},
    "source": [
-    "We can see the first two or three points might be outliers - hardly surprising in this case as the data comes from a very short simulation of a system that begun from the crystal structure conformation and so was almost certainly not equilibrated/relaxed. A \"real\" MM-GBSA calculation needs an equilbrated, well-sampled, trajectory to work with.\n",
+    "We can see that, at least for this trajectory, the predicted binding energy fluctuates considerably, but after about 1 nanosecond of simulation it oscillates around a fairly steady mean value. \n",
     "\n",
     "## Summary\n",
     "\n",
@@ -285,7 +287,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "62274420-6aa7-4bee-9fa7-1ffbf5da2e6b",
+   "id": "7f627f65-12ae-4157-9221-45d3e5e5f4ad",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -307,7 +309,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.19"
+   "version": "3.12.11"
   }
  },
  "nbformat": 4,