Dear Charles,
I have tried to build Cd-S clusters as you suggested before and then generate
PSF file and to calculate structures using anneal.py script with constraints of
NOE, RDC. I got many violations on RDC, Angel, NOE etc. The structures are bot
correct when compared with that obtained from Cyana as well as those similar
protein structures. Please see the attached file for the scripts I used. Any
suggestions?
I appreciate your help!
Hongyan
Dr. Hongyan Li
Department of Chemistry
The University of Hong Kong
Pokfulam Road
Hong Kong
seq="""
MET ASP PRO GLU THR CYS PRO CYS PRO SER GLY GLY SER CYS THR CYS ALA ASP SER CYS
LYS CYS GLU GLY CYS LYS CYS THR SER CYS LYS LYS SER CYS CYS SER CYS CYS PRO ALA
GLU CYS GLU LYS CYS ALA LYS ASP CYS VAL CYS LYS GLY GLY GLU ALA ALA GLU ALA GLU
ALA GLU LYS CYS SER CYS CYS GLN
"""
import psfGen
psfGen.seqToPSF(seq,startResid=1)
xplor.command("""
topo
MASS CD 112.4100 ! cadmium
RESIdue CD {cadmium}
GROUp
ATOM CD TYPE=CD CHARge=2.0 END
END {CD}
! patch to create a Cd-SG cluster
!
!
presidue CDSG
delete atom 2HG end
delete atom 3HG end
delete atom 4HG end
delete atom 5HG end
add bond 1CD 2SG
add bond 1CD 3SG
add bond 1CD 4SG
add bond 1CD 5SG
add angle 2SG 1CD 3SG
add angle 3SG 1CD 4SG
add angle 4SG 1CD 5SG
add angle 2CB 2SG 1CD
add angle 3CB 3SG 1CD
add angle 4CB 4SG 1CD
add angle 5CB 5SG 1CD
end
end
! add Cd atoms
segment
number=401
chain
sequence CD CD CD CD CD CD CD end
end
end
! create 1st cluster
patch CDSG
reference 1=(resid 401)
reference 2=(resid 16)
reference 3=(resid 20)
reference 4=(resid 25)
reference 5=(resid 30)
end
! create 2nd cluster
patch CDSG
reference 1=(resid 402)
reference 2=(resid 8)
reference 3=(resid 14)
reference 4=(resid 16)
reference 5=(resid 27)
end
! create 3rd cluster
patch CDSG
reference 1=(resid 403)
reference 2=(resid 6)
reference 3=(resid 8)
reference 4=(resid 22)
reference 5=(resid 25)
end
! create 4th cluster
patch CDSG
reference 1=(resid 404)
reference 2=(resid 35)
reference 3=(resid 34)
reference 4=(resid 45)
reference 5=(resid 49)
end
! create 5th cluster
patch CDSG
reference 1=(resid 405)
reference 2=(resid 35)
reference 3=(resid 37)
reference 4=(resid 38)
reference 5=(resid 51)
end
! create 6th cluster
patch CDSG
reference 1=(resid 406)
reference 2=(resid 38)
reference 3=(resid 42)
reference 4=(resid 45)
reference 5=(resid 67)
end
! create 7th cluster
patch CDSG
reference 1=(resid 407)
reference 2=(resid 51)
reference 3=(resid 64)
reference 4=(resid 66)
reference 5=(resid 67)
end
write psf output=MT3-all.psf end
""")
xplor.command("""
! new bond parameter
param
bond CD S 1000 2.778
angle S CD S 500 109
angle CT S CD 500 106
angle CD S CD 500 109
! eps sigma eps(1:4) sigma(1:4)
! (kcal/mol) (A)
! ---------------------------------------
nonbonded CD 0.01 1.942 0.01 1.942
end
""")
#read in coordinates
import protocol
protocol.initCoords("MT3-cns-test.pdb")
#delete atoms not in the pdb file
xplor.command("delete sele=(not known) end")
#print out bonds violated by more than 0.1A
xplor.command("print thresh=0.1 bond end")
from xplorPot import XplorPot
#make sure we can calculate bond energy
print XplorPot("BOND").calcEnergy().energy
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints. RDC force constance was increased from 0.001 to 2 Kcal mole-1 Hz-2
#
# 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=100
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
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
xplor.command("""
! new bond parameter
param
bond CD S 1000 2.778
angle S CD S 500 109
angle CT S CD 500 106
angle CD S CD 500 109
! eps sigma eps(1:4) sigma(1:4)
! (kcal/mol) (A)
! ---------------------------------------
nonbonded CD 0.01 1.942 0.01 1.942
end
""")
#
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF(open('protG.seq').read())
protocol.initStruct("MT3-all.psf") # - or from file (open ('filename.psf'))
#
# generate a random extended structure with correct covalent geometry
#
protocol.genExtendedStructure("MT3-all_extended_%d.pdb" % seed)
# or read an existing model
#
#protocol.initCoords("MT3-cns-test.pdb")
#protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
#
# a PotList conatins a list of potential terms. This is used to specify which
# terms are active during refinement.
#
potList = PotList()
# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams
rampedParams=[]
highTempParams=[]
# orientation Tensor - used with the dipolar coupling term
#
oTensor = create_VarTensor("oTensor")
oTensor.setDa(-5.6)
oTensor.setRh(0.271)
# dipolar coupling restraints for protein amide NH.
#
rdc = create_RDCPot(name="rdc",
oTensor=oTensor,
file="MT3-rdc.tbl")
rdc.setPotType('square')
rdc.setThreshold(0.5)
print rdc.info()
potList.append(rdc)
# set up NOE potential
from noePotTools import create_NOEPot
noe = create_NOEPot("noe",file="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 32:68",
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.001
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,2.,"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, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=structLoopAction,
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,rdc],
averageTopFraction=0.2, #report only on best 50% of structs
# averageAccept=accept, #only use structures which pass accept()
averageContext=FinalParams(rampedParams),
averageFilename="ave.pdb", #generate average structure
averageFitSel="name CA",
averageCompSel="not hydro" ).run()
_______________________________________________
Xplor-nih mailing list
[EMAIL PROTECTED]
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
