Activity Coefficients of Aqueous Sodium, Calcium ... - ACS Publications

Subscriber access provided by STEPHEN F AUSTIN STATE UNIV

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Activity Coefficients of Aqueous Sodium, Calcium, and Europium Nitrate Solutions from Osmotic Equilibrium MD Simulations Michael Bley, Magali Duvail, Philippe Guilbaud, and Jean-Francois Dufrêche J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b04950 • Publication Date (Web): 15 Jul 2018 Downloaded from http://pubs.acs.org on July 17, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Activity Coefficients of Aqueous Sodium, Calcium, and Europium Nitrate Solutions from Osmotic Equilibrium MD Simulations Michael Bley,∗,† Magali Duvail,† Philippe Guilbaud,‡ and Jean-Fran¸cois Dufrêche∗,† †ICSM, CEA, CNRS, ENSCM, Univ Montpellier, Marcoule, France ‡CEA, Nuclear Energy Division, Research Department on Mining and Fuel Recycling Processes (SPDS/LILA), BP 17171, F-30207 Bagnols sur Cèze, France E-mail: [email protected]; [email protected] Phone: + 33 4 66 33 92 06

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Osmotic and activity coefficients of three aqueous electrolyte solutions with cations of a similar ionic radius, but different charges are described by Molecular Dynamics with the help of the osmotic equilibrium method using polarizable force fields up to high concentration. Simulations of vapor-liquid interfaces of aqueous solutions of NaNO3 , Ca(NO3 )2 , and Eu(NO3 )3 at different concentrations and at 298.15 K provide timeaveraged number density profiles and consequently the quantity of solvent molecules in the vapor phase. These three cations of similar ionic radii exhibit an increasing amount of water in their first coordination sphere due to their increasing charge. The solvent activity is directly determined by the vapor phase density at different salt concentrations with respect to the vapor phase density of the pure solvent. The obtained densities of the liquid bulk and the osmotic and activity coefficients for the three different nitrate salts are in good agreement with the experimental results. Time-averaged concentration profiles and the interpretation of radial distribution functions are used to explain the role of coordination on the thermodynamic properties of aqueous electrolyte solutions.

Introduction There has been a long history of investigating the relation between the interfacial ion behavior and the consequences for physical chemistry. 1–4 Many phenomena in solution, such as salt solubility, osmotic and activity coefficients, surface tension, critical micellar concentrations, and many more, are determined by the properties of the dissolved ionic species. 5 In late 1880s, Hofmeister 6 and his co-worker Lewith 7 observed that different aqueous salt solutions lead to ion specific results in a systematic series of experiments. Years later the Hofmeister series for cations and anions were introduced. Cations of a high and anions of a low charge density increase the solubility of organic compounds (”salting in”) and decrease the surface tension, whereas cations of a low and anions of high charge density tend

2


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to ”salt out” and increase the surface tension due to structuration effects. 5 Robinson 8 and Stokes 9,10 subsequently developed the concept of ionic hydration explaining the behavior of activity coefficients obtained from isopiestic measurements above the dilute Debye-H¨ uckel (DH) regime. 11 The later developed mean spherical approximation (MSA) theory that allows determining ion hydration properties from activity coefficient measurements. 12–15 In 1942, Robinson et al. 8 measured the osmotic and activity coefficients of bivalent concentrated aqueous nitrate salt solutions and they showed that the activity coefficients of the salts increase with the cation’s charge density. Lowest activity coefficients were observed for aqueous solutions of the biggest and thus softest cation Ba(NO3 )2 , and highest activity coefficents were observed for aqueous Mg(NO3 )2 solutions. The corresponding ionic radii for a coordination number (CN) of 6 nearest neighbors are 1.35 ˚ A (Ba2+ ) and 0.72 ˚ A (Mg2+ ), respectively. 16 Herein, interfacial properties and activity coefficients of differently charged cations with similar radii are studied. The cation charge density increases from NaNO3 over Ca(NO3 )2 to Eu(NO3 )3 , and consequently an increasing activity coefficient in agreement with the Hofmeister series is expected.

Figure 1: Schematic representation of a simulation box presenting an aqueous sodium nitrate solution bulk liquid phase (bE,ini. = 3.0 mol kg−1 ) in equilibrium with the adjacent vapor phase.

These specific ion effects and the interfacial structuring of charges species has been intensively studied by Jungwirth, Tobias, and co-workers by means of molecular dynamics (MD) simulations in the past years. 17–21 In 2001, Jungwirth and Tobias 17 showed that the solute 3



depletion from the interface decreased along the Hofmeister series towards softer anions. An experimental confirmation of these observed ion density distributions at the vapor-liquid interface is hard to probe due to corrugation, disordering, and non-stationary behavior at a subnanometer scale. 5,18,22 Recently, our group proposed a new method for calculating activities in solution by MD simulations of vapor-liquid interfaces. 23 A schematic representation of the simulation boxes used for this approach is depicted in Figure 1. This method has been applied on concentrated aqueous Dy(NO3 )3 solutions 23 and of ethanol-water mixtures. 24 Solvent activities are directly obtained by the time-averaged amount of the solvent molecules present in the gas phase. This multiscale approach, the so-called osmotic equilibrium method, relates atomic and molecular density distributions with thermodynamics parameters of bulk phases. This is a major advantage of this approach compared to previously developed simulation techniques of osmotic coefficients such as the McMillan-Mayer theory, 25–28 the Kirkwood-Buff approach, 29,30 the grand-canonical monte carlo simulations, 31–33 or the osmotic membrane method. 34,35 In this article, the previously developed osmotic equilibrium method is applied on aqueous mono, bi, and trivalent nitrate salt solutions up to high concentrations. Na+ , Ca2+ , and Eu3+ are highly hydrated cations of similar ionic radii of 1.02 ˚ A (CN = 6), 1.06 ˚ A (CN = 7), and 1.07 ˚ A (CN = 8), respectively. 16 Their hydration properties have been part of countless experimental and theoretical studies, thus reliable force field parameters taking into account explicit polarization are available. 36–39 Herein, the trends for most likely first coordination spheres and time-averaged density and concentration profiles of all relevant compounds are directly related with the resulting water activities and solute activity coefficients at different salt concentrations. This is, to the authors best knowledge, the first time that specific ion effects at the interface and corresponding thermodynamic properties of the bulk have been simulated in very good agreement with experimental findings.

4


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Theoretical Methods Osmotic and Activity Coefficients The activity of solvent such as water aW (bE ) at a given salt molality bE is given by the ratio of the partial vapor pressures at a given concentration pW (bE ) and the vapor pressure of pure water p∗W . Considering the water in the vapor phase behaving as an ideal gas 23 allows expressing the solvent activity as a fraction of the number densities of the vapor phase in equilibrium with a mixture (nW (bE )) and the pure liquid (n∗W ), respectively.

aW (bE ) =

∗ ¯W (bE ) nW (bE ) NW (bE )VW N pW (bE ) = = = ∗ ¯∗ p∗W n∗W NW VW (bE ) N W

(1)

∗ ∗ are the are mean amounts of water in the gas phase, VW (bE ) and VW where NW (bE ) and NW

¯W (bE ) and N ¯ ∗ are mean amounts of solvent in gas corresponding gas phase volumes, and N W phase normalized to a reference volume of the vapor phase for a given concentration bE and pure water, respectively. Equation 1 is only valid if the liquid phase is in equilibrium with its gas phase. The osmotic coefficient of water at a given salt concentration φW (bE ) is

φW (bE ) = −

ln aW (bE) νMW bE

(2)

where ν is the amount of ions per formula unit of salt and MW is the molar mass of water. The osmotic coefficient of water is related to the activity coefficient of the salt γE (bE ) through the Gibbs-Duhem relation. 10 The resulting relation between the two thermodynamic properties is given by an integration of the molality b up to a given concentration bE Z

b=bE

ln γE (bE ) = φW (bE ) − 1 + b=0

5

φW (b) − 1 db b


(3)


Page 6 of 34

A good approximation for both ln γE (bE ) and φW (bE ) is the extended Debye-H¨ uckel (DH) equation of jth order. 40 The activity coefficient of the salt is obtained as √ j |z+ z− |Ab kbE X √ βi (bE )i ln γE (bE ) = − + 1 − Bb kbE i=1

(4)

where Ab and Bb are DH parameters on the molality scale, |z+ z− | is the absolute product of the ionic charges, and kbE is the equivalent expression for the ionic strength. The factor k depends on the valency of the electrolyte, thus kNaNO3 = 1, kCa(NO3 )2 = 3, and kEu(NO3 )3 = 6 are used. The sum in the last term adds fitting parameter contributions up to the jth order. The related formula for the osmotic coefficient of water gives √ p p |z+ z− |Ab kbE 1 √ φW (bE ) = 1 − (1 + Bb kbE ) − 2 ln 1 + Bb kbE − (Bb )3 kbE 1 + Bb kbE j X i + βi (bE )i i+1 i=1

(5)

In Equations 4 and 5, the first term is the extended Debye-H¨ uckel limiting law covering the dilute regime for activity and osmotic coefficients, respectively. The second term sums up additional fitting terms up to the jth order that are necessary to describe the activity coefficient at salt concentrations at higher salt concentrations. The parameter Ab is constant according to the equation p e3 2πNA ρW (T ) Ab = (4π0 ) 3/2 kB 3/2 r 3/2 T

3/2

(6)

where ρW (T ) is the mass density of the pure solvent at the temperature T , NA is the Avogardro number, e is the elemental charge, 0 is the vacuum permittivity of pure water, r is the relative permittivity of water and kB is the Boltzmann constant. The second parameter Bb is evaluated as

s Bb = αb Bb∗ = αb

2e2 NA ρW (T ) kB T r 0

where the parameter αb is fitted towards the experimental data.

6


(7)

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Precision of the Method In a previous article about this method, 23 the statistics of evaporation were described by a Poisson model. This approach allows estimating the error for the simulated mean amount of water molecules in the gas phase and for the resulting activities and osmotic coefficients. ¯W (bE ) follows a Assuming that the mean amount of water molecules in the gas phase N Poisson probability distribution leads to

¯W (bE ) = τW λW (bE ) N

(8)

where τW is the mean residence time of water in the gas phase and λW (bE ) is the concentrationdependent evaporation rate of water. τW is obtained through a Maxwell-Boltzmann distribution in the x-direction r τW = lx

πMW 2RT

(9)

where lx is the total length of the gas phase in x direction, MW is the molar mass of water, and R is the ideal gas constant. Combining Equations 1 and 8 yields

aW (bE ) =

¯W (bE ) λW (bE ) N = ∗ ¯ λ∗W NW

(10)

This shows that the activity of water depends on the fraction of the evaporation rates for the mixture λW (bE ) and for pure water λ∗W . The relative error of the mean amount of solvent molecules in the vapor phase δN¯W (bE ) is subsequently calculated

δN¯W (bE ) =

¯W (bE ) ∆N 1 p = ¯W (bE ) N λW (bE )tP

(11)

¯W (bE ) is the total error for the amount of water found in the gas phase and tP is where ∆N the production time of vapor-liquid equilibria in the NVT ensemble. The production time

7



Page 8 of 34

is given by tP = tTot − tEq

(12)

where tTot is the total simulation time and tEq is the time needed to equilibrate the system. The corresponding relative error of the activity of water δaW (bE ) is 1 ∆aW (bE ) =√ δaW (bE ) = aW (bE ) tP

1

1 p +p ∗ λW λW (bE )

!

r =

τW tP

1 1 p +p ¯∗ N N¯W (bE ) W

! (13)

where ∆aW (bE ) is the total error of the water activity. The latter is calculated via r ∆aW (bE ) =

τW tP

¯W (bE ) ¯W (bE ) N N + ¯∗ ¯ ∗ )3/2 N (N W W

p

! (14)

The resulting total error of the osmotic coefficient ∆φW (bE ) is d φW (bE ) δaW (bE ) 1 ∆aW (bE ) = = ∆φW (bE ) = d aW (bE ) νMW bE νMW bE

r

τW tP

1

1 p −p ¯∗ ¯W (bE ) N N W

! (15)

Simulation Details Molecular Dynamics Classical molecular dynamics simulations of aqueous nitrate salt solutions at different concentrations ranging from 0.5 to 11.0 mol kg−1 have been carried out with SANDER14, a module of AMBER14, 41 using explicit polarization in the NPT and NVT ensembles. Periodic boundary conditions were applied to the simulation box in all directions. The equations of motion were numerically integrated using a 1.0 fs time step and the long-range interactions have been calculated using the particle-mesh Ewald method. 42 The systems were equilibrated at 298.15 K and 1 bar (0.1 MPa) for at least 5 ns, and production runs were afterwards conducted for up to 15 ns. Table 1 lists the force field parameters used in the present work. All atomic coordinates were written to the trajectory file every ps. Here, the 8


Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


van der Waals energy is described by a 12-6 Lennard-Jones potential. The water molecules were described by the rigid POL3 model, 43,44 which takes into account the polarization. The Lennard-Jones parameter for sodium are described by Aqvist 36 , whereas the atomic polarizability developed by Dang 37 is used. This parametrisation provides a Na+ coordination in good agreement with experimental 45,46 and simulated data. 36,37,47 For Ca2+ , both Lennard-Jones parameter and atomic polarizability are taken from Dang et al. 38 . The obtained simulated coordination for the calcium cation is in agreement with previous experiments 45,48,49 and simulation 38,50 . The Eu3+ cation is described by the model previously developped by 39 set to reproduce the experimental hydration properties. This polarizable force field for describing the Eu3+ cation in water provides structural and dynamical properties in good agreement with experimental 51,52 and theoretical 53–55 results. The NO− 3 anion is described by the polarizable model defined in a previous study by Duvail et al. 56 , providing structural properties of NO− 3 in aqueous solutions in reasonable agreement with published experimental 57 and calculated data. 19 All three aqueous salt solutions provide simulated densities similar to experiment. 58–60 Table 1: Parameters Used for the MD Simulations

a

atom ii a σii b qi c αd Na+ 0.012 3.328 +1.000 0.250 2+ Ca 0.418 2.913 +2.000 0.470 Eu3+ 0.322 3.271 +3.000 0.823 NNO−3 0.711 3.25 +0.896 0.530 ONO−3 0.879 3.067 -0.632 0.434 HH 2 O +0.365 0.170 OH2 O 0.655 3.204 -0.730 0.528 Energies in kJ mol−1 . b Distances in ˚ A. c Atomic partial charges in e. polarizabilities in ˚ A3 .

d

Atomic

All simulation boxes have been created with the PACKMOL package. 61 The composition of the boxes and other characeristics are tabulated in the Supplementary Information (Tables S1 to S3). The equilibration stage was conducted for 5.0 to 10.0 ns in the NPT canonical

9



ensemble at 298.15 K and 1.0 bar using a τPres. of 0.1 ps, where τPres. is the relaxation time of the Berendsen pressure coupling. 62 For all simulation boxes for NaNO3 and Eu(NO3 )3 , three different equilibrated coordinate files of an identical composition have been used for each salt concentration. In the case of aqueous Ca(NO3 )2 solutions, more different concentrations but only a single trajectory for each concentration have been simulated. The simulation of vaporliquid equilibria of a pure compound and the corresponding amount of solvent molecules in the gas phase has already been discussed in our previous articles. 23,24 Subsequently, simulation boxes have been extended to a finite size of 495.0 ˚ A in x direction to obtain vapor-liquid interfaces leading to a centered equilibrated bulk liquid phase and two adjacent vacuum areas with a width of around 225 ˚ A each. It should be noted that this vacuum space will hereafter be referred as the vapor or gas phase. A schematic representation of the simulation boxes used for the production runs is given in Figure 1. These geometries are used to produce vapor-liquid equilibria in the NVT ensemble for 15.0 ns at 298.15 K. The cut-off radius for all production runs is 12.0 ˚ A.

Data Analysis Data analysis has been performed on the 93 production runs for aqueous salt solutions (30 for each NaNO3 and Ca(NO3 )2 , and 33 for Eu(NO3 )3 ). Equilibrium between the liquid and the vapor phase is obtained after around 8 ns for each production run for each salt. Therefore, only the results of the remaining 7.0 ns are taken into account for calculating the time-averaged amount of evaporated water hNW,i (bE )i at salt concentration bE during the ith run. The determination of the amount of water present in the gas phase for a single step at a given time NW,i (t, bE ) requires the locations of two interfacial yz-planes perpendicular to the x axis. The latter are therefore located at 200 and 295 ˚ A for all trajectories, providing two gas phases of identical volume. All oxygen atoms of water OW within the range from 0 to 200.0 ˚ A and from 295.0 to 495.0 ˚ A are counted for all frame of an individual trajectory. Time-averaged density profiles of the production runs for water and the dissolved ions are 10


Page 10 of 34

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


obtained by calculating the arithmetic mean of a given compound for every 20th frame of the trajectory file within an interval of 1.0 ˚ A in x direction. The time-averaged amount of water in the gas phase hNW,i (bE )i for the i th run at the concentration bE is obtained by ∗ VW hNW,i (bE )i = tP VW,i (bE )

Z

tTot

NW,i (t, bE )dt

(16)

tEq

where NW,i (t, bE ) is the amount of water molecules present in the gas phase at a simulation ∗ time t during the i th simulation run for a given salt concentrationm, and VW and VW,i (bE )

are gas phase volumes of a simulation box for pure water and of the simulation box of the ∗ /VW (bE ) in Equation 16 originates from i th run for a given salt molality. The prefactor VW

the definition of the number densities (Equation 1) and it is usually close to 1. It transforms the expression of hNW,i(bE ) i so that the value corresponds to the same volume of the gas phase as in the case of pure water. The equilibration time as defined in Equation 12 is set to 8 ns, hence for a global simulation time tTot of 15 ns, the analyzed production time tP is 7 ns. A production time of 7 ns corresponds to 7000 analyzed frames for each obtained trajectory. The amount of molecules in the gas phase is obtained for slightly different gas phase volumes for each run i. The volumes of the gas phase differ due to slight different ¯W (bE ) is initial atomic coordinate files. The mean amount of water in the vapor phase N obtained by calculating the arithmetic mean

¯W (bE ) = N

Pi=j

i=1 hNW,i (bE )i

j

(17)

The mean of j slightly different trajectories for a concentration bE allows then calculating the raw water activity with Equation 1 and the evaporation rate λW (bE ) via Equation 8.

11



Results and Discussion Liquid Phase The liquid densities after equilibration in the NPT ensemble is in agreement with the experimental findings. For aqueous sodium nitrate solutions around saturation (bNaNO3 = 11.0 mol kg−1 ), simulation provides a mass density of 1.378 ± 0.007 compared to 1.398 g cm−3 from experiment. 59 Experimental densities of Ca(NO3 )2 solutions at high concentration have been reported by Ewing and Mikovsky 58 up to 20 mol kg−1 . At an initial calcium nitrate concentration of 8.0 mol kg−1 , simulation provides 1.530 ± 0.008 compared to 1.568 g cm−3 from literature. In the case of Eu(NO3 )3 , Rard 60 measured density only in the dilute regime until 1.10 mol kg−1 . Comparing simulated and experimental densities for bEu(NO3 )3 = 1.0 mol kg−1 provides ρSim. = 1.242 ± 0.006 g cm−3 and ρExp. = 1.257 g cm−3 . Using a polynomial extrapolation of the experimental data leads to ρSim. = 2.031 g cm−3 and ρExp. = 1.975 g cm−3 at 6.5 mol kg−1 . Simulated and experimental densities for all concentrations and other characteristics of the simulation are part of the Supplementary Information (Tables S1 to S3). Table 2 provides selected structural properties for the three aqueous electrolyte solutions at different concentrations. It should be noted that the three cations exhibit similar ionic radii. For Na+ , the most probable coordination number in the first coordination sphere is 6 and the ionic radius is 1.02 ˚ A. 16 Here, the usual coordination number for the calcium cation is 7 and for europium 8 nearest neighbors are found in our system. The ionic radii of Ca2+ and Eu3+ are 1.06 and 1.066 ˚ A, respectively. 16 The increase of the coordination number goes hand in hand with the increase of the charge density, and Eu3+ is a way harder ion than Na+ . Thus, the increasing ionic charge should predominantly determine the coordination of these cations. Globally, the amount of water in the ion’s first hydration sphere CNM−OW decreases for an increasing salt concentration, whereas consequentially an increasing amount nitrate anion close to the cation is found. The cation-water distance M-OW remains almost

12


Page 12 of 34

Page 13 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


constant for all concentrations, but slightly increases along the cation series, since Eu3+ is slightly bigger than Na+ . Table 2: Hydration Properties of M(NO3 )x in Aqueous Solution at 298.15 K M(NO3 )x a bE,ini b NaNO3 0.50 3.00 6.00 11.00 Ca(NO3 )2 0.50 3.00 6.00 8.00 Eu(NO3 )3 0.50 3.00 6.00

rM−OW c 2.40 2.39 2.39 2.39 2.43 2.43 2.43 2.42 2.45 2.44 2.44

CNM−OW d 5.48 4.85 3.94 2.87 6.87 6.38 5.36 4.63 8.04 7.96 5.83

r1M−NN e 3.08 3.03 3.02 3.02 3.01 2.99 2.99 2.99 3.02 2.99

r2M−NN f 3.41 3.44 3.42 3.41 3.61 3.59 3.59 3.59 3.69 3.68 3.68

CN1 /CN2g 40.9 36.5 37.6 38.8 10.8 12.7 11.8 10.2 0.7 0.6

h CN1+2 M−NN 0.13 0.86 1.86 3.03 0.07 0.62 1.69 2.45 0.01 0.21 2.04

rNN −NN i 5.75 5.13 4.81 4.70 4.51 4.16 4.03 4.07 4.39 4.04 4.03

rM−M j 3.62 3.66 3.73 3.76 8.45 6.92 6.57 6.43 11.07 8.42 6.70

rOW −OW k 2.74 2.74 2.74 2.74 2.74 2.75 2.78 2.83 2.74 2.85 2.90

a

Salt M(NO3 )x with the cation Mx+ dissolved in water. b Initial M(NO3 )x concentration on a molality scale in mol kg−1 . c First maximum peak of the M-OW RDF in ˚ A. d Number of water molecules in the Mx+ first coordination sphere. e First maximum peak of the M − NN RDF in ˚ A. f Second maximum peak of the M − NN RDF in ˚ A. g Ratio of NO− 3 in the first h and second coordination sphere in %. Number of water molecules in the Mx+ first coordination sphere in mono- and bidentate coordination. i First maximum peak of the NN − NN RDF in ˚ A. j First maximum peak of the M-M RDF in ˚ A. k First maximum peak of the OW − OW RDF in ˚ A. However, the way the nitrate coordinates the cations drastically changes with the charge density of the cation. The lower the charge is and thus lower the amount of neighbors in the first coordination sphere are, the higher the percentage of nitrates is undergoing bidentate coordination. For aqueous NaNO3 solutions, around 40 % of the nitrate in the Na+ first coordination sphere is in bidentate coordination closely to the cation. Ca2+ strongly prefers a mondentate coordination (CN1 /CN2 ) around 10-12 %, whereas almost no bidentate coordination is observed in aqueous Eu(NO3 )3 solutions. The distance between the cation 1 and the nitrogen of the nitrate anion rM−N determined from the first maximum peak of the N

radial distribution function (RDF) is usually found around 3 ˚ Afor the bidentate coordination. 2 On the other hand, the distance for the monodentate coordination rM−N slightly increases N

with the charge of the cation. Figure 2 represents the radial distribution functions and the 13



corresponding coordination numbers for all three aqueous salt solutions at bE = 3.0 mol kg−1 . For Na+ and Ca2+ , two peaks at distances of 3.0 and around 3.5 ˚ A are observed indicating the occurrence of the two coordination states mono- and bidentate. At this salt concentration, the highest amount of nitrate in the first coordination sphere is found for Na+ and the lowest for Eu+ . However, in the second coordination sphere at a distance of around 5.0 ˚ A the inverse is observed. gHrM-NN L NNN HrM-NN L

7 6

NNN HrM-NN L

5

gHrM-NN L

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

NaNO3

4

CaHNO3 L2

3 2

EuHNO3 L3

1 0 0

2

4

6

8

10

12

rM-NN

Figure 2: M-NN radial distribution functions (solid lines) and coordination numbers (dashed lines) obtained for binary aqueous solutions of NaNO3 (orange), Ca(NO3 )2 (green), and Eu(NO3 )3 (gray) solutions at bE,ini = 3.0 mol kg−1 .

The first coordination sphere of all studied cations at high salt concentration (bE = 6.0 mol kg−1 ) is depicted in Figure 3. In the case of sodium nitrate (Figure 3a) it should be noted that the cation-cation distance rM−M is very small compared to the aqueous solutions of calcium and europium nitrate. The Na+ ions are sharing the nearest nitrate anion with their neighbor cations. This results in complex interactions between the cationic and anionic species. NO− 3 coordinates in a monodentate way one cation and at the same time another Na+ is approached by the remaining two oxygen atoms of the nitrate. This is an extremely dynamic process with coordination changes below the picosecond scale. For the bi- and 14


Page 15 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trivalent nitrate salts, the distance between two anions is larger even at high salt molalities. Bidentate coordination is observed for Ca2+ (Figure 3b), but only monodentate coordination occurs for the trivalent europium salt in Figure 3c. The tendency to undergo bidentate coordination is thus less likely in the case of steric constraints due to a higher amount of solvent molecules in the first Eu3+ coordination sphere compared to Ca2+ . This also 2 explains the increasing rM−N distance for the monodentate coordination. Experimental N

structural studies 63 and unconstrained simulation 39,56 of aqueous lanthanoid nitrate solutions showed that the monodentate nitrate coordination is preferred. For concentrated aqueous Ca(NO3 )2 solutions, X-ray diffraction and Raman spectroscopy showed that the monodentate coordination is preferred. 48

Figure 3: Snapshots representing characteristic cation first coordination spheres at 6.0 mol kg−1 for aqueous solution of (a) sodium (orange), (b) calcium (green), and (c) europium nitrate (gray). For sodium, Na+ pairs are observed due to the small rNa−Na distances.

The distance between two nitrate anions decreases for the different cations with the salt concentration and the ionic charges since less water molecules are available for hydration (Table 15



2). Furthermore, the increasing ionic charge results in higher intercationic distances rM−M and for Ca2+ and Eu3+ a higher salt content leads to a decreasing cation-cation distance. The Na+ -N+ distance increases slightly with bE . The ion density also influences the mean distance between the water molecules. The more ions are present in the liquid phase, the less dense is the water network, since the water molecules are required for a proper hydration of the dissolved ionic species. Additional data concerning the structural information of the aqueous electrolyte solutions is available in Supplementary Information (Tables S4 to S9 and Figures S1 to S6).

Vapor-Liquid Interface Since MD simulations for calculating the osmotic equilibrium approach require the presence of a vapor-liquid interface, it is necessary to understand the influence of the latter on the spatial distribution of all species involved. Therefore, time-averaged concentrations profiles have been derived from the MD trajectories for all the concentrations and all the species simulated. This is required since the appearance of a vapor-liquid interface changes the local concentration. Charged species move away from the interface towards the bulk liquid phase and thus the interfacial layer is formed by pure water. The time-averaged concentration profiles for the three different aqueous solutions at a medium concentration of 3.0 mol kg−1 is shown in Figure 4. Water accumulates at the interface and this water rich-zone is followed by a bulk liquid phase containing the dissolved electrolytes. Thus, ions are globally expelled from the interface. Repulsive forces, such as image charges as described in the model of Onsager and Samaras dominate over attractive effects. 64 The distribution for NaNO3 (Figure 4a) and Ca(NO3 )2 (Figure 4b) solutions is almost homogeneous, whereas the Eu(NO3 )3 (Figure 4c) exhibits an inhomogeneous distribution of the charges species due to the higher ion density. This has been previously observed for aqueous solutions of Dy(NO3 )3 23 and this behavior originated in the reduced amount of solvent molecules per ion pair and allowing thus the sharing of the water molecules between different ions. The time-averaged concen16


Page 16 of 34

Page 17 of 34

tration profiles for Na+ and NNO3 − are very similar, but not identical (Figure 4a). Additional Comparative time-averaged molar concentration profiles are provided in Supplementary Information (Figures S7 to S12). 60

Na+

HaL

50

cΑ mol dm-3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ca2+ HbL OH2 O NNO3-

Eu3+ HcL

60 50

40

40

30

30

20

20

10

10

0

dx Å

230 250 270

dx Å

230 250 270

dx Å

0

230 250 270

Figure 4: Time-averaged molar concentrations cα profiles for a species α along the x -axis for aqueous (a) sodium, (b) calcium, and (c) europium nitrate solutions at an initial salt molality bE,ini. = 3.0 mol kg−1 .

Hence, the time-averaged concentration profiles allow determining the distance between the interfacial water layer and the interfacial ions layer for all the ion species. The position of the interfacial layer of a species α is defined as the position dx where cα = 0.3 mol dm−3 . For the three different cations, similiar trends for all concentrations are observed for the distance between the interfacial water layer and the first interfacial ion layer. At low concentrations (bE = 0.5 mol kg−1 ), first ions are usually found at a distance of around 5.0 ˚ A. This distance decreases to 2.0 - 2.5 ˚ A for a molality of 6.0 mol kg−1 for the cations. The increasingly charged cations tend to stay in the bulk liquid keeping their preferred first sphere coordination number, while the softer nitrate anions prefer moving towards the vapor liquid interface. The presence of nitrate close to the interface also depends on the amount of dissolved nitrate ions. The longest distance of 4.5 ˚ A is observed for a low NaNO3 concentration (0.5 mol kg−1 ), whereas the shortest distance of 0.5 ˚ A is consequently observed for the highest Eu(NO3 )3 concentration of 6.0 mol kg−1 . This means that the charge of the cation and thus the equivalent stoichiometric corresponding amount of nitrate anions 17



determine the interfacial composition. This is due to steric constraints and the necessity of hydrating the cations even with interfacial water molecules. A Table listing all these interfacial distances is found in Supplementary Information (Table S10).

Figure 5: Parts of snapshots of the simulation boxes presenting the liquid phase (blue background) in contact with the two vapour phases(green background) for aqueous (a) sodium, (b) calcium, and (c) europium nitrate solutions at an initial salt molality bE,ini. = 3.0 mol kg−1 .

The depletion of charged species at the vapor-liquid interface in x direction is also observed on the corresponding snapshots of Figure 5 for bE = 3.0 mol kg−1 . The accumulation of charges species also changes the concentration of the bulk liquid phase and this needs to be taken into account before determining thermodynamic properties such as osmotic and activity coefficients. Fortunately, the effective bulk liquid densities can be calculated from the number densities in the middle of the time-averaged density profiles. Consequently, this corrected molality bE,liq. is used for further calculations of thermodynamics properties. Globally, the difference between initial composition bE,ini. and the corrected effective molality bE,liq. is small for low salt concentrations. However, for highly concentrated aqueous NaNO3 solutions an important offset to the initial value is found (e.g bE,ini. = 11.0 mol kg−1 and bE,liq. = 24.9 mol kg−1 ) due to the strong Na+ -Na+ interactions (Table 2). Initial and effective bulk liquid concentrations on a molality scale, and also a comparison of the snapshots for different concentrations are part of the Supplementary Information (Tables S1 to S3 and Figure S13 and S14).

18


Page 18 of 34

Page 19 of 34

Osmotic and Activity Coefficients Time-averaged number density profiles of the oxygen atom of water OW and the determination of the quantity of water present in the gas phase are required for calculating the thermodynamic properties from MD simulations. The mean amount of water molecules in ¯ ∗ = 2.29 as liquid phase has previously been determined 23 the gas phase for pure water N W ∗ for a reference gas phase volume VW = 507 nm−3 . Therefore, Equation 16 allows calculating

the normalized average number of water hNW,i (bE )i for a single production run, and thus Equation 17 provides the mean water amount for a given salt concentration. While using Equations 1 and 14, the water activity aW (bE ) and the corresponding total error δaW (bE ) is calculated for all investigated systems. The resulting water activities are plotted against the corrected bulk liquid phase concentrations bE,liq. in Figure 6. 1.2

NaNO3

1.0

aW HbE,liq. L

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.8

CaHNO3 L2

0.6 0.4 0.2

EuHNO3 L3

MD Exp.

0.0 0

2

bE,liq. mol kg 4

6

8

10

-1

Figure 6: Comparison of simulated and experimental water activities aW for aqueous NaNO3 (orange), Ca(NO3 )2 (green), and Eu(NO3 )3 solutions. Dashed lines represent experimental data. 10,65,66

Figure 6 compares simulated data points up to bE,liq. = 10.0 mol kg−1 with the corresponding experimental data (dashed lines). Raw simulations results are in good qualitative and quan19



Page 20 of 34

titative agreement with experimental findings for NaNO3 , 65 Ca(NO3 )3 , 10 and Eu(NO3 )3 66 up to high concentration. Sodium nitrate exhibits the highest water activities for both, simulation and experiment, and decreases only slowly for an increasing salt content. The different curves for the solvent activities of mono-, bi-, and trivalent salt solutions can be explained by Raoult’s law considering an ideal behavior for the system. In the dilute regime, aW (bE ) = xW (bE ) yields aW (bE ) =

1 1 + νMW bE

(18)

and corresponding first order Taylor series equivalent

aW (bE ) = 1 − νMW bE

(19)

where only the stoichiometric amount of ions ν varies for different salt valencies. Equation 19 corresponds to a linear behavior for solvent activities and the higher ν, the steeper the slope and the lower the activities for an increasing concentration bE . The simulated water activities and their total errors for the aqueous NaNO3 solutions are usually between aW = 0.8 and 1.0 along the investigated concentration range. For Ca(NO3 )2 and even more for Eu(NO3 )3 , the decreasing solvent activities are found in very good agreement with experimental data. Calcium nitrate exhibits more data points (16) compared to NaNO3 (10) and Eu(NO3 )3 (11), but only a single production run has been conducted in the case of Ca(NO3 )2 . On the contrary, all simulated water activities for aqueous sodium and europium nitrate solutions are derived from three different production runs and this shows that both approaches lead to reasonable results. Additional data for raw simulated results and the corresponding errors are provided in Supplementary Information (Figure S15 and Tables S11 to S16).

20


Page 21 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Fitting Parameters using Extended Deybe-H¨ uckel Equations. Salt αb a βMD b χ2MD c βExp. d χ2Exp. e NaNO3 2.85 -0.0699 0.351 -0.0135 0.999 Ca(NO3 )2 4.0 0.1777 0.338 0.1639 0.999 Eu(NO3 )3 4.3 0.4145 0.899 0.4038 0.999 a Fitted Debye-H¨ uckel parameter αb in ˚ A. b First-order fitting parameter in kg mol−1 for molecular dynamics results. c Value of the goodness of the fits for molecular dynamics results. d First-order fitting parameter in kg mol−1 for experimental results. e Value of the goodness of the fits for experimental results. These water activity are subsequently used for the calculation of osmotic coefficients φW (bE ) via Equation 2. The resulting osmotic coefficients are then fitted with extended DebyeH¨ uckel equations of the 1st order (Equation 5 and j = 1). The Debye-H¨ uckel parameter Ab and Bb have been calculated for 298.15 K and a dielectric constant of water r = 78.5. In Table 3 the fitting parameters for both experimental and simulated data are listed for the corresponding Extended Deybe-H¨ uckel Equations. All simulated osmotic coefficients φW,MD for bulk liquid concentrations higher than 1.0 mol kg−1 have been used for the fits. The Debye-H¨ uckel parameter αDH determines the shape of the fitted curve for low concentrations and usually determines the local minimum of the curves. NaNO3 has a low influence on the water activity, the osmotic coefficient decreases over the whole investigated molality range. This explains the difference for αDH compared to the rather similiar coefficients for aqueous calcium or europium nitrate solutions. Since the same value for αDH can be used for fitting the experiments and the simulations, the quality of the simulation results are directly compared with experimentally observed behavior via the remaining fitting parameter β. Table 3 shows a good agreement for simulated βMD and experimental βExp. for calcium and europium nitrate solutions, but a discrepancy of a factor 2 between βMD and βExp. for NaNO3 in water. Nevertheless, the qualitative trends are the same for all three species. The estimates of the goodness of the fit (χ2 ) are rather low since the last term in Equation 4 shows linear behavior in the 1st order. This means that the fitted curve represents an arithmetic mean of the simulated osmotic coefficient data. Fitting the simulated data with higher-order fitting

21



terms would result in some zigzag curves not properly representing the qualitative trends of simulation.

EuHNO3 L3

2.2 2.0

MD Exp.

1.8

ΦW HbE,liq. L

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 34

CaHNO3 L2

1.6 1.4 1.2 1.0

NaNO3

0.8 0.6 0

bE,liq. mol kg-1

2

4

6

8

Figure 7: Comparison of experimental 10,65,66 (dashed lines) and fitted simulated (lines) water osmotic coefficients φW for aqueous NaNO3 (orange), Ca(NO3 )2 (green), and Eu(NO3 )3 solutions.

The resulting fitted curves for the osmotic coefficients of all three aqueous salt solutions are depicted in Figure 7. The fitted curves for europium and calcium nitrate solutions are in good agreement with experimental findings, 10,66 whereas simulations underestimate the osmotic coefficients for NaNO3 . Since the osmotic coefficient of water φW and the activity coefficient of the salt γE are related through the Gibbs-Duhem relation, Equation 4 allows calculating activtiy coefficients while using the same fitting parameters αDH and β from Table 3.

22


Page 23 of 34

EuHNO3 L3

1.6 1.4

MD Exp.

1.2

ΓE HbE,liq. L

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CaHNO3 L2

1.0 0.8 0.6

NaNO3

0.4 0.2 0

1

bE,liq. mol kg-1

2

4

6

8

Figure 8: Comparison of simulated (lines) and experimental 10,65,66 (dashed lines) salt activity coefficients γE calculated with Equation 4 for aqueous NaNO3 (orange), Ca(NO3 )2 (green), and Eu(NO3 )3 solutions.

A comparison of fitted experimental and simulated results based on Equation 4 and the parameters given in Table 3 is shown in Figure 8. Only for aqueous NaNO3 solutions, simulations and experiments lead to slightly different results. This different behavior of NaNO3 solutions compared to Eu(NO3 )3 and Ca(NO3 )2 solutions has two possible reasons. The first explanation for this offset is the quality of the force field parameters used for Na+ . Indeed, the chosen current parametrisation for Na+ provides good results for the first coordination sphere, the liquid phase densities, the time-averaged density profiles, and the qualitative trends for the thermodynamic properties in agreement with the trends of the Hofmeister series. 5 However, it appears that the potential for Na+ used here is too attractive. The second aspect in the case of sodium nitrate is the fact that the experimental water activities lie between 0.7 and 1.0 and the change of aW with the salt concentration is rather small. Therefore, for a broad concentration range simulated mean amounts of water ¯W (bE ) are quite close to the reference value N ¯ ∗ . Previous applications of the osmotic N W equilibrium approach on aqueous Dy(NO3 )3 solutions at low ion density (bE < 2.0 mol kg−1 ) 23



showed that the precision is low and error of these activities are higher compared to the results in the concentrated regime. 23 It should be noted that the precision increases for the bi- and trivalent nitrate salts, since the amount of dissolved ions at the same salt molality is increased. The precision of this simulation method is remarkably high for salt molalities above 2.0 mol kg−1 up to saturation. It was recognized that in order to enhance the precision of the method it is possible either to perform longer simulations for every concentration or to keep the same simulation time but for more concentrations. Both procedures provide results in very good agreement together with the experiments (Figure 6). Raw simulated osmotic coefficient of water and their corresponding total errors, the extended Debye-H¨ uckel fit of the simulated data, experimental data points and the results for the activity coefficients are depicted in Supplementary Information (Figures S16 to S21).

Conclusion Molecular dynamics simulations of vapor-liquid equilibria of aqueous mono- (NaNO3 ), bi(Ca(NO3 )2 , and trivalent (Eu(NO3 )3 ) nitrate salt solutions at different salt molalities have been performed using explicit polarization. Structural properties in the first coordination sphere of the three dissolved ions up to high concentration are consistent and are in good agreement with previous experiments and simulations. The obtained simulated densities for the bulk liquid phases are similar compared to experimental measurements. Furthermore, the time-averaged density or concentration profiles taking into account the vapor-liquid interface show that a homogeneous distribution of all species in the bulk liquid phase is granted for all the systems studied. For all three salts, an ion-depleted interfacial layer of water molecules is found which is followed by a the bulk liquid phase. The amount of dissolved ions determines the distance of nitrate anions to the vapor-liquid interface, thus implying that NO− 3 closer to the interface are found for aqueous solutions of calcium and europium nitrate. Furthermore, the higher concentration of charged species leads to a pattern of the

24


Page 24 of 34

Page 25 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molar concentration profiles with respect to the interface mostly for Eu(NO3 )3 and up to a certain extend also for Ca(NO3 )2 at salt molalities above 4.0 mol kg−1 . This behavior originates from the need of proper ion hydration shells and the possibility of sharing neighbored water molecules. The cation coordination and their interfacial density profiles are determined by ion specific effects which is consistent with the literature. 17,18 The polarizable x+ force field parameters for water and the ionic species NO− allow calculating the 3 and M

water activities, the osmotic coefficients of water, and the activity coefficients of the salts in good qualitative and quantitative agreement with experimental data. The experimental and fitted simulated curves for Ca(NO3 )2 and Eu(NO3 )3 match each other, whereas NaNO3 exhibits a slight offset between the two data sets. It is shown that the osmotic equilibrium approach is selective for different dissolved cations and provides salt activity coefficents via the extended Debye-H¨ uckel equation and the Gibbs-Duhem relation in good agreement with experimental studies. It should be noted that this approach for accessing salt activity coefficients in aqueous solutions relies on the directly obtained water activities aW form the MD simulation of the vapor-liquid equilibrium. The results from the extended DH fit of the corresponding osmotic coefficients are an interpretation of the simulated results. Simulating a higher amounts of different salt molalities provides more data points and thus a higher confidence in the outcome of the fit. Since the calculated thermodynamic properties of the bulk liquid and interfacial ion distribution are in agreement with the concept of the Hofmeister series and the literature, 5,18 the simulated trends for salt activities might be another proof for the correctness of their observations regarding ion density profiles. The osmotic equilibrium method provides a powerful tool for accessing the activity coefficients via molecular simulations using a multiscale approach. Using consistent force field parameter for single atoms or molecules allows predicting their coordination and interfacial properties by the means of molecular dynamics. The resulting time-averaged density profiles of vapor-liquid interfaces for aqueous salt solutions allow determining the composition and the properties of an adjacent vapor phase. Therefore, the amount of the water molecules

25



Page 26 of 34

calculated from MD simulations directly leads to the thermodynamic properties in solution. It has been shown that the osmotic equilibrium method is a versatile approach for predicting physico-chemical properties for a huge bandwidth of different systems such as complex organic phases used in separation chemistry or solutions difficult to handle, which are part of the treatment of spent nuclear fuels. In a future work, a more detailed insight in ion-paring kinetics and the comparison with experiment could be assessed, similarly to what was previously perfomed for lanthanide cations binary solutions using the McMillan-Mayer approach by Molina et al. 67 .

Supporting Information Available Simulation details such as compositions of the boxes (Tables S1 to S3), a detailed description of the coordination in the aqueous phase (Tables S4 to S9 and Figures S1 to S6), a comparison of time-averaged concentration profiles of the vapor-liquid interface (Figures S7 to S12), the distance of all species to the vapor-liquid interface (Table S10), the difference between initial and effective bulk liquid phase concentration (Figure S13), a comparison of snapshots of the simulation boxes (Figure S14), all raw water activities (Figure S15), and raw and fitted data for the osmotic and activity coefficients compared with experimental data (Figures S16 to S21 and Tables S11 to S16) are part of the Supplementary Information file. This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement This work was made possible thanks to the high performance computing facilities of TGCC/CCRT and the computing center of CEA Marcoule.

26


Page 27 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Boström, M.; Williams, D. R. M.; Ninham, B. W. Surface Tension of Electrolytes: Specific Ion Effects Explained by Dispersion Forces. Langmuir 2001, 17, 4475–4478. (2) Kunz, W.; Lo Nostro, P.; Ninham, B. W. The Present State of Affairs with Hofmeister Effects. Curr. Opin. Colloid Interface Sci. 2004, 9, 1–18. (3) Boström, M.; Kunz, W.; Ninham, B. W. Hofmeister Effects in Surface Tension of Aqueous Electrolyte Solution. Langmuir 2005, 21, 2619–2623. (4) Parsons, D. F.; Boström, M.; Maceina, T. J.; Salis, A.; Ninham, B. W. Hofmeister Effects: Interplay of Hydration, Nonelectrostatic Potentials, and Ion Size. Langmuir 2011, 13, 12352–12367. (5) Kunz, W., Ed. Specific Ion Effects; World Scientific Publishing: Singapore, 2010. (6) Hofmeister, F. Zur Lehre von der Wirkung der Salze. Arch. f¨ ur Exp. Pathol. und Pharmakologie 1888, 25, 1–30. (7) Lewith, S. Zur Lehre von der Wirkung der Salze. Arch. f¨ ur Exp. Pathol. und Pharmakologie 1887, 24, 1–16. (8) Robinson, R. A.; Wilson, J. M.; Ayling, H. S. The Activity Coefficients of Some Bivalent Metal Nitrates in Aqueous Solution at 25◦ from Isopiestic Vapor Pressure Measurements. J. Am. Chem. Soc. 1942, 64, 1469–1471. (9) Stokes, R. H.; Robinson, R. A. Ionic Hydration and Activity in Electrolyte Solutions. J. Am. Chem. Soc. 1948, 70, 1870–1878. (10) Stokes, R. H.; Robinson, R. A. Electrolyte Solutions - Second Revised Edition; Dover Publications: Mineola, New York, 2002.

27



(11) Debye, P.; H¨ uckel, E. Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. Phys. Z. 1923, 24, 185–206. (12) Triolo, R.; Grigera, J. R.; Blum, L. Simple Electrolytes in the Mean Spherical Approximation. J. Phys. Chem. 1976, 80, 1858–1861. (13) Cartailler, T.; Turq, P.; Blum, L.; Condamine, N. Thermodynamics of Ion Association in the Mean Spherical Approximation. J. Phys. Chem. 1992, 96, 6766–6772. (14) Simonin, J.-P.; Bernard, O.; Blum, L. Real Ionic Solutions in the Mean Spherical Approximation. 3. Osmotic and Activity Coefficients for Associating Electrolytes in the Primitive Model. J. Phys. Chem. B 1998, 102, 4411–4417. (15) Villard, A.; Bernard, O.; Dufrêche, J. F. Non-Additivity of Ionic Radii in Electrolyte Solutions: Hofmeister Effect on Mixtures Modeled by an Associated MSA Model. J. Mol. Liq. 2018, (16) Shannon, R. D. Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides. Acta Crystallogr. Sect. A 1976, 32, 751–767. (17) Jungwirth, P.; Tobias, D. J. Molecular Structure of Salt Solutions: A New View of the Interface with Implications for Heterogeneous Atmospheric Chemistry. J. Phys. Chem. B 2001, 105, 10468–10472. (18) Jungwirth, P.; Tobias, D. J. Specific Ion Effects at the Air / Water Interface. Chem. Rev. 2006, 106, 1259–1281. (19) Dang, L. X.; Chang, T.-M.; Roeselova, M.; Garrett, B. C.; Tobias, D. J. On NO3 –H2 O Interactions in Aqueous Solutions and at Interfaces. J. Chem. Phys. 2006, 124, 066101. (20) Minofar, B.; Jungwirth, P.; Das, M. R.; Kunz, W.; Mahiuddin, S. Propensity of Formate, Acetate, Benzoate, and Phenolate for the Aqueous Solution/Vapor Interface:

28


Page 28 of 34

Page 29 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Surface Tension Measurements and Molecular Dynamics simulations. J. Phys. Chem. C 2007, 111, 8242–8247. (21) Thomas, J. L.; Roeselová, M.; Dang, L. X.; Tobias, D. J. Molecular Dynamics Simulations of the Solution-Air Interface of Aqueous Sodium Nitrate. J. Phys. Chem. A 2007, 111, 3091–3098. (22) Koelsch, P.; Viswanath, P.; Motschmann, H.; Shapovalov, V. L.; Brezesinski, G.; Möhwald, H.; Horinek, D.; Netz, R. R.; Giewekemeyer, K.; Salditt, T. et al. Specific Ion Effect in Physicochemical and Biological Systems: Simulations, Theory and Experiments. Colloids Surfaces A: Physicochem. Eng. Asp. 2007, 303, 110–136. (23) Bley, M.; Duvail, M.; Guilbaud, P.; Dufrêche, J.-F. Simulating Osmotic Equilibria: A New Tool for Calculating Activity Coefficients in Concentrated Aqueous Salt Solutions. J. Phys. Chem. B 2017, 121, 9647–9658. (24) Bley, M.; Duvail, M.; Guilbaud, P.; Penisson, C.; Theisen, J. Molecular Simulation of Binary Phase Diagrams from the Osmotic Equilibrium Method : Vapour Pressure and Activity in Water Ethanol Mixtures. Mol. Phys. 2018, 116, 2009–2021. (25) McMillan, W. G.; Mayer, J. E. The Statistical Thermodynamics of Multicomponent Systems. J. Chem. Phys. 1945, 13, 276–305. (26) Molina, J. J.; Dufrêche, J. F.; Salanne, M.; Bernard, O.; Jardat, M.; Turq, P. Models of Electrolyte Solutions from Molecular Descriptions : The Example of NaCl Solutions. Phys. Rev. E. 2009, 80, 065103. (27) Fyta, M.; Kalcher, I.; Dzubiella, J.; Netz, R. R. Ionic Force Field Optimization based on Single-Ion and Ion-Pair Solvation Properties. 2010, 132, 024911. (28) Nguyen, T. N.; Duvail, M.; Villard, A.; Molina, J. J.; Guilbaud, P.; Dufrêche, J. F. Multi-Scale Modelling of Uranyl Chloride Solutions. J. Chem. Phys. 2015, 142, 024501. 29



(29) Mouˇcka, F.; Nezbeda, I.; Smith, W. R. Chemical Potentials, Activity Coefficients, and Solubility in Aqueous NaCl Solutions: Prediction by Polarizable Force Fields. J. Chem. Theory Comput. 2015, 11, 1756–1764. (30) Yadav, S.; Chandra, A. Preferential Solvation, Ion Pairing, and Dynamics of Concentrated Aqueous Solutions of Divalent Metal Nitrate Salts. J. Chem. Phys. 2017, 147, 244503. (31) Lyubartsev, A. P.; Laaksonen, A. Calculation of Effective Interaction Potentials from Radial Distribution Functions: A Reverse Monte Carlo Approach. Phys. Rev. E. 1995, 52, 3730–3737. (32) Vrbka, L.; Lund, M.; Kalcher, I.; Dzubiella, J.; Netz, R. R.; Kunz, W. Ion-Specific Thermodynamics of Multicomponent Electrolytes: A Hybrid HNC/MD Approach. J. Chem. Phys. 2009, 131, 154109. (33) Smith, W. R.; Moucka, F.; Nezbeda, I. Osmotic Pressure of Aqueous Electrolyte Solutions via Molecular Simulations of Chemical Potentials: Application to NaCl. Fluid Phase Equilib. 2016, 407, 76–83. (34) Luo, Y.; Roux, B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett 2010, 1, 183–189. (35) Kohns, M.; Horsch, M.; Hasse, H. Activity Coefficients from Molecular Simulations using the OPAS Method. J. Chem. Phys. 2017, 147 . (36) Aqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021–8024. (37) Dang, L. X. Computational Study of Ion Binding to the Liquid Interface of Water. J. Phys. Chem. B 2002, 106, 10388–10394.

30


Page 30 of 34

Page 31 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(38) Dang, L. X.; Schenter, G. K.; Fulton, J. L. EXAFS Spectra of the Dilute Solutions of Ca2+ and Sr2+ in Water and Methanol. J. Phys. Chem. B 2003, 107, 14119–14123. (39) Duvail, M.; Guilbaud, P. Understanding the Nitrate Coordination to Eu3+ Ions in Solution by Potential of Mean Force Calculations. Phys. Chem. Chem. Phys. 2011, 13, 5840–5847. (40) Hamer, W. J.; Wu, Y.-C. Osmotic Coefficients and Mean Acitivity Coefficients of Uniunivalent Electrolytes in Water at 25◦ C. J. Phys. Chem. Ref. Data. 1972, 1, 1047–1100. (41) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, III, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H. et al. Amber 14. University of California, San Francisco, 2014. (42) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald : An N log(N ) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. (43) Caldwell, J. W.; Kollman, P. A. Structure and Properties of Neat Liquids Using Nonadditive Molecular Dynamics: Water, Methanol, and N-Methylacetamide. J. Phys. Chem. 1995, 99, 6208–6219. (44) Meng, E. C.; Kollman, P. A. Molecular Dynamics Studies of the Properties of Water around Simple Organic Solutes. J. Phys. Chem. 1996, 100, 11460–11470. (45) Neilson, G.; Enderby, J. Neutron and X-ray Diffraction Studies of Concentrated Aqueous Electrolyte Solutions. Annu. Rep. Prog. Chem., Sect. C Phys. Chem. 1979, 76, 185–220. (46) Skipper, N.; Neilson, G. X-ray and Neutron Diffraction Studies on Concentrated Aqueous Solutions of Sodium Nitrate and Silver Nitrate. J. Phys. Condens. Matter 1989, 1, 4141–4154.

31



(47) Smith, D. E.; Dang, L. X. Computer Simulations of NaCl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757–3766. (48) Smirnov, P.; Yamagami, M.; Wakita, H.; Yamaguchi, T. An X-Ray Diffraction Study on Concentrated Aqueous Calcium Nitrate Solutions at Subzero Temperatures. J. Mol. Liq. 1997, 74, 305–316. (49) Fulton, J. L.; Heald, S. M.; Badyal, Y. S.; Simonson, J. M. Understanding the Effects of Concentration on the Solvation Structure of Ca2+ in Aqueous Solution. I: The Perspective on Local Structure from EXAFS and XANES. J. Phys. Chem. A 2003, 107, 4688–4696. (50) Larentzos, J. P.; Criscenti, L. J. A Molecular Dynamics Study of Alkaline Earth Metal - Chloride Complexation in Aqueous Solution. J. Phys. Chem. B 2008, 112, 14243– 14250. (51) Habenschuss, A.; Spedding, F. H. The Coordination (Hydration) of Rare Earth Ions in Aqueous Chloride Solutions from X-ray Diffraction. III. SmCl3 , EuCl3 , and Series Behavior. J. Chem. Phys. 1980, 73, 442–450. (52) Ishiguro, S.-i.; Umebayashi, Y.; Kato, K.; Takahashi, R.; Kazuhiko, O. Strong and Weak Solvation Steric Effects on Lanthanoid (III) Ions in N,N-Dimethylformamide N,N-Dimethylacetamide Mixtures. J. Chem. Soc., Faraday Trans. 1998, 94, 3607–3612. (53) Clavaguéra, C.; Pollet, R.; Soudan, J. M.; Brenner, V.; Dognon, J. P. Molecular Dynamics Study of the Hydration of Lanthanum(III) and Europium(III) including Many-Body Effects. J. Phys. Chem. B 2005, 109, 7614–7616. (54) Duvail, M.; Spezia, R.; Vitorge, P. A Dynamic Model to explain Hydration Behaviour along the Lanthanide Series. ChemPhysChem 2008, 9, 693–696.

32


Page 32 of 34

Page 33 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(55) Duvail, M.; Vitorge, P.; Spezia, R. Building a Polarizable Pair Interaction Potential for Lanthanoids(III) in Liquid Water: A Molecular Dynamics Study of Structure and Dynamics of the whole Series. J. Chem. Phys. 2009, 130, 104501. (56) Duvail, M.; Ruas, A.; Venault, L.; Moisy, P.; Guilbaud, P. Molecular Dynamics Studies of Concentrated Binary Aqueous Solutions of Lanthanide Salts: Structures and Exchange Dynamics. Inorg. Chem. 2010, 49, 519–530. (57) Caminiti, R.; Licheri, G.; Piccaluga, G.; Pinna, G. On NO− 3 - H2 O Interactions in Aqueous Solutions. J. Chem. Phys. 1978, 68, 1967–1970. (58) Ewing, W. W.; Mikovsky, R. J. Calcium Nitrate. V. Partial Molal Volumes of Water and Calcium Nitrate in Concentrated Solutions. J. Am. Chem. Soc. 1950, 72, 1390–1393. (59) Isono, T. Density, Viscosity, and Electrolytic Conductivity of Concentrated Aqueous Electrolyte Solutions at Several Temperatures. Alkaline-Earth Chlorides, LaCl3 , Na2 SO2 , NaNO3 , NaBr, KNO3 , KBr, Cd(NO3 )2 . J. Chem. Eng. Data 1984, 29, 45–52. (60) Rard, J. A. Osmotic and Activity Coefficients of Aqueous Lanthanum(III) Nitrate and Densities and Apparent Molal Volumes of Aqueous Europium(III) Nitrate at 25 ◦ C. J. Chem. Eng. Data 1987, 32, 92–98. (61) Mart´ınez, L.; Andrade, R.; Birgin, E. G.; Mart´ınez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2010, 30, 2157–2164. (62) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (63) Caminiti, R.; Cucca, P.; D’Andrea, A. Hydration Phenomena in a Concentrated Aque-

33



ous Solution of Ce(NO3 )3 . X-Ray Diffraction and Raman Spectroscopy. Z. Naturforsch. 1983, 38a, 533–539. (64) Onsager, L.; Samaras, N. N. T. The Surface Tension of Debye-H¨ uckel Electrolytes. J. Chem. Phys. 1934, 2, 528–536. (65) Wu, Y. C.; Hamer, W. J. Revised Values of the Osmotic Coefficients and Mean Activity Coefficients of Sodium Nitrate in Water at 25 ◦ C. J. Phys. Chem. Ref. Data. 1980, 9, 513–518. (66) Rard, J. A.; Spedding, F. H. Isopiestic Determination of the Activity Coefficients of Some Aqueous Rare-Earth Electrolyte Solutions at 25 ◦ C. 6. Eu(NO3 )3 , Y(NO3 )3 , and YCl3 . J. Chem. Eng. Data 1982, 27, 454–461. (67) Molina, J. J.; Duvail, M.; Dufrêche, J. F.; Guilbaud, P. Atomistic Description of Binary Lanthanoid Salt Solutions: A Coarse-Graining Approach. J. Phys. Chem. B 2011, 115, 4329–4340.

Graphical TOC Entry

34


Page 34 of 34

Activity Coefficients of Aqueous Sodium, Calcium ... - ACS Publications

Recommend Documents