pKa Calculation Tutorial

pKa Calculation Example

(Modified from http://mccammon.ucsd.edu/~chem215/lectures/lecture10/param/node3.html )

Here are the pKa calculation procedures in step by step detail:

First make a subdirectory for this exercise as follows:
```
# mkdir pka 
# cd pka
```

Second copy (or download) the six necessary files into your working directory (copy will be much faster than download):

# cp /class/pka/* . (<--Be careful of the dot.) 
# ls -l 
-rwxr--r--    1 kwlee    kwlee         752 Apr 25 15:44 pkaS-addH*
-rw-r--r--    1 kwlee    kwlee       31126 Apr 25 15:26 pkaS.dat
-rw-r--r--    1 kwlee    kwlee         839 Apr 25 15:26 pkaS-doinp.inp
-rw-r--r--    1 kwlee    kwlee         361 Apr 25 15:26 pkaS-dops.inp
-rwxr--r--    1 kwlee    kwlee        1215 Apr 25 15:26 pkaS-dosbs*
-rw-r--r--    1 kwlee    kwlee        1477 Apr 25 15:26 pkaS-hbuild.inp

Download/prepare your pdb file for the calculation. I have already downloaded PhoE porin (1pho.pdb) file from the Brookhaven database for you. Copy or rename pkaS.pdb from your pdb file for the next step:
```
# cp 1pho.pdb pkaS.pdb 
```

Now let's edit the PDB file to fix it up for the pKa procedure. Remove all of the header infos, crystal waters,... and finally leave the ATOM CARD(or lines) in the pdb file. The pkaS.pdb is the file you want to make.

# head pkaS.pdb

ATOM      1  N   ALA     1      99.249  40.193  37.394  1.00 34.89
ATOM      2  CA  ALA     1     100.516  40.633  38.037  1.00 34.89
ATOM      3  C   ALA     1     101.298  41.520  37.084  1.00 34.89
ATOM      4  O   ALA     1     101.379  41.233  35.892  1.00 34.89
ATOM      5  CB  ALA     1     101.356  39.431  38.419  1.00 14.52
ATOM      6  N   GLU     2     101.843  42.617  37.600  1.00 35.20
ATOM      7  CA  GLU     2     102.617  43.532  36.772  1.00 35.20
ATOM      8  C   GLU     2     104.058  43.081  36.802  1.00 35.20
ATOM      9  O   GLU     2     104.803  43.428  37.707  1.00 35.20
ATOM     10  CB  GLU     2     102.501  44.962  37.289  1.00 69.96

# tail pkaS.pdb


ATOM   2618  C   PHE   340      97.666  37.450  39.007  1.00 26.37
ATOM   2619  O   PHE   340      96.598  37.524  38.364  1.00 26.37
ATOM   2620  CB  PHE   340      99.107  35.721  40.108  1.00 27.69
ATOM   2621  CG  PHE   340      99.460  35.178  38.760  1.00 27.69
ATOM   2622  CD1 PHE   340      99.039  33.920  38.374  1.00 27.69
ATOM   2623  CD2 PHE   340     100.242  35.916  37.882  1.00 27.69
ATOM   2624  CE1 PHE   340      99.398  33.402  37.126  1.00 27.69
ATOM   2625  CE2 PHE   340     100.601  35.408  36.643  1.00 27.69
ATOM   2626  CZ  PHE   340     100.180  34.150  36.264  1.00 27.69
ATOM   2627  OXT PHE   340      98.654  38.164  38.766  1.00 27.69

There are still some problems with this PDB file that you would encounter when you tried to add hydrogens. Some of the residues have information about alternate positions of some of the side chain atoms. One of those must be selected for each case and the residue name fixed, the other atom positions deleted, and the occupancy set to 1.00. See the below example:

...... lines not shown ......
ATOM    432  CG  LYS    53      -5.378  14.291  11.615  1.00 11.19   
ATOM    433  CD  LYS    53      -6.203  14.696  10.348  1.00 22.57   
ATOM    434  CE  LYS    53      -5.338  15.443   9.380  1.00 27.91  
ATOM    435  NZ  LYS   A53      -5.851  16.739   8.906  0.50 21.67 
ATOM    436  NZ  LYS   B53      -4.770  14.533   8.344  0.50 26.94 
ATOM    437  N   HIS    54      -3.582  16.632  13.804  1.00  6.17  
ATOM    438  CA  HIS    54      -2.193  17.084  13.582  1.00  5.11  
ATOM    439  C   HIS    54      -1.334  16.703  14.783  1.00  6.15   
...... lines not shown ......

...... lines not shown ......
ATOM    432  CG  LYS    53      -5.378  14.291  11.615  1.00 11.19  
ATOM    433  CD  LYS    53      -6.203  14.696  10.348  1.00 22.57  
ATOM    434  CE  LYS    53      -5.338  15.443   9.380  1.00 27.91  
ATOM    435  NZ  LYS    53      -5.851  16.739   8.906  1.00 21.67 
ATOM    437  N   HIS    54      -3.582  16.632  13.804  1.00  6.17  
ATOM    438  CA  HIS    54      -2.193  17.084  13.582  1.00  5.11  
ATOM    439  C   HIS    54      -1.334  16.703  14.783  1.00  6.15   
...... lines not shown ......

After you revised/deleted some atoms(lines) you should renumber the atoms in the pdb file so type the following commands:

% renumber

Enter name of PDB file to renumber ---> pkaS.pdb
Renumbered PDB file will appear in:pkaS.pdb.new
          948 Atoms Processed Successfully

% mv pkaS.pdb.new pkaS.pdb
overwrite pkaS.pdb? (yes/no)[no] : y

Edit the file again and put a HEADER line at the top. You can do this by any editor.


# head pkaS.pdb

HEADER
ATOM      1  N   ALA     1      99.249  40.193  37.394  1.00 34.89
ATOM      2  CA  ALA     1     100.516  40.633  38.037  1.00 34.89
ATOM      3  C   ALA     1     101.298  41.520  37.084  1.00 34.89
ATOM      4  O   ALA     1     101.379  41.233  35.892  1.00 34.89
ATOM      5  CB  ALA     1     101.356  39.431  38.419  1.00 14.52
ATOM      6  N   GLU     2     101.843  42.617  37.600  1.00 35.20
ATOM      7  CA  GLU     2     102.617  43.532  36.772  1.00 35.20
		...... lines not shown ......

Now look to see how many histidines there are and of what type. There is one and is of type HISA (HD1 present in neutral form).
```
# grep 'CA  HIS' pkaS.pdb
ATOM    154  CA  HIS    21      83.821  20.395  40.368  1.00 11.67
 
```
Now look to see how many cysteines there are and of what type (i.e. are they all involved in disulphide bonds). Use insightII or VMD to locate the CYS residues to see. If they are not involved in disulphide bonds, they need to be treated as ionizable.
```
# grep 'CA  CYS' pkaS.pdb
```

Now edit the file: pkaS-doinp.inp (vi pkaS-doinp.inp)...

          NAME of mol 1 file
pkaSH.pdb
          NAME of charge and radius file
pkaS.dat
          NUMBER of polypeptide chains of mol 1
1
          NUMBERS of first residues of each chain
1
          NUMBERS of last residues of each chain
340
          NO of histidins, sites: 1=ND1=HISB; 2=NE2=HISA; 0=CE1 in A or B
1
2
          NO GRIDS; SPACING and dime (max 5 grids)
4
2.5 45 45 45
1.2 15 15 15
0.75 15 15 15
0.25 20 20 20
          MAXIMAL number of iterations for elec
300
          TEMPERATURE in K
293.0
          DIELECTRIC constants of solvent and protein
80.0 20.0
          IONIC strength and radius of ions
150.0 2.0
          NO ADD TITR SITES, SITES DATA:a4,1x,a4,1x,f5.1,1x,i2
0
        CYS to be included, how many? and their res numbers
0
          DATA for new diel. map: nmap and nsph. Put nsph=0 for no_map
1.4 500

If you have two CYS residues not involved in disulphide bonds, you have to do like this:

...... lines not shown ......
          IONIC strength and radius of ions
150.0 2.0
          NO ADD TITR SITES, SITES DATA:a4,1x,a4,1x,f5.1,1x,i2
0
        CYS to be included, how many? and their res numbers
2
9 50
          DATA for new diel. map: nmap and nsph. Put nsph=0 for no_map
1.4 500

Now you are ready to run the pkaS-addH procedure to add hydrogens to all polar atoms in your system such that all residues are in their neutral form.
```
% ./pkaS-addH
setenv UHBD_EXEC /usr/local/uhbd
cp ./pkaS.pdb ./protein.pdb
/usr/local/uhbd/pdb2charmms
/usr/local/uhbd/charmm
/usr/local/uhbd/charmm2pdbs
```
Now look at the resulting pdb file (vi pkaSH.pdb). Look to ensure that all atoms have reasonable coords (i.e. that none have 999.00 as coord values) and that the N- and C-terminal residues have the right numbers of hydrogens. Look at the ASP and GLU residues to see that the COO groups have a hydrogen (COOH), ... Now exit out of vi. Also look near the end of the log file (pkaS-hbuild.log) to see if any errors were reported by CHARMM (i.e. pay special attention to the end of this file to see whether a NORMAL TERMINATION BY NORMAL STOP was attained). Normally you have to delete the 1st line of the pkaSH.pdb file, because it is the copy of the 2nd line.

Now, check to make sure that any non-amino acid groups that you have in your pdb file have parameters in the pkaS.dat file. If not (as is true in this case), edit that file and add some (there are no parameters for a residue named FE):

...... lines not shown ......
ACH  H93   0.092  0.000  0.000  1.200
CA   CA    2.000  0.000  0.000  1.900
SO4  S    -2.000  0.000  0.000  1.900
ZN   ZN    2.000  0.000  0.000  1.900
MN   MN    2.000  0.000  0.000  1.900
HEM  FE    2.000  0.000  0.000  0.740
...... lines not shown ......

(add the FE params as shown)

...... lines not shown ......
ACH  H93   0.092  0.000  0.000  1.200
CA   CA    2.000  0.000  0.000  1.900
SO4  S    -2.000  0.000  0.000  1.900
ZN   ZN    2.000  0.000  0.000  1.900
MN   MN    2.000  0.000  0.000  1.900
FE   FE    2.000  0.000  0.000  1.900
HEM  FE    2.000  0.000  0.000  0.740
...... lines not shown ......

One final file needs to be edited (vi pkaS-dops.inp). Enter a proper header card and set the pH to the desired value and make sure that the other numerical values are set as shown.

                                         header for the run
Porin OmpF at pH 7
                                pH and absolute temperature
7.00 293.0
          number of MC steps, number of initial steps, seed
1000000 100000 23313.0
                                               parameter s1
0.99
                      number and data on additional charges
0

Now run the pKa calculation (this will take ca. 40 minutes so take this time to work on your projects):
./pkaS-dosbs >& pka.log
You will get lots of repeating information going to the pka.log file. Once the calc. has finished, look at the file: pkaS-dops.out for the distribution of ionization states.

***Checking the result files for HEMERYTHRIN (1HMO)***
( HEMERYTHRIN (1HMO) 에 대한 계산 결과보기)
(Sorry for not giving the result files of OmpF. It will be updated soon.)

pkaS-dops.out :

HEME


INPUT DATA:

    Temperature:                        pH:
          293.00                       7.00

 No of MC steps:    No of initial MS steps:
         1000000                     100000

           seed:               parameter s:
        23313.0                       0.990



RESULTS:


Protein charge from all steps      -3.27 +/-     0.00

Protein charge from mean occupancies      -3.29

First protein energy 0.115892D+03 kcal/mol

Minimal free energy found -.124325D+03 kcal/mol

Free energy of the protein -.132142D+03 kcal/mol

Number of accepted states     641111

Number of changes done        460297

Mean final occupancies, charges and pKs of sites

Final pK << pH if occupancy of the site is 0
Final pK >> pH if occupancy of the site is 1

     titratable site    occupancy    charge    pK_app   pK_intr
   1  N    GLYN    1       0.5036     0.504      7.01      7.10
   2  CG   ASP     6       0.0005    -0.999      3.72      4.05
   3  OH   TYR     8       1.0000     0.000     >> pH     10.93
   4  SG   CYS     9       0.9865    -0.013      8.87      8.56
   5  CG   ASP    11       0.0000    -1.000     << pH      3.88
   6  CZ   ARG    15       1.0000     1.000     >> pH     11.36
   7  OH   TYR    18       1.0000     0.000     >> pH      9.52
   8  CG   ASP    22       0.0000    -1.000      2.34      2.74
   9  CG   ASP    23       0.0001    -1.000      3.09      3.55
  10  CD   GLU    24       0.0003    -1.000      3.44      3.29
  11  NE2  HISA   25       0.0000     0.000     << pH     -7.54
  12  NZ   LYS    26       1.0000     1.000     11.40      9.16
  13  CG   ASP    39       0.0011    -0.999      4.04      3.80
  14  CG   ASP    42       0.0007    -0.999      3.82      3.65
  15  NE2  HISA   43       0.1298     0.130      6.17      5.42
  16  CD   GLU    46       0.0003    -1.000      3.44      4.05
  17  CZ   ARG    48       1.0000     1.000     >> pH     10.97
  18  CZ   ARG    49       1.0000     1.000     >> pH     11.30
  19  SG   CYS    50       0.9273    -0.073      8.11      8.69
  20  NZ   LYS    53       0.9993     0.999     10.13      9.04
  21  NE2  HISA   54       0.0000     0.000     << pH     -7.13
  22  CD   GLU    58       0.0000    -1.000     << pH     -4.86
  23  OH   TYR    67       0.9956    -0.004      9.36      9.21
  24  OH   TYR    70       0.9990    -0.001     10.01     10.06
  25  CD   GLU    72       0.0010    -0.999      3.99      4.01
  26  NE2  HISA   73       0.0000     0.000     << pH     -5.03
  27  NZ   LYS    74       0.9998     1.000     10.75      9.97
  28  NZ   LYS    75       1.0000     1.000     >> pH     10.19
  29  NE2  HISA   77       0.0000     0.000     << pH     -6.48
  30  CG   ASP    78       0.0000    -1.000      1.90      3.77
  31  CG   ASP    79       0.0001    -1.000      3.16      4.15
  32  NE2  HISA   82       0.1668     0.167      6.30      6.13
  33  NZ   LYS    83       0.9999     1.000     11.05      9.72
  34  CG   ASP    85       0.0003    -1.000      3.53      3.88
  35  CG   ASP    88       0.0009    -0.999      3.97      3.67
  36  CG   ASP    90       0.0011    -0.999      4.04      3.46
  37  OH   TYR    93       0.9983    -0.002      9.78      9.83
  38  NZ   LYS    95       0.9997     1.000     10.55      9.87
  39  NE2  HISA  101       0.0000     0.000     << pH     -6.43
  40  NZ   LYS   103       0.9964     0.996      9.44      8.66
  41  CG   ASP   106       0.0000    -1.000     << pH     -5.27
  42  NZ   LYS   108       0.9985     0.999      9.83      9.97
  43  OH   TYR   109       0.9994    -0.001     10.22      4.59
  44  CZ   ARG   110       1.0000     1.000     >> pH     11.36
  45  NZ   LYS   112       0.9990     0.999      9.98     10.04
  46  C    ILEC  115       0.0000    -1.000      2.53      3.19

Charges of sites for most probable state
   1  1.  0.  1.  0.  1.  0.  1.  0.  0.  1.  1.  1.  0.  1.  0.  0.  1.  1.  0.  0.
   2 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
   3  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
   4  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  0.  0. -1.  0.  0.  0.  0.  0.
   5 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
   6  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
   7  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
   8 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
   9 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  10 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  11  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  12  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  13 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  14 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  15  0.  0.  0.  0.  1.  1.  0.  0.  1.  1.  0.  0.  0.  1.  0.  1.  0.  0.  0.  0.
  16 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  17  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  18  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  19  0.  0.  0.  0.  0.  0. -1. -1.  0.  0. -1.  0. -1. -1.  0. -1.  0.  0.  0.  0.
  20  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  21  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  22 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  23  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1. -1.  0.
  24  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  25 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  26  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  27  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  28  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  29  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  30 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  31 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  32  0.  0.  1.  1.  0.  0.  0.  0.  1.  1.  1.  0.  1.  0.  0.  0.  0.  0.  0.  0.
  33  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  34 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  35 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  36 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  37  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  38  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  39  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  40  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  0.
  41 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
  42  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  43  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  44  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  45  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  46 -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1. -1.
     -3. -4. -2. -3. -2. -3. -4. -5. -2. -1. -3. -4. -4. -3. -5. -4. -4. -4. -5. -5.

Free electrostatic energies of 20 lowest states
   1    -.124324867925D+03 kcal/mol
   2    -.124302478390D+03 kcal/mol
   3    -.123373072944D+03 kcal/mol
   4    -.123370862409D+03 kcal/mol
   5    -.123194999014D+03 kcal/mol
   6    -.123180065509D+03 kcal/mol
   7    -.122794275556D+03 kcal/mol
   8    -.122755982322D+03 kcal/mol
   9    -.122234359329D+03 kcal/mol
  10    -.122229113834D+03 kcal/mol
  11    -.121871963476D+03 kcal/mol
  12    -.121868415206D+03 kcal/mol
  13    -.121853849241D+03 kcal/mol
  14    -.121847561646D+03 kcal/mol
  15    -.121821614872D+03 kcal/mol
  16    -.121816724441D+03 kcal/mol
  17    -.121333239585D+03 kcal/mol
  18    -.121298564718D+03 kcal/mol
  19    -.121218807284D+03 kcal/mol
  20    -.121120569050D+03 kcal/mol

Note that the two lowest energy states (1 and 2) have almost the same energy. These should ocnsidered to be degenerate states and a close inspection should be made of the difference betweeen these two states.
- Upon closer inspection, we see that the N-terminus is ionized in the 1st state and neutral in the second state.
```
   1  1.  0.
```
- Further, we see that the N-terminus is predicted to be only 50% ionized (i.e. an occupancy of the proton of 0.5036).
```
     titratable site    occupancy    charge    pK_app   pK_intr
   1  N    GLYN    1       0.5036     0.504      7.01      7.10
```
- We would need to look at the enzyme to see where this residue lies to see whether it's ionization state is important for the further studies that we wish to perform. It's very likely not to be important and is probably totally solvent exposed. I would treat it as ionized unless it was close to the catalytic site.

For a better estimate of the predicted pKas, you should look at the file from the hybrid alanysis (hybrid.out), included below:

     21
   ph     mean q     Gion       dGion      CPU(s)   gcut   tcut    Gelec
   -5.00   15.46     -315.27     96.71    0.012464  0.000  0.050     16.47
   -4.00   15.07     -294.23     88.73    0.010796  0.000  0.050     10.39
   -3.00   14.91     -273.52     80.43    0.020011  0.000  0.050      8.32
   -2.00   14.48     -253.17     71.77    0.037045  0.000  0.050      3.34
   -1.00   14.08     -233.48     62.44    0.076233  0.000  0.050     -1.31
    0.00   13.95     -214.13     52.79    0.076418  0.500  0.050     -2.34
    1.00   13.48     -195.08     42.83    0.127944  0.700  0.050     -3.63
    2.00   11.44     -177.62     31.36    0.385696  0.700  0.050     -8.34
    3.00    7.19     -164.53     16.22    0.194681  0.400  0.050    -15.36
    4.00    2.15     -158.19     -0.94    0.214232  0.400  0.050    -20.60
    5.00   -0.56     -157.43    -12.55    0.118689  0.400  0.050    -22.32
    6.00   -1.78     -159.04    -19.29    0.089442  0.400  0.050    -22.69
    7.00   -3.28     -162.55    -22.41    0.103565  0.400  0.050    -23.08
    8.00   -4.51     -167.93    -23.13    0.115720  0.300  0.050    -23.40
    9.00   -6.21     -175.24    -22.55    0.131236  0.300  0.050    -23.39
   10.00   -8.94     -185.59    -19.50    0.071692  0.300  0.050    -22.69
   11.00  -12.67     -200.41    -11.77    0.071817  0.300  0.050    -19.40
   12.00  -16.88     -220.83     -2.95    0.071364  0.300  0.050    -12.81
   13.00  -20.75     -246.89      3.96    0.077185  0.200  0.050     -4.05
   14.00  -23.54     -277.66      7.52    0.013917  0.000  0.050      4.01
   15.00  -24.74     -311.18      8.51    0.008982  0.000  0.050      7.85
     46
 Group    pk(model)  pK(app)     pK(app) - pK(model)
   1      7.500        7.042          -0.458
   2      4.000        3.558          -0.442
   3      9.600       13.357           3.757
   4      8.300        8.951           0.651
   5      4.000        3.081          -0.919
   6     12.000       12.418           0.418
   7      9.600       12.420           2.820
   8      4.000        1.952          -2.048
   9      4.000        2.817          -1.183
  10      4.400        3.382          -1.018
  11      6.300      999.000         992.700
  12     10.400       11.150           0.750
  13      4.000        3.757          -0.243
  14      4.000        3.350          -0.650
  15      6.300        6.188          -0.112
  16      4.400        2.983          -1.417
  17     12.000       13.270           1.270
  18     12.000       14.223           2.223
  19      8.300        8.226          -0.074
  20     10.400       11.236           0.836
  21      6.300      999.000         992.700
  22      4.400      999.000         994.600
  23      9.600        9.575          -0.025
  24      9.600       10.368           0.768
  25      4.400        3.818          -0.582
  26      6.300       -2.102          -8.402
  27     10.400       12.249           1.849
  28     10.400       12.120           1.720
  29      6.300      999.000         992.700
  30      4.000        2.508          -1.492
  31      4.000        2.026          -1.974
  32      6.300        6.318           0.018
  33     10.400       11.978           1.578
  34      4.000        2.522          -1.478
  35      4.000        3.726          -0.274
  36      4.000        3.381          -0.619
  37      9.600       10.377           0.777
  38     10.400       11.130           0.730
  39      6.300      999.000         992.700
  40     10.400        9.566          -0.834
  41      4.000      999.000         995.000
  42     10.400       10.502           0.102
  43      9.600       11.845           2.245
  44     12.000       13.086           1.086
  45     10.400       11.197           0.797
  46      3.800        3.359          -0.441
   -5.00 (pH)
    1.000
    0.000
    0.000
    0.000
    0.000
    1.000
    0.000
    0.000
    0.000
    0.000
    0.079
    1.000
    0.000
    0.000
    1.000
    0.000
    1.000
    1.000
    0.000
    1.000
    0.227
   -0.877
    0.000
    0.000
    0.000
    0.958
    1.000
    1.000
    0.011
    0.000
    0.000
    1.000
    1.000
    0.000
    0.000
    0.000
    0.000
    1.000
    0.030
    1.000
   -0.971
    1.000
    0.000
    1.000
    1.000
    0.000
   -4.00 (pH)
    1.000
    0.000
    0.000
    0.000
    0.000
    1.000
    0.000
    0.000
    0.000
    0.000
    0.014
    1.000
    0.000
    0.000
    1.000
    0.000
    1.000
    1.000
    0.000
    1.000
    0.041
   -0.961
    0.000
    0.000
    0.000
    0.940
    1.000
    1.000
    0.007
    0.000
    0.000
    1.000
    1.000
    0.000
    0.000
    0.000
    0.000
    1.000
    0.015
    1.000
   -0.991
    1.000
    0.000
    1.000
    1.000
    0.000
...... remaining lines not shown ......

Take a look at the hybrid.chart for a more readable form of these results.
Some ionizable residues have predicted pKas that are in the range of 999. This just means that the residue titrates beyond the pH ranges considered. It usually also means that there are titratable groups nearby that are more ionizable and prevent it from ionizing.
Note that this file has ionization state information for every ionizable residue for a wide range of pH values. It is possible to plot a titration curve for whatever residue(s) that you want with a simple conversion program that I have written (/usr/local/pka/hybrid2xmgr). You will be asked for the residue number for the ionizable residue that you want data for. The plotting program: xmgr (/usr/local/pka/xmgr) can be used to make the plot(s).
If you want to continue on to do an electrostatics calculation, then run the following program: /usr/local/pka/dops2uhbd
A file will result with the name: editq.out with commands in it appropriate for using in a UHBD input file for an electrostatics calculation. You should take one like I gave in one of the earlier computational labs on electrostatics calculations
```
edit charge   0.19 atnum     1 end
edit charge  -0.30 atnum    46 end
edit charge  -0.65 atnum    69 end
edit charge  -0.30 atnum   100 end
edit charge   0.80 atnum   143 end
edit charge  -0.65 atnum   182 end
edit charge  -0.30 atnum   217 end
...... lines not shown ......
```

These commands should be used in the following way:

read mol1 file "molecule.pdb" pdb end
set charge radii file "qcharmm.dat" para mine end

edit charge   0.19 atnum 1 end
edit charge  -0.30 atnum    46 end
edit charge  -0.65 atnum    69 end
edit charge  -0.30 atnum   100 end
edit charge   0.80 atnum   143 end
edit charge  -0.65 atnum   182 end
edit charge  -0.30 atnum   217 end
...... lines not shown ......

print apar mol1 end

elec calc mol1
...... lines not shown ......

***Checking the result files for HIV Integrase***
(HIV Integrase에 대한 계산 결과보기)

Homework with HslV (1e94) structure [download ]