- Email: [email protected]
Binding of hydrogen bonding solutes at liquid–vapour interfaces of molecular fluids
Chemical Physics Letters 400 (2004) 515–519 www.elsevier.com/locate/cplett
Binding of hydrogen bonding solutes at liquid–vapour interfaces of molecular fluids Sandip Paul, Amalendu Chandra
*
Department of Chemistry, Indian Institute of Technology, Kanpur 208016, India Received 6 June 2004; in final form 12 October 2004 Available online 21 November 2004
Abstract We have calculated the potential of mean force (PMF) for the transfer of a solute molecule across a liquid–vapour interface for four different systems: (a) one methanol molecule in water, (b) one water molecule in methanol, (c) one acetonitrile molecule in water and (d) one water molecule in acetonitrile by means of constrained molecular dynamics simulations. A minimum of the PMF is found near the Gibbs dividing surface for methanol and acetonitrile solutes although the degree of surface activity is found to be somewhat different due, in part, to varying hydrogen bonding nature of these two solutes. 2004 Elsevier B.V. All rights reserved.
1. Introduction The studies of interfacial systems are extremely important in chemistry, biology, materials science and also in environmental and atmospheric sciences. The uptake of pollutants by water clouds is an important atmospheric phenomenon [1] which involves adsorption at aqueous liquid–vapour interfaces. Theoretical studies, including computer simulations, have been carried out extensively to study the liquid–vapour interfaces of pure liquids, mixtures and and ionic solutions [2–13]. Most of the these studies dealt with structural, dynamical and thermodynamic properties. The transport of solute molecules at liquid–vapour interfaces, which is governed by free energy profiles, has received relatively less attention. There are few studies where free energy profiles have been calculated using ions as solutes [14–17]. Among the studies on transport of molecular solutes, we note the work of Dang and Feller [18] who studied the transfer of a benzene molecule across the liquid–vapour interface of water. They found a large free energy minimum *
Corresponding author. Fax: +91 512 2597436/2590007. E-mail address: [email protected] (A. Chandra).
0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.10.144
near the Gibbs dividing surface (GDS) with an adsorption free energy about 4 kcal/mol [18]. In another study, Pohorille and Benjamin [19] carried out molecular dynamics simulation of phenol at the liquid–vapour interface of water. They reported that the adsorption free energy of phenol at the water–vapour interface is 2.8 ± 0.4 kcal/mol [19]. Phenol is a hydrogen bonding solute and it is logical to consider other hydrogen bonding solutes such as methanol or acetonitrile at water liquid–vapour interfaces and see to what extent the binding of these solutes at the interfaces differs from that of phenol. In the present work, we have carried out a series of constrained molecular dynamics simulations to study the binding of a single methanol or an acetonitrile molecule at the liquid–vapour interface of water and also the binding of a water molecule at methanol and acetonitrile liquid–vapour interfaces. The main focus is to calculate the free energy profiles associated with changes in position of the solute molecule as we move across the liquid–vapour interface of the solvent. The outline of the present Letter is as follows. In Section 2, we discuss the molecular models. The simulation details and computational methodology are presented in
516
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
Section 3. Section 4 deals with the results and discussion and our conclusions are summarized in Section 5.
2. Molecular models The water molecules are characterized by the SPC/E potential where each water molecule is characterized by three interaction sites located on oxygen and hydrogen atoms [20]. The acetonitrile molecules are modeled as a rigid three-site molecule comprising a united atom model for the CH3 group, C and N units in a linear ˚ arrangement. The CH3–C bond distance is 1.458 A ˚ and C–N distance is 1.157 A [21]. For methanol molecule, we have used HI potential model where methanol molecules are modeled as rigid three-site molecule [22]. The CH3 group is modeled as a united atom, CH3–O ˚ , O–H distance is 0.945 A ˚ and bond distance is 1.425 A the COH angle is 108.5 [22]. In these models, the interaction between atomic sites of two different molecules is expressed as " 6 # 12 qa qb rab rab uab ðra ; rb Þ ¼ 4ab ; ð1Þ þ rab rab rab where rab is the distance between the atomic sites a and b and qa is the charge of the ath atom. The LennardJones parameters rab and ab are obtained by using the combination rules rab = (ra + rb)/2 and ab ¼ pffiffiffiffiffiffiffiffi a b . The values of the potential parameters qa, ra and a for water, acetonitrile and methanol are shown in Table 1.
3. Simulation details and computational methodology In this study, we used constrained molecular dynamics simulations to calculate the free energy changes associated with the transport of a solute molecule across the liquid–vapour interface of a solvent. The coordinate for solute transport can be considered as the z-coordinate of Table 1 Values of the Lennard-Jones and electrostatic interaction potential parameters of water, acetonitrile and methanol Atom
r ˚) (A
(kcal/mol)
Charge (e)
Water (SPC/E model)
O H
3.166 –
0.1554 –
–0.8476 +0.4238
Acetonitrile
CH3 C N
3.60 3.40 3.30
0.384 0.1 0.1
+0.269 0.129 –0.398
Methanol (H1 model)
CH3 O H
3.861 3.083 –
01823 0.1758 –
+0.297 0.728 +0.431
e represents the magnitude of electronic charge.
the solute molecule in the liquid–vapour interface normal direction. The Helmholtz free energy difference, DA(zs), between a state where the solute molecular centre of mass located at zs and a reference state where the; centre of mass of solute molecule z0, is [15,16] Z zs DAðzs Þ ¼ Aðzs Þ A0 ¼ fs ðz0s Þ dz0s ; ð2Þ z0
fz ðz0s Þ
where is the total z-component force exerted on the solute molecule located at position z0s and A0 is the free energy when the solute molecule is at z0 which was taken to be in the bulk region. The z-coordinate of the centre of mass of the solute molecule was kept fixed during a given simulation and the average zcomponent force acting on the solute was then calcualated. The free energy was then evaluated by integrating the average forces acting on the solute molecule along the z-direction. For each system, we first carried out a bulk simulation in a cubic box of 500 molecules including solvent and solute (499 solvent molecules and one solute molecule), periodically replicated in all three dimensions. We considered four systems: one methanol molecule in water, one water molecule in methanol, one acetonitrile molecule in water and one water molecule in acetonitrile. The box length L was adjusted according to the experimental density of the bulk solutions at 298 K. After this bulk solution was properly equilibrated, two empty boxes of equal size were added on either side of the original simulation box along the z-dimension and this larger rectangular box was taken as the simulation box in the next phase of the simulation run. That means the simulation cell was chosen to be of dimension L · L · Lz with Lz = 3L and the system was reequilibrated by imposing periodic boundary conditions in all three dimensions. This resulted in a lamella of approximate width L separated by vacuum layers of approximate width 2L. Some of the molecules were found to vapourize to the empty space to form liquid– vapour interfaces on both sides of the lamella. In Table 2, we have included, for each system, the number of solvent and solute molecules and also the lengths of the simulation box along x- and z-dimensions. In all simulations, the long range electrostatic interactions were treated by using the three dimensional Ewald method [23]. The short range Lennard-Jones interactions were calculated by using a spherical cut-off at distance L/2. We employed the quaternion formulation of the equations of rotational motion [23] and, for the integration over time, we adapted the leap-frog algorithm with a time step of 1015 s (1 fs). For the calculation of free energy changes for transfer of a solute molecule across the liquid– vapour interface, the position of the solute molecule is ˚ , with changed from z = 0 (middle of the box) to 30 A ˚ in each simulation run. a position increment of 1.0 A The equilibration period of each simulation run was
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
517
Table 2 The number of solvent Nsolvent and solute Nsolute molecules and the lengths of the simulation box along x (Lix) and z (Lz) directions System
Solvent
Solute
Nsolvent
Nsolvent
Lx ˚) (A
Lz ˚) (A
1 2 3 4
Water Methanol Water Acetonitrile
Methanol Water Acetonitrile Water
499 499 499 499
1 1 1 1
24.70 32.31 24.68 35.25
74.10 96.93 74.04 105.75
200 ps and then the average force on the solute was calculated over a production period of 200 ps.
4. Results and discussion In Fig. 1a, we have shown the free energy profile for transfer of a solute methanol molecule across the liquid– vapour interface of water at 298 K. It is seen from the figure that as the methanol molecule begins to approach from bulk liquid to the interface there is practically no ˚ away from change in the free energy up to about 8 A the middle of the bulk liquid. Then the free energy decreases to a minimum that corresponds to a surface active state in which the methanol molecule adsorbed at the liquid–vapour interface of water. The distance of ˚ away from the location this minimum is about 0.5 A of the Gibbs dividing surface (GDS). The free energy then increases as the methanol molecule leaves the interface and enters the gas phase. When the methanol molecule becomes a gas phase molecule the mean force acting on it is essentially zero and so the free energy is
(a)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2
(b)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2 30
24
18
12
6
0
z [Å] Fig. 1. PMF for transferring (a) a methanol molecule across the water liquid/vapour interface and (b) a water molecule across the methanol liquid/vapour interface.
constant. The free energy of adsorption at interface, which is the free energy difference for a methanol molecule being located in the water liquid–vapour interface minimum relative to the gas phase, is about 6.5 kcal/ mol. The free energy difference for a methanol molecule when it is in the interface and in the bulk liquid is about 2 kcal/mol. This clearly shows that methanol molecule is surface active and it is adsorbed at the liquid–vapour interface of water. The adsorption of methanol molecule at the liquid–vapour interface of water–methanol mixture has also been reported earlier by Matsumoto et al. [10] who calculated the density profiles of methanol molecules rather than the free energy profile. The free energy of solvation is defined as the free energy difference between a molecule being solvated in the bulk versus the molecule in the gas phase. Thus, from Fig. 1a, we can also calculate the free energy of solvation. For methanol solute in water solvent this free energy of solvation is 4.47 kcal/mol which agrees well with the experimental value of 5.10 kcal/mol for solvation of a methanol molecule in pure water [24]. We have also calculated the free energy profile for transfer of a water molecule across the liquid–vapour interface of methanol solvent (Fig. 1b) at the same temperature. In contrast to the results of Fig. 1a, we do not observe any minimum at or near GDS in Fig. 1b. Actually, the DA value is higher when solute water molecule is at the interface than that in bulk liquid methanol. Thus, we can infer that water molecule is not adsorbed at the liquid–vapour interface of methanol, i.e., water molecule is not surface active in methanol. The free energy of solvation in liquid methanol and the free energy of adsorption at interface are found to be 5.90 and 5.5 kcal/mol, respectively. When we transfer an acetonitrile solute across the liquid–vapour interface of water (Fig. 2a), again a minimum of free energy occurs near the Gibbs dividing surface is found. In this case the free energy of adsorption at interface is about 5.5 kcal/mol and the free energy difference between a state when the acetonitrile molecule is present in the interface and in the bulk liquid water is about 1 kcal/mol. This indicates that acetonitrile molecule is also surface active. From the free energy calculation for acetonitrile solute in pure water, we have got a value of4.35 kcal/mol for the free energy of solvation. By comparing the results of Figs. 1 and 2a, we can say that methanol
518
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
(a)
4
GDS
∆A (z) [ kcal/mol]
6
2
vapour
liquid
0 –2
(b)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2 30
24
18
12
6
0
o
z[A] Fig. 2. PMF for transferring (a) a acetonitrile molecule across the water liquid/vapour interface and (b) a water molecule across the acetonitrile liquid/vapour interface.
molecule interact more strongly with the liquid–vapour interface of water than an acetonitrile molecule. The methanol molecule shows about 1 kcal/mol more stabilization free energy at the water–vapour interface than that of an acetonitrile molecule. This varying surface activity of methanol and acetonitrile molecules is believed to be partly due to their different hydrogen bonding properties in liquid water [25]. In Fig. 2b, we have shown the free energy profile for transport of a water molecule across the liquid–vapour interface of acetonitrile. As for water in methanol, no free energy minimum is found at or near the GDS and free energy of adsorption at interface is 3.5 kcal/mol. The free energy of solvation for this system is 4.11 kcal/mol. These values are somewhat higher (less negative) than the corresponding values for water in methanol which can be attributed to the stronger hydrogen bonds between water and methanol molecules as compared to that between water and acetonitrile molecules [25].
5. Conclusion We have presented constrained molecular dynamics simulation to study the free energy changes associated with the transfer of a hydrogen bonding solute molecule across the liquid–vapour interface of a solvent. For the transfer of a methanol molecule across liquid–vapour interface of water, there exists a minimum near the GDS with well depth of about 2.0 kcal/mol but there is no such minimum when we consider the transfer of a water molecule across the liquid–vapour interface of
methanol. The acetonitrile solute in water also shows a minimum minimum near GDS. The well depth for this system is about 1.0 kcal/mol. No minimum of free energy is observed for water at the liquid–vapour interface of acetonitrile solvent. We have also calculated the free energy of solvation for all the four systems studied here. For methanol solute in water solvent the calculated free energy of solvation agrees well with the experimental result. We are not aware of similar experimental results of acetonitrile solute in water or of water solute in acetonitrile or methanol solvents. We note that our calculation of the potential of mean force or the work function for transferring a solute from one location to another in the z-direction is based on a relation (Eq. (2)) that has been widely used in the literature [15,16,18,26,27]. Its derivation can be demonstrated by using the thermodynamic integration method [28]. Recently, Mu¨lders et al. [29] showed that although the thermodynamic force, i.e., the derivative of potential of mean force (PMF) with respect to the reaction coordinate, is not exactly equal to the constraint mechanical force for a single configuration, the equality of the two quantities holds when the emsemble averaging is performed so that the difference of the PMF can be computed by integrating the average force computed in simulation as done in the present work. However, a more recent and more rigorous analysis [30] has shown that the conclusion of Mu¨lders et al. [29] holds when the contribution of the Jacobian of the transformation between the generalized and Cartesian coordinates to the PMF is ignored. In the present problem, strictly speaking, such contribution should be included for an exact evaluation of the PMF difference since the reaction coordinate (z) and the Cartesian coordinates are not independent variables. However, the relative contribution of this Jacobian to the PMF (or free energy) difference as compared to the contribution from average force (Eq. (2)) for macroscopic systems like the ones considered here is still not known. We are not aware of any free energy calculations of macroscopic systems where this additional contribution has been explicitly included. The fact that the calculated free energy difference using Eq. (2) has been found to be in good agreement with available experimental results for many systems [15,18,27] seems to indicate that the relative importance of this additional contribution from Jacobian is probably not very significant but it is yet to be verified through rigorous calculations for macroscopic systems. In the present calculations, the water, methanol and acetonitrile molecules are modelled as non-polarizable molecules with fixed partial charges on various atomic sites. The solvent and solute molecules at liquid–vapour interfaces, however, encounter varying environments and hence the degree of polarization of individual molecules can vary significantly across such interfaces. Thus, the polarization effects can play an important role in determining the equilibrium and dynamical properties
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
of liquid–vapour interfaces as has been shown in some of the recent theoretical studies on aqueous ionic solutions [11,12]. We also note in this context two recent studies on the binding of water and ammonia solutes at aqueous surfaces [27,31] which conclude that the exclusion of many-body polarization effects from the potential models does not interfere with the ability to predict the solvation mechanism at the interface. However, inclusion of such many-body effects improves the quantitative agreement with experimental results. Thus, although the precise role of the many-body polarizability effects in the PMF of hydrogen bonding solutes considered here is yet to be investigated, the results of the non-polarizable models presented here are expected to be valid qualitatively. The present study can be extended in several directions, the most important one would be to include the many-body polarization effects discussed above. There are several methods to account for the polarization effects in molecular dynamics simulations such as the fluctuating charge model [32], charge-on-spring model [33,34] and the self-consistent reaction field method [35]. We are currently developing polarizable models for the water, methanol and acetonitrile molecules using the charge-on-spring model [33,34] and the Ewald boundary condition and the issue of many-body polarisation effects will be investigated once the development of these models are completed. It would also be interesting to study the binding of more complex solute molecules such as long chain surfactants and to investigate the solute adsorption at interfaces by using more rigorous methods such the ab initio molecular dynamics. Work in these directions is in progress.
[2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
[23] [24] [25] [26] [27] [28]
Acknowledgements [29]
We gratefully acknowledge the financial support from the Council of Scientific and Industrial Research (CSIR) and Department of Science and Technology, Government of India.
[30] [31]
References
[33] [34]
[1] R.W. Johnson, G.E. Gordon (Eds.), The Chemistry of Acid Rain: Sources and Atmospheric Processes, American Chemical Sympo-
[32]
[35]
519
sium Series, 1987, vol. 93, American Chemical Society, Washington, DC, 1989, p. 2909. M.A. Wilson, A. Pohorille, L.R. Pratt, Chem. Phys. 129 (1989) 209. M. Matsumoto, Y. Kataoka, J. Chem. Phys. 90 (1989) 2398. E.N. Brodskya, A.I. Rusanov, Mol. Phys. 62 (1987) 251. N.I. Christou, J.S. Whitehouse, D. Nicholson, N.G. Parsonage, Mol. Phys. 55 (1985) 397. I. Benjamin, Chem. Rev. 96 (1996) 1449, and references therein. R.S. Taylor, L.X. Dang, B.C. Garret, J. Phys. Chem. 100 (1996) 11720. J. Alejandre, D.J. Tildesley, G.A. Chapela, J. Chem. Phys. 102 (1995) 4574. T.-M. Chang, K.A. Peterson, L.X. Dang, J. Chem. Phys. 103 (1995) 7502. M. Matsumoto, Y. Takaoka, Y. Kataoka, J. Chem. Phys. 98 (1993) 1464. P. Jungwirth, D.J. Tobias, J. Phys. Chem. B 105 (2001) 10468. P. Jungwirth, D.J. Tobias, J. Phys. Chem. B 106 (2002) 6361. S. Paul, A. Chandra, Chem. Phys. Lett. 373 (2003) 87. I. Benjamin, J. Chem. Phys. 95 (1991) 3698. L.X. Dang, J. Phys. Chem. B 106 (2002) 10388. L.X. Dang, J. Chem. Phys. 119 (2003) 6351. M.A. Wilson, A. Pohorille, J. Chem. Phys. 95 (1991) 6005. L.X. Dang, D. Feller, J. Phys. Chem. B 104 (2000) 4403. A. Pohorille, I. Benjamin, J. Chem. Phys. 94 (1991) 5599. H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. V. Molinero, D. Laria, R. Kapral, J. Chem. Phys. 109 (1998) 6844. M. Haughney, M. Ferrario, I.R. Mcdonald, J. Phys. Chem. 91 (1987) 4934; Mol. Phys. 58 (1986) 849. M.P. Alien, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. A. Ben-Naim, Y. Marcus, J. Chem. Phys. 81 (1984) 2016. S. Chowdhuri, A. Chandra, Chem. Phys. Lett. 373 (2003) 79. M. Hayoum, M. Meyer, P. Turq, J. Phys. Chem. 98 (1994) 6626. L.X. Dang, B.C. Garrett, Chem. Phys. Lett. 385 (2004) 309. D. Frenkel, in: G. Ciccotti, W.G. Hoover (Eds.), Molecular Dynamics Simulations of Statistical Mechanical Systems, NorthHolland, Amsterdam, 1986. T. Mu¨lders, P. Kru¨ger, W. Swegat, J. Schlitter, J. Chem. Phys. 104 (1996) 4869. W.K. den Otter, W.J. Briels, J. Chem. Phys. 109 (1998) 4139. R.S. Taylor, D. Ray, B.C. Garrett, J. Phys. Chem. B 101 (1997) 5473. S.W. Rick, S.J. Stuart, B.J. Berne, J. Chem. Phys. 101 (1994) 6141. T.P. Straatsma, J.A. McCammon, Mol. Simulat. 5 (1990) 181. H. Yu, T. Hansson, W.F. van Gunsteren, J. Chem. Phys. 118 (2003) 221. F.J. Luque, J.M. Bofill, M. Orozco, J. Chem. Phys. 103 (1995) 10183.
Binding of hydrogen bonding solutes at liquid–vapour interfaces of molecular fluids Sandip Paul, Amalendu Chandra
*
Department of Chemistry, Indian Institute of Technology, Kanpur 208016, India Received 6 June 2004; in final form 12 October 2004 Available online 21 November 2004
Abstract We have calculated the potential of mean force (PMF) for the transfer of a solute molecule across a liquid–vapour interface for four different systems: (a) one methanol molecule in water, (b) one water molecule in methanol, (c) one acetonitrile molecule in water and (d) one water molecule in acetonitrile by means of constrained molecular dynamics simulations. A minimum of the PMF is found near the Gibbs dividing surface for methanol and acetonitrile solutes although the degree of surface activity is found to be somewhat different due, in part, to varying hydrogen bonding nature of these two solutes. 2004 Elsevier B.V. All rights reserved.
1. Introduction The studies of interfacial systems are extremely important in chemistry, biology, materials science and also in environmental and atmospheric sciences. The uptake of pollutants by water clouds is an important atmospheric phenomenon [1] which involves adsorption at aqueous liquid–vapour interfaces. Theoretical studies, including computer simulations, have been carried out extensively to study the liquid–vapour interfaces of pure liquids, mixtures and and ionic solutions [2–13]. Most of the these studies dealt with structural, dynamical and thermodynamic properties. The transport of solute molecules at liquid–vapour interfaces, which is governed by free energy profiles, has received relatively less attention. There are few studies where free energy profiles have been calculated using ions as solutes [14–17]. Among the studies on transport of molecular solutes, we note the work of Dang and Feller [18] who studied the transfer of a benzene molecule across the liquid–vapour interface of water. They found a large free energy minimum *
Corresponding author. Fax: +91 512 2597436/2590007. E-mail address: [email protected] (A. Chandra).
0009-2614/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.10.144
near the Gibbs dividing surface (GDS) with an adsorption free energy about 4 kcal/mol [18]. In another study, Pohorille and Benjamin [19] carried out molecular dynamics simulation of phenol at the liquid–vapour interface of water. They reported that the adsorption free energy of phenol at the water–vapour interface is 2.8 ± 0.4 kcal/mol [19]. Phenol is a hydrogen bonding solute and it is logical to consider other hydrogen bonding solutes such as methanol or acetonitrile at water liquid–vapour interfaces and see to what extent the binding of these solutes at the interfaces differs from that of phenol. In the present work, we have carried out a series of constrained molecular dynamics simulations to study the binding of a single methanol or an acetonitrile molecule at the liquid–vapour interface of water and also the binding of a water molecule at methanol and acetonitrile liquid–vapour interfaces. The main focus is to calculate the free energy profiles associated with changes in position of the solute molecule as we move across the liquid–vapour interface of the solvent. The outline of the present Letter is as follows. In Section 2, we discuss the molecular models. The simulation details and computational methodology are presented in
516
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
Section 3. Section 4 deals with the results and discussion and our conclusions are summarized in Section 5.
2. Molecular models The water molecules are characterized by the SPC/E potential where each water molecule is characterized by three interaction sites located on oxygen and hydrogen atoms [20]. The acetonitrile molecules are modeled as a rigid three-site molecule comprising a united atom model for the CH3 group, C and N units in a linear ˚ arrangement. The CH3–C bond distance is 1.458 A ˚ and C–N distance is 1.157 A [21]. For methanol molecule, we have used HI potential model where methanol molecules are modeled as rigid three-site molecule [22]. The CH3 group is modeled as a united atom, CH3–O ˚ , O–H distance is 0.945 A ˚ and bond distance is 1.425 A the COH angle is 108.5 [22]. In these models, the interaction between atomic sites of two different molecules is expressed as " 6 # 12 qa qb rab rab uab ðra ; rb Þ ¼ 4ab ; ð1Þ þ rab rab rab where rab is the distance between the atomic sites a and b and qa is the charge of the ath atom. The LennardJones parameters rab and ab are obtained by using the combination rules rab = (ra + rb)/2 and ab ¼ pffiffiffiffiffiffiffiffi a b . The values of the potential parameters qa, ra and a for water, acetonitrile and methanol are shown in Table 1.
3. Simulation details and computational methodology In this study, we used constrained molecular dynamics simulations to calculate the free energy changes associated with the transport of a solute molecule across the liquid–vapour interface of a solvent. The coordinate for solute transport can be considered as the z-coordinate of Table 1 Values of the Lennard-Jones and electrostatic interaction potential parameters of water, acetonitrile and methanol Atom
r ˚) (A
(kcal/mol)
Charge (e)
Water (SPC/E model)
O H
3.166 –
0.1554 –
–0.8476 +0.4238
Acetonitrile
CH3 C N
3.60 3.40 3.30
0.384 0.1 0.1
+0.269 0.129 –0.398
Methanol (H1 model)
CH3 O H
3.861 3.083 –
01823 0.1758 –
+0.297 0.728 +0.431
e represents the magnitude of electronic charge.
the solute molecule in the liquid–vapour interface normal direction. The Helmholtz free energy difference, DA(zs), between a state where the solute molecular centre of mass located at zs and a reference state where the; centre of mass of solute molecule z0, is [15,16] Z zs DAðzs Þ ¼ Aðzs Þ A0 ¼ fs ðz0s Þ dz0s ; ð2Þ z0
fz ðz0s Þ
where is the total z-component force exerted on the solute molecule located at position z0s and A0 is the free energy when the solute molecule is at z0 which was taken to be in the bulk region. The z-coordinate of the centre of mass of the solute molecule was kept fixed during a given simulation and the average zcomponent force acting on the solute was then calcualated. The free energy was then evaluated by integrating the average forces acting on the solute molecule along the z-direction. For each system, we first carried out a bulk simulation in a cubic box of 500 molecules including solvent and solute (499 solvent molecules and one solute molecule), periodically replicated in all three dimensions. We considered four systems: one methanol molecule in water, one water molecule in methanol, one acetonitrile molecule in water and one water molecule in acetonitrile. The box length L was adjusted according to the experimental density of the bulk solutions at 298 K. After this bulk solution was properly equilibrated, two empty boxes of equal size were added on either side of the original simulation box along the z-dimension and this larger rectangular box was taken as the simulation box in the next phase of the simulation run. That means the simulation cell was chosen to be of dimension L · L · Lz with Lz = 3L and the system was reequilibrated by imposing periodic boundary conditions in all three dimensions. This resulted in a lamella of approximate width L separated by vacuum layers of approximate width 2L. Some of the molecules were found to vapourize to the empty space to form liquid– vapour interfaces on both sides of the lamella. In Table 2, we have included, for each system, the number of solvent and solute molecules and also the lengths of the simulation box along x- and z-dimensions. In all simulations, the long range electrostatic interactions were treated by using the three dimensional Ewald method [23]. The short range Lennard-Jones interactions were calculated by using a spherical cut-off at distance L/2. We employed the quaternion formulation of the equations of rotational motion [23] and, for the integration over time, we adapted the leap-frog algorithm with a time step of 1015 s (1 fs). For the calculation of free energy changes for transfer of a solute molecule across the liquid– vapour interface, the position of the solute molecule is ˚ , with changed from z = 0 (middle of the box) to 30 A ˚ in each simulation run. a position increment of 1.0 A The equilibration period of each simulation run was
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
517
Table 2 The number of solvent Nsolvent and solute Nsolute molecules and the lengths of the simulation box along x (Lix) and z (Lz) directions System
Solvent
Solute
Nsolvent
Nsolvent
Lx ˚) (A
Lz ˚) (A
1 2 3 4
Water Methanol Water Acetonitrile
Methanol Water Acetonitrile Water
499 499 499 499
1 1 1 1
24.70 32.31 24.68 35.25
74.10 96.93 74.04 105.75
200 ps and then the average force on the solute was calculated over a production period of 200 ps.
4. Results and discussion In Fig. 1a, we have shown the free energy profile for transfer of a solute methanol molecule across the liquid– vapour interface of water at 298 K. It is seen from the figure that as the methanol molecule begins to approach from bulk liquid to the interface there is practically no ˚ away from change in the free energy up to about 8 A the middle of the bulk liquid. Then the free energy decreases to a minimum that corresponds to a surface active state in which the methanol molecule adsorbed at the liquid–vapour interface of water. The distance of ˚ away from the location this minimum is about 0.5 A of the Gibbs dividing surface (GDS). The free energy then increases as the methanol molecule leaves the interface and enters the gas phase. When the methanol molecule becomes a gas phase molecule the mean force acting on it is essentially zero and so the free energy is
(a)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2
(b)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2 30
24
18
12
6
0
z [Å] Fig. 1. PMF for transferring (a) a methanol molecule across the water liquid/vapour interface and (b) a water molecule across the methanol liquid/vapour interface.
constant. The free energy of adsorption at interface, which is the free energy difference for a methanol molecule being located in the water liquid–vapour interface minimum relative to the gas phase, is about 6.5 kcal/ mol. The free energy difference for a methanol molecule when it is in the interface and in the bulk liquid is about 2 kcal/mol. This clearly shows that methanol molecule is surface active and it is adsorbed at the liquid–vapour interface of water. The adsorption of methanol molecule at the liquid–vapour interface of water–methanol mixture has also been reported earlier by Matsumoto et al. [10] who calculated the density profiles of methanol molecules rather than the free energy profile. The free energy of solvation is defined as the free energy difference between a molecule being solvated in the bulk versus the molecule in the gas phase. Thus, from Fig. 1a, we can also calculate the free energy of solvation. For methanol solute in water solvent this free energy of solvation is 4.47 kcal/mol which agrees well with the experimental value of 5.10 kcal/mol for solvation of a methanol molecule in pure water [24]. We have also calculated the free energy profile for transfer of a water molecule across the liquid–vapour interface of methanol solvent (Fig. 1b) at the same temperature. In contrast to the results of Fig. 1a, we do not observe any minimum at or near GDS in Fig. 1b. Actually, the DA value is higher when solute water molecule is at the interface than that in bulk liquid methanol. Thus, we can infer that water molecule is not adsorbed at the liquid–vapour interface of methanol, i.e., water molecule is not surface active in methanol. The free energy of solvation in liquid methanol and the free energy of adsorption at interface are found to be 5.90 and 5.5 kcal/mol, respectively. When we transfer an acetonitrile solute across the liquid–vapour interface of water (Fig. 2a), again a minimum of free energy occurs near the Gibbs dividing surface is found. In this case the free energy of adsorption at interface is about 5.5 kcal/mol and the free energy difference between a state when the acetonitrile molecule is present in the interface and in the bulk liquid water is about 1 kcal/mol. This indicates that acetonitrile molecule is also surface active. From the free energy calculation for acetonitrile solute in pure water, we have got a value of4.35 kcal/mol for the free energy of solvation. By comparing the results of Figs. 1 and 2a, we can say that methanol
518
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
(a)
4
GDS
∆A (z) [ kcal/mol]
6
2
vapour
liquid
0 –2
(b)
4
GDS
∆A (z) [ kcal/mol ]
6
2
vapour
liquid
0 –2 30
24
18
12
6
0
o
z[A] Fig. 2. PMF for transferring (a) a acetonitrile molecule across the water liquid/vapour interface and (b) a water molecule across the acetonitrile liquid/vapour interface.
molecule interact more strongly with the liquid–vapour interface of water than an acetonitrile molecule. The methanol molecule shows about 1 kcal/mol more stabilization free energy at the water–vapour interface than that of an acetonitrile molecule. This varying surface activity of methanol and acetonitrile molecules is believed to be partly due to their different hydrogen bonding properties in liquid water [25]. In Fig. 2b, we have shown the free energy profile for transport of a water molecule across the liquid–vapour interface of acetonitrile. As for water in methanol, no free energy minimum is found at or near the GDS and free energy of adsorption at interface is 3.5 kcal/mol. The free energy of solvation for this system is 4.11 kcal/mol. These values are somewhat higher (less negative) than the corresponding values for water in methanol which can be attributed to the stronger hydrogen bonds between water and methanol molecules as compared to that between water and acetonitrile molecules [25].
5. Conclusion We have presented constrained molecular dynamics simulation to study the free energy changes associated with the transfer of a hydrogen bonding solute molecule across the liquid–vapour interface of a solvent. For the transfer of a methanol molecule across liquid–vapour interface of water, there exists a minimum near the GDS with well depth of about 2.0 kcal/mol but there is no such minimum when we consider the transfer of a water molecule across the liquid–vapour interface of
methanol. The acetonitrile solute in water also shows a minimum minimum near GDS. The well depth for this system is about 1.0 kcal/mol. No minimum of free energy is observed for water at the liquid–vapour interface of acetonitrile solvent. We have also calculated the free energy of solvation for all the four systems studied here. For methanol solute in water solvent the calculated free energy of solvation agrees well with the experimental result. We are not aware of similar experimental results of acetonitrile solute in water or of water solute in acetonitrile or methanol solvents. We note that our calculation of the potential of mean force or the work function for transferring a solute from one location to another in the z-direction is based on a relation (Eq. (2)) that has been widely used in the literature [15,16,18,26,27]. Its derivation can be demonstrated by using the thermodynamic integration method [28]. Recently, Mu¨lders et al. [29] showed that although the thermodynamic force, i.e., the derivative of potential of mean force (PMF) with respect to the reaction coordinate, is not exactly equal to the constraint mechanical force for a single configuration, the equality of the two quantities holds when the emsemble averaging is performed so that the difference of the PMF can be computed by integrating the average force computed in simulation as done in the present work. However, a more recent and more rigorous analysis [30] has shown that the conclusion of Mu¨lders et al. [29] holds when the contribution of the Jacobian of the transformation between the generalized and Cartesian coordinates to the PMF is ignored. In the present problem, strictly speaking, such contribution should be included for an exact evaluation of the PMF difference since the reaction coordinate (z) and the Cartesian coordinates are not independent variables. However, the relative contribution of this Jacobian to the PMF (or free energy) difference as compared to the contribution from average force (Eq. (2)) for macroscopic systems like the ones considered here is still not known. We are not aware of any free energy calculations of macroscopic systems where this additional contribution has been explicitly included. The fact that the calculated free energy difference using Eq. (2) has been found to be in good agreement with available experimental results for many systems [15,18,27] seems to indicate that the relative importance of this additional contribution from Jacobian is probably not very significant but it is yet to be verified through rigorous calculations for macroscopic systems. In the present calculations, the water, methanol and acetonitrile molecules are modelled as non-polarizable molecules with fixed partial charges on various atomic sites. The solvent and solute molecules at liquid–vapour interfaces, however, encounter varying environments and hence the degree of polarization of individual molecules can vary significantly across such interfaces. Thus, the polarization effects can play an important role in determining the equilibrium and dynamical properties
S. Paul, A. Chandra / Chemical Physics Letters 400 (2004) 515–519
of liquid–vapour interfaces as has been shown in some of the recent theoretical studies on aqueous ionic solutions [11,12]. We also note in this context two recent studies on the binding of water and ammonia solutes at aqueous surfaces [27,31] which conclude that the exclusion of many-body polarization effects from the potential models does not interfere with the ability to predict the solvation mechanism at the interface. However, inclusion of such many-body effects improves the quantitative agreement with experimental results. Thus, although the precise role of the many-body polarizability effects in the PMF of hydrogen bonding solutes considered here is yet to be investigated, the results of the non-polarizable models presented here are expected to be valid qualitatively. The present study can be extended in several directions, the most important one would be to include the many-body polarization effects discussed above. There are several methods to account for the polarization effects in molecular dynamics simulations such as the fluctuating charge model [32], charge-on-spring model [33,34] and the self-consistent reaction field method [35]. We are currently developing polarizable models for the water, methanol and acetonitrile molecules using the charge-on-spring model [33,34] and the Ewald boundary condition and the issue of many-body polarisation effects will be investigated once the development of these models are completed. It would also be interesting to study the binding of more complex solute molecules such as long chain surfactants and to investigate the solute adsorption at interfaces by using more rigorous methods such the ab initio molecular dynamics. Work in these directions is in progress.
[2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
[23] [24] [25] [26] [27] [28]
Acknowledgements [29]
We gratefully acknowledge the financial support from the Council of Scientific and Industrial Research (CSIR) and Department of Science and Technology, Government of India.
[30] [31]
References
[33] [34]
[1] R.W. Johnson, G.E. Gordon (Eds.), The Chemistry of Acid Rain: Sources and Atmospheric Processes, American Chemical Sympo-
[32]
[35]
519
sium Series, 1987, vol. 93, American Chemical Society, Washington, DC, 1989, p. 2909. M.A. Wilson, A. Pohorille, L.R. Pratt, Chem. Phys. 129 (1989) 209. M. Matsumoto, Y. Kataoka, J. Chem. Phys. 90 (1989) 2398. E.N. Brodskya, A.I. Rusanov, Mol. Phys. 62 (1987) 251. N.I. Christou, J.S. Whitehouse, D. Nicholson, N.G. Parsonage, Mol. Phys. 55 (1985) 397. I. Benjamin, Chem. Rev. 96 (1996) 1449, and references therein. R.S. Taylor, L.X. Dang, B.C. Garret, J. Phys. Chem. 100 (1996) 11720. J. Alejandre, D.J. Tildesley, G.A. Chapela, J. Chem. Phys. 102 (1995) 4574. T.-M. Chang, K.A. Peterson, L.X. Dang, J. Chem. Phys. 103 (1995) 7502. M. Matsumoto, Y. Takaoka, Y. Kataoka, J. Chem. Phys. 98 (1993) 1464. P. Jungwirth, D.J. Tobias, J. Phys. Chem. B 105 (2001) 10468. P. Jungwirth, D.J. Tobias, J. Phys. Chem. B 106 (2002) 6361. S. Paul, A. Chandra, Chem. Phys. Lett. 373 (2003) 87. I. Benjamin, J. Chem. Phys. 95 (1991) 3698. L.X. Dang, J. Phys. Chem. B 106 (2002) 10388. L.X. Dang, J. Chem. Phys. 119 (2003) 6351. M.A. Wilson, A. Pohorille, J. Chem. Phys. 95 (1991) 6005. L.X. Dang, D. Feller, J. Phys. Chem. B 104 (2000) 4403. A. Pohorille, I. Benjamin, J. Chem. Phys. 94 (1991) 5599. H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. V. Molinero, D. Laria, R. Kapral, J. Chem. Phys. 109 (1998) 6844. M. Haughney, M. Ferrario, I.R. Mcdonald, J. Phys. Chem. 91 (1987) 4934; Mol. Phys. 58 (1986) 849. M.P. Alien, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. A. Ben-Naim, Y. Marcus, J. Chem. Phys. 81 (1984) 2016. S. Chowdhuri, A. Chandra, Chem. Phys. Lett. 373 (2003) 79. M. Hayoum, M. Meyer, P. Turq, J. Phys. Chem. 98 (1994) 6626. L.X. Dang, B.C. Garrett, Chem. Phys. Lett. 385 (2004) 309. D. Frenkel, in: G. Ciccotti, W.G. Hoover (Eds.), Molecular Dynamics Simulations of Statistical Mechanical Systems, NorthHolland, Amsterdam, 1986. T. Mu¨lders, P. Kru¨ger, W. Swegat, J. Schlitter, J. Chem. Phys. 104 (1996) 4869. W.K. den Otter, W.J. Briels, J. Chem. Phys. 109 (1998) 4139. R.S. Taylor, D. Ray, B.C. Garrett, J. Phys. Chem. B 101 (1997) 5473. S.W. Rick, S.J. Stuart, B.J. Berne, J. Chem. Phys. 101 (1994) 6141. T.P. Straatsma, J.A. McCammon, Mol. Simulat. 5 (1990) 181. H. Yu, T. Hansson, W.F. van Gunsteren, J. Chem. Phys. 118 (2003) 221. F.J. Luque, J.M. Bofill, M. Orozco, J. Chem. Phys. 103 (1995) 10183.