| hello Charles, I want to do the rigid body conformational sampling using the dock2.py in the attachment. but, 1) I found that all of the 5 domains do not rotate randomly and freely, instead just rotate or translate around the center of mass. Then the sampling may not search the conformational space enough. 2) I do not know how to edit the scripts to let the 5 domains move randomly and freely. Do they can be effected with “randomizeDomainPos”? |
#dock EIN to EIC # repeat numDock times: # randomize domain positions # enforce symmetry # satisfy rdcs # dock using xray term (possibly including Rgyr to collapse) - # allowing only translation # # keep coordinates from best run
# then: attempt to make broken bond by moving only the linker region
# if we can't, print fail and go to next structure
# if we can, do refinement
import os
print 'process id:', os.getpid()
numberOfStructures=120
import protocol
protocol.initRandomSeed(42)
ensembleSize=1
cutoffRfactor=18.22 # need an R-factor better than this to accept a docked struct
numDock=2
pdbTemplate="SCRIPT_STRUCTURE_MEMBER.fit"
protocol.loadPDB("wl-14p-1.pdb",deleteUnknownAtoms=True)
groupedResids = [(1,71), ( 85,152), ( 167,214), ( 221,305), ( 315,412) ]
breakAtoms = ["C5'", "O5'"]
linkerResids = sum( [[( groupedResids[i][1],groupedResids[i+1][0])] for
i in range(len(groupedResids)-1)] , [] )
breakResids = [int(0.5*(pair[0]+pair[1])) for pair in linkerResids]
linkerSel = " or ".join(["resid %d:%d" % pair for pair in linkerResids])
import regularize
regularize.fixupCovalentGeomIVM("not resid %d:%d" % groupedResids[2],
rigidRegions=["resid %d:%d" % pair for
pair in groupedResids[:2]+
groupedResids[3:]],
)
from ensembleSimulation import EnsembleSimulation
esim=EnsembleSimulation('ens',ensembleSize)
from potList import PotList
potList = PotList()
crossTerms=PotList('cross terms')
crossTerms=PotList('cross terms')
from simulationTools import StaticRamp, MultRamp
from simulationTools import InitialParams, FinalParams
rampedParams=[]
highTempParams=[]
media=[]
from varTensorTools import create_VarTensor, addAxisAtoms
from varTensorTools import calcTensor
axisAtoms=addAxisAtoms(resid=999)
ten1=create_VarTensor("medium1",axis=axisAtoms)
media.append(ten1)
ten1.setFreedom("varyDa, varyRh")
#ten1.setFreedom("fixDa, fixRh")
#rdcs=PotList("rdc")
#from rdcPotTools import create_RDCPot
#rdc1 = create_RDCPot("rdc1","free_Mg_A.tbl",ten1)
#calcTensor(ten1) #this to set Da, Rh
#rdc1.addRestraints(open("free_Mg_B.tbl").read())
#rdc1.setShowAllRestraints(1)
#rdcs.append(rdc1)
#potList.append(rdcs)
#rampedParams.append( StaticRamp('for t in media: calcTensor(t)') )
from avePot import AvePot
from xplorPot import XplorPot
#
#from xplorPot import XplorPot
#protocol.initDihedrals("torsion_linker.tbl",
# useDefaults=False)
#dihed=AvePot(XplorPot,"CDIH")
#potList.append(dihed)
#highTempParams.append( StaticRamp("potList['CDIH'].setScale(0.1)") )
#rampedParams.append( MultRamp(0.1,200,"potList['CDIH'].setScale(VALUE)") )
#
selPairs=[(linkerSel,"not PSEUdo")]
for i,pair1 in enumerate(groupedResids):
for pair2 in groupedResids[i+1:]:
selPairs.append( ("resid %d:%d"%pair1,"resid %d:%d"%pair2) )
pass
pass
from repelPotTools import create_RepelPot,initRepel
repel = AvePot(create_RepelPot('repel',
selPairs=selPairs,
selection=AtomSel('not PSEUDO',esim.member())))
potList.append(repel)
#
#from solnXRayPotTools import create_solnXRayPot
#import solnXRayPotTools
#xray=create_solnXRayPot('xray',experiment='EI-Mg-saxswaxs-trunc.dat',
# numPoints=26,
# normalizeIndex=-3,preweighted=False)
#
#xrayCorrect=create_solnXRayPot('xray-c',experiment='EI-Mg-saxswaxs-trunc.dat',
# numPoints=26,
# normalizeIndex=-3,preweighted=False)
#
#solnXRayPotTools.useGlobs(xray)
#
#xray.addRigidRegion(cTerminus)
#xray.addRigidRegion("segid A and (%s)" % nTerminus)
#xray.addRigidRegion("segid B and (%s)" % nTerminus)
#
#
#xray.setNumAngles(50)
#xrayCorrect.setNumAngles(500)
#xray.setScale(400)
#xray.setCmpType("plain")
#potList.append(xray)
#crossTerms.append(xrayCorrect)
#
#print xray.calcEnergy()
#from solnScatPotTools import fitParams
##fitParams(xrayCorrect)
#
##xray.calcGlobCorrect(xrayCorrect.calcd())
##corrects I(q) to higher accuracy, and include solvent contribution corrections
#rampedParams.append( StaticRamp("fitParams(xrayCorrect)") )
#rampedParams.append( StaticRamp("xray.calcGlobCorrect(xrayCorrect.calcd())") )
#
#from residueAffPotTools import create_ResidueAffPot
#contact = create_ResidueAffPot("contact")
#crossTerms.append(contact)
#
#protocol.initRamaDatabase()
#rama= AvePot(XplorPot,'RAMA')
#potList.append(rama)
#rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
# Statistical torsion angle potential
#
from torsionDBPotTools import create_TorsionDBPot
torsionDBPot = create_TorsionDBPot('tDB',system="rna")
potList.append( torsionDBPot )
rampedParams.append( MultRamp(.002,2,"torsionDBPot.setScale(VALUE)") )
bond= AvePot(XplorPot,"BOND")
refBondViolations=bond.violations() # don't wan't more violations than this
print 'refBondViolations:', refBondViolations
print 'initial bond energy:', bond.calcEnergy()
potList.append( bond)
angl= AvePot(XplorPot,"ANGL")
potList.append( angl)
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )
impr= AvePot(XplorPot,"IMPR")
potList.append( impr )
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )
protocol.massSetup()
from ivm import IVM
dynRigid = IVM()
#dynRigid.fix("resid %d:%d" % groupedResids[2])
for pair in groupedResids:
dynRigid.group("resid %d:%d" % pair)
pass
for resid in breakResids:
dynRigid.breakBond("resid %d and (name %s or name %s)" %
(resid,breakAtoms[0],breakAtoms[1]))
pass
protocol.cartesianTopology(dynRigid)
#ten1.setFreedom("fix")
dynTranslate = IVM() # this one fixes the orientation of the a and ab domains
#dynTranslate.fix("resid %d:%d" % groupedResids[2])
for pair in groupedResids:
dynTranslate.hinge("translate","resid %d:%d" % pair)
dynTranslate.group("resid %d:%d" % pair)
pass
for resid in breakResids:
dynTranslate.breakBond("resid %d and (name %s or name %s)" %
(resid,breakAtoms[0],breakAtoms[1]))
pass
protocol.cartesianTopology(dynTranslate)
from simulationTools import AnnealIVM
init_t = 3000
#cool = AnnealIVM(initTemp =init_t,
# finalTemp=25,
# tempStep =25,
# ivm=dynTranslate,
# 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.
"""
minEnergy = 1e30 # big number
startPos = xplor.simulation.atomPosArr()
InitialParams(rampedParams)
k=0
while k < numDock:
xplor.simulation.setAtomPosArr(startPos)
from atomAction import randomizeDomainPos
for pair in groupedResids:
randomizeDomainPos( "resid %d:%d" % pair,
deltaPos=80 )
pass
# generate a new structure with randomized torsion angles
#
from monteCarlo import randomizeTorsions
randomizeTorsions(dynRigid)
# #fit RDCs
# protocol.initMinimize(dynRigid,
# potList=[dSymm,ncs,rdcs],
# numSteps=1000)
#
# dynRigid.run() ; calcTensor(ten1); dynRigid.run(); calcTensor(ten1)
#fix broken bond
# if bond is excluded, linker2 can extend a long distance
# to make an unphysical connection
protocol.initMinimize(dynRigid,
potList=[bond,],
numSteps=1000)
# bond.setScale(1e-5) ; dynRigid.run() ;
# calcTensor(ten1)
# bond.setScale(1e-3) ; dynRigid.run() ;
# calcTensor(ten1)
bond.setScale(1) ; dynRigid.run() ;
# from simulationTools import analyze
# print analyze(rdcs)
# from rdcPotTools import Rfactor_infinite
# rFac = Rfactor_infinite(rdcs[0])
# print 'rFactor: %.2f' % rFac
# if rFac>cutoffRfactor:
# print 'bad R factor. Discarding this configuration'
# continue
# #dock using Rgyr, saxs, vdw, contact?
# protocol.initMinimize(dynTranslate,
# potList=[bond,dSymm,ncs,rdcs,Rgyr],
# numSteps=1000)
#
## xray.setVerbose(2)
## dynTranslate.setVerbose( dynTranslate.printNodeDef )
# dynTranslate.run() ;
# dynTranslate.setVerbose( 0 )
#
# InitialParams(rampedParams) #done for xray-correct fitting
# protocol.initMinimize(dynTranslate,
# potList=[bond,dSymm,ncs,rdcs,Rgyr,xray],
# numSteps=1000)
# dynTranslate.run()
#
InitialParams(rampedParams) #done for xray-correct fitting
#attempt to fix linker geom
try:
protocol.fixupCovalentGeom(sel=linkerSel)
except protocol.CovalentViolation:
pass
print 'number of bond violations:', bond.violations()
if bond.violations()>refBondViolations:
print 'too many bond violations. Discarding this configuration'
continue
# try fixing everything together
protocol.initMinimize(dynRigid,
potList=[bond,angl,impr,repel],
numSteps=1000)
dynRigid.run() ; dynRigid.run()
# InitialParams(rampedParams) # refit xray-correct before scoring
scorePot = PotList()
for p in (bond,repel):
scorePot.append(p)
pass
scoreEnergy = scorePot.calcEnergy()
print 'dock iteration: %d score energy: %.2f' %(k,scoreEnergy)
if scoreEnergy < minEnergy :
minPos = xplor.simulation.atomPosArr()
minEnergy = scoreEnergy
pass
k += 1
pass
xplor.simulation.setAtomPosArr(minPos)
InitialParams( rampedParams )
InitialParams( highTempParams )
#high-temp dynamics
protocol.initDynamics(dynRigid,
potList=potList,
bathTemp=init_t,
initVelocities=1,
finalTime=80, # stops at 80ps or 800 steps
numSteps=800, # whichever comes first
printInterval=100)
dynRigid.setETolerance( init_t/100 ) #used to det. stepsize. default: t/1000
dynRigid.run()
# # initialize parameters for cooling loop
# InitialParams( rampedParams )
#
#
# # initialize integrator for simulated annealing
# #
# protocol.initDynamics(dynTranslate,
# potList=potList,
# numSteps=100, #at each temp: 100 steps or
# finalTime=.2 , # .2ps, whichever is less
# printInterval=100)
#
# # perform simulated annealing
# #
# cool.run()
#
# # refit SAXS for final minimization
# FinalParams( rampedParams )
# final torsion angle minimization
#
protocol.initMinimize(dynRigid)
dynRigid.run()
pass
from simulationTools import StructureLoop
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=pdbTemplate,
doWriteStructures=True,
structLoopAction=calcOneStructure,
genViolationStats=True,
averageRegularize=False,
averagePotList=potList,
averageCrossTerms=crossTerms,
averageFitSel="resid %d:%d" % groupedResids[2],
averageTopFraction=0.5, #report on all structs
averageContext=FinalParams(rampedParams)).run()
Yikan Ph. D candidate TsingHua University School of life Sciences Beijing, China |
_______________________________________________ Xplor-nih mailing list [email protected] https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
