Baixe Study Material on Molecular Structures and Physics from Various Sources e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity! Density Functional Theory and Isodesmic Reaction Based Prediction of Four Stepwise Protonation Constants, as log KH(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, KH, and dissociation, Ka, constants is of special interest to many chemists and life scientists1 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 reported1-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 (Gsol); usually low-level continuum models were employed for the purpose. When the above protocol is used, the absolute pKa value is obtained. To avoid uncertainties related to the solvation energies of either H+ or H3O+ ions, isodesmic reaction (IRn) was incorporated within TC;1,3-7 this protocol of calculation results in relative3 pKa values. Results reported to date predominantly describe the calculations of singly charged molecules, either anions1-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 KH(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 XXXX, xxx, 000 A 10.1021/jp9092964 XXXX American Chemical Society degree the selection of a conformer impacts the accuracy in computed protonation constants, (iv) if the selection of a different protonated form of the reference molecule has an influence on the accuracy of computed protonation constants, (v) which one, structural similarity or the charge on a protonated form of a reference molecule, plays a more important role when the IRn-based methodology described in this work is employed, (vi) if it is possible to predict protonation constants in correct order, as determined from an experiment, (vii) to what extent the selection of level of theory and a basis set has an influence on the accuracy in computed protonation constants, and (viii) how significant is the kind of solvation model employed. These are just a few of the important questions, and this work will address only some of them. We have chosen the ligand nitrilotriacetic acid (NTA) because it is an important derivative of glycine that is widely studied due to its excellent chelating abilities.57 This is a ligand that has enjoyed numerous applications in medicine,58,59 biochemistry,58,59 and industry.60-62 In medicinal and biological studies it was shown that aliphatic amine salts of NTA inhibit the growth of bacteria and fungi and have herbicidal activity.63 NTA has also been used as a transient phytoextraction agent that combines high biodegradability and low phytotoxicity with chelating strength.64 The application of IRn requires a reference molecule and its selection appears to be crucial. Molecules shown in Figure 1 (iminodiacetic acid (IDA), methyliminodiacetic acid (MIDA), ethyliminodiacetic acid (EIDA), propyliminodiacetic acid (PIDA), and hydroxyethyliminodiacetic acid (HIDA)) were used here as they have many structural similarities with the compound of interest (NTA) and their experimental stepwise protonation constants are well-known.57 2. Computational Details All calculations were performed using GAUSSIAN 03, revision D.01,65 on a 64-bit Linux workstation in parallel environment (Opensuse 10.3). Molecular visualizations were accomplished with the aid of GaussView 4.66 Since it is of paramount importance to include diffuse functions for anions,19 both gas-phase and solvent (water, ε ) 78.39) optimizations were performed at the RB3LYP level of theory67 in conjunction with a 6-311+G(d,p) basis set. Full solvent-optimization was performed with the default solvation model provided by Gaussian, i.e., Tomasi’s polarized continuum model (PCM),68-70 and UA0 radii (united atom topological model). We have chosen this model because it generated acceptable results in the Figure 1. Top view of the ligands (in fully protonated forms) discussed in this work. B J. Phys. Chem. A, Vol. xxx, No. xx, XXXX Govender and Cukrowski IRn employed here can be seen as a competition reaction between two ligands for a proton (proton transfer reaction), and for the first protonation constant of NTA it can be written as To investigate the impact the kind and protonated form of the reference molecule has on the theoretically generated protonation constants of NTA, a large number of isodesmic reactions was tested. Here, each reference molecule has three protonation constants (NTA has four), hence for each pair of ligands (L(1) and L(2)) 12 isodesmic reactions (such as eq 1, but involving different protonated forms of the ligands) had to be considered. For simplicity, only the first protonation reaction (PRn) in which each of the two ligands (NTA and a reference molecule) is involved is shown as eqs 2 and 3 A protonation reaction is the reverse of a weak acid dissociation reaction (DRn) and in the case of stepwise reactions the following holds where k ) 1 + m - n, m and n represent the highest dissociation constant (here m ) 4) and an nth consecutive dissociation constant (1 e n e m), respectively, and k applies to a kth consecutive protonation constant, 1 e k e m. Note that the ligand NTA has three acidic groups and only three dissociation constants would be reported when, for example, the TC-based methodology was employed. However, due to the protonation/ deprotonation of the N-atom in NTA, it is of paramount importance to consider also the first protonation constant, log KH(1). From this it follows that the fourth dissociation constant of NTA is linked through eq 4 with the first protonation constant of this ligand. In the case of the reference molecule, the third dissociation constant is linked through eq 4 with the first protonation constant of this reference molecule. (The above is provided here for convenience and to ensure clarity in nomen- clature used in this work.) The relationships between the change in Gibbs energies for protonation and dissociation reactions applicable to eqs 2 and 3 can be written as where ∆GPRn(aq)L(1)(k) and ∆GDRn(aq)L(1)(n) refer to the kth (here k ) 1) stepwise protonation reaction and relevant nth (here n ) 4 for NTA) stepwise dissociation reaction, respectively (the same applies to the reference molecule L(2), but n ) 3). The change in Gibbs energies for each PRn, eqs 2 and 3, can be written, respectively, as The isodesmic reaction of interest (eq 1) can be obtained by subtracting eq 3 from eq 2, and hence by subtracting eq 8 from eq 7 one obtains expressions for the change in Gibbs energy applicable to this isodesmic reaction where the uncertainty related to Gaq(H+) is no longer applicable as this term cancels out (this eliminates any error that might have been introduced by the use of an experimental value for this quantity). Equation 10 was used to calculate ∆GIRn(aq) of IRn (eq 1) from appropriate Gibbs energies obtained for relevant and fully solvent-optimized structures of the ligand NTA and reference molecule L(2). Table S8 in Supporting Information provides the ZPVE-corrected minimum energies, Emin, as well as the Gibbs free energies of NTA and all of the reference molecules studied here. The value for ∆GPRn(aq)L(2) was obtained from the well- known relationship using the reported protonation constants57 (at 20 and 25 °C, µ ) 0.0 and 0.1 M) of the reference molecule L(2). Once ∆GIRn(aq) and ∆GPRn(aq)L(2) have been calculated, the value of ∆GPRn(aq)L(1), which is needed to calculate the protonation constants of NTA from eq 11, was obtained from eq 9. Table 1 provides values for the functions required to calculate protonation constants, calculated and experimental protonation constants of NTA, along with differences between calculated and experimental protonation constants (δ). Values for ∆GPRn(aq)L(2)(k) were calculated from the experimentally available stepwise protonation constants,57 which have been reproduced in Table S9, Supporting Information. Only isodesmic reactions that produced the best results are shown in Table 1; the remaining results are provided in Table S10, Supporting Information. There are several interesting observations one can make from the analysis of data seen in Table 1. (i) All isodesmic reactions seen in Table 1 predicted four protonation constants in correct order, log KH(1) > log KH(2) > log KH(3) > log KH(4), as observed from an experiment. One must realize that the experimental values of the second, third, and fourth protonation constants of NTA differ between each other only by one log unit (or less) and this is a typical error reported in theoretically predicted values of dissociation constants reported to date for organic acids containing only one carboxylic group. Clearly, results obtained here should be seen as satisfac- tory and encouraging for further studies aimed at improvements in accuracy of theoretical predictions. (ii) As one goes from IDA to HIDA, it is seen that the prediction of protonation constants becomes more accurate, with HIDA yielding the most accurate estimates. The observed trend L(1)(aq) + HL(2)(aq) ) HL(1)(aq) + L(2)(aq) ∆GIRn(aq) (1) H+ + L(1) T HL(1) ∆GPRn(aq)L(1) (2) H+ + L(2) T HL(2) ∆GPRn(aq)L(2) (3) log KH (k) ) log 1 Ka (n) ) pKa (n) (4) ∆GPRn(aq)L(1) (1) ) -∆GDRn(aq)L(1) (4) (5) ∆GPRn(aq)L(2) (1) ) -∆GDRn(aq)L(2) (3) (6) ∆GPRn(aq)L(1) ) Gaq(HL(1)) - Gaq(H +) - Gaq(L(1)) (7) ∆GPRn(aq)L(2) ) Gaq(HL(2)) - Gaq(H +) - Gaq(L(2)) (8) ∆GIRn(aq) ) ∆GPRn(aq)L(1) - ∆GPRn(aq)L(2) (9) ∆GIRn(aq) ) Gaq(HL(1)) - Gaq(L(1)) - Gaq(HL(2)) + Gaq(L(2)) (10) ∆G(aq) ) -RT ln K (11) Protonation Constants for Nitrilotriacetic Acid J. Phys. Chem. A, Vol. xxx, No. xx, XXXX E can be linked with an increase in structural similarity between the reference molecule and NTA. (iii) Regardless of the reference molecule studied here, the use of HL- and L2- of L(2) resulted in the best estimates of the first protonation constant of NTA that involves L3- and HL2- forms of NTA. This can be generalized: computing the first protonation constant of studied molecule requires the first protonation constant and components involved in the first protonation reaction of a reference molecule. (iv) The best estimates in the second and third theoretical protonation constants of NTA were always obtained (regardless of the reference molecule involved in IRn) when the second protonation constant of the reference molecule was used. (v) The best prediction of the fourth (highest) protonation constant of NTA always required the highest (third) protonation constant of the reference molecule. Here again one might generalize that, to calculate the highest protonation constant of the studied molecule, one must involve either the same or as high as possible protonation constant of the reference molecule. (vi) The least accurate computed log KH, regardless of the reference molecule used, was always obtained for the first protonation reaction of NTA that involved the most negatively charged forms, that of the studied ligand NTA, L(1)3-, and the reference molecule L(2)2-. (vii) Interestingly, for all reference molecules studied, the smallest errors in the predicted log KH values were obtained for the second and third protonation constants of NTA and all of them might be regarded as acceptable estimates at µ ) 0.0 M. One might rationalize this observation in terms of a charge placed on the studied and reference molecules. In these isodesmic reactions, the HL2-, H2L-, and H3L forms of NTA as well as HL- and H2L forms of the reference molecules were used where the charge varied from -2 to 0. Since good predictions in log KH(2) and log KH(3) of NTA were obtained, regardless of the reference molecule studied here, one might conclude that (a) the more similar the charges on the studied and reference molecule are, the better prediction is achieved, even though significant differences in structures of the reference molecules are present, and (b) the charge on a reference molecule has more significant impact than its structural similarity to NTA, when accuracy in predicted protonation constants is considered. (viii) Unexpectedly, somewhat worse estimates in log KH(4) values were obtained when molecules involved in the isodesmic reactions had 0 and +1 charge; this applies to all systems studied here. This might be a significant observation because many important ligands (among them macrocyclic ligands) do not have carboxylic groups (they are neutral in their fully deprotonated form) and when protonated they have multiple and positive charges on them. Clearly, this requires a dedicated investigation to establish whether the IRn-based protocol described here can be applied successfully for positively charged molecules. If one considers the structural features of each of the reference molecules and that of NTA, it is possible to conclude that the molecule that is structurally most similar to NTA is HIDA. Both NTA and HIDA have one N-donor atom (R3-N:) and three O-donor atoms, with one of them being part of the -OH group (in HIDA) instead of the -COOH group (in NTA). All the other reference molecules have only two O-atoms. The structural similarity of HIDA to NTA correlates well with results seen in Table 1. It also appears that the cavity of the reference molecule, when full energy optimization is performed in solvent, plays a significant role. The values of δ for the first protonation constant of NTA (at µ ) 0.0 M) were 4.81, 3.59, 3.05, and 3.08 log unit when IDA, MIDA, EIDA, and PIDA were used as the reference molecule, respectively. The same trend, the decrease in error with an increase in the cavity of the reference molecule, is seen for the fourth protonation constant, namely -2.19, -1.63, -0.99, and -0.70 log unit, respectively, for the same reference TABLE 1: Comparison of Experimental57 (Exp) and Calculated (from Isodesmic Reactions) Stepwise Protonation Constants of NTA (L(1)), as log KH, Using Protonation Constants of the Reference Molecules at Ionic Strength µ ) 0.0 or 0.1 M and 25 °Ca reaction ∆GIRn(aq) ∆GPRn(aq) for L(1) log KH expb δ ∆GPRn(aq) for L(1) log KH expc δ L(2) ) IDA L(1)3- + HL(2)- ) HL(1)2- + L(2)2- -7.298 -20.654 15.14 10.334 4.81 -20.040 14.69 9.66 5.03 HL(1)2- + H2L(2) ) H2L(1)- + HL(2)- 0.052 -3.822 2.80 2.94 -0.14 -3.522 2.58 2.52 0.06 H2L(1)- + H2L(2) ) H3L(1) + HL(2)- 2.504 -1.371 1.00 2.00d -1.00 -1.071 0.78 1.81 -1.03 H3L(1) + H3L(2)+ ) H4L(1)+ + H2L(2) 4.143 1.620 -1.19 1.00c -2.19 1.729 -1.27 1.00 -2.27 L(2) ) MIDA L(1)3- + HL(2)- ) HL(1)2- + L(2)2- -5.342 -18.998 13.93 10.334 3.59 -18.425 13.51 9.66 3.85 HL(1)2- + H2L(2) ) H2L(1)- + HL(2)- -0.171 -3.705 2.72 2.94 -0.22 -3.336 2.45 2.52 -0.07 H2L(1)- + H2L(2) ) H3L(1) + HL(2)- 2.280 -1.253 0.92 2.00d -1.08 -0.885 0.65 1.81 -1.16 H3L(1) + H3L(2)+ ) H4L(1)+ + H2L(2) 3.458 0.866 -0.63 1.00c -1.63 0.866 -0.63 1.00 -1.63 L(2) ) EIDA L(1)3- + HL(2)- ) HL(1)2- + L(2)2- -4.435 -18.255 13.38 10.334 3.05 -18.010 13.20 9.66 3.54 HL(1)2- + H2L(2) ) H2L(1)- + HL(2)- -1.605 -5.288 3.88 2.94 0.94 -4.633 3.40 2.52 0.88 H2L(1)- + H2L(2) ) H3L(1) + HL(2)- 0.847 -2.836 2.08 2.00d 0.08 -2.182 1.60 1.81 -0.21 H3L(1) + H3L(2)+ ) H4L(1)+ + H2L(2) 2.175 -0.008 0.01 1.00c -0.99 -0.008 0.01 1.00 -0.99 L(2) ) PIDA L(1)3- + HL(2)- ) HL(1)2- + L(2)2- -4.074 -18.303 13.42 10.334 3.08 -17.785 13.04 9.66 3.38 HL(1)2- + H2L(2) ) H2L(1)- + HL(2)- -2.081 -5.478 4.02 2.94 1.08 -5.137 3.77 2.52 1.25 H2L(1)- + H2L(2) ) H3L(1) + HL(2)- 0.370 -3.027 2.22 2.00d 0.22 -2.686 1.97 1.81 0.16 H3L(1) + H3L(2)+ ) H4L(1)+ + H2L(2) 1.089 -0.412 0.30 1.00c -0.70 -1.326 0.97 1.00 -0.03 L(2) ) HIDAc L(1)3- + HL(2)- ) HL(1)2- + L(2)2- -3.535 -15.377 11.27 10.334 0.94 -15.377 11.27 9.66 1.61 HL(1)2- + H2L(2) ) H2L(1)- + HL(2)- -1.398 -4.399 3.22 2.94 0.28 -4.399 3.22 2.52 0.70 H2L(1)- + H2L(2) ) H3L(1) + HL(2)- 1.054 -1.948 1.43 2.00d -0.57 -1.948 1.43 1.81 -0.38 H3L(1) + H3L(2)+ ) H4L(1)+ + H2L(2) 1.601 -0.581 0.43 1.00c -0.57 -0.581 0.43 1.00 -0.57 a All energies are reported in kcal mol-1. b µ ) 0.0 M and 20 °C. c µ ) 0.1 M and 25 °C. d µ ) 0.0 M and 25 °C. F J. Phys. Chem. A, Vol. xxx, No. xx, XXXX Govender and Cukrowski molecules. Also, careful attention needs to be paid to the positioning and presence of atoms, especially the heteroatom, as the additional -OH group present in HIDA, which is not present in IDA, MIDA, EIDA, or PIDA, seems to make a huge difference, as far as prediction of protonation constants is concerned. Additional Test. All four theoretically predicted protonation constants of NTA seen in Table 1 differ from experimental values by less than a single log unit when HIDA was used; results of this accuracy are often referred to as excellent in the literature. On the other hand, the absolute protonation constants obtained from TCs carried significantly larger errors (see Table S7, Supporting Information). One might argue that accuracy obtained from straight continuum model free energies in solution calculation is not directly comparable with that obtained from TC as totally different chemical reactions are involved. Since some solvent structures were not preserved in the gas phase, we performed a single point frequency calculation in the gas phase on the solvent-optimized structures. We selected IDA for the test, as inclusion of this reference molecule in the proton competition reaction (IRn) resulted in the worst results; hence one might expect significant improvement in protonation constants calculations. The test was performed for the reaction where L(1) and L(2) are NTA and the reference molecule IDA, respectively. Our aim here was to find out whether improvement of the Ggas component to a better (higher) level of theory improves the computed protonation constant when the TC-based protocol is implemented. The following protocol was imple- mented. Let Gsol stand for the Gibbs free energy in solution obtained from a single continuum model calculation at B3LYP/ 6-311+G(d,p) in conjunction with the PCM/UA0 solvation model. For each component involved in the above IRn (eq 12) one can compute the free energy of solvation from ∆Gsol ) Gsol - Ggas(1), where Ggas(1) is computed as a single point on the solution-phase geometries at the same B3LYP/6-311+G(d,p) level of theory. Let Ggas(2) stand for the improved estimate obtained at the higher level of theory, in this case also calculated as a single point on the solution-phase geometry, then Gsol(TC), the improved Gibbs free energy obtained from thermodynamic cycle calculations, can be expressed as Gsol(TC) ) Ggas(2) + ∆Gsol. The values of Ggas(2) obtained at the RMP2/6-311+G(d,p) level of theory along with other necessary components are presented in Table 2. From Gsol(TC) values we calculated the change in the Gibbs free energy for reaction 12, ∆Gsol(TC) ) 6.137 kcal mol-1, followed by the fourth protonation constant of NTA, log KH(4) ) -4.50. This is a much worse estimate when compared with the value of -1.19 obtained directly from implementation of IRn (see Table 1). An attempt was made to obtain Ggas(2) values using the G3/ 6-311+G(d,p) level of theory for the single point frequency calculations in the gas phase on the structures optimized at the B3LYP/6-311+G(d,p) level of theory in conjunction with the PCM/UA0 solvation model. Unfortunately, the solvent structures of H3L(1) and H2L(2) have not been preserved. The H-atom has moved away from nitrogen and protonated the -COO- group; see Figure S8, Supporting Information. Conformational Considerations. All protonated forms of NTA and the reference molecules were subjected to the Schrödinger’s Maestro85 conformational analysis in solvent. Generated, on the basis of molecular mechanics/molecular dynamics (MM/MD) principles, lowest energy conformers C-1 are shown in Figures S9-S14, Supporting Information. Table 3a provides energies EC-1 to EC-5 (in kJ/mol) for the five lowest energy conformers of all the protonated forms of NTA. MM/ MD-based SPC in solvent was also performed on the NTA DFT- structures shown in Figure 4; results are shown as ESPC in Table 3a. It is seen that all ESPC are larger than energies of the MM/ MD-generated conformers. Initially, all of the C-1 conformers of NTA were fully DFT-optimized in solvent; results obtained are shown in Table 3b. It was gratifying to note that the differences δG in Table 3b between the relevant structures became almost negligibly small. For all protonated forms of NTA, except H4L+, the value of δG is about (0.1 kcal mol-1. This is equivalent to about (0.07 log unit of the computed protonation constant, a typical experimental uncertainty. A similar procedure was applied to all the protonated forms of IDA, MIDA, EIDA, PIDA, and HIDA (Table S11, Supporting Information). It is seen in Table 3b that only two out of six of the DFT-optimized C-1 conformers have lower energies (HL2- and H2L-) when compared with the energies of the original structures seen in Figure 4. A similar observation applies to energies of the DFT-optimized C-1 conformers of the reference molecules. The lowest Gibbs free energies Gaq were used to calculate protonation constants; results are provided in Table S12, Supporting Information. Because the use of C-1 conformers has not resulted in any significant change in the computed protonation constants and since there was no direct correlation between energies of the DFT-optimized C-1 conformers and original structures, we have decided to use six lowest energy MM/MD-generated conformers of NTA and HIDA for further studies. HIDA was selected because this reference molecule, when used in IRn’s seen in Table 1, generated best estimates in protonation constants. Energy-minimized in solvent conformers C-1 to C-6 of all the protonated forms of NTA and HIDA are shown in Figures S15-S23, Supporting Information. Their ZPVE-corrected ener- gies (Emin) and Gibbs free energies (Gaq) along with the values of the self-constructed (S-c) structures are shown in Table 4; for each protonated form of L(1) and L(2), the lowest Emin and Gaq are printed in bold. Implementation of an extended con- formational analysis resulted in three new lower in energy conformers of NTA (H2L-, H3L, and H4L+) as well as HIDA (L2-, H2L, and H3L+). TABLE 2: Thermochemical Data Used To Calculate the Fourth Protonation Constant of NTA Involving IRn in Eq 12a molecule Gsol, au Ggas(1), au ∆Gsol, au Ggas(2), au Gsol(TC), au H3L(1) -740.332226 -740.252149 -0.080077 -738.373664 -738.453741 H3L(2)+ -512.862287 -512.727315 -0.134972 -511.414007 -511.548979 H4L(1)+ -740.764545 -740.625849 -0.138696 -738.746987 -738.885683 H2L(2) -512.423365 -512.349858 -0.073507 -511.03375 -511.107257 a Ggas(1) values were obtained at the B3LYP/6-311+G(d,p) level of theory. Ggas(2) values were obtained at the RMP2/6-311+G(d,p) level of theory. Gsol values were obtained at the B3LYP/6-311+G(d,p) level of theory in solvent (PCM/UA0). L(1) ) NTA, L(2) ) IDA. For details see the text. H3L(1) + H3L(2) + + H4L(1) + + H2L(2) (12) Protonation Constants for Nitrilotriacetic Acid J. Phys. Chem. A, Vol. xxx, No. xx, XXXX G involves highly charged anions (-3 and -2 for the studied and the reference molecule, respectively), a significantly larger than experimental value (by 3 log units) was obtained. Clearly, more work has to be done to improve this value. Nevertheless, we are convinced that results reported here can be used as a guide in constructing isodesmic reactions that are useful for the theoretical prediction of protonation/dissociation constants, particularly for compounds for which solvent and gas structures differ significantly. Acknowledgment. Financial support of the National Research Foundation of South Africa and the University of Pretoria is highly appreciated. Supporting Information Available: Reported crystal struc- tures of IDA, MIDA, and HIDA; structural matrix of the solvent- optimized H3L* and H3L forms of NTA, H2L form of IDA and MIDA, and H3L form of HIDA; energy-minimized structures in solvent (PCM/UA0) at the RB3LYP/6-311+G(d,p) level of theory for all self-constructed protonated forms of IDA, MIDA, EIDA, PIDA, and HIDA; results obtained from computations involving thermodynamic cycles; ZPVE-corrected minimum and Gibbs free energies of all protonated forms of NTA and reference molecules (IDA, MIDA, EIDA, PIDA, and HIDA) obtained from full energy optimization at the RB3LYP/6- 311+G(d,p) level of theory in conjunction with the PCM/UA0 solvation model; experimental stepwise protonation constants of IDA, MIDA, EIDA, PIDA, and HIDA at µ ) 0.0 and 0.1 M and 25 °C; comparison of experimental and calculated stepwise protonation constants of NTA using protonation constants of the reference molecules IDA, MIDA, EIDA, PIDA, and HIDA; output structures (H3L form of NTA and H2L form of IDA) obtained from the single point frequency calculations at the G3/ 6-311+G(d,p) level of theory on the solvent optimized mol- ecules optimized at the RB3LYP/6-311+G(d,p) level of theory in conjunction with the PCM/UA0 solvation model; energy- minimized structures, at the RB3LYP/6-311+G(d,p) level of theory in conjunction with the PCM/UA0 solvation model, for all protonated forms of C-1 conformers of NTA, IDA, MIDA, EIDA, PIDA, and HIDA; minimum energies of MM/MD- generated conformers in solvent and energies obtained from MM-based SPC performed on the IDA, MIDA, EIDA, PIDA, and HIDA structures; DFT-calculated solvent-optimized energies of all self-constructed protonated forms of IDA, MIDA, EIDA, PIDA, and HIDA and lowest energy MM/MD-generated C-1 conformers; comparison of experimental and calculated stepwise protonation constants of NTA obtained from the lowest energy MM/MD conformers of IDA, MIDA, EIDA, PIDA, and HIDA; energy-minimized structures in solvent (PCM/UA0) at the RB3LYP/6-311+G(d,p) level of theory for MM/MD-generated conformers of the protonated forms of NTA and HIDA. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Namazian, M.; Kalantary-Fotooh, F.; Noorbala, M. R.; Searles, D. J.; Coote, M. C. J. Mol. Struct. (THEOCHEM) 2006, 758, 275. (2) Jacquemin, D.; Perpete, E. A.; Ciofini, I.; Adamo, C. J. Phys. Chem. A 2008, 112, 794. (3) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 5610. (4) Namazian, M.; Heidary, H. J. Mol. Struct. (THEOCHEM) 2003, 620, 257. (5) Namazian, M.; Halvani, S.; Noorbala, M. R. J. Mol. Struct. (THEOCHEM) 2004, 711, 13. (6) Namazian, M.; Halvani, S. J. Chem. Thermodyn. 2006, 38, 1495. (7) Namazian, M.; Zakery, M.; Noorbala, M. R.; Coote, M. L. Chem. Phys. Lett. 2008, 451, 163. (8) Charif, I. E.; Mekelleche, S. M.; Villemin, D.; Mora-Diez, N. J. Mol. Struct. (THEOCHEM) 2007, 818, 1. (9) Liptak, M. D.; Shields, G. C. Int. J. Quantum Chem. 2001, 85, 727. (10) Liptak, M. D.; Shields, G. C. J. Am. Chem. Soc. 2001, 123, 7314. (11) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. J. Am. Chem. Soc. 2002, 124, 6421. (12) Bensen, M. T.; Moser, M. L.; Peterman, D. R.; Dinescu, A. J. Mol. Struct. (THEOCHEM) 2008, 867, 71. (13) Sang-Aroon, W.; Ruangpornvisuti, V. Int. J. Quantum Chem. 2008, 108, 1181. (14) Jang, Y. H.; Goddard, W. A., III; Noyes, K. T.; Sowers, L. C.; Hwang, S.; Chung, D. S. J. Phys. Chem. B 2003, 107, 344. (15) Jitariu, L. C.; Masters, A. J.; Hiller, I. H. J. Chem. Phys. 2004, 121, 7795. (16) Bryantsev, V. S.; Diallo, M. S.; Goddard, W. A., III. J. Phys. Chem. A 2007, 111, 4422. (17) Gao, D.; Svoronos, P.; Wong, P. K.; Maddalena, D.; Hwang, J.; Walker, H. J. Phys. Chem. A 2005, 109, 10776. (18) Murlowska, K.; Sadlej-Sosnowska, N. J. Phys. Chem. A 2005, 109, 5590. (19) Saracino, G. A. A.; Improta, R.; Barone, V. Chem. Phys. Lett. 2003, 373, 411. (20) da Silva, C. O.; da Silva, E. C.; Nascimento, M. A. C. J. Phys. Chem. A 1999, 103, 11194. (21) Li, G.; Cui, Q. J. Phys. Chem. B 2003, 107, 14521. (22) Dahlke, E. E.; Cramer, C. J. J. Phys. Org. Chem. 2003, 16, 336. (23) Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. J. Phys. Chem. A 2003, 107, 9380. (24) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 2493. (25) Pliego, J. R., Jr.; Riveros, J. M. J. Phys. Chem. A 2002, 106, 7434. (26) Klicic, J. J.; Friesner, R. A.; Liu, S.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327. (27) Jang, Y. H.; Sowers, L. C.; Cagin, T.; Goddard, W. A., III. J. Phys. Chem. A 2001, 105, 274. (28) Himo, F.; Eriksson, L. A.; Blomberg, M. R. A.; Siegbahn, P. E. M. Int. J. Quantum Chem. 2002, 76, 714. (29) Himo, F.; Noodleman, L.; Blomberg, M. R. A.; Siegbahn, P. E. M. J. Phys. Chem. A 2002, 106, 8757. (30) Chipman, D. M. J. Phys. Chem. A 2002, 106, 7413. (31) Mujika, J. I.; Mercero, J. M.; Lopez, X. J. Phys. Chem. A 2003, 107, 6099. (32) Scharnagl, C.; Raupp-Kossmann., R. A. J. Phys. Chem. B 2004, 108, 477. (33) Lee, I.; Kim, C. K.; Lee, I. Y.; Kim, C. K. J. Phys. Chem. A 2000, 104, 6332. (34) Maksic, Z. B.; Vianello, R. J. Phys. Chem. A 2002, 106, 419. (35) Magill, A. M.; Cavell, K. J.; Yates, B. F. J. Am. Chem. Soc. 2004, 126, 8717. (36) Li, H.; Hains, A. W.; Everts, J. E.; Robertson, A. D.; Jensen, J. H. J. Phys. Chem. B 2002, 106, 3486. (37) da Silva, C. O.; da Silva, E. C.; Nascimento, M. A. C. J. Phys. Chem. A 2000, 104, 2402. (38) da Silva, C. O.; da Silva, E. C.; Nascimento, M. A. C. Chem. Phys. Lett. 2003, 381, 244. (39) Alvarado-Gonzalez, M.; Orrantia-Borunda, E.; Glossman-Mitnik, D. J. Mol. Struct. (THEOCHEM) 2008, 869, 105. (40) Sulpizi, M.; Sprik, M. Phys. Chem. Chem. Phys. 2008, 10, 5238. (41) Tao, L.; Han, J.; Tao, F. M. J. Phys. Chem. A 2008, 112, 775. (42) Kieseritzky, G.; Knapp, E. W. J. Comput. Chem. 2008, 29, 2575. (43) da Silva, G.; Kennedy, E. M.; Dlugogorski, B. Z. J. Phys. Chem. A 2006, 110, 11371. (44) Ulander, J.; Broo, A. Int. J. Quantum Chem. 2005, 105, 866. (45) Gross, K. C.; Seybold, P. G. Int. J. Quantum Chem. 2000, 80, 1107. (46) Gross, K. C.; Seybold, P. G. Int. J. Quantum Chem. 2001, 85, 569. (47) Pliego, J. R., Jr.; Riveros, J. M. Chem. Phys. Lett. 2000, 332, 597. (48) Pliego, J. R., Jr. Chem. Phys. Lett. 2003, 367, 145. (49) Adam, K. R. J. Phys. Chem. A 2002, 106, 11963. (50) Kim, Y. J.; Streitwieser, A. J. Am. Chem. Soc. 2002, 124, 5757. (51) Schuurmann, G.; Cossi, M.; Barone, V.; Tomasi, J. J. Phys. Chem. A 1998, 102, 6706. (52) Soriano, E.; Cerdan, S.; Ballesteros, P. J. Mol. Struct. (THEOCHEM) 2004, 684, 121. (53) Sadlej-Sosnowska, N. Theor. Chem. Acc. 2007, 118, 281. (54) Govender, K. K.; Cukrowski, I. J. Phys. Chem. A 2009, 113, 3639. (55) Private communication. Gaussian Technical Support. (56) Cramer, C. J. In Essentials of Computational Chemistry; John Wiley & Sons Ltd.: New York, 2002; pp 371-372. (57) NIST Standard Reference Database 46. NIST Critically Selected Stability Constants of Metal Complexes Database; Version 8.0; Smith, R. M., Martell, A. E., Eds.; US Department of Commerce, National Institute of Standards and Technology: Gaithersburg, MD, 2004. J J. Phys. Chem. A, Vol. xxx, No. xx, XXXX Govender and Cukrowski (58) Grigoriev, M. S.; Auwer, C. D.; Meyer, D.; Moisy, P. Acta Crystallogr. 2006, C62, m163. (59) Yu, L. C.; Liu, S. L.; Liang, E. X.; Wen, C. L. J. Coord. Chem. 2007, 60, 2097. (60) Malevich, D.; Wang, Z.; Tremaine, P. R. J. Sol. Chem. 2006, 35, 1303. (61) Vohra, M. S.; Davis, A. P. J. Colloid Interface Sci. 1997, 194, 59. (62) Karhu, J.; Harju, L.; Ivaska, A. Anal. Chim. Acta 1999, 380, 105. (63) Souaya, E. R.; Hanna, W. G.; Ismail, E. H.; Milad, N. E. J. Coord. Chem. 2004, 57, 825. (64) Quartacci, M. F.; Irtelli, B.; Baker, A. J. M.; Navari-Izzo, F. Chemosphere 2007, 68, 1920. (65) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Pettersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; 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.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision D.01; Gaussian, Inc.: Wallingford, CT, 2004. (66) GaussView 4.1.2; Gaussian Inc.; Wallingford, CT, 2004. (67) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (68) Cances, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032. (69) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Chem. Phys. Lett. 1996, 255, 327. (70) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (71) Barone, V.; Cossi, M. J. Phys. Chem. A 1998, 102, 1995. (72) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. J. Comput. Chem. 2003, 24, 669. (73) Allen, F. H. Acta Crystallogr. 2002, B58, 380. (74) Skrzypczak-Jankun, E.; Smith, D. A.; Maluszynska, H. Acta Crystallogr. 1994, C50, 1097. (75) Kaneyoshi, M.; Bond, A.; Jones, W. Acta Crystallogr. 1999, C55, 1260. (76) da Silva, G.; Bozzelli, J. W. J. Phys. Chem. A 2006, 110, 13058. (77) da Silva, G.; Bozzelli, J. W. J. Phys. Chem. C 2007, 111, 5760. (78) Guo, Y.; Gao, H.; Twamley, B.; Shreeve, J. M. AdV. Mater. 2007, 19, 2884. (79) Espinosa-Garcia, J.; Garcia-Bernaldez, J. C. Phys. Chem. Chem. Phys. 2002, 4, 4096. (80) Krossing, I.; Raabe, I. Chem.sEur. J. 2004, 10, 5017. (81) Kutt, A.; Movchun, V.; Rodima, T.; Dansauer, T.; Rusanov, E. B.; Leito, I.; Kaljurand, I.; Koppel, J.; Pihl, V.; Koppel, I.; Ovsjannikov, G.; Toom, L.; Mishima, M.; Medebielle, M.; Lork, E.; Roschenthaler, G. V.; Koppel, I. A.; Kolomeitsev, A. A. J. Org. Chem. 2008, 73, 2607. (82) Nolan, E. M.; Linck, R. G. J. Am. Chem. Soc. 2000, 122, 11497. (83) Peeters, D.; Leroy, G.; Wilante, C. J. Mol. Struct. 1997, 416, 21. (84) Wang, L.; Heard, D. E.; Pilling, M. J.; Seakins, P. J. Phys. Chem. A 2008, 112, 1832. (85) Schrodinger, L. L. C. Schrodinger Maestro; 32nd Floor, Tower 45, 120 West Forty-Fifth Street, New York, 10036, 2003. JP9092964 Protonation Constants for Nitrilotriacetic Acid J. Phys. Chem. A, Vol. xxx, No. xx, XXXX K