Article pubs.acs.org/JPCB
Concentration Dependence of the Dielectric Permittivity, Structure, and Dynamics of Aqueous NaCl Solutions: Comparison between the Drude Oscillator and Electronic Continuum Models Richard Renou,† Minxia Ding,† Haochen Zhu,‡ Anthony Szymczyk,† Patrice Malfreyt,§ and Aziz Ghoufi*,∥ †
Institut des Sciences Chimiques de Rennes, CNRS, UMR 6226, Université de Rennes 1, 263 Avenue du Général Leclerc, 35042 Rennes, France ‡ State Key Laboratory of Pollution Control and Resources Reuse, College of Environmental Science and Engineering, Tongji University, Mingjing Building, 1239 Siping Road, Shanghai 200092, China § Clermont Université, Université Blaise Pascal, Institut de Chimie de Clermont-Ferrand, BP 10448, F-63000 Clermont-Ferrand, France ∥ Institut de Physique de Rennes, CNRS, UMR 6251, Université Rennes 1, 263 Avenue Général Leclerc, 35042 Rennes, France ABSTRACT: We report molecular dynamics simulations of aqueous sodium chloride solutions at T = 298 K and p = 1 bar in order to investigate the salt concentration dependence of the dielectric permittivity, the structure, and the dynamical properties. Different models were applied up to 7 m salt concentration: the Drude oscillator model with a negative Drude particle (SWM4-NDP), the TIP4P/2005-Reif nonpolarizable model, and an electronic continuum polarizable model (MDEC). Both SWM4-NDP and MDEC polarizable models were able to quantitatively reproduce the concentration dependence of the dielectric permittivity of NaCl aqueous solutions. On the contrary, the nonpolarizable TIP4P/2005 water model failed to quantitatively predict this concentration dependence. In contrast with the SWM4-NDP model, the MDEC model was unable to capture the concentration dependence of the structure and the dynamics of NaCl solutions. The SWM4-NDP model proved to be the most efficient polarizable model to reproduce quantitatively the concentration dependence of the dielectric permittivity, the dynamics, and the structure of NaCl solutions.
■
INTRODUCTION Accurate potential models for simulating the ion behavior in aqueous solutions are essential to reproduce the structure and dynamics of salt solutions in a wide range of concentrations. Models used in most molecular simulations have been developed on the basis of single ion properties.1 As a result, the applicability of these models to concentrated aqueous solutions remains challenging. Nevertheless, the development of force fields for aqueous electrolytes remains an active field of research.2−6 The decrease in the dielectric permittivity (ε) of water with increasing NaCl concentration was first observed in 1894.7 This phenomenon has been widely investigated both experimentally7−11 and theoretically. It has been interpreted in terms of dielectric saturation.9,12 Microscopically, the decrease in ε is due to the screening of water dipoles by ions involving a decrease in dipole−dipole correlations.13,14 Two competing phenomena are invoked in the liquid phase: (i) the destructuration of the hydrogen-bonding network leading to an increase in the dipolar fluctuations12,15,16 and (ii) the dielectric saturation due to the orientation of water dipoles around ions, which decreases dipolar fluctuations. It has been shown that the polarizability of ions is a key parameter to understand water structure around ions at liquid−vapor interfaces.17−21 Despite a huge amount of work devoted to © 2014 American Chemical Society
the impact of polarization effects, its influence on the dielectric properties has been much less studied.22,23 Accounting for ion polarizability explicitly can increase prohibitively the simulation time. The classical Drude oscillator24−26 and the MDEC6,27,28 models are attractive tools to account for ion polarizability without leading to a significant increase of the computational cost.6,28 Actually, these models conserve the simple functional form of the pairwise additive interactions while accounting for the electronic polarizability. As a result, these two models can be combined straightforwardly with various treatment methods for longrange electrostatic interactions. Within the Drude oscillator formalism, SWM4-DP,25 SWM4-NDP,26 and SWM629 water models have been proposed. The SWM6 water model29 adds two extra sites compared to the SWM4-DP and SWM4-NDP models25,26 in order to represent oxygen lone pairs. Since SWM6 requires a 2-fold increase in the computational effort,29 the SWM4-NDP model (which corresponds to a reparameterization of the original SWM4-DP model) was preferred in the present work. In the MDEC model, the polarization of Received: December 3, 2013 Revised: March 21, 2014 Published: March 24, 2014 3931
dx.doi.org/10.1021/jp4118419 | J. Phys. Chem. B 2014, 118, 3931−3940
The Journal of Physical Chemistry B
Article
charge model.37−40 In the core−shell model, the polarizable atom is represented by a massive core and a massless shell (Drude particle), connected by a harmonic spring. The core and shell carry different partial charges such that their sum is equal to the charge on the original atom. Thus, this model uses two sites: (i) a Drude particle attached to the core and carrying a charge qD and (ii) the core with a charge of q − qD such that q is the total charge of the atom. Both particles are electrostatically screened from each other and only interact via a harmonic restoring potential Ucore−shell = (1/2)kCS(rC − rS)2, where (rC − rS) is the distance between the core and shell and kCS is the force constant of the harmonic spring which is set to 4180 kJ mol−1 Å−2 for all Drude oscillators.4 An illustration of this model is provided in Figure 1, while the partial charges are
electrolytes is taken into account by scaling their charges by a factor of 1/(εe)1/2, where εe is the part of the relative permittivity which is due to electronic degrees of freedom (εe = 1.78 for water).6,27,28 The charges of the oxygen and hydrogen atoms of water molecules do not need to be rescaled because they have already been effectively scaled when developing common nonpolarizable water force fields by fitting them to experimental properties of water.6,28 The MDEC model can then be considered as a refinement of a nonpolarizable force field. Thus, it was then shown that the MDEC charge scaling procedure results in a better quantification of interactions between ions and nonpolarizable water modes in bulk phases.28 Recent simulations have also shown that the MDEC model improves significantly the interfacial properties of aqueous salt solutions compared with nonpolarizable models.30−32 Indeed, both ion and water structures in the vicinity of the air−water interface have been found in good agreement with experiments.30,31 In this work, we propose to examine the effect of the polarizability of water and ions on the dielectric permittivity, the structure, and the dynamics of ions and water in NaCl solutions of various concentrations. The concentration dependence of the structure and dynamics properties of aqueous NaCl solutions has already been studied by Lyubartsev and Laaksonen33 and Sala et al.23 Lyubartsev and Laaksonen used the flexible nonpolarizable SPC water model while the ions were treated as charged Lennard-Jones particles. They found a weak influence of the salt concentration on the structure, while a strong impact was found on the dynamics. More recently, Sala et al. investigated the dielectric and dynamics properties of aqueous NaCl solutions using the rigid SPC/E water model and RPOL polarizable model. They showed that both models did not reveal significant differences for Na−Cl and Cl−Cl radial distribution functions (RDFs), while significant deviations of RDF in Na−Na pairs were established. In this work, the dynamics of NaCl solutions will be explored by using the Drude and MDEC models. We will evaluate the performance of the classical Drude oscillator and MDEC models to reproduce the salt concentration dependence of the dielectric permittivity. The manuscript is outlined as follows: we first describe the different models, the computational procedures, and the calculation of the dielectric permittivity. Thereafter, we report and discuss results obtained for density, dielectric permittivity, structure, and dynamics of the various NaCl solutions. Finally, we conclude by providing a brief summary of our main results and addressing some interesting perspectives.
Figure 1. Illustration of ion and water models. For water, O represents the oxygen atom, H the hydrogen atom, and M a nonphysical particle. D is the Drude particle.
Table 2. Lennard-Jones Well Depth ε, Size σ, and Partial Charges q for Different Drude Oscillator Polarizable Water Modelsa SWM4-NDP polarizable model26 O H M D a
■
water ions
SWM4-NDP26 NDP4
TIP4P/200542 Reif5
TIP4P/200542 Reif5
1.716 36 0.55733 −1.114 66 −1.716 36
Table 3. Lennard-Jones Well Depth ε, Size σ, and Partial Charges q for Na+ and Cl− Ions Modeled with the Drude Oscillator4
+
MDEC
q (|e|)
106.15 0 0 0
kB is the Boltzmann constant.
Table 1. Water and NaCl Models and Their Respective Notations nonpolarizable (NP)
ε/kB (K)
given in Tables 2 and 3. Polarizable water was modeled from a Drude oscillator model labeled SWM4-NDP.3 SWM4-NDP is a reparameterization of the initial force field SWM4-DP developed by Roux and co-workers.25 The SWM4-NDP polarizable model is a five-site model which differs in its functional form by the sign of the charge of the Drude particle
SIMULATION METHODOLOGY The different models and their respective notations are listed in Table 1. Models and Interactions. Drude Oscillator Model. The polarizability of water molecules and ions was accounted for by the Drude oscillator model.4,26 This model accounts fairly well for induced polarization with limited computational cost in contrast to the induced dipole model34−36 and the fluctuation
Drude
σ (Å) 3.183 95 0 0 0
Na D Cl− D 3932
σ (Å)
ε/kB (K)
q (|e|)
2.9234 0 4.9630 0
15.84 0 36.19 0
+1.688 −0.688 +2.46 −3.46
dx.doi.org/10.1021/jp4118419 | J. Phys. Chem. B 2014, 118, 3931−3940
The Journal of Physical Chemistry B
Article
corresponds to the dipolar polarization. From the fluctuation−dissipation relation,46,47 εbulk can be expressed as
and that of the oxygen atom. In the SWM4-NDP model, the Drude particle D (bounded to the oxygen atom by a harmonic spring25) is negatively charged. The readers are directed to ref 41 for a comprehensive description of the polarizable Drude oscillator model. MDEC Polarizable Model. In the molecular dynamics in electronic continuum (MDEC) model, the polarization of aqueous electrolytes is taken into account by scaling their charges by a factor 1/(εe)1/2, where εe is the part of the relative permittivity due to the electronic degrees of freedom (εe = 1.78 for water).6,27,28 The charges of the oxygen and hydrogen atoms of the water model are not rescaled because they have already been effectively scaled in the development of this model by fitting to experimental properties of water.6,28 The LennardJones parameters of ions correspond to those of the nonpolarizable force field. This apparent simplicity makes this model attractive to account for polarization in a nonpolarizable force field. Nonpolarizable Model. In order to evaluate the ability of the above-mentioned polarizable models to predict the dielectric permittivity, we compared them with a nonpolarizable model for both water molecules and ions. We opted for the nonpolarizable rigid TIP4P/2005 model42 to model water molecules, while Na+ and Cl− ions were modeled from the force field developed by Reif and co-worker.5 TIP4P/2005 and Reif’s models were chosen because they are considered as very efficient to describe both water structure and ion hydration in the bulk phase. The force-field parameters for water and NaCl have been taken from ref 41. Computational Procedure. The equation of motions were solved using the Nosé−Hoover43 algorithm using 0.1 ps for the thermostat relaxation time with a time step of 1 fs. MD simulations were performed with the DL_POLY code,44 in which the core−shell dynamics was controlled by the adiabatic shell method.45 The initial simulation box was a rectangular box of dimensions Lx = Ly = Lz = 20 Å filled with water and ions. The total number of water molecules was fixed to 400, and the number of Na+ and Cl− ions was increased from 1 to 30 depending on the molality. Equilibration and acquisition steps were carried out under constant-NpT conditions with T = 300 K and p = 1 bar. Periodic boundary conditions were applied in the three directions. Standard deviations of the ensemble averages were calculating by breaking the production runs into block averages. The cutoff radius was set to 9 Å. The systems were equilibrated for 2 ns followed by an acquisition phase of 10 ns. Dielectric Permittivity Calculation. The dielectric permittivity of a medium is a macroscopic concept which is defined by the relationship between the total dipolar moment (M) and the electric field inside the medium (E).46 From a microscopic point of view, the polarization of the isotropic bulk medium (Mbulk) is related to the dipole density. The static dielectric constant in the bulk phase (εbulk) can be determined from the analysis of the total dipolar moment of water25,46,47 ε bulk = 1 +
4π ⟨M2⟩ 3ε0VkBT
ε bulk = 1 +
M=
∑ μi i
(2)
where (⟨M2⟩ − ⟨M⟩2) is the excess polarization or polarization density. For an isotropic medium, the ideal limit of zero net dipole moment implies ⟨M⟩ = 0 and then eq 1 is recovered. The dielectric permittivity of electrolytic solutions was calculated by considering Nwater
M=
∑ i
Nions − pairs
Ncore − shell
μi +
∑ i
μi +
∑ i
μi
(3)
As established in the MDEC method, the calculated dielectric permittivity is rescaled by (εe)1/2 = 1.335.28
■
RESULTS AND DISCUSSION NaCl Concentration Dependence of the Density. We report in Figure 2 the densities of the NaCl solutions obtained
Figure 2. Electrolyte solution densities predicted by the different models as a function of the salt concentration. The experimental densities are represented by the solid black line. The statistical fluctuations, corresponding to about 0.5% of the simulated density, are smaller than the size of the different symbols.
with the different models as a function of the salt concentration in the range [0, 4.2] mol kg−1. The comparison with the experimental densities leads to some relevant conclusions being drawn. First, in the range [0, 0.1] mol kg−1, we observe that the densities predicted by the nonpolarizable (NP) model are slightly overestimated with an average deviation of 0.5%, whereas the polarizable Drude oscillator model slightly underestimates the solution density with an average deviation from experiments of 1%. Within the same range of concentrations, the MDEC model better describes the solution density (deviation from experiments around 0.3%). For higher salt concentrations, the NP model deviates significantly from experiments, leading to a maximum deviation of 6%, whereas the Drude model leads to an underestimation of the salt densities with a maximum deviation of 4%. For the highest salt concentrations, the smallest deviation (
