Calculation of pH dependence for small molecules

Question

Is there any possibility to calculate with computational methods the pH dependency of protonation states for small molecules? For example I want to calculate which chemical structure formic acid will have at pH 2 and pH 15. This is a super easy example. I want to calculate it for intermediate products to see if they are stable and part of a possible pathway for some reaction. Does anyone have an idea? I have gaussian software package. Is there maybe a way to have a workaround with explicit water atoms or H+ atoms? We tried it with the acidity for the molecules (MM article), but my supervisor always asks me to calculate the state of the molecule at specific pH values. She wants me to give a few different pHs to the calculation and have a look what structure will be calculated. I have no clue how to do so, even if it makes sense if I am honest.

Answer 1

What quantum chemistry programs can give you is the pKa values of each functional group in your molecule. Using the pKa values, you can calculate the dominant form (dominant protonation state) of your molecule under different pHs from basic analytical chemistry formulas, without the further aid of a quantum chemistry program. As far as I know, no quantum chemistry program can directly give the dominant forms of a molecule under different pHs, but only Gibbs free energies of the relevant species, because it is trivial to calculate the pKa and then the pH dependence of the dominant form of your molecule from the Gibbs free energy data by hand (or, particularly if you have a lot of molecules and/or a lot of pH values to calculate, by the mathematical formula functionality of Excel, which saves you from typing the same formula over and over again on your calculator).

A simple example: suppose that your molecule is a diprotic acid ($H_2 A$), and you have calculated the two pKas of the acid:

$$ \begin{align} H_2 A &= HA^- + H^+ \quad &\rm{pKa_1} \\ HA^- &= A^{2-} + H^+ \quad &\rm{pKa_2} \end{align} $$

First of all, convert the pKa values to Ka values (equilibrium constants), and the pH value to $[H^+]$:

$$ \rm{K_a} = 10^{-\rm{pKa}} \\ [H^+] = 10^{-\rm{pH}} $$

Then, use the definition of the equilibrium constants to calculate the ratios of the different species:

$$ \rm{K_{a1}} = \frac{[H^+][HA^-]}{[H_2 A]} \\ \rm{K_{a2}} = \frac{[H^+][A^{2-}]}{[HA^-]} $$

From these one obtains $[HA^-]/[H_2 A]$ and $[A^{2-}]/[HA^-]$, which can be trivially converted to the relative proportions of each species with respect to all species, $[H_2 A]/([H_2 A]+[HA^-]+[A^{2-}])$, $[HA^-]/([H_2 A]+[HA^-]+[A^{2-}])$ and $[A^{2-}]/([H_2 A]+[HA^-]+[A^{2-}])$. If the first one is the largest, then $H_2 A$ is the dominant form at the given pH. Same for $HA^-$ and $A^{2-}$.

If your supervisor insists on seeing the molecular structure at a given pH, then it is done as follows: now that you have calculated $\rm{pKa_1}$ and $\rm{pKa_2}$, you surely have already optimized the structures of $H_2 A$, $HA^-$ and $A^{2-}$. Suppose also that by the above manual calculations, you found that $HA^-$ is the dominant form of your molecule at the given pH. Then you simply report the optimized structure of $HA^-$. Now you see that, the process is more like to compute all possible structures and then using some calculations to pick the correct one and to discard the wrong ones; the wrong ones have to be calculated anyway because they are involved in the calculations that determine which structure is the correct one. What your supervisor had in mind might be that you tell the program the pH and the program spits out the correct structure without also spitting out the wrong structures, but this is not how things actually work in ab initio calculations.

To stress again: no quantum computation needs to be redone if you want to do the calculation at a second pH. You simply calculate all the pKa values of your molecule, then use the pKa values to manually calculate the dominant form of your molecule under different pHs, and finally pick the structure of that particular protonation state from your calculation results.