- Email: [email protected]
Computational study of anion solvation in nitrobenzene
Chemical Physics Letters 436 (2007) 362–367 www.elsevier.com/locate/cplett
Computational study of anion solvation in nitrobenzene Elvis S. Bo¨es a, Jones de Andrade a, Hubert Stassen a
a,*
, Paulo F.B. Gonc¸alves
b
Grupo de Quı´mica Teo´rica, Instituto de Quı´mica, Universidade Federal do Rio Grande do Sul, Av. Bento Gonc¸alves 9500, 91540-000 Porto Alegre-RS, Brazil b Centro Universita´rio La Salle, Av. Victor Barreto, 2288, 92010-000 Canoas-RS, Brazil Received 18 December 2006 Available online 27 January 2007
Abstract The solvation of anions in nitrobenzene was investigated by the polarizable continuum model (PCM). The PCM was parameterized from structural information obtained by molecular dynamics simulations (MD) of anionic solutions in nitrobenzene. The parameterization was performed against experimental free energies of solvation for 22 anions containing H, C, N, O, F, S, Cl, Br, Se, and I atoms. The calculated Gibbs free energies of solvation present a mean absolute deviation from the experimental data of 2.4 kcal/mol. Hartree–Fock and DFT computations produce equivalent results. Ó 2007 Elsevier B.V. All rights reserved.
1. Introduction Solvation studies of ionic species furnish important details for physical and chemical processes in catalytic, biological, and pharmacological systems [1]. For example, the efficiencies of many important drugs acting as ionic species depend on the ion’s affinities with water (hydrophilicity) and apolar organic phases (lipophilicity) [2]. A quantitative measure for this affinity is provided by the Gibbs free energy of solvation for the species and the Gibbs energy of transfer between the two considered phases. The determination of the Gibbs free energy of transfer of ions between two immiscible solutions is usually performed by experimental electrochemical techniques [3]. Within water/organic phase transfer electrochemical studies, nitrobenzene is one of the most employed solvents due to its low miscibility with water and its dielectric constant ( = 34.82 [5]) enabling dissolvation of electrolytes [4]. Motivated by the relevance of the nitrobenzene solvent in experimental studies on the partition of ions between aqueous and organic phases, we decided to investigate theoret*
Corresponding author. Fax: +55 51 3308 7304. E-mail addresses: [email protected] (E.S. Bo¨es), [email protected] (J.de Andrade), [email protected] (H. Stassen), [email protected] (P.F.B. Gonc¸alves). 0009-2614/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.01.063
ically the solvation of anions in nitrobenzene Only few theoretical studies of the solvation in nitrobenzene are reported in the literature. We cite MD of pure 2-nitrophenyl octyl ether, nitrobenzene, and solutions of K+ and Cl ions in these two solvents [6] and a quantum mechanical approach to the transfer energy of ions in water/nitrobenzene [7]. In this article, we focus on the solvation of anions in nitrobenzene. We propose the use of MD on anion solutions in nitrobenzene to define the solute/solvent interface in these solutions. The PCM [8,9] has been utilized to calculate the Gibbs free energy providing details for the thermodynamics of the solvation process in this solvent. 2. Compilation of experimental data To obtain Gibbs energies of solvation for the anions in nitrobenzene, we combined experimental Gibbs free energies of hydration with the transfer energies for the anions between water and nitrobenzene. Experimental hydration Gibbs energies for the anions were collected from Pliego and Riveros [10], Abraham and Liszi [11], and Marcus [4]. These experimental hydration data were corrected, when necessary, in agreement with the recently established standard value for the Gibbs free energy of solvation for the proton in water, DGsolv ðHþ Þ ¼ 264:0 kcal=mol [12],
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
corresponding to the process of transferring the proton from ideal gas at 1 atm to ideal diluted solution at 1 mol/ l. We adopted a single consistent standard state for defining the Gibbs energy of solvation for all anions according to Ben-Naim, corresponding to the process of transfer of 1 mole of solute per liter of ideal gas to 1 mole/liter of ideal solution [13]. In these standard states, the Gibbs free energy of hydration for the proton is DGsolv ðHþ Þ ¼ 265:9 kcal=mol [10,14]. Transfer free energies between water and nitrobenzene were taken from Refs. [2,15–17]. These data, combined with the Gibbs free energies of hydration, provide the Gibbs free energies of solvation in nitrobenzene. We have collected data for 22 anions dissolved in nitrobenzene. 3. Computational methods The molecular geometries of the 22 anions were optimized in the gas phase, at Hartree–Fock (HF) RHF/631+G(d,p) and Density Functional Theory (DFT) B3LYP/6-31+G(d,p) levels with the GAMESS package [18]. These optimized structures were utilized in the calculations of free energies in nitrobenzene from the PCM. It is wellknown that the use of optimized gas phase structures in condensed phase simulations does not introduce significant errors for solutes with restricted conformational changes [19,20]. In this study, we have used the GEPOL algorithm as implemented in the GAMESS package to create the cavities for the solutes [21]. The GEPOL algorithm generates a cavity for the solute from the positions of the atoms and their radii. Starting from a set of original atomic radii Roi , solvent dependent radii Rcav have been established by adopti ing a factor f to Roi , ¼ f Roi : Rcav i
ð1Þ
The algorithm creates a cavity corresponding to the molecular shape by interlocking spheres of radii Rcav i . Addi˚ were tional smaller spheres with minimum radius of 0.2 A added to smooth the cavity’s surface. The surface of each sphere was divided into 60 discrete elements of area (tesserae). In this work, the original atomic radii Roi from Emsley, modified by the Pisa group [22] were employed as implemented in the GAMESS package. These atomic radii are listed in Table 1. The determination of the suitable scaling factor f from Eq. (1) is one of our main objectives. Therefore, MD of a single anion Cl, Br, I, and CN in nitrobenzene were performed with the GROMACS 3.3 software [23]. Halide anions were described by the OPLS force field [24], CN by the parameters from Ref. [25], and liquid nitrobenzene
363
Table 2 Properties of the liquid solvent nitrobenzene at 298.15 K Property Dielectric constant, [5] ˚ ) [11] Molecular radius (A Density (g cm3) [5]
34.82 2.773 1.198
by the OPLS-AA model [26]. The MD were carried out in the NVT ensemble with 255 nitrobenzene molecules and a single anion arranged in a cubic box with dimensions corresponding to the experimental density (see Table 2) of the solvent at simulation temperature (T = 298 K). Periodic boundary conditions and the minimum image convention were applied. Long range electrostatic interactions were taken into account by the Particle Mesh Ewald method (PME) [27]. Equations of motion were integrated by the leap-frog algorithm [28] with a time step of 1 fs. The bond length of CN anion was fixed by the SHAKE algorithm [30]. The temperature was maintained by the weak coupling Berendsen thermostat with a time constant of 0.1 ps [29]. Data were averaged from production runs of 500 ps preceeded by a stabilization period of 1 ns. The atom–atom pair distribution functions were computed from configurations saved in intervals of 50 time steps. The radial distribution functions (RDFs) for X–H atom pairs, where X represents one of the anion atoms and H the hydrogen atom of the nitrobenzene molecules, have been analyzed in more detail. The distance corresponding to the maximum position in the first peak of these RDFs was utilized to localize the average position of the hydrogen atoms responsible for weak hydrogen bonding between the first shell solvent molecules and the solutes. This structural parameter for the solvent structure around the solutes was used to define the scaling factor from Eq. (1). This methodology of defining the solute–solvent interface in the PCM from MD has been applied before with success in PCM studies of anionic solvation in acetonitrile and DMF [20]. In the PCM formalism, the process of solvation is described by inserting the solute inside a molecular cavity into the dielectric medium representing the solvent. As in our earlier study [20], the free energy of solvation DGsolv is partitioned into the electrostatic energy DGele, the cavitation energy DGcav, and the van der Waals energy DGvdw: DGsolv ¼ DGele þ DGcav þ DGvdW :
ð2Þ
More details about the formalism have been reviewed by Refs. [8,9]. The electrostatic term DGele was calculated from the Integral Equation Formalism version of the PCM, as implemented in the GAMESS package [31,32] at the same
Table 1 ˚ ) for the definition of the solute cavities Atomic radii (in A Element
H
C
N
O
F
S
Cl
Se
Br
I
Radius
1.20
1.70
1.60
1.50
1.35
1.85
1.81
1.90
1.95
2.15
364
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
levels of the gas-phase optimizations, the HF RHF/631+G(d,p) and DFT B3LYP/6-31+G(d,p). Charge density escaping from the cavity boundary is corrected by the charge renormalization procedure of Mennucci and Tomasi [33]. The physical parameters for the nitrobenzene solvent within PCM are summarized in Table 2. The cavitation contribution DGcav was computed from scaled particle theory adapted to molecular cavities as described in Ref. [34]. To computed the van der Waals repulsion–dispersion energy DGvdW, we assume a linear relationship between this term and the solvent accessible surface area of the solute: DGvdW ¼
N X
ni A i :
ð3Þ
i¼1
Ai is the solvent accessible surface area of the ith solute atom with van der Waals parameter ni optimized for this atom [8]. The optimization procedure for ni involve a multiple linear regression procedure minimizing the difference between experimental solvation free energies and the sum of the computed electrostatic and cavitational contributions DGsolv . Thus, the PCM is parameterized from two steps: (i) defining from MD the scaling factors for atomic radii in the cavity creation, and (ii) optimizing the van der Waals parameters for the computation of DGvdW.
Fig. 2. Radial distribution functions for N–H pairs in nitrobenzene solutions.
4. Results and discussion From the MD of the anionic (CN, Cl, Br or I) solutions in nitrobenzene, we computed the RDFs representing solute–solvent atom pair correlations. The most interesting RDFs are those representing the X–H pairs. The H atom of the nitrobenzene molecule possesses a positive charge of q = 0.115e in OPLS-AA [26]. Thus, one might expect sharp first peaks in these RDFs that we have chosen to define the solute–solvent interface in the continuum model. In Figs. 1–5, we have depicted the RDFs for the X–H pairs. Each figure contains three functions corresponding to ortho, meta, and, para hydrogens in the nitro-
Fig. 3. Radial distribution functions for Cl–H pairs in nitrobenzene solutions.
Fig. 4. Radial distribution functions for Br–H pairs in nitrobenzene solutions.
Fig. 1. Radial distribution functions for C–H pairs in nitrobenzene solutions.
benzene ring. For all the X–Hortho pairs, the RDFs present diffuse profiles due to the repulsion between negative charges of the solute and oxygen sites in the nitro group
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
365
Table 3 ˚ 2) optimized for the solvation Van der Waals parameters ni (in kcal/mol A of anions in nitrobenzene
Fig. 5. Radial distribution functions for I–H pairs in nitrobenzene solutions.
of nitrobenzene. On the other hand, RDFs for X–Hpara and XHmeta pairs exhibit well-defined sharp peaks. We observe that the distance corresponding to the maximum in the first peak of the RDFs is slightly shorter for meta hydrogens. Consequently, we utilized the RDFs for meta hydrogens as our standard to establish the solute–solvent interface in the nitrobenzene solutions. For CN, we calculated RDFs for the carbon and the nitrogen sites illustrated in Figs. 1 and 2. The RDF for ˚ the C–Hmeta pair possesses a maximum at rmax = 2.778 A and, for the pair N–Hmeta, a first maximum is located at ˚ . The atomic Emsley radius for carbon is rmax = 2.685 A ˚ , whereas the nitrogen atomic radius is 1.60 A ˚. 1.70 A Therefore, the ratios rmax =Roi are 1.63 and 1.68, respectively. As in our earlier study [20], these ratios are used to define the scaling factor f from Eq. (1). The RDF for the Cl–Hmeta (3) exhibits a first peak with maximum at ˚ . With the chlorine atomic radius, we obtain rmax = 3.086 A the ratio 1.70 defining f. Fig. 4 presents the RDFs for the Br solution. The first maximum in the RDF involving ˚ producing the meta hydrogen is located at rmax = 3.148 A f = 1.61 if the atomic radius for the bromine from Table ˚ ) is used. Finally, in the RDF for the I solution 1 (1.95 A ˚ . The (Fig. 5), we localize a first maximum at rmax = 3.376 A ratio between this distance and the atomic radius of the ˚ ) is f = 1.57. Averaging the five scaling iodine atom (2.15 A factors, we computed an averaged scaling factor of f = 1.64 for anionic solutions in nitrobenzene. In the PCM calculations, we adopted this scaling factor to all the anions investigated in the present study. Next, using the average scaling factor, we computed the electrostatic and cavitation contributions to the Gibbs free energy of solvation for the anions in nitrobenzene. Comparing the sum of the calculated terms with the experimental energy of solvation, we optimized the van der Waals atomic parameters ni of Eq. (3) using the criterium of minimizing differences between experimental and theoretical free energies of solvation. Details of this procedure are described in Ref. [20]. The optimized parameters ni are
Atom
ni (RHF)
ni (DFT)
H C N O F S Cl Se Br I
0.5010 0.0740 0.2499 0.1163 0.3151 0.1625 0.2861 0.6348 0.3780 0.2445
0.5285 0.0710 0.2625 0.1177 0.3282 0.1547 0.2873 0.6050 0.3804 0.2445
shown in the Table 3 for both, the HF and DFT levels of theory. Employing the optimized ni, we computed the van der Waals term of the energy of solvation for the anions according to Eq. (3). Afterwards, DGsolv according to Eq. (2) was computed for the set of 22 anions in nitrobenzene. Tables 4 and 5 summarize these results for the HF and DFT levels, respectively, comparing with the experimental values. The mean absolute deviation between theoretical and experimental data is 2.4 kcal/mol for both, HF and DFT calculations. The mean signed error is 0.6 kcal/mol and 0.7 kcal/mol using HF and DFT levels, respectively. The standard deviation between calculated and experimental values is 3.4 kcal/mol at the HF level and 3.5 kcal/mol for the solvation energies calculated at DFT level Also, it becomes evident from Tables 4 and 5 that the calculations performed at HF or DFT levels produce similar values for the energies of solvation. For small anions, the difference in the electrostatic contribution calculated employing the HF or DFT levels is approximately 0.5 kcal/mol. For most anions in the set, we achieved a very good agreement between calculated and experimental free energies of solvation. No systematic deviations are observed in our results. Larger deviations from experimental data have been computed for HCO 2 , CH3 CO2 , and picrate anions. For the HCO2 anion, we found a deviation of 9.4 kcal/mol at HF level and 9.9 kcal/mol at DFT level. In the case of the CH3 CO 2 anion, the error is 6.5 kcal/mol at HF level and 6.6 kcal/mol at DFT level. For the Pic anion, the observed deviation is negative, 5.4 kcal/mol when the calculations are performed at HF level and 5.5 kcal/mol at DFT level. We believe that the origin of these deviations is the same for the three cases and is related to DGvdW. In our methodology, the van der Waals parameters ni are statistically adjusted to produce the best agreement between the calculated and the experimental free energies of solvation for all the solutes. This method attributes equal statistical weights to all the solutes of the set. With eleven anions in the set containing H, C, and O atoms, the optimization by multivariate fitting favors the minimization of the errors for most of the anions of the subset containing H, C, and O atoms. However, eight of these eleven anions possess larger surfaces and, therefore, the ni
366
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
Table 4 Electrostatic, DGele, cavitational, DGcav, and van der Waals, DGvdW, contributions to the calculated solvation free energies, DGsolv ðcalcÞ, for the anions in nitrobenzene at HF level Anion
DGele
DGcav
DGvdW
DGsolv ðcalcÞ
DGsolv ðexpÞ
Error
54.2 50.3 46.1 56.4 53.8 46.4 46.9 51.8 48.5 48.1 47.7 56.3 54.6 50.1 54.1 49.5 53.5 50.9 44.3 43.7 47.9 38.6
3.9 4.4 5.1 5.7 4.6 6.7 6.9 6.6 7.7 8.1 8.3 6.1 8.7 9.5 11.1 13.7 14.8 12.6 15.2 15.2 15.0 20.2
11.8 18.1 14.2 13.8 7.1 11.6 29.0 8.1 11.7 15.7 10.5 9.3 17.7 18.0 26.2 34.8 28.2 26.5 26.9 26.8 34.9 27.5
62.1 64.0 55.2 64.5 56.3 51.3 69.0 53.3 52.5 55.7 49.9 59.5 63.6 58.6 69.2 70.6 66.9 64.8 56.0 55.3 67.8 45.9
63.6 61.8 55.2 64.3 60.5 51.3 69.0 55.4 49.2 59.9 49.4 68.9 70.1 58.6 69.5 68.9 66.2 66.4 59.1 52.8 65.4 40.5
+1.5 2.2 0.0 0.2 +4.2 0.0 0.0 +2.1 3.3 +4.2 0.5 +9.4 +6.5 0.0 +0.3 1.7 0.7 +1.6 +3.1 2.5 2.4 5.4
Cl Br I N 3 CN SCN SeCN NO 3 ClO 3 BrO 3 ClO 4 HCO 2 CH3 CO 2 F2 CHCO 2 CH3 CH2 CO 2 CH3 CH2 CH2 CO 2 C6 H5 CO2 C6H5O 3-NO2C6H4O 4-NO2C6H4O 2-CH3C6H4O Pic Mean absolute error Mean error RMS error SD
2.4 0.6 3.3 3.4
DGsolv ðexpÞ are experimental data. Units are in kcal/mol.
Table 5 Same as Table 4, but for the DFT computations Anion
DGele
DGcav
DGvdW
DGsolv ðcalcÞ
DGsolv ðexpÞ
Error
Cl Br I N 3 CN SCN SeCN NO 3 ClO 3 BrO 3 ClO 4 HCO 2 CH3 CO 2 F2 CHCO 2 CH3 CH2 CO 2 CH3 CH2 CH2 CO 2 C6 H5 CO2 C6H5O 3-NO2C6H4O 4-NO2C6H4O 2-CH3C6H4O Pic
54.1 50.2 46.1 55.5 53.5 46.4 45.3 51.5 48.4 47.9 47.5 55.6 53.7 49.3 52.9 47.5 52.2 49.7 43.2 42.9 46.6 37.8
3.9 4.4 5.1 5.7 4.7 6.7 7.1 6.7 7.7 8.1 8.5 6.2 8.8 9.6 11.3 13.9 14.9 12.7 15.4 15.4 15.2 20.6
11.8 18.2 14.2 14.7 7.3 11.6 30.8 8.5 11.8 15.9 10.9 9.6 18.6 18.9 27.6 36.9 29.4 27.8 28.1 28.0 36.7 28.8
62.0 64.0 55.2 64.5 56.1 51.3 69.0 53.3 52.5 55.7 49.9 59.0 63.5 58.6 69.2 70.5 66.7 64.8 55.9 55.5 68.1 46.0
63.6 61.8 55.2 64.3 60.5 51.3 69.0 55.4 49.2 59.9 49.4 68.9 70.1 58.6 69.5 68.9 66.2 66.4 59.1 52.8 65.4 40.5
+1.6 2.2 0.0 0.2 +4.4 0.0 0.0 +2.1 3.3 +4.2 0.5 +9.9 +6.6 0.0 +0.3 1.6 0.5 +1.6 +3.2 2.7 2.7 5.5
Mean absolute error Mean error RMS error SD
2.4 0.7 3.4 3.5
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
parameters were optimized taking into account these large surfaces. Consequently, the obtained parameters describe well the large solutes such as CH3 CH2 CO 2 , CH3 CH2 CH2 CO 2 , and C6 H5 CO2 , but might fail for the smaller solutes HCO 2 and CH3 CO2 . With the ni related to the solvent accessible surface area by Eq. (3), we expect an underestimation of the van der Waals contribution for the small solutes. In the case of the rather large Pic anion, the seven oxygen atoms are exposed to the solvent in external positions of this anion. Thus, the large surface areas of these atoms is responsible for the overestimation of the van der Waals energy. We notice that for most of the anions of the set, DGele is approximately 50 kcal/mol. For smaller and simpler anions, from Cl to ClO 4 in Tables 4 and 5, the electrostatic contribution corresponds to 85–95% of the total energy of solvation. For the larger anions, for example those derived from organic acids, the importance of the electrostatic term is lower and the steric terms (DGcav and DGvdW) become more significant for the solvation free energy. For these organic ions, the van der Waals interactions with the solvent molecules play a very important role in the solvation process in nitrobenzene. 5. Conclusion In the present article, we presented a solvation study of anions in nitrobenzene. We have focused on the calculation of the Gibbs free energy of solvation of 22 representative anions in nitrobenzene from an appropriate parameterization of the PCM for calculations at HF and DFT levels. Relevant structural characteristics extracted from MD of a few anionic solutions in nitrobenzene were utilized to establish the suitable parameter of scaling for the atomic radii in the definition of the solute–solvent interface for the calculations with the continuum model. The scaling factor established for solute cavities in nitrobenzene, f = 1.64, is in agreement with the fact that polar aprotic solvents, such as nitrobenzene, are weak hydrogen bond donors towards anionic solutes [35]. The obtained scaling factor possesses physical sense if one compares with other solventes. For example, the study of anion solvation in N,N-dimethylformamide has been described using a scaling factor f = 1.52 [20]. We have parameterized the PCM and calculated the Gibbs free energies of solvation for 22 anions containing up to 18 atoms. Quite general, the theoretical data agree well with experimental findings corroborating that the parameterized PCM model is able to treat anionic solvation in nitrobenzene. Comparing the parameterization for HF and DFT computations, we conclude that both levels of theory produce qualitatively equivalent solvation free energies. Furthermore, the partitioning of the total energy of solvation into the three contributions, cavitation, electrostatic, and van der Waals, permits better understanding about the role of each of these terms and the nature of the interactions between the solvent nitrobenzene and the anio-
367
nic anion solutes. We believe that the parameterized model can be employed in further studies on the solvation of anionic drugs or other anions with biological relevance. Acknowledgements E.S.B. acknowledges a scholarship from the Brazilian agency CAPES. J.d.A. acknowledges a scholarship from the Brazilian agency CNPq. H.S. acknowledges a grant from Brazilian agency CNPq (309211/2003-4. References [1] A.G. Volkov, Liquid Interfaces in Chemical, Biological and Pharmaceutical Applications, Marcel Dekker, New York, 2001. [2] R. Gulaboski et al., J. Phys. Chem. B 108 (2004) 4565. [3] F. Scholz, Sˇ. Komorsky-Lovric´, M. Lovric´, Electrochem. Commun. 2 (2000) 112. [4] Y. Marcus, Ion Properties, Marcel Dekker, New York, 1997. [5] J.A. Riddick, W.B. Bunger, Organic Solvents Physical Properties and Methods of Purification, John Wiley & Sons, New York, 1986. [6] M. Jorge, R. Gulaboski, C.M. Pereira, M.N.D.S. Cordeiro, J. Phys. Chem. B 110 (2006) 12530. [7] T. Osakai, A. Ogata, K. Ebina, J. Phys. Chem. B 101 (1997) 8341. [8] F.J. Luque et al., Phys. Chem. Chem. Phys. 5 (2003) 3827. [9] J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev. 105 (2005) 2999. [10] J.R. Pliego Jr., J.M. Riveros, Phys. Chem. Chem. Phys. 4 (2002) 1622. [11] M.H. Abraham, J. Liszi, J. Chem. Soc. Faraday Trans. 74 (1978) 1604. [12] M.D. Tissandier et al., J. Phys. Chem A 102 (1998) 7787. [13] A. Ben-Naim, J. Phys. Chem. 82 (1978) 792. [14] D.M. Camaioni, C.A. Schwerdtfeger, J. Phys. Chem. A 109 (2005) 10795. [15] R. Gulaboski, K. Riedl, F. Scholz, Phys. Chem. Chem. Phys. 5 (2003) 1284. [16] Sˇ. Komorsky-Lovric´, K. Riedl, R. Gulaboski, V. Mircˇeski, F. Scholz, Langmuir 18 (2002) 8000. [17] Sˇ. Komorsky-Lovric´, K. Riedl, R. Gulaboski, V. Mircˇeski, F. Scholz, Langmuir 19 (2003) 3090. [18] M.W. Schmidt et al., J. Comput. Chem. 14 (1993) 1347. [19] C. Curutchet, A. Bidon-Chanal, I. Soteras, M. Orozco, F.J. Luque, J. Phys. Chem. B 109 (2005) 3565. [20] E.S. Bo¨es, P.R. Livotto, H. Stassen, Chem. Phys. 331 (2006) 142. [21] J.L. Pascual-Ahuir, E. Silla, J. Tomasi, R. Bonaccorsi, J. Comput. Chem. 8 (1987) 778. [22] J. Emsley, The Elements, Clarendon Press, Oxford, 1991. [23] E. Lindahl, B. Hess, D. van der Spoel, J. Mol. Mod. 7 (2001) 306. [24] N.A. McDonald, E.M. Duffy, W.L. Jorgensen, C.J. Swenson, J. Am. Chem. Soc. 120 (1998) 5104. [25] M. Ferrario, I.R. McDonald, M.C.R. Symons, Mol. Phys. 77 (1992) 617. [26] M.L.P. Price, D. Ostrovsky, W.L. Jorgensen, J. Comput. Chem. 22 (2001) 1340. [27] U. Essman, L. Perela, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577. [28] R.W. Hockney, S.P. Goel, J. Comput. Phys. 14 (1974) 148. [29] H.J.C. Berendsen, J.P.M. Postma, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [30] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [31] E. Cance`s, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3032. [32] E. Cance`s, B. Mennucci, J. Math. Chem. 23 (1998) 309. [33] B. Mennucci, J. Tomasi, J. Chem. Phys. 106 (1996) 5151. [34] J. Langlet, P. Claverie, J. Caillet, A. Pullman, J. Phys. Chem. 92 (1988) 1617. [35] A.J. Parker, Chem. Rev. 69 (1969) 1.
Computational study of anion solvation in nitrobenzene Elvis S. Bo¨es a, Jones de Andrade a, Hubert Stassen a
a,*
, Paulo F.B. Gonc¸alves
b
Grupo de Quı´mica Teo´rica, Instituto de Quı´mica, Universidade Federal do Rio Grande do Sul, Av. Bento Gonc¸alves 9500, 91540-000 Porto Alegre-RS, Brazil b Centro Universita´rio La Salle, Av. Victor Barreto, 2288, 92010-000 Canoas-RS, Brazil Received 18 December 2006 Available online 27 January 2007
Abstract The solvation of anions in nitrobenzene was investigated by the polarizable continuum model (PCM). The PCM was parameterized from structural information obtained by molecular dynamics simulations (MD) of anionic solutions in nitrobenzene. The parameterization was performed against experimental free energies of solvation for 22 anions containing H, C, N, O, F, S, Cl, Br, Se, and I atoms. The calculated Gibbs free energies of solvation present a mean absolute deviation from the experimental data of 2.4 kcal/mol. Hartree–Fock and DFT computations produce equivalent results. Ó 2007 Elsevier B.V. All rights reserved.
1. Introduction Solvation studies of ionic species furnish important details for physical and chemical processes in catalytic, biological, and pharmacological systems [1]. For example, the efficiencies of many important drugs acting as ionic species depend on the ion’s affinities with water (hydrophilicity) and apolar organic phases (lipophilicity) [2]. A quantitative measure for this affinity is provided by the Gibbs free energy of solvation for the species and the Gibbs energy of transfer between the two considered phases. The determination of the Gibbs free energy of transfer of ions between two immiscible solutions is usually performed by experimental electrochemical techniques [3]. Within water/organic phase transfer electrochemical studies, nitrobenzene is one of the most employed solvents due to its low miscibility with water and its dielectric constant ( = 34.82 [5]) enabling dissolvation of electrolytes [4]. Motivated by the relevance of the nitrobenzene solvent in experimental studies on the partition of ions between aqueous and organic phases, we decided to investigate theoret*
Corresponding author. Fax: +55 51 3308 7304. E-mail addresses: [email protected] (E.S. Bo¨es), [email protected] (J.de Andrade), [email protected] (H. Stassen), [email protected] (P.F.B. Gonc¸alves). 0009-2614/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.01.063
ically the solvation of anions in nitrobenzene Only few theoretical studies of the solvation in nitrobenzene are reported in the literature. We cite MD of pure 2-nitrophenyl octyl ether, nitrobenzene, and solutions of K+ and Cl ions in these two solvents [6] and a quantum mechanical approach to the transfer energy of ions in water/nitrobenzene [7]. In this article, we focus on the solvation of anions in nitrobenzene. We propose the use of MD on anion solutions in nitrobenzene to define the solute/solvent interface in these solutions. The PCM [8,9] has been utilized to calculate the Gibbs free energy providing details for the thermodynamics of the solvation process in this solvent. 2. Compilation of experimental data To obtain Gibbs energies of solvation for the anions in nitrobenzene, we combined experimental Gibbs free energies of hydration with the transfer energies for the anions between water and nitrobenzene. Experimental hydration Gibbs energies for the anions were collected from Pliego and Riveros [10], Abraham and Liszi [11], and Marcus [4]. These experimental hydration data were corrected, when necessary, in agreement with the recently established standard value for the Gibbs free energy of solvation for the proton in water, DGsolv ðHþ Þ ¼ 264:0 kcal=mol [12],
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
corresponding to the process of transferring the proton from ideal gas at 1 atm to ideal diluted solution at 1 mol/ l. We adopted a single consistent standard state for defining the Gibbs energy of solvation for all anions according to Ben-Naim, corresponding to the process of transfer of 1 mole of solute per liter of ideal gas to 1 mole/liter of ideal solution [13]. In these standard states, the Gibbs free energy of hydration for the proton is DGsolv ðHþ Þ ¼ 265:9 kcal=mol [10,14]. Transfer free energies between water and nitrobenzene were taken from Refs. [2,15–17]. These data, combined with the Gibbs free energies of hydration, provide the Gibbs free energies of solvation in nitrobenzene. We have collected data for 22 anions dissolved in nitrobenzene. 3. Computational methods The molecular geometries of the 22 anions were optimized in the gas phase, at Hartree–Fock (HF) RHF/631+G(d,p) and Density Functional Theory (DFT) B3LYP/6-31+G(d,p) levels with the GAMESS package [18]. These optimized structures were utilized in the calculations of free energies in nitrobenzene from the PCM. It is wellknown that the use of optimized gas phase structures in condensed phase simulations does not introduce significant errors for solutes with restricted conformational changes [19,20]. In this study, we have used the GEPOL algorithm as implemented in the GAMESS package to create the cavities for the solutes [21]. The GEPOL algorithm generates a cavity for the solute from the positions of the atoms and their radii. Starting from a set of original atomic radii Roi , solvent dependent radii Rcav have been established by adopti ing a factor f to Roi , ¼ f Roi : Rcav i
ð1Þ
The algorithm creates a cavity corresponding to the molecular shape by interlocking spheres of radii Rcav i . Addi˚ were tional smaller spheres with minimum radius of 0.2 A added to smooth the cavity’s surface. The surface of each sphere was divided into 60 discrete elements of area (tesserae). In this work, the original atomic radii Roi from Emsley, modified by the Pisa group [22] were employed as implemented in the GAMESS package. These atomic radii are listed in Table 1. The determination of the suitable scaling factor f from Eq. (1) is one of our main objectives. Therefore, MD of a single anion Cl, Br, I, and CN in nitrobenzene were performed with the GROMACS 3.3 software [23]. Halide anions were described by the OPLS force field [24], CN by the parameters from Ref. [25], and liquid nitrobenzene
363
Table 2 Properties of the liquid solvent nitrobenzene at 298.15 K Property Dielectric constant, [5] ˚ ) [11] Molecular radius (A Density (g cm3) [5]
34.82 2.773 1.198
by the OPLS-AA model [26]. The MD were carried out in the NVT ensemble with 255 nitrobenzene molecules and a single anion arranged in a cubic box with dimensions corresponding to the experimental density (see Table 2) of the solvent at simulation temperature (T = 298 K). Periodic boundary conditions and the minimum image convention were applied. Long range electrostatic interactions were taken into account by the Particle Mesh Ewald method (PME) [27]. Equations of motion were integrated by the leap-frog algorithm [28] with a time step of 1 fs. The bond length of CN anion was fixed by the SHAKE algorithm [30]. The temperature was maintained by the weak coupling Berendsen thermostat with a time constant of 0.1 ps [29]. Data were averaged from production runs of 500 ps preceeded by a stabilization period of 1 ns. The atom–atom pair distribution functions were computed from configurations saved in intervals of 50 time steps. The radial distribution functions (RDFs) for X–H atom pairs, where X represents one of the anion atoms and H the hydrogen atom of the nitrobenzene molecules, have been analyzed in more detail. The distance corresponding to the maximum position in the first peak of these RDFs was utilized to localize the average position of the hydrogen atoms responsible for weak hydrogen bonding between the first shell solvent molecules and the solutes. This structural parameter for the solvent structure around the solutes was used to define the scaling factor from Eq. (1). This methodology of defining the solute–solvent interface in the PCM from MD has been applied before with success in PCM studies of anionic solvation in acetonitrile and DMF [20]. In the PCM formalism, the process of solvation is described by inserting the solute inside a molecular cavity into the dielectric medium representing the solvent. As in our earlier study [20], the free energy of solvation DGsolv is partitioned into the electrostatic energy DGele, the cavitation energy DGcav, and the van der Waals energy DGvdw: DGsolv ¼ DGele þ DGcav þ DGvdW :
ð2Þ
More details about the formalism have been reviewed by Refs. [8,9]. The electrostatic term DGele was calculated from the Integral Equation Formalism version of the PCM, as implemented in the GAMESS package [31,32] at the same
Table 1 ˚ ) for the definition of the solute cavities Atomic radii (in A Element
H
C
N
O
F
S
Cl
Se
Br
I
Radius
1.20
1.70
1.60
1.50
1.35
1.85
1.81
1.90
1.95
2.15
364
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
levels of the gas-phase optimizations, the HF RHF/631+G(d,p) and DFT B3LYP/6-31+G(d,p). Charge density escaping from the cavity boundary is corrected by the charge renormalization procedure of Mennucci and Tomasi [33]. The physical parameters for the nitrobenzene solvent within PCM are summarized in Table 2. The cavitation contribution DGcav was computed from scaled particle theory adapted to molecular cavities as described in Ref. [34]. To computed the van der Waals repulsion–dispersion energy DGvdW, we assume a linear relationship between this term and the solvent accessible surface area of the solute: DGvdW ¼
N X
ni A i :
ð3Þ
i¼1
Ai is the solvent accessible surface area of the ith solute atom with van der Waals parameter ni optimized for this atom [8]. The optimization procedure for ni involve a multiple linear regression procedure minimizing the difference between experimental solvation free energies and the sum of the computed electrostatic and cavitational contributions DGsolv . Thus, the PCM is parameterized from two steps: (i) defining from MD the scaling factors for atomic radii in the cavity creation, and (ii) optimizing the van der Waals parameters for the computation of DGvdW.
Fig. 2. Radial distribution functions for N–H pairs in nitrobenzene solutions.
4. Results and discussion From the MD of the anionic (CN, Cl, Br or I) solutions in nitrobenzene, we computed the RDFs representing solute–solvent atom pair correlations. The most interesting RDFs are those representing the X–H pairs. The H atom of the nitrobenzene molecule possesses a positive charge of q = 0.115e in OPLS-AA [26]. Thus, one might expect sharp first peaks in these RDFs that we have chosen to define the solute–solvent interface in the continuum model. In Figs. 1–5, we have depicted the RDFs for the X–H pairs. Each figure contains three functions corresponding to ortho, meta, and, para hydrogens in the nitro-
Fig. 3. Radial distribution functions for Cl–H pairs in nitrobenzene solutions.
Fig. 4. Radial distribution functions for Br–H pairs in nitrobenzene solutions.
Fig. 1. Radial distribution functions for C–H pairs in nitrobenzene solutions.
benzene ring. For all the X–Hortho pairs, the RDFs present diffuse profiles due to the repulsion between negative charges of the solute and oxygen sites in the nitro group
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
365
Table 3 ˚ 2) optimized for the solvation Van der Waals parameters ni (in kcal/mol A of anions in nitrobenzene
Fig. 5. Radial distribution functions for I–H pairs in nitrobenzene solutions.
of nitrobenzene. On the other hand, RDFs for X–Hpara and XHmeta pairs exhibit well-defined sharp peaks. We observe that the distance corresponding to the maximum in the first peak of the RDFs is slightly shorter for meta hydrogens. Consequently, we utilized the RDFs for meta hydrogens as our standard to establish the solute–solvent interface in the nitrobenzene solutions. For CN, we calculated RDFs for the carbon and the nitrogen sites illustrated in Figs. 1 and 2. The RDF for ˚ the C–Hmeta pair possesses a maximum at rmax = 2.778 A and, for the pair N–Hmeta, a first maximum is located at ˚ . The atomic Emsley radius for carbon is rmax = 2.685 A ˚ , whereas the nitrogen atomic radius is 1.60 A ˚. 1.70 A Therefore, the ratios rmax =Roi are 1.63 and 1.68, respectively. As in our earlier study [20], these ratios are used to define the scaling factor f from Eq. (1). The RDF for the Cl–Hmeta (3) exhibits a first peak with maximum at ˚ . With the chlorine atomic radius, we obtain rmax = 3.086 A the ratio 1.70 defining f. Fig. 4 presents the RDFs for the Br solution. The first maximum in the RDF involving ˚ producing the meta hydrogen is located at rmax = 3.148 A f = 1.61 if the atomic radius for the bromine from Table ˚ ) is used. Finally, in the RDF for the I solution 1 (1.95 A ˚ . The (Fig. 5), we localize a first maximum at rmax = 3.376 A ratio between this distance and the atomic radius of the ˚ ) is f = 1.57. Averaging the five scaling iodine atom (2.15 A factors, we computed an averaged scaling factor of f = 1.64 for anionic solutions in nitrobenzene. In the PCM calculations, we adopted this scaling factor to all the anions investigated in the present study. Next, using the average scaling factor, we computed the electrostatic and cavitation contributions to the Gibbs free energy of solvation for the anions in nitrobenzene. Comparing the sum of the calculated terms with the experimental energy of solvation, we optimized the van der Waals atomic parameters ni of Eq. (3) using the criterium of minimizing differences between experimental and theoretical free energies of solvation. Details of this procedure are described in Ref. [20]. The optimized parameters ni are
Atom
ni (RHF)
ni (DFT)
H C N O F S Cl Se Br I
0.5010 0.0740 0.2499 0.1163 0.3151 0.1625 0.2861 0.6348 0.3780 0.2445
0.5285 0.0710 0.2625 0.1177 0.3282 0.1547 0.2873 0.6050 0.3804 0.2445
shown in the Table 3 for both, the HF and DFT levels of theory. Employing the optimized ni, we computed the van der Waals term of the energy of solvation for the anions according to Eq. (3). Afterwards, DGsolv according to Eq. (2) was computed for the set of 22 anions in nitrobenzene. Tables 4 and 5 summarize these results for the HF and DFT levels, respectively, comparing with the experimental values. The mean absolute deviation between theoretical and experimental data is 2.4 kcal/mol for both, HF and DFT calculations. The mean signed error is 0.6 kcal/mol and 0.7 kcal/mol using HF and DFT levels, respectively. The standard deviation between calculated and experimental values is 3.4 kcal/mol at the HF level and 3.5 kcal/mol for the solvation energies calculated at DFT level Also, it becomes evident from Tables 4 and 5 that the calculations performed at HF or DFT levels produce similar values for the energies of solvation. For small anions, the difference in the electrostatic contribution calculated employing the HF or DFT levels is approximately 0.5 kcal/mol. For most anions in the set, we achieved a very good agreement between calculated and experimental free energies of solvation. No systematic deviations are observed in our results. Larger deviations from experimental data have been computed for HCO 2 , CH3 CO2 , and picrate anions. For the HCO2 anion, we found a deviation of 9.4 kcal/mol at HF level and 9.9 kcal/mol at DFT level. In the case of the CH3 CO 2 anion, the error is 6.5 kcal/mol at HF level and 6.6 kcal/mol at DFT level. For the Pic anion, the observed deviation is negative, 5.4 kcal/mol when the calculations are performed at HF level and 5.5 kcal/mol at DFT level. We believe that the origin of these deviations is the same for the three cases and is related to DGvdW. In our methodology, the van der Waals parameters ni are statistically adjusted to produce the best agreement between the calculated and the experimental free energies of solvation for all the solutes. This method attributes equal statistical weights to all the solutes of the set. With eleven anions in the set containing H, C, and O atoms, the optimization by multivariate fitting favors the minimization of the errors for most of the anions of the subset containing H, C, and O atoms. However, eight of these eleven anions possess larger surfaces and, therefore, the ni
366
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
Table 4 Electrostatic, DGele, cavitational, DGcav, and van der Waals, DGvdW, contributions to the calculated solvation free energies, DGsolv ðcalcÞ, for the anions in nitrobenzene at HF level Anion
DGele
DGcav
DGvdW
DGsolv ðcalcÞ
DGsolv ðexpÞ
Error
54.2 50.3 46.1 56.4 53.8 46.4 46.9 51.8 48.5 48.1 47.7 56.3 54.6 50.1 54.1 49.5 53.5 50.9 44.3 43.7 47.9 38.6
3.9 4.4 5.1 5.7 4.6 6.7 6.9 6.6 7.7 8.1 8.3 6.1 8.7 9.5 11.1 13.7 14.8 12.6 15.2 15.2 15.0 20.2
11.8 18.1 14.2 13.8 7.1 11.6 29.0 8.1 11.7 15.7 10.5 9.3 17.7 18.0 26.2 34.8 28.2 26.5 26.9 26.8 34.9 27.5
62.1 64.0 55.2 64.5 56.3 51.3 69.0 53.3 52.5 55.7 49.9 59.5 63.6 58.6 69.2 70.6 66.9 64.8 56.0 55.3 67.8 45.9
63.6 61.8 55.2 64.3 60.5 51.3 69.0 55.4 49.2 59.9 49.4 68.9 70.1 58.6 69.5 68.9 66.2 66.4 59.1 52.8 65.4 40.5
+1.5 2.2 0.0 0.2 +4.2 0.0 0.0 +2.1 3.3 +4.2 0.5 +9.4 +6.5 0.0 +0.3 1.7 0.7 +1.6 +3.1 2.5 2.4 5.4
Cl Br I N 3 CN SCN SeCN NO 3 ClO 3 BrO 3 ClO 4 HCO 2 CH3 CO 2 F2 CHCO 2 CH3 CH2 CO 2 CH3 CH2 CH2 CO 2 C6 H5 CO2 C6H5O 3-NO2C6H4O 4-NO2C6H4O 2-CH3C6H4O Pic Mean absolute error Mean error RMS error SD
2.4 0.6 3.3 3.4
DGsolv ðexpÞ are experimental data. Units are in kcal/mol.
Table 5 Same as Table 4, but for the DFT computations Anion
DGele
DGcav
DGvdW
DGsolv ðcalcÞ
DGsolv ðexpÞ
Error
Cl Br I N 3 CN SCN SeCN NO 3 ClO 3 BrO 3 ClO 4 HCO 2 CH3 CO 2 F2 CHCO 2 CH3 CH2 CO 2 CH3 CH2 CH2 CO 2 C6 H5 CO2 C6H5O 3-NO2C6H4O 4-NO2C6H4O 2-CH3C6H4O Pic
54.1 50.2 46.1 55.5 53.5 46.4 45.3 51.5 48.4 47.9 47.5 55.6 53.7 49.3 52.9 47.5 52.2 49.7 43.2 42.9 46.6 37.8
3.9 4.4 5.1 5.7 4.7 6.7 7.1 6.7 7.7 8.1 8.5 6.2 8.8 9.6 11.3 13.9 14.9 12.7 15.4 15.4 15.2 20.6
11.8 18.2 14.2 14.7 7.3 11.6 30.8 8.5 11.8 15.9 10.9 9.6 18.6 18.9 27.6 36.9 29.4 27.8 28.1 28.0 36.7 28.8
62.0 64.0 55.2 64.5 56.1 51.3 69.0 53.3 52.5 55.7 49.9 59.0 63.5 58.6 69.2 70.5 66.7 64.8 55.9 55.5 68.1 46.0
63.6 61.8 55.2 64.3 60.5 51.3 69.0 55.4 49.2 59.9 49.4 68.9 70.1 58.6 69.5 68.9 66.2 66.4 59.1 52.8 65.4 40.5
+1.6 2.2 0.0 0.2 +4.4 0.0 0.0 +2.1 3.3 +4.2 0.5 +9.9 +6.6 0.0 +0.3 1.6 0.5 +1.6 +3.2 2.7 2.7 5.5
Mean absolute error Mean error RMS error SD
2.4 0.7 3.4 3.5
E.S. Bo¨es et al. / Chemical Physics Letters 436 (2007) 362–367
parameters were optimized taking into account these large surfaces. Consequently, the obtained parameters describe well the large solutes such as CH3 CH2 CO 2 , CH3 CH2 CH2 CO 2 , and C6 H5 CO2 , but might fail for the smaller solutes HCO 2 and CH3 CO2 . With the ni related to the solvent accessible surface area by Eq. (3), we expect an underestimation of the van der Waals contribution for the small solutes. In the case of the rather large Pic anion, the seven oxygen atoms are exposed to the solvent in external positions of this anion. Thus, the large surface areas of these atoms is responsible for the overestimation of the van der Waals energy. We notice that for most of the anions of the set, DGele is approximately 50 kcal/mol. For smaller and simpler anions, from Cl to ClO 4 in Tables 4 and 5, the electrostatic contribution corresponds to 85–95% of the total energy of solvation. For the larger anions, for example those derived from organic acids, the importance of the electrostatic term is lower and the steric terms (DGcav and DGvdW) become more significant for the solvation free energy. For these organic ions, the van der Waals interactions with the solvent molecules play a very important role in the solvation process in nitrobenzene. 5. Conclusion In the present article, we presented a solvation study of anions in nitrobenzene. We have focused on the calculation of the Gibbs free energy of solvation of 22 representative anions in nitrobenzene from an appropriate parameterization of the PCM for calculations at HF and DFT levels. Relevant structural characteristics extracted from MD of a few anionic solutions in nitrobenzene were utilized to establish the suitable parameter of scaling for the atomic radii in the definition of the solute–solvent interface for the calculations with the continuum model. The scaling factor established for solute cavities in nitrobenzene, f = 1.64, is in agreement with the fact that polar aprotic solvents, such as nitrobenzene, are weak hydrogen bond donors towards anionic solutes [35]. The obtained scaling factor possesses physical sense if one compares with other solventes. For example, the study of anion solvation in N,N-dimethylformamide has been described using a scaling factor f = 1.52 [20]. We have parameterized the PCM and calculated the Gibbs free energies of solvation for 22 anions containing up to 18 atoms. Quite general, the theoretical data agree well with experimental findings corroborating that the parameterized PCM model is able to treat anionic solvation in nitrobenzene. Comparing the parameterization for HF and DFT computations, we conclude that both levels of theory produce qualitatively equivalent solvation free energies. Furthermore, the partitioning of the total energy of solvation into the three contributions, cavitation, electrostatic, and van der Waals, permits better understanding about the role of each of these terms and the nature of the interactions between the solvent nitrobenzene and the anio-
367
nic anion solutes. We believe that the parameterized model can be employed in further studies on the solvation of anionic drugs or other anions with biological relevance. Acknowledgements E.S.B. acknowledges a scholarship from the Brazilian agency CAPES. J.d.A. acknowledges a scholarship from the Brazilian agency CNPq. H.S. acknowledges a grant from Brazilian agency CNPq (309211/2003-4. References [1] A.G. Volkov, Liquid Interfaces in Chemical, Biological and Pharmaceutical Applications, Marcel Dekker, New York, 2001. [2] R. Gulaboski et al., J. Phys. Chem. B 108 (2004) 4565. [3] F. Scholz, Sˇ. Komorsky-Lovric´, M. Lovric´, Electrochem. Commun. 2 (2000) 112. [4] Y. Marcus, Ion Properties, Marcel Dekker, New York, 1997. [5] J.A. Riddick, W.B. Bunger, Organic Solvents Physical Properties and Methods of Purification, John Wiley & Sons, New York, 1986. [6] M. Jorge, R. Gulaboski, C.M. Pereira, M.N.D.S. Cordeiro, J. Phys. Chem. B 110 (2006) 12530. [7] T. Osakai, A. Ogata, K. Ebina, J. Phys. Chem. B 101 (1997) 8341. [8] F.J. Luque et al., Phys. Chem. Chem. Phys. 5 (2003) 3827. [9] J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev. 105 (2005) 2999. [10] J.R. Pliego Jr., J.M. Riveros, Phys. Chem. Chem. Phys. 4 (2002) 1622. [11] M.H. Abraham, J. Liszi, J. Chem. Soc. Faraday Trans. 74 (1978) 1604. [12] M.D. Tissandier et al., J. Phys. Chem A 102 (1998) 7787. [13] A. Ben-Naim, J. Phys. Chem. 82 (1978) 792. [14] D.M. Camaioni, C.A. Schwerdtfeger, J. Phys. Chem. A 109 (2005) 10795. [15] R. Gulaboski, K. Riedl, F. Scholz, Phys. Chem. Chem. Phys. 5 (2003) 1284. [16] Sˇ. Komorsky-Lovric´, K. Riedl, R. Gulaboski, V. Mircˇeski, F. Scholz, Langmuir 18 (2002) 8000. [17] Sˇ. Komorsky-Lovric´, K. Riedl, R. Gulaboski, V. Mircˇeski, F. Scholz, Langmuir 19 (2003) 3090. [18] M.W. Schmidt et al., J. Comput. Chem. 14 (1993) 1347. [19] C. Curutchet, A. Bidon-Chanal, I. Soteras, M. Orozco, F.J. Luque, J. Phys. Chem. B 109 (2005) 3565. [20] E.S. Bo¨es, P.R. Livotto, H. Stassen, Chem. Phys. 331 (2006) 142. [21] J.L. Pascual-Ahuir, E. Silla, J. Tomasi, R. Bonaccorsi, J. Comput. Chem. 8 (1987) 778. [22] J. Emsley, The Elements, Clarendon Press, Oxford, 1991. [23] E. Lindahl, B. Hess, D. van der Spoel, J. Mol. Mod. 7 (2001) 306. [24] N.A. McDonald, E.M. Duffy, W.L. Jorgensen, C.J. Swenson, J. Am. Chem. Soc. 120 (1998) 5104. [25] M. Ferrario, I.R. McDonald, M.C.R. Symons, Mol. Phys. 77 (1992) 617. [26] M.L.P. Price, D. Ostrovsky, W.L. Jorgensen, J. Comput. Chem. 22 (2001) 1340. [27] U. Essman, L. Perela, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577. [28] R.W. Hockney, S.P. Goel, J. Comput. Phys. 14 (1974) 148. [29] H.J.C. Berendsen, J.P.M. Postma, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [30] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [31] E. Cance`s, B. Mennucci, J. Tomasi, J. Chem. Phys. 107 (1997) 3032. [32] E. Cance`s, B. Mennucci, J. Math. Chem. 23 (1998) 309. [33] B. Mennucci, J. Tomasi, J. Chem. Phys. 106 (1996) 5151. [34] J. Langlet, P. Claverie, J. Caillet, A. Pullman, J. Phys. Chem. 92 (1988) 1617. [35] A.J. Parker, Chem. Rev. 69 (1969) 1.