We have
submitted a manuscript at Molecular Physics (DOI 10.1080/00268976.2024.2396535).
Abstract:
For 21 carboxylic acids and 14 phenols, the energies of deprotonation were computed at DFT level (B3LYP, B3LYP-D4 and wB97X-D4 functionals) for building linear fits with their experimental pKa values. In general, the agreement with the experiment improves when two to eight explicit water molecules are added to the model, but the set of phenols contains three outliers which mask this trend because these outliers contain polar groups in adjacent positions, and the explicit water molecules, which are added to them, form an arbitrary number of H-bonds between each other. The mean errors of the linear fits are approximately 0.12 pKa. Extension of the basis set from 6-31++G(D,P) to aug-cc-pVTZ does not improve the linear fit accuracy. A similar conclusion can be made about the role of entropy factors: when they are included in the computation (the Gibbs energies are computed), the agreement with the experiment is not improved. Besides that, if the Gibbs energies are computed in standard ‘rigid rotor – harmonic oscillator’ approximation without using a low frequencies threshold, the agreement with the experiment becomes even worse than when the entropies are not used at all, due to the anharmonicity of small vibrational modes.
One more note: one of the conclusions of our work is that there is almost no difference between the results computed by B3LYP-D4/6-31++G(D,P) and B3LYP-D4/aug-cc-pVTZ level; in other words, the error caused by the incompleteness of the 6-31++G(D,P) basis set is nearly zero in comparison to other sources of errors. Here is the correlation between the energies computed by these two methonds:
