Download Notebook View in GitHub Open in Google Colab

Restraining Atom Positions

Sometimes you want to run a simulation in which certain particles have their motion restricted, unable to move too far from where they started. This is called restraining particles (not to be confused with contraining particles, in which they are kept completely fixed, and which is described in another entry).

Position restraints are typically implemented by adding a harmonic force that binds each particle to its initial position. The further it moves, the stronger the force pushing it back. This is generally done by adding a CustomExternalForce to the System.

As an example, let’s copy the beginning of the simulatePdb.py script from the examples directory. It loads a PDB file consisting of villin in water and builds a System from it.

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

pdb = PDBFile('villin.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
                                 nonbondedCutoff=1*nanometer, constraints=HBonds)

Now we will create a CustomExternalForce that binds each particle to a specified location.

restraint = CustomExternalForce('k*periodicdistance(x, y, z, x0, y0, z0)^2')
system.addForce(restraint)
restraint.addGlobalParameter('k', 100.0*kilojoules_per_mole/nanometer)
restraint.addPerParticleParameter('x0')
restraint.addPerParticleParameter('y0')
restraint.addPerParticleParameter('z0')

2

The energy of each particle equals a global constant k multiplied by the square of the distance between the particle’s current position (x, y, z) and a reference position (x0, y0, z0). We specify that x0, y0, and z0 are per-particle parameters, so each particle can have a different reference position. Also note that we compute the distance with the periodicdistance() function. This takes periodic boundary conditions into account. If we were instead simulating a non-periodic system, we would specify the energy expression as 'k*((x-x0)^2+(y-y0)^2+(z-z0)^2)'.

Now we call addParticle() to tell it which particles to apply the restraining force to. You can choose the particles in any way you want, depending on your application. For example, suppose we want to restrain the alpha carbons so the side chains and solvent can equilibrate without large changes to the overall fold. We could write

for atom in pdb.topology.atoms():
    if atom.name == 'CA':
        restraint.addParticle(atom.index, pdb.positions[atom.index])

Now that we have our System ready, we can create a Simulation and run some dynamics. The restraint force will limit the motion of the backbone.

integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.step(10)

Sometimes it is useful to change the target position of the restraint part way through a simulation. This can happen, for example, in steered molecular dynamics where you want to make a particle follow a defined trajectory with time. To do this, simply update the values of the per-particle parameters then call updateParametersInContext() to copy the new values over to the existing Context. The following moves the target position of the first restrained particle by 0.1 nm in the x direction.

index, position = restraint.getParticleParameters(0)
restraint.setParticleParameters(0, index, position+Vec3(0.1, 0, 0))
restraint.updateParametersInContext(simulation.context)