Restraining Dihedrals¶
Sometimes you want to run a simulation in which certain dihedral angles have their motion restricted, unable to move out of a narrow range. You might do this to prevent large changes to the backbone of a protein, or to sample a particular sidechain rotamer. This can be implemented using a PeriodicTorsionForce to add restraint potentials to the dihedrals of interest.
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.
[1]:
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 PeriodicTorsionForce to use as the restraint force. We can add as many torsions as we want to it, each one restraining a different dihedral. The energy of each torsion is given by \(E=k(1+cos(n\theta-\theta_0))\), where \(k\) (the force constant), \(n\) (the periodicity), and \(\theta_0\) (the phase) can be specified separately for each torsion. For this example, we restrain the first residue’s \(\chi_1\) angle, which is defined by atoms 0, 4, 6, and 9. We set the periodicity to 1 since we only want a single minimum in the potential.
[2]:
restraint = PeriodicTorsionForce()
system.addForce(restraint)
restraint.addTorsion(0, 4, 6, 9, 1, 0*radians, 100*kilojoules_per_mole)
[2]:
0
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 first sidechain.
[3]:
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 angle of the restraint part way through a simulation. This can happen, for example, in steered molecular dynamics where you want to make a dihedral 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 changes the target angle of the restraint to 0.1 radians.
[4]:
restraint.setTorsionParameters(0, 0, 4, 6, 9, 1, 0.1*radians, 100*kilojoules_per_mole)
restraint.updateParametersInContext(simulation.context)
