From: Kenno Vanommeslaeghe (kvanomme_at_rx.umaryland.edu)
Date: Fri Oct 25 2013 - 15:05:05 CDT
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
>
>>>
>
This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:49 CST