openmm.py

##########################################################################
# this script was generated by openmm-builder. to customize it further,
# you can save the file to disk and edit it with your favorite editor.
##########################################################################

from __future__ import print_function
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import numpy as np

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')

system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
    nonbondedCutoff=1.0*nanometers, constraints=HAngles, rigidWater=True)
integrator = LangevinIntegrator(300*kelvin, 91.0/picoseconds, 2.0*femtoseconds)

platform = Platform.getPlatformByName('CUDA')
properties = {'CudaDeviceIndex': '0', 'CudaPrecision': 'mixed'}
simulation = Simulation(pdb.topology, system, integrator, platform, properties)
simulation.context.setPositions(pdb.positions)

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

# Generate random initial velocities from Maxwell-Boltzmann distribution.
velocities = Quantity(np.zeros([system.getNumParticles(), 3], np.float32),
                      nanometers / picosecond)
kT = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA * 300*kelvin
for atom_index in range(natoms):
    atom_mass = system.getParticleMass(atom_index)
    # standard deviation of velocity distribution for each coord for this atom
    sigma = sqrt(kT / atom_mass)
    velocities[atom_index, :] = sigma * np.random.normal(size=3)
system.context.setVelocities(velocities)

print('Equilibrating...')
simulation.step(100)

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

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