Density Functional Theory and Isodesmic Reaction Based Prediction of Four Stepwise Protonation Constants, as log K H (n) , for Nitrilotriacetic Acid. The Importance of a Kind and Protonated Form of a Reference Molecule Used Krishna K. Govender and Ignacy Cukrowski* Department of Chemistry, Faculty of Natural Sciences, UniVersity of Pretoria, Lynnwood Road, Hillcrest, Pretoria 0002, South Africa ReceiVed: July 29, 2009; ReVised Manuscript ReceiVed: NoVember 25, 2009 An explicit application of isodesmic reaction (a proton exchange between the studied and structurally similar reference molecule), where the free energy change of the protonation reaction in water was obtained using the free energies in solution from a single continuum model, was used to predict stepwise protonation constants of nitrilotriacetic acid. Calculations were performed at the RB3LYP/6-311+G(d,p) level of theory in conjunction with the PCM-UA0 solvation model. Five reference molecules were investigated. It has been established that one must pay special attention to structural similarities between the studied and reference molecules and selection of a protonated form of the reference molecule. The protonation reactions in which the studied and reference molecule are involved in must be (if possible) of the same order; e.g., the first (or generally nth) protonation reaction of the reference molecule must be used to compute the first (or nth) protonation constant of the studied molecule. The lowest energy conformer must always be used. The first, second, third, and fourth computed protonation constants differed, on average, from experimental values by 3.3, 0.8, 0.2, and 0.2 log units, respectively. It appears that the charge on the reference molecule has more decisive influence on the accuracy of computed protonation constants than its structural differences when compared with the studied molecule. Results reported can be used as a guide in constructing isodesmic reactions useful for the theoretical prediction of protonation constants by use of methodology described in this work. 1. Introduction Knowledge of protonation, K H , and dissociation, K a , constants is of special interest to many chemists and life scientists 1 as they constitute important thermodynamic property of a com- pound that might be of biological, medicinal, or industrial (just to mention few) importance. Although a number of experimental techniques has been developed to measure protonation/dissocia- tion constants under various experimental conditions, many of the chemical species are not easily amenable to a full experi- mental characterization. 2 A number of papers has been reported 1-54 on theoretical prediction of dissociation constants. Most of them employed thermodynamic cycles (TC) to compute the free energies of dissociation reaction. Often, high-level theories were used in the gas-phase calculations where they are known to be accurate. The solution-phase calculations were used to provide the solvation energies (G sol ); usually low-level continuum models were employed for the purpose. When the above protocol is used, the absolute pK a value is obtained. To avoid uncertainties related to the solvation energies of either H + or H 3 O + ions, isodesmic reaction (IRn) was incorporated within TC; 1,3-7 this protocol of calculation results in relative 3 pK a values. Results reported to date predominantly describe the calculations of singly charged molecules, either anions 1-29 (a study of doubly charged anions is very rare) or cations. 30-35 This is most likely due to the fact that (i) it is very difficult for DFT methods to properly describe anions (with multiple negative charges) in the gas phase because in the absence of an external stabilization of the charge (e.g., solvent) DFT methods have a bias toward “over-delocalization” of the charge (one might observe bonds that are longer than expected and significant reduction on the HOMO-LUMO gap) 55 and (ii) computational evaluations of ionic solvation free energies for highly charged anions are inaccurate (these energies are highly dependent on the solvation model used due to different models chosen to generate the “best” electrostatic cavity). 56 Accuracies achieved thus far for computed dissociation constants (for a singly dissociable organic acids) are often within (1.0 log unit, on average, when compared with experimentally available values, but differences of several log units are not uncommon. 3-5,56 Recently we reported the DFT-predicted four stepwise protonation constants, expressed as log K H (n) , for a highly charged molecule nitrilotripropanoic acid (NTPA). 54 An explicit application of an isodesmic reaction involving two structurally similar ligands, where the free energy change of the protonation reaction in solution was obtained using the free energies in solution from a single continuum model, resulted in the average difference between predicted and experimental stepwise pro- tonation constants being (0.5 log unit. This suggested that, in principle, even though serious concerns were expressed, 55,56 accurate determination of stepwise protonation constants for highly negatively charged molecules is possible. In this paper our focus is on parameters that influence the accuracy in predicting four consecutive protonation constants when the IRn-based procedure, as reported by us recently, 54 is employed. It is important to investigate a wide range of polycharged compounds (with negative and positive charges) to establish (i) whether implementation of the protocol can indeed produce consistently good predictions, (ii) how significant the selection of a reference molecule is from the point of view of its structural similarity to the studied compound, (iii) to what * Corresponding author. E-mail: ignacy.cukrowski@up.ac.za. J. Phys. Chem. A 2010, 114, 1868–1878 1868 10.1021/jp9092964  2010 American Chemical Society Published on Web 01/11/2010