Theory and Computer Simulation of Solute Effects on the Surface

Phone: 785-532-5109. Fax: 785-532-6666. ... However, the analysis of computer simulation data on surface adsorption is rather new. We will briefly der...
0 downloads 0 Views 685KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 8975–8984

8975

Theory and Computer Simulation of Solute Effects on the Surface Tension of Liquids Feng Chen and Paul E. Smith* Department of Chemistry, 111 Willard Hall, Kansas State UniVersity, Manhattan, Kansas 66506-3701 ReceiVed: NoVember 20, 2007; ReVised Manuscript ReceiVed: April 15, 2008

A complete description of the thermodynamics of planar mixed solute-solvent interfaces suitable for the analysis of computer simulation data is provided. The approach uses surface probability distributions to characterize the interface regions, coupled with radial distribution functions and the Kirkwood-Buff theory of solutions to characterize the bulk solution properties. The approach is then used to understand the relationship between changes in the surface tension, the degree of surface adsorption or depletion, and the bulk solution properties of two aqueous solute systems. The first, aqueous NaCl solutions, provides an example of a surface excluded solute. The second, aqueous methanol solutions, provides an example of a surface adsorbed solute. The numerical results support the theoretical relationships described here and provide a consistent picture of the thermodynamics of solution interfaces involving any number of components which can be applied to a wide variety of systems. Introduction The distribution of solutes at the solution/vapor interface has important consequences in chemistry and chemical engineering. Consequently, a variety of experimental and theoretical approaches have been used to help understand surface adsorption or exclusion.1–5 It is well established that an increase in the surface tension of a solution due to the addition of a solute indicates exclusion of that solute from the interface region, and vice versa.6 However, atomic level detail concerning the surface distributions has been more difficult to achieve. In principle, molecular simulation can provide such detail. Simulations require an accurate force field and an adequate degree of sampling of the property of interest if the results are to be interpreted with any confidence. Recently, a wealth of information has been provided by simulations of solutes at the water vapor (or vacuum) interface.7–18 We will not discuss all the results here as they have been summarized in several recent reviews.19–22 The main issues arising from a simulation perspective are the observed dependence of the surface structure on the force field used,8,14 and the long times required to obtain precise values for the surface tension and surface structure probability distributions.8,23 In particular, the use of polarizable force fields has generally indicated an increased probability for the location of larger, more polarizable, anions at the surface.8 However, the surface adsorption must be accompanied by a sub surface depletion in order for the Gibbs adsorption equation to be obeyed.15,19 Unfortunately, it has been difficult to fully quantify such structural and thermodynamic changes so that they may be compared with experimental data on surface adsorption. This is a major focus of the current study. Background and Theory Thermodynamics of Surfaces. The general theory of surface adsorption has been established for over a century. However, the analysis of computer simulation data on surface adsorption is rather new. We will briefly derive the Gibbs adsorption isotherm using a simple approach which can be directly related * To whom correspondence should be addressed. Phone: 785-532-5109. Fax: 785-532-6666. E-mail: [email protected].

to the Kirkwood-Buff (KB) theory of solutions, which we will utilize later. KB theory has been used previously to analyze experimental surface adsorption data,24 albeit in a slightly different form. To our knowledge it has not been used to analyze simulation data on the thermodynamics of surface tension changes due to additives. Throughout this article we shall refer to a mixture of a solvent (1), usually water, and any number of additional solutes (2, 3, 4, etc.). Consider the planar surface region of a solution in contact with a vacuum or vapor region containing a negligible number of molecules on one side and bulk solution on the other. The surface lies in the xy plane and is perpendicular to the z axis. The Gibbs-Duhem equations for the surface and the bulk regions at constant T and P are given by

Adγ + N1dµ1 + N2dµ2 + N3dµ3 + · · · ) 0

surface (1a)

n1dµ1 + n2dµ2 + n3dµ3 + · · · ) 0

bulk

(1b)

respectively. Here, A is the surface area of the interface, γ is the surface tension, and µi is the chemical potential of species i. The Ni values represent the number of solvent and solute molecules in the interface region, while the ni values represent the corresponding number in the bulk solution. At equilibrium, the chemical potentials of each species are the same in the surface and bulk regions. The number of molecules in the surface region differs from that in the bulk due to the perturbing effect of the interface which requires a redistribution of solution components in order to maintain equilibrium with the bulk solution. One can use the bulk equation to eliminate the solvent chemical potential change (dµ1) from the surface equation to give

(

Adγ ) - N2 -

) (

)

n2 n3 N dµ - N3 - N1 dµ3 - · · · (2) n1 1 2 n1

or,

dγ ) -Γ2,1dµ2 - Γ3,1dµ3 - · · ·

(3)

where Γi,1 is the usual Gibbs excess surface adsorption (per unit surface area) of each solute i relative to the solvent (1). Taking

10.1021/jp711062a CCC: $40.75  2008 American Chemical Society Published on Web 07/09/2008

8976 J. Phys. Chem. B, Vol. 112, No. 30, 2008

Chen and Smith

derivatives of eq 3 with respect to the number density (Fi ) ni/V) of the solute one finds

( ) ∂γ ∂F2

) -Γ2,1 T,P

( ) ∂µ2 ∂F2

- Γ3,1 T,P

( ) ∂µ3 ∂F2

-· · ·

(4)

∫0Z gi(z)dz

(5)

The distance Z defines the extent of the surface region from some (arbitrary) origin in the vacuum or vapor phase. Beyond Z, the solute and solvent distributions are the same as that in the bulk solution. The excess surface adsorption of solute i per unit surface area is then given by

Γi,1 )

(

)

ni 1 N - N ) Fi A i n1 1

∫0Z [gi(z) - g1(z)]dz ) FiIi,1

(6)

and the integration limit can be formally extended to Z ) ∞, i.e., one does not need to know the extent of the surface region. The integral (Ii,1) over the surface probability distributions can be viewed as a measure of the surface structure. The change in surface tension and the excess adsorption of a species (as defined by eq 6) can be obtained relatively easily from computer simulations.26,27 Even so, we note that no such integration has been performed in recent simulation studies of non volatile solutes.8,10,20,28 The derivatives on the rhs of eq 4 are seemingly more difficult to obtain and have not traditionally been used to analyze computer simulation data on surface adsorption. However, we note that they are properties of the bulk solution. KB Theory of Solutions. A central aspect of the current approach involves the application of KB theory to relate the simulation data obtained for bulk solution properties to the corresponding activity derivatives.29,30 Previous studies have used KB theory to investigate a variety of experimental and theoretical solvation effects.31–33 It is an exact theory of solutions. In particular, it is important to emphasize that KB theory does not involve any approximations or limitations concerning the size or character of the molecules to which it can be applied.32 The theory relates several properties of solution mixtures (containing any number of components) to KB integrals which are defined by29

Gij ) Gji ) 4π

∫

∞

( ) ( )

T,P

which is the same result that has been derived previously.6,25 An expression for Ni can be obtained from the (onedimensional) surface probability distribution functions gi(z),

Ni ) AFi

One property of solutions that will prove particularly useful during the present discussion is the change in solute activity (a2) with solute concentration in a binary solution as defined by

∫

a22 ) β

∂µ2 ∂ln F2

)

T,P

∂ln a2 ∂ln F2

) T,P

1 (8) 1 + F2(G22 - G12)

where β ) 1/RT, and R is the gas constant. To enable the transformation between different concentration scales we will also require the KB expression for the partial molar volume (pmv) of the solvent,29

Vj1 )

1 + F2(G22 - G12) 1 ) η ηa22

(9)

where η ) F1 + F2 + F1F2∆G12 and ∆G12 ) G11 + G22 2G12. All the integrals in eqs 8 and 9 can be obtained from simulations performed at the appropriate bulk solution compositions. The application of KB theory to neutral salt solutions is somewhat more involved.40 We will use the indistinguishable ion approach where the solute number density in eqs 8 and 9 refers to the total ion concentration, and the KB integrals are determined by ignoring the identity of the individual ions.37,40 The consequences of this approach for determining the surface excess or deficit of salts, as provided by eq 6, will be discussed later. Hence, we shall distinguish between the normal molar (or molal) salt concentration (cs or ms), and the total ion concentration (F2 ) n(cs or m2 ) n(ms) for a salt that releases a total of n( ) n+ + n- ions in solution. Therefore, we also have dµs ) n(dµ2, d ln cs ) d ln F2, and d ln a2 ) d ln(y2F2) ) d ln(y(cs), where y is the molar activity coefficient. More details can be found in the literature.37,41,42 Combined Approach for Binary Systems. Here we combine the usual surface adsorption approach with the results provided by the KB theory of solutions to provide a consistent picture of the thermodynamics of surfaces in terms of distribution functions in binary systems-all of which can be obtained from computer simulations. The surface probability distribution functions are normalized by reference to the corresponding bulk distribution. Therefore, one can express the change in surface tension with solute concentration as

( )

β

∂γ ∂F2

) -a22 T,P

Γ2,1 ) -a22 F2

∫0Z [g2(z) - g1(z)]dz ) -a22I2,1

R [gijµVT(r) - 1]r2dr ≈ 4π 0 [gijNPT(r) - 1]r2dr 0

(10)

(7)

where a22 is determined for the corresponding bulk solution, and is usually positive for real solutions (see later). We have used a22 rather than ∂µ2/∂F2 as this removes the inherent singularity in the latter at low solute concentrations. Hence, if the change in surface tension is proportional to the solute concentration, as observed for most salts,43 then the changes in the surface distribution of solute and solvent molecules mirrors the changes in a22. Equation 10 is exact for non volatile solvents and solutes, and all contributions can be obtained from computer simulations. Therefore, an accurate simulated value of a22 will be important for a correct description of surface adsorption. Unfortunately, many common force fields we have tested do not (a priori) reproduce the experimental KB integrals and a22 values for binary solutions.44–48 We note that eq 10 is directly analogous to equations derived for the treatment of solute effects on peptides and proteins.37,49,50 In many cases it is more convenient to express the surface tension changes using other solute concentration scales. To

where gij(r) is the radial distribution function (rdf) in the grand canonical (:VT) ensemble between species i and j as a function of their intermolecular separation r. The negligible approximation in eq 7 is used to enable the determination of KB integrals in closed (NPT) systems,31,34,35 as this is the more typical system to be studied. In this case, the integral is truncated at a distance R beyond which the rdfs are essentially unity. KB integrals can be obtained from experimental data (densities, compressibilities and activities) on solution mixtures,30 or directly from computer simulations.33,36,37 In principle, expressions for the required chemical potential derivatives can be obtained for any number of solute components.29,31 In practice, the expressions become rather cumbersome as the number of components increases and so we will focus on the expressions provided for a binary system of a single solute (2) in a solvent (1). Expressions for ternary systems are available.32,38,39

Solute Effects on the Surface Tension of Liquids

J. Phys. Chem. B, Vol. 112, No. 30, 2008 8977

achieve this one can use the following standard thermodynamic transformations between the molar and molal (m), or mole fraction (x), concentrations in binary solutions,

( ) ∂F2 ∂m2

) F21Vj1 ) F1φ1 and T,P

( ) ∂F2 ∂x2

) F2Vj1

(11)

T,P

j 1 is the partial where φ1 is the volume fraction of the solvent, V molar volume of the solvent, and F ) F1 + F2 is the total number density. Consequently, we have the expression,

( )

β

∂γ ∂x2

)T,P

F2 Γ2,1 F Γ2,1 )η F2 η x2

(12)

for the mole fraction derivative and,

( )

∂γ β ∂m2

F21 Γ2,1 F1 Γ2,1 ))η F η m2 T,P 2

βbcs β∆γ )a22 a22

( )

Γ2,1 ) -βx2

(13)

for the molality derivative. We note that the bulk solute molality is given by m2 ) n2/n1 to within a conversion factor of 1000/ M1 ) 55.51 mol/Kg for water. Analysis of the stability criteria for binary solution mixtures indicates that the value of η must be positive for miscible solutions.31 Consequently, the value of a22 will also be positive for compositions where the pmv of the solvent is positive (see eq 9). This is the case for the majority of water mixtures. Hence, the sign of the surface excess in eqs 10, 12, and 13 is given by the slope of the surface tension against solute concentration using any concentration scale. The conditions for ideal bulk solution behavior depend on the concentration scale being used. On the molar concentration scale an ideal solution (a2 ) F2) is characterized by G22 - G12 ) 0 (or a22 ) 1) for all compositions. Alternatively, on the mole fraction scale one has ∆G12 ) 0 (or η ) F), while ideal behavior on the molality scale is provided by F1∆G12 ) -1 (or η ) F1). The above expressions have several consequences for the analysis of computer simulation data which, to our knowledge, have not been considered in previous studies. First, to fully and accurately describe the thermodynamics and structure of surfaces (as provided by gi) the force fields used must reproduce the change in surface tension with solute concentration, the values of gi (distribution of solute and solvent at the surface), and a22 (or η) a property of the bulk solution. Actually, the above equation is exact and therefore only two of the three need be correct. Second, the integration over g2 - g1 should be performed to a distance at which the integral remains unchanged, i.e., g2 ) g1 ) 1, which may involve many solvation shells away from the surface. The need to include the exact value of a22 or η in experimental studies has been noted previously.51 Some Simple Cases. Here we will investigate several simplified expressions that are obtained in a few limiting cases. These situations provide some useful reference points for real solutions. The first arises when the surface tension change is proportional to the solute molarity - which is common for many salt solutions.43 Therefore, one has ∆γ ) γ - γ1 ) bcs, where γ1 is the surface tension of the pure solvent and b ) (∂γ/∂cs)T,P. Consequently, the surface excess can be obtained from eq 10 which reduces to

Γ2,1 ) -

to the salt (µs ) n+µ+ + n-µ-), whereas the excess is in terms of ions if the indistinguishable ion approximation is used (as is adopted here). The latter is equivalent to using the mean ionic activity coefficient and the total salt concentration in determining a22. The above expression can be reduced further if one assumes or observes ideal behavior for the bulk solution mixture (a22 ) 1). The second type of simplification arises for symmetric ideal (SI) solutions. SI solutions have been discussed in detail by Ben-Naim in the context of KB theory.32 Here, the mole fraction activity coefficients are unity for both species at all compositions in a binary solution. In this case one has ∆G12 ) G11 + G22 2G12 ) 0 and therefore η ) F for all compositions. Therefore, eq 12 reduces to

(14)

The above expression is also valid for salt solutions. In this case the surface excess refers to salt “molecules” if the chemical potential (activity) used in the definition of a22 refers

∂γ ∂x2

(15) T,P

One can also further assume that the surface tension of an SI mixture of 1 and 2 is given by γ ) x1γ1 + x2γ2 where γi is the surface tension of pure i at the same T and P. The above equation then simplifies to give

Γ2,1 ) -βx2(γ2 - γ1)

(16)

which provides another useful reference point for solution mixtures. We will label this type of system as SI2. It is clear that SI2 solution mixtures still give rise to an excess surface adsorption unless the pure solutions also have similar surface tensions. Finally, we note that the limiting slope obtained as the solute concentration tends to zero is provided by

(

β

∂γ ∂ ln y2

)

) -Γo2,1

(17)

T,P,y2f0

with y ) F, x, or m, and where the zero superscript indicates a property of the system at infinite dilution of the solute. Some interesting cases appear in the literature and provide information concerning the structure of the interface at very low solute concentrations as defined by the limiting value of the integral I2,1. For salts displaying a linear dependence of the surface tension on the salt concentration one finds

Io2,1 ) -

βb n(

(18)

In this case the integral refers to the average probability distribution for the ions, i.e., g2 ) (n+g+ + n-g-)/n(. Alternatively, for surface adsorbed solutes an empirical relationship has been developed for low (x2 < 0.01) solute concentrations,52,53

γ ) γ1[1 - B log10(1 + x2 ⁄ C)]

(19)

where B and C are positive constants. Using this expression one can show that at infinite dilution of the solute we have

Io2,1 )

βγ1 B 2.303Fo C

(20)

1

and F1° is the number density of pure water. Hence, knowledge of the parameter b, or B and C, can provide quantitative estimates for the nature of the surface structure. Methods All molecular dynamics simulations were performed using the KBFF force fields for NaCl and methanol,46,54 together with

8978 J. Phys. Chem. B, Vol. 112, No. 30, 2008

Chen and Smith

the SPC/E water model,55 as implemented in the GROMACS program (v3.3.1).56,57 A time-step of 2 fs was used and the methanol bond lengths were constrained using Lincs,58 while the water geometry was constrained using SETTLE.59 A twin range cutoff of 1.0 nm/1.5 nm was employed with a nonbonded pair list update of every 10 steps. Long range electrostatic interactions were evaluated using the PME approach.60 Simulations performed to determine the surface tension and surface adsorption involved approximately 12 000 atoms in a fixed rectangular box of dimensions 4 × 4 × 12 nm3 with the atoms located in a slab geometry occupying the central 4 × 4 × 8 nm3 and generating two surfaces of area A ) 16 nm2. The slab simulations were performed at constant volume and a constant temperature of 300 K using the weak coupling approach.61 The surface tensions were then determined from the pressure tensor elements via the following relationship,

1 1 γ ) Lz Pzz - (Pxx - Pyy) 2 2

[

]

n2 no2 - N2(Z) N1(Z) N1(Z) ) N2(Z) - o n1 n - N (Z) 1

Γ2,1 ) Γ+,1 + Γ-,1 )

(

) (

)

n+ n1 1 N+ - N1 + N- - N1 A n1 A n1

(23) or in terms of the probability distribution functions,

Γ2,1 ) F+

∫0Z [g+(z) - g1(z)]dz + F-∫0Z [g-(z) - g1(z)]dz (24)

and the value of F2 is taken as F+ + F- in eq 6. In all cases the results represent an average over both surfaces within the system, and no evaporation was observed for any of the systems. Preliminary investigations with larger systems (8 × 8 × 12 nm3) do not indicate a significant system size effect on the probability distributions.

(21)

as described previously.62,63 Here, Lz ) 12 nm is the box dimension in the extended (z) direction and PRR are the diagonal elements of the pressure tensor. No long-range dispersion corrections to the surface tension values were included as this is nontrivial for mixed solutions. Fortunately, the contribution to the surface tension of the neglected dispersion interactions beyond 1.5 nm is usually negligible, especially as we are focusing on changes in surface tension. Typically, systems were equilibrated for 5 ns and then simulate for a further 15-30 ns to ensure precise results for both the surface tension and surface adsorption data. The values of a22 for the bulk solutions were determined previously during the initial parametrization procedure.46,54 Both the simulated and experimental data for bulk solutions were taken directly from these earlier studies. The excess adsorption values were determined from the simulations using a simple counting procedure. One can define a local number of solute and solvent molecules which depends on the distance from our arbitrary origin in the vacuum phase, i.e., Ni(Z). The excess adsorption can then be determined as a function of Z. When Z is large enough to approach the bulk solution distribution, the value of Γ2,1 will tend to the constant value required for the current analysis. In order to evaluate this precisely we have used the following expression for the excess adsorption,64

AΓ2,1(Z) ) N2(Z) -

eq 6 becomes

1

(22) where ni° are the initial number of i ions or molecules in the system. The values of ni° and ni differ due to the redistribution of molecules to and from the surface in our finite size system. Therefore, the final bulk solution composition n2/n1 is somewhat different to the initial composition n2°/n1°, and this difference can be significant as N1 can be large. An alternative viewpoint is that any ions or molecules which contribute to the value of Ni should not be considered part of the bulk solution. Using eq 22 one does not have to generate the corresponding rdfs in advance. This can be a problem as the normalization procedure requires the bulk solution concentrations which are unknown if one does not know the extent of the surface region, i.e., Z in eq 5. The bulk solute and solvent concentrations were therefore determined, after equilibration, by averaging over the central 4 nm of the slab. The total excess surface adsorption for salts involves both the distribution of the cations and anions. Hence,

Results To illustrate the relationship between the different terms in eqs 10 and 12 we have performed molecular dynamics simulations to investigate the effect of two different solutes on the surface thermodynamics of water solutions. The first solute was NaCl which produces an increase in the surface tension of the solution from that of pure water and is therefore excluded from the surface. The second solute was methanol which decreases the surface tension and is therefore adsorbed at the surface. In studying these two systems our primary objective was to illustrate and establish the relationships provided by eqs 10 and 12. We have also compared the results to the corresponding experimental data. However, it should be noted that the agreement with experiment is a function of the quality of force field being used. Being exact for non volatile solvents and solutes, eqs 10 and 12 should hold even for force fields which provide a poor description of the solution and/or surface thermodynamics. The above two solutes were also chosen as we have developed corresponding force fields which were specifically designed to reproduce the experimental densities and KB integrals (and thereby a22 and η) for the solution mixtures.46,54 The degree of agreement with experiment one can obtain for a22 and η is provided in Figures 1 and 2. The details of the calculations can be found elsewhere.46,54 The agreement with experiment for the value of a22 in NaCl solutions is very good for all but the highest >4 M) concentrations, particularly in comparison with other force fields.46 The agreement for the methanol η values is less quantitative, although the appropriate trend is clearly reproduced. Furthermore, the NaCl solution densities are in very good agreement with experiment.46 The methanol solution densities are reproduced well at low methanol concentrations but are somewhat lower than experiment for high methanol compositions.54 The simulated and experimental surface tension values for the two solutes are displayed in Table 1 and in Figures 1 and 2. As expected, NaCl solutions displayed an increase in surface tension from that of pure water, while methanol solutions exhibited a decrease. Unfortunately, the simulated surface tension of pure SPC/E water is low (59.7 mN/m, no long-range dispersion correction included) compared to the experimental value of 71.6 mN/m.65 We note that the exact value for SPC/E water has been the subject of several recent studies.23,63,66–68 Our value is in agreement with the most recent findings.69–71 On the other hand, the pure methanol value is slightly too high (26.2 mN/m) compared to the value of 21.6 mN/m observed

Solute Effects on the Surface Tension of Liquids

J. Phys. Chem. B, Vol. 112, No. 30, 2008 8979 TABLE 1: Simulated Surface Tensions of Aqueous Solutions of NaCl and Methanola soln

n1 °

water NaCl

4340 4158 4096 3900 CH3OH 3808 3500 2484 1836 1320 900 546 250 0

n2° cs° (M) cs (M) x2° 152 308 616 192 500 828 1102 1320 1502 1634 1744 2000

0.00 0.99 2.00 4.01

0.00 0.000 1.09 2.21 4.34 0.048 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

x2

γ (mN/m) Tsim (ns)

0.000

59.7 60.7 64.4 67.3 55.4 42.0 38.4 34.7 27.7 29.1 26.5 25.1 26.2

0.036 0.106 0.223 0.350 0.479 0.610 0.731 0.873 1.000

30 20 30 20 20 20 20 20 20 20 20 20 20

a All simulations were performed at 300 K for a total simulation time of Tsim. Symbols are ni°, the total number of solute and water molecules used in the simulations; cs° (x2°), the initial bulk solute concentrations; cs (x2), the final bulk solute concentrations. Typical estimated errors for the surface tension values are (0.3 mN/m.

Figure 1. Thermodynamics of aqueous NaCl solutions. Top: The change in surface tension (∆γ in mN/m) with solute molarity (cs). The symbols represent the raw simulation data and the thin line represents the corresponding fit provided by γ ) 59.3 + 1.93cs. The thick line represents the corresponding experimental data (extrapolated beyond 1M). Bottom: Bulk solution activity derivative (a22) as a function of solute molarity. The symbols represent the raw simulation data and the thin line represents the corresponding fit. The thick line represents the corresponding experimental data.

TABLE 2: Simulated and Experimental Surface Adsorption Data for Aqueous Solutions of NaCl and Methanola η

a22 soln water NaCl

cs

x2

∂γ/∂y

Γ2,1

MD exptl MD exptl MD exptl MD

0.00 0.000 1.00 1.00 1.09 1.03 1.01 2.21 1.19 1.19 4.34 1.55 1.72 0.036 CH3OH 0.106 0.223 0.350

55.1 55.1

54.8 54.8 55.4 53.2

53.2 51.7 53.3 56.4

1.93 1.93 1.93 -129 -81 -50 -30

2.08 2.08 2.08 -215 -117 -54 -27

-0.49 -0.86 -1.30 1.15 2.30 3.40 3.47

exptl -0.54 -0.93 -1.26 1.86 3.12 3.43 3.21

a All data correspond to 300 K. See text for symbol definitions. In the case of ∂γ/∂y, y ) cs for NaCl and y ) x2 for methanol. The values of the relative surface adsorptions correspond to those obtained from the fitted surface tension changes and bulk solution properties (a22 and η). Units are cs (M), η (M), ∂γ/∂y (mN/m/M for NaCl and mN/m for methanol), and Γ2,1 in ions or molecules/nm2.

Figure 2. Thermodynamics of aqueous methanol solutions. Top: The change in surface tension (∆γ in mN/m) with solute mole fraction (x2). The circles represent the raw simulation data and the thin line represents the corresponding fit provided by γ ) 60.3 - 33.80x0.5 - 45.43x + 81.85x2 - 37.20x3. Other symbols represent two corresponding experimental data sets. Bottom: Bulk solution η values (in M) as a function of solute mole fraction. The symbols represent the raw simulation data and the thin line represents the corresponding fit. The thick line represents the corresponding experimental data.

experimentally.65 However, the correct trends in surface tension are observed with increasing solute concentration. The change in surface tension with NaCl concentration is almost quantitatively reproduced by the current simulations. The changes for methanol are reasonable. As a direct consequence of the ions being excluded from the surface, and our use of a relatively small system size, the bulk concentrations of NaCl (cs) have been increased slightly from the initial concentrations (cs°). The

opposite is true for the methanol simulations. The correct bulk solution compositions were therefore determined from the slab simulations by averaging the solute and solvent concentrations over the central 4 nm of the approximately 8 nm slab. These results are provided in Table 1 and used in all subsequent figures. In order to proceed any further one requires fitting equations for both the experimental and simulated surface tensions of both solutes. The functional form and quality of the fit can have a significant effect on the calculated surface adsorption values. In the case of NaCl we have assumed linear behavior of the surface tension with solute molarity for the simulated data, as observed experimentally.43,72 However, we note that the slope of the plots for NaCl vary between experimental studies. For instance, values for ∂γ/∂cs of 1.67,65 1.90,72 and 2.08 mN/m/ M43 have appeared in the literature. We have used the most recent value observed for NaCl concentrations up to 1 M.43 Obtaining an appropriate fitting equation for methanol proved more difficult. In addition, there are two similar sets of experimental data available for fitting.65,73 The best functional form appeared to be ∆γ ) ax0.5 + bx + cx2 + dx3. Even then, the corresponding derivatives at larger solute mole fractions (>0.3) varied widely between the two experimental data sets and were therefore considered unreliable. Using the above fits, together with our previous determinations of a22 and η, the corresponding experimental values of ∂γ/∂y (y ) cs or x2) and Γ2,1 were determined for both solutes. The analogous simulated

8980 J. Phys. Chem. B, Vol. 112, No. 30, 2008

Chen and Smith

Figure 3. Surface adsorption (Γ2,1 in ions/nm2) of aqueous NaCl solution/vacuum interfaces as a function of solute molarity (cs). The symbols correspond to the results obtained after integrating the surface probability distributions (eq 6). The solid line is the expected surface exclusion as determined from eq 10 using the simulated values of the surface tension derivative and a22. The thick line is the corresponding experimental result after extrapolation (dotted line). The straight line corresponds to the surface exclusion expected using the experimental surface tension derivative with a22 ) 1 (an ideal solution).

Figure 4. Surface adsorption (Γ2,1 in molecules/nm2) of aqueous methanol solution/vacuum interfaces as a function of solute mole fraction (x2). The symbols correspond to the results obtained after integrating the surface probability distributions (eq 6). The solid line is the expected surface exclusion as determined from eq 12 using the simulated values of the surface tension derivative and η. The thick line is the corresponding experimental result after averaging over two data sets. The dashed line corresponds to the experimental surface adsorption expected for SI solutions (eq 15), while the straight dotted line corresponds to the experimental surface adsorption expected for SI2 solutions (eq 16).

values corresponding to the bulk solvent compositions studied here were obtained by using the same fitting equations as used for the experimental data, but applied to the simulated surface tension data. In addition, the surface adsorption was determined directly from the trajectories using eqs 5 and 6. The data are provided in Table 2. The surface depletion data for NaCl solutions is displayed in Figure 3. It can be seen that the surface depletion predicted from the surface tension changes and the values of a22 (eq 10), are in complete agreement with the surface depletion obtained directly from eqs 5 and 6. Hence, we have numerical proof of the validity of eq 10 and the associated analysis of the computer simulation data. Furthermore, the simulated and experimental data are in excellent agreement providing a high degree of

confidence in the underlying surface distributions. Assuming an ideal solution (a22 ) 1) is clearly acceptable at low (1 M) salt concentrations, but produces significant errors at higher concentrations. Both the simulated and experimental adsorption values display an initial linear dependence on solute concentration, which then levels off at higher solute concentrations. Interestingly, the curve resembles an inverted Langmuir binding isotherm. The same approach has been applied to the methanol solutions and the results are displayed in Figure 4. Methanol displays the expected positive surface adsorption. Again, the validity of eq 12 is verified by the simulated data obtained at low methanol mole fractions. The simulated data for x2 ) 0.223 displays some disagreement. In our opinion, this is a minor problem which is

Solute Effects on the Surface Tension of Liquids

J. Phys. Chem. B, Vol. 112, No. 30, 2008 8981

Figure 5. Surface properties of an aqueous 2.2 M NaCl solution/ vacuum interface. The origin has been shifted so that g1(z ) 0) ) 0.5 for convenience. Top: Surface probability distributions (gi) for water, sodium, and chloride ions as a function of distance from the interface (z). Middle: Surface adsorption (Γ2,1 in ions/nm2) as a function of integration distance (Z) for all ions (black line and circle), sodium, and chloride ions. The thin dashed line is the expected surface exclusion as determined from eq 10 using the simulated values of the surface tension derivative and a22. The thick dashed line corresponds to the experimental surface exclusion provided by eq 10. The adsorptions observed at Z ) 1 nm were taken as the final simulated values. Bottom: The average salt molality (ms) as a function of distance from the interface (z) obtained from the simulations. The dashed line is the average bulk molality (2.32 m).

Figure 6. Surface properties of an aqueous x2 ) 0.106 methanol solution/vacuum interface. The origin has been shifted so that g1(0) ) 0.5 for convenience. Top: Surface probability distributions (gi) for water and methanol as a function of distance from the interface (z). Middle: Surface adsorption (Γ2,1 in molecules/nm2) of methanol (black line and circle) as a function of integration distance (Z). The thin dashed line is the expected surface adsorption as determined from eq 12 using the simulated values of the surface tension derivative and η. The thick dashed line corresponds to the experimental surface adsorption provided by eq 12. The adsorptions observed at Z ) 1 nm were taken as the final simulated values. Bottom: The average methanol mole fraction (x2) as a function of distance from the interface (z) obtained from the simulations. The dashed line is the average bulk mole fraction (0.106).

probably a consequence of either: (i) errors associated with the fitting procedure or (ii) an inability to accurately model the rather large surface excesses using our relatively small system sizes. Hence, the agreement between experiment and simulation is less satisfactory for this system. The data presented in Table 2 indicate that the disagreement with experiment is essentially due to an incorrect description of the surface tension changes, and not due to differences between the η values. Assuming an ideal solution (η ) F) is clearly acceptable at low (10 ns) are required in order to precisely determine the surface tension changes, and the surface adsorption via eqs 5 and 6. In addition, somewhat larger system sizes are also required to ensure that one has an accurate description of the bulk solution. Hence, our use of nonpolarizable models in the current study. One observes in Figures 5 and 6 that the integrated adsorption values do not strictly remain constant beyond the surface region. This behavior has been observed before during our bulk solution simulations and typically disappears as the system size increases.34 As noted earlier, simulations which provide quantitative agreement with experiment for the surface tensions (derivatives) and the bulk solution values of a22 or η, must also then be in quantitatively agreement with the experimental surface adsorption data (Γ2,1). Unfortunately, the values of Γ2,1 represent an integral over the surface probability distributions gi(z) and therefore it is possible that many different surface distributions can produce the same value of the surface adsorption (see eq 6). This is also true for the bulk solution KB integrals Gij. However, in our previous experience with bulk solution simulations this is generally not the case.45,46,74,75 In fact, it is often difficult to reproduce the density and experimental KB integrals with any parameter set (force field), and we have not observed two different force fields providing the same set of three KB integrals. Hence, we expect the same to be true for the surface probability distributions. The results displayed here for NaCl provide essentially quantitative agreement with experiment for ∆γ, a22, and Γ2,1. This has been achieved using a force field without explicit polarization effects, although the normal combination rules for

Na-O interactions were broken. The current force field was capable of reproducing both bulk solution and interfacial thermodynamic data for a system where changes in polarization would be expected to be significant. While this was not the objective of the present simulations, it does suggest that nonpolarizable force fields can be used with some confidence to study surface distributions as long as a well parametrized force field is used. Other common force fields for NaCl solutions can produce values of a22 which are 3-4 times too low, primarily because the G22 values are too large and positive.46 Hence, as eq 6 is exact, they cannot simultaneously reproduce both the change in surface tension and the correct degree of surface adsorption or exclusion. For instance, a previous NaCl force field comparison provides values of 0.34, 0.31, and 0.82 for a22 at 2 M NaCl.46 The current force field provides a value of 1.08 compared to the experimental result of 1.15. Hence, even if all these force fields reproduced the change in surface tension reasonably accurately, the resulting surface exclusion of ions obtained from eq 10 would be a factor of 3.38, 3.71, 1.40, and 1.06 too large, respectively. Clearly, the use of some of these force fields would result in non negligible errors. This type of behavior (the degree of exclusion being dependent on a22) has also been observed in our previous studies of cavity formation in mixed solvents.76 The NaCl results appear to be consistent with the experimental studies of Raymond and Richmond.4 In particular, the above study concluded that there were ions at the interface, but not at the surface. The simulated surface probability distributions presented here displayed a simple exclusion of both salt ions with very few features. This contrasts somewhat with previous simulation results which suggest a more structured interface.8,10 The reasons for this are undoubtedly due to the different force fields employed in each study. We have confidence in our distributions as we have a force field which reproduces both the changes in surface tensions and the activity derivative (a22) for bulk solutions. Finally, we note that for 1:1 salts one must have Γ+,1 ) Γ-,1 in order to preserve an electrically neutral interface. However, this does not mean that there is no surface potential.72 The latter will depend on the asymmetry between the anion and cation distributions.

Solute Effects on the Surface Tension of Liquids

J. Phys. Chem. B, Vol. 112, No. 30, 2008 8983

Figure 8. Variation in the surface structure (in units of ions or molecules/nm2/M) for aqueous NaCl solutions (top) and methanol solutions (bottom) as a function of composition obtained from the experimental data.

In the theory section we investigated some simple cases where the structure of the surface at low solute concentrations, as characterized by I2,1, could be related to the fitting parameters describing the changes in surface tension. Changes in the surface structure at larger concentrations determined from the experimental data are displayed in Figure 8. If one ignores the initial cusp region, due to the Debye-Hueckel contribution to the activity coefficient, then the structure of the interface for NaCl solutions varies essentially linearly with solute concentration. The trend is toward a less structured interface as the salt concentration is increased. Hence, while the surface exclusion increases continually with concentration (see Figure 3), this is due to an increase in the salt concentration, not an increase in the surface structure. We note that using eq 14 one can write

I2,1 ) -

βb 1 n( a22

(25)

and that an excellent fit to the experimental data between 0.5 and 5.3 M is provided by 1/a22 ) 1.118 - 0.124cs, i.e., the value of G22 - G12 is approximately constant (see eq 8). Hence, the observed linear dependence of the above integral on salt concentration. For low concentrations of surface adsorbed solutes eqs 17 and 19 provide

I2,1 ≈

βγ1 B [1 + x2(F1∆G12 - 1 ⁄ C)] 2.303F C

(26)

As B and C are positive constants, and methanol displays a negative deviation from Raoult’s Law, one finds that ∆G12 is negative for methanol solutions (see Figure 2) and therefore the surface structure decreases with solute concentration. This is the behavior observed in Figure 8 where methanol solutions exhibit a sharp decrease in the surface structure even though the surface adsorption is increasing (see Figure 4). In summary, the structure of the interface decreases with increasing solute concentration for both NaCl and methanol solutions and therefore the most structured interface is observed at low solute concentrations. It should be noted that the approach used here has been developed for non volatile solutes where the vapor phase

contains a negligible number of molecules. In this case, one does not require any definition of where the surface is, or how diffuse it may be. Consequently, it is not necessary to locate the Gibbs dividing surface or its variants. The integration can originate at any z coordinate in the vacuum phase and continues until both surface distribution functions reach their bulk solution values. In practice, this requires the simulation of reasonably large systems to ensure the bulk phase distribution can be simulated accurately and the values of Γ2,1(Z) converge to a constant value. Using eq 22 one can reduce the integration step to a simple counting procedure. Equation 4 can be applied to systems with any number of components and KB theory can be used to provide expressions for ∂µ2/∂F2, ∂µ3/∂F2, etc. in terms of KB integrals.32,38 This illustrates several further advantages of the current approach. First, while the required chemical potential derivatives for ternary (and higher) systems are often difficult to obtain experimentally, they can be obtained relatively easy from simulated systems as all rdfs and corresponding KB integrals are available. Second, the chemical potential derivatives can be viewed as (composition dependent) scaling factors relating the different relative adsorptions (Γ2,1, Γ3,1, etc.) to the change in surface tension. If one does not know the values of these derivatives then it is impossible to say which solute contributes the most to the overall surface tension change even when all the Γi,1 values are known. The approach presented here is similar to that of previous studies of Lennard-Jones-type mixtures where the surface adsorptionisobtainedfromintegratedprobabilitydistributions,26,27,77–79 or other more approximate approaches.18 However, the current approach differs in several aspects. First, it is easily extended to multiple solutes. Second, the bulk solution activities are obtained from bulk solution simulations using KB theory. Previous studies have typically used simulated phase equilibria data,27 which is difficult to obtain for nonvolatile solutes and/ or multicomponent systems. The KB analysis described here can be performed for both volatile and non volatile liquid mixtures. Third, we have provided the appropriate expressions for a variety of concentration scales. The methodolgy described above can be extended to include more volatile solutes using the same approach as outlined by Telo da Gama and Evans for determining the surface adsorption.26 However, this was not necessary for the current applications. Conclusions We have provided an approach for the description of solute effects on the surface tension of planar solution/vapor interfaces that is suitable for the analysis of computer simulations involving any number of solute components. The approach uses surface probability distributions and the KB theory of solutions to provide a complete and consistent thermodynamic and structural description of the interface. It is illustrated that the relationships derived here are observed in simulations of binary solutions where the aqueous solute is either adsorbed or excluded from the surface. Furthermore, it is strongly suggested that simulations of the corresponding bulk solutions be performed in order to determine the required chemical potential (activity) derivatives, as these can significantly affect the observed degree of surface adsorption or exclusion, even when the surface tension changes are accurately reproduced. Unfortunately, one cannot assume that current force fields accurately reproduce these derivatives. Finally, we are currently attempting to modify the above approach to treat curved interfaces, and applications involving systems with multiple solutes are also under investigation.

8984 J. Phys. Chem. B, Vol. 112, No. 30, 2008 Acknowledgment. The project described was supported by Grant No. R01GM079277 from the National Institutes of General Medical Sciences. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of General Medical Sciences or the National Institutes of Health. References and Notes (1) Onsager, L.; Samaras, N. N. T. J. Chem. Phys. 1934, 2, 528–536. (2) Pappenheimer, J. R.; Lepie, M. P.; Wyman, J. J. Am. Chem. Soc. 1936, 58, 1851–1855. (3) Richmond, G. L.; Robinson, J. M.; Shannon, V. L. Prog. Surf. Sci. 1988, 28, 1–70. (4) Raymond, E. A.; Richmond, G. L. J. Phys. Chem. B 2004, 108, 5051–5059. (5) Cacace, M. G.; Landau, E. M.; Ramsden, J. J. Quart. ReV. Biophys. 1997, 30, 241–277. (6) Moore, W. J. Physical Chemistry; Prentice-Hall: Englewood Cliffs, NJ, 1972. (7) Wilson, M. A.; Pohorille, A. J. Chem. Phys. 1991, 95, 6005–6013. (8) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361– 6373. (9) Walker, D. S.; Hore, D. K.; Richmond, G. L. J. Phys. Chem. B 2006, 110, 20451–20459. (10) Bhatt, D.; Newman, J.; Radke, C. J. J. Phys. Chem. B 2004, 108, 9077–9084. (11) Chen, B.; Siepmann, J. I.; Klein, M. L. J. Am. Chem. Soc. 2002, 124, 12232–12237. (12) Dang, L. X.; Chang, T. M. J. Chem. Phys. 1997, 106, 8149–8159. (13) Wick, C. D.; Dang, L. X.; Jungwirth, P. J. Chem. Phys. 2006, 125. (14) Kuo, I. F. W.; Mundy, C. J.; Eggimann, B. L.; McGrath, M. J.; Siepmann, J. I.; Chen, B.; Vieceli, J.; Tobias, D. J. J. Phys. Chem. B 2006, 110, 3738–3746. (15) Brown, E. C.; Mucha, M.; Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2005, 109, 7934–7940. (16) Bahadur, R.; Russell, L. M.; Alavi, S. J. Phys. Chem. B 2007, 111, 11989–11996. (17) Paul, S.; Chandra, A. Chem. Phys. Lett. 2003, 373, 87–93. (18) Paul, S.; Chandra, A. J. Chem. Theory Comput. 2005, 1, 1221– 1231. (19) Vrbka, L.; Mucha, M.; Minofar, B.; Jungwirth, P.; Brown, E. C.; Tobias, D. J. Curr. Opin. Colloid Interface Sci. 2004, 9, 67–73. (20) Jungwirth, P.; Tobias, D. J. Chem. ReV. 2006, 106, 1259–1281. (21) Mucha, M.; Frigato, T.; Levering, L. M.; Allen, H. C.; Tobias, D. J.; Dang, L. X.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 7617–7623. (22) Chang, T. M.; Dang, L. X. Chem. ReV. 2006, 106, 1305–1322. (23) Ismail, A. E.; Grest, G. S.; Stevens, M. J. J. Chem. Phys. 2006, 125, 014702. (24) Tronel-Peyroz, E.; Douillard, J. M.; Bennes, R.; Privat, M. Langmuir 1989, 5, 54–58. (25) Randles, J. E. B. Phys. Chem. Liquids 1977, 7, 107–179. (26) Telo Da Gama, M. M.; Evans, R. Mol. Phys. 1983, 48, 229–250. (27) Mecke, M.; Winkelmann, J.; Fischer, J. J. Chem. Phys. 1999, 110, 1188–1194. (28) Gopalakrishnan, S.; Jungwirth, P.; Tobias, D. J.; Allen, H. C. J. Phys. Chem. B 2005, 109, 8861–8872. (29) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774–777. (30) Ben-Naim, A. J. Chem. Phys. 1977, 67, 4884–4890. (31) Ben-Naim, A. Statistical Thermodynamics for Chemists and Biochemists; Plenum Press: New York, 1992. (32) Ben-Naim, A. Molecular Theory of Solutions; Oxford University Press: New York, 2006. (33) Matteoli, E.; Mansoori, G. A. Fluctuation Theory of Mixtures. Vol 2.; Taylor and Francis: New York, 1990. (34) Chitra, R.; Smith, P. E. J. Chem. Phys. 2001, 114, 426–435. (35) Weerasinghe, S.; Pettitt, M. B. Mol. Phys. 1994, 82, 897–912. (36) Chitra, R.; Smith, P. E. J. Phys. Chem. B 2001, 105, 11513–11522. (37) Smith, P. E. J. Phys. Chem. B 2004, 108, 18716–18724.

Chen and Smith (38) Smith, P. E. Biophys. J. 2006, 91, 849–856. (39) Ruckenstein, E.; Shulgin, I. Fluid Phase Equilib. 2001, 180, 281– 297. (40) Friedman, H.; Ramanathan, P. S. J. Phys. Chem. 1970, 74, 3756– 3765. (41) Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1987, 86, 5110–5116. (42) Chitra, R.; Smith, P. E. J. Phys. Chem. B 2002, 106, 1491–1500. (43) Weissenborn, P. K.; Pugh, R. J. Langmuir 1995, 11, 1422–1426. (44) Chitra, R.; Smith, P. E. J. Chem. Phys. 2001, 115, 5521–5530. (45) Weerasinghe, S.; Smith, P. E. J. Phys. Chem. B 2003, 107, 3891– 3898. (46) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 119, 11342– 12349. (47) Perera, A.; Sokolic, F. J. Chem. Phys. 2004, 121, 11272–11282. (48) Auffinger, P.; Cheatham, T. E.; Vaiana, A. C. J. Chem. Theory Comput. 2007, 3, 1851–1859. (49) Parsegian, V. A.; Rand, R. P.; Rau, D. C. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 3987–3992. (50) Baynes, B. M.; Trout, B. L. J. Phys. Chem. B 2003, 107, 14058– 14067. (51) Strey, R.; Viisanen, Y.; Aratono, M.; Kratohvil, J. P.; Yin, Q.; Friberg, S. E. J. Phys. Chem. B 1999, 103, 9112–9116. (52) Meissner, H. P.; Michaels, A. S. Ind. Eng. Chem. 1949, 41, 2782– 2787. (53) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The properties of gases and liquids; McGraw-Hill: New York, 2001. (54) Weerasinghe, S.; Smith, P. E. J. Phys. Chem. B 2005, 109, 15080– 15086. (55) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (56) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (57) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (58) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463–1472. (59) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952– 962. (60) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (61) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (62) Harris, J. G. J. Phys. Chem. 1992, 96, 5077–5086. (63) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574–4583. (64) Kang, M.; Smith, P. E. Fluid Phase Equilib. 2007, 256, 14–19. (65) Weast, R. C. Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1985. (66) Shi, B.; Sinha, S.; Dhir, V. K. J. Chem. Phys. 2006, 124, 204715. (67) Lu, Y. J.; Wei, B. Appl. Phys. Lett. 2006, 89, 164106. (68) Wemhoff, A. P.; Carey, V. P. Int. J. Thermophys. 2006, 27, 413– 436. (69) Wynveen, A.; Bresme, F. J. Chem. Phys. 2006, 124, 104502. (70) Chen, F.; Smith, P. E. J. Chem. Phys. 2007, 126, 221101. (71) in’t Veld, P. J.; Ismail, A. E.; Grest, G. S. J. Chem. Phys. 2007, 127, 144711. (72) Jarvis, N. L.; Scheiman, M. A. J. Phys. Chem. 1968, 72, 74–78. (73) Efremov, Y. V. Russ. J. Phys. Chem. 1968, 42, 1003–1005. (74) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 118, 10663– 10670. (75) Kang, M.; Smith, P. E. J. Comput. Chem. 2006, 27, 1477–1485. (76) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 118, 5901– 5910. (77) Lee, D. J.; Dagama, M. M. T.; Gubbins, K. E. J. Phys. Chem. 1985, 89, 1514–1519. (78) Salomons, E.; Mareschal, M. J. Phys.-Condens. Matter 1991, 3, 3645–3661. (79) Holcomb, C. D.; Clancy, P.; Zollweg, J. A. Mol. Phys. 1993, 78, 437–459.

JP711062A