From: Revthi Sanker (revthi.sanker1990_at_gmail.com)
Date: Wed Oct 23 2013 - 01:19:40 CDT
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
Has anybody used this script or tell me if this will meet my objective?
Kindly provide your valuable suggestions regarding the same.
Thank you so much.
Yours sincerely.
Revathi.S
M.S. Research Scholar
Indian Institute Of Technology, Madras
India
______________________________________________
International Conference on Bimolecular Simulations and Dynamics
*Official website:*
http://cmsrv.iitm.ac.in/icbsd2013/
_______________________________________________
On Tue, Oct 22, 2013 at 10:26 PM, Kenno Vanommeslaeghe <
kvanomme_at_rx.umaryland.edu> wrote:
> Let me start with trying to answer the last question in Revathi e-mail:
> no, I don't know of any way to accurately calculate a binding free energy
> without having to run simulations for protein (P), ligand (L) and
> protein-ligand (PL) complex separately.
>
> It is hard to predict what kind of errors are introduced by
> post-processing a CHARMM FF trajectory with the AMBER FF like this. I guess
> the worst errors will be canceled out by running all 3 simulations (PL, P,
> L) with the same CHARMM-based methodology prior to switching to AMBER. The
> remaining errors could be bad, or they could be benign compared to the
> approximation of using implicit solvent; I can't easily tell. The problem
> could be avoided altogether by using the CHARMM FF with the AMBER software
> (which is possible); however, the radii used to generate the cavity may be
> optimized for AMBER, which introduces an error that may or may not be
> greater than the one discussed above.
>
> If you have access to the CHARMM software, you can do the MM/PBSA analysis
> and be free of energy-related as well as file conversion issues, which
> would seem to be preferable.
>
>
>
> On 10/22/2013 06:01 AM, Jason Swails wrote:
>
>>
>>
>>
>> On Tue, Oct 22, 2013 at 5:21 AM, Peter Jones <pm-jones_at_bigpond.com
>> <mailto:pm-jones_at_bigpond.com>> wrote:
>>
>> Dear Rethvi,
>>
>> I did this by converting the NAMD-generated dcd trajectory to AMBER
>> format using the software Simulaid. Then you just prepare the
>> appropriate AMBER topology and parameter files, and use the AMBER
>> MM/PBSA protocols to analyse your trajectory. There might be issues
>> here regarding artefacts or inaccuracies due to differences in the
>> forcefields, I couldn't comment on that, but the results I got made
>> very good sense. Converting the trajectory was no laughing matter
>> either, and whichever way you go you're probably going to have to just
>> dig into the errors until you get rid of them,
>>
>>
>> For what it's worth, the AmberTools program that performs MM/PBSA analyses
>> can read DCD files natively, so there's no need to convert them. The
>> tricky part is getting the topology file if you are starting with
>> CHARMM-based files.
>>
>> Hope this helps,
>> Jason
>>
>> --
>> Jason M. Swails
>> BioMaPS,
>> Rutgers University
>> Postdoctoral Researcher
>>
>
>
This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:48 CST