# pka from molecular dynamics

How can one calculate pKa values from first principle calculations, static or dynamics?

I have encountered some works using CPMD where people try to find pKa values from -(i) constrained Molecular Dynamics (MD) and (ii) metadynamics. These are somewhat more intricate ways, I believe. Is there any other simpler ways (theoretical/computational)? Can not we do it from simple MD? Please suggest some references as well, if possible.

• Initially I tried to use the basic formula for pKa, given in wikipedia.. $$pK_a=-\ln \left[ \frac{[A^-][H_3O^+]}{[HA][H_2O]}\right].$$For my case, I used this formula as follows, $$pK_a=-\ln \left[ \frac{69\times 2}{31\times 58}\right]=1.12$$ against experimental values of about 2.0. Here the numbers that I used represent not the molar concentration, but rather the proper number of molecule/ion. How well do you think the result is reliable? Commented Mar 4, 2016 at 18:40
• The calculation you suggest would be valid if the given concentrations were at the system's equilibrium. If you don't intend on directly simulating a grand-canonical ensemble, feel free to accept the (only) answer below.
– Inon
Commented Mar 8, 2016 at 17:58

• That's what I'm trying to do now, as you have suggested! I would not be able to use gromacs as my calculations are based on AIMD simulations. For free energy calculations, as statistical physics says, the formula is $F=-kT \ln Z$, where Z be the partition function. How can I calculate Z from MD trajectory? Commented Mar 4, 2016 at 18:29