Hello Charles
After a few attempts to perform dynamics with rigid bodies under RDC
constraints for my multidomain protein, I have to ask again for your help.
I tried to adapt the dock.py script from the eginiput/diffTens directory.
However, instead of getting 10 different structures, all my structures are
exactly the same. IS there something obvious for you that I did wrong?
Below is my script. I hope you will find some time to give a look at it.
Thanks anyway
H?l?ne D?m?n?
xplor.requireVersion("2.18")
#
# 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(["quick"]) # check for command-line typos
quick=False
for opt in opts:
if opt[0]=="quick": #specify -quick to just test that the script runs
quick=True
pass
pass
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=20
if quick:
numberOfStructures=3
pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(34210) #explicitly set random seed
#
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('protG.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("start1.pdb")
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
potList1 = PotList()
potList0= PotList()
potList2 = 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='2KUI1.pdb',
# cmpSel="not name H*")
# 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,addAxisAtoms
media={}
# medium Da rhombicity
for (medium,Da,Rh,resid) in [
('peg12', 26.6, 0.06,991),
('peg34', 24.9, 0.12,992) ,
('pha12', -11.7, 0.35,995),
('pha34', -7.3, 0.49,997)
]:
oTensor = create_VarTensor(medium, axis=addAxisAtoms(resid=resid))
oTensor.setDaMax(200)
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_toNH
rdcs = PotList('rdc')
for (medium,expt,file, scale) in \
[
('pha12','NH','pha12_nh.tbl' ,1),
('pha34','NH','pha34_nh.tbl' ,1),
('peg12','NH' ,'peg12_nh.tbl' ,1),
('peg34','NH','peg34_nh.tbl' ,1),
]:
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
# rdc.setPotType('square')
rdcs.append(rdc)
pass
potList0.append(rdcs)
potList1.append(rdcs)
potList2.append(rdcs)
rampedParams.append( MultRamp(0.05,5.0, "rdcs.setScale( VALUE )") )
# 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("calcTensorOrientation(media['%s'])" %
medium) )
pass
from simulationTools import analyze
print analyze(potList0)
# gyration volume term
#
# gyration volume term
#
#Rama torsion angle database
#
from xplorPot import XplorPot
protocol.initRamaDatabase()
potList2.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList2['RAMA'].setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList2.append( XplorPot('VDW') )
potList1.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,
onlyCA=1)""") )
potList2.append( XplorPot("BOND") )
potList2.append( XplorPot("ANGL") )
potList2['ANGL'].setThreshold( 5 )
rampedParams.append( MultRamp(0.4,1,"potList2['ANGL'].setScale(VALUE)") )
potList2.append( XplorPot("IMPR") )
potList2['IMPR'].setThreshold( 5 )
rampedParams.append( MultRamp(0.1,1,"potList2['IMPR'].setScale(VALUE)") )
# Give atoms uniform weights, except for the anisotropy axis
#
# IVM setup
# the IVM is used for performing dynamics and minimization in torsion-angle
# space, and in Cartesian space.
#
# 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("fixDa, fixRh") #fix tensor parameters
# varTensorTools.topologySetup(dyn,m) #setup tensor topology
#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()
# reset ivm topology for torsion-angle dynamics
#
from ivm import IVM
from ivm import PublicIVM
dyn_fix = IVM()
dyn_fix.reset()
dyn_fix.group( 'resid 358:421' )
dyn_fix.group( 'resid 426:488' )
dyn_fix.group( 'resid 496:552' )
dyn_fix.group( 'resid 556:625' )
import varTensorTools
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
pass
varTensorTools.topologySetup(dyn_fix,list=media.values())
protocol.torsionTopology(dyn_fix,oTensors=media.values())
dyn_free_sch = IVM()
dyn_free_sch.reset()
dyn_free_sch.group( 'resid 358:421 and \
(name CA or name C or name N or name O or
name HN)' )
dyn_free_sch.group( 'resid 426:488 and \
(name CA or name C or name N or name O or
name HN)')
dyn_free_sch.group( 'resid 496:552 and \
(name CA or name C or name N or name O or name HN)')
dyn_free_sch.group( 'resid 556:625 and \
(name CA or name C or name N or name O or name HN)')
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
pass
varTensorTools.topologySetup(dyn_free_sch,list=media.values())
protocol.torsionTopology(dyn_free_sch,oTensors=media.values())
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3000. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=10,
tempStep =12.5,
ivm=dyn_fix,
rampedParams = rampedParams)
# Give atoms uniform weights, except for the anisotropy axis
#
from atomAction import SetProperty
import varTensorTools
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
AtomSel("all ").apply( SetProperty("fric",10.) )
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
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.
"""
from simulationTools import analyze
print analyze(potList0)
# 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 )
# initial minimization
protocol.initMinimize(dyn_fix, potList=potList0,
printInterval=50)
dyn_fix.run()
dyn_fix.run()
protocol.initMinimize(dyn_fix, potList=potList0,
printInterval=50)
dyn_fix.run()
dyn_fix.run()
protocol.initMinimize(dyn_fix, potList=potList0,
printInterval=50)
dyn_fix.run()
dyn_fix.run()
dyn_fix.run()
# k2=0
# while k2 < 0:
# initialize parameters for cooling loop
InitialParams( rampedParams )
# initialize integrator for simulated annealing
#
protocol.initDynamics(dyn_fix,
potList=potList2,
numSteps=1000, #at each temp: 100 steps or
finalTime=1 , # 1 ps, whichever is less
printInterval=100)
# perform simulated annealing
#
cool.run()
# k2=k2+1
# print k2
# final torsion angle minimization
#
protocol.initMinimize(dyn_free_sch,potList=potList2,
printInterval=50)
dyn_fix.run()
# final all- atom minimization
#
#do analysis and write structure
loopInfo.writeStructure(potList2)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averagePotList=potList2,
averageSortPots=[potList2['BOND'],potList2['ANGL'],potList2['IMPR'],
rdcs],
# averageCrossTerms=refRMSD,
averageTopFraction=0.5, #report only on best 50% of structs
averageAccept=accept, #only use structures which pass
accept()
averageContext=FinalParams(rampedParams),
averageFilename="SCRIPT_ave.pdb", #generate regularized
ave structure
averageFitSel="name CA",
averageCompSel="not resname ANI and not name H*" ).run()
--
Passerelle antivirus CBS
