PBEQ Solver Tutorial

Electrostatic potential and pKa calculations with the PBEQ module

pbeq.inp This input is under the "epot" directory.
Once we read a protein structure, a generic stream file "pbeq.str" is called to calculate the electrostatic potential and solvation energy of the protein. One can choose different options in the stream file:

set EpsP    =   1.0      ! dielectric constant for the protein interior
set EpsW    = 80.0    ! solvent dielectric constant
set Conc    =   0.0      ! salt concentration (M)
set Focus   = true     ! to have a refined calculation focussed
! on the site using a finer grid
set Dcel_c =   1.2     ! the grid spacing in the finite-difference
! (centered on Xcen,Ycen,Zcen)
set Dcel_f =   0.4      ! the grid spacing in the finite-difference
! (centered on Xcen,Ycen,Zcen)
set LEdge   = 10.0    ! distance between a protein atom and a grid
                         ! LEdge*2 for coarse-gird calculations and
! LEdge/2 for fine-grid calculations (see below)
set Options = watr 1.4 reentrant ! Let's use the molecular surface

Then, the center of the protein and number of grid points along XYZ are determined as

coor stat
coor orient norotate
coor stat
set Xcen = ?xave
set Ycen = ?yave
set Zcen = ?zave

calc Nclx_c = int ( ( @LEdge * 4.0 + ?Xmax - ?Xmin ) / @{Dcel_c} )
calc Ncly_c = int ( ( @LEdge * 4.0 + ?Ymax - ?Ymin ) / @{Dcel_c} )
calc Nclz_c = int ( ( @LEdge * 4.0 + ?Zmax - ?Zmin ) / @{Dcel_c} )

if Focus eq true then
   calc Nclx_f = int ( ( @LEdge * 1.0 + ?Xmax - ?Xmin ) / @{Dcel_f} )
   calc Ncly_f = int ( ( @LEdge * 1.0 + ?Ymax - ?Ymin ) / @{Dcel_f} )
   calc Nclz_f = int ( ( @LEdge * 1.0 + ?Zmax - ?Zmin ) / @{Dcel_f} )
endif

To obtain a meaningful result, one have to use an optimized radii for each atom. This is done by

prnlev 0
stream radii_prot_na.str
prnlev 5
scalar wmain statistics select .not. type H* end
define check select (.not type H* ) .and. ( property wmain .eq. 0.0 ) end
if ?nsel ne 0 stop       !some heavy atom have a zero radius

Then, we perform 4 PB calculations: coarse and fine grid calculations for solvent and vacuum;

PBEQ

! in solution
SOLVE Nclx @{nclx_c} Ncly @{ncly_c} Nclz @{nclz_c} .......
if focus eq true SOLVE Nclx @{nclx_f} Ncly @{ncly_f} Nclz @{nclz_f} .....

set Ener80 = ?enpb

! write phi
open write file unit 40 name @pdb.phi80
write PHI kcal unit 40
close unit 40

write phi kcal card xfirst -100.0 xlast 100.0 unit 6

! in vacuum
SOLVE Nclx @{nclx_c} Ncly @{ncly_c} Nclz @{nclz_c} ......
if focus eq true SOLVE Nclx @{nclx_f} Ncly @{ncly_f} Nclz @{nclz_f} ......

set Ener1 = ?enpb

calc dEner = @Ener80 - @Ener1

RESET
END

The above calculations can be executed by

$CHARMMEXEC pdb=3gb1 < pbeq.inp > 3gb1.out

One can repeat the calculations with "pdb=1ppf".

There is a file called "view.com" which will visualize the calculated potential "@pdb.phi80" using the DINO visualization program.

pka_single.inp
pka_multiple.inp These inputs are under the "pka" directory. In this example, we will compute pKa values of Asp7, Glu10, Glu19, Asp27, and Glu43 of OMTKY3 using two PDB structures from NMR (PDB:1OMU) and X-ray (PDB:1PPF). We will examine the influence of protein interior dielectric constant as well as treatment of titration sites on the computed pKa. "pka_single.inp" consider one site at a time, while "pka_multiple.inp" treat the multiple (5 in this example) sites all together. When multiple sites are treated at the same time, all other sites are neutralized during the PB calculations of a specific site. The "single site" file is executed by

$CHARMMEXEC pdb=1ppf epsp=20 resid=7 < pka_multiple.inp > 1ppf_pka.out,resid7_epsp20

This will produce "1ppf_pka.dat,resid7_epsp20" which contains

site pKa_intr pKa_mod pKa_shift pKa_Born pKa_back dG_Born dG_back
1 3.023    3.86    -0.836    0.313 -1.150    -0.431    1.579

pKa_intr(insic) = pKa_mod(el) + pKa_shift
pKa_shift = pKa_Born + pKa_back

The "multiple sites" file is executed by

$CHARMMEXEC pdb=1omu epsp=20 < pka_multiple.inp > 1omu_pka.out,epsp20

This will produce "1omu_pka.dat,epsp20" which contains

site pKa_intr pKa_mod pKa_shift pKa_Born pKa_back dG_Born dG_back

1 2.683 3.86 -1.177 0.265    -1.442       -0.363 1.979
2 2.706 4.07 -1.363 0.292    -1.656    -0.401 2.273
3 2.102 4.07 -1.967 0.014    -1.982    -0.020 2.720
4 2.794 3.86 -1.065 0.585    -1.651       -0.804 2.266
5 3.489 4.07 -0.580 0.107 -0.688 -0.147 0.944

The input files call "pka_smeared.str" which sequentially calls "pka_charge.str", "pka_shift.str", and "pka_inter.str".

Experimental values are
Asp7 : < 2.7
Glu10 : 4.2
Glu19 : 3.2
Asp27 : < 2.3
Glu43 : 4.8