* mmpbsa.inp : run MM-PBSA analysis on dynamics trajectories
* Dynamics were run in explicit solvent. A desolvated trajectory
* was produced with the merge coor command (as used below)
* Entropy is calculated from quasiharmonic analysis of heavy atoms.
*

bomlev -2

set smi 3
set seg c
set mod b
set cycle 1
set runs 1
set dcd dyna/dcd

set 1 ~wardjm/charmm ! Directory holding rtf and prm files
set 2 @1/martin ! Working directory
set crd @2/crd ! Directory holding crd files
set dcd @2/@dcd ! Directory holding trajectory files

! Read in the topology and parameter files
open read unit 1 form name @1/newtop.rtf
read rtf card unit 1
close unit 1

open read unit 1 form name @1/newpar.prm
read para card unit 1
close unit 1

! Load the psf: protein plus solvent box
open read unit 2 form name @2/@smi.@seg.@mod.solv.psf
read psf card unit 2
close unit 2

! Read in initial coordinate set
open read unit 2 form name @crd/@smi.@seg.@mod.eq2.crd
read coor card unit 2
close unit 2

coor copy comp

delete atom sele resn tip3 end
define solute sele .not. resn tip3 end

! Open a unit to hold the trajectory file(s) -- # of files = @runs
set open 20

label openloop
open read unit @open unform name -
@dcd/@smi.@seg.@mod.@cycle.solu.dcd

incr cycle by 1
incr open by 1

if @cycle le @runs goto openloop

! Get trajectory file parameters
traj iread 20
traj readset begin ?start ! My file = 100100
set skip ?skip ! = 100
set nfile ?nfile ! = 3000
calc nframes = @nfile * @runs ! total number of snapshots in analysis

! MM-PBSA Analysis
! For Emm, dGpb, dGasa calculation, analyze 10 ps snapshots
! For quasiharmonic analysis, use full trajectory with 100 fs skip

traj iread 20 nunits @runs begin 110000 skip 10000

set a 10
set b 10
calc z = 300 * @runs

set dcel 0.40 ! Grid spacing variable for PB equaition solution

! MM-PBSA equiation
! G = + + - TdS
set emm 0 ! Solute internal energy
set inte 0 ! ENER - VDW - ELEC
set exte 0 ! VDW + ELEC
set dgpb 0 ! solvation free energy from finite diffrence PB eq
set dgasa 0 ! nonpolar contribution to solvation free energy
set n 0 ! snapshot counter

nbonds ! set up nonbond list (defaults from parameter set used)

label trajloop
traj read

! Emm
energy
calc emm = @emm + ?ener
calc exte = @exte + ?elec + ?vdw
calc inte = @inte + ( ?ener - ?elec - ?vdw )

! Gasa Gasa = 0.00542 * ASA + 0.92
coor surf rprobe 1.4
calc dgasa = @dgasa + ( 0.00542 * ?area + 0.92 )

! Gpb Finite difference PB eq
! Establish grid dimensions
coor stat

set x ?xmin
let x = abs @x
let x = maxi @x ?xmax
calc x = 2 * ( 5 + @x )

set y ?ymin
let y = abs @y
let y = maxi @y ?ymaxset z ?zmin
calc y = 2 * ( 5 + @y )

set z ?zmin
let z = abs @z
let z = maxi @z ?zmax
calc z = 2 * ( 5 + @z )

calc x = @x / @dcel
calc xcel = int( @x )
if @xcel lt @x incr xcel by 1
incr xcel by 1

calc y = @y / @dcel
calc ycel = int( @y )
if @ycel lt @y incr ycel by 1
incr ycel by 1

calc z = @z / @dcel
calc zcel = int( @z )
if @zcel lt @z incr zcel by 1
incr zcel by 1


pbeq
scalar wmain = radius

solve dcel @dcel nclx @xcel ncly @ycel nclz @zcel maxit 1000 -
epsw 80.0 epsp 4.0 temp 298 watr 1.4 intbp
set ener80 = ?enpb

solve dcel @dcel nclx @xcel ncly @ycel nclz @zcel maxit 1000 -
epsw 1.0 epsp 4.0 temp 298 watr 1.4 intbp
set ener1 = ?enpb

calc dgpb = @dgpb + (@ener80 - @ener1)

reset
end

incr n by 1
incr a by @b
if @a le @z goto trajloop

! Find ensemble average values
calc avgemm = @emm / @n
calc avginte = @inte / @n
calc avgexte = @exte / @n
calc avgdgpb = @dgpb / @n
calc avgdgasa = @dgasa / @n

! Generate desolvated trajectory w/o hydrogens for quasiharmonics
! orient wrt mass weighted rms
open write unit 10 unform name @dcd/@smi.@seg.@mod.-h.dcd

merge coor firstu 20 nunit @runs outputu 10 -
begin @begin skip @skip nfile @nframes -
sele .not. hydrogen end -
orie mass sele .not. hydrogen end

open read unit 20 unform name @dcd/@smi.@seg.@mod.-h.dcd

delete atom sele hydrogen end

! Calculate average structure and copy to comparison set
coor dyna firstu 20 nunit 1 begin @begin skip @skip -
sele .not. hydrogen end orie mass sele .not. hydrogen end

coor copy comp

! Calculate number of vectors to use for vibrational analysis
set nmod ?natom
calc nmod = @nmod * 3

vibran nmod @nmod

quasi temp 298 firstu 20 nunit 1 begin @begin skip @skip

thermo temp 298
set stot ?stot

end

calc tds = 298 * @stot
calc dG = @avgemm + @avgdgpb + @avgdgasa - @tds

open write unit 15 form name @2/dat/@smi.@seg.@mod.@cycle.mmpbsa.dat
write title unit 15
* MM-PBSA results for @smi.@seg.@mod
*************************************************************
* Emm: Total @emm Avg @avgemm
* Internal @inte @avginte
* External @exte @avgexte
* dgpb: Total @dgpb Avg @avgdgpb
* dgasa: Total @dgasa Avg @avgdgasa
* stot = @stot Tds = @tds
**************************
* dG = @dG
*

stop



Joshua Ward Graduate Student Purdue University Department of Medicinal Chemistry and Pharmacology