# Calculation of pH in Molecular Dynamics simulation?

I have a system that has 5M aluminum hydroxide (10 molecules) ion [Al(OH)4]- and 5M NaOH (10 molecules) in a box of side length 15 angstroms each. I have 113 water molecules in the system, so that the density of water is consistent with 1g/cc. I need a way to find the pH of the simulation. I tried calculating it by counting the H+ ions, but that gave me a value of ~1.2, which did not really make sense.

• +1 but have you tried using "activity" instead of "concentration"? Commented Mar 8 at 5:03
• I am trying to match the simulations with experiments. So I need to know the pH, because the experiment was conducted at a certain pH. Commented Mar 8 at 12:02
• How would you count H+ ions if there isn't any? pH can be inferred from how many hydroxide ions are there. Commented Mar 8 at 12:56
• What kind of force field are you using? Or is this ab initio? Typical molecular dynamics simulations have fixed protonation states that you set initially. Another option is constant-pH MD, where you set the pH and the protonation states fluctuate. If you are using a reactive force field (ReaxFF?), it might be possible to observe spontaneous deprotonation/protonation, but the number of hydronium ions may be too small to observe. Commented Mar 8 at 14:26
• I'm using a modified implementation of ReaxFF. Commented Mar 8 at 14:40

What matters is the number of "free" H+ ions. A quick estimation using the autoionization constant of water shows that the H+ concentration is on the order of $$10^{-15}$$ mol/L. Obviously, it is essentially impossible to see even one H+ ion in a simulation of merely 100 water molecules. If you do see at least one H+, please check if you have used a wrong criterion for determining whether a hydrogen atom is H+, or if your MD simulation setup has a severe problem (for example you have added the wrong number of Na+ ions).

In any case, if you do not need very high accuracy, say if your result only needs to be accurate to about $$\pm$$1, then this is an elementary analytical chemistry exercise. You can simply look up the experimental pKb of [Al(OH)4]- and solve for the equilibrium concentration of OH-; then, using the autoionization constant of water, you can easily find the concentration of H+ and therefore the pH. For higher accuracy, you may need to consider the activity coefficients (https://en.wikipedia.org/wiki/Activity_coefficient), which can be rather accurately calculated by the ionic strength, the latter of which can be calculated easily from the components of the solution. Still, this is something you can calculate by hand within one hour. By contrast, MD simulations probably have larger error bars than what can be obtained by the above procedure (due to method and statistical errors), and is much more involved, so I don't recommend using MD to solve this problem.

• As a quick estimate -- 1L = 1e24 nm^3, so 1M is about 0.6 particles (or ions or molecules) per cubic nanometer. Commented Mar 8 at 9:55
• I know the OP started the phrasing of "free H+ ions". But typically we are taught that free hydrogen ions dont exist in water -- cuz if it did, it'd bind to a water molecule to form hydronium. So it's all hydronium ions. Is hydronium ion what you mean when you wrote free hydrogen ions? Commented Mar 8 at 15:25
• I do not have free H+ ions in the box. What I have are 10 OH- ions from NaOH and 10 [Al(OH)4]- ions in the box. Commented Mar 8 at 17:30
• @Argyll the "hydronium ion" is useful for chemical reasoning (HCl doesn't spontaneously ionise -- it transfers a proton to water). But at molecular level H+ is solvated by anywhere from 2 to 20 water molecules at once -- so "H3O+" isn't fully accurate either, even though it's a good model for chem ed. See J. Chem. Educ. 2011, 88, 7, 875 and related papers for more info. Commented Mar 9 at 1:15
• @ShernRenTee: I see. Thank you for the reference. Commented Mar 9 at 1:57

Here's some useful intuition for solution molecular dynamics:

One litre contains $$10^{24}$$ cubic nanometers. One mole contains $$6 \times 10^{23}$$ particles. So 1M of volumetric density corresponds to a number density of 0.6 per cubic nanometers.

The self-dissociation constant of water is $$K_w = 10^{-14} [\mathrm{M}^2]$$. So at equilibrium and pH 7, both $$\mathrm{H}^+$$ and $$\mathrm{OH}^-$$ will have a concentration of $${10^{-7}}$$ M.

By the above calculation, this is 60 ions of each type per cubic micron.

As such, it is not feasible to directly simulate water dissociation to probe pH. Rather, the process is modelled using grand-canonical Monte Carlo exchange with a reservoir (essentially "chemostatting", in analogy to barostatting and thermostatting). For example this recent paper chemostats on the combination of acid/base dissociation and electrolyte ion insertion/removal to reach feasible reaction probabilities.

ReaxFF results would certainly not be trustworthy until proven otherwise, since how could they possibly calibrate parameters to exactly achieve sub-micromolar concentrations of hydronium and hydroxide in pure water?