Calculating pKa's using Gaussian

I've been learning to use Gaussian recently and I found a question here that seemed like a good opportunity for me to practice. I know that I can obtain $K_a$ from $$\Delta G=-RT\ln(K_a)$$ using the general reaction $$\ce{HA(aq) + H_2O(aq)->A-(aq) +H_3O+(aq)}$$ and taking the sum of thermal and electronic free energy for each compound as their value of $G$ However, I'm not confident that the method I'm using in Gaussian is correct.

I performed an optimization/frequency calculations on all these molecules using the APFD functional, with basis set 6-311G+(2d,p) and SMD solvation using water.

After working out my calculation, I obtained $pK_a=26.10$. My main question comes down to: did I use a sufficient level of theory, at least approximately, to capture the effect of solvation? I've looked at the Gaussian white book for thermochemistry as well as a few journal articles (1 2) and couldn't get a clear sense of if I had the right level of theory or if my calculations were completely off-base.

Another finer point might be: are there any fairly similar experimental values to which I could compare my calculation to see if I'm in the right ballpark?

• Such a small model for water protonation is generally way off, and you will not improve it by eg better basis set. – Greg Mar 13 '17 at 23:08
• @Greg So would it require better theory or a better solvation model(like explicit solvation) to make this work? If it winds up being too burdensome a calculation I might just leave it, but I would be interested to see any resources that point to a sufficient model size for this type of problem. – Tyberius Mar 13 '17 at 23:52
• In my experience, the solvation is the critical part that often fails when heterolytic dissociation is studied. – Greg Mar 14 '17 at 0:14
• Are you just calculating the free energy change of the dissociation using implicit solvation? One thing that you could do is to use a more suitable thermodynamic cycle. You can optimise the structures at a lower level of theory and then calculate electronic energies in the gas phase and in implicit solvent at a higher level to get better solvation free energies while using the experimental value of the solvation energy of the proton (-265.9 kcal/mol). – Snakeysnakes Mar 14 '17 at 10:20
• The approach itself is flawed. There is no H3O+ in water. – AMT Mar 14 '17 at 13:13

This is unfortunately not a real answer to your question as such, but I'd like to provide some pointers why this approach is a lot more complicated than it seems at first. It might very well be (currently) impossible.
That is probably also a reason why you don't find a lot of publications on the theory of it and so far no black-box method has been developed. Solvation is still one of the major problems in computational chemistry, hence the use of implicit solvents and counter-charges. Such methods naturally come with the implicit disclaimer that these are approximate corrections only.
The most elaborate method I know so far is COSMO-RS, which still has large error margins, at least that is the my current state of knowledge.

I think for calculations you would at least need to run molecular dynamics simulations with many explicit solvent molecules to describe the equilibrium. This is important, because the proton which dissociates will be caught in the bulk of water molecules and a highly complex network of hydrogen bonds. This is very well demonstrated by the Grotthuss mechanism, which shows how easily bond lengths can change. I would say it is impossible to describe this complex system with a single quantum chemical calculation.

Furthermore, DFT is notoriously bad at describing non-covalent interactions, and including dispersion corrections will only get you so far. I have not checked recently, but I believe there are quite a few publications discussing various sizes of water clusters on this level. I have neither used nor heard of AFPD before, so it would be wise to thoroughly search for literature (benchmarks) about its performance on these clusters. Suffice to say, that one particular functional will not give you a reliable result as you have nothing to compare it to. Especially in such complicated cases you have to calibrate your results.

Another problem here is that you are the describing heterolytic bond dissociation. More often than not this may lead to additional problems. I don't always have confidence in my results when I treat simple donor-acceptor complexes where there is no charge separation and all singlet species. A supermolecular approach would probably a better choice, but could heavily suffer from a BSSE. However, it might be worth considering $$\ce{H2O\bond{...}HA <=> H2OH\bond{...}A}.$$

I think the key for calculating acid dissociation constants is treating the whole system as one ensemble.

Recall the common definition of the equilibrium constant for the reaction \begin{align}\ce{AH (aq) &<=> A- (aq) + H+ (aq)},& K &= \frac{a(\ce{A-})\,a(\ce{H+})}{a(\ce{HA})}.\end{align} We conveniently ignore water here, because it is abundant and we assume ideal solution. However, when your system is too small, you will get a large error for $K$. This is most likely due to your solvated proton, which will be described worst with solvation models. Explicitly you are writing $$\ce{HA (aq) + H2O (aq) -> A− (aq) + H3O+ (aq)},$$ which is not really what you are calculating. Instead you are (more ore less) using already approximations to treat thermal and solvent corrections $$\ce{HA (g, th, solv) + H2O (g, th, solv) -> A− (g, th, solv) + H3O+ (g, th, solv)}.$$ And this doesn't even consider that you are calculating two states in an equilibrium and you will get the equilibrium constant (using activities, not concentrations) for only these two states (which might not even be the ground sates), $$K = \exp\left\{-\frac{\Delta G}{RT}\right\}.$$ You have also not accounted for the possible conformational space of $\ce{HA}$/$\ce{A-}$, which might become very important as soon as you have ethyl moieties, or even cyclohexyl moieties.

These are just some issues that immediately come to mind, there might very well be more underlying insufficiencies. I have looked into this a couple of years ago and came to the conclusion not to bother with it any more.

Assessing acid strengths does not need to be at the level of what experimental chemists understand. If they really need the $\ce{pKa}$ they should run the appropriate measurements and obtain reliable results instead of educated guessing. If they want you to classify whether an acid is stronger or weaker, use a rigorous approach to obtain the conformational space for all species and compare them with isodesmic equations. This will have the advantage that insufficiencies of the methodology will likely cancel each other.

Another purely theoretical approach would be calculating gas-phase proton affinities. From this you can judge if a compound behaves stronger or weaker. Keep in mind that this approach is better suited for bases or nucleophiles. If you have the resources, you can also evaluate Lewis acidities, $$\ce{Base + BX3 <=> Base-BX3}.$$

This should give you at least a good idea how on a relative level your compounds behave.

You can also use a first-principle derived value for your proton, which will likely cause the largest error. A quick check turned up a few promising results. (For example: Chang-Guo Zhan and David A. Dixon. J. Phys. Chem. A 2001, 105 (51), 11534–11540.) You can then go ahead and calculate know compounds and try to find a correlation between $\ce{pKa}$ and $\Delta G$. This way of "calibration" might lead to some significant error compensation.

I would also recommend to look at COSMO-RS in more detail. There might be a lot of development that I have missed and they might have found a way to obtain more reliable results. As far as I remember, the foundation of this was solid and a good step into easier capturing solvent effects without the need to resort to full blown MD simulations.

• Thank you! This is definitely helpful and gives me a nice range of ways that I could go forward with the calculation or look at the problem from a different angle. I'm might try using something other than DFT (I defaulted to using it since it was what is used in general in book I'm using to learn Gaussview). I also think looking at the supermolecular equilibrium could be promising. – Tyberius Mar 15 '17 at 20:18
• @Tyberius In my opinion GaussView is not a programme worth learning as it is only a visualisation tool of the results, not even a good one. Learn how Gaussian operates, how an input needs to be structured, what different things get reported in the output, how you can analyse those results further. When you do this, it's easier to switch to other calculation packages and analysis tools. I have to reiterate, without a proper size assessment and MD treatment, you won't get reasonable results. – Martin - マーチン Mar 16 '17 at 0:01
• I misphrased that, it is mainly for Gaussian and includes the specific structure of the input/output files, though it gives instructions for how to view results or set up jobs within Gaussview. But yeah I can see that finding a pKa for an individual compound might be outside of the scope of strictly Gaussian calculations without supplementing with MD, experiment, or some calibration of with a series of similar compounds. – Tyberius Mar 16 '17 at 15:59
• @Tyberius Recently the topic came up on CCL (ccl.net/cgi-bin/ccl/message-new?2020+04+22+001), and Marcel Swart shared this link: comp.chem.umn.edu/solvation, which is a good resource I thought I share with you. – Martin - マーチン Apr 30 at 22:01