example_tip5p_nvt_nve.py

from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

pdb = app.PDBFile('water256_bulk_tip5p.pdb')
#pdb = app.PDBFile('a.pdb')
#pdb = app.PDBFile('water14_cluster.pdb')
boxDim = 19.3996888399961804/10.
temperature = 300
boxSize = (boxDim, boxDim, boxDim) * unit.nanometer
pdb.topology.setUnitCellDimensions(boxSize)
forcefield = app.ForceField('amber10.xml', 'tip5p.xml')

################################## NVT

system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=0.9*unit.nanometers, constraints=None, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.VerletIntegrator(0.2*unit.femtoseconds)

system.addForce(mm.AndersenThermostat(temperature, 1./unit.picoseconds))

# platform = mm.Platform.getPlatformByName('Reference')
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'double'}
simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.computeVirtualSites()


simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)

simulation.reporters.append(app.StateDataReporter("tip5p_nvt.log", 50, step=True, time=True,
    potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True,
    progress=True, remainingTime=True, speed=True, totalSteps=500000,
    separator='\t'))

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

final_nvt_state = simulation.context.getState(getVelocities=True, getPositions=True)
positions = final_nvt_state.getPositions()
velocities = final_nvt_state.getVelocities()
################################## NVE

system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=0.9*unit.nanometers, constraints=None, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.VerletIntegrator(0.2*unit.femtoseconds)

simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(positions)
simulation.context.computeVirtualSites()
simulation.context.setVelocities(velocities)

simulation.reporters.append(app.StateDataReporter("tip5p_nve.log", 50, step=True, time=True,
    potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True,
    progress=True, remainingTime=True, speed=True, totalSteps=500000,
    separator='\t'))

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