Hello,
I am trying to determine the structure of a peptide from NOE constraints.
For this purpose, I used the following scripts which did not include the
terms corresponding to the dihedral angles of the force field, i.e., it did
not include the term "DIHE". Thus the structures obtained for one of the
compounds presented values of Phi and Psi that are outside the Ramachandran
diagram.
I tried to solve this by including the term "DIHE". When I did that many
NOE restraints gave problems, and by reviewing them and weakening some of
them I got structures that do get into the allowed zones of the
Ramachandran map but are very different from each other.
Is there any problem with using the DIHE term? How to solve the problem of
the dihedral angles or the high RMSD in the structures?
Thank you in advance,
Rafael RodrÃguez
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
xplor.requireVersion("2.24")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2013/07/10
#
# this checks for typos on the command-line. User-customized arguments can
# also be specified.
#
xplor.parseArguments()
# filename for output structures. This string must contain the STRUCTURE
# literal so that each calculated structure has a unique name. The SCRIPT
# literal is replaced by this filename (or stdin if redirected using <),
# but it is optional.
#
inFilename = "out_sa/anneal300iso_*.sa" #structures to analyze
outFilename = "out_refine/SCRIPT_STRUCTURE.sa" #fit output structures
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#Inicializar par'ametros m'ios
#protocol.parameters['protein']="parallhdg.pro"
protocol.initParams(("/nfs/home/rafaelr/protein-3.2-mod.par"))
protocol.initStruct("CIGB-300-iso-21.psf")
command = xplor.command
# read an existing model
#
protocol.initCoords("out_sa/anneal300iso_ave.pdb")
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="out_sa/anneal300iso_ave.pdb",
psf=("CIGB-300-iso-21.psf"),
cmpSel="not name H*")
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"restr.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
pot.setPotType("soft") # if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
rampedParams.append( MultRamp(2,30, "noe.setScale( VALUE )") )
## set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","jna_CI.tbl",
# A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
## Set up dihedral angles
from xplorPot import XplorPot
#protocol.initDihedrals("dihe_CI.tbl",
# 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)") )
# hbdb - distance/angle bb hbond term
#
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
# Statistical torsion angle potential
#
#from torsionDBPotTools import create_TorsionDBPot
#torsionDBPot = create_TorsionDBPot('tDB')
#potList.append( torsionDBPot )
#rampedParams.append( MultRamp(.002,2,"torsionDBPot.setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond(nbxmod=4)") )
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=0)""") )
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()
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
protocol.torsionTopology(dyn)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
protocol.cartesianTopology(minc)
# 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=25,
tempStep =12.5,
ivm=dyn,
rampedParams = rampedParams)
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>1:
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()>0:
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
#
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
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,
numSteps=1000,
printInterval=50
)
minc.run()
#do analysis and write structure when this routine is finished.
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(structLoopAction=calcOneStructure,
pdbFilesIn=inFilename,
pdbTemplate=outFilename,
# numStructures=numberOfStructures,
calcMissingStructs=True, #calculate only missing structures
doWriteStructures=True, #analyze and write coords after calc
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH'],potList['VDW']], #
averageCrossTerms=refRMSD,
averageTopFraction=0.1, #report only on best 5% of structs
averageAccept=accept, #only use structures which pass accept()
averageContext=FinalParams(rampedParams),
averageFilename="out_refine/SCRIPT_ave.pdb", #generate
regularized ave structure
averageFitSel="name CA or name N or name C or name O ",
averageCompSel="not name H*" ).run()
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
xplor.requireVersion("2.38")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2013/07/10
#
# this checks for typos on the command-line. User-customized arguments can
# also be specified.
#
xplor.parseArguments()
# filename for output structures. This string must contain the STRUCTURE
# literal so that each calculated structure has a unique name. The SCRIPT
# literal is replaced by this filename (or stdin if redirected using <),
# but it is optional.
#
inFilename = "out_refine/refine300iso*.sa" #structures to analyze
outFilename = "out_minimiz/SCRIPT_STRUCTURE.min" #fit output structures
#numberOfStructures=100 #usually you want to create at least 20 desactivado
si se lee
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#Inicializar par'ametros m'ios
#protocol.parameters['protein']="charmm22/parallh22x.pro"
protocol.initParams(("protein-3.2-mod.par"))
protocol.initStruct("CIGB-300-iso-21.psf")
command = xplor.command
# read an existing model
#
protocol.initCoords("out_refine/refine300iso_ave.pdb")
#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 StaticRamp, InitialParams
minimParams=[]
# 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='out_refine/refine300iso_ave.pdb',
psf=("CIGB-300-iso-21.psf"),
cmpSel="not name H*")
# hbdb - distance/angle bb hbond term
#
from xplorPot import XplorPot
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
# Statistical torsion angle potential
#
#from torsionDBPotTools import create_TorsionDBPot
#torsionDBPot = create_TorsionDBPot('tDB')
#potList.append( torsionDBPot )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList.append( XplorPot('VDW') )
potList.append( XplorPot("BOND") )
potList.append( XplorPot("ANGL") )
potList['ANGL'].setThreshold( 5 )
potList.append( XplorPot("IMPR") )
potList['IMPR'].setThreshold( 5 )
# 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
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
protocol.cartesianTopology(minc)
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
# final all- atom minimization
#
protocol.initMinimize(minc,
numSteps=15000,
potList=potList,
printInterval=500,
dEPred=10
)
minc.run()
#do analysis and write structure when this routine is finished.
pass
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"restr.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
# pot.setPotType("soft") # if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
## set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","jna_CI.tbl",
# A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
# Set up dihedral angles
from xplorPot import XplorPot
#protocol.initDihedrals("dihe_CI.tbl",
# useDefaults=False # by default, symmetric sidechain
# restraints are included
# )
potList.append( XplorPot('CDIH') )
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>0:
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()>0:
return False
return True
from simulationTools import StructureLoop, FinalParams
StructureLoop(structLoopAction=calcOneStructure,
pdbFilesIn=inFilename,
pdbTemplate=outFilename,
# numStructures=numberOfStructures,
calcMissingStructs=True, #calculate only missing structures
doWriteStructures=True, #analyze and write coords after calc
genViolationStats=1,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH'],potList['VDW']],
#potList['CDIH'],
averageCrossTerms=refRMSD,
averageTopFraction=0.1, #report only on best 10% of structs
averageAccept=accept, #only use structures which pass accept()
# averageContext=FinalParams(rampedParams),
averageFilename="out_minimiz/SCRIPT_ave.pdb", #generate
regularized ave structure
averageFitSel="name CA",
averageCompSel="not name H*" ).run()
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
xplor.requireVersion("2.24")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, J-coupling restraints.
#
# this script performs annealing from an extended structure.
# It is faster than the original anneal.py
#
# CDS 2009/07/24
#
# this checks for typos on the command-line. User-customized arguments can
# also be specified.
#
xplor.parseArguments()
# filename for output structures. This string must contain the STRUCTURE
# literal so that each calculated structure has a unique name. The SCRIPT
# literal is replaced by this filename (or stdin if redirected using <),
# but it is optional.
#
outFilename = "out_sa/SCRIPT_STRUCTURE.sa"
numberOfStructures=200 #usually you want to create at least 20
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed() #set random seed - by time
#command = xplor.command('parameter @TOPPAR:parallhdg.pro
@TOPPAR:aldrin_par_1.0.pro end')
command = xplor.command #averiguar para que esta aqui
#Inicializar par'ametros m'ios
#protocol.parameters['protein']="parallhdg.pro"
protocol.initParams(("/nfs/home/rafaelr/protein-3.2-mod.par"))
#Inicializar estructura
protocol.initStruct("CIGB-300-iso-21.psf")
# generate random extended initial structure with correct covalent geometry
#
protocol.genExtendedStructure()
#Inicializar coordenadas
#protocol.initCoords("template_VAD315.pdb")
#
# 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='modelo.pdb',psf='AVV036.psf',
# cmpSel="not name H*")
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('all',1,"restr.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
pot.setPotType("soft") # if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
rampedParams.append( MultRamp(2,30, "noe.setScale( VALUE )") )
## set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","jna_CI.tbl",
# A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
## Set up dihedral angles
from xplorPot import XplorPot
#dihedralRestraintFilename="dihe_CI.tbl"
#protocol.initDihedrals(dihedralRestraintFilename,
# 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 )
# hbdb - hbond database-based term
# CHEQUEAR SI SE INCLUYE O NO
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
#New torsion angle database potential
#from torsionDBPotTools import create_TorsionDBPot
#torsionDB = create_TorsionDBPot('torsionDB')
#potList.append( torsionDB )
#rampedParams.append( MultRamp(.002,2,"torsionDB.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,
onlyCA=0)""") )
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, configure bath/molecule friction coeff.
#
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()
# configure ivm topology for torsion-angle dynamics
#
protocol.torsionTopology(dyn)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
protocol.cartesianTopology(minc)
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3500. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =12.5,
ivm=dyn,
rampedParams = rampedParams)
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
# generate a new structure with randomized torsion angles
#
from monteCarlo import randomizeTorsions
randomizeTorsions(dyn)
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
# protocol.fixupCovalentGeom(maxIters=1000,useVDW=1)
# set torsion angles from restraints
#
from torsionTools import setTorsionsFromTable
# setTorsionsFromTable(dihedralRestraintFilename)
protocol.writePDB(loopInfo.filename()+".init")
# calc. initial tensor orientation
#
# 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
#
protocol.initDynamics(dyn,
potList=potList, # potential terms to use
bathTemp=init_t,
initVelocities=1,
finalTime=100, # stops at 100ps or 1000 steps
numSteps=1000, # whichever comes first
# finalTime=1000, # stops at 800ps or 8000 steps
# numSteps=50000, # whichever comes first
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
# numSteps=1000, #at each temp: 1000 steps or
# finalTime=.4 , # .4ps, whichever is less
printInterval=100)
# perform simulated annealing
#
cool.run()
# final torsion angle minimization
#
protocol.initMinimize(dyn,
printInterval=50) #adicionar numSteps=1000
dyn.run()
# final all- atomic degrees of freedom minimization
#
protocol.initMinimize(minc,
potList=potList,
dEPred=10
)
minc.run()
#do analysis and write structure when this function retunrs
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
doWriteStructures=True,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=True,
averageTopFraction=0.1, #report stats on best 10% of structs
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH'],potList['VDW']],
averageContext=FinalParams(rampedParams),
averageFitSel=("name CA or name C or name N or name O"),
#criterio de selecci'on
# averageCrossTerms=refRMSD,
averageFilename="out_sa/SCRIPT_ave.pdb",
averagePotList=potList).run()
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
