2099
J . Phys. Chem. 1990, 94, 2099-2105
Monte Carlo Studies on Aqueous Solutions of Methylamine and Acetonlttile. Hydration of sp8 and sp Nitrogen William J . Dunn III* and Peter I. Nagyt College of Pharmacy, P.O. Box 6998, The University of Illinois at Chicago, Chicago, Illinois 60680 (Received: January 25, 1989; In Final Form: September 14, 1989)
Monte Carlo simulations have been carried out on dilute aqueous solutions of methylamine and acetonitrile to model the hydration of the sp3 and sp nitrogens. In general, the calculated enthalpy of solvation, relative free energy, and partial molar volume compared favorably with experimental data. Solution structures were considered as relevant to analysis of the hydrogen bonding around the nitrogen atom. There are approximately three water molecules in the first hydration shell around the amino group of methylamine forming two hydrogen bonds to the solute. Acetonitrile has 1.3-1.4 water molecules around the nitrogen atom in the first shell with all of them hydrogen bonding. No hydrogen bonds were found in the methyl regions that are hydrated by 19 water molecules in the first shell for both solutes.
Introduction Since theoretical modeling of the structure of pure water by Stillinger and Rahman' and Lie et al.,2 investigations of the structures of both pure liquid and dilute aqueous solutions have become of central importance. By use of Monte Carlo3 or molecular dynamics4 methods, particular interest has been paid to the hydration of hydrophobic solute^,^ dilute aqueous solutions of small, polar, neutral species,+I2 and organic and inorganic ions.l3-I6 Recently, the investigationswere extended to large, single molecules of biological interest," as well as to enzyme-inhibitor c o m p l e x e ~ . ' ~ JSome ~ calculations used a very powerful procedure known as the free energy perturbation method or statistical perturbation theory, SPT?O to determine free energy differences for some simple solutes in aqueous solution.2'*z2 Despite the biochemical interest of many of the above calculations, none of them deal directly with the solution structure of neutral sp3 and sp nitrogen containing compounds. The NH, group is an important biochemical functional group that is found in peptides and proteins. It is responsible for a number of properties of biopolymers. The cyano group is rarely present in natural molecules. Some synthetic drugs, however, contain this functional group. An investigation of the different hydration potentials of molecules containing nitrogen in two different valence states is of interest purely from a theoretical standpoint. In our previous approximations to build up the first hydration shell of polar molecules, simple methods have been tried. These consider the solvent accessible surface areaz3or the molecular electrostatic potential of solute.z4 They fail, however, to take into consideration the influence of the thermal motion on the equilibrium hydration structure. The more sophisticated quantum chemical supermolecule approach2s has a similar shortcoming. Thus, we report here statistical calculations for neutral compounds with nitrogen atoms in sp3 and sp hybridization states. These investigations supplement previous ones dealing almost exclusively with nitrogen in its sp2 hybridization state.9J0J2*22b Method and Calculations Monte Carlo calculations using periodic boundary conditions were carried out for N,P,T ensembles containing 264 water molecules and one C H 3 N H 2or CH3CN solute molecule. Parameters P and T were l atm and 25 OC, respectively. The software used was provided by Prof. W. L. Jorgensen and adapted in this laboratory to run on an IBM 3081 computer. The details of such calculations are well d e ~ c r i b e d . ~ ~ JHere l ~ J ~only " some computational aspects will be discussed. Calculations using different computers and computational parametersz6 have resulted in somewhat different values for the *To whom correspondence should be
addressed.
'On leave from Chemical Works of Gedeon Richter Ltd., Budapest, Hungary.
0022-3654/90/2094-2099$02.50/0
TABLE I: Parameters for the Potential Function and Geometric Values" u RIA e/(kcal/mol) RIA anale/dep.
MeNH2 Me N
H
0.25 0.20 -1.05 -0.90 0.40 0.35
3.800 3.775 3.300 3.250 0.0 0.0
0.170 0.207 0.170 0.170 0.0 0.0
Me-N
1.456 1.474 1.033 1.011
N-H Me-N-H
121.6 110.5 104.4 106.2
H-N-H MeCN Me
C N
0.15 0.28 -0.43
3.775 3.650 3.200
0.207 0.150 0.170
Me-C C-N
Me-C-N
1.458 1.157 180.0
"OPLS parameters for methylamine: I (upper), ref 27a; I1 (lower), ref 27b. Geometry: I (upper), ref 27a; I1 (lower), ref 28b,c. Case A: parameter I, geometry I. Case B: parameter 11, geometry 11. Case C: parameter I, geometry 11. OPLS parameters for acetonitrile: ref 27c. Geometry: ref 28a.
average interaction energy, density, heat capacity, etc., for liquid water. To obtain correct reorganization energy in the simulations ( I ) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 61, 4973. (2) Lie, G. C.; Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,64, 2314. (3) (a) Valleau, J. P.; Whittington, S.G. In Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 5, pp 137-168. (b) Valleau, J. P.; Torrie, G. M. In ref 3a, pp 169-194. (4) (a) Kushick, J.; Berne, B. J. In Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 5. (b) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press:
Cambridge, 1987.
( 5 ) Pangali, C.; Rao, M.; Berne, B. J. J . Chem. Phys. 1979, 71, 2982.
(6) (a),Owicki, J. C.; Scheraga, H. A. J . Am. Chem. SOC.1977,99,7413. (b) Swaminathan, S.; Harrison, S. W.; Beveridae, - D. L. J . Am. Chem. SOC. i978,100,5705. (7) (a) Okazaki, S.; Nakanishi, K.; Touhara, H.; Watanabe, N.; Adachi, Y. J. Chem. Phys. 1981,74,5863. (b) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J . Phys. Chem. 1985,89, 3470. (8) (a) Linse, P.; Karlstrom, G.; Jonsson, B. J . Am. Chem. Soc. 1984, 106, 4096. (b) Ravishanker, G.; Mehrotra, P. K.; Mezei, M.; Beveridge, D. L. J . Am. Chem. SOC.1984, 106, 4102. (9) (a) Marchese, F. T.; Mehrotra, P. K.; Beveridge, D. L. J . Phys. Chem. 1984, 88, 5692. (b) Jorgensen, W. L.; Swenson, C. J. J . A m . Chem. SOC. 1985, 107, 1489.
0 1990 American Chemical Society
2100
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990
Dunn and Nagy
TABLE 11: Results of the Monte Carlo Simulations for Aqueous Solutions of H20, CH3NH2,and CH3CNa-b H20
Ess
runs 4100K
ETOT -2659.0 (3.3)
Esx -19.9 (0.4)
-2639.1 (3.3)
4000K 6000K 6000K 4000K 6000K 3100K 4600K
-2670.7 -2667.0 -2661.1 -2662.7 -2665.5 -2671.4 -2672.2
-15.5 -15.3 -14.2 -17.2 -17.3 -16.4 -15.8
-2655.2 -2651.7 -2646.9 -2645.5 -2648.2 -2655.0 -2656.4
4Ess 9.9 (4.5)
AH -10.6 (4.5)
AH,,, -10.5'
AV
A Vexp
CH,NH, A
B C CH,CN
(2.2) (2.1) (2.0) (2.4) (2.9) (2.5)
(0.3) (0.3) (0.2) (0.3) (0.3) (0.3)
(2.2) (2.1) (2.0) (2.4) (2.9) (2.5)
3.8 (4.0) 7.3 (3.9) 12.1 (3.9) 13.5 10.8 (4.1) 4.0 (4.4) 2.6 (4.1)
-12.3 -8.6 -2.7 -4.3 -7.1 -13.0 -13.8
(4.0) (3.9) (3.9) (4.1) (4.4) (4.2)
-10.6d
51 58 54 54 39
(17) (14) (12) (14)
41.7'
1 1 5 (12)
-8.3''
95 (12)
49.4f
OEnergies in kcal/mol and volumes in cm3/mol. bError bars are in parentheses for ETOT,Esx, and AV. Values in parentheses for Ess, hEss, and AH are calculated as square root of the sum of squares of error bars for the corresponding terms in eqs 1 and 2. cReference 26a. dReference 39. eReference 32. 'References 37 and 38. $Reference 33. for solutions, calculations were carried out for the pure water system. Intermolecular energies were calculated by using the 12-6-1 OPLS and the TIP4P model for water.26a,b A cutoff radius of 8.5 A was used for both solvent-solvent and solute-solvent interactions. No cutoff corrections2" were applied since it was postulated that the correction for the solvent-solvent interaction is largely removed by subtraction in calculating the reorganization energy. The Lennard-Jones parameters and atomic charges (Table I) were taken from Impey et (set I) and ( I O ) (a) Kuharski, R. A.; Rossky, P. J. J . Am. Chem.Soc. 1984,106,5786. (b) Cristinziano, P.: Lelj. F.: Amodeo, P.; Barone, V. Cfiem. Phys. L e x 1987, 140, 401. ( I I ) (a) Jorgensen, W. L.; Madura, J. D. J . Am. Chem. SOC.1983, 105, 1407. (b) Alagona, G.; Tani, A. J . Mol. Struct. (THEOCHEM) 1988, 166, 375. (12) (a) Cieplak, P.; Kollman, P. A. J . Am. Chem. SOC.1988, 110, 3734. (b) Beveridge, D. L.; Maye, P. V.; Jayaram, B.; Ravishanker, G.; Mezei, M. J . Biomol. Srrucr. Dyn. 1984, 2, 261. (13) (a) Alagona, G.; Ghio, C.; Kollman, P. A. J . Am. Cfiem.SOC.1985, 107,2229. (b) Jayaram, B.: Mezei, M.; Beveridge. D. L. J . Am. Cfiem. SOC. 1988. 110. 1691. (14) Jorgensen, W. L.; Buckner, J. K.; Huston, S. E.; Rossky, P. J. J . Am. Chem. SOC.1987, 109, 1891. (15) (a) Mezei, M.; Mehrotra, P. K.; Beveridge, D. L. J . Biomol. Struct. Dyn. 1984, 2, I . (b) Alagona, G.: Ghio, C.; Kollman, P. A. J . Am. Chem. SOC.1986. 108. 185. (c) Joraensen. W. L.; Gao. J. J. Phys. Chem. 1986, 90, 2174. (16) (a) Chandrasekhar,J.; Spellmeyer, D. C.; Jorgensen, W. L. J . Am. Cfiem. SOC.1984, 106, 903. (b) Elliott, R. J.; Goodfellow, J. M. J . Theor. Biol. 1987, 127,403. (c) Elliott, R. J.; Goodfellow, J. M. J. Theor. Biol. 1987, 128, 121. (17) (a) Corongiu, G.; Clementi, E. Biopolymers 1982, 21, 763. (b) Jorgensen, W. L.; Tirado-Rives,J. J . Am. Chem. SOC.1988, 110, 1657. (18) (a) Ahlstrom, P.; Teleman, 0.;Jonsson, B.; Forsen, S. J . Am. Chem. SOC.1987, 109. 1541. (b) Clementi, E.; Chin, S.; Logan, D. Isr. J . Cfiem. 1986, 27, 127. (c) Tapia, 0.;Eklund, H.; Branden, C. I. In Steric Aspects of Biomolecular Inferacfions;Naray-Szabo. G., Simon, K., Eds.; CRC Press: Boca Raton, FL, 1987; pp 159-180. (19) (a) Lybrand, T. P.; McCammon, J. A,; Wipff, G. Proc. Natl. Acad. Sci. U.S.A. 1986,83,833. (b) Wong, C. F.; McCammon, J. A. J . Am. Chem. SOC.1986, 108, 3830. (c) Bash, P. A,; Singh, U. C.; Brown, F. K.; Langridge, R.; Kollman, P. A. Science 1987, 235, 574. (20) Zwanzig, R. W. J . Chem. fhys. 1954, 22, 1420. (21) (a) Jorgensen, W. L.; Ravimohan, C. J. Chem. Phys. 1985,83, 3050. (b) Mezei, M.; Mehrotra, P. K.; Beveridge, D. L. J . Am. Chem. SOC.1985, 107, 2239. (c) Jorgensen, W. L.; Buckner, J. K. J. Phys. Chem. 1987, 91, 6083. (d) Jorgensen, W. L.; Gao, J. J . A m . Cfiem. SOC.1988, 110, 4212. (22) (a) Singh, U.C.; Brown, F. K.; Bash, P. A,; Kollman, P. A. J . Am. Chem. SOC.1987, 109, 1607. (b) Cieplak, P.; Bash, P.; Singh, U.C.; Kollman, P. A . J . Am. Cfiem. SOC.1987, 109, 6283. (23) (a) Dunn,W. J., 111; Grigoras, S.; Koehler, M. G. J . Med. Cfiem. 1987, 30, 1121. (b) Koehler, M. G.; Grigoras, S.; D u m , W. J. 111. Quanr. Struct.-Act. Relot. 1988, 7, 150. (24) (a) Nagy, P. J . Mol. Srruct. (THEOCHEM) 1988, 181, 361. (b) Nagy, P. J . Comput.-Aided Mol. Des. 1988, 2, 65. (c) Alagona, G.; Ghio, C.; Nagy, P. J . Mol. Srrucr. (THEOCHEM) 1989, 187, 219. (25) Pullman, A. In Quantum Theory of Chemical Reactions; Daudel, R., Pullman, A., Salem, L., Veillard, A., Eds.; Reidel: Dordrecht, 1980 Vol. 11. pp 1-32. (26) (a) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (b) Jorgensen, W. L.; Madura. J. D. Mol. fhys. 1985,56, 1381. ( c ) Jorgensen, W. L. Chem. fhys. Lett. 1982, 92,405. (d) Jorgensen, W . L. J . Am. Cfiem.Soc. 1981,103,335. (e) Madura, J. D.; Pettitt, B. M. Mol. f h y s . 1988, 64, 325. (27) (a) Impey, R. W.; Sprik, M.; Klein, M. L. J . Am. Chem. SOC.1987, 109, 5900. (b) Jorgensen, W. L.; Briggs, J. M . J . Am. Chem. Soc. 1989, I l l , 4190. (c) Jorgensen. W. L.; Briggs, J. M. Mol. Pfiys. 1988, 63, 547.
Jorgensen and B r i g g (set ~ ~ 11) ~ ~for CH3NH2and from Jorgensen The united atom approximation for and B r i g g for ~ ~CH3CN. ~~ the methyl group was used for both solutes. The molecular geometry of acetonitrile corresponds to that determined by microwave spectroscopy.28a Geometry I for methylamine was fitted in ref 27a to obtain optimal results for liquid methylamine in a molecular dynamics calculation. Geometry I1 corresponds to the experimental structure.28bcThere were three different calculations for methylamine, A (geometry I, set I), B (geometry 11, set 11), and C (geometry 11, set I), with results in Table 11. Preferential sampling due to Owicki and Scheraga30a and JorgensenMbwas used to enhance the rate of convergence for dilute solutions. The probability of a new configuration generated by translating and rotating a single solvent molecule was proportional to 1/(R2 + c ) . R is the distance between the solute and solvent molecules, and c is a constant set here equal to 150. Generation of a new configuration by translating and rotating the solute was attempted on every 50th step. A volume change was tried after 1000-1500 configurations. The proportion of the accepted new configurations as compared to all those generated was 41-42%. We used lOOK blocks of simulations to obtain final error bars based on the separate averages for blocks. In some parts of this work random sampling was used, and the results were compared with those obtained by preferential sampling. Preferential sampling was used in the equilibration phase with lOOOK configurations for acetonitrile and with 3000K configurations for methylamine. There were 3 1OOK configurations taken for acetonitrile and 4000K configurations for methylamine in the averaging phases. Continuing these simulations with random sampling, 1500K further configurations were considered for acetonitrile and 3000K for methylamine. Thus, altogether 4600K and 6000K configurations were taken for acetonitrile and methylamine (case A), respectively, in the averaging phases. This gives a ratio of 2:1 for the number of configurations with preferential and random sampling. Calculations for methylamine using only preferential sampling are presented in cases B and C. Relative free energy calculations between methylamine and acetonitrile were carried out in cases A and C. A linear coupling parameter, A, was with a change of 0.125. Reference states corresponded to X = 0.0 (methylamine), 0.125,0.375,0.625,0.875, and 1.0 (acetonitrile). An equilibration phase of lOOOK configurations was followed by 2000K configurations in the averaging phase. Our previous results29showed a fast convergence of the relative free energy with no considerable change after 1500K configurations. Results and Discussion Thermodynamics. Thermodynamic data calculated for the dilute aqueous solutions of CH,NH2 and CH3CN are given in Table 11. Results obtained by us for a system with 263 water (28) (a) Costain, C. C. J . Chem. ffiys. 1958, 29, 864. (b) Nishikawa, T. J . Phys. Soc. Jpn. 1957, 12,668. (c) Fink, W. H.; Allen, L. C. J . Cfiem. Phys. 1967, 46, 2276.
(29) Nagy, P. I.; Dunn, W. J., 111; Nicholas, J. B. J . Cfiem. f f i y s . 1989. 91, 3707. (30) (a) Owicki, J. C.; Scheraga, H. A. Chem. Phys. Lett. 1977, 47,600. (b) Jorgensen, W. L. J . Pfiys. Cfiem. 1983, 87, 5304.
Aqueous Solutions of Methylamine and Acetonitrile
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2101
where Esx and Ess are the solute-solvent and solvent-solvent interaction energies in solution. The solvent reorganization energy is AEss = E= - Eso, the difference of the solventsolvent energies for the solution and pure solvent. The latter value was calculated here as -2659.0 kcal for 264 solvent molecules. The solvation enthalpy for water is -10.6 kcal, in excellent agreement with the experimental value of -10.51 kcal (Table 11). The calculated solvent reorganization energy around a H 2 0 solute is 9.9 kcal, and the solutesolvent interaction energy is -19.9 kcal. The corresponding values were reported as 10 and -20 kcal.lla Results obtained by preferential sampling followed by random sampling are shown for CH3NH,(A) and CH3CN. The error bars (given in parentheses in Table 11) for ET,, are of similar magnitude within both pairs. The Esxvalues are well-determined, having an error of only 0.3 kcal. Therefore, the solvent-solvent energies have almost the same uncertainties as the total energies (see eq 1). The error bar for Esso is 3.3 kcal obtained from ETOT for pure water and is -4 kcal for the solvent-solvent reorganization energy. In conclusion, there is no significant difference within the pairs of values in Table 11. The uncertainties are of f 0 . 3 kcal for Esx and f2-4 kcal for ETOT, Ess, and AEss. Calculations B and C for methylamine using only preferential sampling result in the same uncertainties after considering 6000K configurations. The Esx values above are less negative than that for methanol, -18 kcal,IIass and the present value of -19.9 kcal for the interaction energy of the H 2 0 solute with the solvent. These values indicate weaker solutesolvent interactions in the solutions of N-containing molecules as compared to 0-containing ones. This finding emphasizes the importance of considering a large number of solvent molecules in investigating solvent effects since ab initio calculations with large basis sets31 have found the H,N.-HOH interaction energy to be more negative than that for the water dimer. The calculated solvation enthalpy for methylamine is 7.1-8.6 kcal, differing by 2.0-3.5 kcal from the experimental value in cases A and C. Parametrization B results in a less acceptable value. The enthalpy of solvation for acetonitrile differs by 4.7-5.5 kcal from the experimental value. The order of the calculated mean AH values is reversed when compared to experiment. This may be a weak point of our calculations. However, since the difference between the experimental solvation enthalpies of CH3NH2 and CH3CN is only 2.3 kcal, which is within the statistical uncertainty of our calculations, the order of the limiting values cannot be decided at this level. The calculated final volumes are from 7915 f 12 to 7948 f 15 A3 for the methylamine and 8009 f 11 A3 for the acetonitrile solutions. Taking a volume of 7851 f 17 A3 calculated for pure water, the partial molar volumes, AV, are from 38.5 to 58.4 cm3 for methylamine and 95.1 cm3 for acetonitrile. The values for methylamine compare favorably with the experimental value of 41.7 cm3,32while the calculated partial molar volume for acetonitrile is overestimated by 93%, in comparison with the experimental value of 49.4 cm3.33 The convergence in calculating the
partial molar volume is rather slow.29 The good agreement between the experimental and calculated values of AV for methylamine in case C is attributed to the long simulation, using only preferential sampling to enhance the convergence rate. The effect of the longer simulation is reflected also in the case of acetonitrile in comparing the AVvalues after 3100K and 4600K configurations. Before discussing the results of the free energy calculations, one should compare the dipole moments of the solutes due to the different parametrization. In case A for methylamine, the gasphase value of 1.3 1 D34was maintained by the balance of a large bond polarity and the rather flat pyramidal structure of the amino group. In case B the bond polarities decreased, but due to the near-tetrahedral bond angles, the overall dipole moment increased to 1.68 D. Enhancing the bond polarities and preserving the geometry in case C, the dipole moment reached 1.94 D. Table I1 shows that the highest dipole moment with the set of larger bond polarities results in the most negative solute-solvent interaction energy (case C). However, the larger bond polarity seems to be more important than the higher dipole moment when comparing the Esx values for A and B. The formula of ZwanzigZ0shows that the sign of the relative free energy depends on the sign of the difference for the Esx's in question. The relative solvation free energy, AAG = AG(methylamine) - AG(acetonitrile), is -0.7 kcal from experimental values.35 This suggests a more negative solutesolvent interaction energy for methylamine than for acetonitrile. Our calculated value for methylamine is -15.3 kcal in case A, close to that for acetonitrile of -15.8 kcal. Esx is considerably more negative in case C and less negative in case B (Table 11). Thus, we calculated AAG only for cases A and C. In case B also the AH value was surprisingly poor, indicating an inappropriate parametrization. Reliability of the acetonitrile parameters was supported indirectly by our preliminary calculations for the transformation of methanol , ~ ~the to a ~ e t o n i t r i l e : the ~ ~ experimental value is 1.2 k ~ a l and calculated value is 1.6 kcal. The free energy difference between methylamine and acetonitrile solvation was calculated as +0.74 kcal in case A and -1.93 kcal in case C. In the latter case the sign is correct, but the value is in error by some 1.2 kcal. In both cases we used the same OPLS parameter set of Impey et al.27a The qualitatively different values show the high sensitivity of the relative free energy results to the molecular geometry of methylamine. An intermediate geometry between I and I1 in Table I may give good agreement. The correct sign for AAG in case C, as well as the good partial molar volume obtained with it, indicates a geometry near the experimental one to be preferable for Monte Carlo simulations of the aqueous solution of methylamine. In this paper we are mainly interested in the solution structures in the vicinities of the solutes. Neither the radial nor the energy distribution functions have been considerably modified in going from case A to C (see next section). Thus, we conclude that there is no need to optimize the methylamine geometry in order to obtain relevant structural information about the solution structure. Solution Structure. Analysis of the radial distribution functions provides information about the solution structure. These functions are very similar in cases A and C for methylamine, and so the results of them will be discussed together. Figure 1 shows the 0-0 radial distribution function for water and the solutions of CH3NH2and CH3CN. There is very close coincidence of the first maxima (2.75-2.77 h;) and minima (3.40 A). The second maxima (4.38-4.50 A) are also very similar.
(31) (a) Alagona, G.; Ghio, C.; Cammi, R.; Tomasi, J. In Molecules in Physics, Chemistry and Biology; Maruani, J., Ed.; Reidel: Dordrecht, 1988; Vol. 2, pp 507-547. (b) Szczesniak, M. M.; Scheiner, S.J . Chem. Phys. 1986, 84,6328. (c) Latajka, Z.; Scheiner, S.J . Comput. Chem. 1987, 8, 663. (d) Del Bene, J. E. Int. J . Quantum Chem., Quantum Biol. Symp. 1987, No. 27, 14. (32) (a), Cabani, S.;Conti, G.; Lepori, L.J . Phys. Chem. 1974, 78, 1030. (b) Shahidi, F.; Farrell, P. G.; Edward, J. T. J . Chem. SOC.,Faraday Trans. I 1977,715. (33) Janelli, L.; Pansini, M.; Jalenti, R. J. Chem. Eng. Data 1984, 29, 266.
(34) Weast, R. C., Ed. CRC Handbook of Chemistry and Physics, 61st ed.; CRC Press: Boca Raton, FL, 1980-81. (35) Pearson, R. G. J . Am. Chem. SOC.1986, 108,6109, and references therein. (36) Nagy, P. I.; Dunn, W. J., 111. Presented at the XXII Midwest Theoretical Chemistry Conference, Indianapolis, May 1989. (37) An, X.-W.; Mansson, M. J . Chem. Thermodyn. 1983, IS, 287. (38) (a) Morcom, K. W.; Smith, R. W. J . Chem. Thermodyn. 1969, I , 503. (b) Tomkjns, R. P. T.; Turner, P. J. J . Chem. Thermodyn. 1977,9, 707. (39) Felsing, W. A,; Wohlford, P. H. J . Am. Chem. SOC.1932, 54, 1442.
molecules as solvent and one molecule as solute are also given. Accepting the usual approximations for N,P, T ensemb l e ~ ~ , ~ ~ , and ~ ~ disregarding , I l ~ , l ~ ~ the , ~ change in the internal energy of the solute by solvation, the total intermolecular energy (ETOT) and heat of solvation ( A H ) are given as
2102
The Journal of Physical Chemistry, Vol. 94, No. 5. 1990
Dunn and Nagy
CH3CN
0 1 2.0
I
4.0
I
5.0
4.0
3.0
3.0
2.0
,
1
6.0
Figure 3. N-.O radial distribution functions in aqueous solution for
R /A
CH,NH2 (top) and for CH3CN (bottom).
'
Figure 1. 0-0 radial distribution functions for the solvent in solution of H 2 0 (solid line), CH3NH2(+), and CHJN (m).
1
o'2
CH3NH2
d
I
1.0
-30
-20
0
-10
Bonding energy lkcal/moll
2.0
3.0
I
4.0
5.0
6.0
Figure 4. N-H(water) radial distribution functions in aqueous solution for CHJNH2 (top) and for CH$N
(bottom).
Figure 2. Solvent-solvent bonding energy distribution (a). Solute-solvent energy distribution in aqueous solution: CH3NH2C (b), CH3NH2 A (c),and CHJN (d).
The solvent-solvent and solute-solvent bonding energy distribution functions are shown in Figure 2. The former is unaltered in going from pure water to solution. The solute-solvent interaction energy curves are higher and narrower than the solventsolvent ones at their peak values. The peak of the distribution function in case C for methylamine moved toward a more negative value as compared to case A, in agreement with the more negative value for Esx. The solute-solvent radial distribution functions are shown in Figures 3-5. The N-0 function for C H 3 N H 2(Fi ure 3, upper curve) has distinct maxima at 2.80-2.85 and 4.25 , a shoulder at 3.10 A, and minima at 3.30-3.40 and 5.40 A. Integrating up to the first minimum as the boundary of the first hydration shell gives a coordination number of 3.0-3.2 for the oxygen of water in this interval. The N.-H(water) radial distribution function (Figure 4, upper curve) has the first maximum at 1.85 A and minimum at 2.45 A. Integration up to the latter value gives a coordination number of 1.1-1.3 for the water hydrogens around the amine nitrogen. Radial distribution functions referring to possible hydrogenbonding sites with acceptor waters are shown in Figures 5 and 6. The (N)H--O function (Figure 5) has a peak and shoulders
1
1.0
2.0
3.0
4.0
5.0
R
1;'
Figure 5. (N)H-O radial distribution function for CH3NH2in aqueous solution.
in the 1.60-2.50-A region and a broader and higher peak at 3.15-3.30 A. Integrating up to the first minimum indicates that
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990 2103
Aqueous Solutions of Methylamine and Acetonitrile
r\
CH3NH2
0
2.0
3.0
4.0
5.0
R
/a
6.0
Figure 6. Methyl-0 radial distribution functions in aqueous solution for CH3NH2and CH3CN.
each amine hydrogen is associated with 0.6-0.7 water molecules. Considering their largely overlapping spheres (see Figure 8), this corresponds to 1 water molecule being around the hydrogens. The CH3.-0 radial distribution function for methylamine (Figure 6) indicates a loosely structured first shell between 2.9 and 5.5 A, containing 21-23 water molecules. The N - 0 radial distribution function for CH3CN (Figure 3, lower curve) does not show a well-defined peak. There is only a shoulder which does not change during the course of simulation. Integrating up to 3.12 A on this curve results in a coordination number of 1.3-1.4 with water oxygens in the first hydration sphere. The Ne-H(water) radial distribution function has its first maximum at 1.95-2.00 8, (Figure 4, lower curve). The coordination number of the water hydrogens around the nitrile nitrogen is 1.3-1.4, obtained from integrating the radial distribution function up to its first minimum at 2.50 A. The CH3.-0 radial distribution function for the acetonitrile has the first maximum at 3.6 A. Integration to the first minimum at 5.5 A gives -21 water molecules around the methyl group. The second maxima of the Ne-0 and N-.H(water) radial distribution functions for methylamine are near 4.25 and 3.25 A, respectively. The functions calculated in these regions for contributions from water molecules in the second hydration shell and from those in the vicinity of the methyl site have broad bands. Second bands of the N-0 and N--H radial distribution functions for acetonitrile are very flat and poorly defined. On the contrary, the (N)H-.O function for methylamine has a characteristic second band with a maximum value considerably greater than that of the first one. The solute-solvent energy pair distribution functions are given in Figure 7 . Methylamine has a well-defined maximum in the interval -7.8 to -4.1 kcal and a definite, but poorly defined, shoulder between -4.1 and -2.4 kcal. Integrals to -4.1 and -2.4 kcal give values of 0.9 and 2.0, respectively, for set A and values of 1.1 and 2.1 for set C. These results indicate a 1:l ratio of water molecules for both kinds of hydrogen bond. Acetonitrile has a smooth shoulder between -4 and -3 kcal; integration up to -2.5 gives a value of 1.4 for the number of hydrogen bonds. From the results discussed above it is possible to suggest a model of the first hydration shells for methylamine and acetonitrile. More direct information may be obtained by statistical analysis of coordinates for the solvent molecules, but relevant, although indirect, information can also be obtained from analysis of radial and energy pair distribution functions. From the N - 0 radial distribution function of methylamine, there are three water molecules in the first hydration shell of the amino group. This sphere has a radius of 3.3-3.4 A (see Figure 3). Based on the energy pair distribution function, -2 water molecules are hydrogen bonded to the solute. Calculations for
-
-
Interaction energy
lkcal/moll
Figure 7. Solute-solvent energy-pair distribution for CH3NH2C (a), CH3NH2A (b), and CH3CN (c) in aqueous solution. TABLE 111: OPLS Interaction Energies of the CH3NH2-H,O and CH3CN-H20 Dimers at Selected Arrangements‘
H ‘
H, C H ~ - C E N N - - - H -0
CHs-
C
180.0
-6.51
2.00 2.30 2.05
180.0 180.0 158.4
-3.67 -1.38 -3.58
2.30
107.0
-1.04
2.96 3.12
2.00 2.16
180.0 180.0
-4.60 -4.20
2.97
2.50
110.1
-3.28
3.05
2.09
180.0
-1.39
3.03 3.33 3.03
H
CH,-NI~-*H \H,--O’
1.89
sN
I O.
H
HI,,. N-CH,.--0
I
H
NEC-CH3.--0
/H
3.80
180.0
-0.50
3.60
180.0
+2.71
‘H /H
H ‘
Distances in A, angles in deg, and energies in kcal/mol. the pair-interaction energy between water and methylamine in case A (Table 111) show that the energy of a linear N-H-0 dimer at a separation of 2.85 A is -6.5 kcal, compatible with the strong hydrogen-bond energies in the first region of the energy pair distribution function in Figure 7. Thus, one of the water molecules forms a N.-H-0 hydrogen bond to the solute (see 0,in Figure 8). The difference of the minimal N-H and N-eO distances is about 1 A (see Figures 3 and 4), which is nearly equal to the 0-H bond length in the water molecule. These separations suggest a linear hydrogen bond. The explanation is supported by the position of the second peak of the N-H radial distribution function (Figure 4). Its location is about 3.25 A, near the theoretical value of 3.12 A for a linear hydrogen bond. Further support comes from the CH,--H(water) radial distribution function that has a maximum at 2.60-2.80 A (see for distances Figure 8) and contains 1 .l-1.3
2104
Dunn and Nagy
The Journal of Physical Chemistry, Vol. 94, No. 5, 1990
/
A '
'"/
Figure 8. Scheme for the first hydration shell of the amino roup of CH3NH2in aqueous solution. Spheres with radii of 2.5 and 2.9 around the nitrogen and methyl group, respectively, are excluded for water oxygen; sphere of radius of 1.5 8, around nitrogen is excluded volume for water hydrogens as well.
f
water hydrogens within a sphere of 2.85 A. Based on these values, a linear hydrogen bond would have a CH3-N-0 angle of 101-1 16'. which is near the direction of the nitrogen lone pair. The second hydrogen-bonded water molecule in the first hydration shell (0,in Figure 8) is not as tightly bound as the first. Calculations for linear and bent N-H-0 dimers (Table 111) gave interaction energies between -3.7 and -1.0 kcal. The more negative values are highly compatible with the region of the weaker hydrogen bonds. The third water is not hydrogen bonded to the solute but may form a bridge between the two other waters (0, in Figure 8). One of its hydrogen atoms points toward the nitrogen, contributing to the coordination number of 1.1-1.3 for nitrogen with hydrogen. No symmetrical arrangement of O2and O3around the amine hydrogens was postulated since it should give a value of > I for the number of the N-He-0 hydrogen bonds, contrary to our finding. A bifurcated hydrogen-bond structure, however, was examined. The data in Table 111 show that this arrangement is unfavorable. Thus, a linear N-H--0 hydrogen bond being formed and broken around both amine hydrogens is more favorable, resulting in a solvation of each hydrogen on an average of 50%. Our hydration model is much simpler for CH3CN. Table 111 shows that the energies of both the collinear and bifurcated dimers are compatible with the hydrogen-bond energies (-5.5 to -2.5 kcal), but bifurcation corresponds to a H(water) coordination number twice the value of that for O(water). By integration of the radial distribution functions to the end of the plateau at 3.12 8, (Figure 3) and to the first minimum (Figure 4), the value of 1.3-1.4 was obtained for both coordination numbers around the nitrogen. Thus, only a linear model is consistent here. Since this number of water molecules equals the value of waters having an interaction energy more negative than the limit of -2.5 kcal, all water molecules in this sphere are considered H bonded. The energy of an arrangement with the 0 - H bond perpendicular to the axes of acetonitrile is -1.4 kcal, which is out of the H-bond energy interval. Furthermore, a stable arrangement in the perpendicular orientation would require a considerable value of the C-H(water) radial distribution function at a distance of 2.4 A. The corresponding G ( R ) value is very small here, indicating that there is no water in this orientation. Since the number of H-bonded water molecules is greater than 1 in the first shell, all of the N-H-0 bonds cannot be simultaneously collinear to the C-N bond. A straight hydrogen bond with C N O angle differing considerably from 180° may, however, be also preferable. Jorgensen and B r i g g ~ found ~ ' ~ an optimal value of 166O for the C N O angle with a b initio calculations. According to our model (Figure 9), there are two water molecules spending on average 65-70% of their time (corresponding to 1.3-1.4 molecules) within a sphere of 3.12 %I around the nitrogen, forming nearly linear hydrogen bonds. The maxi-
Figure 9. Scheme of the first hydration shell of the cyano group of CH,CN in aqueous solution. For the meaning of spheres with radii of 1.5, 2.5, and 2.9 A, see Figure 8.
mum of the radial density function for N-H(water) is at 2.1 A, in good accordance with an 0 location of 3.06 A from N. Calculations for linear dimers where the water is located at the edge of the methyl group (Table 111) show that, at distances corresponding to the maxima of the CH3-0 radial distribution functions, the interaction energy is slightly negative (CH3NH2) or strongly positive (CH3CN). This indicates that the interaction with water at the methyl edge should not contribute to the hydrogen-bond energy; thus, only molecules belonging to the nitrogen sphere are hydrogen bonded to the solute. The first hydration shell of the methyl group contains 21-23 water molecules at both compounds. This value, however, includes water molecules hydrating the nitrogenic site, as well, since the first shells of the polar groups are embedded in those of the methyl groups (see Figures 8 and 9). A correct analysis should involve investigations of the proximity introduced by Beveridge et aLgaqm A rough estimation of the coordination number of the methyl group may, however, be given by subtracting the number of the water molecules considered to belong to the shell of the nitrogen atom from the value obtained for the first shell of the methyl group. The values calculated in this way are 18-20 and 19 for methylamine and acetonitrile, respectively. Upon comparison of the present results for methylamine with those for dilute aqueous methanol,"a!29 considerable similarities are found. The first maxima of the X-O(water) (X = 0(methanol), N ) radial distribution functions are in the interval of 2.75-2.85 A. The first minima, considered as the boundary of the first hydration shell, are in the range 3.30-3.40 A. There are 3.0-3.2 water molecules around the polar sites. The number of the water molecules hydrogen bonded to them is 2.0-2.4. The first shell of the methyl group contains 18-20 water molecules at methylamine and 19-20 waters for methanol. ConcI usions
Calculated enthalpies of solvation of gaseous C H 3 N H z and CH3CN in water agree with the experimental values within 2.0-5.5 kcal. The relative free energy differs from the experimental value by 1.2 kcal due to a lack of geometry optimization of the methylamine solute. Good agreement with the experimental partial molar volume has been achieved for methylamine. This quantity is overestimated, however, for acetonitrile in a shorter simulation. The first hydration shell of the amino group in methylamine (with a radius of -3.3 A) contains three water molecules with two of them hydrogen bonded to the solute. The third water molecule appears to form a bridge between the other two. The boundary of the first hydration shell of CH3CN is not as structured. The first shell with a radius of 3.12 8, contains 1.4 water molecules with linear N-H-O hydrogen bonds to the solute. The coordination number of the methyl group with water is 18-20 for methylamine and 19 for acetonitrile. There is no hydrogen bond between the water molecules and the methyl group (40) Mehrotra, P.K.:Beveridge, D.L. J . Am. Chem. SOC.1980,102,4287.
J . Phys. Chem. 1990, 94, 2105-2113 at these solutes since the pair-interaction energy is hardly negative even in favorable orientations.
Acknowledgment. The authors are grateful to Prof. W. L. Jorgensen for making possible the use of his Monte Carlo program
2105
and to Dr. G. Glen for valuable discussions. Also, the authors acknowledge the technical assistance of Dr. R. Goldstein of the Academic Data Network of the University of Illinois at Chicago. Registry No. MeNH2, 74-89-5; MeCN, 75-05-8.
Monte Carlo Simulations of an Electric Double Layer Bo Svensson,* Bo Jonsson, and C. E. Woodward Physical Chemistry 2, Chemical Centre, POB 124, S-22100 Lund, Sweden (Received: February 6, 1989; In Final Form: July 20, 1989)
Monte Carlo simulations of an electric double layer in an electrolyte solution have been performed in the canonical ensemble. The chemical potential of the electrolyte was evaluated with a modified Widom method. The effect of salt concentration, surface charge density, and ionic size on the ion distributions and on thermodynamic quantities was investigated. Particular emphasis was put on the evaluation of the osmotic pressure, and the results were compared with predictions from the Poisson-Boltzmann equation. The Poisson-Boltzmann results are in good agreement with simulation for systems containing monovalent counterions. In the presence of divalent counterions, the approximate theory fails and there is a qualitative difference between Poisson-Boltzmann and the exact results. For a 1:2 salt at sufficiently high concentration, there appears to be a negative net osmotic pressure as a result of the different correlation effects in the bulk and in the double layer.
1. Introduction Electrolyte solutions containing highly charged aggregates occur in numerous chemical and biochemical systems. Though the chemical composition and size of aggregates may differ widely, the long-ranged Coulombic interaction often gives these systems several common properties. These electrical properties are usually discussed in terms of the double layer. A double layer is formed as a result of the strong attraction between an aggregate and its counterions and repulsion of its co-ions. This causes an accumulation of opposite charge close to the aggregate surface, which is counteracted by the loss in entropy. The result is a spatially diffuse charge distribution adjacent to the aggregate surface. When two such aggregates approach, their double layers interact, and, depending on the conditions, they may appear to either attract or repel one another. One very important system where the double-layer interaction plays a fundamental role is in the stability of charged colloids.’** In fact, a repulsive double-layer interaction is the main stabilizing factor here. The double-layer interaction is mediated entirely by electrostatic interactions and can be modulated in a number of different ways, e.g., by varying the ambient salt concentration or by changing the salt valency. One may also be able to influence the charge of the colloid, eg., in AgI sols the surface charge density can be modulated by the addition of iodide ions. Another system of considerable technical importance, where double-layer interactions play an important role, is solutions containing charged surfactants. Surfactant molecules are known to form a number of aggregates of different geometry (micelles, liquid crystals, etc.), and the thermodynamic stability of these phases can be rationalized on the basis of electrostatic interactiom3 The theoretical foundation for understanding of all these types of systems has long since been based on the so-called DLVO theory.’V2 The DLVO theory for a solution of like-charged aggregates has two main ingredients. The first is a double layer repulsion between particles, which is described by the PoissonBoltzmann (PB) equation, within the framework of the primitive model. The second is the attractive van der Waals interaction, ( I ) Derjaguin, B. V.; Landau, L. Acta Phys. Chim. URSS 1941, 14, 633. (2) Verwey, E. J. W.; Overbeek, J. Th. G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (3) Jonsson, Be.; Wennerstrom, H. J . Colloid Interface Sci. 1981, 80,482.
0022-3654/90/2094-2105$02.50/0
arising from a difference in dielectric permittivity between the aggregates and the solvent. We will focus only on the double-layer interaction. In previous work, we used Monte Carlo (MC) simulations to test the accuracy of the PB theory for interacting double layers. We treated the case where the only small ion species were point counterions.“ The mean field approximation of the PB equation was shown to be excellent for weakly coupled systems, e.g., monovalent counterions in a high dielectric ~olvent.~JOn the other hand, in the presence of divalent counterions, the PB equation failed to even qualitatively describe the nature of the interaction between charged aggregates6 This is despite the fact that the calculated counterion profiles appeared to be in good agree~nent.~ This failure of the PB theory is due to its complete neglect of ionic correlations between the mobile cunterions-the mean field approximation assumes the pair correlation function for small ions to be unity, Le., Similar conclusions emerged from nonuniform integral equation theories such as the anisotropic hypernetted chain (AHNC) equation.’ In these earlier simulations, two like-charged aggregates were modeled as charged, parallel walls, which contained between them the neutralizing ionic solution. For this system, ionic correlations are manifested in two separable contributions to the net osmotic pressure between the walls. These are most easily described by dividing the region between the charged walls into two ‘double layers” with an imaginary, centrally placed, parallel plane. The first contribution is the repulsive entropic term, which is proportional to the counterion density at the midplane. It is in fact the only term that appears in the PB theory. The simulated value to this term is significantly smaller than the result from the PB approximation, due mainly to correlations between ions on the same side of the midplane (intra-double-layer correlations). The second contribution to the pressure is due to ionic correlations (4) Jonsson, Bo; Wennerstrom, H.; Halle, B. J . Phys. Chem. 1980, 84, 2179. ( 5 ) Torrie, G. M.; Valleau, J. P. Chem. Phys. Lett. 1979, 65, 343. (6) Guldbrand, L.; JBnsson, Bo,Wennerstrom, H.; Linse, P. J . Chem. Phys. 1984, 80,2221. (7) Kjellander, R.; Marcelja, S.Chem. Phys. Lett. 1984, 112, 49.
0 1990 American Chemical Society
