Hi Charles,
Thank you for your reply. I did modify the
'newrefine.py'(attached file) but got the following errors:
InternalDynamics::step:
large timestep detected. Halving.
%atoms "ALT2-35 -CYSP-SG2 " and
"ALT3-35 -CYSP-SG2 " only 0.01 A apart
InternalDynamics::step:
large timestep detected. Halving.
InternalDynamics::step: large
timestep detected. Halving.
InternalDynamics::step: large timestep
detected. Halving.
InternalDynamics::step: large timestep detected.
Halving.
InternalDynamics::step: large timestep detected. Halving.
InternalDynamics::step:
large timestep detected. Halving.
InternalDynamics::step: large
timestep detected. Halving.
InternalDynamics::step: large timestep
detected. Halving.
InternalDynamics::step: large timestep detected.
Halving.
InternalDynamics::step: large timestep detected. Halving.
InternalDynamics::step:
large timestep detected. Halving.
InternalDynamics::step: large
timestep detected. Halving.
InternalDynamics::step: stepsize too
small.
Traceback (most recent call last):
File "<string>",
line 2, in <module>
File
"/opt/xplor-nih-2.26/python/trace.py", line 180, in run
exec cmd
in dict, dict
File "<string>", line 1, in <module>
File "newRefine.py", line 634, in <module>
averagePotList=potList).run()
File
"/opt/xplor-nih-2.26/python/simulationTools.py", line 281, in run
s.structLoopAction(s)
File "newRefine.py", line 569, in
calcOneStructure
dyn.run()
File
"/opt/xplor-nih-2.26/python/ivm.py", line 465, in run
(done,s.stepsize_) = PublicIVM.step(s,s.stepsize());
File
"/opt/xplor-nih-2.26/python/wrappers/publicIVM.py", line 136, in step
def step(self, *args, **kwargs): return _publicIVM.PublicIVM_step(self,
*args, **kwargs)
_publicIVM.IVMError: IVM error: stepsize too small
PyInterp::command:
error executing: >execfile('newRefine.py')<
HEAP: maximum
use= 22720953 current use= 19890309
X-PLOR: total CPU time=
19.0100 s
X-PLOR: entry time at 19:24:50 7-Sep-10
X-PLOR: exit
time at 19:25:13 7-Sep-10
Also, I am not quite understand
the clock atom and tauc part. How are the initial values (resid 777-783;
taucapp=2.9,3.0) chosen?
I want to calculate form extended
structures, and the modified 'anneal.py'(attached file) did work. But in
the output file, there are VDW violations between conformers of the
same cysp like these:
viol% amount index restraint name
---------------------------------------------------------------------------
100.0
0.96 0 ( ALT1 7 CYSP SG1 ) ( ALT2 7 CYSP SG2 )
100.0
1.01 1 ( ALT1 7 CYSP SG1 ) ( ALT3 7 CYSP SG2 )
100.0
1.02 2 ( ALT1 7 CYSP SG2 ) ( ALT2 7 CYSP SG1 )
100.0
1.85 3 ( ALT1 7 CYSP SG2 ) ( ALT2 7 CYSP SG2 )
100.0
1.18 4 ( ALT1 7 CYSP SG2 ) ( ALT2 7 CYSP CL1 )
Are
these normal? I use MOLMOL to calculate the distance between 'OS1' and
'HN' atoms. The calculated ones are larger than the max distance in the
input '.tbl' files, but they were not reported in the '.viols' files.
Thank
you very much for your help!
xplor.requireVersion("2.14.4")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# this script performs annealing from an extended structure
#
# CDS 2005/05/10
#
# this checks for typos on the command-line. User-customized arguments can
# also be specified.
#
xplor.parseArguments()
# filename for output structures. This string must contain the STRUCTURE
# literal so that each calculated structure has a unique name. The SCRIPT
# literal is replaced by this filename (or stdin if redirected using <),
# but it is optional.
#
outFilename = "tatA_STRUCTURE.sa"
#numberOfStructures=1 #usually you want to create at least 20
numberOfStructures=100 #usually you want to create at least 20
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed() #set random seed - by time
nconf = 3
command = xplor.command
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('xh.seq')
# make cys-MTSL substitutions
#
#import os
#try:
# os.unlink(successFile)
#except:
# pass
xplor.parseArguments() # check for typos on the command-line
import protocol
protocol.initTopology(('protein'))
protocol.initParams(('protein'))
#import psfGen
#psfGen.seqToPSF(seq,seqType='prot',startResid=1)
protocol.initStruct("../7Inputs/generate.psf")
protocol.initCoords("../7Inputs/random3Alts1.pdb")
# atoms in this selection will not be moved at all
#
#fixGeom=AtomSel('''(not (name ca or name n or name c or
# name o or name hn or name ha))
# and (resid 10 or resid 35 or resid 33))''')
# generate a random extended structure with correct covalent geometry
# saves the generated structure in the indicated file for faster startup
# next time.
#
protocol.genExtendedStructure()
#
# 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
#
#from posDiffPotTools import create_PosDiffPot
#refRMSD = create_PosDiffPot("refRMSD","name CA or name C or name N",
# pdbFile='g_xray.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' ,'rdc_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 *= ( 9.9 / rdc.oTensor.Da(0) )**2
# rdc.setScale(scale)
# rdcs.append(rdc)
# pass
#potList.append(rdcs)
#rampedParams.append( MultRamp(0.01,1.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,"../9Restraints/noe_cysp_N.tbl"),('pre',1,"../9Restraints/pre_cysp.tbl"),('hbond',1,"../9Restraints/hbond.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("../9Restraints/dihe.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 )
# radius of gyration term
#
#protocol.initCollapse("resid 1:56",
# Rtarget=10.16)
#potList.append( XplorPot('COLL') )
# hbda - distance/angle bb hbond term
#
#protocol.initHBDA('hbond.tbl')
#potList.append( XplorPot('HBDA') )
from ivm import IVM
dyn = IVM()
# group MTSL
for res in ([ 10, 35, 33, 7]) :
for i in range(0,nconf+1) :
dyn.group("""segid ALT%-d and ((resid %d and (name *S* or name SG* or name CL1 or name HL*)))""" % (i, res))
pass
pass
def setConstraints():
#number of conformers (support up to 10 conformers)
scaleVdw = 1.0 / float( nconf )
xplor.command("""
constraints
! interaction (all) (all)
inter = (segid " ") (segid " ")
inter = (segid "ALT1") (segid "ALT1")
inter = (segid "ALT2") (segid "ALT2")
inter = (segid "ALT3") (segid "ALT3")
inter = (segid "ALT4") (segid "ALT4")
inter = (segid "ALT5") (segid "ALT5")
inter = (segid "ALT6") (segid "ALT6")
inter = (segid "ALT7") (segid "ALT7")
inter = (segid "ALT8") (segid "ALT8")
inter = (segid "ALT9") (segid "ALT9")
inter = (segid "ALT0") (segid "ALT0")
weights * 1 end
inter = (segid " ") (segid "ALT1")
inter = (segid " ") (segid "ALT2")
inter = (segid " ") (segid "ALT3")
inter = (segid " ") (segid "ALT4")
inter = (segid " ") (segid "ALT5")
inter = (segid " ") (segid "ALT6")
inter = (segid " ") (segid "ALT7")
inter = (segid " ") (segid "ALT8")
inter = (segid " ") (segid "ALT9")
inter = (segid " ") (segid "ALT0")
weights * 1 vdw %f end
end""" % scaleVdw )
return
#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( StaticRamp("setConstraints()") )
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()
# initialize ivm topology for torsion-angle dynamics
#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())
protocol.torsionTopology(dyn)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
for res in ([ 10, 35, 33, 7]) :
for i in range(0,nconf+1) :
minc.group("""(segid ALT%-d and ((resid %d and (name *S* or name SG* or name CL1 or name HL*))))""" % (i, res))
pass
pass
#for m in media.values():
# m.setFreedom("varyDa, varyRh") #allow all tensor parameters float here
# pass
#protocol.cartesianTopology(minc,oTensors=media.values())
protocol.cartesianTopology(minc)
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3500. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =12.5,
ivm=dyn,
rampedParams = rampedParams)
#cart_cool is for optional cartesian-space cooling
cart_cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =12.5,
ivm=minc,
rampedParams = rampedParams)
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
# generate a new structure with randomized torsion angles
#
from monteCarlo import randomizeTorsions
randomizeTorsions(dyn)
# calc. initial tensor orientation
#
#from varTensorTools import calcTensorOrientation
#for medium in media.values():
# calcTensorOrientation(medium)
# pass
# 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=800, # stops at 800ps or 8000 steps
numSteps=8000, # 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()
# optional cooling in Cartesian coordinates
#
protocol.initDynamics(minc,
potList=potList,
numSteps=100, #at each temp: 100 steps or
finalTime=.4 , # .2ps, whichever is less
printInterval=100)
#cart_cool.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,
averageTopFraction=1, #report stats on best 50% of structs
averageContext=FinalParams(rampedParams),
# averageCrossTerms=refRMSD,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
# noe,rdcs,potList['CDIH']],
noe,potList['CDIH']],
averagePotList=potList).run()
xplor.requireVersion("2.20")
# this checks for typos on the command-line. User-customized arguments can
# also be specified.
#
(opts,args) = xplor.parseArguments(["quick"])
quick=False
for opt in opts:
if opt[0]=="quick": #specify -quick to just test that the script runs
quick=True
pass
pass
# filename for output structures. This string must contain the STRUCTURE
# literal so that each calculated structure has a unique name. The SCRIPT
# literal is replaced by this filename (or stdin if redirected using <),
# but it is optional.
#
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=2 if quick else 4
# protocol module has many high-level helper functions.
#
import protocol, xplor
protocol.initRandomSeed() #set random seed - by time
# parameters
protocol.initTopology(("protein"))
protocol.initParams(("protein"))
# psf
protocol.initStruct("../7Inputs/generate.psf")
# coords
protocol.initCoords("../7Inputs/random3Alts1.pdb")
#xplor.simulation.deleteAtoms("resname ANI")
import protocol
#protocol.fixupCovalentGeom(useVDW=True,verbose=1,
# angleTol=2.5)
protocol.genExtendedStructure()
nconf = 3
import protocol, xplor
#
# 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()
#crossTerms = 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
#
#from posDiffPotTools import create_PosDiffPot
#refRMSD = create_PosDiffPot("refRMSD","name CA or name C or name N",
# pdbFile='g_xray.pdb',
# cmpSel="not name H*")
#crossTerms.append(refRMSD)
# 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 [ ('m1', 8.6, 0.61),
# ('m2', 8.5, 0.6 ),
# ('m3', 8.0, 0.59), ]:
# 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 \
# [('m1','NH' , 'inputs/dip_nh.tbl' ,1),
# ('m1','CAHA','inputs/dip_caha.tbl' ,1),
# ('m1','HNC', 'inputs/dip_hnc.tbl' ,0.108),
# ('m1','NCO', 'inputs/dip_nc.tbl' ,0.05),
# ('m2','HNp', 'inputs/17dip_protnh.tbl' ,1),
# ('m2','HNd', 'inputs/junk_dna_nh.tbl' ,1),
# ('m2','CHd',('inputs/final_dna_sug-new.tbl',
# 'inputs/final_dna_aro.tbl' ),1),
# ('m3','CHg', 'inputs/CAH_gabriel.tbl' ,1),
# ('m3','HHp', 'inputs/wt_HH_dipy-new.tbl' ,0),
# ('m3','HHn', 'inputs/HH_dip_negx-new.tbl' ,0),
# ]:
# rdc = create_RDCPot("%s"%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 *= ( 9.9 / rdc.oTensor.Da(0) )**2
# rdc.setScale(scale)
# rdcs.append(rdc)
# pass
#for name in ('HHp', 'HHn'):
# pot = rdcs[name]
# scale_toNH(pot)
# pot.setGyroA( 10.46847258 )
# pot.setUseSign(False)
# pot.setPotType('square') #scales rdc values and vary distance dependence
# pass
#from varTensorTools import calcTensorOrientation
#for medium in media.values():
# calcTensorOrientation(medium)
# pass
#potList.append(rdcs)
#rampedParams.append( MultRamp(0.01,1.0, "rdcs.setScale( VALUE )") )
#
# PRE restraints
#
from prePotTools import create_PREPot
pre = PotList('pre')
# add clock pseudo-atoms
#
#from prePotTools import addClockAtoms
#for resid in range(777,783):
# addClockAtoms(resid)
# pass
#
# for CYSP10
#
cysp10 = PotList('cysp10')
for (name,file,assignType) in [
("Cysp10hn","pre_cysp_10.tbl","normal")
]:
p = create_PREPot(name,"../9Restraints/"+file,assignType,600)
cysp10.append(p)
pass
pre.append(cysp10)
#
# for CYSP35
#
cysp35 = PotList('cysp35')
for (name,file,assignType) in [
("Cysp35hn","pre_cysp_35.tbl","normal")
]:
p = create_PREPot(name,"../9Restraints/"+file,assignType,600)
cysp35.append(p)
pass
pre.append(cysp35)
#
# for CYSP33
#
cysp33 = PotList('cysp33')
for (name,file,assignType) in [
("Cysp33hn","pre_cysp_33.tbl","normal")
]:
p = create_PREPot(name,"../9Restraints/"+file,assignType,600)
cysp33.append(p)
pass
pre.append(cysp33)
#
# for CYSP7
#
cysp7 = PotList('cysp7')
for (name,file,assignType) in [
("Cysp7hn","pre_cysp_7.tbl","normal")
]:
p = create_PREPot(name,"../9Restraints/"+file,assignType,600)
cysp7.append(p)
pass
pre.append(cysp7)
#
def runSBmode():
print 'configuring SB mode'
from prePotTools import setupSBmode
from simulationTools import flattenPotList
for p in flattenPotList(pre): setupSBmode(p)
#apparent tauc
tcappCysp10=2.9
tcappCysp35=3.0
for p in list(cysp10)+list(cysp7):
p.setTcType("fix")
p.setTauC( tcappCysp10 )
pass
for p in list(cysp33)+list(cysp35):
p.setTcType("fix")
p.setTauC( tcappCysp35 )
pass
# for p in list(es2at750)+list(es3at750):
# p.setTcType("opt")
# p.setTauC( tcappEs2 )
# p.setTcMin( tcappEs2 )
# p.setTcMax( tcappEs2 * 1.6 )
# for p in es4at800:
# p.setTcType("opt")
# p.setTauC( tcappEs4 )
# p.setTcMin( tcappEs4 )
# p.setTcMax( tcappEs4 * 1.6 )
# pass
return
bbPre = PotList('bbPre')
for p in (cysp10['Cysp10hn'],
cysp7['Cysp7hn'],
cysp33['Cysp33hn'],
cysp35['Cysp35hn']):
bbPre.append(p)
pass
#crossTerms.append(bbPre)
#scPre = PotList('scPre')
#for p in (es2at500['Es2sc500'], es2at500['Es2ns500'], es2at750['Es2sc750'],
# es3at500['Es3sc500'], es3at500['Es3ns500'],
# es3at750['Es3sc750'], es3at750['Es3ns750'],
# es4at500['Es4sc500'], es4at500['Es4ns500'],
# es4at800['Es4sc800'], es4at800['Es4ns800']):
# p.setScale(0.5)
# scPre.append(p)
# pass
#crossTerms.append(scPre)
def runSBMFmode():
print 'configuring SBMF mode'
from prePotTools import setupSBMFmode
from simulationTools import flattenPotList
for p in flattenPotList(pre): setupSBMFmode(p)
pre.calcEnergy()
# tcappEs2=2.9
# tcappEs4=3.0
# maxtcEs2 = tcappEs2 / es2at500['Es2hn500'].aveS2()
# maxtcEs4 = tcappEs4 / es4at500['Es4hn500'].aveS2()
# for p in list(es2at500)+list(es3at500):
# p.setTcType("opt")
# p.setTcMin( tcappEs2 )
# p.setTcMax( maxtcEs2 )
# pass
# for p in es4at500 :
# p.setTcType("opt")
# p.setTcMin( tcappEs4 )
# p.setTcMax( maxtcEs4 )
# pass
# for p in list(es2at750)+list(es3at750) :
# p.setTcType("opt")
# p.setTcMin( tcappEs2 )
# p.setTcMax( 1.6*maxtcEs2 )
# pass
# for p in es4at800 :
# p.setTcType("opt")
# p.setTcMin( tcappEs4 )
# p.setTcMax( 1.6*maxtcEs4 )
# pass
return
potList.append(pre)
rampedParams.append( MultRamp(0.05,1.0, "pre.setScale( VALUE )") )
rampedParams.append( StaticRamp("runSBMFmode()") )
highTempParams.append( StaticRamp("runSBmode()") )
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('noe',1,["../9Restraints/noe_cysp_N.tbl",
# "inputs/dna_trial1-new.tbl", # modified version of dna_trial1.tbl
# "inputs/endnoe_dnax.tbl",
# "inputs/pp_trial.tbl",
# "inputs/trial_inter3.tbl",
# "inputs/6_22nh_inter.tbl",
# "inputs/inter_hb.tbl",
# "inputs/ends_bdnazzz-new.tbl", # modified version of ends_bdnazzz.tbl
]),
('hbond',1,["../9Restraints/hbond.tbl",
])
]:
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 )") )
from xplorPot import XplorPot
# carbon 13 chemical shifts
#protocol.initCarb13("inputs/c13_shifts.tbl")
#potList.append(XplorPot("CARB"))
# set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","inputs/hnha_xplr.tbl",
# A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
#x_jCoup = create_JCoupPot("xjcoup","inputs/dna_pcoup.tbl",
# A=15.3,B=-6.2,C=1.5,phase=120.0)
#crossTerms.append(x_jCoup)
# Set up dihedral angles
protocol.initDihedrals(["../9Restraints/dihe.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 )
# radius of gyration term
#
#protocol.initCollapse("resid 4:81 or resid 105:112 or resid 117:124",
# Rtarget=15,scale=0.25)
#potList.append( XplorPot('COLL') )
#Rama torsion angle database
#
#protocol.initRamaDatabase('nucleic')
protocol.initRamaDatabase('protein')
potList.append( XplorPot('RAMA') )
ramaProtClasses=['phi_psi_ala', 'phi_psi_arg', 'phi_psi_asn', 'phi_psi_asp',
'phi_psi_chi1_arg', 'phi_psi_chi1_asn', 'phi_psi_chi1_asp',
'phi_psi_chi1_cyo', 'phi_psi_chi1_cyr', 'phi_psi_chi1_gln',
'phi_psi_chi1_glu', 'phi_psi_chi1_his', 'phi_psi_chi1_ile',
'phi_psi_chi1_leu', 'phi_psi_chi1_lys', 'phi_psi_chi1_met',
'phi_psi_chi1_phe', 'phi_psi_chi1_ser', 'phi_psi_chi1_thr',
'phi_psi_chi1_trp', 'phi_psi_chi1_tyr', 'phi_psi_chi1_val',
'phi_psi_csp', 'phi_psi_cyo', 'phi_psi_cyr', 'phi_psi_gln',
'phi_psi_glu', 'phi_psi_gly', 'phi_psi_his', 'phi_psi_ile',
'phi_psi_leu', 'phi_psi_lys', 'phi_psi_met', 'phi_psi_phe',
'phi_psi_pro', 'phi_psi_ser', 'phi_psi_thr', 'phi_psi_trp',
'phi_psi_tyr', 'phi_psi_val', 'phi_psi_xcp', 'phi_psi_xpr',
'chi1_chi2_arg', 'chi1_chi2_asn', 'chi1_chi2_asp',
'chi1_chi2_chi3_arg', 'chi1_chi2_chi3_gln',
'chi1_chi2_chi3_glu', 'chi1_chi2_chi3_lys',
'chi1_chi2_chi3_met', 'chi1_chi2_cyo', 'chi1_chi2_gln',
'chi1_chi2_glu', 'chi1_chi2_his', 'chi1_chi2_ile',
'chi1_chi2_leu', 'chi1_chi2_lys', 'chi1_chi2_met',
'chi1_chi2_phe', 'chi1_chi2_trp', 'chi1_chi2_tyr',
'chi2_chi3_arg', 'chi2_chi3_gln', 'chi2_chi3_glu',
'chi2_chi3_lys', 'chi2_chi3_met', 'chi4_arg', 'chi4_lys']
#FIXME: be able to tune nucleic fc separate from protein
rampedParams.append( MultRamp(1,1,"potList['RAMA'].setScale(VALUE)") )
rampedParams.append( MultRamp(.002,1,
"""for pclass in ramaProtClasses:
xplor.command('rama class %s force %f end ' % (pclass, VALUE))""") )
#initialize the aa-aa positional database
#xplor.command("@inputs/bbpd.setup")
#potList.append( XplorPot('ORIE') )
#planarity restraints
#xplor.command("restraints plane @inputs/dnaPlanBp.tbl end")
#potList.append(XplorPot("plan"))
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond()") )
rampedParams.append( StaticRamp("setConstraints()") )
rampedParams.append( MultRamp(0.9,0.78,
"xplor.command('param nbonds repel VALUE end end')") )
rampedParams.append( MultRamp(.004,4,
"xplor.command('param nbonds rcon VALUE end end')") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
rcon=1,
tolerance=45,
repel=1.2,
onlyCA=1)""") )
#highTempParams.append(
# StaticRamp(
# xplor.command("""constraints
# inter = (not (name ca or name p or name mn)) (all)
# weights * 1 angl 0.4 impr 0.1 vdw 0 elec 0 end
# inter = (name ca or name p or name mn) (name ca or name p or name mn)
# weights * 1 angl 0.4 impr 0.1 vdw 1.0 end
# end""") ) )
def setConstraints():
#number of conformers (support up to 10 conformers)
scaleVdw = 1.0 / float( nconf )
xplor.command("""
constraints
! interaction (all) (all)
inter = (segid " ") (segid " ")
inter = (segid "ALT1") (segid "ALT1")
inter = (segid "ALT2") (segid "ALT2")
inter = (segid "ALT3") (segid "ALT3")
inter = (segid "ALT4") (segid "ALT4")
inter = (segid "ALT5") (segid "ALT5")
inter = (segid "ALT6") (segid "ALT6")
inter = (segid "ALT7") (segid "ALT7")
inter = (segid "ALT8") (segid "ALT8")
inter = (segid "ALT9") (segid "ALT9")
inter = (segid "ALT0") (segid "ALT0")
weights * 1 end
inter = (segid " ") (segid "ALT1")
inter = (segid " ") (segid "ALT2")
inter = (segid " ") (segid "ALT3")
inter = (segid " ") (segid "ALT4")
inter = (segid " ") (segid "ALT5")
inter = (segid " ") (segid "ALT6")
inter = (segid " ") (segid "ALT7")
inter = (segid " ") (segid "ALT8")
inter = (segid " ") (segid "ALT9")
inter = (segid " ") (segid "ALT0")
weights * 1 vdw %f end
end""" % scaleVdw )
return
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("resname TAU ").apply(SetProperty("mass", 300.0))
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()
# group MTSL
for res in ([ 10, 35, 33, 7]) :
for i in range(0,nconf+1) :
dyn.group("""segid ALT%-d and ((resid %d and (name *S* or name SG* or name CL1 or name HL*)))""" % (i, res))
pass
pass
# allow fo tau_c optimization
from simulationTools import flattenPotList
for p in flattenPotList(pre):
p.setTcType("opt")
pass
# initialize ivm topology for torsion-angle dynamics
#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)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
# group EDTA-Mn
for res in ([ 10, 35, 33, 7]) :
for i in range(0,nconf+1) :
minc.group("""(segid ALT%-d and ((resid %d and (name *S* or name SG* or name CL1 or name HL*))))""" % (i, res))
pass
pass
#for m in media.values():
# m.setFreedom("varyDa, varyRh") #allow all tensor parameters float here
# pass
#protocol.cartesianTopology(minc)
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3000. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =100 if quick else 12.5,
ivm=dyn,
rampedParams = rampedParams)
def calcOneStructure(loopInfo):
# """ this function calculates a single structure, performs analysis on the
# structure, and then writes out a pdb file, with remarks.
# """
# calc. initial tensor orientation
#
# from varTensorTools import calcTensorOrientation
# for medium in media.values():
# calcTensorOrientation(medium)
# pass
# 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=800, # stops at 800ps or 8000 steps
numSteps=5 if quick else 8000, # 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()
# optional cooling in Cartesian coordinates
#
protocol.initDynamics(minc,
potList=potList,
numSteps=100, #at each temp: 100 steps or
finalTime=.4 , # .2ps, whichever is less
printInterval=100)
protocol.initMinimize(minc,
potList=potList,
dEPred=10)
minc.run()
# print close-contact atoms to a .bumps file
# this XPLOR command will honor the contraints interaction setup
# so that overlaps the copies of the tag(s) are not counted. The
# output in the .viols will report these overlaps.
xplor.command("""
distance
from=(not PSEUDO) to=(known)
cuton=0 cutoff=1.5 disp=print
output=%s
end
""" % (loopInfo.filename() + '.bumps') )
#do analysis and write structure
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averageTopFraction=1.0, #report stats on best 50% of structs
averageContext=FinalParams(rampedParams),
# averageCrossTerms=crossTerms,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH']],
averagePotList=potList).run()
_______________________________________________
Xplor-nih mailing list
[email protected]
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
