Dear Charles and xplor-users,
I am trying to use ensemble calculation for unstructured proteins. However,
the output structures are too compact than expected, so I did the following
test:
The simulated annealing python input script from the PRE distribution was
used, but removed all experimental restraint terms. Only "angl", "vdw",
"bond", "impr", and "rama" were left.
And I added or modified the following lines in proper positions of this
script for ensemble calculation:
########################################
from ensembleSimulation import EnsembleSimulation
ensembleSize=1
esim = EnsembleSimulation("ensemble",ensembleSize)
hitemp_potList = PotList()
xpTerms = ["BOND", "ANGL", "IMPR", "VDW", "RAMA"]
for xn in xpTerms :
hitemp_potList.add( AvePot(XplorPot, xn) )
pass
dyn = IVM(esim)
########################################
If the ensemble size equals to 1, the results are quite reasonable: I
calculated 500 times (numStructures=500) and analyzed their average radius
of gyration, their sizes are widely distributed with a averaged value close
to expected random coil size(in my case, 76 residues with radius of
gyration around 40 A).
However, if the ensemble size is set to 7, the output structures are much
compact, the radius of gyration is only around 26A. Something wrong with
the ensemble calling? Thanks in advanced!!
Best wishes,
Jie-rong
Here is the whole script I used for ensemble calculation.
########################
#
#
command = xplor.command
from jCoupPot import JCoupPot
from noePot import NOEPot
from xplor import select
from xplorPot import XplorPot
from rdcPotTools import *
from pdbTool import *
from atomAction import *
from selectTools import *
from simulationTools import *
from ivm import IVM
import string
import protocol
from ensembleSimulation import EnsembleSimulation
simWorld.setRandomSeed( 5764 )
command("""
structure @ubq1conf.psf end
param @./TopPar/parallhdg_new.pro
@./TopPar/parnah1er1_mod_new.inp
@./TopPar/par_axis_3.pro
@./TopPar/edtaThymine.par
@./TopPar/manganese.par
end
set echo = on message = on end
""")
protocol.initParams(('./TopPar/cysp.par'))
command("""
coor @random1conf1.pdb
coor copy end
""")
ensembleSize=7
esim = EnsembleSimulation("ensemble",ensembleSize)
tmps = open("./TopPar/proShortRange.list").read()
proShortRange = tmps.split(None)
#
#Rama torsion angle databasse
#
protocol.initRamaDatabase()
#
# setup parameters for nbond term.
#
command("""
parameters ! just repulsive term
nbonds
atom
nbxmod 3 wmin 0.01 cutnb 4.5 tolerance 0.5
rexp 2 irex 2 ctonnb=2.99 ctofnb=3.
repel= 0.5 rcon= 1.0 ! changed later
end
end
""")
#
# annealing setting
#
ini_rad = 0.9 ; fin_rad = 0.78
ini_con = 0.004 ; fin_con = 4.0
ini_ang = 0.4 ; fin_ang = 1.0
ini_imp = 0.1 ; fin_imp = 1.0
ini_rama = 0.002 ; fin_rama = 1.0
hot_cdi = 10.0 ; cool_cdi = 200.0
hitemp_potList = PotList()
xpTerms = ["BOND", "ANGL", "IMPR", "VDW", "RAMA"]
for xn in xpTerms :
hitemp_potList.add( AvePot(XplorPot, xn) )
pass
#set mass and fric
AtomSel("not (resname ANI or resname TAU)").apply(SetProperty("mass",
100.0))
AtomSel("resname ANI or resname TAU ").apply(SetProperty("mass", 300.0))
AtomSel("all").apply(SetProperty("fric", 10.0))
#
# IVM setup
#
dyn = IVM(esim)
dyn.reset() ; dyn.potList().removeAll()
#break ribose rings
IVM_breakRiboses(dyn)
#break proline rings
IVM_breakProlines(dyn)
IVM_groupRigidSidechain(dyn)
dyn.autoTorsion()
def setConstraints(k_ang,k_imp):
nconf = 1
scaleVdw = 1.0 / float( nconf )
command("""
constraints
inter = (segid "C1") (segid "C1")
weights * 1 angl %f impr %f vdw %f end
end""" % ( k_ang, k_imp, scaleVdw) )
return
def checkAndCorrect(dyn, pene, ene, tmstp, ni, maxEner):
# successful or not?
#
badIter = 0
if ene > maxEner :
badIter = 1
print "energy is too high (%7.1f)" % ene
if ene > pene and ene/pene > 10.0 :
badIter = 1
print "energy is higher than ten times of previous value (current:
%7.1f, previous: %7.1f)" % (ene, pene)
if tmstp < 1.0e-7 : #oritinal 1.0e-4
print "timestep is too short (%e ps)" % tmstp
badIter = 1
# for successful case
# - copy the coordinates to comparison set
#
if badIter == 0 :
command("""
coor copy end ! copy good coordinates
""")
return 1
# for unsuccessful case
# - copy back the previous coordinates
# - randomize initial velocities
command("""
coor swap end ! call previous good coordinates
coor copy end ! <- to erase bad one
""")
if ni > 3 :
print "Too many corrections. Exiting programs"
exit(1)
bath = dyn.bathTemp()
randomizeVelocities( bath )
dyn.resetReuse() #necessary to get randomized velocities
return 0
def dynamicsWithCheck(s, timestep, bath, finalTime, maxcycles):
s.calcEnergy()
ok = 0
ni = 0
while ok == 0 :
pene = s.Epotential()
s.setNumSteps(maxcycles)
s.setStepsize( timestep )
s.setBathTemp( bath )
s.setFinalTime( finalTime )
s.setPrintInterval( 500 ) #original 10
s.setETolerance( bath/1000 )
s.setResetCMInterval(10)
s.setResponseTime(20)
#dyn.setVerbose(0xffff)
dyn.setMaxDeltaE(10000)
#dyn.setAdjustStepsize(0)
s.run()
ene = s.Epotential()
tmstp = s.stepsize()
ok = checkAndCorrect(s, pene, ene, tmstp, ni, 100000)
ni = ni + 1
pass
return
def setProRamaForce(k_rama):
trace.suspend()
for cn in (proShortRange):
command("rama class %s force %f end " % (cn, k_rama))
pass
trace.resume()
print "(setProRamaForce) current k_rama: ", k_rama
return
def structLoopAction(loopInfo):
#
# read initial structure
#
command("""
coor init end
coor @random1conf1.pdb
coor copy end
""")
dyn.init()
#
# regularization of covalent geometry
#
command("""
flags exclude * include bond angl impr end
mini powell nstep 1000 nprint 50 end
""")
#
# high temp dynamics with vdw
#
# potential terms to use
dyn.setPotList( hitemp_potList)
setProRamaForce( ini_rama )
command("restraints dihed scale %f end" % hot_cdi )
# VDW setting for high temperature
command("""
parameters
nbonds
atom
nbxmod 3
wmin = 0.01 cutnb = 100 tolerance 45
repel= 1.2 rexp = 2 irex = 2 rcon=1.0
end
end
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
""")
# For PRE-calc in SB mode
# dynamics
timestep = 0.002
bath = 3000.0
dyn.setStepType("pc6")
randomizeVelocities( bath )
dyn.resetReuse() #necessary to get randomized velocities
dyn.calcEnergy()
for i_high in range(0,10) :
print " high temp dynamics -- subcycle No. %d" % i_high
dynamicsWithCheck(dyn, timestep, bath, 1.0, 500)
#print pre.showRestraints(0)
pass
#
# cooling
#
#number of cycles
numSteps = 240
global k_ang, k_imp
# VDW setting for cooling
command("""
parameters
nbonds
atom
nbxmod 3
wmin 0.01
cutnb 4.5
tolerance 0.5
rexp 2 irex 2
ctonnb=2.99 ctofnb=3.
repel=0.5 ! changed later
rcon=1.0 ! changed later
end
end
""")
setConstraints(ini_ang, ini_imp)
# setting force constants and etc.
k_rama = MultRamp(ini_rama,fin_rama,"setProRamaForce( VALUE )")
k_ang = MultRamp(ini_ang,fin_ang)
k_imp = MultRamp(ini_imp,fin_imp)
radius = MultRamp(ini_rad,fin_rad, "command('param nbonds repel VALUE
end end')")
k_vdw = MultRamp(ini_con,fin_con, "command('param nbonds rcon VALUE end
end')")
rampedParams = [k_rama,radius,k_vdw,k_ang,k_imp]
for k in rampedParams :
k.init(numSteps)
pass
command("restraints dihed scale %f end" % cool_cdi )
#dynamics
initTemp = 3000.0
finalTemp = 25.0 #original 25
coolTimestep = 0.002
tempStep = (initTemp-finalTemp)/float(numSteps)
finalTime = 1.0 #original 0.2
#maxCycle = 100 #<==original value
maxCycle = 500
bathTemp = initTemp;
for i_cool in range(0,numSteps):
bathTemp -= tempStep
for k in rampedParams :
k.update()
pass
setConstraints(k_ang.value(), k_imp.value())
print "cooling cycle - No.%d-" % i_cool
print "current bathTemp = %7.2f" % bathTemp
dynamicsWithCheck(dyn, coolTimestep, bathTemp, finalTime, maxCycle)
pass
#
# Powell minimization in torsion angle space
#
dyn.setStepType("powell")
dyn.setNumSteps(20000)
dyn.setMaxCalls(20000)
dyn.setETolerance(1.0e-7)
dyn.setGTolerance(0.01)
dyn.setDEpred(1.0)
dyn.run()
#
# Powell minimization in Cartesian space
#
xplor.simulation.potList().removeAll()
command("""
flags exclude *
!include bond angle impr vdw rama orie plan noe cdih coup carb
sani xdip coll scripting
include bond angle impr vdw rama
end
mini powell nstep 2000 nprint 10 end
""")
dyn.resetReuse()
dyn.init()
dyn.calcEnergy()
outPDB_template = "./RandomConfNOPRE/%d_MEMBER_%d.pdb" % (ensembleSize,
(loopInfo.count+1))
outFile = PDBTool( outPDB_template )
en = {}
eAll = dyn.Epotential()
rms = {}
violations = {}
for (term, threshold) in (("BOND",0.05),
("ANGL", 5.0),
("IMPR", 5.0),
("NOE", 0.5),
("CDIH", 5.0)
):
result = command("print thresh %f %s"% (threshold,term),
("RESULT","VIOLATIONS"))
rms[term] = float(result[0])
violations[term] = int(result[1])
pass
outFile.addRemark("OVERALL ENERGY: %7.2f " % eAll)
bars = "----------------------------------------------------------"
outFile.addRemark(bars)
outFile.addRemark("RMS")
outFile.addRemark("bond: %7.4f, angl: %7.4f, impr: %7.4f" %
(rms["BOND"], rms["ANGL"], rms["IMPR"]))
outFile.write()
pass
StructureLoop(numStructures=50,startStructure=0 ,
structLoopAction=structLoopAction).run()
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
http://dcb.cit.nih.gov/pipermail/xplor-nih/attachments/20081106/ac126daa/attachment-0001.html
