Hi Charles,
Thanks for the reply, but even when commenting out the rigid regions angl,
impr, and cdih blowup. So I went back to using the the dna_refi as a
template for subsequent refinement. I think that is it working ok, in the
attached script I am reading in a structure calculated with the RDCs from
the upper and lower stems and applying only the iminos from the entire RNA
but I do not understand why ANI ends up in the middle of the molecule. And
when I include multiple sets of RDC data sets, why the ANI for each data
set are overlapped when the tensors differ (should they not be distinct?).
Is there a way to move ANI farther from the molecule? Or is it something in
the script? Or is it the version of xplor (I am using 2.33)?
On Tue, May 26, 2015 at 7:05 AM, Charles Schwieters <[email protected]>
wrote:
>
> Hello--
>
> >
> > 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?
> >
>
> I probably am not understanding, but is it sufficient to simply
> comment-out all of your rigidRegions? I must be missing something.
>
> Charles
>
xplor.requireVersion("2.24")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2005/05/10
#
(opts,args) = xplor.parseArguments(["slow"]) # check for command-line typos
quick=True
for opt in opts:
if opt[0]=="slow": #specify -slow for full calc
quick=False
pass
pass
outFilename = "./REF2/SCRIPT_STRUCTURE.sa"
numberOfStructures=200
if quick:
numberOfStructures=4
pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#
# annealing settings
#
command = xplor.command
protocol.initParams("nucleic")
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('RNA1.seq')
#protocol.initStruct("g_new.psf") # - or from file
# 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())
# or read an existing model
#
protocol.loadPDB("../../RDC2_Refine/REF1/rdc_refine_16.sa")
xplor.simulation.deleteAtoms("not known")
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
#
# 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*")
#} {rm 2Apr}
# 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 [ #('ph', -50.00, 0.093),
#('a', -38.32, 0.154),
('b', -35.7, 0.083)
#('c', -27.36, 0.50)
]:
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_toCH
rdcs = PotList('rdc')
for (medium,expt,file, scale) in \
[#('ph','NH_3' ,'../../cons/SL2RDCNH1.inp' ,1),
#('ph','CH_3','../../cons/SL2RDCCH1.inp' ,1),
#('a','NH_4','../../cons/SL1RDCNH2.inp' ,1),
#('a','CH_4','../../cons/SL1RDCCH2.inp' ,1),
('b','NH_2','../../cons/FLRDCNH3.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(3) # in Hz
rdcs.append(rdc)
pass
potList.append(rdcs)
rampedParams.append( MultRamp(0.005,0.25, "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
#planarity restraints
xplor.command("@../../cons/FLPLANE2.inp")
potList.append(AvePot(XplorPot("plan",xplor.simulation)) )
#initialize the aa-aa positional database
xplor.command("@../../cons/FL.setup")
potList.append(XplorPot("orie"))
rampedParams.append( StaticRamp("potList['ORIE'].setScale(0.2)") )
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"../../cons/FLRES27.inp"),
('hb',1,"../../cons/FLHBOND4.inp")
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
#pot.setPotType("soft") #soft - if you think there may be bad NOEs
pot.setScale(scale)
pot.setThreshold(0.4)
noe.append(pot)
#rampedParams.append( MultRamp(2,50, "noe.setScale( VALUE )") ) #was 50 2Apr
rampedParams.append( StaticRamp("noe.setScale( 50 )") )
# 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)
# {rm 2Apr}
# Set up dihedral angles
from xplorPot import XplorPot
protocol.initDihedrals("../../cons/FLTOR7.inp",
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 ) #5 degrees is the default value, though
# gyration volume term
#
# gyration volume term
#
#{from gyrPotTools import create_GyrPot
#gyr = create_GyrPot("Vgyr",
# "resid 1:56") # selection should exclude disordered tails
#potList.append(gyr)
#rampedParams.append( MultRamp(.002,1,"gyr.setScale(VALUE)") )
#}{rm 2Apr}
# hbda - distance/angle bb hbond term
#
#{protocol.initHBDA('hbda.tbl')
#potList.append( XplorPot('HBDA') )
#}{rm 2Apr}
#Rama torsion angle database
#
protocol.initRamaDatabase('nucleic')
potList.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].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(0.9,0.8,
"command('param nbonds repel VALUE end end')") )
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,
# selStr="nameP")""") ) #included ,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()
# initially minimize in Cartesian space with only the covalent constraints.
# Note that bonds, angles and many impropers can't change with the
# internal torsion-angle dynamics
# breaks bonds topologically - doesn't change force field
#
#dyn.potList().add( XplorPot("BOND") )
#dyn.potList().add( XplorPot("ANGL") )
#dyn.potList().add( XplorPot("IMPR") )
#
#dyn.breakAllBondsIn("not resname ANI")
#import varTensorTools
#for m in media.values():
# m.setFreedom("fix") #fix tensor parameters
# varTensorTools.topologySetup(dyn,m) #setup tensor topology
#
#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
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
protocol.torsionTopology(dyn)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
#for m in 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 = 700. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=25, # {was 12}
tempStep =25, #{was 12.5}
ivm=dyn,
rampedParams = rampedParams)
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>0:
return False
if potList['rdc'].rms()>2.0: #was 1.2 23Apr this might be tightened some
return False
if potList['CDIH'].violations()>0:
return False
if potList['BOND'].violations()>0:
return False
if potList['ANGL'].violations()>0:
return False
if potList['IMPR'].violations()>1:
return False
return True
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
# 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
#
# high temperature bit - using only P-P nonbonded terms
# protocol.initNBond(repel=1.2,
# cutnb=100,
# tolerance=45,
# selStr="name P")
protocol.initDynamics(dyn,
potList=potList, # potential terms to use
bathTemp=init_t,
initVelocities=1,
finalTime=10, # stops at 10ps or 5000 steps
numSteps=5000, # whichever comes first; was 20 and 10000
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()
# final all- atom minimization
#
protocol.initMinimize(minc,
potList=potList,
dEPred=10)
minc.run()
#do analysis and write structure
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,rdcs,potList['CDIH'],potList['RAMA']],
#averageCrossTerms=refRMSD,
averageTopFraction=0.5, #report only on best 50% of structs
averageAccept=accept, #only use structures which pass accept()
averageContext=FinalParams(rampedParams),
averageFilename="REF2/SCRIPT_ave.pdb", #generate regularized ave structure
averageFitSel="not name H* and not resname ANI",
averageCompSel="not name H* and not resname ANI" ).run()
_______________________________________________
Xplor-nih mailing list
[email protected]
https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
