From 1f93f528efe694f31ce234098eae9440f22306d2 Mon Sep 17 00:00:00 2001 From: "S. Farr" Date: Thu, 28 Dec 2023 22:11:48 +0100 Subject: [PATCH] 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",