Download Notebook View in GitHub Open in Google Colab

Computing Interaction Energies

Sometimes you want to compute the interaction energy between two molecules, or between two sets of molecules such as solvent and solute. This is complicated by the fact that many Force classes can compute nonbonded interactions, and different ones must be handled differently. We consider three cases.

1. NonbondedForce is the most common class used for nonbonded interactions. It does not have an option to directly calculate interaction energies, only the total energy of the whole system. Instead we can perform three energy evaluations: one for the two molecules together, and one for each of the molecules individually to get its internal energy. Subtracting gives the interaction energy.

2. CustomNonbondedForce is also often used to compute nonbonded interactions. It supports “interaction groups”, which can be used to compute only the interaction energy between two groups of particles.

3. Some interactions are not pairwise, such as implicit solvent or polarizable force fields. The interaction between two particles depends on many other particles, including ones in other molecules. In these cases, the concept of an “interaction energy” is not well defined. You must consider carefully what quantities you actually want to calculate and why. We do not consider this case further here.

Let’s start by loading a file for the villin headpiece in water, and modelling it with the CHARMM36 force field.

from openmm import *
from openmm.app import *
from openmm.unit import *

pdb = PDBFile('villin.pdb')
forcefield = ForceField('charmm36.xml', 'charmm36/water.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME)

We will compute a solute-solvent interaction energy. We need to identify the two sets of atoms whose interaction we want. In this file, the solvent consists of water and chloride ions. We can select the atoms by residue name.

solvent = set([a.index for a in pdb.topology.atoms() if a.residue.name in ('HOH', 'Cl')])
protein = set([a.index for a in pdb.topology.atoms() if a.index not in solvent])

Now we will modify the Forces in the System as follows.

• For NonbondedForce objects, we will add parameter offsets that can be used to “zero out” the parameters of particles in each set, causing them to not interact. We also zero out the exceptions, since they are used for bonds within a single molecule, not for interactions between molecules.

• For CustomNonbondedForce objects, we add interaction groups to compute just the solute-solvent interaction energy.

• We also sort the Force objects into different force groups so we can evaluate them separately: group 0 for NonbondedForce, group 1 for CustomNonbondedForce, and group 2 for everything else.

for force in system.getForces():
    if isinstance(force, NonbondedForce):
        force.setForceGroup(0)
        force.addGlobalParameter("solute_scale", 1)
        force.addGlobalParameter("solvent_scale", 1)
        for i in range(force.getNumParticles()):
            charge, sigma, epsilon = force.getParticleParameters(i)
            # Set the parameters to be 0 when the corresponding parameter is 0,
            # and to have their normal values when it is 1.
            param = "solute_scale" if i in protein else "solvent_scale"
            force.setParticleParameters(i, 0, 0, 0)
            force.addParticleParameterOffset(param, i, charge, sigma, epsilon)
        for i in range(force.getNumExceptions()):
            p1, p2, chargeProd, sigma, epsilon = force.getExceptionParameters(i)
            force.setExceptionParameters(i, p1, p2, 0, 0, 0)
    elif isinstance(force, CustomNonbondedForce):
        force.setForceGroup(1)
        force.addInteractionGroup(protein, solvent)
    else:
        force.setForceGroup(2)

Let’s create a Context for performing calculations. The integrator is not important, since we will only be performing single point energy evaluations.

integrator = VerletIntegrator(0.001*picosecond)
context = Context(system, integrator)
context.setPositions(pdb.positions)

CHARMM36 uses NonbondedForce for Coulomb interactions and CustomNonbondedForce for Lennard-Jones interactions. To compute the Coulomb interaction energy, we evaluate group 0 three times to subtract the internal energy of each set from the total energy.

def coulomb_energy(solute_scale, solvent_scale):
    context.setParameter("solute_scale", solute_scale)
    context.setParameter("solvent_scale", solvent_scale)
    return context.getState(getEnergy=True, groups={0}).getPotentialEnergy()

total_coulomb = coulomb_energy(1, 1)
solute_coulomb = coulomb_energy(1, 0)
solvent_coulomb = coulomb_energy(0, 1)
print(total_coulomb - solute_coulomb - solvent_coulomb)

-6345.739608764648 kJ/mol

The Lennard-Jones interaction energy is much simpler. We just evaluate group 1.

print(context.getState(getEnergy=True, groups={1}).getPotentialEnergy())

-472.0870351791382 kJ/mol

Other force fields may divide up the energy differently. For example, they may use a single NonbondedForce to compute both Coulomb and Lennard-Jones interactions. In that case, the energy computed from NonbondedForce alone represents the total interaction energy.

If you still want to separate the Coulomb and Lennard-Jones interactions in that case, it can be done by defining separate parameters for the two. In this example we decompose the interaction energy for Amber14, which uses a single NonbondedForce for all nonbonded interactions.

forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME)
for force in system.getForces():
    if isinstance(force, NonbondedForce):
        force.setForceGroup(0)
        force.addGlobalParameter("solute_coulomb_scale", 1)
        force.addGlobalParameter("solute_lj_scale", 1)
        force.addGlobalParameter("solvent_coulomb_scale", 1)
        force.addGlobalParameter("solvent_lj_scale", 1)
        for i in range(force.getNumParticles()):
            charge, sigma, epsilon = force.getParticleParameters(i)
            force.setParticleParameters(i, 0, 0, 0)
            if i in protein:
                force.addParticleParameterOffset("solute_coulomb_scale", i, charge, 0, 0)
                force.addParticleParameterOffset("solute_lj_scale", i, 0, sigma, epsilon)
            else:
                force.addParticleParameterOffset("solvent_coulomb_scale", i, charge, 0, 0)
                force.addParticleParameterOffset("solvent_lj_scale", i, 0, sigma, epsilon)
        for i in range(force.getNumExceptions()):
            p1, p2, chargeProd, sigma, epsilon = force.getExceptionParameters(i)
            force.setExceptionParameters(i, p1, p2, 0, 0, 0)
    else:
        force.setForceGroup(2)

Now we can evaluate the interaction energies as before, by subtracting internal energies from the total energy.

integrator = VerletIntegrator(0.001*picosecond)
context = Context(system, integrator)
context.setPositions(pdb.positions)

def energy(solute_coulomb_scale, solute_lj_scale, solvent_coulomb_scale, solvent_lj_scale):
    context.setParameter("solute_coulomb_scale", solute_coulomb_scale)
    context.setParameter("solute_lj_scale", solute_lj_scale)
    context.setParameter("solvent_coulomb_scale", solvent_coulomb_scale)
    context.setParameter("solvent_lj_scale", solvent_lj_scale)
    return context.getState(getEnergy=True, groups={0}).getPotentialEnergy()

total_coulomb = energy(1, 0, 1, 0)
solute_coulomb = energy(1, 0, 0, 0)
solvent_coulomb = energy(0, 0, 1, 0)
total_lj = energy(0, 1, 0, 1)
solute_lj = energy(0, 1, 0, 0)
solvent_lj = energy(0, 0, 0, 1)
print('Coulomb interaction energy:', total_coulomb - solute_coulomb - solvent_coulomb)
print('LJ interaction energy:', total_lj - solute_lj - solvent_lj)

Coulomb interaction energy: -5638.60320864596 kJ/mol
LJ interaction energy: 220.66263641439218 kJ/mol