* 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   =   1           ! Number of sites for pKa calculations (titratable sites)
define pKaSite1 select resid @resid end   

open write card unit 10 name @pdb_pka.dat,resid@resid_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