Analyzing Energy Contributions¶
It is sometimes useful to decompose the total energy of a system into separate components: bonds, angles, nonbonded interactions, etc. In OpenMM, a force field is made up of Force objects, one for each type of interaction. Decomposing the energy means identifying the contribution of each Force object.
Let’s start by loading a PDB file and building a System.
[1]:
from openmm.app import *
from openmm import *
from openmm.unit import *
pdb = PDBFile('ala_ala_ala.pdb')
forcefield = ForceField('amber14-all.xml')
system = forcefield.createSystem(pdb.topology)
OpenMM does not have a way to directly query the energy of a Force object. Instead, it lets you query the energy of a force group. We therefore need each Force object to be in a different group.
[2]:
for i, f in enumerate(system.getForces()):
f.setForceGroup(i)
Like any modification to a System or the Forces it contains, this must be done before you create a Simulation. Now we can go ahead and create it.
[3]:
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
We use Context.getState() to query any kind of state information: positions, velocities, forces, energy, etc. It has an optional groups argument which can take a set of force groups. When this is provided, the forces and energy returned in the State reflect only the specified groups. We can take advantage of this to get the energy from each group (which corresponds to a single Force, since we put each one into a different group).
[4]:
for i, f in enumerate(system.getForces()):
state = simulation.context.getState(getEnergy=True, groups={i})
print(f.getName(), state.getPotentialEnergy())
HarmonicBondForce 18.088199615478516 kJ/mol
HarmonicAngleForce 70.43384552001953 kJ/mol
PeriodicTorsionForce 84.32029724121094 kJ/mol
NonbondedForce -130.988525390625 kJ/mol
CMMotionRemover 0.0 kJ/mol
