Molecular Dynamics Study of GasLiquid Aqueous Sodium Halide

J. Phys. Chem. C 2007, 111, 721-737

721

Molecular Dynamics Study of Gas-Liquid Aqueous Sodium Halide Interfaces. I. Flexible and Polarizable Molecular Modeling and Interfacial Properties Tatsuya Ishiyama and Akihiro Morita* Department of Computational Molecular Science, Institute for Molecular Science, Myodaiji, Okazaki 444-8585, Japan ReceiVed: August 11, 2006; In Final Form: October 14, 2006

Gas-liquid interfacial structures of NaCl and NaI aqueous solutions are investigated via molecular dynamics simulations using a flexible and polarizable water model we have developed. The new five-site model of water aims at suitably describing interfacial properties, including vibrational sum frequency spectroscopy where both flexibility and polarization are crucial. The performance of the water model is systematically examined and demonstrated by a number of properties of bulk and interface, including density, vaporization energy, dipole moment, diffusion coefficient, radial distribution function, infrared and Raman spectra of the O-H stretching region, surface potential, and surface excess of ions. The orientational structure of surface water is investigated in detail in connection with the issue of surface solvation of anions. These investigations will be utilized to analyze the sum frequency generation spectra in relation to the orientational structure at the molecular level.

1. Introduction Ion species play crucial roles for the structure of aqueous interfaces in many fields of chemistry, such as interfaces at electrodes in electrochemistry, structure of proteins and membranes, and corrosion processes in environmental chemistry. The effects of ions on the gas-liquid interface structure of aqueous salt solutions have drawn great and renewed interest for the past decade, mainly in connection with the atmospheric chemistry. A number of gas-liquid heterogeneous kinetic measurements imply that uptake reactions of gaseous species with ions in aqueous solutions are accelerated at the interface,1,2 though the mechanism is not very well understood. Surface concentration of chloride anions in NaCl aqueous solution is suggested as a key parameter in tropospheric heterogeneous chemical processes.3,4 It had long been believed that surfaces of aqueous salt solutions are devoid of ions. It is well-known that the surface tensions of aqueous solutions of most inorganic salts are greater than that of pure water, which indicates, through the Gibbs adsorption equation from the thermodynamic point of view, a negative surface excess of ions.5 This fact is interpreted in the continuum dielectric theory that a point charge is repelled from the interface of two different dielectric media due to the effect of the so-called image force.6,7 However, recent molecular dynamics (MD) simulation studies pose challenging results against the above traditional picture of aqueous interfaces. Jungwirth and Tobias showed, by MD simulations of sodium halide aqueous solutions with the point dipole force field model, that the surface concentration of large and highly polarizable halide anions (e.g., Br- and I-) at the air-solution interface is relatively high compared to that in the bulk region.8-10 The enhanced concentration of anions at the interfaces was also demonstrated by other polarizable force fields11 and free energy * To whom correspondence should be addressed. E-mail: amorita@ ims.ac.jp. Also at Graduate University for Advanced Studies, Myodaiji, Okazaki 444-8585, Japan.

calculations.12-14 The interface structure turned out to be sensitive to the force field, as a former MD simulation using a nonpolarizable molecular model did not predict the surface enhancement of anions, though anions are more weakly repelled from the surface than cations.15,16 The predictions by MD calculations have invoked experimental studies to detect the gas-liquid interface structure at a molecular level.17 It is however a challenging issue partly because surface-sensitive techniques applicable to liquid interfaces are quite limited compared to those for the solid interfaces. Most of the experimental studies on the aqueous salt surfaces aim at detecting ions by surface-sensitive techniques. Petersen and Saykally employed second harmonic generation to detect anion species at the aqueous interfaces.18-20 The high pressure X-ray photoelectron spectroscopy experiment by Ghosal et al. also observed the surface enhancement of iodide anion,21 whereas Weber et al. carried out the valence-band photoemission experiment of NaI aqueous solution and interpreted their results as a negative surface excess of ions.22 In contrast to the above attempts to detect the surface ions, sum frequency generation (SFG) spectroscopy can probe the interfacial structure of water by a monolayer sensitivity. SFG is a versatile and powerful experimental technique to study liquid interfaces with a number of applications, as discussed in ref 23 and references therein. Since orientational structure of water is sensitive to the electric field, if it is formed, SFG provides valuable information on the double layer of ions, presumably associated with the segregation of anions. The capabilities of SFG as a surface probe will be augmented by close collaboration with theoretical calculations. Therefore, we have developed theoretical methods to compute the SFG spectroscopy by molecular dynamics simulations, and first applied them to the gas-liquid water surface with considerable success.24,25 In the present and following papers, we extend the treatments and apply them to the analysis of the electrolyte aqueous solutions. For this purpose, first, we discuss the development of a flexible and polarizable (FP) water model

10.1021/jp065191s CCC: $37.00 © 2007 American Chemical Society Published on Web 12/10/2006

722 J. Phys. Chem. C, Vol. 111, No. 2, 2007 suitable to the SFG computation, a key ingredient for the SFG analysis. The molecular model has to allow internal vibration for use of the time-dependent formulation of SFG and also polarization to discuss the optical response. In our previous simulations of a pure water surface, however, the treatment of polarization was rather incomplete, where nonpolarizable models of water were used to generate MD trajectories for convenience of MD calculations. Consequently, the effect of polarization was not incorporated in the structure and dynamics of the interface, while it was accounted for by the calculations of instantaneous dipole and polarizability. This approximate treatment of polarization is no more valid for the electrolyte solutions, since the polarizable interaction is essential to examine the surface segregation of anions. In the past decade, development of FP water models has been an attractive avenue for further improvements of empirical force fields.26-33 The performances of these FP models have been investigated in the bulk liquid simulations with particular focus on their structure, energetics, thermodynamics, transport dynamics, dielectrics, spectroscopy, and second virial coefficient. Although the FP models are expected to be more elaborate than other models including only internal vibration or polarization, the implementation of these FP models encountered several practical problems that need to be overcome. First, early FP models failed to reproduce the tetrahedral structure for liquid water, characterized by the second peak of the oxygen-oxygen radial distribution function.26,27 This problem is resolved by introducing an additional interaction site in the present study, as suggested by ref 34 (see sections 2 and 3). The second problem is related to general polarizable models with the point charge/polarizability approximation, which breaks down at a short internuclear separation and causes the so-called “polarization catastrophe”. This problem becomes more troublesome when the flexibility is incorporated in the models. The proposed remedy for this problem is to apply a screening function for short internuclear Coulombic interaction or smeared charge distribution, with various functional forms, such as simple shield type,26 Thole type,33,35 Gaussian type,27,32 and Slater type.29 In the present study, we use the Thole type35 or its alternative type36 to prevent the second problem (see section 2). The third problem is about the equilibrium geometry of water molecules. Most of the classical flexible water models yield smaller H-O-H angles of a water molecule in the condensed phase than in the gas phase, whereas a recent experiment37 and ab initio calculations38,39 have shown the opposite trend. While this is a common problem to most of the general flexible water models, the polarization effect magnifies this problem. Burnham and Xantheas resolved this problem by elaborating the geometry dependence of polarization, which properly results in an increased H-O-H angle in water clusters.30 Ren and Ponder, on the other hand, empirically adjusted the equilibrium H-O-H angle of an isolated molecule so as to reproduce the proper geometry in the condensed phase.33 The former remedy for the third problem is physically more reasonable, though the latter works as well to practically describe water properties in the condensed phase. We employ the latter idea in the present study because of its simplicity (see sections 2 and 3). The FP model of water in the present work is motivated by the interfacial properties, including vibrational SFG spectroscopy in particular, of aqueous systems. To our knowledge, the present study is the first attempt to apply a FP force field to the interfaces of pure water and aqueous solutions. The main purpose of the present study (part I) is to construct a FP force field of water that well reproduces not only bulk but also

Ishiyama and Morita interfacial properties of pure water and aqueous electrolyte solutions. The model is utilized to study the SFG spectroscopy for these interfaces in the following paper (part II).23 The next section describes the details of the FP water force field model and molecular dynamics simulation based on the model. To check the performance of the present model, we examine several bulk properties by the molecular dynamics simulation and compare those with available experiments in section 3. In section 4, the interfacial properties of pure water and aqueous sodium halide solutions are discussed with the present model, including density profiles, water orientational order, induced dipole moment profiles, and surface potential difference between pure water and the aqueous solutions, and surface excess of ions. Then, the concluding remarks follow in section 5. 2. Molecular Models and Simulation 2.1. Force Fields. Using the FP molecular model, the total potential energy of the MD system is written as the sum of the intramolecular vibrational energy, Uintra, the intermolecular pairwise additive energy consisting of the van der Waals and Coulomb interactions, Upair, and the nonadditive energy associated with the induced polarization, Upol:

Utotal ) Uintra + Upair + Upol

(1)

In the present simulations, the electronic polarization is represented by an induced point dipole located at the molecular center of mass. Geometry. The present FP force field of water consists of five interaction sites, that is, three atomic sites of hydrogen (H) and oxygen (O), the center of mass (G), and an in-plane auxiliary site (M). The equilibrium geometry of the three atomic sites is eq the same as that of the SPC model,40 req OH ) 1.0 Å and θHOH ) 109.5°. The center of mass (G) is determined by the three atomic sites with their standard atomic weights,41 and the auxiliary site (M) is located along the H-O-H bisector toward the H atoms with a distance of rOM ) 0.21 Å. The locations of the sites G and M are constrained by the above conditions during the intramolecular vibration. A Lennard-Jones core is placed on the oxygen (O) site, while three point charges are on the two H sites and the site M. A single point dipole is located on the center of mass (G) to represent the electronic polarization. Intramolecular Terms. The potential function Uintra for the internal vibration of water consists of the harmonic and anharmonic O-H stretching vibration terms, harmonic bending term, and these coupling terms: 6

Uintra )

∑i n)2 ∑kn(ri,1n + ri,2n) + (kθ/2)ri,32 + [

krθri,3(ri,1 + ri,2) + krr′γi,1ri,2] (2) where ri,1 and ri,2 are the displacements of the two O-H bond lengths from the equilibrium O-H separation and ri,3 is that from the equilibrium H-H separation of the ith molecule. The parameters in eq 2 are tabulated in Table 1; they are empirically determined to reproduce the experimental IR, Raman, and SFG spectra of condensed water, as discussed later in section 3.3 and the following paper.23 We found that the higher order anharmonicity terms in the O-H stretching in eq 2 are needed to reproduce the experimental red shifts in the O-H stretching vibrations, though these terms are fully empirical in nature. The bending and the coupling terms are followed by Marti et al.,42

Flexible and Polarizable Molecular Modeling

J. Phys. Chem. C, Vol. 111, No. 2, 2007 723

TABLE 1: Intramolecular Potential Parameters of a Water Molecule force constants

unit

k2 ) 4.359 k3 ) -5.295 k4 ) 11.118 k5 ) -409.669 k6 ) 1389.472 kθ ) 1.697 krθ ) -1.280 krr′ ) 1.055

[mdyn Å-1] [mdyn Å-2] [mdyn Å-3] [mdyn Å-4] [mdyn Å-5] [mdyn Å-1] [mdyn Å-1] [mdyn Å-1]

The induced dipole, µi, is defined by the following selfconsistency condition:

µi ) riEi ) ri[E0i -

σ (Å)

(kJ mol-1)

O H M G

3.170 0.000 0.000 0.000

Water 0.765 0.000 0.000 0.000 Ionsa

Na+ ClI-

2.583 4.401 5.167

0.4184 0.4184 0.4184

Tij )

0.000 +0.5245 -1.0490 0.000

eq 8

+1 -1 -1

0.250 3.250 6.900

The σ and values are from refs 65 and 94, and the R values are from refs 95-97.

where the original value of krr′ ) 0.88 mdyn Å-1 is adjusted to the present FP model. Intermolecular Pairwise Terms. For the intermolecular pairwise part, Upair, we employ the usual Lennard-Jones potentials, ULJ, and Coulombic potentials, UC:

Upair ) ULJ + UC )

∑ ∑ i>j a,b

[ {( ) ( ) } σab

4ab

12

-

rai,bj

σab

6

rai,bj

+

]

(1) fai,bj qaqb

rai,bj

(3)

where rai,bj is the separation between site a of the ith molecule and site b of the jth molecule, ab and σab the Lennard-Jones potential parameters between sites a and b, and qa the partial charge at site a. The Lennard-Jones parameters between different sites are evaluated by the Lorentz-Berthelot mixing rules: σab ) (σa + σb)/2 and ab ) xab. The parameters used in the present study are listed in Table 2. The partial charges of water sites are consistent with the dipole moment of an isolated water (1) molecule, 1.85 D.43 fai,bj in eq 3 is a short-range damping factor for Coulomb interactions, as described later. Polarization Terms. The polarization energy associated with the induced dipoles, Upol, can be calculated by44

Upol ) -

1 2

∑i

µi E0i

(4)

where µi is the induced dipole moment on the ith molecule which is self-consistently determined. E0i is the electrostatic field exerted on site G of the ith molecule by the surrounding permanent charges: site

E0i

)

∑∑

j(*i) b

(2) fi,bj qb

ri,bj3

ri,bj

rij3

I-

3wij4 rij5

rijrij

(7)

where I is the identity matrix and wij represents a short-range damping factor for the dipole tensor. The polarizability tensor of a water molecule has an anisotropic form that depends on both the water orientation and the internal conformation, which was determined by quantum chemical calculations:25

a

site

4wij3 - 3wij4

R (Å3)

q (e)

(6)

where ri is the polarizability tensor of the ith molecule and Ei is the total electrostatic field acting on site G of the ith molecule. Tij is the dipole-dipole interaction tensor:

TABLE 2: Intermolecular Potential Parameters of Water and Ions site

∑ Tijµj]

j(*i)

(5)

(2) where fi,bj is the second short-range damping factor for the electric field acting on site G.

[

† † ri,1Ai,1 + Ai,2 ri,2Ai,2 Ri ) Ai,1

D1 + D4ri,k 0 0 D + D r 0 ri,k ) 2 5 i,k 0 D3 + D6ri,k + D7ri,k2 0 0

]

(8a)

(8b)

where the suffix k ()1,2) is associated with each O-H bond, k. Equation 8 represents the water polarizability as the sum of contributions from the two O-H bonds. Ai,k is the direction cosine matrix to rotate the space fixed coordinates to the local coordinates for the O-H bond, k, where z is along the O-H bond and y is normal to the molecular plane. ri,k is the bond polarizability for the O-H bond, k, defined in the local frame by eq 8b, where the Di values are cited from ref 25. On the other hand, the polarizability tensors of ions are isotropic, and their values are tabulated in Table 2. (1) Damping Factors. In eqs 3, 5, and 7, damping factors fai,bj , (2) fi,bj , and wij are involved to avoid the so-called “polarization catastrophe” at short intermolecular distance,35 as discussed in section 1. In particular, flexible force field models tend to facilitate the unphysical events, as the deformation of molecular geometry could help accelerate the divergence of polarization. To avoid this problem, we introduce damping functions optimized for the Coulomb interaction part, electric field, and (1) dipole tensor, respectively. The damping factor fai,bj for the Coulomb terms takes the functional form proposed by Thole:35 (1) (r) ) fai,bj

{

(r/s)4 - 2(r/s)3 + 2(r/s) (r < s) (r > s) 1

(9)

where r ) rai,bj is the distance between site a of the ith molecule and site b of the jth molecule. The scaling distance, s, for sites a and b in eq 9 is optimized to the present FP model as sab ) 2.2(RaRb)1/6, where Ra is a volume parameter for site a. For the sites of water, RO ) 0.862 Å3 and RH ) 0.514 Å3 are employed,35 and for the ions, the polarizability values in Table 2 are employed. On the other hand, for the electric field and dipole tensor, the damping functions developed by Levy and co-workers36 are used: (2) (r) ) 1 - exp(-γ(r/s′)n) fi,bj

(10)

724 J. Phys. Chem. C, Vol. 111, No. 2, 2007

wij(r) )

{

r/s′ + (1/m)[1 - (r/s′)m] (r < s′) (r > s′) 1

Ishiyama and Morita TABLE 3: Bulk Properties of Pure Water

(11)

where γ ) 0.85, n ) 3, and m ) 32 are employed.36 r ) ri,bj in eq 10 is the distance between the center of mass (site G) of the ith molecule and site b of the jth molecule, and r ) rij in eq 11 is the center-of-mass distance between the ith and jth molecules. The scaling distance, s′, in eqs 10 and 11 is 1.662(RaRb)1/6,36 where the volume parameter for the center of mass of water is RH2O ) 1.44 Å3.45 For a periodic system, which is treated by the Ewald summation as described in section 2.2, we note that the shortrange damping functions are applied to the real part of Ewald terms. The details are discussed in the Appendix. 2.2. MD Simulation. Molecular dynamics calculations for pure water, NaCl, and NaI aqueous solutions were performed using the force fields of water and ions described in section 2.1. The MD computations are executed using a slab geometry of liquid in three-dimensional periodic boundary conditions. The periodic cell is rectangular with dimensions of Lx × Ly × Lz ) 30 Å × 30 Å × 150 Å, where the gas-liquid interfaces are normal to the z axis. The number of molecules per each cell is 1000 for the pure water system, 1000 (water) + 20 (anion) + 20 (cation) for 1.1 M aqueous solution systems, and 1000 (water) + 40 (anion) + 40 (cation) for 2.1 M systems. For the correction of the long-range electrostatic forces, the Ewald summation method is employed for partial charges and dipoles.46,47 The Ewald separation parameter is set to 0.242 Å-1, and the real and reciprocal-space cutoffs are 13 Å and 1.47 Å-1, respectively, the values of which are the optimized parameters to maximize computational performance.48 The details of the Ewald summation for the dipolar systems including the shortrange damping functions are presented in the Appendix. The van der Waals interactions and the real part of the electrostatic interactions were calculated with a cutoff length of 13 Å, and the Verlet neighbor list49 was employed for these calculations with a shell thickness of 2 Å. Newton’s equations of motion are integrated numerically using the velocity Verlet algorithm49 with a time step of 0.61 fs. At each time step, eq 6 for the self-consistent induced polarization is solved iteratively until the convergence is obtained, where the average root-mean-square differences of the induced dipoles for successive iterations,

x∑Ni |µi(new)-µi(old)|2/N, is below a threshold of 2.5 × 10-10

D. The molecules are initially placed with random displacements and orientations from cubic lattice points to form the slab geometry, and the initial velocities are randomized according to the Maxwell-Boltzmann distribution at a temperature of 298 K. The systems of 256 replicas are first equilibrated in parallel for 100 ps each, at the constant temperature using the Berendsen thermostat50 with a coupling constant of 0.4 ps. Then, the statistical samplings are taken in parallel for a total of 30 ns under the NVE ensemble, with the average kinetic temperature of 298 K. The calculations and analysis procedures of bulk or interface properties are discussed in sections 3 and 4. The program code for the MD calculations is written by us. The MD calculations were performed on the Hitachi HA8000 parallel machine at the Research Center for Computational Science, Okazaki, Japan. 3. Bulk Properties To assess the reliability of the present force field model, we evaluate representative properties related to the bulk water and aqueous solutions. Some of the bulk properties in this

F (g/cm3) U/N (kJ/mol) D (109 m2/s) p (Debye)

this model

experiment

0.995 -41.4 2.59 2.66d

0.997a -41.5b 2.30c 2.9 ( 0.6e

a Reference 98. b Reference 99. c Reference 100. d p ) µperm + µ ) 2.11 + 0.55 (D). e Reference 101.

section are used as reference data for optimizing the force field parameters in section 2.1, such as density, vaporization energy, site-site radial distribution functions, and vibrational spectra. 3.1. Pure Water. Table 3 summarizes the bulk density, F, average potential energy per molecule (vaporization energy), U/N, diffusion coefficient, D, average total [permanent (µperm) + induced (µ)] dipole moment per molecule, p. The diffusion coefficient, D, is derived from the mean-squared displacement at long times:

D ) lim tf∞

〈

∑|ri(t) - ri(0)|2 N i

1 1 6t

N

〉

(12)

where ri(t) is the center-of-mass position of the ith water molecule at time t. Table 3 shows that the calculated values of F and U/N well reproduce the experimental ones, as these values are used as a part of reference data for the potential parameters. The present model also reasonably reproduces the experimental values of D and p, though the force field parameters are not optimized for D or p. The calculated values in Table 3 are in fairly good agreement with those of other force fields commonly used (see, for example, Table 2 of ref 51 for mutual comparison). Regarding the equilibrium molecular geometry in liquid water, the average bond length and angle are, respectively, 1.026 Å and 104.1°, which differ from those of the isolated water model by +0.026 Å and -5.4°. In particular, the noticeable decrease in the bond angle is noted as a shortcoming of a flexible water model, as discussed in section 1. In the present model, the decrease in the angle is empirically offset by the tetrahedral angle (109.5°) of the isolated model. The result of the average angle in liquid is thereby consistent with the recent neutron diffraction experiment ranging from 103.7 to 105.4°.52 Radial Distribution Functions. Figure 1 shows the calculated O-O, O-H, and H-H radial distribution functions, gOO, gOH, and gHH, in pure water, where the experimental results at 25 °C by Soper and Phillips53 are also displayed for comparison. The present calculations well reproduce the experimental peak positions and heights of the three site-site radial distribution functions (RDFs), in particular the second peak of gOO at rOO ∼ 4.5 Å indicative of the tetrahedral coordination of the hydrogen bond network and the marginal third peak of gOO at rOO ∼ 6.5 Å. The second peak is quite sensitive to the position of the additional site M, where the negative partial charge is located.34 We confirmed that a FP model with rOM ) 0 Å, such as PSPC,44 results in the absence of the peak at the tetrahedral coordination.26 The position of site M is optimized to best describe the three site-site RDFs. The hydration number in the first hydration shell, Nh, is calculated by integrating gOO up to the first minimum position, R1 ()3.34 Å for gOO):

Nh ) nw

∫0R g(r)4πr2 dr 1

(13)



Figure 1. (a) O-O, (b) O-H, and (c) H-H radial distribution functions in pure water. The dashed lines indicate experimental results.53 (Reproduced with permission. Copyright 1986 Elsevier.)

TABLE 4: Bulk Properties of Aqueous Solutions Density, F (g/cm3) calcd exptla

1.1 M NaCl

2.1 M NaCl

1.1 M NaI

2.1 M NaI

1.042 1.040

1.077 1.080

1.113 1.118

1.201 1.227

Diffusion Coefficients, D (109 m2/s) 1.1 M NaCl

2.1 M NaCl

1.1 M NaI

2.1 M NaI

calcd exptla,b

2.11

H2O 1.74 2.01

2.22 2.24

1.96 2.22

calcd exptla,b

0.89 1.23

Cation 0.73 1.13

0.98 1.26

0.86 1.24

1.18 1.77

Anion 0.98 1.61

1.44 1.76

1.26 1.61

calcd exptla,b

a Reference 100. b The experimental diffusion coefficients are reported at 1.0 and 2.0 M.

where nw is the bulk number density of water. The calculated Nh value is 5.08, which is in reasonable agreement with the experimental value 5.2 ( 0.4.54 3.2. Aqueous Solutions. Next, we discuss the bulk properties of NaCl and NaI aqueous solutions of 1.1 and 2.1 M concentrations. Table 4 compares the calculated densities, F, and selfdiffusion coefficients, D, of constituent species to available experimental values. The calculated densities, F, very well reproduce the experiments, including the dependence of solute and concentration. On the other hand, the self-diffusion coefficients, D, of water and ions agree with the experimental values in the order of magnitude, and the concentration dependence of D is properly elucidated. To assess the accuracy of D, the present model is comparable in performance to other nonpolarizable force field models reported to date.55-58 Here, we note that the diffusion behavior of solute ions strongly depends on the treatment of long-range interaction. Properly incorporating

Figure 2. (a) O-O, (b) O-cation, and (c) O-anion radial distribution functions in NaCl and NaI aqueous solutions.

the long-range Coulombic interaction by the Ewald summation is essential for an accurate description of the ionic mobility,59 and other treatments using cutoff procedures may cause unexpected artifacts.60,61 In our preliminary simulations employing cutoff treatments based on the spherical truncation with a tapering function at the boundary62 or the minimum-image convention, we experienced exceedingly small D values of solute ions. Radial Distribution Functions. Next, we discuss RDFs in aqueous solutions. Figure 2 shows the site-site RDFs around the water oxygen site O, that is, O-O, O-cation (Na+), and O-anion (Cl- or I-). The O-O RDF, gOO, in Figure 2a is quite analogous with that of pure water in Figure 1a, implying little perturbation by the presence of ions on gOO within the present concentration range. As shown in Figure 2b and c, the RDFs between water oxygen and ions show little concentration dependence of ions and counterion effect. This tendency at room temperature is in contrast to the aqueous solutions at higher temperatures where ion pairing is facilitated, and hence, the O-Na+ RDF shows strong concentration dependence.63,64 The O-ion RDFs in Figure 2b and c can derive the hydration numbers, Nh, around the ion species, Na+, Cl-, and I-, by using eq 13. Table 5 shows the calculated results of Nh for the ions,


Ishiyama and Morita

TABLE 5: Calculated Hydration Numbers, Nh, of Ions 1.1 M NaCl 2.1 M NaCl 1.1 M NaI 2.1 M NaI ref 65 Na+-water Cl--water I--water a

6.1 7.1

6.2 7.3

6.3

6.2

7.4

7.3

a

5.8 7.2 7.9

MD at infinite dilution with a rigid nonpolarizable model.

ones, indicating that the SSIP states are more stable in the present simulations. There are other MD studies reporting ionion RDFs or PMFs, which are in a qualitative discrepancy each other on the relative stability of the CIP and SSIP states.57,66,67 Lyubartsev and Laaksonen reported that the CIP state is more stable than the SSIP in NaCl aqueous solution with a flexible but nonpolarizable force field.57 Guardia et al. concluded that SSIP is preferred in NaCl solution using a flexible force field model,66 while Smith and Dang examined the NaCl PMF in an infinite dilution system with a polarizable force field, reporting a more stable CIP state.67 The relative stability strongly depends on the force field employed, since it is governed by the competition of stability between the cation-anion direct interaction and the hydration around each ion. Experimental investigations associated with this problem will be needed to construct a more reliable ion-ion interaction potential. The present FP model tends to prefer the SSIP state, presumably because the direct anion-cation Coulomb interaction is smeared at a short internuclear separation, as discussed in section 2, at which the validity of the point charge approximation is questionable.35 On the other hand, the cation-cation and anion-anion RDFs shown in parts b and c of Figure 3, respectively, are qualitatively similar to those by Lyubartsev and Laaksonen,57 implying that these RDFs are less sensitive to the force field. 3.3. IR and Raman Spectra. The line shapes of infrared (IR) and Raman spectra are described with the dipole moment, M, and polarizability, A, of the whole system, respectively.69 Accordingly, we briefly mention the calculations of A and M; the details are described in ref 25. The system dipole moment, M, and polarizability, A, are expressed as

M(t) )

({rai}) + µi] ∑i pi(t) ) ∑i [µperm i

(15)

∑i reffi (t)

(16)

A(t) )

Figure 3. (a) Cation-anion, (b) cation-cation, and (c) anion-anion radial distribution functions in the aqueous solutions.

which are in fair agreement with previous MD results at infinite dilution with a rigid nonpolarizable model.65 Ion-ion interaction in aqueous solutions, particularly cationanion association, has been investigated either by RDFs or equivalently by the potentials of mean force (PMF, W(r)).63,64,66-68 The latter is related to the former through

W(r) ) -kBT ln g(r)

where pi(t) is the dipole moment of the ith molecule at time t, consisting of the permanent term, µperm ({rai}), and the induced i term, µi. The analytical form of the former is presented in ref 25 as a function of the intramolecular configuration {rai}, and the latter is calculated by eq 6. reff i (t) in eq 16 is the effective polarizability defined as ri(t) fi(t), where fi(t) denotes the local field correction factor.25 Note that eqs 15 and 16 will be used also in the SFG calculations in the following paper.23 The line shape function for the IR absorption, IIR(ω), is defined via

IIR(ω) ≡

1 2π

∫-∞∞dt e-iωt〈M(0)‚M(t)〉

(17)

where the time correlation function is given quantum mechanically. The absorption coefficient, Rabs(ω), is69

(14)

where kB is the Boltzmann constant and T is the temperature. Figure 3 shows the cation-anion, cation-cation, and anionanion RDFs calculated in aqueous solutions. In the cationanion RDFs in panel a, two noticeable peaks exist beside a third broad maximum. The first peak at r ∼ 3 Å represents the socalled contact ion pair (CIP) configuration, and the second peak at r ∼ 5 Å, the solvent-separated ion pair (SSIP).67 This panel shows that the SSIP peaks are considerably higher than the CIP

Rabs(ω) )

4π2ω(1 - e-pω/kBT) IR I (ω) 3Vpcn(ω)

(18)

where V is the volume, p is the Planck constant divided by 2π, c is the speed of light, n(ω) is the (real part of) refractive index. The classical analogue of the line shape function, IIR cl (ω), is defined in the same way and is related to eq 17 with the harmonic correction and symmetrized factors:70


IIR cl (ω) )

1 2π


∫-∞∞dt e-iωt〈M(0)‚M(t)〉cl )

[

( )]

pω pω coth 2kBT 2kBT

-1

1 + e-pω/kBT IR I (ω) (19) 2

where 〈 〉cl denotes the classical correlation function. The absorption coefficient, Rabs(ω), is thus expressed using IclIR(ω) as

Rabs(ω) )

4π2ω2 ω2 IR IclIR(ω) ∼ I (ω) 3Vcn(ω)kBT n(ω) cl

(20)

On the other hand, the classical line shape functions of polarized and depolarized Raman scattering are introduced analogously:

IRaman pol,cl (ω) ≡ IRaman depol,cl(ω) ≡

∫0∞dt e-iωt〈Ah (0) Ah (t)〉cl

(21)

∫0∞dt e-iωt〈Tr(B(0)‚B(t))〉cl

(22)

1 2π

1 2π

Figure 4. Calculated IR spectrum, Rabs(ω), of pure water. The inset shows the experimental IR spectrum.71 (Reproduced with permission. Copyright 1996 Society for Applied Spectroscopy.)

where A h ) (1/3) Tr A is the spherical part of the polarizability of the entire system and B ) A - A h I is the anisotropic part. The differential cross sections for polarized and depolarized Raman scattering, (d2σpol/dω dΩ) and (d2σdepol/dω dΩ), are

( ) (

)

( ) (

)

4 d2σpol 1 2πn pω ) I Raman(ω) (ω0 - ω) dω dΩ c kBT 1 - e-pω/kBT pol,cl (23) 2

4

d σdepol 1 2πn pω ) I Raman (ω) (ω0 - ω) dω dΩ c kBT 1 - e-pω/kBT depol,cl (24)

where ω0 is the frequency of the incident radiation. The IR spectrum of pure water, Rabs(ω), is evaluated via eq 20 with calibration by the experimental values of n(ω).71 The calculated IR spectrum is shown in Figure 4 and compared with the experimental spectrum by Bertie and Lan.71 The peak position at about 3400 cm-1 and the spectral shape of the calculated Rabs(ω) are in good agreement with the experimental ones, though the present calculation yields a slightly smaller relative intensity at the high and low frequency tail regions. While the spectral shape in Figure 4 is fairly well reproduced, we also note that the absolute value of the absorption coefficient is underestimated by a few orders of magnitude, due to the accuracy of the “quantum correction factor”.72 Thus, we focus on the analysis of spectral shape in the O-H stretching region. The Raman spectra are calculated via eqs 23 and 24, where the incident light wavelength, 2πc/ω0, is set to be 488 nm and the dispersion of n is neglected. To facilitate comparison with experimental scattering of perpendicular polarization, the depolarized scattering is scaled by 1/10, since (d2σ/dω dΩ)⊥ ) (1/10)(d2σdepol/dω dΩ).69 Figure 5 shows the calculated Raman spectra in comparison with the experimental results by Brooker et al.73 The calculated Raman spectra reasonably well reproduce the relative intensity of experimental polarized to depolarized spectra and the line shapes, except for the apparently weak intensity of the calculated polarized spectrum at about 3200 cm-1. The experimental band at about 3200 cm-1 presumably corresponds to the calculated band component blue-shifted to the maximum at 3350 cm-1. It is difficult to optimize the empirical force field to the polarized spectrum in this frequency region while retaining the current accuracy for the other properties. The intensity of this region is known to be

Figure 5. Calculated polarized (red line) and depolarized (blue line) Raman spectra, (d2σpol/dω dΩ) and (d2σ/dω dΩ)⊥, of pure water. The inset shows the experimental Raman spectra.73 (Reproduced with permission. Copyright 1989 John Wiley & Sons.)

emphasized in ice, and collective excitation of intermolecular O-H coupled vibration is suggested to be important.74 While the present FP model captures the spectral features of IR and Raman, further accurate description of water molecules with a large red shift in O-H vibration should require a more sophisticated modeling of intermolecular interactions and/or perturbation on the molecular polarizability. 4. Interfacial Properties In this section, we apply the FP model to gas-liquid interfacial properties of water and aqueous solutions to investigate interfacial structure in detail. This information will be valuable to interpret vibrational sum frequency spectra in molecular detail with the help of MD simulation.23 The present MD simulation employed the slab geometry of liquid, as discussed in section 2.2, and therefore, the gas-liquid interface lies on both sides of the slab. The two interface structures are sampled independently and averaged in the statistical sense. 4.1. Surface Density Profiles. Figure 6 shows the density profiles of constituent molecules and ions near the gas-liquid interfaces. In this figure and hereafter, we use a local z coordinate, zˆ ) z - zGibbs, for convenience of describing gasliquid interfacial structure, where zˆ ) 0 is located at the Gibbs dividing surface for solvent water and zˆ > 0 (zˆ < 0) denotes the gas (liquid) side. The density profiles of solvent water at the aqueous surfaces in Figure 6 are well represented with the following function:

F(zˆ) )

a zˆ 1 - tanh 2 0.455δ

[

(

)]

(25)


Figure 6. Calculated density profiles at the gas-liquid interfaces of pure water (green dashed lines in panels a and b), NaCl (panel a), and NaI (panel b) aqueous solutions. The ordinate is normalized with the bulk density of each component.

where a is the bulk liquid density and δ is called the 10-90 thickness. These parameters are determined by a least-squares method,75 and the calculated 10-90 thickness parameters, δ10-90, are summarized in Table 6. These values are close to 3 Å for all of the aqueous interfaces, which are in good agreement with the ones found by the SPC/E model76 or by other commonly used water models.77 For the NaCl solutions shown in Figure 6a, the density profiles of Na+ and Cl- have relatively small peaks with normalized heights, F/Fbulk ) 1.2-1.4 at zˆ ∼ -5 Å. In these normalized profiles, the concentration dependence of NaCl is found to be minor. These profiles are in accord with the conventional picture of aqueous surfaces devoid of ions, though we also observe a layer slightly enriched with ions about 5 Å below the Gibbs dividing surface. We also note that the profiles of Na+ and Cl- are analogous each other, while a little difference in the peak locations indicates a weak electric double layer structure between Na+ and I- layers. These results are qualitatively consistent with previous MD studies by Jungwirth and co-workers with rigid polarizable force field models,8,9,78 though they reported the anions penetrate the interface slightly more. Figure 6b shows the density profiles of NaI aqueous solutions, which clearly indicates the strong surface enrichment of I- at about zˆ ) -2 to -3 Å. Na+ also shows a comparable peak in the density profile at a different position at zˆ ) -5 to -6 Å. The peaks in I- and Na+ form apparent electric double layers with a separation of about 3-4 Å. These behaviors support the previously reported results using rigid polarizable force field models.8,9,11 Regarding the concentration dependence, the normalized peak heights, F/Fbulk, for both Na+ and I- decrease with ion concentration, presumably because of the ion-ion

Ishiyama and Morita electrostatic repulsion. Comparing the present simulation with the previous works in a quantitative sense, the present simulation predicts a rather weak double layer structure in a somewhat deeper region. The present peak value of F/Fbulk for I- is about 1.75 in 1.1 M NaI solution, which tends to be lower than previously reported values, 2.5-3.5.8,9,11 The maximum of Idensity is observed in Figure 6b at a distance about 2.5 Å below the Gibbs dividing surface toward the liquid, whereas the maximum of I- is located at about 1.0 Å below the Gibbs dividing surface in other previous studies.8,9,11 These quantitative discrepancies are attributed to the difference in the force fields. For our preliminary simulations employing other force fields, we found that one of the important differences in the force fields between the present and the previous studies is in the shortrange damping treatment for the electrostatic interaction, particularly for the induced dipolar interaction, as discussed in section 2. Although the polarizable force fields commonly predict the electric double layer structure of the segregated Na+ and I- layers, the details of the ion density profiles are sensitive to the force fields. The present FP force field well reproduces the experimental surface properties reflecting the local structure, such as the features of surface potential and the SFG spectra, though more details of the force fields should be further examined. 4.2. Orientational Order of Water Molecules. Orientational structure of water at liquid interfaces is straightforwardly investigated by MD simulation with a desirable depth resolution. The water orientation is defined with two coordinates, θdip and φ, as schematically shown in Figure 7. θdip is an angle between the water permanent dipole vector and the surface normal vector, and φ is a dihedral angle between the molecular plane and the plane containing the surface normal and the permanent dipole. These coordinates of molecular orientation are described using the vector of permanent dipole, which is determined by the partial charges and their variable geometry. The orientation of the induced dipole, on the other hand, will be discussed in section 4.3. We first show the average profile of cos θdip in Figure 8 to illustrate the perturbation mechanism of ions on the water orientation. Then, the full probability distribution of water orientation79,80 (θdip, φ) is discussed as a function of the vertical depth, zˆ, in Figure 9. Figure 8 shows 〈cos θdip(zˆ)〉 of water molecules in pure water, NaCl, and NaI aqueous solutions. (Note that the gas region zˆ J 5 Å is noisy due to scarce sampling.) It is noteworthy in pure water that 〈cos θdip〉 changes its sign near the Gibbs dividing surface. The negative 〈cos θdip〉 in the liquid side means that the water molecules direct their dipoles toward the liquid phase, and vice versa in the vapor side.81,82 In NaCl aqueous solutions denoted with blue lines in Figure 8, the 〈cos θdip(zˆ)〉 profile of water orientation is hardly perturbed except for a slightly positive region at about zˆ ∼ -5 Å. This orientational change is attributed to the weak electric double layers of Na+ and Cl-, which are much weaker but essentially exert the same perturbation as those of the NaI case discussed below. In NaI aqueous solution denoted with red lines in Figure 8, the perturbation is more conspicuous. There exists a large positive 〈cos θdip〉 region in the liquid side of the Gibbs dividing surface. This positive region indicates that the water molecules direct their dipoles toward the vapor phase, apparently because the electric double layers of Na+ and I-, as seen in Figure 6b, create an electric field between the layers, which in turn orders water dipoles toward the vapor phase. The reversed orientation of water dipoles is one of the remarkable consequences of the solute perturbation.83,84



TABLE 6: Calculated 10-90 Thickness, δ10-90, and Surface Potential, O, at the Interfaces of Pure Water and Aqueous Solutions δ10-90 (Å) φ (mV)

pure water

1.1 M NaCl

2.1 M NaCl

1.1 M NaI

2.1 M NaI

3.14 -421.2 ( 0.8

2.98 -420.2 ( 3.5

2.89 -415.9 ( 2.7

3.29 -390.5 ( 2.0

3.24 -376.0 ( 3.5

Figure 7. Definition of orientational angles (a) θdip and (b) φ for a surface water molecule.

Figure 8. Average distribution of 〈cos θdip(zˆ)〉 for pure water and NaCl and NaI aqueous solutions.

Figure 9 illustrates the full orientational distribution in the two-dimensional space (θdip, φ) with a depth resolution at 1.5 Å spacing. In this figure, we display the anisotropic part

P′(cos θdip, φ; zˆ) ≡ P(cos θdip, φ; zˆ) - Pisotropic

(26)

where Pisotropic is the isotropic probability, Pisotropic ) 1/(4π). Accordingly, the positive P′ regions, colored blue in Figure 9, denote higher probabilities than the isotropic average, while the negative P′ regions, colored red, denote lower probabilities than the average. Typical configurations with high probabilities are labeled a-f, and they are illustrated in Figure 10. Figure 9 displays φ in the range 0-180°, though φ < 90° and φ > 90° should be symmetric. In the vapor side of the Gibbs dividing surface, 0.0 Å < zˆ < 3.0 Å, for pure water, 2.1 M NaCl, and 2.1 M NaI aqueous solution surfaces, two major orientations, labeled a and b, are commonly observed in Figure 9 and illustrated in Figure 10. In the a state, a water molecule directs one OH bond almost vertically to the vapor phase and the other to the liquid phase. The a state is thought to be characteristic of the dangling OH and relatively emphasized in the outermost layer (1.5 Å < zˆ < 3.0 Å). In the b state, the molecular plane is close to parallel to the interface but with both OH bonds slightly toward the liquid phase. This configuration is relatively significant in the vicinity of the Gibbs dividing surface (0 Å < zˆ < 1.5 Å). In the liquid side of the Gibbs dividing surface (-1.5 Å < zˆ < 0.0 Å), the c state, with the molecular plane parallel to the

interface, is dominant commonly in pure water and aqueous solutions. This state is modified with the dipole toward the vapor phase to a slight extent (d state) in the region of -3.0 Å < zˆ < -1.5 Å. The b, c, and d states are commonly featured with cos θdip ∼ 0 and φ ∼ 90°, which means that the molecular plane is nearly lying on the interface. Such an orientational configuration is regarded as dominant at the gas-water interface.24,81 We also notice a small variation in the maximum location of cos θdip among the three states b-d. As depicted in Figure 10, the maximum of cos θdip is slightly positive for b in the vapor side, close to zero for c, and slightly negative for d in the liquid side. In the region of -3.0 Å < zˆ < -1.5 Å, another orientational state, e, emerges in pure water and NaCl solution systems, whereas this e state is missing in the NaI solution. The e state indicates a water molecule with one OH bond almost vertical to the liquid phase. These two states d and e in one monolayer region below the Gibbs dividing surface may stem from the hydrogen bonding network, whereas the network in NaI solution may be somewhat perturbed due to the strong enhancement of iodide anion concentration or the electric double layer formation. In the deeper region of -6.0 Å < zˆ < -3.0 Å, the water orientation is close to isotropic in pure water and NaCl solution, whereas the water dipoles in the NaI solution retain strong orientation toward the vapor phase (f state). This orientation is apparently due to the electric double layers of Na+ and I-. (We note in Figure 6 that the maximum densities of Na+ and I- are located at about zˆ ∼ - 6 and -2 Å, respectively, in 2.1 M NaI solution.) To summarize the above results, water orientation is hardly perturbed by ions in the vapor side of the Gibbs dividing surface (0.0 Å < zˆ < 3.0 Å), where a common state a is observed among pure water and the electrolyte solutions. Near the Gibbs dividing surface, the b, c, and d states, with the molecular plane nearly parallel to the interface, are dominant commonly among the three systems. However, the strong enhancement of anion concentration in the NaI system gradually disturbs the hydrogen bonding network in the liquid side of the Gibbs dividing surface (-3.0 Å < zˆ < 0.0 Å), and therefore, the NaI system loses the e state which is observed in pure water and NaCl solution. The most characteristic orientation of NaI solution emerges in the region zˆ < -3 Å, where the water molecules direct their dipoles toward the vapor phase in response to the electric double layer formation. 4.3. Induced Dipole Moment Profiles. In section 4.2, we have elucidated the orientational order of the water permanent dipole moment. On the other hand, the objective of this section is to investigate the induced dipole moments of water molecules and ions. Figure 11 displays the average induced dipole, 〈|µ|〉 (panel a), and its z component, 〈µz〉 (panel b), of a water molecule as a function of the vertical depth, zˆ. In the bulk liquid, the average induced dipole takes a finite value, 〈|µ|bulk〉 ∼ 0.6 D, owing to the reaction field, while |µ| should vanish in the gas phase. Panel a shows the gas-liquid transition behavior of the induced dipole |µ|. This transition behavior of |µ| poses a considerable difference to that of the density profile, F, in Figure 6, though both quantities decay in the vapor side. The water density, F, changes at around the Gibbs dividing surface relatively abruptly


Ishiyama and Morita

Figure 9. Orientational probability distributions of water molecules, P′(cos θdip, φ; z), for pure water (left column), 2.1 M NaCl (center), and 2.1 M NaI (right) aqueous solutions. The blue regions refer to higher probabilities than average (P′ > 0), while the red regions, to lower probabilities (P′ < 0). The major orientations labeled a-f are illustrated in Figure 10.

by eq 25 with a 10-90 thickness of δ10-90 ∼ 3 Å, whereas the induced dipole of water decays significantly slower for the vapor side. The difference in the transition behaviors is manifested at zˆ ) 3 Å, for example, where 〈|µ|〉/〈|µ|bulk〉 = 0.6 while F/Fbulk = 0.03. A water molecule in the vapor side, which is connected to the liquid by a limited number of hydrogen bonds, can retain a substantial magnitude of the induced dipole. Regarding the z component of the induced dipole moment, µz, on the other hand, Figure 11b shows a remarkable maximum in µz at zˆ ∼ -4 Å, where the substantial electric field is present due to the double layers of ions. However, the ratio of µz to |µ| in panels a and b, which roughly corresponds to 〈cos θdip〉 in Figure 8a, is much smaller than unity, implying that the induced dipoles of water are not fully ordered on average along the surface normal even in the electric double layer region. Figure 12 shows the average induced dipole, |µ|, and its z component, µz, of anions per molecule. The polarization of the cation (Na+) is omitted here because of its negligible role compared to the anions. Figure 12 clearly indicates a largely

induced dipole moment near the interface, particularly for the iodide anion (red lines), with a maximum location in the vapor side of the Gibbs dividing surface. The average induced dipole of ions vanishes at zˆ < 5 Å in the isotropic bulk region. An important feature in Figure 12 is that the magnitude of µz is comparable to |µ| near the interface, indicating that the anions direct their induced dipoles almost toward the surface normal from the liquid phase to the vapor phase. It is also interesting to see that the induced dipole moment of an anion in the 1.1 M system is larger than that in the 2.1 M system. This is explained by the fact that the anion-anion dipolar interaction near the interface tends to destabilize the induced dipoles along the surface normal.23,84 4.4. Surface Potential. The surface potential reflects the electrostatic distribution at the interface, and therefore is considered a good indicator of molecular orientation at the interface. While it has been a challenging issue for a long time to determine the absolute value of the surface potential at the gas-vapor interface of water, the effects of inorganic salts on



Figure 10. Schematic of the major water orientations a-f depicted from the results of Figure 9.

Figure 12. Average induced dipole moment profiles of an anion: (a) absolute, |µ|; (b) z component, µz.

∆φ ) φ(z1) - φ(z0) ) -

∫zz dz′ Ez(z′) 1

(27)

0

where Ez(z) is the normal component of the electric field at z. Ez(z) is related to the charge density, Fq(z), and the induced dipole density, Fµz(z), as

0Ez(z) + Fµz(z) )

∫zzdz′ Fq(z′)

(28)

0

where 0 is the dielectric permittivity of free space. From eqs 27 and 28, the surface potential, ∆φ, is given as

∆φ ) ∆φq + ∆φµz

(29)

where

∆φq ) -

Figure 11. Average induced dipole moment profiles of a water molecule: (a) absolute, |µ|; (b) z component, µz.

the surface potential have been well examined experimentally.5,85 Here, we calculate the surface potential of the aqueous solutions to assess the performance of the present force field on these interfacial properties. First, we briefly formulate the surface potential, ∆φ, which is an electrostatic potential difference between the vapor and the bulk liquid. In the case of a plane interface, the reference points z0 and z1 are put far from the interface in the vapor phase and liquid phase, respectively. Then, ∆φ can be defined as12,14,82

1 0

∫zz dz ∫zzdz′ Fq(z′), 1

0

0

∆φµz )

1 0

∫zz dz Fµ (z) 1

0

z

(30)

In the polarizable MD simulations, Fq(z) and Fµz(z) can be easily evaluated. Figure 13a shows the average net (positive + negative) charge density distributions, Fq(z), for pure water and electrolyte solutions. All of the profiles in panel a commonly exhibit a positive region of Fq(z) in the vapor side of the Gibbs dividing surface and a conspicuous negative region in the adjacent layer in the liquid side, implying a polarization of partial charge distribution of water near the Gibbs dividing surface.82 Besides these common features among pure water and electrolyte solutions, influence of the segregated ion distributions is obviously seen in the charge distribution, Fq(z), of the NaI solutions, particularly the 2.1 M solution denoted by a solid red line. The negative region at about zˆ ∼ -2 Å becomes more pronounced due to the I- layer, and another positive region


Ishiyama and Morita

Figure 13. (a) Charge density, Fq(zˆ), and (b) induced dipole density, Fµz(zˆ), for pure water and NaCl and NaI aqueous solutions.

emerges at about zˆ ∼ -5 Å due to the Na+ layer. Panel b indicates that the net induced dipole density exhibits a large maximum at about zˆ ∼ -3 Å for the electrolyte solutions, particularly for the NaI solutions, in response to the electric double layer formation. Using Fq(z) and Fµz(z) in Figure 13, the electrostatic potential difference between vapor and liquid (surface potential) can be calculated by eqs 29 and 30. The results of ∆φ and their components (∆φq and ∆φµz) are shown in Figure 14 as a function of z1. Note that ∆φ(z1) coincides with the surface potential, ∆φ, when z1 is set far from the interface in the liquid (zˆ ) z1 f -∞). In the pure water (panel a), the asymptotic value of ∆φ(z1) is estimated to be -0.421 V, which is in accord with previous calculations with various force fields in the range -0.9 to 0 V (including quadrupole contribution).82 In panel a, the surface potential, ∆φ, of pure water is dominated by ∆φq (red line), along with a slight positive contribution of ∆φµz (blue line). This implies that the surface potential is mostly governed by the ∆φq value of pure water. The sign of ∆φµz appears opposite to the calculation by Dang and Chang,12 although the reported value of ∆φ ) -0.50 ( 0.001 V is nearly consistent with the present result. The sign of ∆φµz is determined by integrating the Fµz(z) profile in Figure 13b, as a consequence of cancellation of positive and negative regions. ∆φµz appears to be more sensitive to the force field than ∆φq, though the former is of secondary importance in the total ∆φ value for pure water. In the case of aqueous solutions in Figure 13b and c, however, ∆φµz plays an equally important role as ∆φq in the surface potential. From zˆ ) 0 Å to -5 Å, ∆φq (red lines) rises toward the liquid phase due to the electric double layer of ions, whereas ∆φµz (blue lines) drops so as to cancel out the above contribution of ∆φq, resulting in a nearly flat total surface potential, ∆φ (black lines). A similar conclusion is also reached in a recent simulation by Wick et al.86 As a consequence of cancellation,

Figure 14. Calculated surface potential, ∆φ(z1), via eq 27 and the components ∆φq(z1) and ∆φµz(z1) via eq 30 as a function of zˆ ) z1 for (a) pure water and (b) NaCl and (c) NaI solutions.

the surface potential shows a little change between aqueous solutions and pure water. The calculated difference, ∆φ, for the aqueous solutions to that for the pure water, ∆φ(solution) ∆φ(water), are plotted in comparison with the experiment values85 in Figure 15. The present simulation with the FP model reproduces quite well the experimental results of ∆φ(solution) - ∆φ(water). 4.5. Surface Excess. It has been argued whether the surface excess estimated from experiments is in accord with the density profiles obtained by MD simulation. In principle, the Gibbs surface excess, Γi, for i species can be evaluated from the calculated density profile, Fi(zˆ), as10,87

Γi )

∫z0[Fi(zˆ) - Fi,liq]dzˆ + ∫0z [Fi(zˆ) - Fi,gas]dzˆ 0

(31)

1

where Fi,gas and Fi,liq are the bulk gas and liquid densities for i species. From the definition of the Gibbs dividing surface, Γi vanishes for the solvent water and Fi,gas ) 0 for the nonvolatile ions. For 1:1 electrolyte solutions such as NaCl and NaI, the electroneutrality condition results in Γ+ ) Γ-(≡Γ) for cations and anions.88 In the present study, we calculated Γ from the density profile (see Figure 6) via eq 31 in Table 7, where the



Figure 15. Comparison of the calculated difference in the surface potential between aqueous solutions and pure water, ∆φ(solution) ∆φ(water), with the experiments at 20 °C.85

< 1 for I- is apparently observed in -9 Å j zˆ j -5 Å. This cancellation mechanism certainly plays an important role in the calculation of Γ. Besides this cancellation, the present force field model tends to produce a negative Γ value, partly because the double layer structure of ions is formed somewhat to the liquid side of the Gibbs dividing surface, as discussed in section 4.1. However, we notice that considerable discrepancies remain between the calculated Γ values and experimental ones in Table 7. The surface excess, Γ, appears to be very sensitive to the density profiles and thus a difficult quantity to reproduce. Radke and co-workers examined the surface tension and the surface excess for NaF, NaCl, and NaBr aqueous solutions by nonpolarizable MD simulations88,93 and reported that those quantities deviate from the experimental ones particularly at high concentration.93 In their simulations, the surface excess for NaCl solutions changes from negative to positive with increasing salt concentration and that for NaBr solutions shows a somewhat oscillating behavior,93 which is clearly inconsistent with the experimental measurements. As was pointed out by them, the system size dependence of Γ may also be an important factor to be examined. To reproduce Γ by MD simulation is an important subject for future studies. 5. Concluding Remarks

Figure 16. Ion surface excess, Γexp, derived from the experimental surface tension at 25 °C of NaCl and NaI solutions91 and activity coefficients.92 The solid lines denote approximated surface excess, -(x/ RT) dγ/dx, where x is the mole fraction of salt.

error range is estimated from numerical discrepancy of Γi for cations and anions in the two sides of the liquid slab. The surface excess of ions, Γ, is related to the concentration dependence of the surface tension, γ, by the following thermodynamic condition:89

Γexp ) -

a dγ RT da

(32)

where a is the activity of ions. To distinguish the two definitions, the surface excess derived via eq 32 is labeled Γexp.90 Figure 16 displays the experimental Γexp values derived by the available data of surface tension91 and activity coefficients.92 As stated in section 1, the surface tension measurements lead to Γexp < 0 in Figure 16, owing to the fact that the surface tension increases with ion concentration. The calculated Γ values by eq 31 are compared to the experimental Γexp values by eq 32 at the same concentrations in Table 7. The calculated Γ values tend to predict negative surface excess except for the case of the 1.1 M NaI solution, though the large error range makes it hard to determine the sign of Γ. Vrbka et al. argued87 that the negative Γ could be compatible with the surface segregation of ions, since the accumulated ion concentration is canceled by the presence of an ion deficient region beneath the ion rich layer. For the 2.1 M NaI solution in Figure 6b, for example, the region of F/Fbulk

In the present work, we developed a flexible and polarizable water model, particularly for describing the interfacial properties of aqueous systems. The internal vibration is incorporated for purposes of properly describing intramolecular O-H stretching vibrations in the condensed phase, which are probed by interfacial sum frequency generation spectroscopy as well as the infrared and Raman spectroscopies. While the effects of polarizability are well recognized on the interface modeling of electrolyte aqueous solutions, simultaneous effects of internal vibration could be expected to have additional impact on interfacial properties. The present model of water was thoroughly examined by calculating representative properties associated with bulk liquid and the gas-liquid interface of pure water and NaCl and NaI aqueous solutions. These properties include density, vaporization energy, dipole moment, diffusion coefficient, radial distribution function, infrared and Raman spectra of the O-H stretching region, surface potential, and surface excess of ions. The calculated results demonstrate that the present model has equivalent performance to commonly used polarizable molecular models in describing bulk and interface properties not directly pertinent to intramolecular vibrations. The infrared and Raman spectra for intramolecular vibrations are reasonably described by classical MD simulation with the present model, which empirically takes account of anharmonicity and frequency shift in the condensed phase. One of the noticeable features of the present model is relatively weaker electrostatic interaction of ions in a shortrange region than other rigid polarizable models. This feature is manifested in the calculated ion-ion radial distribution functions in the bulk solutions or ion density profiles near the gas-liquid interface. This difference is mainly attributed to the short-range damping treatment rather than a direct consequence of the internal vibrations. The molecular flexibility rather augments practical importance of the short-range damping treatment to avoid the unphysical divergence of polarization. While the short-range correction is invoked to account for the well-known shortcomings of conventional interaction site models with the point charge/polarizability approximation, the


Ishiyama and Morita

TABLE 7: Calculated and Experimental Surface Excess of Ions, Γ Γ (nm-2) calcd Γexp (nm-2) exptla a

1.1 M NaCl

2.1 M NaCl

1.1 M NaI

2.1 M NaI

-0.101 ( 0.086 -0.461

-0.087 ( 0.075 -0.758

0.077 ( 0.219 -0.255

-0.038 ( 0.066 -0.396

At 20 °C.91

correction is in itself of empirical nature. The reliability of the correction will be an important issue to be further examined. In connection to this issue, the calculated density profile of iodide anion for NaI solution interfaces exhibits a maximum slightly to the liquid side of the Gibbs dividing surface. This implies that the water structure of the topmost layer is not significantly perturbed by the ions, irrespective of the electric double layer formation of Na+ and I- near the interface. Another feature worth noting is in the calculation of ion surface excess, which the present model apparently fails to reproduce. Although the present model tends to predict a small or negative surface excess of ions, its magnitude severely underestimates the experimental results. Further investigation of surface excess is a remaining crucial issue to understand the gas-liquid interfacial structure of electrolyte solutions. The present model of water and detailed information on the interface structure are utilized in the following paper, which mainly discusses the relation between the vibrational sum frequency generation spectra and microscopic interface structure. Acknowledgment. The authors would like to thank Prof. Susumu Okazaki, Dr. Noriyuki Yoshii, and Dr. Tateki Ishida for valuable discussions. This work was supported by the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan. Appendix: Ewald Sum and Short-Range Damping Treatment

UC ) Upol ) -

∑i

µiE0i

+

1

qaiqbj ∑∑ r a,b

2 i,j(*i)

∑

2 i,j(*i)

µiTijµj +

1 2

1 2

∑i

µiE0i

The derivatives of these energies with respect to site c of molecule k, which yield the force terms acting on the site, are written as

FCck ) -∇ckUC ) qckE0ck

(35)

) -∇ckUpol )

∑i µi(∇ckE0i ) - ∑ µi(∇ckTik)µk + 2 (∇ckrk)EkEk i(*k)

(36a)

5

[A1ik + A2ik]

(36c)

where

A1ik ) (µiµk)rik + (rikµi)µi + (rikµk)µk A2ik ) -

5 (rikµi)(rikµk)rik rik2

(37a) (37b)

In eqs 36b and c, mck is the mass of site c of molecule k and Mk k ) ∑on mc is the mass of molecule k. The center of mass is c k k located at rk ) ∑on mckrck/∑on mck, which leads to ∇ck ) c c (mck/Mk)∇k. Next, we consider a three-dimensional periodic system. Here, the short-range damping is not employed, assuming f (1) ) f (2) ) w ) 1 in eqs 3, 5, and 7. Implementation of the damping functions is discussed later. For the periodic system, the terms of Coulomb energy, UC, electric field, E0i (E0ai), dipole tensor, Tij (Tai,j), and its derivative, ∇kTik, in eqs 33-36 are replaced with the following expressions:

1 2

′∑qai qbj φ(rai,bjn) ∑n ∑ i,j a,b

(38a)

∑n ∑j ′∑b qbj∇iφ(ri,bjn)

(38b)

∑n ∑j ′∑b qbj∇aiφ(rai,bjn)

(38c)

E0i ) E0ai ) -

∇kTik ) (34)

3

(36b)

k i(*k) rik

i(*k)

Tai,j )

∑i EiriEi )

1

mck

(33)

ai,bj

k i(*k)

∑ µi(∇ckTik)µk ) M ∑

Tij )

1

-

Fpol ck

i(*k)

UC )

In this appendix, we detail the present computation for Coulombic and induced dipolar interactions with use of the Ewald sum and the short-range damping treatment. Energy and Force. First, we give general expressions of energy and force terms before introducing Ewald and damping functions. The Coulomb and polarization energies, UC and Upol, are presented as

1

mck

∑i µi(∇ckE0i ) ) - ∑ qckTck,iµi + M ∑ ∑a qaiTai,kµk

∑n ∇i∇jφ(rijn)

(38d)

∑n ∇ai∇jφ(rai,jn)

(38e)

∑n ∇k∇i∇kφ(rikn)

(38f)

where φ(r) ) 1/r and rai,bjn ) |rai,bjn| ) |rai - rbj + (Lxnx, Lyny, Lznz)|. n ) (nx, ny, nz) is the integer vector specifying lattice. The prime sign in the summation denotes that the term with i ) j and n ) 0 is omitted. Now, we introduce two sets of auxiliary functions, Bi and Ci, for use of later expressions:

2R -(Rr)2 r e xπ

(39a)

4R3 -(Rr)2 3 e r 3xπ

(39b)

B1(r) ) erfc(Rr) + B2(r) ) B1(r) +


8R5 -(Rr)2 5 B3(r) ) B2(r) + r e 15xπ 2R -(Rr)2 e C1(r) ) erf(Rr) r xπ C2(r) ) C1(r) -

4R3 -(Rr)2 3 e r 3xπ

J. Phys. Chem. C, Vol. 111, No. 2, 2007 735 where

[

(39c)

∑n ′

Treal ij )

(40a)

rijn3

4π

Trec ij )

(40b)

Urec C

)

2π V

1

∑∑ ∑ 2 n i,j a,b ′

erfc(Rrai,bjn) qaiqbj

rai,bjn

∑ F(|G|)∑ ∑qaiqbj cos[G(rai - rbj)] G*0 i,j a,b Uself C

)Umask C

R

)-

∑i ∑a qai

2

xπ

1

∑∑ ∑ 2 i a

b(*a)

real ) Tai,j

E0,rec ) i

4π V

)

(46a)

(46b)

4R3 Iδij 3xπ

(46c)

∑n

[

I

rai,jn3

B1(rai,jn) - 3

]

rai,jnrai,jn rai,jn5

B2(rai,jn)

(41d)

rec ) Tai,j

4π V

mask )Tai,j

[

∑ GGF(|G|) cos[G(rai - rj)]

ri,bjn3

B1(ri,bjn)

∑ GF(|G|)∑j ∑b qbj sin[G(ri - rbj)] G*0 )0 E0,self i

I 3

rai,j

C1(rai,j) - 3

]

rai,jrai,j

C2(rai,j) δij (47c)

rai,j5

In eqs 46c and 47c, δij is Kronecker’s delta. For the derivative of the dipole tensor, ∇ckTik ) (mck/Mk)∇kTik is included in eq 36a or 36c as the form of µi(∇kTik)µk, which consists of the following real, reciprocal, and self terms:

µi(∇kTreal ik )µk ) -

)E0,mask i

(42a) µi(∇kTrec ik )µk )

∑n

3

rikn5

[A1iknB2(rikn) + A2iknB3(rikn)] (48a)

qairi,ai

∑a

ri,ai3

C1(ri,ai)

(42d)

For the electric field E0ai, the analogous expressions for E0,real ai , 0,self E0,rec , and E are derived from eqs 42a-c, respectively, by ai ai , is given replacing the suffix i with ai. The mask term, E0,mask ai as

)E0,mask ai

∑

b(*a)

qbirai,bi rai,bi3

C1(rai,bi)

4π V

∑ G(µiG)(µkG)F(|G|) sin[G(ri - rk)]

G*0

(48b)

µi(∇kTikself)µk ) 0

(42b) (42c)

(47b)

G*0

qbjri,bjn

∑n ∑j ′∑b

(47a)

(41b)

In eq 41b, F(R) ) exp[-R2/(4R2)]/R2. G ) 2π(hx/Lx, hy/Ly, hz/ Lz) is the reciprocal lattice vector, and (hx, hy, hz) is the integer vector specifying the reciprocal lattice. R is called the Ewald separation parameter. The mask term, Umask C , is required for a molecular system to account for the self-interaction among the intramolecular sites. For the electric field E0i ,

E0,real i

]

and

(41c)

rai,bi

B2(rijn)

rijn5

G*0

Tself ij )

(41a)

erf(Rrai,bi) qaiqbi

B1(rijn) - 3

∑ GGF(|G|) cos[G(ri - rj)]

V

Ewald Expressions. The Ewald expression for the UC potential is presented as the sum of the real part, Ureal C , self mask reciprocal part, Urec C , self part, UC , and mask part, UC :

Ureal C )

rijnrijn

I

(48c)

Note that the dipole tensor Tij in the above Ewald expressions formally includes the self-interaction (i ) j), and thus, the summation relevant to the dipole tensor should be modified accordingly. This means that ∑j(*i) in eq 6 should be replaced with ∑j using the Ewald summation. Damping Functions. Then, we apply the short-range damping functions, f (1), f (2), and w in eqs 3, 5, and 7, to the above Ewald expressions. The damping functions in the short range are invoked to avoid divergence of polarization. For this purpose, it is sufficient and convenient to apply them to the real part of the Ewald expressions. In the present study, therefore, eqs 41a, 42a, 46a, and 48a are modified as

(43) ) Ureal/ C

For the dipole tensors Tij and Tai,j, rec self Tij ) Treal ij + Tij + Tij

(44)

real rec mask + Tai,j + Tai,j Tai,j ) Tai,j

(45)

1

(1) ′∑ fai,bj (rai,bjn)qaiqbj ∑ ∑ 2 n i,j a,b

) E0,real/ i

∑n ∑j ′∑b fi,bj(2)(ri,bjn)

erfc(Rrai,bjn) rai,bjn

(49)

qjri,bjn ri,bjn3

B1(ri,bjn)

(50)


Ishiyama and Morita

Treal/ ) ij

∑n

[

I

rijn

3

B1(rijn)(4wijn3 - 3wijn4) - 3

µi(∇kTikreal/)µk )

∑n

[

(rijnrijn) rijn5

B2(rijn)wijn4

]

(51) -3

rikn5

{A1iknB2(rikn)(4wikn3 - 3wikn4) +

A2iknB3(rikn)wikn4} + 12(∇kwikn) wikn3) -

{

(µiµk) rikn3

(µirikn)(µkrikn) rikn5

B1(rikn)(wikn2 -

B2(rikn)wikn3

}]

(52)

(1) (2) where fai,bj (rai,bjn), fi,bj (ri,bjn), and wijn ) wij(rijn) are defined in eqs 3, 5, and 7, respectively.

References and Notes (1) Finlayson-Pitts, B. J.; Pitts, J. N. J. Chemistry of the Upper and Lower Atmosphere; Academic Press: San Diego, CA, 2000. (2) Hu, J. H.; Shi, Q.; Davidovits, P.; Worsnop, D. R.; Zahniser, M. S.; Kolb, C. E. J. Phys. Chem. 1995, 99, 8768. (3) Seinfeld, J. H. Science 2000, 288, 285. (4) Knipping, E. M.; Lakin, M. J.; Foster, K. L.; Jungwirth, P.; Tobias, D. J.; Gerber, R. B.; Dabdub, D.; Finlayson-Pitts, B. J. Science 2000, 288, 301. (5) Randles, J. E. B. Phys. Chem. Liq. 1977, 7, 107. (6) Landau, L. D.; Lifshitz, E. M. Electrodynamics of continuous media; Course of Theoretical Physics, Vol. 8; Pergamon Press: New York, 1960. (7) Onsager, L.; Samaras, N. N. T. J. Chem. Phys. 1934, 2, 528. (8) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2001, 105, 10468. (9) Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2002, 106, 6361. (10) Jungwirth, P.; Tobias, D. J. Chem. ReV. 2006, 106, 1259. (11) Archontis, G.; Leontidis, E.; Andreou, G. J. Phys. Chem. B 2005, 109, 17957. (12) Dang, L. X.; Chang, T. M. J. Phys. Chem. B 2002, 106, 235. (13) Dang, L. X. J. Phys. Chem. B 2002, 106, 10388. (14) Chang, T. M.; Dang, L. X. Chem. ReV. 2006, 106, 1305. (15) Benjamin, I. J. Chem. Phys. 1991, 95, 3698. (16) Wilson, M. A.; Pohorille, A. J. Chem. Phys. 1991, 95, 6005. (17) Petersen, P. B.; Saykally, R. J. Anuu. ReV. Phys. Chem. 2006, 57, 333. (18) Petersen, P. B.; Johnson, J. C.; Knutsen, K. P.; Saykally, R. J. Chem. Phys. Lett. 2004, 397, 46. (19) Petersen, P. B.; Saykally, R. J. Chem. Phys. Lett. 2004, 397, 51. (20) Petersen, P. B.; Saykally, R. J. J. Am. Chem. Soc. 2005, 127, 15446. (21) Ghosal, S.; Hemminger, J. C.; Bluhm, H.; Mun, B. S.; Hebenstreit, E. L. D.; Ketteler, G.; Ogletree, D. F.; Requejo, F. G.; Salmeron, M. Science 2005, 307, 563. (22) Weber, R.; Winter, B.; Schmidt, P. M.; Widdra, W.; Hertel, I. V.; Dittmer, M.; Faubel, M. J. Phys. Chem. B 2004, 108, 4729. (23) Ishiyama, T.; Morita, A. J. Phys. Chem. C, 2007, 111, 738. (24) Morita, A.; Hynes, J. T. Chem. Phys. 2000, 258, 371. (25) Morita, A.; Hynes, J. T. J. Phys. Chem. B 2002, 106, 673. (26) Wallqvist, A. Chem. Phys. 1990, 148, 439. (27) Zhu, S. B.; Singh, S.; Robinson, G. W. J. Chem. Phys. 1991, 95, 2791. (28) Corongiu, G. Int. J. Quantum Chem. 1992, 42, 1209. (29) Saint-Martin, H.; Hernandez-Cobos, J.; Bernal-Uruchurtu, M. I.; Ortega-Blake, I.; Berendsen, H. J. C. J. Chem. Phys. 2000, 113, 10899. (30) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (31) Iuchi, S.; Morita, A.; Kato, S. J. Phys. Chem. B 2002, 106, 3466. (32) Jeon, J.; Lefohn, A. E.; Voth, G. A. J. Chem. Phys. 2003, 118, 7504. (33) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (34) Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 1996, 105, 8274. (35) Thole, B. T. Chem. Phys. 1981, 59, 341. (36) Bernardo, D. N.; Ding, Y.; Krogh-Jespersen, K.; Levy, R. M. J. Phys. Chem. 1994, 98, 4180. (37) Ichikawa, K.; Kameda, Y.; Yamaguchi, T.; Wakita, H.; Misawa, M. Mol. Phys. 1991, 73, 79. (38) Moriarty, N. W.; Karlstro¨m, G. J. Chem. Phys. 1997, 106, 6470.

(39) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572. (40) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (41) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 86th ed.; CRC Press: Boca Raton, FL, 2005. (42) Marti, J.; Padro, J. A.; Guardia, E. J. Mol. Liq. 1994, 62, 17. (43) Dyke, T. R.; Muenter, J. S. J. Chem. Phys. 1973, 59, 3125. (44) Ahlstro¨m, P.; Wallqvist, A.; Engstro¨m, S.; Bojo¨nsson, B. Mol. Phys. 1989, 68, 563. (45) Stillinger, F. H. In The Liquid State of Matter; Fluids, Simple and Complex; Montroll, E. W., Lebowitz, J. L., Eds.; North-Holland: Amsterdam, The Netherlands, 1982. (46) Toukmaji, A.; Sagui, C.; Board, J.; Darden, T. J. Chem. Phys. 2000, 113, 10913. (47) Aguado, A.; Madden, P. A. J. Chem. Phys. 2003, 119, 7471. (48) Fincham, D. Mol. Simul. 1994, 13, 1. (49) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, U.K., 1987. (50) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (51) Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. J. Chem. Phys. 2003, 119, 5185. (52) Dore, J. C.; Garawi, M.; Bellissent-Funel, M. C. Mol. Phys. 2004, 102, 2015. (53) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47. (54) Soper, A. K. J. Chem. Phys. 1994, 101, 6888. (55) Chowdhuri, S.; Chandra, A. J. Chem. Phys. 2001, 115, 3732. (56) Paul, S.; Chandra, A. Chem. Phys. Lett. 2003, 373, 87. (57) Lyubartsev, A. P.; Laaksonen, A. J. Phys. Chem. 1996, 100, 16410. (58) Uchida, H.; Matsuoka, M. Fluid Phase Equilib. 2004, 219, 49. (59) Perera, L.; Essmann, U.; Berkowitz, M. L. J. Chem. Phys. 1995, 102, 450. (60) Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. (61) Yonetani, Y. Chem. Phys. Lett. 2005, 406, 49. (62) Steinhauser, O. Mol. Phys. 1982, 45, 335. (63) Dang, L. X. J. Chem. Phys. 1992, 97, 1919. (64) Koneshan, S.; Rasaiah, J. C. J. Chem. Phys. 2000, 113, 8125. (65) Koneshan, S.; Rasaiah, J. C.; Lynden-Bell, R. M.; Lee, S. H. J. Phys. Chem. B 1998, 102, 4193. (66) Guardia, E.; Rey, R.; Padro, J. A. Chem. Phys. 1991, 155, 187. (67) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (68) Peslherbe, G. H.; Ladanyi, B. M.; Hynes, J. T. J. Phys. Chem. A 2000, 104, 4533. (69) McQuarrie, D. A. Statistical Mechanics; Harper Collins Publishers: New York, 1976. (70) Bader, J. S.; Berne, B. J. J. Chem. Phys. 1994, 100, 8539. (71) Bertie, J. E.; Lan, Z. Appl. Spectrosc. 1996, 50, 1047. (72) Egorov, S. A.; Everitt, K. F.; Skinner, J. L. J. Phys. Chem. A 1999, 103, 9494. (73) Brooker, M. H.; Hancock, G.; Rice, B. C.; Shapter, J. J. Raman Spectrosc. 1989, 20, 683. (74) Buch, V. J. Phys. Chem. B 2005, 109, 17771. (75) Press, W. H.; Teukolski, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran77: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: New York, 1992. (76) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (77) Garrett, B. C.; Schenter, G. K.; Morita, A. Chem. ReV. 2006, 106, 1355. (78) 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. (79) Jedlovszky, P.; Vincze, A.; Horvai, G. J. Chem. Phys. 2002, 117, 2271. (80) Pa´rtay, L.; Jedlovszky, P.; Vincze, A.; Horvai, G. J. Phys. Chem. B 2005, 109, 20493. (81) Matsumoto, M.; Kataoka, Y. J. Chem. Phys. 1988, 88, 3233. (82) Sokhan, V. P.; Tildesley, D. J. Mol. Phys. 1997, 92, 625. (83) Brown, E. C.; Mucha, M.; Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2005, 109, 7934. (84) Ishiyama, T.; Morita, A. Chem. Phys. Lett., 2006, 431, 78. (85) Jarvis, N. L.; Scheiman, M. A. J. Phys. Chem. 1968, 72, 74. (86) Wick, C. D.; Dang, L. X.; Jungwirth, P. J. Chem. Phys. 2006, 125, 024706. (87) Vrbka, L.; Mucha, M.; Minofar, B.; Jungwirth, P.; Brown, E. C.; Tobias, D. J. Curr. Opin. Colloid Interface Sci. 2004, 9, 67. (88) Bhatt, D.; Newman, J.; Radke, C. J. J. Phys. Chem. B 2004, 108, 9077. (89) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces, 6th ed.; Wiley: New York, 1997. (90) A factor of 1/2 is not included in eq 32, as Γexp means the surface excess of each ion component.89

Flexible and Polarizable Molecular Modeling (91) Washburn, E. W., Ed. International Critical Tables; McGrawHill: New York, 1928. (92) Robinson, R. A.; Stokes, R. H. Electrolyte Solutions, 2nd ed.; Butterworth: London, 1959. (93) Bhatt, D.; Chee, R.; Newman, J.; Radke, C. J. Curr. Opin. Colloid Interface Sci. 2004, 9, 145. (94) Chowdhuri, S.; Chandra, A. J. Chem. Phys. 2003, 118, 9719. (95) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1991, 95, 1954. (96) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1994, 100, 3085.

J. Phys. Chem. C, Vol. 111, No. 2, 2007 737 (97) Markovich, G.; Perera, L.; Berkowitz, M. L.; Cheshnovsky, O. J. Chem. Phys. 1996, 105, 2675. (98) Wagner, W.; Pruss, A. J. Phys. Chem. Ref. Data 2002, 31, 387. (99) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys 1983, 79, 926. (100) Oki, M., et al., Ed. Kagaku Binran (in Japanese), 5th ed.; Maruzen: Tokyo, 2004. (101) Badyai, Y. S.; Saboungi, M. L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, A. K. J. Chem. Phys. 2000, 112, 9206.

Molecular Dynamics Study of GasLiquid Aqueous Sodium Halide

Recommend Documents