**From:** Kenno Vanommeslaeghe (*kvanomme_at_rx.umaryland.edu*)

**Date:** Fri Oct 25 2013 - 15:05:05 CDT

**Next message:**Roy Fernando: "Mysterious slow down in parallel"**Previous message:**Joseph Farran: "Re: Fwd: Re: [Checkpoint] BLCR and NAMD"**In reply to:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Next in thread:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Reply:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

Good stuff! For the psfgen issue, I think I might have a possible solution:

- CHARMM version 36b1 and newer appear to contain code to read XPLOR psf

files, so you can read the psfgen psf and thereby sidestep the atom

ordering issue.

- If you're stuck with an older CHARMM version, the attached shell script

should be able to convert an xplor psf into a CHARMM psf that can likewise

be read into CHARMM (though I haven't thoroughly tested it).

One question for Hannes: what was the reason for removing the

quasiharmonic analysis from the original script at

http://www.charmm.org/ubbthreads/ubbthreads.php?ubb=showflat&Number=1118

?

On 10/23/2013 03:39 AM, Hannes Loeffler wrote:

*> 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
*

*>
*

*>>>
*

*>
*

- text/plain attachment: psf_xplor2charmm

**Next message:**Roy Fernando: "Mysterious slow down in parallel"**Previous message:**Joseph Farran: "Re: Fwd: Re: [Checkpoint] BLCR and NAMD"**In reply to:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Next in thread:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Reply:**Hannes Loeffler: "Re: MMPBSA-like analysis"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

*
This archive was generated by hypermail 2.1.6
: Tue Dec 31 2013 - 23:23:53 CST
*