From: Hannes Loeffler (Hannes.Loeffler_at_stfc.ac.uk)
Date: Wed Oct 23 2013 - 02:39:29 CDT
Hi,
well, I am the author of this script and it certainly served me well
for my work published in DOI:10.1002/prot.24339. You may in particular
be interested in checking out the supplement where I did some tests
with CHARMM and AMBER on Ras-Raf.
A few notes though: The script below uses apbs through the iapbs
interface although you could of course try the other PB solvers in
CHARMM as well. The problem for me was that I couldn't get iapbs work
together with CHARMM on a 64-bit machine we had so I ended up only doing
the force field energy part with CHARMM and I had to modify the code to
write out PQR files. I then fed those PQR files "by hand" to standalone
apbs.
Another problem is that you will have to redo the topology files with
CHARMM. CHARMM assumes the PDB file in "standard" PDB order but
psfgen, at least in the version 1.4 I was using, did some odd swapping
of patched residues (disulfide bridges, C-terminal cap). This means
that the topology and the original NAMD DCD where out of sync with
respect to atom order. So I also had to reorder the DCD with catdcd
(part of VMD).
Hope that helps,
Hannes.
On Wed, 23 Oct 2013 11:49:40 +0530
Revthi Sanker <revthi.sanker1990_at_gmail.com> wrote:
> Dear all,
> Thank you so much for your suggestions. I am using the same dcd files
> while I am still working on the topology file conversion. Meanwhile I
> found this script:
>
> http://www.stfc.ac.uk/cse/25247.aspx
>
> ** mmpbsa.inp: run MM-PBSA analysis on trajectories*
> * only non-entropic contributions are computed
> *-
> * IMPORTANT: DEFInes must be modified (after PSF read) to match your
> system *-
> * usage: charmm mol=(complex|receptor|ligand) data=file < mm-pbsa.inp
> *-
> * written in 2010 by Hannes H. Loeffler, STFC Daresbury, UK
> *
>
>
> bomlev -2
> wrnlev 0
> prnlev 2
>
>
> set usage "charmm mol=(complex|receptor|ligand) data=file <
> mm-pbsa.inp"
>
> if @?mol eq 0 then
> echo "Error: mol parameter not set"
> echo @usage
> stop
> endif
>
> if @mol ne complex then
> if @mol ne receptor then
> if @mol ne ligand then
> echo "Error: mol parameter not one of complex, receptor or
> ligand" echo @usage
> stop
> endif
> endif
> endif
>
> if @?data eq 0 then
> echo "Error: no data file specified"
> echo @usage
> stop
> endif
>
>
> !== variables to be set by user ==
> ! modify the DEFINEs below too!
>
> set traj "protein.dcd" ! input trajectory
> set fframe 1 ! first frame
> set lframe 0 ! last frame: if < 0 then last frame from
> trajectory
> set skip 5
>
> set psf "complex.psf"
> set crd "complex.crd"
>
> set T 300.0 ! temperature
> set grdx 0.5 ! grid spacings for APBS
> set grdy 0.5
> set grdz 0.5
>
> ! ion parameters
> set ionq1 0 ! ion charge
> set ionc1 0.0 ! concentration in mol/l
> set ionr1 0.0 ! ion radius in Angstrom
> set ionq2 0
> set ionc2 0.0
> set ionr2 0.0
>
> set epsprot 1.0 ! dielectric constant of the protein
> set epssol 80.0 ! dielectric constant of the bulk solvent
> set epsvac 1.0 ! dielectric constant of the vacuum
> set rsol 1.4 ! solvent probe radius
>
> set gammaP 0.00542 ! for Gasa (non-polar dGsolv)
> set betaP 0.92 ! for Gasa
>
>
> !== read in the topology and parameter files ==
> read rtf card name $CHM_HOME/toppar/top_all27_prot_na.rtf
> read param card name $CHM_HOME/toppar/par_all27_prot_na.prm
>
> !== load the structure file ==
> open read unit 2 form name @psf
> read psf card unit 2
> close unit 2
>
> !!! DEFIne here the parts of the system !!!
> define complex select segid A .or. segid B end
> define receptor select segid A end
> define ligand select segid B end
>
> !== read in initial coordinate set ==
> open read unit 2 form name @crd
> read coor card unit 2
> close unit 2
>
> ! copy coordinates to comparison set
> coor copy comp
>
>
> !== setup trajectory for reading: we assume single trajectory for the
> moment ==
> set trajU 20
>
> ! determine last frame
> open read unit @trajU unform name @traj
> traj query unit @trajU
> close unit @trajU
>
> calc trlen = ?nfile - @fframe
>
> if @lframe le 0 then
> calc lframe = @fframe + @skip * int( @trlen / @skip )
> endif
>
> if @lframe le @fframe then
> echo "Error: last frame must be larger than first frame"
> stop
> endif
>
> open read unit @trajU unform name @traj
> traj first @trajU nunits 1 begin @fframe stop @lframe skip @skip
>
>
> !== initialise variables ==
> set datU 15
>
> set n 0 ! frame counter
> calc nmax = int( ( @lframe - @fframe ) / @skip )
>
> ! MM-PBSA equation: G = Egas + Gsolv - TSmm
> set avgEgas 0 ! solute internal energy
> set varEgas 0
> set avgEele 0
> set varEele 0
> set avgEvdW 0
> set varEvdW 0
> set avgEint 0
> set varEint 0
> set avgGasa 0 ! nonpolar contribution to solvation free energy
> set varGasa 0
> set avgGPB 0 ! solvation free energy from finite difference PBEQ
> set varGPB 0
>
> ! data file
> open write card unit @datU name @data
> write title unit @datU
> *#frame Egas Eele EvdW Eint Gasa
> GPB *
>
> !== BEGIN reading the trajectory ==
> label trajloop
> incr n by 1
>
> echo " MM-PBSA> step " @n
>
> traj read
>
> delete atoms select .not. @mol end
>
> coor orient norot ! center coordinates
>
> ! Egas
> energy cutnb 999.0 ctofnb 997.0 ctonnb 995.0 ctexnb 999.0
>
> set Egas = ?ener
> calc delta = @Egas - @avgEgas
> calc avgEgas = @avgEgas + @delta / @n
> calc varEgas = @varEgas + @delta * ( @Egas - @avgEgas )
>
> set Eele = ?elec
> calc delta = @Eele - @avgEele
> calc avgEele = @avgEele + @delta / @n
> calc varEele = @varEele + @delta * ( @Eele - @avgEele )
>
> set Evdw = ?vdw
> calc delta = @Evdw - @avgEvdw
> calc avgEvdw = @avgEvdw + @delta / @n
> calc varEvdw = @varEvdw + @delta * ( @Evdw - @avgEvdw )
>
> calc Eint = ?ener - ?elec - ?vdw
> calc delta = @Eint - @avgEint
> calc avgEint = @avgEint + @delta / @n
> calc varEint = @varEint + @delta * ( @Eint - @avgEint )
>
> ! Gasa
> coor surface rprobe @rsol
> calc Gasa = @gammaP * ?area + @betaP
> calc delta = @Gasa - @avgGasa
> calc avgGasa = @avgGasa + @delta / @n
> calc varGasa = @varGasa + @delta * ( @Gasa - @avgGasa )
>
> != GPB finite difference PBEQ =
> ! establish grid dimensions
> coor stat
>
> != BEGIN Poisson-Boltzmann equation =
> pbeq
> !scalar wmain = radius ! set wmain to vdW radii
> stream radius.str ! Benoit Roux's parameter
>
> apbs mgauto lpbe grdx @grdx grdy @grdy grdz @grdz -
> sdie @epssol pdie @epsprot srad @rsol temp @T -
> ionq1 @ionq1 ionc1 @ionc1 ionr1 @ionr1 -
> ionq2 @ionq2 ionc2 @ionc2 ionr2 @ionr2 -
> srfm 2 bcfl 1 swin 0.3 chgm 1 cmeth 1 calcene 1 debug 0
> set pbext = ?enpb
>
> apbs mgauto lpbe grdx @grdx grdy @grdy grdz @grdz -
> sdie @epsvac pdie @epsprot srad @rsol temp @T -
> srfm 2 bcfl 1 swin 0.3 chgm 1 cmeth 1 calcene 1 debug 0
> set pbvac = ?enpb
>
> calc GPB = @pbext - @pbvac
> calc delta = @GPB - @avgGPB
> calc avgGPB = @avgGPB + @delta / @n
> calc varGPB = @varGPB + @delta * ( @GPB - @avgGPB )
>
> reset
> end
> != END Poisson-Bolztmann equation =
>
> ! write data
> trim n from 1 to 10
> trim Egas from 1 to 10
> trim Eele from 1 to 10
> trim EvdW from 1 to 10
> trim Eint from 1 to 10
> trim Gasa from 1 to 10
> trim GPB from 1 to 10
>
> write title unit @datU
> *@n @Egas @Eele @EvdW @Eint @Gasa @GPB
> *
>
> ! reload the structure file
> open read unit 2 form name @psf
> read psf card unit 2
> close unit 2
> if @n lt @nmax goto trajloop
> !== END reading the trajectory ==
>
> !== calculate average values ==
> trim avgEgas from 1 to 12
> trim avgEele from 1 to 12
> trim avgEvdW from 1 to 12
> trim avgEint from 1 to 12
> trim avgGPB from 1 to 12
> trim avgGasa from 1 to 12
>
> calc n = @n - 1
>
> if @n gt 1 then
> calc varEgas = sqrt( @varEgas / @n )
> calc varEele = sqrt( @varEele / @n )
> calc varEvdW = sqrt( @varEvdW / @n )
> calc varEint = sqrt( @varEint / @n )
> calc varGPB = sqrt( @varGPB / @n )
> calc varGasa = sqrt( @varGasa / @n )
> endif
>
> trim varEgas from 1 to 12
> trim varEele from 1 to 12
> trim varEvdW from 1 to 12
> trim varEint from 1 to 12
> trim varGPB from 1 to 12
> trim varGasa from 1 to 12
>
> calc G = @avgEgas + @avgGPB + @avgGasa
>
> incr n by 1
>
> write title unit @datU
> *#
> *# MM-PBSA results for @mol (@n frames)
> *#
> *# MEAN STD
> *# ============ ============
> *# Eele @avgEele @varEele
> *# EvdW @avgEvdW @varEvdW
> *# Eint @avgEint @varEint
> *# Egas @avgEgas @varEgas
> *# Gasa @avgGasa @varGasa
> *# GPB @avgGPB @varGPB
> *#
> *# G(no_entropy) = Egas + Gasa + GPB = @G
> *#
> *#
> *# input files
> *# trajectory: @traj (first frame #@fframe, last frame #@lframe,
> skip @skip)
> *# PSF: @psf
> *# CRD: @crd
> *#
> *# key parameters
> *# T = @T K
> *# Gasa: gammaP = @gammaP, betaP = @betaP
> *# PBEQ: grid spacings = (@grdx, @grdy, @grdz),
> *# eps_Solv = @epssol, eps_Vac = @epsvac, eps_Prot = @epsprot,
> *# ion1_qcr = (@ionq1, @ionc1, @ionr1), ion2_qcr = (@ionq2,
> @ionc2, @ionr2),
> *# rsol = @rsol
> *#
> *
>
> close unit @datU
>
> stop
> >
-- Scanned by iCritical.
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:23:52 CST