Modeling Flexible Molecules in Solution: A pK_a Case Study

Continuum solvation models have been incredibly successful for the computationally efficient study of chemical reactions in solution. However, their development and application has generally been on focused on investigations of small, rigid molecules. Additional factors must be considered when studying large, flexible and multiply ionizable species. These include whether the use of thermocycle or entirely solution-phase approaches are more appropriate for the calculation of solution-phase free energies, which metrics can be used to reliably identify the conformation(s) adopted by flexible molecules in solution, and how errors due to inaccuracies in the prediction of low energy vibrational frequencies can be avoided. Here we explore these issues using the calculation of pK_as for a diverse set of amine-containing species as a case study. We show that thermocycle-based approaches should only be applied where there are relatively small structural changes between the gas- and solution-phase molecular geometries, and that these methods are generally not appropriate for conformational searching. Using gas- or solution-phase energies or gas-phase free energies can also lead to errors in the identification of the most stable molecular conformation(s). Scaling of low energy vibrational modes (i.e., use of the quasi-harmonic oscillator approximation) is helpful, however care must be taken to ensure modes that change as part of the reaction are not disregarded. Entirely solution-phase approaches to the Gibbs free energy and hence pK_a calculations were found to yield accurate pK_a values for the amine test set studied when each charged site is complexed with an explicit water molecule and a proton exchange scheme is applied with an appropriately chosen reference acid.