Is it possible to calculate the pka value of different H atoms on one molecule?

Question

I like to calculate pka values for two different H atoms in one molecule. Is it possible to do so with computational methods? I have gaussian and orca as software to use. I searched in the internet for articles about it. I often found QSAR. But I don't want to do QSAR calculation, because the can only perform pka values for the hole molecule in one. So I can get not different steps. And may do additional calculations with the geometry. So has anyone an idea if it is possible with DFT calculations or something similar?

Answer 1

Yes. Simply calculate the pKa from the definition:

$$ \rm{pK_a} = \lg \frac{[AH]}{[A^-][H^+]} \tag{1} $$

The equilibrium constant in the logarithm is computed from the Gibbs free energy of the protonation reaction $A^- + H^+ \to AH$:

$$ \frac{[AH]}{[A^-][H^+]} = e^{-\Delta G/RT} \tag{2} $$ $$ \Delta G = G(AH_{aq}) - G(A^-_{aq}) - G(H^+_{aq}) \tag{3} $$

Since in aqueous solution, the three species (especially $A^-$ and $H^+$) interact strongly with water molecules, it may be necessary to include some explicit water molecules, i.e. instead of Eq. (3), we have (note the subscript $l$; we will return to this issue later. Also note that the formula does not require that $m-n-k \ge 0$)

$$ \Delta G = G((AH(H_2O)_m)_{aq}) - G((A^-(H_2O)_n)_{aq}) - G((H^+(H_2O)_k)_{aq}) \\ - (m-n-k)G((H_2 O)_{l}) \tag{4} $$

Some researchers use the following, slightly different form, where instead of calculating the Gibbs free energy of one water molecule and multiply by $(m-n-k)$, they calculate a single water cluster with $(m-n-k)$ molecules. This has advantages and disadvantages, and I'm not sure if this is a more common (or more justified) approach:

$$ \Delta G = G((AH(H_2O)_m)_{aq}) - G((A^-(H_2O)_n)_{aq}) - G((H^+(H_2O)_k)_{aq}) \\ - G(((H_2 O)_{(m-n-k)})_{l}) \tag{5} $$

(this obviously assumes that $m-n-k \ge 0$. When $m-n-k < 0$, one should use $$ \Delta G = G((AH(H_2O)_m)_{aq}) - G((A^-(H_2O)_n)_{aq}) - G((H^+(H_2O)_k)_{aq}) + G(((H_2 O)_{(n+k-m)})_{l}) $$ instead.)

Frequently (but not always), $k=4$ is used. The values of $m$ and $n$ depend on the structure of the acid $AH$; normally, they are chosen as the minimum numbers necessary to describe all strong hydrogen bonds between the solute and the water. Here "strong hydrogen bonds" can be loosely understood as hydrogen bonds that are significantly stronger than the hydrogen bonds between water molecules.

Now the only remaining problem is the calculation of Gibbs free energies involved in Eq. (4) or Eq. (5). There are lots of tutorials for this, like http://gaussian.com/videos/ (for Gaussian) and https://www.orcasoftware.de/tutorials_orca/prop/thermo.html (for ORCA). I would just stress some points that are missed or not sufficiently highlighted by these tutorials:

1. A configurational search is recommended before the geometry optimization, in order to guarantee that you are working on representative, low-energy molecular configurations, and is absolutely necessary if you have explicit water molecules (otherwise you cannot determine the correct orientations or even locations of the water molecules). Instead of just picking the lowest energy (or lowest Gibbs free energy) conformer, an even better alternative is to pick many low energy conformers, do Gibbs free energy calculations on all of them and take the Boltzmann average of the results. Configurational searches are best done by external programs like CREST, Molclus or ABCluster, at semiempirical (preferentially GFN2-xTB) levels of theory.
2. It is recommended to use a solvation model during the geometry optimization. The SMD model is the most accurate among common implicit solvation models, but some people say that it is less numerically stable than PCM and CPCM, which can prevent geometry convergence. Thus, the PCM (for Gaussian) and CPCM (for ORCA) solvation models can be recommended for the geometry optimization.
3. The dispersion correction must be added during the geometry optimization, frequency calculation and the single point energy calculation, unless the method already fully takes dispersion into account (e.g. M06-2X, wB97M-V, DLPNO-CCSD(T)).
4. After the frequency calculation, the user must check for imaginary frequencies, and if there are any, eliminate them. Otherwise the Gibbs free energies can have errors up to 2~3 kcal/mol (see my answer here Is it right to neglect very small imaginary frequencies?), corresponding to a pKa error of 1~2 units, for each imaginary frequency.
5. For the single point energy calculation, it is recommended to go to at least the DLPNO-CCSD(T)/cc-pVTZ level of theory. If possible, the TightPNO keyword and/or the F12 correction may be added to recover more correlation. If DLPNO-CCSD(T)/cc-pVTZ is not affordable, a double hybrid functional (e.g. PWPB95-D3) plus a triple-zeta basis set (e.g. def2-TZVPP) may be recommended. These may sound like an overkill for the average Gaussian user who is accustomed to a B3LYP/6-311+G(d,p) level single point calculation, but since even these high-level single point energy calculations are not much more expensive than the opt+freq calculations (at least when done by ORCA), there is little reason to not use these high-level methods.
6. Solvation effects are absolutely essential at the single point energy level. The accurate SMD method is in general recommended. It can be directly combined with DLPNO-CCSD(T) and double hybrids, although in the former case some people may prefer to do the DLPNO-CCSD(T) calculation in vacuum and add the solvation free energy calculated at e.g. the B3LYP/def2-TZVP level, citing that the SMD method is not thoroughly tested in conjunction with coupled cluster methods. Still another approach is to calculate the solvation free energy at the M05-2X/6-31G* (or the very similar M06-2X/6-31G*) level, which is justified by the fact that the SMD method was parameterized at the M05-2X/6-31G* level. There seems to be no clear consensus as to which of the three methods is the best.
7. The Gibbs free energies, as obtained by adding the Gibbs free energy corrections to the single point energies, must be further manually corrected for concentration effects. This is because the Gibbs free energies generated by most quantum chemistry programs are calculated at 1 atm (i.e. the standard state of ideal gases), which is a different state than the standard state of aqueous solutions (1 mol/L) and liquid water (55.6 mol/L). This point is frequently overlooked, especially the fact that the standard state of water in an aqueous solution is liquid water, not a "1 mol/L aqueous solution of water". To correct a Gibbs free energy from concentration $c_1$ to $c_2$, the following equation is used:

$$ G_2 - G_1 = RT\ln \frac{c_2}{c_1} \tag{6} $$

Finally, note that even when one is only interested in the first deprotonation, the above procedure is still much more reliable than QSAR, and is recommended whenever the computational resources permit. The only advantage of QSAR is that it is extremely fast to calculate, and is suitable for real-time estimation of pKa or the large-scale screening of the pKas of thousands or millions of molecules.