CHARMM Development Project Forums CHARMM Community Script Archive MM-PBSA script

* 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

hallo..

some clarifications and questions about this script

in the mm-pbsa script...let y = maxi @y ?ymaxset z ?zmin, -- seems ?ymaxset is not a substitution parameter. is this ?ymax set z ?zmin ?
------------------------------------------------------

the deleted water residues (delete atom sele resn tip3 end) is causing atom miss match with the trajectory loop

(output file)
#######################################################
CHARMM>

CHARMM> coor copy comp
SELECTED COORDINATES COPIED TO THE COMPARISON SET.


CHARMM>

CHARMM> delete atom sele resn tip3 end
SELRPN> 55719 atoms have been selected out of 60231

Message from MAPIC: Atom numbers are changed.

Message from MAPIC: 18573 residues deleted.

Message from MAPIC: 2 segments deleted.
DELTIC:55719 bonds deleted
DELTIC:18573 angles deleted
DELTIC:18573 acceptors deleted
PSFSUM> PSF modified: NONBOND lists and IMAGE atoms cleared.
PSFSUM> Summary of the structure file counters :
Number of segments = 5 Number of residues = 301
Number of atoms = 4512 Number of groups = 1276
Number of bonds = 4503 Number of angles = 8166
Number of dihedrals = 12012 Number of impropers = 598
Number of HB acceptors = 479 Number of HB donors = 473
Number of NB exclusions = 0 Total charge = -0.00000

CHARMM> define solute sele .not. resn tip3 end
SELRPN> 4512 atoms have been selected out of 4512

CHARMM>

CHARMM> open unit 51 read unform name /big/konidala/mpijobs/mpijob1.trj
VOPEN> Attempting to open::/big/konidala/mpijobs/mpijob1.trj::
OPNLGU> Unit 51 opened for READONLY access to /big/konidala/mpijobs/mpijob1.trj

CHARMM> open unit 52 read unform name /big/konidala/mpijobs/mpijob1_1ns.trj
VOPEN> Attempting to open::/big/konidala/mpijobs/mpijob1_1ns.trj::
OPNLGU> Unit 52 opened for READONLY access to /big/konidala/mpijobs/mpijob1_1ns.trj

CHARMM>

CHARMM> ! Get trajectory file parameters
CHARMM> !traj iread 51
CHARMM>

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

CHARMM> traj iread 51 nunits 2 skip 500
TRAJ: INITIATING READ OF A TRAJECTORY, OPTIONS;
FIRSTU = 51 NUNIT = 2 SKIP = 500

CHARMM>

CHARMM> set a 10
Parameter: A <- "10"

CHARMM> set b 10
Parameter: B <- "10"

CHARMM> !calc z = 300 * 2
CHARMM>

CHARMM> set dcel 0.40 ! Grid spacing variable for PB equaition solution
Parameter: DCEL <- "0.40"

CHARMM>

CHARMM> ! MM-PBSA equiation
CHARMM> ! G = + + - TdS
CHARMM> set emm 0 ! Solute internal energy
Parameter: EMM <- "0"

CHARMM> set inte 0 ! ENER - VDW - ELEC
Parameter: INTE <- "0"

CHARMM> set exte 0 ! VDW + ELEC
Parameter: EXTE <- "0"

CHARMM> set dgpb 0 ! solvation free energy from finite diffrence PB eq
Parameter: DGPB <- "0"

CHARMM> set dgasa 0 ! nonpolar contribution to solvation free energy
Parameter: DGASA <- "0"

CHARMM> set n 0 ! snapshot counter
Parameter: N <- "0"

CHARMM>

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

NONBOND OPTION FLAGS:
ELEC VDW ATOMs CDIElec SHIFt VATOm VSWItch
BYGRoup NOEXtnd NOEWald
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.000
WMIN = 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS = 1.000
NBXMOD = 5
There are 0 atom pairs and 0 atom exclusions.
There are 0 group pairs and 0 group exclusions.
<MAKINB> with mode 5 found 12669 exclusions and 11780 interactions(1-4)
<MAKGRP> found 3459 group exclusions.
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR 1296530 ATOM PAIRS AND 0 GROUP PAIRS

General atom nonbond list generation found:
1207447 ATOM PAIRS WERE FOUND FOR ATOM LIST
67169 GROUP PAIRS REQUIRED ATOM SEARCHES


CHARMM>

CHARMM> label trajloop

CHARMM> traj read

READING TRAJECTORY FROM UNIT 51
NUMBER OF COORDINATE SETS IN FILE: 1120
NUMBER OF PREVIOUS DYNAMICS STEPS: 40500
FREQUENCY FOR SAVING COORDINATES: 500
NUMBER OF STEPS FOR CREATION RUN: 560000

***** ERROR ***** NO. OF ATOMS DO NOT MATCH 4512 60231
*** LEVEL -4 WARNING *** BOMLEV IS -2
BOMLEV HAS BEEN SATISFIED. TERMINATING.

########################################################
---------------------------------------------------------

commenting the line.. delete atom sele resn tip3 end in the script and considering the whole system, the job terminates with the warning from VEHEAP storage request. eventhough i am using a higher heap size of 40480000.


########################################################
CHARMM>

CHARMM> pbeq

Calculations with the Poisson-Boltzmann Equation


PBEQ> scalar wmain = radius

PBEQ>

PBEQ> solve dcel @dcel nclx @xcel ncly @ycel nclz @zcel maxit 1000 -
PBEQ> epsw 80.0 epsp 4.0 temp 298 watr 1.4 intbp
Parameter: DCEL -> "0.40"
Parameter: XCEL -> "252"
Parameter: YCEL -> "252"
Parameter: ZCEL -> "220"

Calculation with 60231 atoms

LINEARIZED PBEQ SOLVER: Successive OverRelaxation (SOR) method

ITERATION PARAMETERs
Maximum # iterations (MAXITS) = 1000
Tolerance of convergence (DEPS) =.200E-05
Mixing factor (LAMBDA,DOME) = 1.000

CHARGE DISTRIBUTION METHOD: the trilinear interpolation

BOUNDARY POTENTIAL CALCULATION METHOD
The Debye-Huckel approximation for half number of boundary points along 1d-axis and
potential of the rest will be interpolated from nearest grid points

PHYSICAL PARAMETERs
Solvent probe radius (WATR) = 1.400 [Angs]
Ion exclusion radius (Stern layer) = 0.000 [Angs]
Solvent dielectric constant (EPSW) = 80.000
Protein dielectric constant (EPSP) = 4.000
Salt concentration (CONC) = 0.000 [moles]/[liter]
Temperature (TEMP) = 298.000 [K]
Debye-Huckel factor (KAPPA2) = 0.000 [1/Angs**2]

NUMBER OF GRID POINTS: 253 253 221
Box in X from -50.400 to 50.400
Box in Y from -50.400 to 50.400
Box in Z from -44.000 to 44.000
VEHEAP> Expanding heap size by14155776 words.
VEHEAP> Expanding heap size by14155776 words.
VEHEAP> Error: MALLOC returned ZERO.

***** LEVEL -4 WARNING FROM <VEHEAP> *****
***** COULD NOT SATISFY STORAGE REQUEST.
******************************************
BOMLEV ( -2) IS REACHED - TERMINATING. WRNLEV IS 5

###########################################################
-----------------------------------------------------------

will this PBEQ method work for randomly distributed lipids and ions around the protein without deleting them for the free energy calculations?


thanks for your comments in advance.

cheers

I have a question about the script of PBEQ

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

The first one is for protein in solvent surrounding. so we set the parameters of epsw is 80.0 for water molecule
and epsp is 4.0 for protein inside, water probe radius is 1.4 angs.

and the second one is for protein in vacuum surrounding. we set the parameters of epsw is 1.0 for vacuum. Do we need
to set the probe radius is 1.4 angstr?
Or maybe we can set the parameter of watr is 0.0?


Thanks~