In reference 1, they calculated pKa of blocked Aspartic acid and Histidine, here I am doing the same thing, the purpose is to find out if the pKa is different from residue in protein. Follwing is an example of aspartic acid(positive amino acids:ARG HIS LYS does not need add_h.inp, negative ones ASP GLU need, TYR does not need to add H).

To prepare mono amino acid,
charmm name=asp < build_amino.inp(outputs:data/asp.crd|pdb|psf )
add H:
charmm in=asp <add_h.inp(outputs:data/asp_h.crd|pdb|psf)
Run a constant pH simulation with implicit solvent:
mpirun -np 4 charmm in=asp <hphmd.inp(1ns, around 7 minutes. outputs:data/aspartic.ph3.7.dcd|ene|rst|lamb)(set struct =data/@in if ARG HIS LYS TYR)
Calculate unprotonated fraction (S) by CptSX.pl:
CptSX.pl asp.ph3.7.lamb 10000 50000 3.7(outputs:*.sx)

CptSX.pl asp.ph3.7.lamb 1 1000 3.7(outputs:*.sx)
Calculate pKa values by fitting to Henderson-Hasselbalch (HH) equation by cptpKa.pl:
cptpKa.pl asp asp.ph3.7.lamb-10000-50000.sx(outputs:*.pka)

cptpKa.pl asp asp.ph3.7.lamb-1-1000.sx(outputs:*.pka)
cat asp.ph3.7.lamb-10000-50000.pka
1 3.56 1.0

Following are pKa's for other residues:

Let's run replica exchange pH simulation:
mpirun -np 14 charmm  <hphmd_phrex_gb.inp
outputs:data/asp.ph*.rst/dcd/ene/lamb_IREP, asp.phrex.log/out_IREP.
The *.out_IREP files are output from each replica, the *.log_IREP files contain the information about exchange attempts (ACCEPT or REJECT), and the energy of each replica.
Now we calculate S for each replica:
CptSX.pl asp.ph1(2,3,4,5).lamb_0(1,2,3,4) 5 1000 1(2,3,4,5)

or

for i in $(seq 1 14); do j=`expr $i-1|bc`; CptSX.pl asp.ph$i.lamb_$j 5 1000 $j; done
Then calculate pKa:
cptpKa.pl asp asp.ph1.lamb_0-10000-50000.sx asp.ph2.lamb_1-10000-50000.sx asp.ph3.lamb_2-10000-50000.sx asp.ph4.lamb_3-10000-50000.sx asp.ph5.lamb_4-10000-50000.sx

or

string=

for i in $(ls *.sx); do string+=" $i"; done

cptpKa.pl asp $string

cat asp.pka
1 4.19 0.89 1.00

Following are pKa's for other residues:

Now let's try explicit solvent model,

charmm in=asp <solvate.inp(outputs:aspartic_h_solvated.pdb/psf/crd/str)
Add ions to the system,
charmm in=asp conc=0.100 <add_ions_hphmd.inp(outputs:asp_h_ions.pdb/crd/str/prm)
Now we combine the structures,edit minimize.inp, and mark off "join SOLV CLA renumber", since the system only has SOD.
charmm in=asp < minimize.inp (outputs:asp_h.solv.ion.min.psf/pdb/crd/pbc.)
The system should looks like this:

Let's run replica exchange pH simulation on this system:
We still use the same hphmd_phrex.inp, but we need to use periodic boundary condition and turn on Image, add hybrid to gbsw. Don't forget to change struct=data/@in_h.solv.ion.min
mpirun -np 14 charmm < hphmd_phrex_ex.inp (pH from 1 to 14, 1.0ns)
pKa=3.61,the titration curve is:

 

Following are the pKa's of 6 titratable residues from 1ns 14 replica exchange simulation

Dimer Aminoacids

Following are the results from dimers of amino acids.

edit built_amino.inp, "read sequ @resname 2"(change from 1 to 2)

Following are the results from dimers of 14 replica GB

Dimer explicit:

For Tyrosine, add following to phmd.in:

TYR 10.1 A B 2.0 (A and B need to be trained)
CZ 0.11 0.40
OH -0.54 -0.40
HH 0.43 0.00 1.0 0.0

For ARG, add following to phmd.in
ARG 12.1 0.0 0.0 0.0
CD 0.20 0.00
HD1 0.09 0.09
HD2 0.09 0.09
NE -0.70 -0.70
HE 0.44 0.43
CZ 0.64 0.40
NH1 -0.80 -0.62
HH11 0.46 0.31
HH12 0.46 0.00 1.0 0.0
NH2 -0.80 -0.62
HH21 0.46 0.31
HH22 0.46 0.31

To find parameter A and B, set A=0, B=0,2.0->0.0. run:

charmm < hphmd_derive.inp

outputs: *.lamb under data/

gnuplot--> stats "*.lamb" using 3

TYR (second column is A=0,B=0,2.0.not too much different)

0.2 7.7521 7.74
0.3 10.8119 10.8124
0.4 13.0275 12.9593
0.5 14.2578 14.2270
0.6 14.5236 14.5245
0.7854 12.7851 12.7876
0.9 10.7390 10.6879
1.0 8.7023 8.6971
1.1 6.5819 6.5819
1.2 4.7294 4.7294
1.3 3.1584 3.1584
1.4 1.8290 1.8290
A= -7.76075 B=1.32243

ARG
0.2 76.3238
0.3 107.2953
0.4 130.2527
0.5 144.2557
0.6 148.9769
0.7854 136.3344
0.9 118.4906
1.0 99.6261
1.1 74.9732
1.2 59.9357
1.3 41.7497
1.4 25.2668
A=-65.8466 B=1.53207

LYS

0.2 43.1129
0.3 58.6399
0.4 67.2173
0.5 68.9941
0.6 63.5223
0.7854 40.6659
0.9 22.8726
1.0 8.8289
1.1 -3.3910
1.2 -10.4622
1.3 -12.7697
1.4 -10.5828
A=-76.1272 B=0.767769

References:
1)Toward the accurate first-principles prediction of ionization equilibria in proteins. Jana Khandogin,Brooks III 2006.