* File name : pka_pb.inp
* Calculation of pKa 
*

set charmmdir = $CHARMMDATA

!Read topologie
open read card unit 10 name @charmmdir/top_all27_prot_lipid.rtf
read  rtf card unit 10

!Read parameters
open read card unit 20 name @charmmdir/par_all27_prot_lipid.prm
read para card unit 20

!Read Sequence
open unit 10 read form name @pdb.pdb
read sequence pdb unit 10
generate PRO0 first NTER last CTER setup
close unit 10

open read card unit 10 name @pdb.pdb
read coor pdb  unit 10 resid
close unit 10

ic param
ic build
hbuild

set  Singleresidue   =  false  ! Only the selected residue is used as a model compound.
                               ! Otherwise, C and O of previous residue, and 
                               ! N, HN, CA, and HA of next resiude are also included in 
                               ! a model compound (see, pka_charge.str)
                               ! if the next residue is PROLine, we have more (see top file)

set  npKsite   =   5           ! Number of sites for pKa calculations (titratable sites)
define pKaSite1 select resid  7 end   
define pKaSite2 select resid 10 end   
define pKaSite3 select resid 19 end   
define pKaSite4 select resid 27 end   
define pKaSite5 select resid 43 end   

open write card unit 10 name @pdb_pka.dat,epsp@epsp
write title unit 10
* site pKa_intr pKa_mod pKa_shift pKa_born pKa_back dG_born dG_back
*

stream pka_smeared.str

stop