Okay. My apologies. The custom forces are loaded from an additional file along with some of the other settings. I will include them below. Almost all of these were obtained through CHARMM-GUI.
inputs are given in a separate file:
Code: Select all
nstep = 500000 # Number of steps to run
dt = 0.002 # Time-step (ps)
nstout = 1000 # Writing output frequency (steps)
nstdcd = 10000 # Writing coordinates trajectory frequency (steps)
coulomb = PME # Electrostatic cut-off method
ewaldTol = 0.0005 # Ewald error tolerance
vdw = Force-switch # vdW cut-off method
r_on = 1.0 # Switch-on distance (nm)
r_off = 1.2 # Switch-off distance (nm)
temp = 310.15 # Temperature (K)
fric_coeff = 1 # Friction coefficient for Langevin dynamics
pcouple = yes # Turn on/off pressure coupling
p_ref = 1.0 # Pressure (Pref or Pxx, Pyy, Pzz; bar)
p_type = membrane # MonteCarloBarostat type
p_XYMode = XYIsotropic # For MonteCarloMembraneBarostat
p_ZMode = ZFree # For MonteCarloMembraneBarostat
p_tens = 0.0 # Sulface tension for MonteCarloMembraneBarostat (dyne/cm)
p_freq = 100 # Pressure coupling frequency (steps)
cons = HBonds # Constraints mehtod
rest = yes # Turn on/off restraints
fc_mpos = 0
fc_x = 1
divide = 1
Then read by another file:
omm_readinputs
Code: Select all
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
class _OpenMMReadInputs():
def __init__(self):
self.mini_nstep = 0 # Number of steps for minimization
self.mini_Tol = 1.0 # Minimization energy tolerance
self.gen_vel = 'no' # Generate initial velocities
self.gen_temp = 300.0 # Temperature for generating initial velocities (K)
self.gen_seed = None # Seed for generating initial velocities
self.nstep = 0 # Number of steps to run
self.dt = 0.002 # Time-step (ps)
self.nstout = 100 # Writing output frequency (steps)
self.nstdcd = 0 # Wrtiing coordinates trajectory frequency (steps)
self.coulomb = PME # Electrostatic cut-off method
self.ewald_Tol = 0.0005 # Ewald error tolerance
self.vdw = 'Force-switch' # vdW cut-off method
self.r_on = 1.0 # Switch-on distance (nm)
self.r_off = 1.2 # Switch-off distance (nm)
self.temp = 300.0 # Temperature (K)
self.fric_coeff = 1 # Friction coefficient for Langevin dynamics
self.pcouple = 'no' # Turn on/off pressure coupling
self.p_ref = 1.0 # Pressure (Pref or Pxx, Pyy, Pzz; bar)
self.p_type = 'membrane' # MonteCarloBarotat type
self.p_scale = True, True, True # For MonteCarloAnisotropicBarostat
self.p_XYMode = MonteCarloMembraneBarostat.XYIsotropic # For MonteCarloMembraneBarostat
self.p_ZMode = MonteCarloMembraneBarostat.ZFree # For MonteCarloMembraneBarostat
self.p_tens = 0.0 # Sulface tension for MonteCarloMembraneBarostat (dyne/cm)
self.p_freq = 15 # Pressure coupling frequency (steps)
self.cons = HBonds # Constraints method
self.rest = 'no' # Turn on/off restraints
self.fc_x = 0.0 # shear force in x direction (kJ/mol/nm^2)
self.fc_bb = 0.0 # Positional restraint force constant for protein backbone (kJ/mol/nm^2)
self.fc_sc = 0.0 # Positional restraint force constant for protein side-chain (kJ/mol/nm^2)
self.fc_mpos = 0.0 # Positional restraint force constant for micelle lipids (kJ/mol/nm^2)
self.fc_lpos = 0.0 # Positional restraint force constant for lipids (kJ/mol/nm^2)
self.fc_ldih = 0.0 # Dihedral restraint force constant for lipids (kJ/mol/rad^2)
self.fc_cdih = 0.0 # Dihedral restraint force constant for carbohydrates (kJ/mol/rad^2)
self.divide = 0.0 # restraint force for the middle of the box.
def read(self, inputFile):
for line in open(inputFile, 'r'):
if line.find('#') >= 0: line = line.split('#')[0]
line = line.strip()
if len(line) > 0:
segments = line.split('=')
input_param = segments[0].upper().strip()
try: input_value = segments[1].strip()
except: input_value = None
if input_value:
if input_param == 'MINI_NSTEP': self.mini_nstep = int(input_value)
if input_param == 'MINI_TOL': self.mini_Tol = float(input_value)
if input_param == 'GEN_VEL':
if input_value.upper() == 'YES': self.gen_vel = 'yes'
if input_value.upper() == 'NO': self.gen_vel = 'no'
if input_param == 'GEN_TEMP': self.gen_temp = float(input_value)
if input_param == 'GEN_SEED': self.gen_seed = int(input_value)
if input_param == 'NSTEP': self.nstep = int(input_value)
if input_param == 'DT': self.dt = float(input_value)
if input_param == 'NSTOUT': self.nstout = int(input_value)
if input_param == 'NSTDCD': self.nstdcd = int(input_value)
if input_param == 'COULOMB':
if input_value.upper() == 'NOCUTOFF': self.coulomb = NoCutoff
if input_value.upper() == 'CUTOFFNONPERIODIC': self.coulomb = CutoffNonPeriodic
if input_value.upper() == 'CUTOFFPERIODIC': self.coulomb = CutoffPeriodic
if input_value.upper() == 'EWALD': self.coulomb = Ewald
if input_value.upper() == 'PME': self.coulomb = PME
if input_param == 'EWALD_TOL': self.ewald_Tol = float(input_value)
if input_param == 'VDW':
if input_value.upper() == 'FORCE-SWITCH': self.vdw = 'Force-switch'
if input_param == 'R_ON': self.r_on = float(input_value)
if input_param == 'R_OFF': self.r_off = float(input_value)
if input_param == 'TEMP': self.temp = float(input_value)
if input_param == 'FRIC_COEFF': self.fric_coeff = float(input_value)
if input_param == 'PCOUPLE':
if input_value.upper() == 'YES': self.pcouple = 'yes'
if input_value.upper() == 'NO': self.pcouple = 'no'
if input_param == 'P_REF':
if input_value.find(',') < 0:
self.p_ref = float(input_value)
else:
Pxx = float(input_value.split(',')[0])
Pyy = float(input_value.split(',')[1])
Pzz = float(input_value.split(',')[2])
self.p_ref = Pxx, Pyy, Pzz
if input_param == 'P_TYPE':
if input_value.upper() == 'ISOTROPIC': self.p_type = 'isotropic'
if input_value.upper() == 'MEMBRANE': self.p_type = 'membrane'
if input_value.upper() == 'ANISOTROPIC': self.p_type = 'anisotropic'
if input_param == 'P_SCALE':
scaleX = True
scaleY = True
scaleZ = True
if input_value.upper().find('X') < 0: scaleX = False
if input_value.upper().find('Y') < 0: scaleY = False
if input_value.upper().find('Z') < 0: scaleZ = False
self.p_scale = scaleX, scaleY, scaleZ
if input_param == 'P_XYMODE':
if input_value.upper() == 'XYISOTROPIC': self.p_XYMode = MonteCarloMembraneBarostat.XYIsotropic
if input_value.upper() == 'XYANISOTROPIC': self.p_XYMode = MonteCarloMembraneBarostat.XYAnisotropic
if input_param == 'P_ZMODE':
if input_value.upper() == 'ZFREE': self.p_ZMode = MonteCarloMembraneBarostat.ZFree
if input_value.upper() == 'ZFIXED': self.p_ZMode = MonteCarloMembraneBarostat.ZFixed
if input_value.upper() == 'CONSTANTVOLUME': self.p_ZMode = MonteCarloMembraneBarostat.ConstantVolume
if input_param == 'P_TENS': self.p_tens = float(input_value)
if input_param == 'P_FREQ': self.p_freq = int(input_value)
if input_param == 'CONS':
if input_value.upper() == 'NONE': self.cons = None
if input_value.upper() == 'HBONDS': self.cons = HBonds
if input_value.upper() == 'ALLBONDS': self.cons = AllBonds
if input_value.upper() == 'HANGLES': self.cons = HAngles
if input_param == 'REST':
if input_value.upper() == 'YES': self.rest = 'yes'
if input_value.upper() == 'NO': self.rest = 'no'
if input_param == 'FC_X': self.fc_x = float(input_value)
if input_param == 'FC_BB': self.fc_bb = float(input_value)
if input_param == 'FC_BB': self.fc_bb = float(input_value)
if input_param == 'FC_SC': self.fc_sc = float(input_value)
if input_param == 'FC_MPOS': self.fc_mpos = float(input_value)
if input_param == 'FC_LPOS': self.fc_lpos = float(input_value)
if input_param == 'FC_LDIH': self.fc_ldih = float(input_value)
if input_param == 'FC_CDIH': self.fc_cdih = float(input_value)
if input_param == 'DIVIDE': self.divide = float(input_value)
return self
def read_inputs(inputFile):
return _OpenMMReadInputs().read(inputFile)
We read the parameters of the system:
omm_readparams
Code: Select all
import os
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
def read_psf(filename):
psf = CharmmPsfFile(filename)
return psf
def read_crd(filename):
crd = CharmmCrdFile(filename)
return crd
def read_charmm_rst(filename):
charmm_rst = CharmmRstFile(filename)
for i, line in enumerate(charmm_rst.header):
line = line.strip()
words = line.split()
if len(line) != 0:
if words[0] == 'CRYSTAL' or words[0] == '!CRYSTAL':
line1 = charmm_rst.header[i+1]
line2 = charmm_rst.header[i+2]
boxlx = Vec3(float(line1.split()[0].replace("D", "E")), 0.0, 0.0)
boxly = Vec3(0.0, float(line1.split()[2].replace("D", "E")), 0.0)
boxlz = Vec3(0.0, 0.0, float(line2.split()[2].replace("D", "E")))
box = (boxlx, boxly, boxlz)
break
positions = charmm_rst.positionsold
new_positions = []
for position in positions:
oldx = position[0]/angstrom
oldy = position[1]/angstrom
oldz = position[2]/angstrom
newx = oldx + boxlx[0]/2.0
newy = oldy + boxly[1]/2.0
newz = oldz + boxlz[2]/2.0
new_positions.append(Vec3(newx, newy, newz))
charmm_rst.box = Quantity(box, angstroms)
charmm_rst.positions = Quantity(new_positions, angstroms)
return charmm_rst
def read_params(filename):
extlist = ['rtf', 'prm', 'str']
parFiles = ()
for line in open(filename, 'r'):
if '!' in line: line = line.split('!')[0]
parfile = line.strip()
if len(parfile) != 0:
ext = parfile.lower().split('.')[-1]
if not ext in extlist: continue
parFiles += ( parfile, )
params = CharmmParameterSet( *parFiles )
return params
def read_box(psf, filename):
for line in open(filename, 'r'):
segments = line.split('=')
if segments[0].strip() == "BOXLX": boxlx = float(segments[1])
if segments[0].strip() == "BOXLY": boxly = float(segments[1])
if segments[0].strip() == "BOXLZ": boxlz = float(segments[1])
psf.setBox(boxlx*angstroms, boxly*angstroms, boxlz*angstroms)
return psf
def gen_box(psf, crd):
coords = crd.positions
min_crds = [coords[0][0], coords[0][1], coords[0][2]]
max_crds = [coords[0][0], coords[0][1], coords[0][2]]
for coord in coords:
min_crds[0] = min(min_crds[0], coord[0])
min_crds[1] = min(min_crds[1], coord[1])
min_crds[2] = min(min_crds[2], coord[2])
max_crds[0] = max(max_crds[0], coord[0])
max_crds[1] = max(max_crds[1], coord[1])
max_crds[2] = max(max_crds[2], coord[2])
boxlx = max_crds[0]-min_crds[0]
boxly = max_crds[1]-min_crds[1]
boxlz = max_crds[2]-min_crds[2]
psf.setBox(boxlx, boxly, boxlz)
return psf
The file that handles PBC:
omm_rewrap
Code: Select all
from __future__ import print_function
import os
from math import *
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
def rewrap(simulation):
bonds = simulation.topology.bonds()
positions = simulation.context.getState(getPositions=True).getPositions()
box = simulation.context.getState().getPeriodicBoxVectors()
boxlx = box[0][0]/angstrom
boxly = box[1][1]/angstrom
boxlz = box[2][2]/angstrom
min_crds = [positions[0][0]/angstrom, positions[0][1]/angstrom, positions[0][2]/angstrom]
max_crds = [positions[0][0]/angstrom, positions[0][1]/angstrom, positions[0][2]/angstrom]
for position in positions:
min_crds[0] = min(min_crds[0], position[0]/angstrom)
min_crds[1] = min(min_crds[1], position[1]/angstrom)
min_crds[2] = min(min_crds[2], position[2]/angstrom)
max_crds[0] = max(max_crds[0], position[0]/angstrom)
max_crds[1] = max(max_crds[1], position[1]/angstrom)
max_crds[2] = max(max_crds[2], position[2]/angstrom)
xcen = (max_crds[0] + min_crds[0]) / 2.0
ycen = (max_crds[1] + min_crds[1]) / 2.0
zcen = (max_crds[2] + min_crds[2]) / 2.0
for bond in bonds:
atom1 = bond[0]
atom2 = bond[1]
atom1id = atom1.index
atom2id = atom2.index
res1 = atom1.residue
res2 = atom2.residue
x1, y1, z1 = positions[atom1id]
x2, y2, z2 = positions[atom2id]
dx = fabs(x1/angstrom - x2/angstrom)
dy = fabs(y1/angstrom - y2/angstrom)
dz = fabs(z1/angstrom - z2/angstrom)
if dx > boxlx/2 or dy > boxly/2 or dz > boxlz/2:
for atom in res2.atoms():
oldx = positions[atom.index][0]/angstrom
oldy = positions[atom.index][1]/angstrom
oldz = positions[atom.index][2]/angstrom
if dx > boxlx/2.0:
if oldx < xcen: newx = oldx + boxlx
else: newx = oldx - boxlx
else:
newx = oldx
if dy > boxly/2.0:
if oldy < ycen: newy = oldy + boxly
else: newy = oldy - boxly
else:
newy = oldy
if dz > boxlz/2.0:
if oldz < zcen: newz = oldz + boxlz
else: newz = oldz - boxlz
else:
newz = oldz
new_position = Vec3(newx, newy, newz)
positions[atom.index] = Quantity(new_position, angstroms)
simulation.context.setPositions(positions)
return simulation
omm_restraints:
Code: Select all
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
def restraints(system, psf, crd, inputs):
boxlx = system.getDefaultPeriodicBoxVectors()[0][0].value_in_unit(nanometers)
boxly = system.getDefaultPeriodicBoxVectors()[1][1].value_in_unit(nanometers)
boxlz = system.getDefaultPeriodicBoxVectors()[2][2].value_in_unit(nanometers)
if inputs.fc_x > 0:
# position restraints for protein
shearForce = CustomExternalForce('fc_x*x')
shearForce.addPerParticleParameter('fc_x')
for residue in psf.residue_list:
if residue.resname == 'TIP3':
for atom in residue.atoms:
shearForce.addParticle(atom.idx, [inputs.fc_x])
system.addForce(shearForce)
if inputs.divide > 0:
# we divide the box in two by applying a strong force on the atoms in the middle of the box.
holdForce = CustomExternalForce('fc_divide*(px^2+py^2+pz^2); \
px=min(dx, boxlx-dx); \
py=min(dy, boxly-dy); \
pz=min(dz, boxlz-dz); \
dx=abs(x-x0); \
dy=abs(y-y0); \
dz=abs(z-z0);')
holdForce.addGlobalParameter('boxlx', boxlx)
holdForce.addGlobalParameter('boxly', boxly)
holdForce.addGlobalParameter('boxlz', boxlz)
holdForce.addPerParticleParameter('fc_divide')
holdForce.addPerParticleParameter('x0')
holdForce.addPerParticleParameter('y0')
holdForce.addPerParticleParameter('z0')
# now we check every atom and if the atom is in a certain z range, we add it.
for line in open('restraints/water_pos.txt', 'r'):
segments = line.strip().split()
atom1 = int(segments[0])
xpos = crd.positions[atom1].value_in_unit(angstroms)[0]/10
ypos = crd.positions[atom1].value_in_unit(angstroms)[1]/10
zpos = crd.positions[atom1].value_in_unit(angstroms)[2]/10
holdForce.addParticle(atom1, [100, xpos, ypos, zpos])
system.addForce(holdForce)
return system
omm_barostat:
Code: Select all
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
def barostat(system, inputs):
if inputs.p_type == 'isotropic':
barostat = MonteCarloBarostat( inputs.p_ref*bar, inputs.temp*kelvin )
if inputs.p_type == 'membrane':
inputs.p_tens = inputs.p_tens*10.0
barostat = MonteCarloMembraneBarostat( inputs.p_ref*bar, inputs.p_tens*bar*nanometers, inputs.temp*kelvin, inputs.p_XYMode, inputs.p_ZMode, inputs.p_freq )
if inputs.p_type == 'anisotropic':
barostat = MonteCarloAnisotropicBarostat( inputs.p_ref*bar, inputs.temp*kelvin, inputs.p_scale[0], inputs.p_scale[1], inputs.p_scale[2], inputs.p_freq )
system.addForce(barostat)
return system
The file to write force and velocity:
omm_reporters
Code: Select all
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
class ForceReporter(object):
def __init__(self, file, reportInterval):
self._out = open(file, 'w')
self._reportInterval = reportInterval
def __del__(self):
self._out.close()
def describeNextReport(self, simulation):
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, False, False, True, False, None)
def report(self, simulation, state):
forces = state.getForces().value_in_unit(kilojoules/mole/nanometer)
for f in forces:
self._out.write('%g %g %g\n' % (f[0], f[1], f[2]))
class VelReporter(object):
def __init__(self, file, reportInterval):
self._out = open(file, 'w')
self._reportInterval = reportInterval
def __del__(self):
self._out.close()
def describeNextReport(self, simulation):
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, False, True, False, False, None)
def report(self, simulation, state):
velocities = state.getVelocities().value_in_unit(nanometer/picosecond)
for v in velocities:
self._out.write('%g %g %g\n' % (v[0], v[1], v[2]))
vfswitching file:
omm_vfswitch
Code: Select all
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
def vfswitch(system, psf, inputs):
r_on = inputs.r_on
r_off = inputs.r_off
# custom nonbonded force for force-switch
chknbfix = False
for force in system.getForces():
if isinstance(force, NonbondedForce):
nonbonded = force
if isinstance(force, CustomNonbondedForce):
nbfix = force
chknbfix = True
# vfswitch
vfswitch = CustomNonbondedForce('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
+step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
-step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
cr6 = ccnbb*ofdif3*rjunk3; \
cr12 = ccnba*ofdif6*rjunk6; \
rjunk3 = r3-recof3; \
rjunk6 = tr6-recof6; \
r3 = r1*tr2; \
r1 = sqrt(tr2); \
tr6 = tr2 * tr2 * tr2; \
tr2 = 1.0/s2; \
s2 = r*r; \
ccnbb = 4.0*epsilon*sigma^6; \
ccnba = 4.0*epsilon*sigma^12; \
sigma = sigma1+sigma2; \
epsilon = epsilon1*epsilon2; \
onoff3 = recof3/on3; \
onoff6 = recof6/on6; \
ofdif3 = off3/(off3 - on3); \
ofdif6 = off6/(off6 - on6); \
recof3 = 1.0/off3; \
on6 = on3*on3; \
on3 = c2onnb*Ron; \
recof6 = 1.0/off6; \
off6 = off3*off3; \
off3 = c2ofnb*Roff; \
c2ofnb = Roff*Roff; \
c2onnb = Ron*Ron; \
Ron = %f; \
Roff = %f;' % (r_on, r_off) )
vfswitch.addPerParticleParameter('sigma')
vfswitch.addPerParticleParameter('epsilon')
vfswitch.setNonbondedMethod(vfswitch.CutoffPeriodic)
vfswitch.setCutoffDistance(nonbonded.getCutoffDistance())
for i in range(nonbonded.getNumParticles()):
chg, sig, eps = nonbonded.getParticleParameters(i)
nonbonded.setParticleParameters(i, chg, 0.0, 0.0) # zero-out LJ
sig = sig*0.5
eps = eps**0.5
vfswitch.addParticle([sig, eps])
for i in range(nonbonded.getNumExceptions()):
atom1, atom2 = nonbonded.getExceptionParameters(i)[:2]
vfswitch.addExclusion(atom1, atom2)
vfswitch.setForceGroup(psf.NONBONDED_FORCE_GROUP)
system.addForce(vfswitch)
# vfswitch14
vfswitch14 = CustomBondForce('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
+step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
-step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
cr6 = ccnbb*ofdif3*rjunk3; \
cr12 = ccnba*ofdif6*rjunk6; \
rjunk3 = r3-recof3; \
rjunk6 = tr6-recof6; \
r3 = r1*tr2; \
r1 = sqrt(tr2); \
tr6 = tr2 * tr2 * tr2; \
tr2 = 1.0/s2; \
s2 = r*r; \
ccnbb = 4.0*epsilon*sigma^6; \
ccnba = 4.0*epsilon*sigma^12; \
onoff3 = recof3/on3; \
onoff6 = recof6/on6; \
ofdif3 = off3/(off3 - on3); \
ofdif6 = off6/(off6 - on6); \
recof3 = 1.0/off3; \
on6 = on3*on3; \
on3 = c2onnb*Ron; \
recof6 = 1.0/off6; \
off6 = off3*off3; \
off3 = c2ofnb*Roff; \
c2ofnb = Roff*Roff; \
c2onnb = Ron*Ron; \
Ron = %f; \
Roff = %f;' % (r_on, r_off) )
vfswitch14.addPerBondParameter('sigma')
vfswitch14.addPerBondParameter('epsilon')
for i in range(nonbonded.getNumExceptions()):
atom1, atom2, chg, sig, eps = nonbonded.getExceptionParameters(i)
nonbonded.setExceptionParameters(i, atom1, atom2, chg, 0.0, 0.0) # zero-out LJ14
vfswitch14.addBond(atom1, atom2, [sig, eps])
system.addForce(vfswitch14)
# vfswitch_NBFIX
if chknbfix:
nbfix.setEnergyFunction('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
+step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
-step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
cr6 = ccnbb*ofdif3*rjunk3; \
cr12 = ccnba*ofdif6*rjunk6; \
rjunk3 = r3-recof3; \
rjunk6 = tr6-recof6; \
r3 = r1*tr2; \
r1 = sqrt(tr2); \
tr6 = tr2 * tr2 * tr2; \
tr2 = 1.0/s2; \
s2 = r*r; \
ccnbb = bcoef(type1, type2); \
ccnba = acoef(type1, type2)^2; \
onoff3 = recof3/on3; \
onoff6 = recof6/on6; \
ofdif3 = off3/(off3 - on3); \
ofdif6 = off6/(off6 - on6); \
recof3 = 1.0/off3; \
on6 = on3*on3; \
on3 = c2onnb*Ron; \
recof6 = 1.0/off6; \
off6 = off3*off3; \
off3 = c2ofnb*Roff; \
c2ofnb = Roff*Roff; \
c2onnb = Ron*Ron; \
Ron = %f; \
Roff = %f;' % (r_on, r_off) )
return system
Then, I run the simulation with this file:
openmm_run
Code: Select all
from __future__ import print_function
import argparse
import sys
import os
from omm_readinputs import *
from omm_readparams import *
from omm_vfswitch import *
from omm_barostat import *
from omm_restraints import *
from omm_rewrap import *
from omm_reporters import *
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *
parser = argparse.ArgumentParser()
parser.add_argument('-i', dest='inpfile', help='Input parameter file', required=True)
parser.add_argument('-p', dest='psffile', help='Input CHARMM PSF file', required=True)
parser.add_argument('-c', dest='crdfile', help='Input CHARMM CRD file', required=True)
parser.add_argument('-t', dest='toppar', help='Input CHARMM-GUI toppar stream file', required=True)
parser.add_argument('-b', dest='sysinfo', help='Input CHARMM-GUI sysinfo stream file (optional)', default=None)
parser.add_argument('-icrst', metavar='RSTFILE', dest='icrst', help='Input CHARMM RST file (optional)', default=None)
parser.add_argument('-irst', metavar='RSTFILE', dest='irst', help='Input restart file (optional)', default=None)
parser.add_argument('-ichk', metavar='CHKFILE', dest='ichk', help='Input checkpoint file (optional)', default=None)
parser.add_argument('-opdb', metavar='PDBFILE', dest='opdb', help='Output PDB file (optional)', default=None)
parser.add_argument('-orst', metavar='RSTFILE', dest='orst', help='Output restart file (optional)', default=None)
parser.add_argument('-ochk', metavar='CHKFILE', dest='ochk', help='Output checkpoint file (optional)', default=None)
parser.add_argument('-odcd', metavar='DCDFILE', dest='odcd', help='Output trajectory file (optional)', default=None)
parser.add_argument('-rewrap', dest='rewrap', help='Re-wrap the coordinates in a molecular basis (optional)', action='store_true', default=False)
args = parser.parse_args()
# Load parameters
print("Loading parameters")
inputs = read_inputs(args.inpfile)
params = read_params(args.toppar)
psf = read_psf(args.psffile)
crd = read_crd(args.crdfile)
if args.sysinfo:
psf = read_box(psf, args.sysinfo)
else:
psf = gen_box(psf, crd)
vel_file = args.odcd[:-3]+"txt"
# Build system
system = psf.createSystem(params, nonbondedMethod=inputs.coulomb,
nonbondedCutoff=inputs.r_off*nanometers,
constraints=inputs.cons,
ewaldErrorTolerance=inputs.ewald_Tol)
if inputs.vdw == 'Force-switch': system = vfswitch(system, psf, inputs)
if inputs.pcouple == 'yes': system = barostat(system, inputs)
if inputs.rest == 'yes': system = restraints(system, psf, crd, inputs)
integrator = LangevinIntegrator(inputs.temp*kelvin, inputs.fric_coeff/picosecond, inputs.dt*picoseconds)
exit
# Set platform
platform = Platform.getPlatformByName('CUDA')
prop = dict(CudaPrecision='single', DisablePmeStream='true')
# Build simulation context
simulation = Simulation(psf.topology, system, integrator, platform, prop)
simulation.context.setPositions(crd.positions)
if args.icrst:
charmm_rst = read_charmm_rst(args.icrst)
simulation.context.setPositions(charmm_rst.positions)
simulation.context.setVelocities(charmm_rst.velocities)
simulation.context.setPeriodicBoxVectors(charmm_rst.box[0], charmm_rst.box[1], charmm_rst.box[2])
if args.irst:
with open(args.irst, 'r') as f:
simulation.context.setState(XmlSerializer.deserialize(f.read()))
if args.ichk:
with open(args.ichk, 'rb') as f:
simulation.context.loadCheckpoint(f.read())
# Re-wrap
if args.rewrap:
simulation = rewrap(simulation)
# Calculate initial system energy
print("\nInitial system energy")
print(simulation.context.getState(getEnergy=True).getPotentialEnergy())
# Energy minimization
if inputs.mini_nstep > 0:
print("\nEnergy minimization: %s steps" % inputs.mini_nstep)
simulation.minimizeEnergy(tolerance=inputs.mini_Tol*kilojoule/mole, maxIterations=inputs.mini_nstep)
print(simulation.context.getState(getEnergy=True).getPotentialEnergy())
# Generate initial velocities
if inputs.gen_vel == 'yes':
print("\nGenerate initial velocities")
if inputs.gen_seed:
simulation.context.setVelocitiesToTemperature(inputs.gen_temp, inputs.gen_seed)
else:
simulation.context.setVelocitiesToTemperature(inputs.gen_temp)
# set up velocity reporter
velocity_dis = simulation.context.getState( getVelocities=True )
velocity_reporter = VelReporter(vel_file, inputs.nstdcd)
# Production
if inputs.nstep > 0:
print("\nMD run: %s steps" % inputs.nstep)
if inputs.nstdcd > 0:
if not args.odcd: args.odcd = 'output.dcd'
simulation.reporters.append(DCDReporter(args.odcd, inputs.nstdcd))
simulation.reporters.append(velocity_reporter)
simulation.reporters.append(
StateDataReporter(sys.stdout, inputs.nstout, step=True, time=True, potentialEnergy=True, temperature=True, progress=True,
remainingTime=True, speed=True, totalSteps=inputs.nstep, separator='\t')
)
simulation.step(inputs.nstep)
# Write restart file
if not (args.orst or args.ochk): args.orst = 'output.rst'
if args.orst:
state = simulation.context.getState( getPositions=True )
with open(args.orst, 'w') as f:
f.write(XmlSerializer.serialize(state))
if args.ochk:
with open(args.ochk, 'wb') as f:
f.write(simulation.context.createCheckpoint())
if args.opdb:
crd = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(psf.topology, crd, open(args.opdb, 'w'))