Re: [apbs-users] Titration Curves using the single site approximation

[Murga, L. <l....@no...>]

Just in case, I am including here the UHBD script that we would use in this case:

read mol 1 file "pkaSH.pdb"        pdb end  ====> coordinates of the protein
read mol 2 file "sitesinpr.pdb"    pdb end  =====> coordinates of the atoms of ionizable residues
set charge radii file "pkaS.dat"         para mine end

 edit all charge 0. end
 edit charge 1. atnum     1 end

read grid epsi binary file  "coarse.epsi" end
read grid epsj binary file  "coarse.epsj" end
read grid epsk binary file  "coarse.epsk" end

 elec calc mol 1
 bcflag 2 solver 0
 reps
 nmap      1.4
 nsph     500
 maxits     300
 temp    293.0
 sdie   80.00
 pdie   20.00
 ionstr   150.000
 rion    2.00
end

print elec phizero all end
print elec phisave all end

 elec calc mol 1
 bcflag 4 solver 0
 spacing   1.20 dime    15    15    15

 gcenter   24.704  20.926  27.944
 nmap      1.4
 nsph     500
 maxits     300
 temp    293.0
 sdie   80.00
 pdie   20.00
 ionstr   150.000
 rion    2.00
end
print elec phisave all end

 elec calc mol 1
 bcflag 4 solver 0
 spacing   0.75 dime    15    15    15

 gcenter   24.704  20.926  27.944
 nmap      1.4
 nsph     500
 samsrf
 maxits     300
 temp    293.0
 sdie   80.00
 pdie   20.00
 ionstr   150.000
 rion    2.00
end
print elec phisave all end

 elec calc mol 1
 bcflag 4 solver 0
 spacing   0.25 dime    20    20    20

 gcenter   24.704  20.926  27.944
 nmap      1.4
 nsph     500
 samsrf
 maxits     300
 temp    293.0
 sdie   80.00
 pdie   20.00
 ionstr   150.000
 rion    2.00
end

Thanks!


Leonel

________________________________
From: Murga, Leonel <l....@no...>
Sent: Tuesday, June 27, 2017 3:46:21 PM
To: apb...@li...
Subject: [apbs-users] Titration Curves using the single site approximation


Hello all,


For years our group has been using UHBD to calculate titration curves using the single site approximation. Under this approximation, a +1 charge is located in a particular atom of an ionizable group of a protein (for example at the C of the COOH group of an aspartate or the N of the Lysine) while all the other charges are set to zero; then, the potentials that this charge creates at the other ionizable groups in the protein are calculated via the PB equation. These are sometimes called interaction potential. The process is carried on until all the interaction potentials created by placing at +1 charge on each ionizable group of the protein have been calculated. The potentials are then used for the calculation of titration curves using one of several methodologies available.


For large proteins the use of a fine enough grid (~0.25 Ang) is not feasible because of the computational cost it represents. Even more so in the above situation in which the PBE has to be solved once for each ionizable residue in the protein. For these cases UHBD uses focusing. In essence, the PB equation is solved once using a coarse grid (~2.0 to 2.5 Ang) centered at the protein center in an initial step. In additional calculations, increasingly finer grids are placed and centered at the atom where the +1 charge was located and the PB equation is solved. This process is repeated until the finest grid is reached. The electrostatic potentials are then calculated at each atomic position in the protein including where the +1 charge is located.


We have tried to implement these algorithm in APBS without success. Here is a section of the input script for APBS for calculating the interaction potentials for the first ionizable group in a protein (this would actually corresponds to the N atom of the N-terminus):

read
    mol pqr pkaSH.pqr ====> protein in its neutral state
    mol pqr pkaSH_Mod.pqr =====> protein with all charges set to zero except a +1 charge on the first ionizable
    mol pqr chargeaa_Mod.pqr  =======> this is the isolated ionizable used in a later calculation that is not included here.
end
elec name prot-coarse
    mg-manual
    dime 129 129 129
    grid 1.00 1.00 1.00
    gcent mol 2
    mol 2
    lpbe
    bcfl sdh
    ion charge 1 conc 0.150 radius 2.0
    ion charge -1 conc 0.150 radius 2.0
    pdie 20.0000
    sdie 80.0000
    srfm smol
    chgm spl2
    sdens 10.00
    srad 1.40
    swin 0.30
    temp 298.15
    calcforce no
end
elec name prot-fine1
    mg-manual
    dime 97 97 97
    grid 0.5 0.5 0.5
    gcent 24.704 20.926 27.944  ======> These are the coordinates of the N atom of the first ionizable group
    mol 2
    lpbe
    bcfl focus
    ion charge 1 conc 0.150 radius 2.0
    ion charge -1 conc 0.150 radius 2.0
    pdie 20.0000
    sdie 80.0000
    srfm smol
    chgm spl2
    sdens 10.00
    srad 1.40
    swin 0.30
    temp 298.15
    calcforce no
end
elec name prot-fine2
    mg-manual
    dime 65 65 65
    grid 0.25 0.25 0.25
    gcent 24.704 20.926 27.944
    mol 2
    lpbe
    bcfl focus
    ion charge 1 conc 0.150 radius 2.0
    ion charge -1 conc 0.150 radius 2.0
    pdie 20.0000
    sdie 80.0000
    srfm smol
    chgm spl2
    sdens 10.00
    srad 1.40
    swin 0.30
    temp 298.15
    calcenergy total
    calcforce no
    write pot dx pot_mulgr_pr
end
The pot_mulgr_pr file is written without any problems.
We then use the multivalue program in the utils folder asking it to calculate the potentials at all the atomic positions and here is what we get (this is only a small section):

2.622700e+01,1.271200e+01,2.357900e+01,-1.649562e+18
2.608400e+01,1.166200e+01,2.258300e+01,NaN
2.535400e+01,1.046600e+01,2.321300e+01,NaN
2.533700e+01,1.031000e+01,2.443700e+01,NaN
2.711200e+01,1.272300e+01,2.407800e+01,-1.460145e+18
2.480500e+01,9.537000e+00,2.245400e+01,NaN
2.405200e+01,8.386000e+00,2.297400e+01,NaN
2.489700e+01,7.502000e+00,2.384900e+01,NaN
2.450400e+01,7.220000e+00,2.497200e+01,NaN
2.350600e+01,7.549000e+00,2.179300e+01,NaN
2.274100e+01,6.293000e+00,2.218200e+01,NaN
2.224200e+01,5.559000e+00,2.093100e+01,NaN
2.331900e+01,5.176000e+00,2.003700e+01,NaN
2.393100e+01,3.984000e+00,2.008300e+01,NaN
2.362200e+01,3.034000e+00,2.096100e+01,NaN
2.489500e+01,3.751000e+00,1.919900e+01,NaN
2.475500e+01,9.746000e+00,2.147200e+01,NaN
2.345400e+01,5.790000e+00,1.926600e+01,NaN
2.296200e+01,3.182000e+00,2.172500e+01,NaN
2.411300e+01,2.152000e+00,2.095500e+01,NaN
2.535300e+01,2.849000e+00,1.918600e+01,NaN

Almost every potential comes as NaN.


In UHBD this process works without any problems.

What am I doing wrong?

Thanks for your help!

Leonel Murga