Hello,
I'm trying to run a version of the anneal.py from de gb1 example. But it,
gives me the following error.
~/Documents/sim_hepcidine$ pyXplor anneal.py
segid: minResid: 1 maxResid: 24
reading directional HB PMF - 310 right
reading directional HB PMF - 310 left
reading directional HB PMF - a N-t & cent
reading directional HB PMF - a C-t & isol
reading directional HB PMF - b anti cent
reading directional HB PMF - b anti long
reading directional HB PMF - b para cent
reading directional HB PMF - b para edge
reading directional HB PMF - lr isolated
reading directional HB PMF - g turn right
reading directional HB PMF - g turn left
calculating directional HB PMF forces
i/i+3i/i+4 : directional forces done
reading linearity HB PMF - i/i+3
reading linearity HB PMF - i/i+4
reading linearity HB PMF - lr
reading linearity HB PMF - i/i+2
calculating linearity HB PMF forces
linearity HB PMF forces done
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BOMLEV= 0 reached. Program execution will be terminated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine DIE called . Terminating
Here is the script I'm using. I also made a psf file from another script so
I don't know if ti is a topology problem.
Thanks in advance
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
xplor.requireVersion("2.34")
#
# 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 2009/07/24
#
# 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 = "out_anneal/SCRIPT_STRUCTURE.sa"
numberOfStructures=200 #usually you want to create at least 20
# if quick:
# numberOfStructures=3
# pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed() #set random seed - by time
command = xplor.command
protocol.initParams(("protein-4.0-mod.par"))
protocol.initStruct("hepcidin_free.psf")
# generate PSF data from sequence and initialize the correct parameters.
#
# from psfGen import seqToPSF
# seqToPSF('protG.seq')
# # 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())
# #
# # 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' ,'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 *= ( 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 list(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,"noe_restr.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
dihedralRestraintFilename="dihe_restr.tbl"
protocol.initDihedrals(dihedralRestraintFilename,
#useDefaults=False # by default, symmetric sidechain
# restraints are included
)
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('hbda.tbl')
# potList.append( XplorPot('HBDA') )
# hbdb - distance/angle bb hbond term
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
# New torsion angle database potential
#
from torsionDBPotTools import create_TorsionDBPot
torsionDB = create_TorsionDBPot('torsionDB', system='protein')
potList.append( torsionDB )
rampedParams.append( MultRamp(.002,2,"torsionDB.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(.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
#
protocol.massSetup()
# 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 list(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)
# for m in list(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 = 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)
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
# calc. initial tensor orientation
#
# from varTensorTools import calcTensorOrientation
# for medium in list(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 when function returns
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
doWriteStructures=True,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=True,
averageTopFraction=0.1, #report stats on best 50% of structs
averageContext=FinalParams(rampedParams),
averageFitSel=("name CA or name C or name N or name O"),
# averageCrossTerms=refRMSD,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH']],
averageFilename="out_sa/SCRIPT_ave.pdb",
averagePotList=potList).run()
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
