arXiv:1001.2648v1 [physics.chem-ph] 15 Jan 2010 Models of electrolyte solutions from molecular descriptions: The example of NaCl solutions John Jairo Molina 1 , 2 , $^{3}$, ∗ Jean-Fran¸cois Dufrˆ eche 1 , 2 , $^{3}$, † Mathieu Salanne 1 , $^{2}$, Olivier Bernard1 , $^{2}$, Marie Jardat1 , $^{2}$, and Pierre Turq1 , 2 1 UPMC-Universit´e Paris 06, UMR 7195, PECSA, F-75005 Paris, France 2 CNRS, UMR 7195, PECSA, F-75005 Paris, France 3 Institut de Chimie S´eparative de Marcoule (ICSM), UMR 5257 CEA–CNRS–Universit´e Montpellier 2, Site de Marcoule, Bˆ atiment 426, BP 17171, 30207 Bagnols-sur-C`eze Cedex, France We present a method to derive implicit solvent models of electrolyte solutions from all-atom descriptions; providing analytical expressions of the thermodynamic and structural properties of the ions consistent with the underlying explicit solvent representation. Effective potentials between ions in solution are calculated to perform perturbation theory calculations, in order to derive the best possible description in terms of charged hard spheres. Applying this method to NaCl solutions yields excellent agreement with the all-atom model, provided ion association is taken into account. Since the pioneering works of Debye, H¨ uckel, and Onsager, electrolyte solutions have been commonly described by continuous solvent models, for which the McMillan-Mayer theory [1] provides a rigorous statistical-mechanical foundation. Within that level of description, simple phenomenological models such as the primitive model (PM), for which the ions are assimi- lated to charged hard spheres [2], can lead to explicit formulas for the thermodynamic and structural proper- ties (e.g., with the help of the mean spherical approxima- tion (MSA) [3] or the binding MSA (BIMSA) [4]). These models are the most practical to use [5], since they allow for a direct link between the experimental measurements and the microscopic parameters of the system. Never- theless, they ignore the molecular structure of the sol- vent. Consequently, they cannot properly account for the complex specific effects of the ions, which appear in numerous biological, chemical, and physical interfacial phenomena [6, 7], without further developments. An alternative procedure consists in carrying out molecular simulations, where both the solvent and solute are treated explicitly. After a rigorous averaging over the solvent configurations, a coarse-grained description of the ions, which still includes the effect of the solvent structure, can be obtained [8–11]. However, this set of methods is purely numeric; they do not provide any an- alytical expression for thermodynamic quantities. They are therefore restricted to simple geometries [12, 13] (bulk solutions or planar interfaces). The description of com- plex systems, such as porous or electrochemical materi- als, is still based on continuous solvent models [14]. In this letter we present a method aimed at bridging the gap between analytical and numerical approaches. It is based on the application of liquid perturbation theory (LPT) [15] to effective ion-ion potentials extracted from $^{∗}$Electronic address: john.molina@etu.upmc.fr $^{†}$Electronic address: jean-francois.dufreche@upmc.fr molecular dynamics (MD) results. Different approxima- tions of the PM are employed for the case of NaCl elec- trolyte solutions: a two component model (MSA2), that only takes free ions into account, and two different three component models (MSA3 and BIMSA3), which include a third species (the contact ion pair). As we proceed to show, LPT allows us to select the best simple model which accurately accounts for the thermodynamics and the physical-chemistry of the system. The first stage consists in calculating the McMillan- Mayer effective ion-ion interaction potentials V eff ij ( r ), by inverting the radial distribution functions (RDF) g$_{ij}$ ( r ) obtained by MD. The simulations were carried out on a box of 2000 water molecules and 48 NaCl pairs us- ing the same interaction potentials as in reference [16]. This setup corresponds to a concentration of 0 . $_{64 mol l}$− $^{1}$. NPT ensemble sampling at standard pressure and tem- perature was enforced, with a time step of 1 fs and a pressure bath coupling constant of 1 ps. An equilibration run of 0.25 ns was followed by a production run of 0.6 ns for five different initial configurations. The averages of the resulting RDF were then used for the potential inver- sion via the HNC closure [15]. These effective potentials are assumed to be concentration independent and will be used for simulations at all concentrations. Subtracting the long-range Coulombic potential V LR ij ( r ) (which depends on the dielectric constant of the solvent) from V eff ij ( r ), we obtain the short-range contri- bution V SR ij ( r ) to the effective potentials. These are given in Fig. 1 (species 1 and 2 refer to Na + and Cl − free ions, respectively). All the short-range potentials exhibit os- cillations corresponding to the solvent layering between the ions, but this effect is particularly important for the cation-anion interaction: a considerable potential barrier ( ≳ 2 k$_{B}$T ) separates the first two attractive wells. To serve as a reference, Monte Carlo (MC) simulations were performed with these effective potentials; a comparison between MD and MC RDF is also provided in Fig. 1. The excellent agreement between both sets of RDF validates the HNC inversion procedure [17], and allows us to com-