Density Functional Theory Calculation of pK a ’s of Thiols in Aqueous Solution Using Explicit Water Molecules and the Polarizable Continuum Model Bishnu Thapa and H. Bernhard Schlegel* Chemistry Department, Wayne State University, Detroit, Michigan 48202, United States * S Supporting Information ABSTRACT: The pK a ’s of substituted thiols are important for understanding their properties and reactivities in applications in chemistry, biochemistry, and material chemistry. For a collection of 175 different density functionals and the SMD implicit solvation model, the average errors in the calculated pK a ’s of methanethiol and ethanethiol are almost 10 pK a units higher than for imidazole. A test set of 45 substituted thiols with pK a ’s ranging from 4 to 12 has been used to assess the performance of 8 functionals with 3 different basis sets. As expected, the basis set needs to include polarization functions on the hydrogens and diffuse functions on the heavy atoms. Solvent cavity scaling was ineffective in correcting the errors in the calculated pK a ’s. Inclusion of an explicit water molecule that is hydrogen bonded with the H of the thiol group (in neutral) or S - (in thiolates) lowers error by an average of 3.5 pK a units. With one explicit water and the SMD solvation model, pK a ’s calculated with the M06-2X, PBEPBE, BP86, and LC-BLYP functionals are found to deviate from the experimental values by about 1.5-2.0 pK a units whereas pK a ’s with the B3LYP, ωB97XD and PBEVWN5 functionals are still in error by more than 3 pK a units. The inclusion of three explicit water molecules lowers the calculated pK a further by about 4.5 pK a units. With the B3LYP and ωB97XD functionals, the calculated pK a ’s are within one unit of the experimental values whereas most other functionals used in this study underestimate the pK a ’s. This study shows that the ωB97XD functional with the 6-31+G(d,p) and 6-311++G(d,p) basis sets, and the SMD solvation model with three explicit water molecules hydrogen bonded to the sulfur produces the best result for the test set (average error -0.11 ± 0.50 and +0.15 ± 0.58, respectively). The B3LYP functional also performs well (average error -1.11 ± 0.82 and -0.78 ± 0.79, respectively). ■ INTRODUCTION Substituted thiols have a wide variety of uses and applications in chemistry, biochemistry, and material chemistry. In biochem- istry, for example, thiols are known for antioxidant properties such as radical quenching. 1,2 In cell redox buffers, their role is to regulate the protein thiol/disulfide composition. The disulfide bonds are important in maintaining the structural stability of soluble proteins. 3,4 Some interesting examples from material science include the use of substituted benzenethiols in molecular electronics, surface-enhanced Raman spectroscopy, and quantum electronic tunneling between plasmonic nano- particle resonators. 5-11 Understanding of the properties and reactivities of thiols as a function of pH requires a reliable set of measured or calculated acid dissociation constants. The experimental determination of pK a ’s is not always easy because of problems such as interference from other solutes in the complex substrate environment, difficulties in isolation of specific residues, complexity due to the solvent system, etc. Hence, there is always a need to calculate pK a ’s using quantum chemical techniques. The calculation of pK a ’s is the subject of a number of recent reviews. 12-14 The pK a for a molecule is obtained from the solution phase free energy of the deprotonation reaction, AH ⇋ A - +H + . The quality of calculated pK a ’s depends on the accuracy of the computed deprotonation energies and the reliability of the estimated solvation energies. Early studies showed that some implicit solvation models can lead to large errors in the calculated pK a ’s. 12-14 However, these errors are often systematic for a given functional group, and suitable estimates of pK a ’s can be obtained from linear correlations between calculated solvation free energies or pK a ’s and known experimental values. Friesner and co-workers developed a protocol for predicting pK a ’s for a wide range of functional groups that involved a linear correlation between experimental pK a ’s and raw pK a ’s computed by using free energies from density functional calculations and a continuum solvation model with radii optimized for each functional group. 15 Zhang Received: May 18, 2016 Revised: June 20, 2016 Published: June 21, 2016 Article pubs.acs.org/JPCA © 2016 American Chemical Society 5726 DOI: 10.1021/acs.jpca.6b05040 J. Phys. Chem. A 2016, 120, 5726-5735