Hello all, I am doing the ensemble simulation against SAXS data, and found that when I set different values of xray_scale, the results are very differnt. when I set my scale as 4000, I get the results as: Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 5867317.88(36155.84) 15.416( 0.034) 7344.7(135.8) 19621 SelNBPot contact 0.00( 0.00) SolnScatPot xray 3819.14(2231.67) 0.549( 0.147) XplorPot ANGL 383549.55(12825.12) 9.295( 0.157) 4057.5( 65.7) 9713 XplorPot BOND 104110.35(6190.68) 0.079( 0.002) 826.1( 23.0) 5499 XplorPot IMPR 5392936.62(25468.00) 51.739( 0.122) 2380.9( 41.9) 4409 XplorPot RAMA -21738.24(429.25) XplorPot VDW 4640.46(1444.40) 80.1( 11.7) Cross Validated Terms Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 0.95( 0.56) 0.549( 0.147) 0.0( 0.0) SolnScatPot xray-c 0.95( 0.56) 0.549( 0.147) Statistics for Additional Quantities SolnScatPot Chi^2 xray 0.320( 0.187) But when I set the scale as 400 the results are: Results for the top 5 (of 10) structures Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 5771888.93(25598.60) 15.565( 0.072) 7538.3( 37.4) 19621 SelNBPot contact 0.00( 0.00) SolnScatPot xray 3552.58(430.47) 1.718( 0.106) XplorPot ANGL 372003.14(12908.97) 9.154( 0.158) 4162.5( 16.3) 9713 XplorPot BOND 108314.08(7276.59) 0.081( 0.003) 875.9( 13.0) 5499 XplorPot IMPR 5302939.72(30393.59) 51.305( 0.147) 2411.7( 13.9) 4409 XplorPot RAMA -20293.64(361.01) XplorPot VDW 5373.05(1789.52) 88.2( 19.1) Cross Validated Terms Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 8.88( 1.08) 1.718( 0.106) 0.0( 0.0) SolnScatPot xray-c 8.88( 1.08) 1.718( 0.106) Statistics for Additional Quantities SolnScatPot Chi^2 xray 2.981( 0.361) when the xray_scale=40 Results for the top 5 (of 10) structures Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 5874492.26(51971.20) 15.878( 0.092) 7349.1(159.0) 19621 SelNBPot contact 0.00( 0.00) SolnScatPot xray 656.99( 21.02) 2.340( 0.037) XplorPot ANGL 384002.99(12819.76) 9.301( 0.155) 4063.4( 93.4) 9713 XplorPot BOND 101573.48(8302.84) 0.078( 0.003) 821.1( 27.8) 5499 XplorPot IMPR 5404253.65(51091.47) 51.793( 0.245) 2383.0( 24.3) 4409 XplorPot RAMA -21204.30(1217.53) XplorPot VDW 5209.46(2849.19) 81.6( 31.4) Cross Validated Terms Average values for potential terms: type name Energy(dev) RMSD(dev) viols(dev) Num PotList TOTAL 16.42( 0.53) 2.340( 0.037) 0.0( 0.0) SolnScatPot xray-c 16.42( 0.53) 2.340( 0.037) Statistics for Additional Quantities SolnScatPot Chi^2 xray 5.512( 0.176) 1) And I got three different chi^2. But I do not know which one is better? or none of them values are not so good? 2) what is the meaning of “dev” in Energy(dev)? Did the value of the “xray energy” term above(3819.14,3552.58,656.99 ) had been scaled by the “xray_scale” in the script? and how? 3) did "nonbonded and torsionDB in this case” means VDW and RAMA? to get the reasonable or good results, how can i set the value of “xray_sacle”(the xray energy term) to match “nonbonded and torsionDB” terms? 4) Last but most important: when i checked the output pdbs, I found the output does not change much compared to my initial structure. But the chi^2 of the xplor-nih is very small in the output file.(but when i caculate the chi^2 by Crysol the chi^2 is very big). I do not know why. 5) I have attached my scripts bellow. My SAXS data has been extrapolate and qmax=0.3. Yikan |
xplor.requireVersion("2.19")
import protocol
protocol.initRandomSeed()
protocol.topology['nucleic']="nucleic-2.0.top"
protocol.parameters['nucleic']="nucleic-2.0.par"
protocol.initParams("nucleic")
xray_scale = 400
numberOfStructures=10
fn_pdb = "Ddb12-2np-qr-noh.pdb"
fn_saxsData='YP1_ddb12.dat'
ensembleSize = 1
#by default, all ensemble members have equal weights, but you can change this:
#
ensWeights=[1.0/ensembleSize for i in range(ensembleSize)]
####ensWeights=[0.8,0.1,0.1]
startStructure=0
outFilename = "SCRIPT_%d_STRUCTURE_MEMBER.pdb" % numberOfStructures
rigidRegions=[
'resid 1:82', 'resid 89:172',
]
breakResids=[83,88]
#
# read in the PSF and initial PDB files
#
import psfGen
####psfGen.pdbToPSF("fix_HAB.pdb")
#
# starting coords
#
#protocol.initCoords("tectoRNA.pdb")
#protocol.loadPDB("fsr0_s0_102.pdb")
#start from a good structure...
protocol.loadPDB(fn_pdb)
#protocol.loadPDB("ref-11-2-2011b_100_0-relocated-d783.pdb")
#protocol.writePDB("/usr/home/schwitrs/t0.pdb")
xplor.simulation.deleteAtoms("not known")
#import regularize
#regularize.fixupCovalentGeomIVM(translateRegions=rigidRegions)
#protocol.covalentMinimize("""not (resid 1:76 or resid 117:156 or resid 173:283 or resid 319:388 or resid 410:617 or resid 769:834 or resid 634:684 or resid 699:750)""")
#protocol.writePDB("/usr/home/schwitrs/t1.pdb")
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
#protocol.initCollapse("resid 1:835",
# Rtarget=38.2)
#potList.append( XplorPot('COLL') )
#SAXS term
from solnXRayPotTools import create_solnXRayPot
import solnXRayPotTools
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(50)
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())",
stride=100))
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)
# t/a db
protocol.initRamaDatabase('nucleic')
rama= AvePot(XplorPot,'RAMA')
potList.append(rama)
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
# database contact term
from residueAffPotTools import create_ResidueAffPot
contact = AvePot(create_ResidueAffPot("contact"))
potList.append(contact)
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 region in rigidRegions:
sel = "(%s) and (name O3' or name P or name O5')" % region
sel = "(%s)" % region
dyn.group(sel)
#dyn.hinge('translate',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) and (name O3' or name P or name O5')" % region
sel = "(%s)" % region
minc.group(sel)
#minc.hinge('translate',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)
init_t = 3000
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=DynMin(dyn),
ivm=dyn,
rampedParams = rampedParams)
from simulationTools import testGradient
#testGradient(potList,eachTerm=1)
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.writePDB("/usr/home/schwitrs/t2.pdb")
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.writePDB("/usr/home/schwitrs/t3.pdb")
#
# torsion angle minimization
#
protocol.initMinimize(dyn)
dyn.run()
# protocol.writePDB("/usr/home/schwitrs/t4.pdb")
##
##all atom minimization
##
minc.run()
# protocol.writePDB("/usr/home/schwitrs/t5.pdb")
#### dyn.group('resid 1:76 or resid 117:156')
#### dyn.group('resid 173:283')
#### dyn.group('resid 319:388 or resid 410:617 or resid 769:834')
#### dyn.group('resid 634:684')
#### dyn.group('resid 699:750')
#### minc.hinge('full','resid 1:76 or resid 117:156')
#### minc.hinge('full','resid 173:283')
#### minc.hinge('full','resid 319:388 or resid 410:617 or resid 769:834')
#### minc.hinge('full','resid 634:684')
#### minc.hinge('full','resid 699:750')
#
# 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
