#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# CDS 5/14/03
# modified 7/2/03
#
# this checks for typos on the command-line. User-customized arguments could
# also be specified.
#
xplor.parseArguments()
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=2
seed=3421
simWorld.setRandomSeed(seed) #set random seed
command = xplor.command
#
# import necessary modules
#
from xplorPot import XplorPot
from rdcPotTools import create_RDCPot
from varTensorTools import create_VarTensor
import varTensorTools
from ivm import IVM
from potList import PotList
import protocol
#
# generate PSF data from sequence and initialize the correct parameters.
#
from psfGen import seqToPSF
seqToPSF(open('protG.seq').read())
#
# generate a random extended structure with correct covalent geometry
#
protocol.genExtendedStructure("gb1_extended_%d.pdb" % seed)
#
# a PotList conatins a list of potential terms. This is used to specify which
# terms are active during refinement.
#
potList = PotList()
# orientation Tensor - used with the dipolar coupling term
#
oTensor = create_VarTensor("oTensor")
oTensor.setDa(-9.9)
oTensor.setRh(0.23)
# dipolar coupling restraints for protein amide NH.
#
rdc = create_RDCPot(name="rdc",
oTensor=oTensor,
file="bicelles_new_nh.tbl")
rdc.setThreshold(0.5)
print rdc.info()
potList.append(rdc)
# set up NOE potential
from noePotTools import create_NOEPot
noe = create_NOEPot("noe",file="protG_NOE.tbl")
noe.setPotType( "soft" )
noe.setRSwitch( 0.5 )
noe.setAsympSlope( 1. )
noe.setSoftExp(1.)
noe.setThreshold(0.5)
print noe.info()
potList.append(noe)
# set up J coupling
from jCoupPotTools import create_JCoupPot
jCoup = create_JCoupPot("jcoup",file="jna_coup.tbl",
A=6.98,B=-1.38,C=1.72,phase=-60)
jCoup.setThreshold(1.0)
print jCoup.info()
potList.append(jCoup)
# Set up dihedral angles
protocol.initDihedrals("protG_dih.tbl",
scale=5, #initial force constant
useDefaults=0)
potList.append( XplorPot('CDIH') )
# radius of gyration term - not used here.
#
#protocol.initCollapse("resid 4:20",
# scale=25.0,
# Rtarget=14.0)
#potList.append( XplorPot('COLL') )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
protocol.initNBond(nbxmod=4, # Can use 4 here, due to IVM dynamic
repel=0.5) # initial effective atom radius
potList.append( XplorPot('VDW') )
#
# annealing settings
#
init_t = 3500. # Need high temp and slow annealing to converge
final_t = 100
cool_steps = 15000 # slow annealing
#
# initial force constant settings
#
ini_rad = 0.4
ini_con = 0.004
ini_ang = 0.4
ini_imp = 0.4 # was 0.1
ini_noe = 1 # was 150
ini_sani = 0.01
potList.add( XplorPot("BOND") )
potList.add( XplorPot("ANGL") )
potList.add( XplorPot("IMPR") )
#
# potential terms used for high-temp dynamics
#
hitemp_potList = PotList()
hitemp_potList.add( XplorPot("BOND") )
hitemp_potList.add( XplorPot("ANGL") )
hitemp_potList.add( XplorPot("IMPR") )
hitemp_potList.add( XplorPot("CDIH") )
hitemp_potList.add( noe )
hitemp_potList.add( jCoup )
hitemp_potList.add( rdc )
# Give atoms uniform weights, except for the anisotropy axis
from atomAction import SetProperty
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
AtomSel("resname ANI ").apply( SetProperty("mass",300.) )
AtomSel("all ").apply( SetProperty("fric",10.) )
#
# IVM setup
#
dyn = IVM()
# 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")
oTensor.setFreedom("fix")
varTensorTools.topologySetup(dyn,oTensor)
protocol.initMinimize(dyn,numSteps=1000)
dyn.run()
#
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
dyn.potList().removeAll()
oTensor.setFreedom("fixDa, fixRh")
varTensorTools.topologySetup(dyn,oTensor)
protocol.torsionTopology(dyn)
#
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
from selectTools import IVM_groupRigidSidechain
IVM_groupRigidSidechain(minc)
minc.breakAllBondsIn("not resname ANI")
oTensor.setFreedom("varyDa, varyRh")
varTensorTools.topologySetup(minc,oTensor)
def setConstraints(k_ang,k_imp):
command("""
constraints
interaction (not resname ANI) (not resname ANI)
weights * 1 angl %f impr %f
end
end""" % ( k_ang, k_imp) )
return
def structLoopAction(loopInfo):
#
# this function calculates a single structure.
#
#
# set some high-temp force constants
#
noe.setScale( 20 ) #use large scale factor initially
rdc.setScale( ini_sani )
command("parameters nbonds repel %f end end" % ini_rad)
command("parameters nbonds rcon %f end end" % ini_con)
setConstraints(ini_ang, ini_imp)
# Initial weight--modified later
command("restraints dihedral scale=5. end")
#
# high temp dynamics
# note no Van der Waals term
#
init_t = 3500
ini_timestep = 0.010
bath = init_t
timestep = ini_timestep
protocol.initDynamics(dyn,
potList=hitemp_potList, # potential terms to use
bathTemp=bath,
initVelocities=1,
finalTime=20,
printInterval=100)
dyn.setETolerance( bath/100 ) #used to det. stepsize. default: bath/1000
dyn.run()
# increase dihedral term
command("restraints dihedral scale=200. end")
#
# cooling and ramping parameters
#
global k_ang, k_imp
# MultRamp ramps the scale factors over the specified range for the
# simulated annealing.
from simulationTools import MultRamp
k_noe = MultRamp(ini_noe,30.,"noe.setScale( VALUE )")
k_rdc = MultRamp(ini_sani,1.,"rdc.setScale( VALUE )")
radius = MultRamp(ini_rad,0.8,
"command('param nbonds repel VALUE end end')")
k_vdw = MultRamp(ini_con,4, "command('param nbonds rcon VALUE end end')")
k_ang = MultRamp(ini_ang,1.0)
k_imp = MultRamp(ini_imp,1.0)
protocol.initDynamics(dyn,
potList=potList,
bathTemp=bath,
initVelocities=1,
finalTime=2,
printInterval=100,
stepsize=timestep)
from simulationTools import AnnealIVM
AnnealIVM(initTemp =init_t,
finalTemp=100,
tempStep =25,
ivm=dyn,
rampedParams = [k_noe,k_rdc,radius,k_vdw,k_ang,k_imp],
extraCommands= lambda notUsed: \
setConstraints(k_ang.value(),
k_imp.value())).run()
#
# final torsion angle minimization
#
protocol.initMinimize(dyn,
printInterval=50)
dyn.run()
#
# final all atom minimization
#
protocol.initMinimize(minc,
potList=potList,
printInterval=100)
minc.run()
#
# analyze and write out structure
#
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=structLoopAction,
genViolationStats=1,
averagePotList=potList).run()