There are 2 parts to this question.
1.Background is: we have managed to make the default scripts run and create the
right thioether linkages and Dha/Dhb residues, but we have two peptides of 28
and 20 aa that have lysine sidechain protons all the way up to 0.9 and 0.99
ppm-referenced to DSS (and there are aromatic residues in the peptide as well).
This is highly suggestive that the lysine HD/HG should stack on an aromatic
sidechain, right? in at least one of the structures, the Lys sidechain could
easily be rotated to come in proximity with the aromatic ring but the
calculation results don't place it there, there are NOEs to support this, but
the calculations don't favor it, in fact such a structure is nowhere in the
lowest 20/100. How do I implement this?
2. I've tried mimicking the csPotTools.p y script, but the calculation won't
run at all saying "cannot convert 4 letter residue to 3 letter" (because we
defined cyclized aminobutyrine/alanine as D-Abu or Abu or D-Ala, ALA depending
on the stereochemistry-see below). I have a feeling even if they were defined
as 3 letter, camShift or whatever it is that does this wouldn't recognize them.
I apologize if I don't know some of the inner workings I should, but I can't
find a comprehensive explanation on how to do this anywhere.
See below for the specific error:
k=create_CSPot("ProcA21_csPot.tbl",cs )
Error converting String DABU into amino acid index: not a valid 3-letter code
Warning: Did not recognize the following atomName: CG
Full info: [DABU][CG][47][3]
Warning: Did not recognize the following atomName: HG1
Full info: [DABU][HG1][48][3]
Warning: Did not recognize the following atomName: HG2
Full info: [DABU][HG2][49][3]
Warning: Did not recognize the following atomName: HG3
Full info: [DABU][HG3][50][3]
Error converting String DALA into amino acid index: not a valid 3-letter code
Error converting String DABU into amino acid index: not a valid 3-letter code
Warning: Did not recognize the following atomName: CG
Full info: [DABU][CG][146][12]
Warning: Did not recognize the following atomName: HG1
Full info: [DABU][HG1][147][12]
Warning: Did not recognize the following atomName: HG2
Full info: [DABU][HG2][148][12]
Warning: Did not recognize the following atomName: HG3
Full info: [DABU][HG3][149][12]
Error converting String DABU into amino acid index: not a valid 3-letter code
Warning: Did not recognize the following atomName: CG
Full info: [DABU][CG][237][18]
Warning: Did not recognize the following atomName: HG1
Full info: [DABU][HG1][238][18]
Warning: Did not recognize the following atomName: HG2
Full info: [DABU][HG2][239][18]
Warning: Did not recognize the following atomName: HG3
Full info: [DABU][HG3][240][18]
CamShift ERROR: Side chain 0 (CYS) incomplete
Could not find atom(s) HG1
CamShift ERROR: Side chain 1 (CYS) incomplete
Could not find atom(s) HG1
CamShift ERROR: Side chain 3 (DABU) incomplete
Could not find atom(s) HB1, HB3
CamShift ERROR: Side chain 9 (DALA) incomplete
Could not find atom(s) HB1
CamShift ERROR: Side chain 12 (DABU) incomplete
Could not find atom(s) HB1, HB3
CamShift ERROR: Side chain 17 (CYS) incomplete
Could not find atom(s) HG1
CamShift ERROR: Side chain 18 (DABU) incomplete
Could not find atom(s) HB1, HB3
CamShift ERROR: Side chain 26 (CYS) incomplete
Could not find atom(s) HG1
CamShift ERROR: Dihedral angle 2 of residue 9 (DALA) incomplete
Chemical shift calculation for this residue likely to be
incorrect
refine_ProcA21.py(711): from atomSel import AtomSel
refine_ProcA21.py(712): k.setDomainSel(AtomSel("not PSEUDO",esim))
Traceback (most recent call last):
File "<string>", line 2, in <module>
File "/Users/silviabobeica/xplor-nih-2.41.1/python/trace.py", line 180, in run
exec cmd in dict, dict
File "<string>", line 1, in <module>
File "refine_ProcA21.py", line 712, in <module>
k.setDomainSel(AtomSel("not PSEUDO",esim))
NameError: name 'esim' is not defined
PyInterp::command: error executing: >execfile('refine_ProcA21.py')<
/usr/local/bin/xplor: line 548: 4987 Abort trap: 6
These are the shifts of one of these lysines.
assign ( resid 17 and name HA ) 3.871 0.003
assign ( resid 17 and name HB1 ) 1.770 0.006
assign ( resid 17 and name HB2 ) 1.672 0.005
assign ( resid 17 and name HD# ) 1.567 0.008
assign ( resid 17 and name HE# ) 2.948 0.003
assign ( resid 17 and name HG1 ) 0.995 0.001
assign ( resid 17 and name HG2 ) 0.904 0.010
assign ( resid 17 and name HN ) 7.324 0.002
I'm going to paste below all the modifications we made to xplor to make these
bonds:
Can you please help us input those lysine 0.9ppm shifts to make the calculation
results relate with these observations by telling us where in the script and
what needs to be added?
.par file changes and additions
! the bond lengths were also reference to the Turpin 2014 paper, 5/14/2017
bond CXX HA $kbon 1.10 ! defined for DHB and DHA
LZ 12/14/2016
bond CXX CT $kbon 1.509 ! defined for DHB LZ
12/14/2016
bond CXX CXX $kbon 1.352 ! defined for DHB and DHA
LZ 12/14/2016
bond CXX NH1 $kbon 1.410 ! defined for DHB and DHA
LZ 12/14/2016
bond CXX C $kbon 1.530 ! defined for DHB and DHA
LZ 12/14/2016
! based on Turpin paper: RCS Adv 2014, 4, 48621-48631
! New CHARMM force field parameters for dehydrated amnio acid resdies, the key
to Lantibiotic Molecular
! dynamic simulation
! added 5/14/2017, L Zhu
angle CXX CXX C 60.00 120.0 ! defined for DHA and DHB
angle CXX CXX NH1 80.00 128.0 ! defined for DHA and DHB
angle C CXX NH1 80.00 110.0 ! defined for DHA and DHB
angle CXX CXX HA 45.00 120.5 ! defined for DHA and DHB
angle CXX C NH1 80.00 116.5 ! defined for DHA and DHB
angle CXX C O 80.00 122.5 ! defined for DHA and DHB
angle CXX C N 80.0 119.0 ! defined for DHB-Pro link
angle CXX NH1 C 50.00 120.0 ! defined for DHA and DHB
angle CXX NH1 H 34.00 117.0 ! defined for DHA and DHB
angle CXX CXX CT 43.50 126.5 ! defined for DHB
angle HA CXX CT 43.50 118.0 ! defined for DHB
angle HA CXX HA 43.50 118.0 ! defined for DHA
angle CXX CT HA 80.0 109.5 ! defined for DHB
improper C CXX +NH1 O $kpla 0 180.0 ! defined for DHA
and DHB LZ 12/15/2016
improper C CXX N O $kpla 0 180.0 ! defined for DHA
and DHB -pro LZ 12/15/2016
improper CXX C NH1 CT $kpx 0 180.0 ! defined for DHA and
DHB, backbone to next resid
improper CXX C N CT $kback 0 180.0 ! NH peptide
planarity, DHB link to RPO
improper HA HA CXX HA $kchi 0 -66.514 ! methyl
improper O C NH1 CXX $kpx 0 0.0 ! CO peptide
planarity, front
improper CT C NH1 CXX $kpx 0 180.0 ! new trans peptide
bond
improper H NH1 C CXX $kpx 0 0.0 ! NH peptide planarity
! based on Turpin's paper, 5/15/2017
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dihedral C CXX CXX HA 3.2000 2 180.0 ! DHA and DHB
dihedral NH1 CXX CXX HA 4.0000 2 180.0 ! DHA and DHB: average
from 3 for DHB, 5 for DHA
dihedral NH1 CXX CXX CT 3.5000 2 180.0 ! DHB
dihedral NH1 CXX CXX CT 1.7000 1 180.0 ! DHB
dihedral NH1 CXX CXX CT 0.5000 3 180.0 ! DHB
dihedral CXX CXX CT HA 0.3000 3 180.0 ! DHB
dihedral CXX CXX C O 0.5700 1 0.0 ! DHB and DHA
dihedral CXX CXX C O 0.5900 2 0.0 ! DHB and DHA
dihedral CXX CXX C O 2.9600 3 0.0 ! DHB and DHA
dihedral CXX CXX C NH1 2.1000 1 180.0 ! DHB and DHA
dihedral CXX CXX C NH1 2.8800 2 0.0 ! DHB and DHA
dihedral CXX CXX C NH1 2.0000 3 0.0 ! DHB and DHA
dihedral HA CXX CXX C 3.2000 2 180.0 ! DHB and DHA
dihedral CT CXX CXX C 3.5000 2 180.0 ! DHB
dihedral CT CXX CXX C 1.7000 1 0.0 ! DHB
dihedral CT CXX CXX C 0.5000 3 0.0 ! DHB
dihedral CXX NH1 C CT 0.5000 3 0.0 ! DHB
dihedral C CXX NH1 C 0.5000 3 0.0 ! DHB
dihedral H NH1 CXX CXX 2.5000 2 180.0 ! DHB and DHA
dihedral CXX CXX NH1 H 2.5000 2 180.0 ! DHA and DHB
dihedral C CXX NH1 H 2.5000 2 180.0 ! DHA and DHB
dihedral CXX C NH1 H 2.5000 2 180.0 ! DHA and DHB
dihedral CXX CXX NH1 C 0.5000 1 180.0 ! DHA and DHB
dihedral CXX CXX NH1 C 0.5000 2 0.0 ! DHA and DHB
dihedral CXX CXX NH1 C 0.2000 3 180.0 ! DHA and DHB
! from Turpin paper 5/23/2017
improper C CXX NH1 O 120.000 0 0.0 ! DHA and DHB CO
peptide planarity, self
improper CXX CXX C NH1 72.000 0 0.0 ! CA chirality, self
dihedral CT CT S CT 0.2400 1 180.0 ! chi1 - chi4 defined
for ALAS and ABUS
dihedral CT CT S CT 0.3700 3 0.0 ! chi1 - chi4 defined
for ALAS and ABUS
dihedral CT S CT CT 0.2400 1 180.0 ! chi1 - chi4 defined
for ALAS and ABUS
dihedral CT S CT CT 0.3700 3 0.0 ! chi1 - chi4 defined
for ALAS and ABUS
dihedral CT S CT HA 0.2800 3 0.0 ! chi1 - chi4 defined
for ALAS and ABUS
dihedral CT CT CT NH1 $kdih 3 0.0 ! chi1 - chi4, defined
for ABU, Zhu, 05-02-2017
NONBonded CXX 0.1450 3.2072 0.1450 3.2072 ! LZ
_______________________
.top file changes/additions:
residue DALA ! D-ALA, added by zhu, L 1/2/2015 based on ALA
group
atom N type=NH1 charge=-0.36 end
atom HN type=H charge= 0.26 end
group
atom CA type=CT charge= 0.00 end
atom HA type=HA charge= 0.10 end
group
atom CB type=CT charge=-0.30 end
atom HB1 type=HA charge= 0.10 end
atom HB2 type=HA charge= 0.10 end
atom HB3 type=HA charge= 0.10 end
group
atom C type=C charge= 0.48 end
atom O type=O charge=-0.48 end
bond N HN
bond N CA bond CA HA
bond CA CB bond CB HB1 bond CB HB2 bond CB HB3
bond CA C
bond C O
improper HA C N CB ! change from L-aa: a N c b, to D-aa: a C N b, stereo CA
improper HB1 HB2 CA HB3 !stereo CB
end
residue ABU ! define ABU, Zhu, L based on VAL, 1/2/2015
group
atom N type=NH1 charge=-0.36 end
atom HN type=H charge= 0.26 end
group
atom CA type=CT charge= 0.00 end
atom HA type=HA charge= 0.10 end
group
atom CB type=CT charge=-0.20 end
atom HB1 type=HA charge= 0.10 end
atom HB2 type=HA charge= 0.10 end ! add HB2 proton
group
atom CG type=CT charge=-0.30 end
atom HG1 type=HA charge= 0.10 end
atom HG2 type=HA charge= 0.10 end
atom HG3 type=HA charge= 0.10 end
group
atom C type=C charge= 0.48 end
atom O type=O charge=-0.48 end
bond N HN
bond N CA bond CA HA
bond CA CB bond CB HB1 bond CB HB2
bond CB CG bond CG HG1 bond CG HG2 bond CG HG3
bond CA C
bond C O
improper HA N C CB !stereo CA
improper HB1 HB2 CA CG !stereo CB, HA HA CT CT: -78.874 methylene
improper HG1 HG2 CB HG3 !stereo CG, HA HA CT HA: -66.514 methyl
dihedral CG CB CA N ! dihedral X CT CT X kdih 3 0.0 chi1-chi4, where
3=n, multiplicity, 0.0=phase shift, aleady defined in par file as X CT CT X
end
residue DABU ! define DABU, Zhu, L/Silvia based on ABU, 10/28/2016
group
atom N type=NH1 charge=-0.36 end
atom HN type=H charge= 0.26 end
group
atom CA type=CT charge= 0.00 end
atom HA type=HA charge= 0.10 end
group
atom CB type=CT charge=-0.20 end
atom HB1 type=HA charge= 0.10 end
atom HB2 type=HA charge= 0.10 end ! add HB2 proton
group
atom CG type=CT charge=-0.30 end
atom HG1 type=HA charge= 0.10 end
atom HG2 type=HA charge= 0.10 end
atom HG3 type=HA charge= 0.10 end
group
atom C type=C charge= 0.48 end
atom O type=O charge=-0.48 end
bond N HN
bond N CA bond CA HA
bond CA CB bond CB HB1 bond CB HB2
bond CB CG bond CG HG1 bond CG HG2 bond CG HG3
bond CA C
bond C O
improper HA C N CB !stereo CA ! change stero from L-ABU to D-ABU: HA
N C CB
improper HB1 HB2 CA CG !stereo CB
improper HG1 HG2 CB HG3 !stereo CG
dihedral CG CB CA N
end
residue DHB ! define DHB, Zhu, L based on ABU, 1/7/2015
group
atom N type=NH1 charge=-0.36 end
atom HN type=H charge= 0.26 end
atom CA type=CXX charge= 0.10 end ! charge changed, 5/4/2017 to make
zero charge of HN-N-CA
group
atom CB type=CXX charge=-0.10 end ! charge changed!, 5/4/2017 to make
zero of CB-HB
atom HB type=HA charge= 0.10 end ! charge changed!
group
atom CG type=CT charge=-0.30 end
atom HG1 type=HA charge= 0.10 end
atom HG2 type=HA charge= 0.10 end
atom HG3 type=HA charge= 0.10 end
group
atom C type=C charge= 0.48 end
atom O type=O charge=-0.48 end
bond N HN
bond N CA
bond CA CB bond CB HB
bond CB CG bond CG HG1 bond CG HG2 bond CG HG3
bond CA C
bond C O
improper HG1 HG2 CB HG3 !stereo CG
improper CA CB C N !stereo CA, Zhu 5/10/2017, based on Turpin's paper
end
residue DHA ! define DHA, Zhu, L based on DHB, 1/7/2015L
group
atom N type=NH1 charge=-0.36 end
atom HN type=H charge= 0.26 end
atom CA type=CXX charge= 0.10 end ! charge changed, 5/4/2017
group
atom CB type=CXX charge=-0.20 end ! charge changed! 5/4/2017
atom HB1 type=HA charge= 0.10 end ! charge changed!
atom HB2 type=HA charge= 0.10 end ! charge changed!
group
atom C type=C charge= 0.48 end
atom O type=O charge=-0.48 end
bond N HN
bond N CA
bond CA CB bond CB HB1 bond CB HB2
bond CA C
bond C O
improper CA CB C N !stereo CA, Zhu 5/10/2017, based on Turpin's paper
end
presidue ABUS ! S-link between HB1 of ABU and CYS ...ABU - S -
CYS...
! added by Zhu, L based on DISU 1/2/2015
group
delete atom 1HB1 end
! modify atom 1CB charge= -0.10 end check charge? good, only one HB
protons left
modify atom 1CB charge= 0.10 end ! to make the new bond charge to zero
group
delete atom 2HG end
! modify atom 2CB charge= 0.20 end
! modify atom 2SG charge=-0.20 end
modify atom 2SG charge=-0.10 end ! make the net charge to zero for
the new bond
add bond 1CB 2SG
add angle 1CA 1CB 2SG
add angle 1CB 2SG 2CB
add angle 1HB2 1CB 2SG ! 5/9/2017
add dihedral 1CA 1CB 2SG 2CB ! Zhu 5/9/2017
add dihedral 1CB 2SG 2CB 2CA
end
presidue ALAS ! S-link bet. S link between ALA and CYS with a
single -C-S-C-... added by Zhu 1/7/2015
group
delete atom 1HB1 end
! modify atom 1CB charge= -0.20 end two HB protons left, charge is the
only differnce from ABUS
modify atom 1CB charge= 0.20 end ! to make the new bond to charge zero
group
delete atom 2HG end
! modify atom 2CB charge= 0.20 end
! modify atom 2SG charge=-0.20 end
modify atom 2SG charge=-0.20 end ! to make the new bond to charge zero
add bond 1CB 2SG
add angle 1CB 2SG 2CB
add angle 1CA 1CB 2SG
add angle 1HB2 1CB 2SG
add angle 1HB3 1CB 2SG ! 5/9/2017
add dihedral 1CA 1CB 2SG 2CB ! Zhu 5/9/2017
add dihedral 1CB 2SG 2CB 2CA
_________
This is our script:
# 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.pdb"
numberOfStructures=100
if quick:
numberOfStructures=3
pass
# protocol module has many high-level helper functions.
#
import protocol
# explicitly set random seed
protocol.initRandomSeed()
# protocol.initRandomSeed() # set random seed - by time
#
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
# seqToPSF('cylll.seq')
# protocol.initStruct("ave_chantal.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())
#protocol.genExtendedStructure()
# or read an existing model
#
protocol.loadPDB("ProcA21_template.pdb")
xplor.simulation.deleteAtoms("not known")
xplor.command("""
patch ABUS reference=1=( residue 4 ) reference=2=( residue 1 ) end
patch ALAS reference=1=( residue 10 ) reference=2=( residue 2 ) end
patch ABUS reference=1=( residue 13 ) reference=2=( residue 18 ) end
patch ABUS reference=1=( residue 19 ) reference=2=( residue 27 ) end
""")
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='ProcA21.pdb',
# 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,"ProcA21_noe.tbl"),('hbond',1,"ProcA21_hbond.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",file="jna_coup.tbl",
A=6.98,B=-1.38,C=1.72,phase=-60.0)
jCoup.setThreshold(1.0),
#jCoup.potType='square'
jCoup.potType='soft'
jCoup.setScale(10)
potList.append(jCoup)
print jCoup.info()
# potList[jCoup].setThreshold( 5 )
# Set up dihedral angles
from xplorPot import XplorPot
dihedralRestraintFilename="ProcA21_dihed.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 ) #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)") )
## # hbda - distance/angle bb hbond term
## #
## protocol.initHBDA('hbda.tbl')
## potList.append( XplorPot('HBDA') )
## or hbdb - hbond databas-based term
## protocol.initHBDB()
## potList.append( XplorPot('HBDB') )
#Rama torsion angle database
#
protocol.initRamaDatabase()
potList.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
## or use new torsion angle database potential: based on tietz file
#
## 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=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
#
from atomAction import SetProperty
## import varTensorTools
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
## varTensorTools.massSetup(media.values(),300)
AtomSel("all ").apply( SetProperty("fric",10.) )
# 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 = 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()>1000:
return False
if potList['VDW'].violations()>100: #this might be tightened some
return False
if potList['CDIH'].violations()>100:
return False
if potList['BOND'].violations()>100:
return False
if potList['ANGL'].violations()>100:
return False
if potList['IMPR'].violations()>100:
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)
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,potList['CDIH']],
averageTopFraction=0.2, #report only on best 50% of structs
averageAccept=accept, #only use structures which pass accept()
averageContext=FinalParams(rampedParams),
averageFilename="ProcA21_ave.pdb", #generate regularized ave
structure
averageFitSel="name CA",
averageCompSel="not resname ANI and not name H*" ).run()
Thank you so much for your help!
Silvia Bobeica
Grad. student
van der Donk lab
########################################################################
To unsubscribe from the XPLOR-NIH list, click the following link:
http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
