Assessment of Methodology and Chemical Group Dependences in the

Subscriber access provided by The Library | University of Bath

Article

Assessment of Methodology and Chemical Group Dependences in the Calculation of the pKa for Several Chemical Groups Toru Matsui, Yasuteru Shigeta, and Kenji Morihashi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00587 • Publication Date (Web): 01 Sep 2017 Downloaded from http://pubs.acs.org on September 3, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Assessment of Methodology and Chemical Group Dependences in the Calculation of the pKa for Several Chemical Groups Toru Matsui a, *, Yasuteru Shigeta b, and Kenji Morihashi a a

Department of Chemistry, Graduate School of Pure and Applied Sciences, University of Tsukuba, 1-1-

1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan. b

Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-

8577, Japan. * Corresponding author. E-mail: [email protected]. FAX +81-29-853-7392.

(ABSTRACT) We have investigated the dependencies of various computational methods in the calculation of acid dissociation constants (pKa values) of certain chemical groups found in protonatable amino acids based on our previous scheme [Matsui et al. Phys. Chem. Chem. Phys. 2012, 14, 4181-4187.]. By changing the quantum chemical (QC) method (Hartree-Fock (HF) and perturbation theory, and composite methods, or exchange-correlation functionals in density functional theory (DFT)), basis sets, solvation models, and the cavities used in the solvent models, we have exhaustively tested about 2,200 combinations to find the best combination for pKa estimation among them. Of the tested parameters, the choice of the basis set and cavity is the most crucial to reproduce experimental values compared to other factors. Concerning the basis set, the inclusion of diffuse functions is quite important for carboxyl, thiol, and phenol groups judging from the mean absolute errors (MAEs) measured from the experimental values. Of the cavity models, between the Pauling, Klamt, and the universal force field (UFF) definitions, the UFF defined cavity is the best choice, resulting in the smallest MAEs. Concerning the QC methods, hybrid DFTs and range-separated DFTs always provide better results than pure DFTs and 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 46

HF. As a result, we found that LC-ϖPBE/6-31+G(d) with PCM-SMD/UFF provides the best pKa estimation with a MAE of within 0.15 pKa units.

1. Introduction The acid dissociation constant (pKa value) determines whether a protonatable/deprotonatable chemical group in a protein is protonated at a certain pH, thus it is very important factor in both fields of chemistry and biology. Particularly for synthetic chemists, it is useful to know the theoretically predicted pKa values of target compound before synthesis for the rational design of functional materials. For both experimental and theoretical biochemists, information concerning the pKa values of the amino acids on a protein surface is of great importance for the preparation of buffer solutions to regulate the pH of experiments, and the theoretical estimation of the pKa is necessary to perform molecular dynamics (MD) simulations of proteins. In addition, the pKa value is of interest in pharmacology. For example, the order of the deprotonation process in polyprotic compounds such as salicylic acid, histamine, and dopamine is very relevant concerning their pharmacological action in the brain. Therefore, the development of a reliable pKa estimation scheme is an important issue in the quantum chemistry.1, 2 So far, many theoretical studies have reported their own computational schemes to evaluate the pKa values of specific compounds. These methods can be mainly categorized into three groups: (1) using MD simulations, (2) using quantum mechanics/molecular mechanics (QM/MM) MD simulations, and (3) using static QM model calculations with a continuum model to describe solvent effects. In MD simulations of proteins, PROPKA

3, 4

is a widely used program to evaluate the pKa of specific amino

acids, and is often used before performing MD simulations in the literatures. However, PROPKA sometimes gives unrealistic pKa values depending on the quality of the X-ray crystal structure of the chosen protein, and is not applicable to substrates, co-factors, and small molecules other than protein amino acids. QM-based methods are required to model the pKa more generally with high accuracy. The 2 ACS Paragon Plus Environment

Page 3 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


most realistic and accurate method to consider the effects of explicit solvent within these three categories is the QM/MM scheme. However, obtaining pKa values using QM/MM methods requires significant computational resources because of the evaluation of the Gibbs energy differences between HA and A– and because free energy calculations require knowledge of the many configurations of the water environemnt in the QM/MM simulation. Instead, many theoretical papers have used Car-Parinello MD simulations (CPMD), over a 1-ps MD simulation to evaluate the free energy difference of molecules, for example, the pKa value of a molecule in the condenced phase, even for small molecules such as amino acids.5-7 However, using these methods to simulate realistic cases is still time-consuming. Thus, a simple computational scheme is necessary to estimate the pKa values. The easiest alternative approach to avoid the sampling problems that arise from all-atom simulations is to carry out static QM calculations combined with an implicit solvent model.

8-17

This

method is simpler because only the geometry optimizations and free energy calculations based on the harmonic vibrational frequencies for the protonated and deprotonated states of a given target molecule are necessary to obtain the Gibbs energy difference in solution. The combinations of QM methods with polarizable continuum models (PCM) is popular, so many theoretical works using PCMs have adopted a thermo-chemical cycle scheme with a constant Gibbs energy of a proton, G(H+) to estimate the deprotonation free energy in a solution. Although some computational works have attempted to estimate methodology dependent G(H+) values via a proton exchange reaction (A– + H3O+ + (n-1)H2O-> HA + nH2O or A– + nH2O -> HA + OH- + (n-1)H2O) to remedy this problem, 18-20 the obtained values strongly depends on the number of water molecules, n. Moreover, when the choice of the basis set is poor, for example, when diffuse functions are not included, the energy difference between OH– and H2O is clearly different from the expected value (with respect to values calculated using a more complete basis set) even if solvent effects are taken into account as shown in our previous study.21 Alternatively, many theoretical works have employed a linear correlation between the computational free energy differences estimated by using thermochemical dynamics and experimental pKa values. Although this approach 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 46

resembles to our scheme in referring to the experimental pKa values, it overlooks G(H+). Whichever of these methodologies is chosen, methodology dependence is an inevitable problem in computing pKa values by QM calculations.22 Although the fundamental biochemistry textbooks always list the pKa values of polar amino acids, few computational studies have succeeded in reliably reproducing these pKa values because of the problems mentioned above. To avoid these problems, we have proposed the “appropriate treatment for pKa and proton’s energy with benchmark set” (AKB) scheme by using QM methods with implicit solvation models.23, 24 This scheme is also based on a linear regression curve between the Gibbs energy difference of HA and A– and the experimental pKa values of reference molecules to estimate pKa values of a target molecule, which is similar to methods used in analytical chemistry. In previous studies, we reported that we could estimate the pKa values of protonatable side chains in amino acids or small oligopeptides with high accuracy (within an error of 0.2–0.5 pKa units). This computational scheme also has an advantage in that it can theoretically estimate the optimal G(H+) for each chemical group depending on the method employed. Although G(H+) is, in principle a constant value in the actual experiments, the PCM may cause the over-/under-estimation of the deprotonation energies because it implicitly treats the solvent molecules around the solute so that reorganization of solvent molecules around the solute upon the deprotonation is not perfectly described. This computational scheme sometimes fails to describe highly charged systems because the PCM stabilizes the charges in aqueous solution excessively, as discussed later. Another problem in this scheme is the need to determine a linear regression curve not only for each chemical group but also for each computational method. (as shown seen in recent works25–27). Therefore, it is necessary to accumulate further data for pKa calculations. Although the methodology dependence of the pKa values is troublesome and has often been discussed in the field of computational chemistry, to date, there have been few reports of its systematic analysis, except for Ref. 28. In our previous studies, several computational methods were adopted, and, thus, the method dependencies of the computational results were not fully discussed. To find which computational method is reliable for 4 ACS Paragon Plus Environment

Page 5 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the computation of pKa values, we must assess the dependencies on the basis sets, the exchangecorrelation functionals used in density functional theory (DFT), and the solvent models in the PCM of the pKa values. By computing many combinations of these computational methods, here, we report the methodology dependencies in computing accurate pKa values and G(H+) for specific chemical groups. Section 2 gives a brief explanation of how we estimated the pKa values and G(H+). In Section 3, we numerically demonstrate the methodology dependencies. In Section 4, we compare the results estimated from the potential energy to those estimated from the Gibbs energy. In Section 5, we discuss which computational method is the best for practical use. Finally, we give a summary in Section 6.

2. Computational Details We introduce the method for computing the pKa values for the side chains of polar amino acids and for halogen-substituted pyrimidine bases used in our recent works. Here, we briefly review the computational method. The pKa value is obtained from the Gibbs energy difference ∆G(solv) of a deprotonation reaction (HA → A– + H+) in solution, i.e.,

pK a =

∆G (solv ) , ( ln10) RT

(1)

where R and T denote the gas constant and the absolute temperature, respectively. In this study, the temperature was fixed at 298.15 K (25 ºC). Recent PCM results enable us to compute the Gibbs energy in an implicit solvent. We assumed that the Gibbs energy difference in solution is

( ) ( )

∆G (solv ) = G A- + G H + − G ( HA ) ,

(2)

where G(A–) and G(HA) can be estimated directly from free energy calculations in solution from vibrational frequency analyses using a combination of a QM method and the PCM; however, one cannot directly estimate G(H+) from any quantum chemical calculation owing to the absence of an electron in H+, so a constant G(H+) value is generally taken from the literature. However, the obtained results are 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 46

not sufficiently accurate, as mentioned before. Moreover, even if one parametrizes G(H+) to fit the pKa of a given compound, Equation (1) is still worse for predicting the pKa judging from a correlation plot between the calculated ∆G(solv) and the experimental pKa values.23 In our AKB scheme, we introduce an effective scaling factor, s, to be multiplied together with the Gibbs energy difference between the reactant and the product. Then, we have

{ ( )

}

( )

sG H + s∆G (solv ) s G A − G ( HA ) . pK a = = + ( ln10) RT ( ln10) RT ( ln10) RT

(3)

The scaling factor s serves not only as the activity coefficient of the deprotonation reaction in solution but also as an error correction for a computational method employed. Because the difference between the computed and experimental pKa values depends on the computational methods, such as the exchange-correlation functionals in DFT, the basis sets, and the cavity model used in the PCM method, the error in Δ‫(ܩ‬solv) is expected to be systematic within similar reactions, i.e., the same chemical group. Thus, here, we regard G(H+) and s as constants, both of which depend on these approximation levels. Consequently, Equation (3) becomes

pK a = k∆G0 + C0 , s k= ( ln10) RT

and

C0 =

(4a)

( )

sG H +

( ln10) RT

,

(4b)

where ∆G0 is the Gibbs energy difference of G(A–)–G(HA). Equation (4) results in an apparent linear correlation between the ∆G0 and pKa values. Therefore, it is possible to compute an unknown pKa from a theoretical Gibbs energy difference by fitting appropriate values of k and C0 from ∆G0 and the experimental pKa values of known reference molecules. We have computed ∆G0 for several reference molecules and performed a least-square fitting with the corresponding experimental pKa values (see SI for details). Once these parameters are obtained by the least-square fitting, we can estimate G(H+) in aqueous solution from the slope and the abscissa as G(H+) = C0/k.

(5) 6

ACS Paragon Plus Environment

Page 7 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Thus, the AKB scheme has an advantage over other schemes, such as MD, QM/MM MD, and thermodynamic cycle methods, in obtaining G(H+) for each chemical group and each computational method. Thus, the AKB scheme is very effective for molecules whose optimized geometries and vibrational frequencies are easily available. Actually, many papers have reported that the linear fitting procedure enables the reproduction of the experimental pKa values of many compounds.29–36 In principle, the resultant scaling factor and G(H+) might be the same regardless of the chemical compound only when one adopts highly accurate QM treatment (for example, basis set limit results using coupled cluster methods, composite methods, and so on) with a reliable solvation environment. However, according to our previous studies using Hartree–Fock (HF), DFT, and second-order Møller–Plesset perturbation theory (MP2) calculations with limited-size basis sets, the obtained results depend strongly on the chemical groups, so k and C0 must be estimated for each chemical group. The G(H+) value obtained by this scheme seems to be variable because of the model without any solvent molecule. Some researches showed that the explicit solvent model gives the constant G(H+) value.37, 38 Such a group dependency remains a serious problem in our scheme, which must be solved when we discuss about the protonation reactions. We used 55 chemical compounds as an experimental data set, including 12 carboxyl acids, 9 primary amines, 9 thiols, 10 phenols, 8 imidazoles, and 7 alcohols in the polar amino acid residues, which are solvated in aqueous solution.39–50 The details (names of compounds and their experimental pKa values) are shown in the Supporting Information. All calculations were carried out using Gaussian 09 D. 01

51

except for Section 3.3. Throughout this study, the grid for DFT calculations was set to “ultrafine” (a pruned (99, 590) grid), and the numbers of d and f functions in the basis functions were set to 5 and 7, respectively.

3. Benchmark Calculations

7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 46

It is necessary to choose four options to carry out the solvent model calculations: (1) the DFT exchange-correlation functional in DFT calculations or the electron correlation method in ab initio methods, (2) the basis set, (3) the solvent model, and (4) the cavity radius. We employed 20 functionals/methods (BLYP, B3LYP,52 BHandHLYP,53 LC-BLYP,54 CAM-B3LYP,55 B2PLYP,56 PBE, PBE1PBE (namely PBE0),57 LC-ϖPBE,58, 59 B97D,60, 61 ϖB97XD,62 PW91, B3PW91, M06, M06-L, M06-2X,63, 64 MN12SX,65 TPSSh,66, 67 HF, and MP2), 12 kinds of the well-known basis sets (6-31G, 631G(d), 6-31G(d,p), 6-31+G(d), 6-31+G(d,p), 6-31++G(d,p), 6-311++G(d,p), cc-pVDZ, aug-cc-pVDZ, double zeta with the Stuttgart Dresden effective core potential (ECP) (SDD),68 LanL2DZ,69,

70

and

SVP). We also employed the “model chemistry calculations” provided in Gaussian (CBS-4M, CBSQB3,71–73 G3MP2, G3, G3B3,74–76 and G4MP277), which give good agreement with the experimental values shown in our previous studies. For each combination, three kinds of the solvent models (IEFPCM (PCM),78 C-PCM,79–81 and PCM-SMD (SMD)82) and three kinds of parameters for the cavity radii (universal force field (UFF), Klamt, and Pauling) were adopted to model the solvent effects. The notation for a combination of these options is represented as “method/basis with model/cavity” for simplicity. In this paper, the accuracy of a given combination is quantitatively discussed in terms of the average of mean absolute errors (MAE) for six chemical groups, noted as “average MAE.” At first, the results of 2,214 combinations of ∆G0 calculations are plotted with the experimental pKa values, as shown in Figure 1. Note here that the LC-ϖPBE/6-31+G(d) with PCM-SMD/UFF method provides the best average MAE, 0.15 pKa units, of the 2,214 computational methods. As shown in Figures 1a–1c, two compounds with different functional groups (for example, thiol and phenol) but similar pKa values do not always yield similar deprotonation energies. Let us consider the cases of thiol, amine, and phenol groups. As shown in Figure 1b, the plot for thiols is similar to that for phenols, whereas that for amines differs from the thiol and phenol plot. On the other hand, as shown in Figure 1c, the plot for thiols almost overlaps that for amines. With ab initio electron correlation calculations, the scaling factor depends only on the proton donor (-NH, -OH, and -SH), as shown in Figure 1d, which 8 ACS Paragon Plus Environment

Page 9 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shows that the slope of amines is almost the same as that of imidazoles. The details of these topics will be discussed later. This trend is also observed when one adopts the G3 method, which is a quantum chemistry composite method. Interestingly, the value of G(H+) determined using G3 with SMD/Pauling was almost independent of the chemical group, as shown in Figure 1f (see the horizontal intercepts for six chemical groups in Figure 1f). As shown in this figure, the linear regression curves for -OH (-SH) and -COOH groups (carboxylic acids, thiols, phenols, and alcohol) are in good agreement with each other, and that is also true for -NH groups (amines and imidazoles), although those slopes are slightly different. Therefore, it is expected that a highly accurate method along with the proper choice of solvation models may give the same slope and intercept, as mentioned in the previous section. Nevertheless, this is sometimes impractical because of the high computational demand.

3.1. Basis set dependence Because the computational methods/theories exhibit similar trends in the linear correlation between the calculated ∆G(solv) and experimental p‫ܭ‬௔ values, the “B3LYP with C-PCM/UFF” method was used for simplicity to extract the basis set dependence of the predicted p‫ܭ‬௔ values. The MAE for each basis set is plotted in Figure 2, and the basis set dependencies on the slopes and intercepts are shown in Table 1. Let us discuss the individual basis set dependencies in detail. The double-zeta basis sets with diffuse functions and polarization functions provide good results, within 0.2 pKa units, in the average MAEs. Surprisingly, the addition of diffuse and polarization functions to hydrogen atoms does not significantly affect the pKa estimation. The MAEs obtained at the 6-31+G(d) are similar to those obtained using the 6-31+G(d,p) and 6-31++G(d,p) basis sets. Although the MAEs of the calculations using a triple-zeta basis set are improved in comparison to those using double-zeta basis sets, the inclusion of diffuse and polarization functions may further improve the MAEs, and the 6-31+G(d) basis set is highly recommended for our method because of the balance between computational cost and accuracy. This 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 46

trend is also true for the most combinations of methods, except for MP2. The MP2 calculations require the aug-cc-pVDZ basis set rather than the Pople-type and ECP basis sets. The basis set dependence of MP2 is discussed later. In the case of the OH group, some basis sets without diffuse functions give better results than those with diffuse functions. However, the pKa for H2O is not close to the experimental value (15.74) without the use of diffuse functions, where the deprotonation energies of H2O are clearly overestimated, and the computational pKa values range from 19 to 22 (the experimental value is 15.74). In contrast, by using diffuse functions, the computational pKa become similar to the experimental value (14.9–16.1). This is probably because the electronic structure of a single OH- cannot be adequately described with a small basis set because the negative charge is delocalized in the deprotonated state and diffuse functions are necessary to describe the electronic structure of the anionic state. As a result, the MAE becomes larger when diffuse functions are not included. On the other hand, the MAEs, which range from 0.26 to 0.32, of the amine and imidazole groups weakly depend on the basis set, except for the amines calculated with the 6-31+G(d) basis set. In the case of these groups, the deprotonation process changes the charge of the compound from +1 to 0. It is unnecessary to describe the anionic state in these compounds. In the case of thiols, large all-electron basis sets are required to achieve an error within 0.2 pKa units. Oddly, ECPs also give good results, especially for the thiols. Because LanL2DZ places an ECP on the sulfur atom, the efficiency of the ECP can be seen in the difference between the LanL2DZ and SDD ECP. However, based on our results, there is little difference between them. Therefore, it is effective to use “all LanL2DZ” calculations for the prediction of thiol pKa values. However, the performance of the ECP on the imidazoles is not as good as for other chemical groups.

3.2. Electron correlation effects in ab initio methods and exchange-correlation functionals in DFT In this subsection, the exchange-correlation functional dependence is discussed. For simplicity, the “6-31+G(d) with C-PCM/UFF” method was used if the method is not explicitly stated. In Table 2, the 10 ACS Paragon Plus Environment

Page 11 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


scaling factors and estimated G(H+) for each chemical group are listed, and Figure 3 illustrates MAEs for twenty different levels of theory. Overall, the scaling factor does not depend on the computational method but depends strongly on the chemical groups. HF and MP2. As expected, MP2 generally gives better results than HF in almost all chemical groups if we use a large basis set, although the results are comparable for the 6-31+G(d) basis set, as shown in Figure 3. MP2 calculations with the 6-31+G(d) give minimum errors for the imidazole groups, whereas DFT calculations do not give equivalent results. However, the MP2 calculations give worse results for the carboxyl and amine groups than the DFT calculations. The best result of the MP2 calculations carried out in this study was obtained with the “MP2/aug-cc-pVDZ with SMD-PCM/UFF” method, which yields an MAE of 0.16 pKa units for all compounds. MP2 calculations are strongly basis set dependent, and, thus, this result clearly shows the dependence too, suggesting that larger basis sets give better results, especially for MP2 calculations. However, the calculation of the Gibbs energy using MP2 is computationally expensive. In this regard, DFT is more suitable for the estimation of the pKa values because of the balance between computational cost and accuracy. In the case of the HF calculations, the results for the thiols are the worst among all the calculations, probably because of the lack of dynamical electronic correlation effects in describing the electronic structures of heavier elements. Pure and hybrid-DFT functionals. The average MAEs determined using pure DFT functionals (BLYP, PBE, and PW91), whose MAEs are no less than 0.2 pKa units, are similar to those obtained by HF. The MAEs of the carboxyl and thiol groups determined using any pure DFT functional are smaller than those determined by HF, whereas those of the phenol and alcohol groups are larger. This error probably originates from the fact that pure DFT functionals and HF are poor at describing the electronic structures of charged aromatic molecules. A similar trend (in scaling factor and G(H+)) can be observed for BLYP, PW91, and PBE, except for the calculation of the amine group by BLYP. The result for M06L is a little different from the other pure DFT functionals because it is a meta-GGA functional, as



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 46

discussed later. As seen in the results using B97D, we estimated a large G(H+) in the phenol group compared with the other pure DFTs. For hybrid DFT functionals, the percentage of added HF exchange might determine the average MAE. In this case, the HF exchange seems to compensate the errors by using pure DFTs. The addition of 20– 25% exchange provides the best performance, as found in B3LYP, B3PW91, and PBE0 functionals. The calculation “B3PW91/6-31+G(d) with PCM-SMD/UFF” gave the second-best result, with an average MAE of 0.15 pKa units. When the exchange ratio is increased to 50%, the results are very similar to those obtained by HF. The calculations using B2PLYP, which is a double hybrid functional, are expected to be better than those using B3LYP because of the inclusion of the exact MP2 correlation term. However, the average MAE of B2PLYP is a little worse than that of B3LYP, and a strong basis set dependence was also found in the case of B2PLYP, as for MP2. Therefore, the trend is almost the same as MP2. Because the computational cost of the B2PLYP calculation is higher than those of both MP2 and B3LYP in Gaussian 09, B2PLYP is not recommended for the pKa calculations. Meta-GGA (Minnesota family). In this study, four of Truhlar’s Minnesota exchange-correlation functionals (M06, M06-L, M06-2X, and MN12SX) were examined as typical examples of meta-GGA functionals. Some papers have pointed out that the Minnesota family with regular grids sometimes fail to describe potential energy surfaces smoothly.83, 84 However, when using an ultrafine grid, this problem is resolved. The structural convergence for the thiol groups is not good using the Minnesota family. The M06 and M06-L functionals are categorized as pure functionals, and these give similar results to those of the pure GGA functionals discussed previously. On including HF exchange, the MAE becomes smaller, as confirmed with M06-2X or MN12SX. In particular, M06-2X gives a rather good result for the amine and imidazole groups. Among the Minnesota family, “MN12SX or M06-2X/6-31+G(d) with C-PCM/UFF” are also recommended for the calculation of pKa.


Page 13 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Long-range corrected DFT functionals. The LC-BLYP, CAM-B3LYP, LC-ϖPBE, and ϖB97X-D exchange-correlation functionals are “range-separated” (RS) or "long-range corrected” (LC) DFT functionals. It has been reported that LC-DFT with dispersion correction provides better results for the interaction energies than MP2.85, 86 As expected, the LC-DFT functionals typically gave the best results among all the DFT functionals used. In addition, the LC-DFT functional adequately describes the charge separation and partly includes the self-interaction energy correction. Therefore, LC-DFTs are expected to be promising methods for pKa prediction. Indeed, LC-BLYP provides the best MAE for the carboxyl and phenol groups among all calculation methods, and LC-ϖPBE gives the best averaged MAE among the DFT functionals. As reported in the next subsection, we employed LC-ϖPBE/6-31+G(d) to assess the solvent model dependence. Composite model chemistry. According to our previous study, composite methods provide very accurate results for the alcohol group compared with the gold standard CCSD(T,AE)/aug-cc-pVTZ calculations.21 In this subsection, we describe a series of composite model calculations that we applied to the six chemical groups to show the limitation of using linear regression curves. The MAEs for the six composite model calculations are shown in Figure 4 to clarify the solvent and cavity model dependencies. The average MAEs of CBS-4M, CBS-QB3, and G4MP2 methods are not as good as those from the HF calculations, at best 0.2 pKa units. The G3 family (G3, G3B3, and G3MP2) gave slightly better results than former composite model methods. As a result, G3B3 gives the best results of the six composite model methods. Unfortunately, the composite model methods do not always give the best result. For example, G3B3 with PCM-SMD/UFF, which gave the best result among 54 combinations of the composite models, is worse than the best of the DFT calculations. The largest error is found in the calculation of the imidazole group, as is the case with any other calculations. Therefore, the correction for the imidazole groups should be considered in further refinement for accurate pKa prediction within our scheme.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 46

3.3. Solvent model For the PCM, many theories have been proposed: conductor-like PCM (C-PCM), integral equation formalism variant PCM (IEF-PCM), and the PCM-solvation model based on density (SMD), and these have been implemented in many QM programs. Independent of the solvent models, there are many cavity models that are used to place a given molecule in the solvent continuum. In this study, three cavity models (the universal force field model (UFF), Klamt model, and Pauling model) were examined. Our previous studies have shown that the united atom models for Kohn-Sham DFT calculations (UAKS) model also is useful in computing the pKa values. However, because UAKS often fails in the geometry optimization of several compounds, we did not use this scheme. Thus, this study examines nine combinations of “solvent model/cavity model” at the same level of theory (“theory/basis set”), where the LC-DFT and composite model chemistry methods were adopted. Note here that the combination of “IEF-PCM with UFF” is the default choice in self-consistent reaction field calculations in Gaussian09. Model dependency. Figure 5 shows the MAEs for each "solvent model/cavity model" with “LCϖPBE/6-31+G(d)”, and Table 3 lists the obtained scaling factor and G(H+). The solvent model does not affect these values significantly but does affect the magnitude of the MAEs. The average MAEs are in the order of PCM-SMD < C-PCM < IEF-PCM, although the MAEs strongly depend on the chemical groups. On the other hand, both the scaling factor and G(H+) depends on the choice of the cavity model. Particularly, G(H+) of the amine group obtained by the Pauling model differs that from the UFF model by more than 80 kJ/mol. This remarkable cavity model dependence also appears in the calculation of imidazole, phenol, and alcohol groups, while no apparent changes were observed in calculations of the carboxyl and thiol pKa values. Although the Pauling model gives the large scaling factor of 0.8 for the imidazole and amine groups, it does not give accurate pKa values, whereas the MAEs for the amine and imidazole groups are over 0.5 pKa units with the Pauling model. The Klamt model exhibits just a


Page 15 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


moderate trend between the UFF and Pauling models. The accuracy of the cavity models is found to be in the order of UFF > Klamt > Pauling. Cavity Radii Parameter and Gaussian version dependency. In many studies which employ PCM as a solvation model, the default parameters determined by Gaussian are frequently used. On the other hand, the parameter α, which determines the cavity radii using in PCM, is important for the solvation Gibbs energy. Some studies investigated the dependencies of α on pKa values and found that α=1.2 gives better results than α=1.1 (default).87 The version of Gaussian is also important in such a benchmark study as investigaed for influences to solvation energy by Takano and Houk.20 Here, we again investigated the alpha dependencies with several versions of Gaussian, i.e., Gaussian 09 E.0188 and Gaussian 16 A.03. 89 For this purpose, M06-2X/6-31+G(d,p) with PCM-SMD/UFF is employed. The α dependency with several versions of Gaussian is presented in Figure 6. As the version is newer, the MAE becomes smaller in this case. The latest Gaussian gives the least MAE (0.16 pH units in 56 compounds). In any versions, the default value (α=1.1) gives the least MAE in several chemical groups. Although the default value works well in this case, the selection of it can influence the result to some extent.

4. Obtaining the pKa from the Potential Energy Differences In our previous study,24 we proposed a simple method of pKa estimation using the energy of the optimized geometry (SKE) scheme, which employs the potential energy difference between the protonated and deprotonated compounds rather than the difference in the Gibbs energies. The Gibbs energy of a molecule X, G(X), is divided into three parts; (i) the potential energy at the optimized geometry ESCF(X), (ii) the zero-point energy EZPE(X), and (iii) the thermochemical term ∆G0-> 298K(X):

G ( X ) = ESCF ( X) + EZPE ( X) + ∆G0→298K ( X) .

(6)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 46

The SKE scheme assumes that (ii) and (iii) depend on the chemical group and are approximated as constant values for each chemical group. When the constant parts are summed into another constant C’0, the pKa can be rewritten as

pK a =

{ ( )

}+

s' E A - − E ( HA )

( ln10) RT

s'Grest = k ' ∆E0 + C '0 , ( ln10) RT

(7)

where Grest is a sum of residual parts of free energy components defined as

( )

( )

( )

Grest = G H + + EZPE ( HA ) + ∆G0→298K ( HA ) − EZPE A - − ∆G0→298K A - .

(8)

Although G(H+) is disregarded within this scheme, all we require for the pKa calculation with the SKE scheme is the energy difference between the protonated and deprotonated states. This treatment considerably reduced the computational costs because it is not necessary to evaluate a Hessian matrix (second derivative with respect to the nuclear coordinates) for the estimation of vibrational frequencies, i.e., zero-point energy and entropy corrections. Now, we compare the SKE scheme with the AKB scheme (using the Gibbs energy) as shown in Figure 7, where the dotted lines indicate that the MAE provided by the former is the same as that of the latter. We used the B3LYP functional and explored 117 kinds of the other combinations (thirteen basis sets, three solvation models, and three cavity models). Judging from Figures 7a–7f, the SKE scheme yielded better MAEs than the AKB scheme for many computational methods, except for those concerning the imidazole group. As shown in Figure 7g, the SKE scheme sometimes exceeds the AKB scheme in some combinations of the methods. On the other hand, the AKB scheme gives the best average MAE of 0.15 pKa units, while the SKE scheme yields a worse average MAE of 0.18 pKa units. Therefore, the SKE scheme provides the pKa value within the reasonable error despite the low computational cost and, thus, serves as a convenient alternative to the AKB scheme.

5. Discussion 16 ACS Paragon Plus Environment

Page 17 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This section discusses which computational method is the best for pKa estimation. Table 4 lists the best method (the lowest MAE) for each chemical group. In some groups, computation with a small basis set provided the best result. We excluded these results because such a basis set does not accurately evaluate G(H+) or the pKa of H2O in the alcohol group. The “Best method” noted in the table shows the best results of 2,214 combinations of the computational methods for each group. When the accuracy of a computational method is evaluated as the average over six chemical groups, the calculation “LCϖPBE/6-31+G(d) with SMD/UFF” noted as the “Recommended method” in Table 4 provides the least average MAE 0.14 pKa units. Each chemical group shows a unique trend, in other words, characteristics distinguishing them from the other groups. The carboxyl group seems to depend strongly on the basis set, whether the basis set contains diffuse functions or not. On the other hand, the MAE of the phenol and alcohol groups does not depend on either the basis set or the computational model. However, the scaling factor of the phenol group is smaller than the other two groups (the carboxyl and alcohol groups), and G(H+) is underestimated. In the cases of the alcohol and phenol groups, the post-HF (MP2) and composite model chemistry methods provide better results than those from DFT and HF calculations. The recommended calculation presents the MAE of 0.03 pKa units for the alcohol group. This result is rather better than the gold standard calculations, CCSD(T)/CBS-limit, presented in our previous study.35 Actually, the “G3B3 with CPCM/Klamt” calculations gives the smallest MAE (0.06 pH units) for the phenol group. The trend for the thiol groups is similar to that for the carboxyl groups because the solvent (cavity) model dependency is not strong. However, unlike for the carboxy group, the best result is obtained with BHandHLYP, which is in-between pure DFT and HF. In the case of the amine and imidazole groups, the cavity model dependence is quite strong, as described in Section 3.3. For these groups, the “Pauling” cavity model should not be used. Because the basis set dependence is weak for the amine group, the basis sets without diffuse and polarization functions are sometimes better than those with them. However, for consistency with the calculations of 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 46

other chemical groups, it is desirable to include diffuse and polarization functions, even though there is little change on their inclusion. When diffuse and polarization functions are included, range-separated DFT calculations reduce the MAE for amine group pKa values. The imidazole group is the most complicated group among the six chemical groups. The basis set dependence exhibits similar trend to that of the amine. On the other hand, the error becomes larger when the diffuse function is included in the calculation. As noted before, the MP2/aug-cc-pVDZ calculation is the best for the imidazoles, though MP2 does not perform well in the calculations involving carboxylic acids and amines. In total, the number of the computational methods whose “average MAE” was smaller than 0.2 pKa units is quite limited. Among the 2,214 combinations, there are only five combinations, ((1) “LCϖPBE/6-31+G(d) with SMD/UFF,” (2) “B3LYP/6-31+G(d,p) with SMD/UFF,” (3) “B3PW91/631+G(d) with C-PCM/UFF,” (4) “ϖB97X-D/6-31+G(d) with PCM/UFF,” and (5) “M06-2X/6-31+G(d) with C-PCM/UFF”), whose “average MAE” is about 0.15 pKa units. Considering that the MAE determined by the composite model chemistry is over 0.17 pKa units, these calculations are very effective for pKa calculations on large molecules. Overall, the significance in the computational dependencies on the pKa calculations decrease in importance from (a) the cavity model, (b) the basis set, (c) the computational method (DFT functional, post-SCF, and composite model chemistry method), to (d) the continuum model. In particular, for (a), reliable calculations can be obtained by using the UFF cavity model. Finally, we applied these methods mentioned above to determine the pKa values of the polar amino acids listed in Table 5. The recommended method reproduced the pKa value of histidine (imidazole) well but overestimated the pKa values of cysteine (thiol), tyrosine (phenol), and lysine (amine). As discussed in our previous studies, this overestimation is due to the localized charge concentration, i.e., the amino acids are solvated in their zwitterionic forms.22 To improve these results, we must consider a larger molecular system such as tripeptides or use a counterion to neutralize the system.


Page 19 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Summary We have investigated several methodology dependencies in the estimation of pKa for specific chemical groups in polar amino acids. The dominant factor is the choice of cavity model. The basis set is also an important factor for many chemical groups. Not only the pKa values but also G(H+) are dependent on the chemical identity and are not constant, except for certain calculations. Among the DFT functionals and MP2 calculations, the hybrid and long-range corrected DFT functionals are reliable for the estimation of pKa because the errors by such calculations are comparable to those obtained by very accurate composite model calculations, i.e., G3 and/or G3B3. The best calculation is “LC-ϖPBE/631+G(d) with PCM-SMD/UFF,” which gave an MAE of 0.14 pKa units for 55 chemical compounds. There are many combinations of calculation methods other than those chosen in this study, but our dataset will be a great help for the pKa estimation, not only for polar amino acids but also for many other chemical compounds.

ACKNOWLEDGMENTS This study was supported by Grants-in-Aid for Innovative Areas "3D-active site science, Dynamical ordering, and Photosynergetics"(Nos. JP26102525, JP25104716 and JP26107004), and for Scientific Research (B) (Nos. 23350064, and 17H03034) from the Japanese Society for the Promotion of Science (JSPS). Some calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS).

SUPPORTING INFORMATION



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 46

Supporting information includes the detailed data of the experimental data for pKa values and all calculated results (2,214 combinations of computational methods). These materials are available on the World Wide Web at http://pubs.acs.org. FIGURE CAPTIONS Figure 1: Computed deprotonation Gibbs energies (ΔG) vs. experimental pKa values. (a) BLYP/LanL2DZ with C-PCM/Klamt, (b) LC-ϖPBE/6-31+G(d) with PCM-SMD/UFF, which provides the best results for total MAE, (c) HF/6-31G(d) with IEF-PCM/UFF, (d) MP2/aug-cc-pVDZ with PCMSMD/UFF, (e) G3 with PCM-SMD/Klamt, and (f) G3 with SMD/Pauling.

Figure 2: Basis dependence on mean absolute error (MAE) for each chemical compound. 631G(631G), 6-31G(d) (denoted “d”), 6-31G(d,p) (dp), 6-31+G(d) (pd), 6-31+G(d,p) (pdp), 631++G(d,p) (ppdp), 6-311++G(d,p) (6311), cc-pVDZ (cVDZ), aug-cc-pVDZ (acDZ), SDD(SDD), LanL2DZ (L2D), and SVP (SVP). The computational method, except for basis set, is fixed to B3LYP with PCM-SMD/Klamt.

Figure 3: Electronic structure calculation level dependence on the mean absolute error (MAE). 631+G(d) with CPCM/UFF. The used functionals (aberrations appeared in the figure) are BLYP (BLY), B3LYP (B3L), BHandHLYP (BHL), B2PYLP (B2P), LC-BLYP (LCB), CAM-B3LYP (CAM), HF (HF), MP2 (MP2), PBE (PBE), PBE0 (PB0), LC-ϖPBE (LCw), PW91 (PW91), B3PW91(B3P), TPSSh (TSh), B97D (97D), ϖB97XD (wXD), M06 (M06), M06-L (M6L), M06-2X (M62), and MN12SX (MNS).


Page 21 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: Composite model dependence on the mean absolute error (MAE). (a) CBS-4M, (b) CBSQB3, (c) G3, (d) G3B3, (e) G3MP2, and (f) G4MP2. Concerning the solvation model, CP, IP, and SM denote C-PCM, IEF-PCM (PCM), and PCM-SMD, respectively. For cavity radii, UF, KL, and PL denote UFF, Klamt, and Pauling, respectively. Figure 5: PCM and cavity model dependencies on the mean absolute error (MAE). LC-ϖPBE/631+G(d). The notation is the same as Figure 4. Figure 6: The dependency of the scaling factor (Alpha-value) or version of Gaussian in PCM on MAE in each chemical compound at M06-2X/6-31+G(d,p) level of theory with SMD-PCM/UFF. (a) Gaussian 09, D.01, (b) Gaussian 09, E.01, and (c) Gaussian16, A.03. The default value of alpha is 1.1 in any verions of Gaussian listed here. Figure 7: Mean absolute error (MAE) in each chemical compound. The horizontal axis is for MAE from AKB scheme (Gibbs energy), and the vertical axis is for that from SKE scheme (potential energy). (a) carboxyl, (b) amine, (c) thiol, (d) phenol, (e) imidazole, (f) alcohol, and (g) overall MAE. The line in the figure shows that the MAE provided by two schemes are the same. We employed the B3LYP functional for the calculations.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 46

Figure 1 Matsui et al.

Figure 2 Matsui et al. 22 ACS Paragon Plus Environment

Page 23 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60





1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 46



Page 25 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60






1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 46



Page 27 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


TABLE 1 Scaling factor and G(H+) (in kJ/mol) in each chemical group (B3LYP with PCM-SMD/Klamt). Thiol Phenol Imidazole Alcohol Carboxyl Amine s

G(H+)

s

G(H+)

s

G(H+)

G(H+)

s

s

G(H+)

s

G(H+)

6-31G

0.343 -1129.20

0.592

-1102.64 0.403

-1057.44 0.211

-989.91

0.584 -1129.45

0.300

-1048.33

6-31G(d)

0.386 -1156.79

0.638

-1095.88 0.566

-1130.97 0.290

-1066.25 0.580 -1110.14

0.350

-1098.05

6-31G(d,p)

0.390 -1171.68

0.643

-1103.74 0.561

-1137.84 0.250

-1053.21 0.578 -1119.41

0.352

-1113.14

6-31+G(d)

0.429 -1102.38

0.728

-1088.29 0.549

-1101.30 0.285

-1014.15 0.582 -1086.16

0.429

-1066.97

6-31+G(d,p)

0.429 -1116.17

0.668

-1086.04 0.569

-1112.26 0.285

-1028.41 0.584 -1095.11

0.431

-1082.51

6-31++G(d,p)

0.429 -1116.33

0.671

-1086.40 0.559

-1110.77 0.287

-1029.99 0.584 -1095.26

0.433

-1083.45

6-311++G(d,p)

0.416 -1114.88

0.667

-1083.23 0.569

-1112.85 0.289

-1028.73 0.590 -1093.44

0.439

-1086.24

cc-pVDZ

0.424 -1170.07

0.622

-1093.21 0.608

-1156.25 0.263

-1049.47 0.550 -1105.83

0.389

-1125.49

aug-cc-pVDZ

0.431 -1117.86

0.673

-1082.32 0.574

-1110.83 0.286

-1028.05 0.629 -1095.83

0.442

-1086.56

SDD

0.383 -1110.66

0.600

-1096.45 0.380

-1063.40 0.223

-979.69

0.579 -1116.22

0.307

-1027.22

LanL2DZ

0.378 -1109.72

0.600

-1096.56 0.400

-1064.20 0.223

-979.41

0.599 -1118.25

0.337

-1051.45

SVP

0.460 -1169.66

0.652

-1097.77 0.563

-1170.04 0.255

-1042.40 0.564 -1110.25

0.436

-1146.53

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 28 of 46

TABLE 2 Scaling factor and G(H+) (in kJ/mol) in each chemical group (6-31+G(d) with C-PCM/UFF) COOH Amine Thiol Phenol Imidazole Alcohol s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

BLYP

0.348

-1091.43 0.798

-1105.63

0.587 -1108.62

0.258

-985.70

0.559 -1069.90 0.385

-1044.13

B3LYP

0.359

-1100.78

0.551

-1039.33 0.602 -1112.60

0.299

-1023.40 0.567 -1075.50 0.397

-1063.68

BHandHLYP

0.359

-1108.60

0.565

-1045.76 0.533 -1101.62

0.321

-1047.38 0.549 -1079.16 0.393

-1076.27

B2PLYP

0.372

-1097.10 0.553

-1037.48 0.556 -1100.40

0.306

-1023.08 0.578 -1069.01 0.404

-1063.92

LC-BLYP

0.387

-1091.99 0.568

-1029.60 0.602 -1099.22 0.330

-1036.43 0.571 -1067.98 0.415

-1067.09

CAM-B3LYP

0.365

-1096.94 0.585

-1041.09 0.608 -1107.77

0.318

-1033.89 0.570 -1072.27 0.408

-1067.72

HF

0.352

-1119.51

0.574

-1055.07 0.512 -1096.12 0.335

-1067.97 0.541 -1087.07 0.384

-1087.88

MP2

0.348

-1085.95 0.546

-1034.68 0.561 -1095.96 0.340

-1036.52 0.603 -1059.09 0.414

-1064.25

PBE

0.366

-1095.52 0.545

-1034.39 0.602 -1107.23

0.255

-983.61

0.561 -1068.00 0.391

-1049.15

PBE0

0.364

-1105.69

0.584

-1048.39 0.610 -1111.94

0.269

-1007.88 0.554 -1074.15 0.400

-1070.81

LC-wPBE

0.387

-1112.12

0.575

-1049.64 0.620 -1120.48

0.327

-1051.40 0.571 -1083.24 0.418

-1087.49

PW91

0.353

-1094.26 0.543

-1034.54 0.598 -1108.09

0.258

-986.39

0.560 -1069.21 0.389

-1049.06

B3PW91

0.354

-1107.26

-1052.20 0.604 -1116.67

0.304

-1031.77 0.553 -1078.57 0.395

-1070.36

TPSSh

0.356

-1106.97 0.545

-1047.10 0.597 -1118.45 0.259

-1003.81 0.556 -1080.16 0.387

-1065.19

0.583

28


Page 29 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


B97-D

0.355

-1108.46 0.545

-1049.36 0.571 -1120.04 0.257

-1001.90 0.557 -1083.69 0.385

-1061.90

wB97XD

0.355

-1108.58 0.602

-1058.79 0.613 -1116.93 0.320

-1046.93 0.559 -1083.80 0.409

-1081.48

M06

0.393

-1105.92 0.556

-1033.49 0.631 -1111.83 0.267

-1001.79 0.558 -1065.81 0.387

-1059.21

M06L

0.364

-1114.06 0.572

-1046.01 0.678 -1126.13 0.246

-993.10

0.579 -1075.81 0.376

-1059.63

M06-2X

0.376

-1105.59 0.585

-1039.31 0.591 -1098.36 0.337

-1044.64 0.561 -1064.36 0.401

-1070.63

MN12SX

0.366

-1098.02 0.578

-1035.86 0.658 -1102.84 0.290

-1013.09 0.527 -1058.35 0.383

-1051.14

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 30 of 46

TABLE 3 Scaling factor and G(H+) (in kJ/mol) in each chemical group (LC-wPBE/6-31+G(d)) Carboxyl Amine Thiol Phenol Imidazole Alcohol s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

s

G(H+)

CPCM/UFF

0.387 -1112.12

0.575 -1049.64

0.620 -1120.48

0.327 -1051.40

0.571 -1083.24

0.418 -1087.49

CPCM/Klamt

0.458 -1112.16

0.695 -1099.48

0.620 -1120.42

0.384 -1069.98

0.591 -1103.38

0.471 -1096.22

CPCM/Pauling

0.555 -1111.54

0.839 -1118.20

0.646 -1117.00

0.499 -1092.10

0.733 -1123.14

0.554 -1104.40

PCM/UFF

0.385 -1112.29

0.581 -1050.83

0.614 -1120.03

0.340 -1056.87

0.575 -1083.67

0.417 -1086.83

PCM/Klamt

0.450 -1112.05

0.688 -1098.70

0.614 -1119.93

0.383 -1069.75

0.598 -1104.07

0.460 -1091.80

PCM/Pauling

0.505 -1107.84

0.826 -1117.12

0.648 -1117.38

0.482 -1089.35

0.625 -1117.32

0.554 -1104.66

SMD/UFF

0.394 -1112.43

0.564 -1039.74

0.575 -1110.24

0.342 -1058.39

0.564 -1072.90

0.421 -1087.23

SMD/Klamt

0.464 -1112.66

0.684 -1088.86

0.583 -1111.48

0.328 -1048.59

0.602 -1094.53

0.441 -1082.86

SMD/Pauling

0.559 -1110.76

0.865 -1111.97

0.672 -1117.71

0.392 -1065.29

0.668 -1110.58

0.521 -1094.01

30


Page 31 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


TABLE 4 Parameters and results of recommended method and the best computational method. Best Recommended Method

Basis

Model

Radii

Carboxyl

PW91

6-31+G(d)

PCM-SMD Klamt

0.0809

-89.02

0.08

0.0691

-76.83

0.19

Amine

wB97XD

6-31+G(d,p)

PCM-SMD UFF

0.0988

-103.58

0.10

0.0989

-102.85

0.10

Thiol

BHandHLYP

6-31++G(d,p)

IEF-PCM

UFF

0.1082

-121.84

0.12

0.1008

-111.92

0.16

Phenol

LC-BLYP

6-311++G(d,p)

CPCM

UFF

0.0619

-65.32

0.06

0.0600

-63.50

0.13

Imidazole

MP2

6-31++G(d,p)

PCM-SMD Klamt

0.1042

-113.20

0.13

0.0989

-106.09

0.23

OH

MP2

aug-cc-pVDZ

IEF-PCM

0.1084

-121.31

0.03

0.0738

-80.20

0.05

Pauling

k

C0

MAE

k

C0

MAE

31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 32 of 46

TABLE 5 Computed deprotonation energy (in kJ/mol) and estimated pKa values of each amino acid with “Best” and “Recommended” computational methods. Best

Recommended

Exp.

∆G

pKa

∆G

pKa

pKa

Asp

1159.84

3.68

1170.03

3.98

3.86

Glu

1166.71

4.14

1175.15

4.33

4.25

Lys

1162.66

11.28

1153.94

11.30

10.79

Cys

1207.26

8.82

1198.61

8.91

8.33

Tyr

1227.37

10.61

1232.55

10.45

10.05

His

1149.77

6.60

1135.26

6.17

6.04

H2O

1271.50

15.62

1296.41

15.43

15.74

32


Page 33 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES (1) Corona-Avendaño, S.; Alarcón-Angeles, G.; Rosquete-Pina, G. A.; Rojas-Hernández, A.; Gutierrez, A.; Ramírez-Silva, M. T.; Romero-Romo, M.; Palomar-Pardavé, M. New Insights on the Nature of the Chemical Species Involved during the Process of Dopamine Deprotonation in Aqueous Solution: Theoretical and Experimental Study. J. Phys. Chem. B 2007, 111, 1640-1647. (2) Baba, T.; Matsui, T.; Kamiya, K.; Nakano, M.; Shigeta, Y. A Density Functional Study on pKa of Small Polyprotic Molecule. Int. J. Quantum Chem. 2014, 114, 1128-1134. (3) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theor. Comput. 2011, 7, 525-537. (4) Søndergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J. Chem. Theor. Comput. 2011, 7, 2284-2295. (5) Cheng, J.; Sulpizi, M.; Sprik, M. Redox Potentials and pKa for Benzoquinone from Density Functional Theory Based Molecular Dynamics. J. Chem. Phys. 2009, 131, 154504. (6) Mangold, M.; Rolland, L.; Costanzo, F.; Sprik, M.; Sulpizi, M.; Blumberger, J. Absolute pKa Values and Solvation Structure of Amino Acids from Density Functional Based Molecular Dynamics Simulation. J. Chem. Theory Comput. 2011, 7, 1951-1961. (7) Chen, Y.-L.; Doltsinis, N. L.; Hider, R.C.; Barlow, D.J. Prediction of Absolute Hydroxyl pKa Values for 3-Hydroxypyridin-4-ones. J. Phys. Chem. Lett. 2012, 3, 2980-2985.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 46

(8) Liptak, M. D.; Shields, G. C. Accurate pKa Calculations for Carboxylic Acids Using Complete Basis Set and Gaussian-n Models Combined with CPCM Continuum Solvation Methods. J. Am. Chem. Soc. 2001, 123, 7314-7319. (9) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. Absolute pKa Determinations for Substituted Phenols. J. Am. Chem. Soc. 2002, 124, 6421-6427. (10)

Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion−

Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066-16081. (11)

Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Single-Ion Solvation Free Energies and the Normal

Hydrogen Electrode Potential in Methanol, Acetonitrile, and Dimethyl Sulfoxide. J. Phys. Chem. B 2007, 111, 408-422. (12)

Verdolino, V.; Cammi, R.; Munk, B. H.; Schlegel, H. B. Calculation of pKa Values of

Nucleobases and the Guanine Oxidation Products Guanidinohydantoin and Spiroiminodihydantoin using Density Functional Theory and a Polarizable Continuum Model. J. Phys. Chem. B 2008, 112, 16860-16873. (13)

Ho, J.; Coote, M. L. pKa Calculation of Some Biologically Important Carbon Acids - An

Assessment of Contemporary Theoretical Procedures. J. Chem. Theory Comput. 2009, 5, 295-306. (14)

Ho, J.; Coote, M. L. A Universal Approach for Continuum Solvent pKa Calculations: Are We

There Yet? Theor. Chem. Acc. 2010, 125, 3-21. (15)

Ho, J. Are Thermodynamic Cycles Necessary for Continuum Solvent Calculation of pKas and

Reduction Potentials? Phys. Chem. Chem. Phys. 2015, 17, 2859-2868.


Page 35 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(16)

Ho, J. Predicting pKa in Implicit Solvents: Current Status and Future Directions. Aust. J. Chem.

2014, 67, 1441−1460. (17)

Jensen, J. H. Predicting accurate absolute binding energies in aqueous solution: thermodynamic

considerations for electronic structure methods. Phys. Chem. Chem. Phys. 2015, 17, 12441-12451. (18)

Pliego, J. R., Jr.; Riveros, J. M. The Cluster−Continuum Model for the Calculation of the

Solvation Free Energy of Ionic Species. J. Phys. Chem. A 2001, 105, 7241-7247. (19)

Pliego, J. R., Jr.; Riveros, J. M. Theoretical Calculation of pKa Using the Cluster−Continuum

Model. J. Phys. Chem. A 2002, 106, 7434-7439. (20)

Takano, Y.; Houk, K. N. Benchmarking the Conductor-like Polarizable Continuum Model

(CPCM) for Aqueous Solvation Free Energies of Neutral and Ionic Organic Molecules. J. Chem. Theory Comput. 2005, 1, 70-77. (21)

Matsui, T.; Kitagawa, Y.; Okumura, M.; Shigeta, Y. Sakaki, S. Consistent Scheme for Computing

Standard Hydrogen Electrode and Redox Potentials. J. Comp. Chem. 2013, 34, 21-26. (22)

Lee, T.-B.; McKee, M. L. Dependence of pKa on Solute Cavity for Diprotic and Triprotic Acids.

Phys. Chem. Chem. Phys. 2011, 13, 10258-10269. (23)

Matsui, T.; Oshiyama, A.; Shigeta, Y. A Simple Scheme for Estimating the pKa Values of 5-

substituted Uracils. Chem. Phys. Lett. 2011, 502, 248-252. (24)

Matsui, T.; Baba, T.; Kamiya, K.; Shigeta, Y. An Accurate Density Functional Theory Based

Estimation of pKa Values of Polar Residues Combined with Experimental Data: from Amino Acids to Minimal Proteins. Phys. Chem. Chem. Phys. 2012, 14, 4181-4187.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(25)

Page 36 of 46

Zhang, S.; Baker, J.; Pulay, P. A Reliable and Efficient First Principles-Based Method for

Predicting pKa Values. 1. Methodology. J. Phys. Chem. A 2010, 114, 425-431. (26)

Zhang, S.; Baker, J.; Pulay, P. A Reliable and Efficient First Principles-Based Method for

Predicting pKa Values. 2. Organic Acids. J. Phys. Chem. A 2010, 114, 432-442. (27)

Miguel, E. L. M.; Silva, P. L. Pliego, J. R. Theoretical Prediction of pKa in Methanol: Testing

SM8 and SMD Models for Carboxylic Acids, Phenols, and Amines. J. Phys.Chem. B 2014, 118. 5730-5739. (28)

Galano, A.; Pérez-González, A.; Castañeda-Arriaga, R.; Muñoz-Rugeles, L.; Mendoza-

Sarmiento, G.; Romero-˜Silva, A.; Ibarra-Escutia, A.; Rebollar-Zepeda, A. M.; Léon-Carmona, J. R.; Hernández-Olivares, M. A.; Alvarez-Idaboy, J. R. Empirically Fitted Parameters for Calculating pKa Values with Small Deviations from Experiments Using a Simple Computational Strategy. J. Chem. Inf. Model. 2016, 56, 1714-1724. (29)

Derbel, N.; Clarot, I.; Mourer, M.; Regnouf-de-Vains J.-B.; Manuel F. Ruiz-López, M. F.

Intramolecular Interactions versus Hydration Effects on p-Guanidinoethyl-phenol Structure and pKa Values. J. Phys. Chem. A 2012, 116, 9404-9411. (30)

Casasnovas, R.; Ortega-Castro, J.; Frau, J.; Donoso, J.; Muñoz, F. Theoretical pKa Calculations

with Continuum Model Solvents, Alternative Protocols to Thermodynamic Cycles. Int. J. Quantum Chem. 2014, 114, 1350-1363. (31)

Rebollar-Zepeda, A. M.; Galano, A. Quantum Mechanical Based Approaches for Predicting pKa

Values of Carboxylic Acids: Evaluating the Performance of Different Strategies. RSC Adv. 2016, 6, 112057-112064.


Page 37 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(32)

Rossini E.; Netz, R. R.; Knapp, E.-W. Computing pKa Values in Different Solvents by

Electrostatic Transformation. J. Chem. Theory Compt. 2016, 12, 3360-3367. (33)

Leon-Carmona, J. R.; Galano, A.; Alvarez-Idaboy, J. R. Deprotonation Routes of Anthocyanidins

in Aqueous Solution, pKa Values, and Speciation under Physiological Conditions. RSC Adv. 2016, 6, 53421-53429. (34)

Maekawa, S.; Matsui, T.; Hirao, K.; Shigeta, Y. Theoretical Study on Reaction Mechanisms of

Nitrite Reduction by Copper Nitrite Complexes: Toward Understanding and Controlling Possible Mechanisms of Copper Nitrite Reductase. J. Phys. Chem. B 2015, 119, 5392-5403. (35)

Matsui, T.; Kitagawa, Y.; Okumura, M.; Shigeta, Y. Accurate Standard Hydrogen Electrode

Potential and Applications to the Redox Potentials of Vitamin C and NAD/NADH. J. Phys. Chem. A 2015, 119, 369-376. (36)

Marenich, A. V.; Ho, J.; Coote, M. L.; Cramer, C. J.; Truhlar, D. G. Computational

electrochemistry: prediction of liquid-phase reduction potentials. Phys. Chem. Chem. Phys. 2014, 16, 15068-15106. (37)

Thapa, B.; Schlegel, H. B. Theoretical Calculation of pKa’s of Selenols in Aqueous Solution

Using an Implicit Solvation Model and Explicit Water Molecules. J. Phys. Chem. A 2016, 120, 8916-8922. (38)

Ishikawa, A.; Nakai, H. Quantum chemical approach for condensed-phase thermochemistry (III):

Accurate evaluation of proton hydration energy and standard hydrogen electrode potential. Chem. Phys. Lett. 2016, 650, 159-164. (39)

Hall, H. K. Jr. Correlation of the Base Strengths of Amines. J. Am. Chem. Soc. 1957, 79, 5441-

5444. 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(40)

Page 38 of 46

Brown, H. C. et al., in Braude, E. A. and Nachod, F. C. Determination of Organic Structures by

Physical Methods, Academic Press, New York, 1955. (41)

Franzen, V. Beziehungen Zwischen Konstitution und Katalytischer Aktivität von Thiolaminen

Bei der Katalyse der Intramolekularen Cannizzaro-Reaktion. Chem. Ber. 1957, 90, 623-633. (42)

Dippy, J. F. J.; Hughes, S. R. C.; Rozanski, A. The Dissociation Constants of Some

Symmetrically Disubstituted Succinic Acids. J. Chem. Soc. 1959, 2492-2498. (43)

Edsall, J. T.; Wyman J. Biophysical Chemistry, Academic Press, Inc., New York, 1958.

(44)

Bruice, T. C.; Schmir, G. L. Imidazole Catalysis. II. The Reaction of Substituted Imidazoles with

Phenyl Acetates in Aqueous Solution. J. Am. Chem. Soc. 1958, 80, 148-156. (45)

Gawron, O.; Duggan, M.; Grelechi, C. J. Manometric Determination of Dissociation Constants

of Phenols. Anal. Chem. 1952, 24, 969-970. (46)

Fickling, M. M.; Fischer, A.; Mann, B. R.; Packer, J.; Vaughan, J. Hammett Substituent

Constants for Electron-withdrawing Substituents: Dissociation of Phenols, Anilinium Ions and Dimethylanilinium Ions. J. Am. Chem. Soc. 1959, 81, 4226-4230. (47)

Mukherjee, L. M.; Grunwald, E. Physical Properties and Hydrogen-bonding in the System

Ethanol–2,2,2-Trifluoro-ethanol. J. Phys. Chem. 1958, 62, 1311-1314. (48)

Ballinger, P.; Long, F. A. Acid Ionization Constants of Alcohols. I. Trifluoroethanol in the

Solvents H2O and D2O1. J. Am. Chem. Soc. 1959, 81, 1050-1053. (49)

Ballinger, P.; Long, F. A. Acid Ionization Constants of Alcohols. II. Acidities of Some

Substituted Methanols and Related Compounds. J. Am. Chem. Soc. 1960, 82, 795-798.


Page 39 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50)

Kreevoy, M. M.; Harper, E. T.; Duvall, R. E. Wilgus, H. S. III; Ditsch, L. T. Inductive Effects on

the Acid Dissociation Constants of Mercaptans. J. Am. Chem. Soc. 1960, 82, 4899-4902. (51)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.;

Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (52)

Becke, A.D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem.

Phys. 1993, 98, 5648-5652. (53)

Becke, A.D. A New Mixing of Hartree–Fock and Local Density-functional Theories. J. Chem.

Phys. 1993, 98, 1372-1377. (54)

Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A Long-Range Correction Scheme for Generalized-

Gradient-Approximation Exchange Functionals. J. Chem. Phys. 2001, 115, 3540-3544. (55)

Yanai, T.; Tew, D.; Handy, N. A New Hybrid Exchange-correlation Functional Using the

Coulomb-attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51-57. (56)

Grimme, S. Semiempirical Hybrid Density Functional with Perturbative Second-order

Correlation. J. Chem. Phys. 2006, 124, 034108. 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(57)

Page 40 of 46

Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable

Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158-6170. (58)

Vydrov, O. A.; Scuseria, G. E. Assessment of a Long-range Corrected Hybrid Functional. J.

Chem. Phys. 2006, 125, 234109. (59)

Vydrov, O. A.; Heyd, J.; Krukau, A.; Scuseria, G. E. Importance of Short-range versus Long-

range Hartree-Fock Exchange for the Performance of Hybrid Density Functionals. J. Chem. Phys. 2006, 125, 074106. (60)

Becke, A. D. Density-functional Thermochemistry. V. Systematic Optimization of Exchange-

correlation Functionals. J. Chem. Phys. 1997, 107, 8554-8560. (61)

Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-range

Dispersion Correction. J. Comp. Chem. 2006, 27, 1787-1799. (62)

Chai, J.-D.; Head-Gordon, M. Long-range Corrected Hybrid Density Functionals with Damped

Atom–atom dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615-6620. (63)

Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-group Thermochemistry,

Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125, 194101. (64)

Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group

Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215-241. (65)

Peverati, R.; Truhlar, D. G. Screened-exchange Density Functionals with Broad Accuracy for

Chemistry and Solid-state Physics. Phys. Chem. Chem. Phys. 2012, 14, 16187-16191. 40 ACS Paragon Plus Environment

Page 41 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(66)

Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional

Ladder: Nonempirical Meta-generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. (67)

Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative Assessment of a New

Nonempirical Density Functional: Molecules and Hydrogen-bonded Complexes. J. Chem. Phys. 2003, 119, 12129-12137. (68)

Fuentealba, P.; Preuss, H.; Stoll, H.; Szentpály, L.v. A Proper Account of Core-polarization with

Pseudopotentials: Single Valence-electron Alkali Compounds. Chem. Phys. Lett. 1982, 89, 418-422. (69)

Hay, P. J.; Wadt, W. R. Ab initio Effective Core Potentials for Molecular Calculations. Potentials

for K to Au Including the Outermost Core Orbitals. J. Chem. Phys. 1985, 82, 299-310. (70)

Wadt, W. R.; Hay, P. J. Ab initio Effective Core Potentials for Molecular Calculations. Potentials

for Main Group Elements Na to Bi. J. Chem. Phys. 1985, 82, 284-298. (71)

Ochterski, J. W.; Petersson, G. A.; Montgomery, J. A. Jr. A Complete Basis Set Model

Chemistry. V. Extensions to Six or More Heavy Atoms. J. Chem. Phys. 1996, 104, 2598-2619. (72)

Montgomery, J. A. Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. A Complete Basis Set

Model Chemistry. VI. Use of Density Functional Geometries and Frequencies. J. Chem. Phys. 1999, 110, 2822-2827. (73)

Montgomery, J. A. Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. A Complete Basis Set

Model Chemistry. VII. Use of the Minimum Population Localization Method. J. Chem. Phys. 2000, 112, 6532-6542.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(74)

Page 42 of 46

Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. Gaussian-3 (G3)

Theory for Molecules Containing First and Second-Row Atoms. J. Chem. Phys. 1998, 109, 77647776. (75)

Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-3 Theory using Density

Functional Geometries and Zero-Point Energies. J. Chem. Phys. 1999, 110, 7650-7657. (76)

Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Rassolov, V.; Pople, J. A. A. Gaussian-3 Theory

Using Reduced Moller-Plesset Orders. J. Chem. Phys. 1999, 110, 4703-4709. (77)

Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 theory. J. Chem. Phys. 2007, 126,

084108. (78)

Klamt, A.; Schüürmann, G. COSMO: a New Approach to Dielectric Screening in Solvents with

Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc., Perkin Trans 2, 1993, 799-805. (79)

Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in

Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995-2001. (80)

Cancès, M. T.; Mennucci, B.; Tomasi, J. A New Integral Equation Formalism for the Polarizable

Continuum Model: Theoretical Background and Applications to Isotropic and Anisotropic Dielectrics. J. Chem. Phys. 1997, 107, 3032-3041. (81)

Mennucci, B.; Tomasi, J. Continuum Solvation Models: A New Approach to the Problem of

Solute’s Charge Distribution and Cavity Boundaries. J. Chem. Phys. 1997, 106, 5151-5158. (82)

Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute

Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378-6396. 42 ACS Paragon Plus Environment

Page 43 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(83)

Johnson, E. R.; Becke, A. D.; Sherrill, C. D.; DiLabio, G. A. Oscillations in Meta-generalized-

gradient Approximation Potential Energy Surfaces for Dispersion-bound Complexes. J. Chem. Phys. 2009, 131, 034111. (84)

Wheeler, S. E.; Houk, K. N. Integration Grid Errors for Meta-GGA-Predicted Reaction Energies:

Origin of Grid Errors for the M06 Suite of Functionals. J. Chem. Theory Comput. 2010, 6, 395-404. (85)

Sato, T.; Tsuneda, T.; Hirao, K. A Density-functional Study on π-aromatic Interaction: Benzene

Dimer and Naphthalene Dimer. J. Chem. Phys. 2005, 123, 104307. (86)

Matsui, T.; Sato, T.; Shigeta, Y.; Hirao, K. Sequence-dependent Proton-transfer Reaction in

Stacked GC pair II: The Origin of Stabilities of Proton-transfer Products. Chem. Phys. Lett. 2009, 478, 238-242. (87)

Klamt, A.; Mennucci, B.; Tomasi, J.; Barone, V.; Curutchet, C.; Orozco, M.; Luque. F. J. On the

Performance of Continuum Solvation Methods. A Comment on “Universal Approaches to Solvation Modeling”. Acc. Chem. Res. 2009, 42, 289-292. (88)


Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 46

A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision E. 01; Gaussian, Inc.: Wallingford, CT, 2009. (89)


Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16, revision A. 03; Gaussian, Inc.: Wallingford, CT, 2016.


Page 45 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents 258x254mm (200 x 200 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Contents (Small file) 129x127mm (72 x 72 DPI)


Page 46 of 46

Assessment of Methodology and Chemical Group Dependences in the

Recommend Documents