FULL PAPER

WWW.C-CHEM.ORG

Aqueous Acidities of Primary Benzenesulfonamides: Quantum Chemical Predictions Based on Density Functional Theory and SMD KeRstutis Aidas,*[a,b] Kiril Lanevskij,[a] Rytis Kubilius,[a] Liutauras Juska,[a] Daumantas Petkevicˇius,[b] and Pranas Japertas[a,c] Aqueous pKa of selected primary benzenesulfonamides are predicted in a systematic manner using density functional theory methods and the SMD solvent model together with direct and proton exchange thermodynamic cycles. Some test calculations were also performed using high-level composite CBS-QB3 approach. The direct scheme generally does not yield a satisfactory agreement between calculated and measured acidities due to a severe overestimation of the Gibbs free energy changes of the gas-phase deprotonation reaction by the used exchange-correlation functionals. The relative pKa values calculated using proton exchange method compare to experimental data very well in both qualitative and quantita-

tive terms, with a mean absolute error of about 0.4 pKa units. To achieve this accuracy, we find it mandatory to perform geometry optimization of the neutral and anionic species in the gas and solution phases separately, because different conformations are stabilized in these two cases. We have attempted to evaluate the effect of the conformer-averaged free energies in the pKa predictions, and the general conclusion is that this procedure is highly too costly as compared C 2015 with the very small improvement we have gained. V Wiley Periodicals, Inc.

Introduction

In drug discovery, it is highly necessary to have theoretical methods capable of predicting aqueous acidities of candidate compounds as precisely as possible. In case of sulfonamidebased drugs, the need for high precision is further emphasized by the fact that many compounds bearing a sulfonamide moiety tend to ionize around physiological pH. In this range, even a small change of pKa can lead to a significantly altered distribution of ionic forms and, therefore, drug’s activity. Achieving the desired accuracy of predictions still poses a significant challenge for empirical quantitative structure–property relationship (QSPR) approaches that are commonly used for this task.[15,16] The main reason is that such methods are heavily dependent on availability of a sufficiently large training set containing high quality experimental data, whereas ionization of sulfonamides has not been studied well enough to meet this requirement. Given the wide variation range of pKa of the sulfonamide group, there are not yet enough data to comprehensively describe the influence of the chemical environment of the ionizable center, as well as the modulating effects of

Compounds bearing a sulfonamide moiety are highly relevant from a pharmacological perspective. A number of representatives of this class of compounds are clinically used drugs generally exhibiting antiviral and antibacterial activity.[1] Examples of sulfonamide drugs include antimicrobials (sulfacetamide, sulfadiazine), diuretics (acetazolamide, hydrochlorothiazide), anticonvulsants (sultiame, topiramate), antidiabetics (glimepiride, tolbutamide), and other medications.[2,3] Over the last decade, sulfonamides have been actively investigated as promising inhibitors of carbonic anhydrases.[4,5] These metalloenzymes were found to exist in 15 different isoforms in the human body, and their enhanced activity has been linked to various diseases. In particular, carbonic anhydrases IX and XII have been associated with various forms of human cancer.[6–8] Because the differences between different isoforms are rather small, potential inhibitors ought to be highly selective, and various classes of sulfonamides have been demonstrated to possess this important feature.[6,9,10] The protonation state of a drug molecule has a profound influence on numerous physiological processes, such as the rate of intestinal absorption,[11] permeability across the bloodbrain barrier,[12] and further distribution within the body.[13] It has been long recognized that anionic forms of sulfonamides are generally more biologically active as compared with the neutral species.[14] The inhibitory activity of sulfonamides against carbonic anhydrases is precisely due to the anionic sulfonamides displacing solvent molecules in the active site of the enzyme, and the nitrogen atom of the sulfonamide moiety binds directly to the Zn(II) ion.[5,6] 2158

Journal of Computational Chemistry 2015, 36, 2158–2167

DOI: 10.1002/jcc.23998

[a] K. Aidas, K. Lanevskij, R. Kubilius, L. Juska, P. Japertas Vs˛I “Aukstieji algoritmai”, A. Mickevicˇiaus g. 29, LT-08117, Vilnius, Lithuania E-mail: [email protected] [b] K. Aidas, D. Petkevicˇius Department of General Physics and Spectroscopy, Faculty of Physics, _ Vilnius University, Sauletekio al. 9, LT-10222 Vilnius, Lithuania [c] P. Japertas ACD/Labs, Inc., 8 King Street East, Suite 107, Toronto, Ontario, Canada M5C 1B5 Contract grant sponsor: Research Council of Lithuania; Contract grant number: MIP-026/2013 C 2015 Wiley Periodicals, Inc. V

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

various substituent patterns. Another issue comes from discrepancies in experimental data published in different sources, as it is not uncommon to encounter values differing by more than 1 pKa unit reported for the same compound.[17] In this work, we focused on developing a quantum chemical approach that would enable us to gain more insight into aqueous dissociation of sulfonamide-based drugs. Our primary aim was to derive a computational procedure geared to accurately calculate acid pKa values of a prototypical set of primary benzenesulfonamides that was basically compiled from a single source[14] to ensure consistence of experimental data. Our approach is based on electronic structure methods rooted in density functional theory (DFT), continuum solvation model and thermodynamic cycles. The thermodynamic cycle-based methods present an effective way to indirectly evaluate the condensedphase free energy change of the acid dissociation reaction, a quantity which solely determines the acidity constant.[18–21] We will consider two methods from the family of thermodynamic cycles—the direct method and the so-called proton exchange scheme. The direct method is by far the simplest of all thermodynamic cycles. However, even though gas-phase acidities can in principle be evaluated with the accuracy of several kJ/mol using, for example, different composite[22,23] or high-level Wn methods,[24] accurate prediction of solvation free energies, especially those of charged compounds, is still a severe issue limiting the accuracy of the direct thermodynamic cycle.[18,20] Here, these two quantities are combined with the experimental-theoretical value of the proton solvation free energy, meaning that the success of the direct method relies heavily on error cancelation. Even though the longstanding debate[18,20,25–28] on the correct value of the proton solvation free energy has converged to a consensus (we refer to the chapter of Alongi and Shields[20] for an excellent discussion of this topic), the direct scheme is nevertheless prone to systematic errors. To deal with these issues, there are two approaches to be generally adopted. One can rely either on QSPRs where pKa is evaluated using certain descriptors typically based on gasphase quantum chemical calculation,[29–32] or on more complex thermodynamic cycles.[18–20,25,33–38] One of the latter methods is the so-called proton exchange method which involves an isodesmic reaction where the proton migrates from the dissociating acid molecule to the conjugate base of another (reference) compound. Here, an explicit reference to the free proton is absent, thus avoiding the need for any related experimental data. In addition, an ionic compound is present on both sides of the isodesmic reaction, and this leads to a potential cancelation of errors in computing the differences in the solvation free energies of neutral and, particularly, ionic species. This cancelation of errors is more likely in the case where the reference compound is structurally similar—at least in the vicinity of the ionizable moiety—to the compound under study. Because this method allows evaluating the relative pKa value with respect to that of the reference only, the accuracy of the predicted pKa thus also depends on how reliable the pKa of the reference is. In this work, we rely on a universal solvation model based on solute electronic density, SMD, which has been

parameterized on an extensive set of mainly organic molecules[39] and thoroughly benchmarked,[40] to compute solvation free energies of benzenesulfonamides. Compounds containing the sulfonamide moiety have not, however, been included in the training set of SMD, and the performance of SMD to predict pKa values of this important class of pharmaceutical molecules is thus a study of its own interest. The sulfonamides we consider in this work are medium-sized molecules possessing a degree of conformational flexibility (vide infra). Proper conformational analysis is thus mandatory to identify low-energy conformers, yet certain constraints have to be imposed on the used level of theory to maintain cost-effectiveness of our scheme. We will thus rely on the B3LYP/6-31G* model chemistry[41–43] which was used in the parameterization of SMD. In addition to testing two different thermodynamic cycles, we will evaluate the effect of different exchange-correlation functional, M05-2X,[44] in combination with one-electron basis sets of different quality on the accuracy of pKa predictions. We will also inspect whether reasonable pKa values of benzenesulfonamides can be expected using free energies of the lowest-energy conformers only, or whether proper conformeraveraged free energies are required. The structural features and gas-phase acidities of various sulfonamide classes have been a subject of several previous theoretical studies.[45–49] Quantum chemical calculations of aqueous acidities of sulfonamides have been reported,[50–52] but the results have not been analyzed in depth. To the best of our knowledge, this is the first study where aqueous acidities of sulfonamides are predicted and analyzed in a systematic manner.

Computational Details The acidity constant, pKa , is a measure of the thermodynamic feasibility of acid dissociation reaction in solution. The standard-state free energy change of the acid dissociation reaction, DGaq , is associated with the aqueous pKa according to[18,20] pKa ¼

DGaq ; RTln10

(1)

where R and T are the universal gas constant and temperature, respectively. Using the direct thermodynamic cycle depicted in Figure 1, the free energy change of the acid dissociation reaction at the standard state of 1 mol/l is expressed as follows: 



DGaq ¼ DGg 1DDGs 1DG

!

(2)

: 

The first term on the right-hand side of eq. (2), DGg , is the Gibbs free energy change of the acid dissociation reaction in the gas phase for the standard state of 1 atm. The term DDGs is the difference between the solvation free energies of the products and the reactant, and is given according to DDGs ¼ DGs ðA- Þ1DGs ðH1 Þ2DGs ðAHÞ:

(3)

Here, the terms DGs ðAHÞ; DGs ðA- Þ and DGs ðH1 Þ are the solvation free energies of the dissociating acid, its conjugate base and the proton, respectively, at the standard state of 1 mol/l. Journal of Computational Chemistry 2015, 36, 2158–2167

2159

FULL PAPER

WWW.C-CHEM.ORG

Figure 1. The direct thermodynamic cycle. 

The last term on the right-hand side of eq. (2), DG ! , accounts for the conversion between the standard states of 1 atm and 1 mol/l, and is evaluated[20] as  !

DG

5 RTln24:4654

(4)

At the temperature of 298.15 K, this term amounts to 7.9 kJ/mol. An isodesmic reaction is considered in the proton exchange thermodynamic cycle as shown in Figure 2, where an auxiliary compound, BH, denoted “reference,” adopts the proton from the dissociating acid, AH. This scheme allows for the evaluation of the relative pKa of the acid AH with respect to that of the reference according to pKa ðAHÞ ¼

DGaq 1pKa ðBHÞ: RTln10

(5)

In the latter equation, DGaq is the standard-state free energy change of the isodesmic reaction in aqueous solution given as 

DGaq ¼ DGg 1DDGs : 









DGg ¼ Gg ðA2 Þ1Gg ðBHÞ2Gg ðAHÞ2Gg ðB2 Þ;

(7)



where Gg ðXÞ is the gas-phase free energy of compound X at the standard state of 1 atm. The second term on the righthand side of eq. (6) is again the difference in the solvation free energies of the products and reactants, and reads DDGs ¼ DGs ðA2 Þ1DGs ðBHÞ2DGs ðAHÞ2DGs ðB2 Þ:

(8)

In case of both direct and proton exchange methods, the gas-phase Gibbs free energy change of the considered reaction was estimated by evaluating the Gibbs free energy of each involved species which includes the terms of electronic, Eelec, and zero-point vibrational energies, EZPV, thermal contributions to enthalpy, Hcorr, and entropy, S, according to Gg ¼ Eelec 1EZPV 1Hcorr 2TS:

Gave

 # X Gi 2G0 ¼ G0 2RT ln 11 exp 2 ; RT i6¼0

(10)

Journal of Computational Chemistry 2015, 36, 2158–2167

where G0 is the lowest free energy among those of all conformers, and the summation runs over all conformers excluding the most stable. The conformer-averaged gas-phase free energy is evaluated using the Gibbs free energies of every gas-phase conformer computed by eq. (9). The conformeraveraged solvation free energy is calculated in this work as follows: DGs;ave ¼ Gaq;ave 2Eelec;ave :

(11)

The first term on the right-hand side of eq. (11), Gaq;ave , is the conformer-averaged liquid-phase free energy, and the averaging is performed over the conformations found in aqueous solution. The second term, Eelec;ave , is the conformeraveraged electronic energy, computed analogously using eq. (10), but the averaging is performed over the conformations found in vacuum which in principle may differ from those found in solution. Only those conformations with energies differing from G0 (or Eelec;0 in the case when Eelec;ave is computed)

(9)

The zero-point energy, entropy, and thermal contributions to enthalpy are estimated using standard ideal gas rigid-rotor harmonic-oscillator approximation.[18,53] The gas-phase free energy of the proton is 226.3 kJ/mol.[18] 2160

"

(6)

In contrast to eq. (2), the term DG ! is absent in eq. (6), because the number of reactants is now equal to that of products,[20] see Figure 2. The free energy change of the isodesmic reaction in the gas phase is expressed according to 

The solution-phase free energies are calculated using the SMD model[39] which belongs to the class of implicit solvation models. The solvation free energy is evaluated as the difference between the solution-phase free energy and the electronic energy of the compound computed in vacuum, as this approach is generally adopted in the parameterization of the implicit solvation models aimed to provide accurate solvation energies. A value of 21112.5 kJ/mol[18,20] is used for the solvation free energy of the proton in this work, as the same value was actually utilized to compile the set of solvation free energies for the parameterization of the SMD model.[39] The SMD model was formulated for geometry-optimized molecules in the gas phase, meaning that the relaxation effects are neglected in this scheme. For single-conformer molecules or in the cases where the same single conformer dominates in both gas and solution phases, this is a sound approximation. In this work, we have chosen to account for the geometry relaxation by optimizing the geometry of all species in both gas phase and aqueous solution separately, the latter being modeled via SMD. In the case, where several conformations of comparable free energies are present, conformer-averaged free energies should in principle be used while calculating the solution phase free energy changes according to eqs. (2), (3), and (6)–(8). The Boltzmann average of the free energy is generally evaluated according to[18,20]

Figure 2. The proton exchange thermodynamic cycle.

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

N1-methylsulfanilamide (2), N1-o-tolylsulfanilamide (3), N1-mtolylsulfanilamide (4), N1-p-aminophenylsulfanilamide (5), N1hydroxyethylsulfanilamide (6), and N1-furfurylsulfanilamide (7). To predict aqueous pKa of these molecules, we have first undertaken an extensive conformational analysis aimed to identify lowest-energy conformers of neutral and deprotonated species in gas and solution phases. In addition, we attempted to find all conformers similar in energy, to assess the effect of using conformer-averaged free energies to predict acidities rather than those of single most stable conformers. For this purpose, we have relied on the B3LYP exchangecorrelation functional and the Pople type 6-31G* basis set as well as SMD model for solvent effects. The ranking of the gasphase conformations is based on their electronic energies. We have focused only on the conformations of similar energies, so that the energy differences do not exceed 5 kJ/mol. Molecular geometries of gas- and solution-phase conformations of all compounds 1–7 in their neutral and deprotonated forms considered in this work are collected in the Supporting Information, including all relevant thermochemical data for each molecule as well as conformer-averaged energies. Conformational analysis Figure 3. A set of benzenesulfonamides under study.

by no more than 5 kJ/mol significantly contribute to the averaged energy.[54] All electronic structure calculations in this work were conducted using Gaussian 09 program.[55] To utilize SMD model consistently, we have performed calculations of solvation free energies at the same levels of theory that were utilized in the parameterization of SMD. This includes the B3LYP and M05-2X exchange correlation functionals[41,44] along with the 6-31G*, 6311G**, or cc-pVTZ basis sets.[42,43,56] We have utilized the pruned (99,590) grid for the evaluation of the two-electron integrals as implemented in Gaussian 09. In several instances, CBS-QB3 method[23] was used to accurately predict gas-phase free energies of the considered species. The geometry optimization in both gas and solution phases was performed using tight convergence criterion which was always followed by the calculation of harmonic vibrational modes. This procedure was applied to verify whether the true minimum on the potential energy surface was found and, in the case of the gas phase, to evaluate zero-point energy, as well as thermal contribution to the gas-phase Gibbs free energy. Normally, unscaled vibrational frequencies were utilized. To account for the vibrational anharmonicity, we have also considered the uniform scaling of the B3LYP/6-31G*-based harmonic frequencies by a factor of 0.9614 as derived by Scott and Radom.[57] Solutionphase calculations were performed using the default settings of SMD as implemented in Gaussian 09. All pKa values were computed at the temperature of 298.15 K.

Results and Discussion Our selection of benzenesulfonamides considered in this work is presented in Figure 3 and includes benzenesulfonamide (1),

The low-energy conformations stemming from the benzenesulfonamide moiety in its neutral and deprotonated forms along with the corresponding labels are presented in Figure 4. The neutral form of benzenesulfonamide, 1, possesses two conformations of Cs symmetry in the gas phase, namely the eclipsed and the staggered conformations, denoted SN1 and SN2 in Figure 4. The B3LYP/6-31G* approach indicates the eclipsed conformation to be more stable, as the electronic energy difference between the two is around 5 kJ/mol, whereas the difference in the Gibbs free energies is even larger and amounts to around 10 kJ/mol. These findings are in line with the results of the gas electron diffraction measurements[58] and previous quantum chemical calculations,[46–48,50,51,58] which indicate a clear preference of eclipsed over the staggered conformation. In aqueous solution, only the staggered conformation survives as the eclipsed conformation turns out to be a saddle point

Figure 4. Conformations of the sulfonamide moiety in its neutral and deprotonated forms as well as their labels. Different atom types are represented by the following colors: blue, nitrogen; yellow, sulfur; red, oxygen; brown, carbon; white, hydrogen.

Journal of Computational Chemistry 2015, 36, 2158–2167

2161

FULL PAPER

WWW.C-CHEM.ORG

Figure 5. Unique conformations of 3 for a given conformation of the sulfonamide moiety. Different atom types are represented by the following colors: blue, nitrogen; yellow, sulfur; red, oxygen; brown, carbon; white, hydrogen.

on the B3LYP/6-31G*/SMD-based potential energy surface. The staggered conformation is stabilized in the condensed phase due to its much larger permanent dipole moment, which was calculated to be 6.08 D in gas phase, compared with 3.66 D for the eclipsed conformation. Two unique conformations were found for the deprotonated form of 1 in the gas phase. They correspond to the SD1 and SD2 conformations shown in Figure 4, and are virtually iso-energetic with a minor difference in the electronic energies of around 0.1 kJ/mol slightly preferring the SD2 conformation. The atomic arrangement of the deprotonated sulfonamide moiety with the proper dihedral HANASAC angle of around 908 has been identified as a global minimum previously.[45] These two conformations are also stabilized in aqueous solution, and the difference of 0.4 kJ/mol between their free energies is again very small, but now the SD1 conformation is slightly more stable. Therefore, the SD1 and SD2 conformers were considered in conformational averaging of the free energies of the deprotonated forms of 1. Attachment of the substituted amino group at the para position of the benzene ring in the benzenesulfonamide breaks the Cs symmetry characteristic to the neutral form of 1. This has an immediate consequence that the eclipsed conformation is no longer the minimum on the B3LYP/6-31G*-based potential energy surface, both in vacuum and solution. Instead, a new minimum emerges, denoted SN3 in Figure 4. We will refer to this conformation of the sulfonamide moiety as “twisted” further in the text. In vacuum, the twisted conformation is always more stable than the staggered conformation for all considered compounds (2–7) with electronic energy differences typically being around 4.7 kJ/mol. Thereby, in all cases we consider only the twisted conformations of the sulfonamide moiety for the conformer-averages of the gas-phase (free) energies. In the condensed phase, the staggered conformation is found to be the most stable again, but—in contrast to benzenesulfonamide 1—the twisted conformations are found to coexist as well. The free energy differences between the two forms are typically around 3 kJ/mol, and thus both forms have to be included in the conformational averaging of the free energies according to eq. (10). Two types of confor2162

Journal of Computational Chemistry 2015, 36, 2158–2167

mations of the deprotonated sulfonamide moiety, SD1 and SD2, are found to have similar electronic energies in the gas phase, and the electronic energy differences between the two does not exceed 1 kJ/mol. In the condensed phase, the same two conformations predominate and possess similar free energies again. In addition, a third conformational form, denoted SD3 in Figure 4, is found to be stabilized in solution. The free energy difference between the SD1 and SD3 is typically around 1.8–2.4 kJ/mol, therefore, we have always included this conformation into the averaging of free energies. We have to note, however, that the SD3 conformational form was not always found, as it evolves down to conformation SD1 quite often. In addition, several other conformations of the deprotonated sulfonamide moiety have been found in solution. In particular, they differ from the conformations depicted in Figure 4 by the location of the hydrogen atom attached to nitrogen, which now faces toward the viewer in the orientation shown in Figure 4. In all cases, the free energies of the conformations of this type were estimated to be at least 5 kJ/mol higher compared with the most stable conformation. Consequently, we have always neglected their presence when calculating the conformer-averaged free energies. Conformations of the compound 2 arise solely from the conformational flexibility of the sulfonamide moiety because the methylamino group at the para position of the benzene ring is virtually planar. However, the attachment of this group breaks the Cs symmetry present in 1, so that the number of unique conformations for 2 increases. The neutral form possesses three unique conformations, two twisted and one staggered in both gas and solution phases. For the deprotonated form of 2, there are four conformations in vacuum, two of SD1 and two of SD2 type, and five in condensed phase, including one extra SD3 type conformation. For compounds 3, 4, 6, and 7, additional conformational degrees of freedom originate from the conformational flexibility of the substituent attached to the nitrogen atom of the amino group. There are two unique conformational states of the tolyl moiety in 3 and 4, see Figure 5. For 3, the conformational state A depicted in Figure 5 is always considerably lower in energy than state B. Therefore, we have discarded all conformers of 3 related to the conformational state B of the tolyl moiety. This leaves us with 6 unique conformations of the neutral form of 3 in both vacuum and solution, whereas 8 unique conformations of SD1 and SD2 type in vacuum are supplemented by 2 SD3 conformations in the condensed phase for the deprotonated form of 3. In contrast, both conformations A and B shown in Figure 5 have very similar energies for a given conformation of the sulfonamide moiety in 4 due to the altered steric effects as the methyl group is now attached to the meta position of the benzene ring. The number of conformations considered in the averaging of the energies of 4 is thus virtually twice that of 3. Rotation of the aminophenyl moiety around the C–N bond does not introduce any additional conformations of 5. However, there is a difference in the orientation of this group with respect to the benzene ring of the benzenesulfonamide found for the deprotonated form in vacuum as compared with all other cases. As shown in Figure 6, the aminophenyl group WWW.CHEMISTRYVIEWS.COM

WWW.C-CHEM.ORG

Figure 6. Conformation of 5 in a) neutral form in both gas and solution phases as well as in deprotonated form in solution phase, and b) in deprotonated form in gas phase. Different atom types are represented by the following colors: blue, nitrogen; yellow, sulfur; red, oxygen; brown, carbon; white, hydrogen. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

possesses an in-plane orientation, A, in the neutral form and in the deprotonated form in solution, but an out-of-plane conformation, B, was found to be stabilized for the deprotonated form in vacuum. Similar out-of-plane orientation of the tolyl moiety is found for the deprotonated forms of 3 and 4 in vacuum as well. Compounds 6 and 7 exhibit the richest conformational variety among the considered molecules. The low-energy conformations of the hydroxyethyl moiety in 6 are shown in Figure 7. The A and B states (see Figure 7) are found to be the preferred conformations of this moiety for the neutral form of 6

Figure 7. Low energy conformations of 6 for a given conformation of the sulfonamide moiety. Different atom types are represented by the following colors: blue, nitrogen; yellow, sulfur; red, oxygen; brown, carbon; white, hydrogen. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FULL PAPER

in vacuum, the electronic energy of the former being lower by about 1.0–1.5 kJ/mol for a given conformation of the sulfonamide group. Both states A and B are apparently stabilized by the favorable electrostatic interaction between the oxygen atom of the hydroxy group and the hydrogen atom of the amino group. In solution, the A configuration remains the most stable, yet a new type of conformation, C, emerges with free energy very similar to that of B. In the conformational averaging of the liquid-phase free energies of neutral form of 6, we have included 6 unique twisted and staggered conformations corresponding to the A conformational state, as well as 4 staggered conformations corresponding to the B and C states of the hydroxyethyl moiety. The most stable anionic forms of 6 in vacuum, D, E, and F as depicted in Figure 7, are quite different from A, B, or C. However, these configurations are also seemingly stabilized by electrostatic interactions, the only difference is that now the hydrogen atom of the hydroxy group is oriented toward the nitrogen atom of the amino group. We have considered a total of 16 conformations of D and E types owing to the 8 unique low-energy SD1 and SD2 conformations of the sulfonamide moiety for each of the two cases. In addition, four conformations of F type were included in the averaging because the other four were found to evolve to the corresponding conformations of the D type. In aqueous solution, the A, B, and C conformations prevail for the deprotonated form of 6. We have included 9 unique conformations of the A type, as 1 SD3 type conformer was found along with 8 SD1 and SD2 conformations, and also 8 conformations of B and C types (4 SD1 type conformations for each respective state of the hydroxyethyl group), yielding a total of 17 unique conformations. The low-energy conformations of the furfuryl moiety of compound 7 are shown in Figure 8. For the neutral form in vacuum, the A type conformation is found to be the most stable, which is stabilized by around 1.4 and 2.4 kJ/mol against the conformations of type C and B, respectively. The conformational averaging of the gas-phase energies of neutral forms of 7 includes 12 conformations in total, all of which correspond to the twisted conformation of the sulfonamide group. In the condensed phase, the neutral form possesses 4 different lowenergy conformations corresponding to types C, D, A, and E, as shown in Figure 8. Both staggered and twisted conformations corresponding to C and D configurations of the furfuryl moiety were considered in the averaging of the free energies together with staggered conformations of A and E types, yielding a total of 16 conformations to be averaged over. The staggered conformation of the sulfonamide and the C state of the furfuryl group comprise the most stable conformation of neutral compound 7 in solution. The deprotonated form has as many as 40 unique conformations of similar energies, including A, B, C, E, and F states of the furfuryl group as shown in Figure 8. Here, SD1 conformation of the sulfonamide group and B conformation of the furfuryl moiety are identified as the global minimum on the potential energy surface. For the deprotonated form in solution, we have found 37 conformations of similar energies, and here the furfuryl moiety favors A, C, D, and E conformational states. The most stable conformation corresponds to the Journal of Computational Chemistry 2015, 36, 2158–2167

2163

FULL PAPER

WWW.C-CHEM.ORG

Figure 8. Low energy conformations of 7 for a given conformation of the sulfonamide moiety. Different atom types are represented by the following colors: blue, nitrogen; yellow, sulfur; red, oxygen; brown, carbon; white, hydrogen. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

C state of the furfuryl group and SD1 configuration of the sulfonamide moiety. Aqueous acidities The aqueous pKa of benzenesulfonamides 1–7 predicted using the direct thermodynamic cycle along with a listing of methods and strategies used for the evaluation of the gas-phase free energy changes of the dissociation reaction and solvation free energies are given in Table 1. In addition to the B3LYP predictions, we have also performed computations of acidities using the M05-2X exchange-correlation functional in combination with three different basis sets. We have assumed that the most stable conformations identified during the B3LYP/6-31G* conformational analysis will remain the most stable when

using the M05-2X functional. We have generally found that the B3LYP/6-31G*-based conformations of the benzenesulfonamides were also minima on the M05-2X potential energy surfaces. The only exception we have observed is related to the SD2 conformational state of the sulfonamide moiety in the deprotonated species in vacuum which was found to evolve to the SD1 type during the M05-2X/cc-pVTZ-based geometry optimization. In addition, we have obtained an imaginary vibrational frequency for the deprotonated forms of 1 in solution using the M05-2X functional together with 6-311G** as well as cc-pVTZ basis sets. Apparently, an even finer grid for the computation of the two-electron integrals is required here. Because these imaginary frequencies are small in magnitude and they are generally not required to compute solvation free energy, we have still considered the corresponding liquidphase geometries of 1 as minima, and the results of the computed pKa values are included in Table 1. As one can see in Table 1, most of the applied approaches significantly overestimate the pKa values of sulfonamides as compared with experimental data. In an attempt to investigate the reason of the failure of DFT-based predictions, we have performed highly accurate CBS-QB3 calculations of the free energy changes of the gas-phase dissociation reaction for selected sulfonamides. The results of the CBS-QB3 calculations are presented in Table 2 along with the B3LYP and M05-2Xbased results for the gas-phase free energy changes, as well as the differences in the solvation free energies estimated according to eq. (3). It is immediately evident from Table 2 that the B3LYP/6-31G* approach overestimates the gas-phase free energy change by around 50 kJ/mol as compared with the corresponding CBS-QB3 result. Previously, a very good agreement between CBS-QB3 and B3LYP results for gas-phase acidities has been reported.[47] The discrepancy found here can be attributed to the much smaller basis set used in our DFT calculations, whereas the 6-3111G** basis was utilized in the computations by Remko.[47] Indeed, the M05-2X functional in combination with the 6-31G* basis set overestimates this quantity by around 37 kJ/mol, but augmenting the basis set with the diffuse functions reduces this discrepancy by about

Table 1. Aqueous pK a values of primary benzenesulfonamides 1–7 calculated using the direct method. 

DGg B3LYP/6-31G* CBS-QB3 B3LYP/6-31G*[a] B3LYP/6-31G* B3LYP/6-31G*[c] B3LYP/6-31G*[d] M05-2X/6-31G* M05-2X/6-311G** M05-2X/cc-pVTZ Exp.[e]

DGs

1

2

3

4

5

6

7

B3LYP/6-31G* B3LYP/6-31G* B3LYP/6-31G* B3LYP/6-31G*[b] B3LYP/6-31G* B3LYP/6-31G*[d] M05-2X/6-31G* M05-2X/6-311G** M05-2X/cc-pVTZ

20.86 11.91 21.11 20.24 26.38 20.57 17.70 13.09 14.94 10.10

22.19 13.44 22.44 22.22 28.10 21.85 19.32 15.51 16.74 10.77

21.54 – 21.80 21.85 27.53 21.14 18.92 14.89 16.09 9.96

21.48 – 21.74 21.98 27.52 21.03 18.14 14.34 15.27 9.74

21.86 – 22.11 22.68 27.90 21.49 19.51 15.10 16.43 10.22

21.87 13.16 22.12 23.93 27.96 21.64 18.74 14.94 16.36 10.92

21.83 13.01 22.08 22.62 27.82 21.32 19.06 14.89 16.42 10.88



The levels of theory used in the calculation of free energy change of the acid dissociation reaction in the gas phase, DGg , and the free energy of solvation, DGs , are indicated in the first and second columns, respectively. The pKa calculations were performed using (free) energies of single lowest energy  conformers, unless specified otherwise. [a] DGg computed using scaled harmonic vibrational frequencies, scaling factor of 0.9614 was used.[57] [b] DGs  calculated for gas-phase optimized lowest-energy molecular geometries. [c] DGg approximated by the difference in the electronic energies, DEelec , that  is zero-point energy and thermal contributions are here neglected. [d] DGg and DGs calculated using conformer-averaged (free) energies. [e] Ref. [14] for 2–7, Ref. [60] for 1.

2164

Journal of Computational Chemistry 2015, 36, 2158–2167

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 2. Free energy change of the acid dissociation reaction in the gas  phase, DGg , and the difference in the solvation free energies, DDGs , as defined in eq. (3), calculated using different model chemistries and the SMD for lowest energy conformers of sulfonamides 2, 6, and 7 as found using the B3LYP/6-31G* level. All values are given in kJ/mol. 2 Method CBS-QB3 B3LYP/6-31G* M05-2X/ 6-31G* M05-2X/ 6-311G** M05-2X/ cc-pVTZ



6

DGg

DDGs

1411.5 1461.4 1448.9



7 

DGg

DDGs

DGg

DDGs

– 21341.4 21344.8

1400.4 1450.1 1436.2

– 21333.2 21337.2

1406.1 1456.4 1443.6

– 21339.7 21342.7

1421.4

21340.4

1409.3

21331.9

1414.4

21337.4

1428.1

21340.1

1419.0

21333.4

1421.6

21335.8

27 kJ/mol. The use of the extensive cc-pVTZ basis does not lead to further improvement, possibly due to the lack of the diffuse functions in this basis set. The differences in the solvation free energies, conversely, are almost insensitive to the level of theory used, and they vary in the range of 4–7 kJ/mol among different model chemistries. We are thus left with the conclusion that overprediction of pKa values shown in Table 1 is due to the overestimated free energy changes of the gasphase dissociation reaction of the benzenesulfonamides as computed using DFT methods. There is a clear correlation between the data provided in  Tables 1 and 2, that is, the more the value of the DGg is overestimated, the larger is the discrepancy between the predicted and experimental acidities. The best quantitative agreement between the calculated and measured pKa is thus seen when CBS-QB3 gas-phase free energy changes are utilized. We have used different approaches to calculate the acidities using the B3LYP/6-31G* level of theory. First, we have used scaled vibrational frequencies to compute zero-point energy as well as enthalpic and entropic contributions to gas-phase Gibbs free energies. The obtained pKa values are seen to be systematically increased in magnitude by 0.25–0.26 pKa units as compared with the results where unscaled frequencies were

considered. Second, we have neglected the geometry relaxation effects in solution by evaluating the solvation free energies on gas-phase optimized geometries of sulfonamides. As shown in Table 1, this approach leads to very small changes of the computed acidities, especially when compared with the large discrepancies between B3LYP results and experimental pKa values. The effect of using the conformer-averaged free energies instead of the energies of single conformers is also very small. Conversely, the zero-point vibrational energy and thermal corrections significantly contribute to the absolute values of predicted acidities, and these contributions cannot be neglected when using the direct thermodynamic cycle. The aqueous acidities computed using the proton exchange method are presented in Table 3. Here, compound 2 and its experimental pKa value were chosen as a reference. In contrast to the direct method-based predictions, relative pKa values compare to the experimental results very well, especially for compounds 3–7. Using the B3LYP/6-31G* level of theory, the predicted acidities deviate from the experimental results by around 0.1–0.3 pKa units for compounds 3, 4, and 5, by around 0.5 pKa units for 6 and 7, and by around 0.7 pKa units for 1. The mean absolute error, MAE, of the pKa prediction for these six compounds is thus 0.4 pKa units. In particular, we also observe a rather correct qualitative picture, as the predicted acidities now quite closely follow the trend of experimental results for compounds 3–7 despite the fact that measured pKa values are distributed within a very narrow range of 1 pKa unit, compound 1 being the most notable exception. When CBS-QB3-based gas-phase free energies are used, virtually no improvement in the predicted pKa values of the considered compounds was found, indicating that the error cancelation is effective already at the DFT level. The scaling of the vibrational frequencies is seen to have no effect on the relative acidity constants, because the effect is extremely systematic and thus tends to fully cancel out. We observe that the prediction quality deteriorates when solvation free energies are estimated using gas-phase optimized geometries in the SMD calculations of the liquid-phase free energies. This

Table 3. Aqueous pK a values of primary benzenesulfonamides 1 and 3–7 calculated using the proton exchange method taking 2 with pK a of 10.77 as a reference. 

DGg B3LYP/6-31G* CBS-QB3 B3LYP/6-31G*[a] B3LYP/6-31G* B3LYP/6-31G*[c] B3LYP/6-31G*[d] M05-2X/6-31G* M05-2X/6-311G** M05-2X/cc-pVTZ Exp.[e]

DGs

1

3

4

5

6

7

MAE

B3LYP/6-31G* B3LYP/6-31G* B3LYP/6-31G* B3LYP/6-31G*[b] B3LYP/6-31G* B3LYP/6-31G*[d] M05-2X/6-31G* M05-2X/6-311G** M05-2X/cc-pVTZ

9.44 9.23 9.44 8.79 9.06 9.49 9.15 8.35 8.97 10.10

10.13 – 10.13 10.40 10.21 10.07 10.37 10.14 10.12 9.96

10.06 – 10.06 10.53 10.20 9.94 9.59 9.61 9.29 9.74

10.44 – 10.44 11.23 10.57 10.41 10.96 10.36 10.46 10.22

10.45 10.48 10.45 12.49 10.63 10.56 10.19 10.21 10.39 10.92

10.41 10.34 10.41 11.17 10.49 10.24 10.51 10.14 10.44 10.88

0.39 0.62 0.39 0.90 0.46 0.35 0.56 0.61 0.49



The levels of theory used in the calculation of free energy change of the acid dissociation reaction in the gas phase, DGg , and the free energy of solvation, DGs , are indicated in the first and second columns, respectively. The pKa calculations were performed using (free) energies of single lowest energy  conformers, unless specified otherwise. Mean absolute error, MAE, with respect to the experimental data is given in the last column. [a] DGg computed using scaled harmonic vibrational frequencies, scaling factor of 0.9614 was used.[57] [b] DGs calculated for gas-phase optimized lowest-energy molecular  geometries. [c] DGg approximated by the difference in the electronic energies, DEelec , that is zero-point energy and thermal contributions are here  neglected. [d] DGg and DGs calculated using conformer-averaged (free) energies. [e] Ref. [14] for 2–7, Ref. [60] for 1.

Journal of Computational Chemistry 2015, 36, 2158–2167

2165

FULL PAPER

WWW.C-CHEM.ORG

finding originates from the fact that different conformations of sulfonamides are found to be the most stable in gas and solution phases. Even in the cases of 1 and 2, where the conformational variability is attributed to the sulfonamide moiety alone, different conformations of the neutral form, eclipsed and staggered, are stabilized in the gas and solution phases, respectively. Geometry relaxation of the deprotonated form was also found to contribute considerably to the solvation free energy, as again different conformations, SD1 and SD2, were identified as the most stable in the two phases. Consequently, we find it mandatory to perform geometry optimizations of sulfonamides in the gas and solution phases separately. Such a procedure was first found to be critical for accurate predictions of aqueous acidities for a family of substituted phenols by Liptak et al.[59] Even though zero-point vibrational energy and thermal corrections have a large influence on the absolute magnitude of pKa predicted using the direct method, relative pKa values are much less sensitive to this approximation as seen in Table 3, and thus, costly calculations of the Hessian can in principle be omitted. However, caution is still necessary, as neglecting these contributions leads to a decrease in pKa prediction accuracy for 1 by 0.5 pKa units, indicating that the magnitude of this effect is likely dependent on the nature of the reference compound. In particular, the high accuracy of pKa predictions is partly due to the fact that zero-point energy and thermal contributions cancel out between similar sulfanilamides, but this cancelation is apparently not so successful in the case of benzenesulfonamide 1. Even though 1 may be seen as a natural reference compound to be used in the proton exchange method, our results thus indicate that caution has to be exercised when relative pKa values of primary benzenesulfonamides against 1 are computed. Higher degree of similarity between the benzenesulfonamide compounds under scrutiny and the reference compound, likely beyond that of the benzenesulfonamide moiety alone, may be generally required to achieve the accuracy similar to that seen in Table 3 for sulfanilamides with 2 used as reference. As one can see in Table 3, the effect of using conformeraveraged energies for pKa predictions is indeed small for all considered benzenesulfonamides, and this daunting procedure is to a large extent unnecessary. The effect is especially insignificant for compounds 3–5. Even though inclusion of this procedure could be expected to improve predicted pKa values of 6 and 7, owing to the conformational richness of these molecules, a slight improvement was only observed for 6, but not for 7. Using the M05-2X functional leads to a slightly worse overall pKa prediction accuracy compared with the B3LYP predictions using the 6-31G* basis set, although the agreement between experimental and calculated values is still very satisfactory, and the prediction error is below 1 pKa unit for most compounds. The inclusion of the diffuse functions in the basis set generally enables us to obtain slightly better results for 3–6, but the improvement is definitely not so dramatic as was the case with the direct method.

Summary In this work, we attempted to gain more insight into the dissociation of primary benzenesulfonamides—a class of molecules with 2166

Journal of Computational Chemistry 2015, 36, 2158–2167

important pharmaceutical properties—using quantum chemical approaches and thermodynamic cycle methods in a systematic manner. We have mainly relied on the electronic structure methods rooted in DFT and the SMD model to account for solvation. The direct and the proton exchange schemes from the family of thermodynamic cycle methods were considered for pKa predictions. We have performed a thorough conformational analysis of neutral and deprotonated forms of seven different benzenesulfonamide compounds in gas and aqueous solution phases. The conformational analysis is based on the B3LYP functional and the 6-31G* basis set along with SMD. For most of the considered molecules, different gas- and solution-phase conformers were identified as global minima on the respective potential energy surfaces. The absolute magnitudes of pKa values predicted using the direct thermodynamic cycle significantly deviate from the experimental measurements. The cause of this disagreement is related to the failure of DFT methods to predict gas-phase free energy changes of the acid dissociation reaction accurately, as the utilized DFT methods in combination with various basis sets overestimate this quantity by tens of kJ/mol as compared with the CBS-QB3 results. The relative pKa values predicted using the proton exchange method compare to the experimental data very well: the deviations are mostly well below 1 pKa unit, and qualitative trends are correctly reproduced as well. To achieve such accuracy, we find it mandatory to perform geometry optimization of the neutral and deprotonated species in the gas and solution phases separately, because different conformations are stabilized in these two cases. However, while using N1-methylsulfanilamide as a reference, high quantitative accuracy of predicted relative pKa values was seen for structurally more similar substituted sulfanilamides, whereas the result for benzenesulfonamide was found to deviate from experimental value somewhat more. This in turn indicates that testing is to be made if benzenesulfonamide—which indeed may seem as a natural choice—is considered as reference compound for primary benzensulfonamides. We have investigated the effect of the conformer-averaged free energies on the accuracy of pKa predictions, and although some improvement was evident, this procedure does not pay off in general. As a result, we recommend predicting pKa values of benzenesulfonamides using simply the free energies of single reasonable conformers. In addition, the effects of zero-point energy and thermal contributions to the free energy cancel out in the relative pKa predictions using the proton exchange scheme, making the costly Hessian-based calculation of vibrational modes unnecessary for this purpose.

Acknowledgment Computations were performed on resources at the High Perform_ ance Computing Center “HPC Sauletekis” at Faculty of Physics, Vilnius University. Keywords: aqueous acidity  primary benzenesulfonamides  SMD  density functional theory  thermodynamic cycle WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

How to cite this article: K. Aidas, K. Lanevskij, R. Kubilius, L. Juska, D. Petkevicˇius, P. Japertas. J. Comput. Chem. 2015, 36, 2158–2167. DOI: 10.1002/jcc.23998

]

Additional Supporting Information may be found in the online version of this article.

[1] C. T. Supuran, A. Innocenti, A. Mastrolorenzo, A. Scozzafava, Mini Rev. Med. Chem. 2004, 4, 189. [2] T. H. Maren, Annu. Rev. Pharmacol. Toxicol. 1976, 16, 309. [3] F. Carta, A. Scozzafava, C. T. Supuran, Expert Opin. Ther. Pat. 2012, 22, 747. [4] C. T. Supuran, F. Briganti, S. Tilli, W. R. Chegwidden, A. Scozzafava, Bioorg. Med. Chem. 2001, 9, 703. [5] V. Alterio, A. D. Fiore, K. D’Ambrosio, C. T. Supuran, G. D. Simone, Chem. Rev. 2012, 112, 4421. [6] M. Aggarwal, B. Kondeti, R. McKenna, Bioorg. Med. Chem. 2013, 21, 1526. [7] K. Nordfors, J. Haapasalo, M. Korja, A. Niemel€a, J. Laine, A.-K. Parkkila, S. Pastorekova, J. Pastorek, A. Waheed, W. S. Sly, S. Parkkila, H. Haapasalo, BMC Cancer 2010, 10, 148. [8] M.-H. Chien, T.-H. Ying, Y.-H. Hsiegh, C.-H. Lin, C.-H. Shih, L.-H. Wei, S. F. Yang, Oral Oncol. 2012, 48, 417.  _ E. Capkauskait _ J. Gylyte, _ M. Kisonaite, _ S. Tumkevicˇius, D. [9] A. Zubriene, e, Matulis, J. Enzyme Inhib. Med. Chem. 2014, 29, 124.  _ A. Zubriene, _ A. Smirnov, J. Torresan, M. Kisonaite, _ J. [10] E. Capkauskait e, _ J. Gylyte, _ V. Michailoviene, _ V. Jogaite, _ E. Manakova, S. Kazokaite, Grazulis, S. Tumkevicˇius, D. Matulis, Bioorg. Med. Chem. 2013, 21, 6937. [11] D. P. Reynolds, K. Lanevskij, P. Japertas, R. Didziapetris, A. Petrauskas, J. Pharm. Sci. 2009, 98, 4039. [12] K. Lanevskij, P. Japertas, R. Didziapetris, A. Petrauskas, J. Pharm. Sci. 2009, 98, 122. [13] T. Rodgers, M. Rowland, Pharm. Res. 2007, 24, 918. [14] P. H. Bell, R. O. Roblin, Jr., J. Am. Chem. Soc. 1942, 64, 2905. [15] C. Liao, M. C. Nicklaus, J. Chem. Inf. Model. 2009, 49, 2801. [16] J. Manchester, G. Walkup, O. Rivin, Z. You, J. Chem. Inf. Model. 2010, 50, 565. [17] F. Milletti, L. Storchi, L. Goracci, S. Bendels, B. Wagner, M. Kansy, G. Cruciani, Eur. J. Med. Chem. 2010, 45, 4270. [18] J. Ho, M. L. Coote, WIREs Comput. Mol. Sci. 2011, 1, 649. [19] J. Ho, M. L. Coote, Theor. Chem. Acc. 2010, 125, 3. [20] K. S. Alongi, G. C. Shields, In Annual Reports in Computational Chemistry, Vol. 6, Chapter 8; R. A. Wheeler, Ed.; Elsevier B. V.: Amsterdam, 2010; pp. 113–138. [21] D. M. Chipman, J. Phys. Chem. A 2002, 106, 7413. [22] L. A. Curtiss, K. Raghavachari, P. C. Redfern, V. Rassolov, J. A. Pople, J. Chem. Phys. 1998, 109, 7764. [23] J. A. Montgomery, Jr., M. J. Frisch, J. W. Ochterski, G. A. Petersson, J. Chem. Phys. 1999, 110, 2822. [24] A. D. Boese, M. Oren, O. Atasoylu, J. M. L. Martin, M. Kallay, J. Gauss, J. Chem. Phys. 2004, 120, 4129. [25] M. D. Liptak, G. C. Shields, J. Am. Chem. Soc. 2001, 123, 7314. [26] M. W. Palascak, G. C. Shields, J. Phys. Chem. A 2004, 108, 3692. [27] M. D. Tissandier, K. A. Cowen, W. Y. Feng, E. Gundlach, M. H. Cohen, A. D. Earhart, J. V. Coe, T. R. Tuttle, Jr., J. Phys. Chem. A 1998, 102, 7787. [28] C. P. Kelly, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. B 2006, 110, 16066. [29] E. Soriano, S. Cerdan, P. Ballesteros, J. Mol. Struct. (Theochem) 2004, 684, 121. [30] P. G. Seybold, W. C. Kreye, Int. J. Quantum Chem. 2012, 112, 3769. [31] K. C. Gross, P. G. Seybold, C. M. Hadad, Int. J. Quantum Chem. 2002, 90, 445.

[32] F. Eckert, A. Klamt, J. Comput. Chem. 2006, 27, 11. [33] S. Zhang, J. Comput. Chem. 2012, 33, 517. [34] R. Casasnovas, D. Fernandez, J. Ortega-Castro, J. Frau, J. Donoso, F. Mu~ noz, Theor. Chem. Acc. 2011, 130, 1. [35] V. S. Bryantsev, M. S. Diallo, W. A. Goddard, III, J. Phys. Chem. B 2008, 112, 9709. [36] A. M. Toth, M. D. Liptak, D. L. Phillips, G. C. Shields, J. Chem. Phys. 2001, 114, 4595. [37] C. P. Kelly, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. A 2006, 110, 2493. [38] J. Ho, M. L. Coote, J. Chem. Theory Comput. 2009, 5, 295. [39] A. V. Marenich, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. B 2009, 113, 6378. [40] A. V. Marenich, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. B 2009, 113, 4538. [41] A. D. Becke, J. Chem. Phys. 1993, 98, 5648. [42] W. J. Hehre, R. J. Ditchfield, J. A. Pople, J. Chem. Phys. 1972, 56, 2257. [43] P. C. Hariharan, J. A. Pople, Theor. Chim. Acta 1973, 28, 213. [44] Y. Zhao, N. E. Schultz, D. G. Truhlar, J. Chem. Theory Comput. 2006, 2, 364. [45] M. A. Murcko, Theor. Chem. Acc. 1997, 96, 56. [46] P. V. Bharatam, A. A. Gupta, D. Kaur, Tetrahedron 2002, 58, 1759. [47] M. Remko, J. Phys. Chem. A 2003, 107, 720. [48] J. R. B. Gomes, P. Gomes, Tetrahedron 2005, 61, 2705. , J. Kozısek, J. Mol. Struct. 2014, 1059, 124. [49] M. Remko, P. Herich, F. Gregan [50] M. Remko, C.-W. von der Lieth, Bioorg. Med. Chem. 2004, 12, 5395. [51] M. Remko, J. Mol. Struct. (Theochem) 2010, 944, 34. [52] J. J. Klicic´, R. A. Friesner, S.-Y. Liu, W. C. Guida, J. Phys. Chem. A 2002, 106, 1327. [53] F. Jensen, Introduction to Computational Chemistry, 2nd ed.; John Wiley: Chichester, 2007. [54] J. H. Jensen, H. Li, In Encyclopedia of Inorganic Chemistry; R. B. King, Ed.; Wiley: Amsterdam, 2009. [55] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009. [56] R. A. Kendall, T. H. Dunning, R. J. Harrison, J. Chem. Phys. 1992, 96, 6796. [57] A. P. Scott, L. Radom, J. Phys. Chem. 1996, 100, 16502. [58] V. Petrov, V. Petrova, G. V. Girichev, H. Oberhammer, N. I. Giricheva, S. Ivanov, J. Org. Chem. 2006, 71, 2952. [59] M. D. Liptak, K. C. Gross, P. G. Seybold, S. Feldgus, G. C. Shields, J. Am. Chem. Soc. 2002, 124, 6421. [60] V. M. Krishnamurthy, B. R. Bohall, C.-Y. Kim, D. T. Moustakas, D. W. Christianson, G. M. Whitesides, Chem. Asian J. 2007, 2, 94.

Received: 13 March 2015 Revised: 1 June 2015 Accepted: 12 June 2015 Published online on 7 July 2015

Journal of Computational Chemistry 2015, 36, 2158–2167

2167