TJ's script for running Osama's j-coupling simulations.

run_jcpl_sim.py

#!/usr/bin/env python

"""
Usage: ./run_sim.py <pdb> <cuda_device> <output_name>
"""

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


print(sys.argv)
print("Running %s on device %s" % tuple(sys.argv[1:3]))
input_file  = sys.argv[1]
output_file = sys.argv[3]

# ---------------------------------------------------------------------------- #
# Parameters

pdb = app.PDBFile(input_file)
timestep = 2.5 * unit.femtoseconds
sim_len  = 100 * unit.nanoseconds

ionic_strengths = { '1BU5' : 60.00,  # mM NaCl
                    '1IGD' : 308.0,  # mM NaCl
                    '1UBQ' : 30.00 } # mM NaCl
                    
pHs =             { '1BU5' : 7.0,  
                    '1IGD' : 6.5,  
                    '1UBQ' : 4.7 }
                    
forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml')

# ---------------------------------------------------------------------------- #



if input_file.find('1BU5') != -1:
    system_id = '1BU5'
elif input_file.find('1IGD') != -1:
    system_id = '1IGD'
elif input_file.find('1UBQ') != -1:
    system_id = '1UBQ'
else:
    raise Exception('TJ, use this script only for Osamas project!')

# --- build the simulation box

modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield, pH = pHs[system_id])
modeller.addSolvent(forcefield, padding=0.75*unit.nanometers)
modeller.addSolvent(forcefield, ionicStrength = ionic_strengths[system_id] * unit.micromolar)


# --- build the system

number_of_steps = int( round(sim_len / timestep) )

system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME, 
    nonbondedCutoff=0.7*unit.nanometers, vdwCutoff=0.9*unit.nanometers,
    constraints=app.HBonds, rigidWater=True, ewaldErrorTolerance=0.0005)
    
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds, timestep)
integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'} #, 'CudaDeviceIndex': sys.argv[2]}
simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)

# save the built system to a PDB for later
positions = simulation.context.getState(getPositions=True).getPositions()
app.PDBFile.writeFile(simulation.topology, positions, open(system_id + '_wbox.pdb', 'w'))

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

simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating for 1ns...')
simulation.step( int( round(1.0 * unit.nanoseconds / timestep) ) )

interval = int( round(25.0 * unit.picoseconds / timestep) )

print('Saving output to: %s' % output_file)
simulation.reporters.append(app.DCDReporter(output_file, interval))
simulation.reporters.append(app.StateDataReporter(stdout, interval, step=True, 
                            potentialEnergy=True, temperature=True))

print('Running production for %d steps...' % number_of_steps)
simulation.step(number_of_steps)

print('Finished simulation!')