Stability of hydroxylamine isomers in aqueous solution: Ab initio study using continuum, cluster-continuum and Shells Theory of Solvation

Chemical Physics Letters 518 (2011) 61–64

Contents lists available at SciVerse ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Stability of hydroxylamine isomers in aqueous solution: Ab initio study using continuum, cluster-continuum and Shells Theory of Solvation Guilherme Ferreira de Lima a, Josefredo R. Pliego Jr. b,⇑, Hélio Anderson Duarte a a b

Grupo de Pesquisa em Química Inorgânica Teórica (GPQIT), Departamento de Química – ICEx – UFMG, 31.270-901 Belo Horizonte, MG, Brazil Departamento de Ciências Naturais, Universidade Federal de São João del-Rei, 36301-160 São João del-Rei, MG, Brazil

a r t i c l e

i n f o

Article history: Received 4 October 2011 In final form 1 November 2011 Available online 7 November 2011

a b s t r a c t The relative stability of the NH2OH and NH3O isomers of hydroxylamine in aqueous solution was investigated at ab initio CCSD(T)/TZVPP+diff level for gas phase energy and Shells Theory of Solvation (STS) method for solvation contribution. Our calculations point out that the zwitterionic form (NH3O) is only 3.5 kcal mol1 less stable than the normal isomer (NH2OH). The solvation contribution obtained by the SM8 and SM8AD models are in close agreement with the STS calculations, whereas the PCM (or CPCM) and cluster-continuum methods need small cavities for predicting the correct stability value. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Hydroxylamine is an interesting chemical species with increased reactivity toward organic subtracts due to the called a-effect [1]. The origin of this effect has been widely investigated for anionic nucleophiles whereas for neutral species the problem is less explored [1–6]. In order to understand its solution chemistry, it is very important to know what chemical species is present and their potential role in the reaction. In fact, an important problem in chemistry is the speciation of molecules. In the case of hydroxylamine, there are two possibilities: the normal isomer, NH2OH, and the zwitterionic form, NH3O. The normal form has long been considered as the most stable species in liquid phase. However, recent reports have suggested the zwitterionic form is relatively stable and could be present in a significant proportion in aqueous solution [7–9]. As a consequence, this species has been considered as the active one in reactions involving hydroxylamine [10]. In addition, Mazera et al. have reported experimental [11] and theoretical [12] studies suggesting that both the isomers play an important role in the reaction of hydroxylamine with carboxylic esters in aqueous solution. Although the isomerization from normal to zwitterionic form is a very simple reaction, this equilibrium has been considered difficult to describe through pure continuum solvation models [9]. In fact, an early report by Kirby et al. [9] has predicted an uncertain around 15 kcal mol1 in the solvation free energy difference for the zwitterionic and normal forms when using different solvation models. For example, using PCM/B3LYP/6-311++G(d,p) level of the-

⇑ Corresponding author. E-mail address: [email protected] (J.R. Pliego Jr.). 0009-2614/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2011.11.001

ory, the normal isomer is more stable by 11.2 kcal/mol in free energy. A cluster model with six water molecules plus the continuum has reversed the stability although this kind of calculation does not provide a true free energy and no conclusion can be reached. More computations using Monte Carlo simulation with a classical forcefield and gas phase computations at MP4//B3LYP/6-311++G(d,p) level has predicted the normal isomer more stable for 8.8 kcal mol1, a small difference from pure continuum solvation model [9]. However, using atomic charges obtained from PCM calculations instead of gas phase charges, the solvation free energy obtained from Monte Carlo method predicts the normal isomer becomes less stable for 2 kcal mol1. These differences are very significant and point out how difficult is describing the zwitterionic isomer in solution [9]. The problem was reexamined by Kirby and co-workers [7] using structure–reactivity correlations and those authors have predicted the normal isomer is the most stable species by only 0.9 kcal mol1. Fernandez et al. [8] has also examined the problem and using CBS-QB3 and CPCM computed data to predict pKa of different species, they have concluded the normal isomer is more stable for 2.2 kcal mol1 than the zwitterionic form. The aim of this Letter is to use ab initio methods and to apply three levels of solvation free energy calculations to describe the stability of hydroxylamine isomers. The first approach corresponds to the popular continuum solvation models, like the different implementations of the polarizable continuum model (PCM) [13–18] and the recently developed SM8 [19], SM8AD [20] and SMD [21] methods. The second approach is the cluster-continuum model, a static hybrid discrete/continuum method [22,23]. The third approach is the linear response discrete/continuum Shells Theory of Solvation, a dynamical hybrid discrete/continuum approach recently developed [24,25].

62

G.F. de Lima et al. / Chemical Physics Letters 518 (2011) 61–64

2. Ab initio calculations Geometry optimization and harmonic frequency calculations were carried out at PCM/X3LYP/6-31+G(d) level of theory. For this PCM method [17,15], we have used the IEF formalism [14] and the following atomic radii: N (1.60), O (1.50), H (1.20), scale factor of 1.20 and the dielectric constant of water. Single point energy calculations were done at MP4 and CCSD(T) levels using the Ahlrichs’s Def2-TZVPP basis set [26] extended with diffuse sp functions on N and O atoms with exponents 0.05 and 0.06, respectively. We have named this basis set as TZVPP+diff. All the PCM and CCSD(T) calculations were done with the GAMESS program [27] and the MP4 calculations were done with the FIREFLY program [28]. More solvation free energy computations were done using different PCM cavities. Some computations with the GAMESS have used the above atomic radii with scale factor of 1.20, 1.10 and 1.00 in conjunction with the X3LYP/6-31+G(d) method for electronic density. We have named these calculations as PCM. Because the widespread use of the GAUSSIAN program for predicting chemical reactivity in solution, we have also tested other usual cavities [18] like UAHF, UA0, UAKS and Pauling available in GAUSSIAN 03 [29] using the CPCM/B3LYP/6-31+G(d,p) level of theory. The SMx family of solvation models of Cramer and Truhlar was extensively parametrized for aqueous solution and have presented a very good performance when compared with different continuum models [30–32]. We have used two recent versions, the SM8 [19] and SM8AD [20] for predicting the relative stability of hydroxylamine isomers in aqueous solution. These calculations were done using the M06-2X functional [33] with the 6–31G(d) basis set and CM4M charges with the GAMESSPLUS program [34]. Other model of the SMx family is the SMD, which combines the cavity–dispersion–solvent structure (CDS) terms for non-electrostatic with PCM for continuum electrostatic interaction. For these calculations, we have tested the effect of using the functionals M06-2X and X3LYP, as well as the effect of diffuse functions. All the SMD computations were done with a recent version of the GAMESS program [27]. The cluster-continuum model was also used for predicting the water solvation effect [22]. One and two explicit water molecules were included in the solvation shell of the zwitterionic form of hydroxylamine. This method considers the solvation process as:

Figure 1. Thermodynamic cycle used in the cluster-continuum model.

[24,25]. In this method, the solute and solvent configurations are generated through a molecular dynamics simulation and a number n of solvent molecules are chosen to interact with the solute. The remained solvent molecules are represented through a continuum solvation model. The solvation free energy of solute A is calculated through the equation:

DGsolv ðAÞ ¼

1 < U AS1 > þ < DGsolv ðAS1 Þ  DGsolv ððAÞS1 Þ > 2

where A is the solute and S the solvent. The first term is the clustering free energy in the gas-phase (1 mol/L standard state), the second and third terms are the solvation free energy of the cluster and of the solvent molecule, respectively. The last term is a correction because the different concentration of the solvent molecules in the gas-phase (1 mol/L standard state) and solution phase (55.5 mol/ L). The above equation was derived using chemical potentials [22] and it is equivalent to the thermodynamic cycle presented in the Figure 1. The calculations for the cluster-continuum method were done at MP4/TZVPP+diff//PCM/X3LYP/6-31+G(d) level for gas-phase and PCM/X3LYP/6-31+G(d) method for solvation with the above defined atomic radii and scale factor of 1.00, 1.10 and 1.20.

where the first term is the average solute–solvent interaction involving the solute A and n explicit solvent molecules of the inner solvation shell. The second term is the solvation free energy of the solute plus the inner shell by the dielectric continuum and the third term is similar to second term excluding the solute and adding a cavity to avoid the dielectric continuum fulfill the solute cavity. The <. . .> brackets mean thermal average. Molecular dynamics simulations for each hydroxylamine tautomer were carried out using a cubic box (10  10  10 Å3) with 32 water molecules and periodic boundary conditions. DC-SCC-DFTB quantum method [35–37] implemented in the deMon-nano [38] was used to describe the potential energy surface of each system. Initially we performed 500 ps of molecular dynamics simulation to equilibrate the systems. During the equilibrations the Berendsen thermostat were used to give an average temperature of 300 K. After the equilibration the thermostat was removed and a production run of 500 ps of molecular dynamics simulation were performed to produce the trajectories. The time step during the production run was 0.5 fs and the deviation of the average temperature was estimated to be 12 K. After the molecular dynamics simulation 2500 structures equally spaced were selected with the hydroxylamine and the 14 water molecules nearest to the center of mass of the solute considering both the oxygen and the hydrogen atoms of water molecule. In our previous work we showed that the 14 explicit water molecules is a reasonable number to calculate the solvation free energy [25]. The interaction between the solute and the 14 water molecules were determined using the density functional theory with the PBE exchange/correlation functional and the 6-31+G(d,p) basis set as it is implemented in GAUSSIAN03 package. The interaction energy was estimated as follows: (i) First one single point calculation of the hydroxylamine with the 14 water molecules was performed. (ii) After, the energy of the solute in the same geometry of the first calculation was calculated. (iii) Finally, the energy of the 14 water molecules in the same geometry of the first calculation was calculated; however, without the solute. The long range effects in the solvation free energy were estimated using the PCM method using the HF/6-31+G(d,p) wave function. The cavities were created using the atomic radii internally storage in GAMESS with the scale factor of 2.0. In the solvation of n solvent molecules without the solute, we have added a Ne atom to the center of mass of the hydroxylamine to fulfill the hydroxylamine cavity. All the calculations using the PCM method were performed in the GAMESS package.

3. Molecular dynamics Shells Theory of Solvation calculations

4. Results and discussion

The solvation free energy was also calculated using the Shells Theory of Solvation (STS) in the linear response approximation

The results of all calculations are presented in Table 1. Our calculations show that in the gas-phase, the normal hydroxylamine

NH3 OðgÞ þ n H2 OðsolÞ ! NH3 OðH2 OÞnðsolÞ And the solvation free energy is calculated through the equation:

DGsolv ðAÞ ¼ DGclust ðAðSÞn Þ þ DGsolv ðAðSÞn Þ  nDGsolv ðSÞ  nRT ln½SðsolÞ 

63

G.F. de Lima et al. / Chemical Physics Letters 518 (2011) 61–64

isomer is more stable than the zwitterionic form by 27.6 kcal mol1 (DG), in good agreement with previous reports. We have also found that MP2, MP4 and CCSD(T) calculations predicts close relative energies. Thus, the zwitterionic form should not play any role in the gas-phase reactions of hydroxylamine. On the other hand, in aqueous solution the zwitterionic form is much more solvated. The most sophisticated method used in the Letter is the Shells Theory of Solvation (STS), where full quantum mechanical solute–solvent interaction energies involving 14 water molecules selected from molecular dynamics trajectories were considered. The remained solvation shell was described through the PCM method. In this approach, the zwitterionic isomer is stabilized (DDGsolv) by 24.0 kcal mol1, resulting that the normal isomer is more stable for only 3.5 kcal mol1 than the zwitterion. We have considered this value is the best one calculated in this Letter. When the CPCM methods with the UAHF, UAKS and UA0 cavities are used, the solvation contribution to the free energy favors the zwitterionic form by 12.8–14.5 kcal mol1, resulting that the greater stability of the normal isomer is in the range of 13.1–14.8 kcal mol1. However, using the smallest cavity, the Pauling radii, the stabilization of the zwitterion increases even more, becoming only 4.4 kcal mol1 less stable than the normal isomer. Similar trends can be observed in the PCM calculations with the GAMESS cavities with different scale factors. For example, for scale factor of 1.20, the normal isomer is more stable than the zwitterion for 12.0 kcal mol1. When the scale factor was set to 1.00, the difference decreased to 5.9 kcal mol1, close to the CPCM/Pauling

results. This somewhat arbitrary definition of cavities is undesirable and it is a weakness of continuum solvation models. An optimal cavity should be found for each solvent and be able to treat both cations and anions simultaneously. However, this goal has been very difficult [39,40]. On the other hand, the performance of the SMx family of continuum models of Cramer and Truhlar was extensively tested for ions and neutral species, with a good performance. Thus, we have included solvation calculations using the SM8, SM8AD and SMD methods. The SM8 and SM8AD methods predict an even higher stabilization of the zwitterion than PCM with small cavities, indicating that the normal isomer is more stable than the zwitterion form by only 3.8 and 4.0 kcal mol1, respectively. These values are in close agreement with the STS calculations. In the case of SMD, the performance is worse. We have tested the effect of using different functionals (M06-2X and X3LYP) and the effect of diffuse functions. All the three SMD computations are close each other and have performance similar to PCM with scale factor of 1.10. These differences in continuum models are related to the size of cavity for electrostatic contribution. Indeed, the SM8 method has smaller cavities and it is in better agreement with the STS calculations. The cavities for SMD are greater, resulting in a poor performance of this method. A similar observation on the effect of the cavities for ions was done in a theoretical study of pKa of alcohols using continuum model [41]. The cluster-continuum model with up to two water molecules was also applied for this problem. We have found that just one explicit water molecule produces the lowest solvation free energy for the zwitterion with scale factor of 1.20 and 1.10 in the PCM

Table 1 Thermodynamic data for hydroxylamine isomerization.a NH2OH ? NH3O

CPCM-UAHF/B3LYP/6-31+G(d,p) CPCM-UA0/B3LYP/6-31+G(d,p) CPCM-UAKS/B3LYP/6-31+G(d,p) CPCM-Pauling/B3LYP/6-31+G(d,p) SM8/M06-2X/6–31G(d) SM8AD/M06-2X/6–31G(d) SMD/M06-2X/6–31G(d) SMD/M06-2X/6-31+G(d) SMD/X3LYP/6-31+G(d) PCM (1.20) /X3LYP/6-31+G(d) PCM (1.10) /X3LYP/6-31+G(d) PCM (1.00) /X3LYP/6-31+G(d) Cluster-continumm (PCM 1.20) Cluster-continumm (PCM 1.10) Cluster-continumm (PCM 1.00) STS

MP2b

MP4b

CCSD(T)b

DG g c

26.27 DGsolvd NH2OH

26.05 DGsolvd NH3O

26.16 DDGsolve

27.57 DGsolf

8.96 7.57 8.96 9.57 6.07 5.73 4.60 5.37 5.16 6.16 8.22 10.91 6.16 8.22 10.91 8.93

21.72 21.57 23.43 32.74 29.81 29.28 21.94 24.18 23.45 21.76 26.71 32.55 23.71 27.45 32.15 32.98

12.76 14.00 14.47 23.17 23.75 23.54 17.34 18.81 18.29 15.60 18.49 21.64 17.55 19.23 21.24 24.05

14.81 13.57 13.10 4.40 3.83 4.03 10.23 8.76 9.28 11.97 9.08 5.93 10.02 8.34 6.33 3.52

Literature data PCM/B3LYP/6-311++G(d,p)g MCh MC(PCM)i Structure–reactivity corr.j CPCM-Pauling/CBS-QB3k a b c d e f g h i j k l

DG g c 25.0 27.60 27.60 – 25.5

DDGsolve 13.8 18.76 29.65 – 23.3

Units of kcal/mol. Ab initio calculations with the TZVPP+diff basis set. Gas-phase free energy change for the NH2OH ? NH3O process using CCSD(T)/TZVPP+diff energies and PCM/X3LYP/6-31+G(d) frequencies. Solvation free energy. Solvation free energy difference for the NH2OH ? NH3O process. Solution phase free energy change for the NH2OH ? NH3O process using the gas-phase CCSD(T)/TZVPP+diff energies and PCM/X3LYP/6-31+G(d) frequencies. Ref. [9]. Monte Carlo free energy perturbation calculation and MP4//B3LYP/6-311++G(d,p) method for gas-phase. Ref. [9]. Monte Carlo free energy perturbation calculation with polarized charges and MP4//B3LYP/6-311++G(d,p) method for gas-phase. Ref. [9]. Estimation based on structure–reactivity correlation. Ref. [7]. CPCM computation with gas-phase CBS-QB3 energies. Ref. [8]. Solution phase free energy change for the NH2OH ? NH3O process.

DGsoll 11.2 8.84 2.05 0.90 2.2

64

G.F. de Lima et al. / Chemical Physics Letters 518 (2011) 61–64

computations. These calculations predict that the normal isomer is more stable for 10.0 and 8.3 kcal mol1, respectively. When using the scale factor of 1.00, the relative stability decrease to 6.3 kcal mol1, against 5.9 kcal mol1 found for pure PCM computations. These results indicate that for this cavity, the extra water is not needed and point out that even a discrete/continuum quasi-chemical approach is not able to predict accurately the stability of zwitterionic isomer using the usual scale factor of 1.20 for PCM. Comparing the continuum methods with the STS results, we can note that the popular cavities in GAUSSIAN03 (UAHF, UAKS, UA0) and the default in GAMESS (scale factor of 1.20) have poor performance, predicting the normal isomer much more stable than the zwitterion. The poor performance of PCM methods for this problem can be attributed to the highly polar zwitterionic structure, which has similarity with ionic species. Many studies have pointed out that accurate calculations of the solvation free energy of ions is hard for continuum solvation models [40,42]. Therefore, zwitterionic species are also a challenge for continuum models and only small cavities are able to predict the real stability of these species in aqueous solution. As a consequence, any theoretical study involving hydroxylamine must take care of the solvation model used. Otherwise, the calculations could lead to wrong prediction of the zwitterion reactivity. 5. Conclusions The present high level theoretical calculations predicts the NH3O isomer is enough stable to be kinetically relevant and should play a key role in many reactions of hydroxylamine in aqueous solution. We have also found the SM8 and SM8AD models are in close agreement with the STS calculations, the SMD method has a poor performance and PCM (or CPCM) method needs small cavities to be reliable for this system. Acknowledgments The authors thank the Brazilian Research Council (CNPq) and the Fundação de Amparo a Pesquisa do Estado de Minas Gerais (FAPEMIG) for the support. References [1] Y.S. Simanenko, A.F. Popov, T.M. Prokopeva, E.A. Karpichev, V.A. Savelova, I.P. Suprun, C.A. Bunton, Russ. J. Org. Chem. 38 (2002) 1341.

[2] S.M. Villano, N. Eyet, W.C. Lineberger, V.M. Bierbaum, J. Am. Chem. Soc. 131 (2009) 8227. [3] A.M. McAnoy, M.R.L. Paine, S.J. Blanksby, Org. Biomol. Chem. 6 (2008) 2316. [4] Y. Ren, H. Yamataka, Chem. Eur. J. 13 (2007) 677. [5] Y. Ren, H. Yamataka, J. Org. Chem. 72 (2007) 5660. [6] E. Buncel, I. Um, Tetrahedron 60 (2004) 7801. [7] A.J. Kirby et al., Chem. Commun. 46 (2010) 1302. [8] M.I. Fernandez, M. Canle, M.V. Garcia, J.A. Santaballa, Chem. Phys. Lett. 490 (2010) 159. [9] A.J. Kirby, J.E. Davies, T.A.S. Brandão, P.F. Silva, W.R. Rocha, F. Nome, J. Am. Chem. Soc. 128 (2006) 12374. [10] A.J. Kirby, D.W. Tondo, M. Medeiros, B.S. Souza, J.P. Priebe, M.F. Lima, F. Nome, J. Am. Chem. Soc. 131 (2009) 2023. [11] D.J. Mazera, J.C. Gesser, J.R. Pliego Jr., ARKIVOC (2007) 199. [12] D.J. Mazera, J.R. Pliego Jr, G.I. Almerindo, J.C. Gesser, J. Braz. Chem. Soc. 22 (2011) 2165. [13] B. Mennucci, E. Cances, J. Tomasi, J. Phys. Chem. B 101 (1997) 10506. [14] E. Cances, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3032. [15] M. Cossi, V. Barone, R. Cammi, J. Tomasi, Chem. Phys. Lett. 255 (1996) 327. [16] S. Miertus, J. Tomasi, Chem. Phys. 65 (1982) 239. [17] S. Miertus, E. Scrocco, J. Tomasi, Chem. Phys. 55 (1981) 117. [18] V. Barone, M. Cossi, J. Tomasi, J. Chem. Phys. 107 (1997) 3210. [19] A.V. Marenich, R.M. Olson, C.P. Kelly, C.J. Cramer, D.G. Truhlar, J. Chem. Theory Comput. 3 (2007) 2011. [20] A.V. Marenich, C.J. Cramer, D.G. Truhlar, J. Chem. Theory Comput. 5 (2009) 2447. [21] A.V. Marenich, C.J. Cramer, D.G. Truhlar, J. Phys. Chem. B 113 (2009) 6378. [22] J.R. Pliego Jr, J.M. Riveros, J. Phys. Chem. A 105 (2001) 7241. [23] E. Westphal, J.R. Pliego Jr, J. Chem. Phys. 123 (2005) 074508. [24] J.R. Pliego Jr, Theor. Chem. Acc. 128 (2011) 275. [25] G.F. de Lima, H.A. Duarte, J.R. Pliego Jr, J. Phys. Chem. B 114 (2010) 15941. [26] F. Weigend, R. Ahlrichs, Phys. Chem. Chem. Phys. 7 (2005) 3297. [27] M.W. Schmidt et al., J. Comput. Chem. 14 (1993) 1347. [28] A.A. Granovsky, PC GAMESS/FIREFLY version 7.1.F. . [29] M.J. Frisch et al., GAUSSIAN 03, Revision B.04, GAUSSIAN, Inc., Pittsburgh, PA, 2003. [30] C.J. Cramer, D.G. Truhlar, Acc. Chem. Res. 41 (2008) 760. [31] A. Klamt, B. Mennucci, J. Tomasi, V. Barone, C. Curutchet, M. Orozco, F.J. Luque, Acc. Chem. Res. 42 (2009) 489. [32] C.J. Cramer, D.G. Truhlar, Acc. Chem. Res. 42 (2009) 493. [33] Y. Zhao, D.G. Truhlar, Theor. Chem. Acc. 120 (2008) 215. [34] M. Higashi, et al., GAMESSPLUS – version 2010-2, University of Minnesota, Minneapolis, 2010, based on the General Atomic and Molecular Electronic Structure System (GAMESS) as described in Ref. [27]. [35] A.F. Oliveira, G. Seifert, T. Heine, H.A. Duarte, J. Braz. Chem. Soc. 20 (2009) 1193. [36] M. Elstner et al., Phys. Rev. B 58 (1998) 7260. [37] L. Zhechkov, T. Heine, S. Patchkovskii, G. Seifert, H.A. Duarte, J. Chem. Theory Comput. 1 (2005) 841. [38] T. Heine, et al., deMon, deMon-nano edn. deMon-Sofware, Mexico, Bremen, 2009. [39] J.R. Pliego Jr, J.M. Riveros, Chem. Phys. Lett. 355 (2002) 543. [40] C. Curutchet, A. Bidon-Chanal, I. Soteras, M. Orozco, F.J. Luque, J. Phys. Chem. B 109 (2005) 3565. [41] I. Tunon, E. Silla, J. Pascual-Ahuir, J. Am. Chem. Soc. 115 (1993) 2226. [42] J.M. Ho, M.L. Coote, Theor. Chem. Acc. 125 (2010) 3.