Example 1: Bovine Pancreatic Trypsin Inhibitor¶
This is a simple example showing how grand can be used to simulate GCMC/MD for a protein solvated in water. The majority of this script below is composed of standard OpenMM functions, with a few additional parts to execute grand canonical sampling. These additional functions are described here.
The
grand.utils.add_ghosts()function is used to add some additional waters to the system topology, so these can be used for GCMC insertion moves - until inserted, these waters are non-interacting.The
ref_atomsis used to choose the set of atoms, for which the mean coordinate is used to define the centre of the GCMC sphere - this should be chosen carefully (along with the sphere radius) to ensure the region of interest is well covered.The
grand.samplers.StandardGCMCSphereSamplerobject contains all of the necessary variables for carrying out GCMC moves, and the arguments given should be self-explanatory.The
gcmc_mover.initialise()function must be executed before starting the simulation, as this feeds some context-specific variables togcmc_moverand ensures that the ghosts are switched off.The
gcmc_mover.deleteWatersInGCMCSphere()removes any waters present in the GCMC sphere at the beginning of the simulation, so that the water sampling will be less biased by the initial water locations.The
gcmc_mover.move(simulation.context, 200)function executes a number of GCMC moves at a given point. For reasons of efficiency, it is best to carry these out in blocks of at least ~20 moves.By running
gcmc_mover.report(), a simulation frame is written out and the log file is updated.
The remaining functions are used to process the trajectory for visualisation and analysis:
grand.utils.shift_ghost_waters()translates the ghost waters far from the simulation box, so that they will not be confused with interacting waters.grand.utils.recentre_traj()is used to recentre the trajectory on a particular atom. However, this can be expensive, so if this atom does not get close to the edges of the simulation cell, this is not necessary.grand.utils.align_traj()is used to align the protein with respect to the initial frame (or a reference structure, via thereferenceargument).grand.utils.write_sphere_trajwrites out a PDB trajectory of the GCMC sphere, which may be helpful for visualisation.
The documentation for these functions can be found in the “grand package” section. The full script is included below.
"""
Description
-----------
Example script of how to run GCMC/MD in OpenMM for a BPTI system
Note that this simulation is only an example, and is not long enough
to see equilibrated behaviour
Marley Samways
"""
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
from openmmtools.integrators import BAOABIntegrator
import grand
# Load in PDB file
pdb = PDBFile('bpti-equil.pdb')
# Add ghost water molecules, which can be inserted
pdb.topology, pdb.positions, ghosts = grand.utils.add_ghosts(pdb.topology,
pdb.positions,
n=5,
pdb='bpti-ghosts.pdb')
# Create system
ff = ForceField('amber14-all.xml', 'amber14/tip3p.xml')
system = ff.createSystem(pdb.topology,
nonbondedMethod=PME,
nonbondedCutoff=12.0*angstroms,
switchDistance=10.0*angstroms,
constraints=HBonds)
# Define atoms around which the GCMC sphere is based
ref_atoms = [{'name': 'CA', 'resname': 'TYR', 'resid': '10'},
{'name': 'CA', 'resname': 'ASN', 'resid': '43'}]
gcmc_mover = grand.samplers.StandardGCMCSphereSampler(system=system,
topology=pdb.topology,
temperature=300*kelvin,
referenceAtoms=ref_atoms,
sphereRadius=4.2*angstroms,
log='bpti-gcmc.log',
dcd='bpti-raw.dcd',
rst='bpti-rst.rst7',
overwrite=False)
# BAOAB Langevin integrator
integrator = BAOABIntegrator(300*kelvin, 1.0/picosecond, 0.002*picoseconds)
platform = Platform.getPlatformByName('CUDA')
platform.setPropertyDefaultValue('Precision', 'mixed')
simulation = Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(300*kelvin)
simulation.context.setPeriodicBoxVectors(*pdb.topology.getPeriodicBoxVectors())
# Switch off ghost waters and those in sphere (to start fresh)
gcmc_mover.initialise(simulation.context, ghosts)
gcmc_mover.deleteWatersInGCMCSphere()
# Equilibrate water distribution - 10k moves over 5 ps
print("Equilibration...")
for i in range(50):
# Carry out 200 moves every 100 fs
gcmc_mover.move(simulation.context, 200)
simulation.step(50)
print("{}/{} equilibration GCMC moves accepted. N = {}".format(gcmc_mover.n_accepted,
gcmc_mover.n_moves,
gcmc_mover.N))
# Add StateDataReporter for production
simulation.reporters.append(StateDataReporter(stdout,
1000,
step=True,
potentialEnergy=True,
temperature=True,
volume=True))
# Reset GCMC statistics
gcmc_mover.reset()
# Run simulation - 5k moves over 50 ps
print("\nProduction")
for i in range(50):
# Carry out 100 GCMC moves per 1 ps of MD
simulation.step(500)
gcmc_mover.move(simulation.context, 100)
# Write data out
gcmc_mover.report(simulation)
#
# Need to process the trajectory for visualisation
#
# Shift ghost waters outside the simulation cell
trj = grand.utils.shift_ghost_waters(ghost_file='gcmc-ghost-wats.txt',
topology='bpti-ghosts.pdb',
trajectory='bpti-raw.dcd')
# Centre the trajectory on a particular residue
trj = grand.utils.recentre_traj(t=trj, resname='TYR', resid=10)
# Align the trajectory to the protein
grand.utils.align_traj(t=trj, output='bpti-gcmc.dcd')
# Write out a PDB trajectory of the GCMC sphere
grand.utils.write_sphere_traj(radius=4.2,
ref_atoms=ref_atoms,
topology='bpti-ghosts.pdb',
trajectory='bpti-gcmc.dcd',
output='gcmc_sphere.pdb',
initial_frame=True)