Hi Charles,
I am finding that I do not know how to relax the rigid body, could you make
some recommendations? The read in structure is refined using the dna_refi
refine.py including the RDCs of the stems so the starting structure going
into the rigid body calculation is pretty good but lacks the global
alignment and consensus fold. My last thought, at the moment I am using
two refinement scripts, but I wonder should I combine the refine and rigid
body minimization into a single script?
On Fri, May 22, 2015 at 7:22 AM, Charles Schwieters <[email protected]>
wrote:
>
> Hello--
>
> >
> > could you let me know if the script I included is a good starting point?
> > Removing the SAXS restraints and including the RDCs?
> >
>
> Aside from not including RDCs it *seems* to be fine.
>
> Charles
>
# writen by Jinbu Wang
#
xplor.requireVersion("2.19")
init_t = 1500
numberOfStructures=1 #was 100
ensembleSize = 3
fn_pdb_init = "../RDC_Refine27/REF/rdc_refine.sa"
#xray_scale = 400
#by default, all ensemble members have equal weights, but you can change this:
#
ensWeights=[1.0/ensembleSize for i in range(ensembleSize)]
startStructure=0
outFilename = "SCRIPT_%d_STRUCTURE_MEMBER.pdb" % numberOfStructures
### RRE Restraints
HBOND = True
NOE = True
DIHEDRAL = True
ORIENT = True
#fn_pre = "rre"
#fn_tbl_dih = "%s_dih.tbl" % fn_pre
fn_tbl_dih = "../cons/FL_TOR7.inp"
#fn_tbl_plan = "%s_plan.tbl" % fn_pre
fn_tbl_plan = "@../cons/FL_PLANE2.inp"
#fn_tbl_hpp = "%s_hpp.tbl" % fn_pre
#fn_tbl_dpp = "%s_dpp.tbl" % fn_pre
#fn_tbl_hb = "%s_hb.tbl" % fn_pre
fn_tbl_hb = "../cons/FL_HBOND4.inp"
#fn_tbl_pp = "%s_pp.tbl" % fn_pre
#fn_tbl_ext = "%s_ext.tbl" % fn_pre
#fn_tbl_noe = "%s_noe.tbl" % fn_pre
fn_tbl_noe = "../cons/FL_RES27.inp"
noe_rsts = [
#("DPP", fn_tbl_dpp, "hard", 50),
#("HPP", fn_tbl_hpp, "hard", 50),
("NHB", fn_tbl_hb, "hard", 50),
#("PP", fn_tbl_pp, "hard", 50),
#("EXT", fn_tbl_ext, "hard", 50),
("NOE", fn_tbl_noe, "hard", 50),
]
rigidRegions=[
'resid 1:14 or resid 55:68',
'resid 21:49',
##'resid 7:15 or resid 219:225',
##'resid 16:19 or resid 212:215',
##'resid 20:23 or resid 207:210',
##'resid 28:31 or resid 200:202',
##'resid 32:36 or resid 194:197',
##'resid 38:44 or resid 97:103',
##'resid 48:49 or resid 68:69',
##'resid 50:56 or resid 60:66 or resid 78:83 or resid 89:94',
##'resid 105:111 or resid 148:154',
##'resid 119:123 or resid 142:146',
##'resid 127:129 or resid 135:137',
##'resid 163:171 or resid 180:188',
]
##breakResids=[226, 6, 218, 15, 211, 19, 206, 27, 199, 31, 193, 37, 96,
## 47, 67, 49, 59,77, 88, 104, 147, 118, 141, 126, 134, 162, 179]
breakResids=[15, 16, 17, 18, 19, 20, 37, 38, 39, 40]
import protocol
protocol.initRandomSeed()
protocol.topology['nucleic'] = "nucleic-2.0.top"
protocol.parameters['nucleic'] = "nucleic-2.0.par"
protocol.initParams("nucleic")
# read in the PSF and initial PDB files
#
import psfGen
protocol.loadPDB(fn_pdb_init)
xplor.simulation.deleteAtoms("not known")
protocol.fixupCovalentGeom(maxIters=100,useVDW=1) #added 20150523
from ensembleSimulation import EnsembleSimulation
esim=EnsembleSimulation('ens',ensembleSize)
esim.setAveType('sum')
# list of potential terms used in refinement
from potList import PotList
potList = PotList()
crossTerms=PotList('cross terms')
# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams
rampedParams=[]
#from xplorPot import XplorPot
#from avePot import AvePot
#from solnXRayPotTools import create_solnXRayPot
#import solnXRayPotTools
#
#fn_saxsData="rre_nsw23c.dat"
#xray=create_solnXRayPot('xray',experiment=fn_saxsData,
# normalizeIndex=-3,preweighted=False)
#
#xrayCorrect=create_solnXRayPot('xray-c',experiment=fn_saxsData,
# normalizeIndex=-3,preweighted=False)
#
#solnXRayPotTools.useGlobs(xray)
#
#xray.setNumAngles(50)
#xrayCorrect.setNumAngles(500)
#xray.setScale(xray_scale)
#xray.setCmpType("plain")
#xray.setEnsWeights( ensWeights )
#xrayCorrect.setEnsWeights( ensWeights )
#
#potList.append(xray)
#crossTerms.append(xrayCorrect)
#
#print xray.calcEnergy()
#from solnScatPotTools import fitParams
#correct I(q) to higher accuracy, and include solvent contribution corrections
#stride=10 specifies that this fit is performed every 100th temperature during
#simulated annealing. During refinement, infrequent solvent correction updates
#are ok because the structure doesn't change too much.
#
#rampedParams.append(
# StaticRamp(
# "fitParams(xrayCorrect);xray.calcGlobCorrect(xrayCorrect.calcd())" ))
from varTensorTools import create_VarTensor
media={}
# medium Da rhombicity
for (medium,Da,Rh) in [ #('ph', -50.00, 0.093),
#('a', -38.32, 0.154),
('b', -50.00, 0.011)
#('c', -27.36, 0.50)
]:
oTensor = create_VarTensor(medium)
oTensor.setDa(Da)
oTensor.setRh(Rh)
media[medium] = oTensor
pass
from rdcPotTools import create_RDCPot, scale_toNH
rdcs = PotList('rdc')
for (medium,expt,file, scale) in \
[#('ph','NH_3' ,'../cons/Mod3_RDCNH1.inp' ,1),
#('ph','CH_3','../cons/Mod3_RDCCH1.inp' ,1),
#('a','NH_4','../cons/Mod4_RDCNH2.inp' ,1),
#('a','CH_4','../cons/Mod4_RDCCH2.inp' ,1),
('b','NH_2','../cons/FL_RDCNH3.inp' ,1)
# ('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.01,1.0, "rdcs.setScale( VALUE )") ) #was 0.05,5.0 2Apr
# calc. initial tensor orientation
# and setup tensor calculation during simulated annealing
#
from varTensorTools import calcTensorOrientation, calcTensor
for medium in media.keys():
calcTensorOrientation(media[medium])
rampedParams.append( StaticRamp("calcTensor(media['%s'])" % medium) )
pass
from avePot import AvePot
from xplorPot import XplorPot
rampedParams.append( StaticRamp("""protocol.initNBond()""") )
rampedParams.append( MultRamp(0.9,0.8,
"xplor.command('param nbonds repel VALUE end end')") )
rampedParams.append( MultRamp(.004,4,
"xplor.command('param nbonds rcon VALUE end end')") )
#highTempParams.append( StaticRamp("""protocol.initNBond(rcon=0.004,
# repel=1.2,
# onlyCA=1)""") )
#
nbond=AvePot(XplorPot,"VDW")
potList.append(nbond)
protocol.initRamaDatabase('nucleic')
rama= AvePot(XplorPot,'RAMA')
potList.append(rama)
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
### new POT
if fn_tbl_dih:
protocol.initDihedrals(fn_tbl_dih,scale=50)
potList.append(AvePot(XplorPot, "cdih"))
###planarity restraints
if fn_tbl_plan:
xplor.command("""
restraints plane
@%s
end
""" % fn_tbl_plan)
potList.append(AvePot(XplorPot,"plan"))
if ORIENT:
#initialize the aa-aa positional database
xplor.command("SET ECHO=OFF END")
xplor.command("SET messages NONE END")
xplor.command("@../cons/FL.setup")
potList.append(AvePot(XplorPot,"orie"))
rampedParams.append( StaticRamp("potList['ORIE'].setScale(0.2)") )
if HBOND:
potList.append(AvePot(XplorPot,'HBON'))
potList['HBON'].setScale(5)
if NOE:
noePots = PotList("noe")
from noePotTools import create_NOEPot
for pot_name, fn_tbl, pot_type, pot_scale in noe_rsts:
if fn_tbl is None: continue
noe = create_NOEPot(pot_name, fn_tbl)
noe.setPotType(pot_type)
noe.setScale(pot_scale)
noePots.append(noe)
potList.append(noePots)
###
for name in ("bond","angl","impr"):
potList.append( AvePot(XplorPot,name) )
pass
rampedParams.append( MultRamp(0.4,1.0,"potList['ANGL'].setScale(VALUE)"))
rampedParams.append( MultRamp(0.1,1.0,"potList['IMPR'].setScale(VALUE)"))
# Give atoms uniform weights, except for the anisotropy axis
#
protocol.massSetup()
from ivm import IVM
dyn = IVM()
protocol.initDynamics(dyn,potList=potList)
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
for region in rigidRegions:
sel = "(%s)" % region
dyn.group(sel)
dyn.hinge('full',sel)
for resid in breakResids:
dyn.breakBond("(resid %d and name O3') or (resid %d and name P)"%(resid,
resid+1))
protocol.torsionTopology(dyn)
##
rigidNucleicSelections = {}
for (res,sel) in \
(("CYT","""name n1 or name c6 or name c5 or name c4 or name n3 or name c2 OR
NAME N4 or name h5 or name h6 or name o2"""),
("GUA","""name n9 or name c4 or name n3 or name c2 or name n1 or name n2 or
name c6 or name c5 or name n7 or name c8 or name o6
or name h1 or name h8"""),
("ADE","""name n9 or name c4 or name n3 or name c2 or name n1 or
name c6 or name c5 or name n7 or name c8 or name n6
or name h2 or name h8"""),
#add n6?
("TED","name n1 or name c6 or name c5 or name c4 or name n3 or name c2"),
("THY","""name n1 or name c6 or name c5 or name c4 or name cm
or name n3 or name c2 or name o4 or name o2 or name h3 or name h6"""),
("URI","""name n1 or name c6 or name c5 or name c4
or name n3 or name c2 or name o4 or name o2 or name h3 or name h5
or name h6""")):
rigidNucleicSelections[res] = sel
minc = IVM()
for region in rigidRegions:
sel = "(%s)" % region
minc.group(sel)
minc.hinge('full',sel)
for resid in breakResids:
minc.breakBond("(resid %d and name O3') or (resid %d and name P)"%(resid,
resid+1))
protocol.cartesianTopology(minc)
protocol.initMinimize(minc,potList=potList)
annealTime=0.2 #time to integrate at a given annealing temp.
annealSteps=0 #max num steps to be taken during a single annealing temp.
from simulationTools import AnnealIVM
anneal= AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =25,
ivm=dyn,
rampedParams = rampedParams)
from simulationTools import testGradient
from pdbTool import PDBTool
def calcOneStructure(loopInfo):
from atomAction import randomizeDomainPos
# initialize parameters for high temp dynamics.
InitialParams( rampedParams )
# high temperature bit
protocol.initDynamics(dyn,
initVelocities=1,
bathTemp=init_t,
potList=potList,
numSteps=5000,
finalTime=10)
dyn.run()
protocol.initNBond() #reset to include all atoms
# 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
#
anneal.run()
protocol.initMinimize(dyn)
dyn.run()
##
##all atom minimization
##
minc.run()
#
# perform analysis and write structure
loopInfo.writeStructure(potList,crossTerms)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
startStructure=startStructure,
structLoopAction=calcOneStructure,
pdbTemplate=outFilename,
genViolationStats=1,
averageTopFraction=0.5,
averagePotList=potList,
averageCrossTerms=crossTerms,
averageContext=FinalParams(rampedParams),
#averageFilename="ave.pdb",
averageFitSel="not hydro and not resname ANI",
averageCompSel="not hydro and not resname ANI" ).run()
_______________________________________________
Xplor-nih mailing list
[email protected]
https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
