Reporting Minimization¶
You can report the status of simulation.minimizeEnergy() using a MinimizationReporter. Note that MinimizationReporter was introduced in OpenMM 8.1.
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 for more explanation.
First we will create a test system then we will show an example.
[1]:
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
pdb = PDBFile('villin.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
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.
[2]:
# The class can have any name but it must subclass MinimizationReporter.
class MyMinimizationReporter(MinimizationReporter):
# within the class you can declare variables that persist throughout the
# minimization
energies = [] # array to record progress
# you must override the report method and it must have this signature.
def report(self, iteration, x, grad, args):
'''
the report method is called every iteration of the minimization.
Args:
iteration (int): The index of the current iteration. This refers
to the current call to the L-BFGS optimizer.
Each time the minimizer increases the restraint strength,
the iteration index is reset to 0.
x (array-like): The current particle positions in flattened order:
the three coordinates of the first particle,
then the three coordinates of the second particle, etc.
grad (array-like): The current gradient of the objective function
(potential energy plus restraint energy) with
respect to the particle coordinates, in flattened order.
args (dict): Additional statistics described above about the current state of minimization.
In particular:
“system energy”: the current potential energy of the system
“restraint energy”: the energy of the harmonic restraints
“restraint strength”: the force constant of the restraints (in kJ/mol/nm^2)
“max constraint error”: the maximum relative error in the length of any constraint
Returns:
bool : Specify if minimization should be stopped.
'''
# Within the report method you write the code you want to be executed at
# each iteration of the minimization.
# In this example we get the current energy, print it to the screen, and save it to an array.
current_energy = args['system energy']
if iteration % 100 == 0: # only print to screen every 100 iterations for clarity of webpage display
print(current_energy)
self.energies.append(current_energy)
# The report method must return a bool specifying if minimization should be stopped.
# You can use this functionality for early termination.
return False
We now create an instance of the reporter and minimize the system with the reporter attached.
[3]:
reporter = MyMinimizationReporter()
simulation.minimizeEnergy(reporter=reporter)
import matplotlib.pyplot as plt
plt.plot(reporter.energies)
plt.ylabel("System energy (kJ/mol)")
plt.xlabel("Minimization iteration")
plt.show()
-142752.4829820002
-172600.61286814426
-176350.9121980728
-178082.39053243893
-179001.32075730423
-179449.9535824655
-179787.81424348694
-130436.52349621359
-153731.81033183183
-158893.09683325
-161371.83576248403
-162798.59830031695
-163750.27338761103
-164583.5929359414
-165111.86466812095
-165685.33011257192
-166095.94886591853
-166419.86401987408
-166714.80609204687
-166905.94190024622
-167103.17940697435
-167268.17701494784
-167420.67852697457
-167555.70245464446
-167657.1597443788
-167756.66600105583
-167831.22170631718
-167907.68236496928
-167963.4909734041
-167999.72042089768
-168059.21264844976
-168105.44811788748
-167198.26465962586
-166862.8159779532
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:
The objective function being minimized not actually the same as the potential energy.
The objective function and the potential energy can both increase between iterations.
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.
All work is provisional. It is possible for the minimizer to throw it out and start over.
