# 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.

Thanks in advance.

• 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? Mar 4 '16 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
Mar 8 '16 at 17:58

## 1 Answer

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.

• 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? Mar 4 '16 at 18:29
• If you can't be bothered to look up a 'smart' way, you can do this directly by binning position and momentum space and build a histogram from real-space trajectories (gotten from a MD simulation, of course).
– Inon
Mar 4 '16 at 21:21