script.py

# This script was generated by OpenMM-Setup on 2020-07-28.

from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *

#try:
#    platform = Platform.getPlatformByName('OpenCL')
#    platformProperties = {'Precision': 'mixed'}
#except:



platform = Platform.getPlatformByName('CPU')
platformProperties = {'Threads': 5}

# Input Files

pdb = PDBFile('1wa3_BA-05_modl-processed.pdb')
forcefield = ForceField('charmm36.xml', 'charmm36/water.xml')

# System Configuration

nonbondedMethod = PME
nonbondedCutoff = 1*nanometers
ewaldErrorTolerance = 1.0e-5
constraints = HBonds
rigidWater = True
constraintTolerance = 1.0e-5

# Integration Options

dt = 2*picoseconds
temperature = 303*kelvin
friction = 1/picosecond
pressure = 1*atmospheres
barostatInterval = 100

# Simulation Options

steps = 1000000
equilibrationSteps = 10000


dcdReporter = DCDReporter('out.dcd', 0.01*nanosecond)
dataReporter = StateDataReporter('out.log', 0.01*nanosecond, totalSteps=steps,
    step=True, time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True, volume=True, density=True, separator='\t')

# Prepare the Simulation

print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
#simulation = Simulation(topology, system, integrator, platform=platform, platformProperties=platformProperties)
simulation = Simulation(topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(positions)

# Minimize and Equilibrate

print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')

from timeit import default_timer as timer
start = timer()
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)
end = timer()


print(f'PLATFORM: {platform.getName()}')
ns_run = dt*equilibrationSteps/nanosecond
real_time_days = (end-start)/3600/24
ns_per_day = ns_run/real_time_days

print(f'Running {ns_per_day:2.2f} ns/day')

# Simulate

print('Simulating...')
#simulation.reporters.append(dcdReporter)
#simulation.reporters.append(dataReporter)
#simulation.currentStep = 0
#simulation.step(steps)