Hi everyone,
I am having trouble using GB potential through the python interface. It
raises a segmentation fault when I start minimization or MD of an IVM.
Strangely, this does not occur using the standard xplor interface (at least it
does not crash). I wrote a script inspired from test/gborn.inp which is
described below, it is just a test prior to a more sophisticated script. I run
xplor-nih 2.25 on a x86-64 RHEL 5.3 station.
If anyone as an idea to help, it would be welcome, because I am not
very at ease with native xplor-language, I prefer the python interface.
Thank you for your attention,
Olivier Serve
#testDyn.py
import protocol
import regularize
from xplor import command
from ivm import IVM
from xplor import select
protocol.initRandomSeed(3432)
#Load amber force field topology / parameters
xplor.command("topology @TOPPAR:topamber.inp @TOPPAR:amberpatches.pro end")
xplor.command("parameter @TOPPAR:paramber.gb.inp end")
xplor.command("parameter ANGLE HP CT N 50.0 109.50 end ! AA
general changed based on NMA nmodes")
from potList import PotList
potList = PotList()
from xplorPot import XplorPot
potList.append( XplorPot('VDW') )
bonds=XplorPot("BOND")
potList.append(bonds)
potList.append( XplorPot("ANGL") )
potList['ANGL'].setThreshold( 5 )
potList.append( XplorPot("IMPR") )
potList['IMPR'].setThreshold( 5 )
potList.append( XplorPot("GBIN") )
potList.append( XplorPot("GBSE") )
#Read Sequence File to create structure for amber topology
xplor.command("""
segment
name="A"
chain
first NMET TAIL nil * END
first NTER TAIL + * END
link PPG1 HEAD - * TAIL + GLY end
link PPG2 HEAD - GLY TAIL + * end
link PPGG HEAD - GLY TAIL + GLY end
link PEPP HEAD - * TAIL + PRO end
link PPGP HEAD - GLY TAIL + PRO end
link PEPT HEAD - * TAIL + * end
last CGLY HEAD nil GLY end
last CTER HEAD - * end
sequence @1ubq.seq end
end
end
""")
protocol.initCoords("refine_0.sa")
regularize.addUnknownAtoms()
#Load volumes for GB calculations
xplor.command("@TOPPAR:volumes.amber")
#Selection of atoms, only sidechains and final loop allowed to be modified
resiSel="1:71"
bbStr='(name C or name CA or name N or name O or name H or name HA)'
sideChainsSel=select('resid '+resiSel+' and '+bbStr)
protocol.initNBond()
protocol.initRamaDatabase()
#NBond parameters extracted from gborn.inp
xplor.command("vector do (charge = charge + .2306) (name ht%)")
xplor.command("vector do (rmsd = rmsd * 0.9) (all) ! reduce volumes by 10%")
xplor.command("parameter reduce selection=(all) overwrite=true mode=average end
end")
xplor.command("""parameter
nbonds
atom cdie trunc
e14fac=0.8333333
cutnb 500. ctonnb 480. ctofnb 490
tolerance=100.0
nbxmod 5 vswitch
wmin=1.0
end
end""")
xplor.command("parameters nbonds EPS=1.0 WEPS=80. offset=0.09 GBHCT end end")
#Using python xplor interface
myDyn=IVM()
myMin=IVM()
myDyn.setGroupList(sideChainsSel)
myMin.setGroupList(sideChainsSel)
protocol.torsionTopology(myDyn)
protocol.initDynamics(myDyn, potList=potList, bathTemp=303, initVelocities=1,
finalTime=5, printInterval=500)
myDyn.setETolerance(3)
protocol.cartesianTopology(myMin)
protocol.initMinimize(myMin, potList=potList, numSteps=100, printInterval=25)
myMin.run()
myDyn.run()
#Using native xplor interface (Inspired from gborn.inp)
xplor.command("""
vector do (rmsd = rmsd * 0.9) (all) ! reduce volumes by 10%
flags include gbse gbin end
parameter reduce selection=(all) overwrite=true mode=average end end
energy end
! weakly restrain Calpha position
vector do (harmonic = 0.5) (name ca)
coor disp=reference @refine_0_xplor.sa
constraints harmonic end
flags include harmonic end
!-----------------
! Now run dynamics
vector do (vx = maxwell(250)) (all)
vector do (vy = maxwell(250)) (all)
vector do (vz = maxwell(250)) (all)
!-------------------------------------
! relax initial bad contacts
minimize rigid group = (resid 1:71 and (name C or name CA or name N or name O
or name H or name HA)) nstep=10 drop=10 end
minimize rigid group = (resid 1:71 and (name C or name CA or name N or name O
or name H or name HA)) nstep=10 drop=10 end
!--------------
! let it go
dynamics internal
tbath=303
group = (resid 1:71 and (name C or name CA or name N or name O or name H
or name HA))
auto torsion
nstep=500 timest=0.001 {ps} ! 5 ps dynamics
nprint=25 ! statistics output
end
""")
protocol.writePDB("testDyn.pdb")
