Hello all and Charles,
I am new to XPLOR-NIH (As well as computation chemistry in general). I have
three questions. I have been attempting to run an experiment using the
refine.py script. The type of experiment I am attempting to run is thus: Begin
from a starting structure which is believed to be within 2-4 Angstroms RMSD
from the actual structure, and do a structure refinement using limited distance
restraints, torsion angle restraints, and chemical shift anisotropy restraints.
Is this a reasonable experiment to attempt with XPLOR-NIH? The script I have
currently been using is posted below (I have lowered the starting temperature
and step, and attempted to comment out parts I do not believe I need.) Are
there any glaring errors within? Lastly, my protein flies apart - is this being
caused by very limited restraints; would there be any way around that issue?
xplor.requireVersion("2.18")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2005/05/10
# Modified for use AB July 2009
xplor.parseArguments() # check for typos on the command-line
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=3
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('protG.seq')
#protocol.initStruct("g_new.psf") # - or from file
# generate a random extended structure with correct covalent geometry
# saves the generated structure in the indicated file for faster startup
# next time.
#
#protocol.genExtendedStructure("gb1_extended_%d.pdb" %
# protocol.initialRandomSeed())
# or read an existing model
#
protocol.loadPDB("Thiodecoy.pdb")
xplor.simulation.deleteAtoms("not known")
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
#
# a PotList contains a list of potential terms. This is used to specify which
# terms are active during refinement.
#
from potList import PotList
potList = PotList()
# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams
rampedParams=[]
highTempParams=[]
# compare atomic Cartesian rmsd with a reference structure
# backbone and heavy atom RMSDs will be printed in the output
# structure files
#
# I would like to put this back in later -AB
from posDiffPotTools import create_PosDiffPot
refRMSD = create_PosDiffPot("refRMSD","name CA",
pdbFile='ThioXRAY.pdb')
# cmpSel="not name H*")
# orientation Tensor - used with the dipolar coupling term
# one for each medium
# For each medium, specify a name, and initial values of Da, Rh.
#
from varTensorTools import create_VarTensor
media={}
# medium Da rhombicity
for (medium,Da,Rh) in [ ('t', -6.5, 0.62),
('b', -9.9, 0.23) ]:
oTensor = create_VarTensor(medium)
oTensor.setDa(Da)
oTensor.setRh(Rh)
media[medium] = oTensor
pass
# dipolar coupling restraints for protein amide NH.
#
# collect all RDCs in the rdcs PotList
#
# RDC scaling. Three possible contributions.
# 1) gamma_A * gamma_B / r_AB^3 prefactor. So that the same Da can be used
# for different expts. in the same medium. Sometimes the data is
# prescaled so that this is not needed. scale_toNH() is used for this.
# Note that if the expt. data has been prescaled, the values for rdc rmsd
# reported in the output will relative to the scaled values- not the expt.
# values.
# 2) expt. error scaling. Used here. A scale factor equal to 1/err^2
# (relative to that for NH) is used.
# 3) sometimes the reciprocal of the Da^2 is used if there is a large
# spread in Da values. Not used here.
#
# from rdcPotTools import create_RDCPot, scale_toNH
# rdcs = PotList('rdc')
# for (medium,expt,file, scale) in \
# [('t','NH' ,'tmv107_nh.tbl' ,1),
# ('t','NCO','tmv107_nc.tbl' ,.05),
# ('t','HNC','tmv107_hnc.tbl' ,.108),
# ('b','NH' ,'bicelles_new_nh.tbl' ,1),
# ('b','NCO','bicelles_new_nc.tbl' ,.05),
# ('b','HNC','bicelles_new_hnc.tbl',.108)
# ]:
# rdc = create_RDCPot("%s_%s"%(medium,expt),file,media[medium])
#1) scale prefactor relative to NH
# see python/rdcPotTools.py for exact calculation
# scale_toNH(rdc) - not needed for these datasets -
# but non-NH reported rmsd values will be wrong.
#3) Da rescaling factor (separate multiplicative factor)
# scale *= ( 1. / rdc.oTensor.Da(0) )**2
# rdc.setScale(scale)
# rdc.setShowAllRestraints(1) #all restraints are printed during analysis
# rdc.setThreshold(1.5) # in Hz
# rdcs.append(rdc)
# pass
# potList.append(rdcs)
# rampedParams.append( MultRamp(0.05,5.0, "rdcs.setScale( VALUE )") )
# calc. initial tensor orientation
#
from varTensorTools import calcTensorOrientation
for medium in media.values():
calcTensorOrientation(medium)
pass
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"rstrtsThio74_108.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
# pot.setPotType("soft") - if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
rampedParams.append( MultRamp(2,30, "noe.setScale( VALUE )") )
# set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","jna_coup.tbl",
# A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
# Set up dihedral angles
from xplorPot import XplorPot
protocol.initDihedrals("TorsionThio74_108.tbl",
useDefaults=0)
potList.append( XplorPot('CDIH') )
highTempParams.append( StaticRamp("potList['CDIH'].setScale(10)") )
rampedParams.append( StaticRamp("potList['CDIH'].setScale(200)") )
# set custom values of threshold values for violation calculation
#
potList['CDIH'].setThreshold( 5 ) #5 degrees is the default value, though
# gyration volume term
#
# gyration volume term
#
from gyrPotTools import create_GyrPot
gyr = create_GyrPot("Vgyr",
"resid 1:56") # selection should exclude disordered tails
potList.append(gyr)
rampedParams.append( MultRamp(.002,1,"gyr.setScale(VALUE)") )
# hbda - distance/angle bb hbond term
#
# protocol.initHBDA('hbda.tbl')
# potList.append( XplorPot('HBDA') )
#Rama torsion angle database
#
protocol.initRamaDatabase()
potList.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond()") )
rampedParams.append( MultRamp(0.9,0.8,
"command('param nbonds repel VALUE end end')") )
rampedParams.append( MultRamp(.004,4,
"command('param nbonds rcon VALUE end end')") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
rcon=0.004,
tolerance=45,
repel=1.2,
onlyCA=1)""") )
potList.append( XplorPot("BOND") )
potList.append( XplorPot("ANGL") )
potList['ANGL'].setThreshold( 5 )
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )
potList.append( XplorPot("IMPR") )
potList['IMPR'].setThreshold( 5 )
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )
# Give atoms uniform weights, except for the anisotropy axis
#
from atomAction import SetProperty
import varTensorTools
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
varTensorTools.massSetup(media.values(),300)
AtomSel("all ").apply( SetProperty("fric",10.) )
# IVM setup
# the IVM is used for performing dynamics and minimization in torsion-angle
# space, and in Cartesian space.
#
from ivm import IVM
dyn = IVM()
ivm = IVM()
ivm
### Script appending from
http://nmr.cit.nih.gov/pipermail/xplor-nih/2008-April/000970.html
### This should create the trajectory in the VMD viewer
from vmdInter import VMDInter
from vmdInter import VMDTraj
vmd=VMDInter()
vmdStruct=vmd.makeObj("Thio74_108"); vmdStruct.bonds( AtomSel("all") )
# vmdStruct=vmd.makeObj("dynStruct"); vmdStruct.bonds( AtomSel("not hydro") )
vmdTraj=VMDTraj(vmdStruct); vmdTraj.saveInterval=5
dyn.setTrajectory( vmdTraj )
# initially minimize in Cartesian space with only the covalent constraints.
# Note that bonds, angles and many impropers can't change with the
# internal torsion-angle dynamics
# breaks bonds topologically - doesn't change force field
#
#dyn.potList().add( XplorPot("BOND") )
#dyn.potList().add( XplorPot("ANGL") )
#dyn.potList().add( XplorPot("IMPR") )
#
#dyn.breakAllBondsIn("not resname ANI")
#import varTensorTools
#for m in media.values():
# m.setFreedom("fix") #fix tensor parameters
# varTensorTools.topologySetup(dyn,m) #setup tensor topology
#
#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
for m in media.values():
# m.setFreedom("fixDa, fixRh") #fix tensor Rh, Da, vary orientation
m.setFreedom("varyDa, varyRh") #vary tensor Rh, Da, vary orientation
protocol.torsionTopology(dyn,oTensors=media.values())
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
for m in media.values():
m.setFreedom("varyDa, varyRh") #allow all tensor parameters float here
pass
protocol.cartesianTopology(minc,oTensors=media.values())
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 1000. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=1,
tempStep =10.0,
ivm=dyn,
rampedParams = rampedParams)
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>0:
return False
if potList['rdc'].rms()>1.2: #this might be tightened some
return False
if potList['CDIH'].violations()>0:
return False
if potList['BOND'].violations()>0:
return False
if potList['ANGL'].violations()>0:
return False
if potList['IMPR'].violations()>1:
return False
return True
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
# initialize parameters for high temp dynamics.
InitialParams( rampedParams )
# high-temp dynamics setup - only need to specify parameters which
# differfrom initial values in rampedParams
InitialParams( highTempParams )
# high temp dynamics
#
protocol.initDynamics(dyn,
potList=potList, # potential terms to use
bathTemp=init_t,
initVelocities=1,
finalTime=10, # stops at 10ps or 5000 steps
numSteps=5000, # whichever comes first
printInterval=100)
dyn.setETolerance( init_t/100 ) #used to det. stepsize. default: t/1000
dyn.run()
# initialize parameters for cooling loop
InitialParams( rampedParams )
# initialize integrator for simulated annealing
#
protocol.initDynamics(dyn,
potList=potList,
numSteps=100, #at each temp: 100 steps or
finalTime=.2 , # .2ps, whichever is less
printInterval=100)
# perform simulated annealing
#
cool.run()
# final torsion angle minimization
#
protocol.initMinimize(dyn,
printInterval=50)
dyn.run()
# final all- atom minimization
#
protocol.initMinimize(minc,
potList=potList,
dEPred=10)
minc.run()
#do analysis and write structure
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR']],
# noe,rdcs,potList['CDIH']],
averageCrossTerms=refRMSD,).run()
# averageTopFraction=0.5, #report only on best 50% of structs
#averageAccept=accept, #only use structures which pass accept()
# averageContext=FinalParams(rampedParams),
# averageFilename = "SCRIPT_ave.pdb", #generate regularized ave structure
# averageFitSel="name CA",
# averageCompSel="not resname ANI and not name H*" ).run()
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
http://dcb.cit.nih.gov/pipermail/xplor-nih/attachments/20090726/cb37bd38/attachment-0001.html
