pka from molecular dynamics
Have you considered a free-energy calculation?

That is, calculate the free-energy of the dissociated species and of the neutral one, and calculate the pKa via the formula for activity.

There was this GROMACS tutorial, but it seems the link is broken.