From 6ae06e43ea383a892b3bf9e1f04079b49619b917 Mon Sep 17 00:00:00 2001 From: Stephen Farr Date: Thu, 2 Nov 2023 14:51:34 +0000 Subject: [PATCH 1/4] Add minimizationReporter example --- cookbook.md | 1 + notebooks/cookbook/report_minimization.ipynb | 129 +++++++++++++++++++ 2 files changed, 130 insertions(+) create mode 100644 notebooks/cookbook/report_minimization.ipynb diff --git a/cookbook.md b/cookbook.md index ad2372a..37997e7 100644 --- a/cookbook.md +++ b/cookbook.md @@ -51,6 +51,7 @@ caption: Analysis and System Inspection notebooks/cookbook/Analyzing Energy Contributions notebooks/cookbook/Computing Interaction Energies notebooks/cookbook/Querying and Modifying Charges and Other Parameters.ipynb +notebooks/cookbook/report_minimization.ipynb ::: diff --git a/notebooks/cookbook/report_minimization.ipynb b/notebooks/cookbook/report_minimization.ipynb new file mode 100644 index 0000000..d233381 --- /dev/null +++ b/notebooks/cookbook/report_minimization.ipynb @@ -0,0 +1,129 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e174828a", + "metadata": {}, + "source": [ + "## Reporting Minimization\n", + "\n", + "You can report the status of [`simulation.minimizeEnergy()`](http://docs.openmm.org/development/api-python/generated/openmm.app.simulation.Simulation.html?#openmm.app.simulation.Simulation.minimizeEnergy) using a [`MinimizerReporter`](http://docs.openmm.org/development/api-python/generated/openmm.openmm.MinimizationReporter.html).\n", + "\n", + "The syntax for doing this is somewhat different to typical OpenMM functionality. Read the [API documentation](http://docs.openmm.org/development/api-python/generated/openmm.openmm.MinimizationReporter.html) for more explanation.\n", + "\n", + "First we will create a test system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3592a17b", + "metadata": {}, + "outputs": [], + "source": [ + "from openmm.app import *\n", + "from openmm import *\n", + "from openmm.unit import *\n", + "from sys import stdout\n", + "\n", + "pdb = PDBFile('villin.pdb')\n", + "forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')\n", + "system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,\n", + " nonbondedCutoff=1*nanometer, constraints=HBonds)\n", + "integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)\n", + "simulation = Simulation(pdb.topology, system, integrator)\n", + "simulation.context.setPositions(pdb.positions)" + ] + }, + { + "cell_type": "markdown", + "id": "2e7b659d", + "metadata": {}, + "source": [ + "We then create a `MinimizationReporter` class that will print the current energy to the screen every 100 steps and save the energies to an array which we can plot later" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8692b40e", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a subclass of MinimizationReporter\n", + "class Reporter(MinimizationReporter):\n", + " interval = 100 # report interval\n", + " energies = [] # array to record progress\n", + " def report(self, iteration, x, grad, args):\n", + " # print current system energy to screen \n", + " if iteration % self.interval == 0:\n", + " print(iteration, args['system energy'])\n", + "\n", + " # save energy at each iteration to an array we can use later\n", + " self.energies.append(args['system energy'])\n", + "\n", + " # The report method must return a bool specifying if minimization should be stopped. \n", + " # You can use this functionality for early termination.\n", + " return False" + ] + }, + { + "cell_type": "markdown", + "id": "9bbc621a", + "metadata": {}, + "source": [ + "We now create an instance of the reporter and minimize the system with the reporter attached." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6be58c3b", + "metadata": {}, + "outputs": [], + "source": [ + "# Create an instance of our reporter\n", + "reporter = Reporter()\n", + "\n", + "# Perform local energy minimization \n", + "simulation.minimizeEnergy(reporter=reporter)\n", + "\n", + "# plot the minimization progress\n", + "import matplotlib as plt\n", + "plt.plot(Reporter.energies)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8" + }, + "tags": [ + "barostat", + "thermostat", + "application layer" + ], + "vscode": { + "interpreter": { + "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From af53b2c9e023761eb7d8c4c9d7802a0c1a7e4a50 Mon Sep 17 00:00:00 2001 From: Stephen Farr Date: Thu, 2 Nov 2023 15:01:34 +0000 Subject: [PATCH 2/4] Fix matplotlib import --- notebooks/cookbook/report_minimization.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/notebooks/cookbook/report_minimization.ipynb b/notebooks/cookbook/report_minimization.ipynb index d233381..6b369b6 100644 --- a/notebooks/cookbook/report_minimization.ipynb +++ b/notebooks/cookbook/report_minimization.ipynb @@ -89,7 +89,7 @@ "simulation.minimizeEnergy(reporter=reporter)\n", "\n", "# plot the minimization progress\n", - "import matplotlib as plt\n", + "import matplotlib.pyplot as plt\n", "plt.plot(Reporter.energies)\n", "plt.show()" ] @@ -111,7 +111,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.8" + "version": "3.12.0" }, "tags": [ "barostat", From 1f93f528efe694f31ce234098eae9440f22306d2 Mon Sep 17 00:00:00 2001 From: "S. Farr" Date: Thu, 28 Dec 2023 22:11:48 +0100 Subject: [PATCH 3/4] update minimization reporter --- notebooks/cookbook/report_minimization.ipynb | 87 ++++++++++++++++---- 1 file changed, 69 insertions(+), 18 deletions(-) diff --git a/notebooks/cookbook/report_minimization.ipynb b/notebooks/cookbook/report_minimization.ipynb index 6b369b6..e8d3b77 100644 --- a/notebooks/cookbook/report_minimization.ipynb +++ b/notebooks/cookbook/report_minimization.ipynb @@ -7,11 +7,11 @@ "source": [ "## Reporting Minimization\n", "\n", - "You can report the status of [`simulation.minimizeEnergy()`](http://docs.openmm.org/development/api-python/generated/openmm.app.simulation.Simulation.html?#openmm.app.simulation.Simulation.minimizeEnergy) using a [`MinimizerReporter`](http://docs.openmm.org/development/api-python/generated/openmm.openmm.MinimizationReporter.html).\n", + "You can report the status of [simulation.minimizeEnergy()](http://docs.openmm.org/latest/api-python/generated/openmm.app.simulation.Simulation.html?#openmm.app.simulation.Simulation.minimizeEnergy) using a [MinimizationReporter](http://docs.openmm.org/development/api-python/latest/openmm.openmm.MinimizationReporter.html). Note that `MinimizationReporter` was introduced in OpenMM 8.1.\n", "\n", - "The syntax for doing this is somewhat different to typical OpenMM functionality. Read the [API documentation](http://docs.openmm.org/development/api-python/generated/openmm.openmm.MinimizationReporter.html) for more explanation.\n", + "The syntax for doing this is somewhat different to typical OpenMM functionality. You will need to define a subclass of `MinimizationReporter` that has a `report()` method. Within this report method you can write code to take any action you want. The report method is called every iteration of the minimizer. Read the [API documentation](http://docs.openmm.org/latest/api-python/generated/openmm.openmm.MinimizationReporter.html) for more explanation.\n", "\n", - "First we will create a test system." + "First we will create a test system then we will show an example." ] }, { @@ -40,7 +40,7 @@ "id": "2e7b659d", "metadata": {}, "source": [ - "We then create a `MinimizationReporter` class that will print the current energy to the screen every 100 steps and save the energies to an array which we can plot later" + "Below is an example which prints the current energy to the screen and saves the energy to an array for plotting. The comments explain each part." ] }, { @@ -50,17 +50,56 @@ "metadata": {}, "outputs": [], "source": [ - "# Define a subclass of MinimizationReporter\n", - "class Reporter(MinimizationReporter):\n", - " interval = 100 # report interval\n", + "\n", + "# The class can have any name but it must subclass MinimizationReporter.\n", + "class MyMinimizationReporter(MinimizationReporter):\n", + "\n", + " # within the class you can declare variables that persist throughout the\n", + " # minimization\n", + "\n", " energies = [] # array to record progress\n", + "\n", + " # you must override the report method and it must have this signature.\n", " def report(self, iteration, x, grad, args):\n", - " # print current system energy to screen \n", - " if iteration % self.interval == 0:\n", - " print(iteration, args['system energy'])\n", + " '''\n", + " the report method is called every iteration of the minimization.\n", + " \n", + " Args:\n", + " iteration (int): The index of the current iteration. This refers \n", + " to the current call to the L-BFGS optimizer.\n", + " Each time the minimizer increases the restraint strength, \n", + " the iteration index is reset to 0.\n", + "\n", + " x (array-like): The current particle positions in flattened order: \n", + " the three coordinates of the first particle, \n", + " then the three coordinates of the second particle, etc.\n", + "\n", + " grad (array-like): The current gradient of the objective function \n", + " (potential energy plus restraint energy) with \n", + " respect to the particle coordinates, in flattened order.\n", + "\n", + " args (dict): Additional statistics described above about the current state of minimization. \n", + " In particular:\n", + " “system energy”: the current potential energy of the system\n", + " “restraint energy”: the energy of the harmonic restraints\n", + " “restraint strength”: the force constant of the restraints (in kJ/mol/nm^2)\n", + " “max constraint error”: the maximum relative error in the length of any constraint\n", + "\n", + " Returns:\n", + " bool : Specify if minimization should be stopped.\n", + " '''\n", + "\n", + " # Within the report method you write the code you want to be executed at \n", + " # each iteration of the minimization.\n", + " # In this example we get the current energy, print it to the screen, and save it to an array. \n", "\n", - " # save energy at each iteration to an array we can use later\n", - " self.energies.append(args['system energy'])\n", + " current_energy = args['system energy']\n", + "\n", + " if iteration % 100 == 0: # only print to screen every 100 iterations for clarity of webpage display\n", + " print(current_energy)\n", + "\n", + " \n", + " self.energies.append(current_energy)\n", "\n", " # The report method must return a bool specifying if minimization should be stopped. \n", " # You can use this functionality for early termination.\n", @@ -82,17 +121,29 @@ "metadata": {}, "outputs": [], "source": [ - "# Create an instance of our reporter\n", - "reporter = Reporter()\n", + "reporter = MyMinimizationReporter()\n", "\n", - "# Perform local energy minimization \n", "simulation.minimizeEnergy(reporter=reporter)\n", "\n", - "# plot the minimization progress\n", "import matplotlib.pyplot as plt\n", - "plt.plot(Reporter.energies)\n", + "plt.plot(reporter.energies)\n", + "plt.ylabel(\"System energy (kJ/mol)\")\n", + "plt.xlabel(\"Minimization iteration\")\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "id": "ccc8979c", + "metadata": {}, + "source": [ + "You will notice that the energy does not change continuously and does not always decrease. This is because the L-BFGS algorithm used by the minimizer does not support constraints. The minimizer therefore replaces all constraints with harmonic restraints, then performs unconstrained minimization of a combined objective function that is the sum of the system’s potential energy and the restraint energy. Once minimization completes, it checks whether all constraints are satisfied to an acceptable tolerance. It not, it increases the strength of the harmonic restraints and performs additional minimization. If the error in constrained distances is especially large, it may choose to throw out all work that has been done so far and start over with stronger restraints. This has several important consequences:\n", + "\n", + "- The objective function being minimized not actually the same as the potential energy. \n", + "- The objective function and the potential energy can both increase between iterations. \n", + "- The total number of iterations performed could be larger than the number specified by the maxIterations argument, if that many iterations leaves unacceptable constraint errors. \n", + "- All work is provisional. It is possible for the minimizer to throw it out and start over. " + ] } ], "metadata": { @@ -111,7 +162,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.0" + "version": "3.12.1" }, "tags": [ "barostat", From 08764b3ba8809e89c8d28aab70c34a45d5da8ec2 Mon Sep 17 00:00:00 2001 From: "S. Farr" Date: Tue, 2 Jan 2024 17:05:28 +0100 Subject: [PATCH 4/4] fix link --- notebooks/cookbook/report_minimization.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/cookbook/report_minimization.ipynb b/notebooks/cookbook/report_minimization.ipynb index e8d3b77..9b155eb 100644 --- a/notebooks/cookbook/report_minimization.ipynb +++ b/notebooks/cookbook/report_minimization.ipynb @@ -7,7 +7,7 @@ "source": [ "## Reporting Minimization\n", "\n", - "You can report the status of [simulation.minimizeEnergy()](http://docs.openmm.org/latest/api-python/generated/openmm.app.simulation.Simulation.html?#openmm.app.simulation.Simulation.minimizeEnergy) using a [MinimizationReporter](http://docs.openmm.org/development/api-python/latest/openmm.openmm.MinimizationReporter.html). Note that `MinimizationReporter` was introduced in OpenMM 8.1.\n", + "You can report the status of [simulation.minimizeEnergy()](http://docs.openmm.org/latest/api-python/generated/openmm.app.simulation.Simulation.html?#openmm.app.simulation.Simulation.minimizeEnergy) using a [MinimizationReporter](http://docs.openmm.org/latest/api-python/generated/openmm.openmm.MinimizationReporter.html). Note that `MinimizationReporter` was introduced in OpenMM 8.1.\n", "\n", "The syntax for doing this is somewhat different to typical OpenMM functionality. You will need to define a subclass of `MinimizationReporter` that has a `report()` method. Within this report method you can write code to take any action you want. The report method is called every iteration of the minimizer. Read the [API documentation](http://docs.openmm.org/latest/api-python/generated/openmm.openmm.MinimizationReporter.html) for more explanation.\n", "\n",