Constraining Atom Positions¶
Sometimes you want to run a simulation in which certain particles are kept completely fixed, unable to move at all. This is called constraining particles (not to be confused with restraining particles, which is described in another entry).
The easiest way to constrain particles is to set their masses to zero. When a particle has zero mass, it is treated fixed. This means
Integrators will never change their positions.
The local energy minimizer will never change their positions.
They do not contribute to the kinetic energy of 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.
[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 set the masses of some particles to 0. You can choose any set of particles you want, depending on your application. For example, suppose we want to equilibrate the water while keeping the protein fixed. Since the protein makes up the first 582 particles, we could write
[2]:
for i in range(582):
system.setParticleMass(i, 0*amu)
Alternatively we could do it based on information from the Topology, such as the names of atoms, residues, or chains. For example, this sets the mass of every particle that is not water to 0.
[3]:
for atom in pdb.topology.atoms():
if atom.residue.name != 'HOH':
system.setParticleMass(atom.index, 0*amu)
Now that we have our System ready, we can create a Simulation and run a few steps.
[4]:
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.step(10)
Let’s query the minimized positions and see how the positions of a few atoms have changed.
[5]:
state = simulation.context.getState(getPositions=True)
positions = state.getPositions()
print('These positions should be the same (at least up to the limits of single precision):')
print(pdb.positions[0])
print(positions[0])
print()
print('These positions should be different:')
print(pdb.positions[600])
print(positions[600])
These positions should be the same (at least up to the limits of single precision):
Vec3(x=2.516, y=1.4160000000000001, z=1.9440000000000002) nm
Vec3(x=2.5160000324249268, y=1.4160000085830688, z=1.944000005722046) nm
These positions should be different:
Vec3(x=1.7550000000000001, y=2.556, z=3.3930000000000002) nm
Vec3(x=1.7894033193588257, y=2.570143699645996, z=3.349152088165283) nm
