Fetch and simulate PDB in openmm.

do_everything.py

from sys import stdout
import simtk.openmm as mm
import simtk.openmm.app as app
import simtk.unit as u
import urllib

def write_pdb(pdb_id):
  url = 'http://www.rcsb.org/pdb/files/%s.pdb' % pdb_id
  pdb_text = urllib.urlopen(url).read()
  file_handle = open("%s.pdb" % pdb_id, 'w')
  file_handle.writelines(pdb_text)
  file_handle.close()

pdb_id = "1vii"
pdb_filename = "%s.pdb" % pdb_id
write_pdb(pdb_id)

pdb = app.PDBFile(pdb_filename)


temperature = 300 * u.kelvin
ionic_strength = 0.01 * u.molar

forcefield = app.ForceField('amber99sbnmr.xml', 'tip3p.xml')

modeller = app.modeller.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
modeller.addSolvent(forcefield, padding=0.8 * u.nanometers, ionicStrength=ionic_strength)


system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME, 
    nonbondedCutoff=0.95*u.nanometers, constraints=app.HAngles, rigidWater=True, 
    ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(temperature, 1.0 / u.picoseconds, 2.0*u.femtoseconds)
integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'single'}
simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)

print('Minimizing...')
simulation.minimizeEnergy()

simulation.context.setVelocitiesToTemperature(temperature)
print('Equilibrating...')
simulation.step(10000)

simulation.reporters.append(app.DCDReporter('output.dcd', 10000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, 
    potentialEnergy=True, temperature=True))

print('Running Production...')
simulation.step(10000)
print('Done!')