Simulating Osmotic Equilibria: A New Tool for ... - ACS Publications

Subscriber access provided by the University of Exeter

Article

Simulating Osmotic Equilibria: A New Tool to Calculate Activity Coefficients in Concentrated Aqueous Salt Solutions Michael Bley, Magali Duvail, Philippe Guilbaud, and Jean-Francois Dufrêche J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b04011 • Publication Date (Web): 19 Sep 2017 Downloaded from http://pubs.acs.org on October 1, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 38

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

Simulating Osmotic Equilibria: A New Tool to Calculate Activity Coefficients in Concentrated Aqueous Salt Solutions Michael Bley,∗,† Magali Duvail,† Philippe Guilbaud,‡ and Jean-Fran¸cois Dufrêche∗,† †Institut de Chimie Séparative de Marcoule (ICSM), UMR 5257 CEA - Université Montpellier - CNRS - ENSCM, BP 17171, 30207 Bagnols-sur-Cèze, 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]

1 ACS Paragon Plus Environment


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 Herein, a new theoretical method is presented for predicting osmotic equilibria and activities, where a bulk liquid and its corresponding vapor phase are simulated by means of molecular dynamics using explicit polarization. Calculated time-averaged number density profiles provide the amount of evaporated molecules present in the vapor phase and consequently the vapor phase density. The activity of the solvent and the corresponding osmotic coefficient are determined by the vapor density at different solute concentrations with respect to the reference vapor density of the pure solvent. With the extended Debye-H¨ uckel equation for the activity coefficient and the corresponding Gibbs-Duhem relation in addition the activity coefficients of the solutes are calculated by fitting the osmotic coefficients. A simple model based on the combination of Poisson processes and Maxwell-Boltzmann velocity distributions is introduced to interpret statistical phenomena observed during the simulations, which are related to evaporation and recondensation. This method is applied on aqueous dysprosium nitrate (Dy(NO3 )3 ) solutions at different concentrations. The obtained densities of the liquid bulk, osmotic and activity coefficients are in good agreement with the experimental results for concentrated and saturated solutions. Density profiles of the liquid-vapor interface at different concentrations provide a detailed insight on the spatial distribution of all compounds.

Introduction Systems containing aqueous electrolyte solutions are found in a wide range of industrial, geological, and biological processes, and the understanding of their nature plays a huge role for various applications such as corrosion prevention and hydrometallurgical separation techniques. 1,2 Dysprosium is found in neutron-absorbing control rods in nuclear power plants and as a substitute for neodymium in permanent magnets for electric cars motor and wind 2 ACS Paragon Plus Environment

Page 2 of 38

Page 3 of 38

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


turbine generators and thus plays a major role in the transition process towards carbon-free energy technologies. 3,4 The development of efficient recycling routes for such products allows reducing toxic waste and to overcome possible shortages of Dy in the near future. 5,6 Understanding solvent extraction based purification and recycling approaches requires detailed information on different thermodynamic equilibrium properties of aqueous phases containing electrolytes. The efficiency of liquid-liquid extraction systems at equilibrium conditions can be predicted by using the mass action law. 7–9 In the case of nuclear waste management, the separation of actinoid An3+ and lanthanoid Ln3+ salts mostly in nitric solution requires activity coefficients of all the compounds involved up to high concentrations. Detailed experimental studies on the osmotic and activity coefficients of aqueous solutions for various lanthanoid nitrates from low concentrations until saturation have been conducted by Rard and Spedding using isopiestic measurements. 10–12 For actinoids, only few experimental data is available due to their radioactivity and their complex chemistry. 2,13–15 Since the fundamental work of Debye and H¨ uckel 16,17 (DH) in the 1920s, several different phenomenological improvements have been introduced to predict osmotic and activity coefficients of aqueous electrolyte solutions beyond the limiting law of higher concentrations. In 1972, Hamer and Wu introduced a concept using several empirical terms to surpass the limitation of the original DH equation for concentrated electrolyte solutions. 18 Brønsted and later Guggenheim introduced the Specific ion Interaction Theory (SIT) to estimate the single-ion activity coefficients at higher concentrations. 19,20 Pitzer evolved this theory and provided the pragmatic concept of the Pitzer equations, which are still widely used in science and industry. 21–23 Methods based on statistical mechanics, such as the Mean Spherical Approximation (MSA) theory 24,25 and Concept of Hypernetted Chains (HNC), 26,27 use integral equations to predict thermodynamics properties. Most of those models are based on continous solvent model as they directly extend DH theory. They are justified by the McMillan-Mayer (MM) theory. 28,29 The latter shows that they correspond to a double step averaging process (i) average over the solvent molecules, which gives a continuous solvent



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

model for the solvents, and (ii) average over the solute molecules. Similarly, several approaches have been proposed to calculate osmotic coefficients in solutions from molecular simulation. These methods to calculate solute and solvent activities of aqueous solutions are described in Figure 1. The MacMillan-Mayer method consists in calculating the MM potential of mean forces between the ions (first stage of the MM theory). In a second stage the configurations of the solute molecules are averaged out on order to obtain the osmotic pressure and the corresponding activity coefficients. 30–33 This method is especially interesting as it automatically satisfies the DH limiting law as long as the MM potential is coulombic at long distances. Nevertheless, it is restricted to relatively dilute solutions because N-body terms in the MM potentials are generally ignored. Another advantage of the MM approach is that equilibrium constants K ◦ (T ) can be calculated because they are direct integrals of the MM potentials. 34 The second method, the Kirkwood-Buff approach, 35–38 requires the radial distribution functions g(r) in the grand canonical ensemble. Then an integration scheme has to be performed in order to get the thermodynamic properties. This intergration step is difficult because of the high precision required and the computational cost to calculate a radial distribution function g(r) in the grand canonical ensemble. Due to this bottleneck, to our knowledge despite its recent interest, the Kirkwood-Buff theory has not been applied in practice to derive activities in solution. Weerasinghe and Smith simulated activity derivatives of aqueous alkali halide solutions using Kirkwood-Buff integrals. 39,40 Another way to obtain activities in solutions is based on grand canonical Monte-Carlo simulations (MC). 30,41,42 Those MC simulations allow the calculation of individual ion activities and the prediction of solubilities at a low noise. This method can be considered as a very flexible approach, but it requires a suitable MC code and has to be conducted in the grand canonical ensemble. A similar MC methodology is the one of Gibbs ensemble, were phase equilibria in fluids without the presence of an interface are simulated. 43,44 The vapor pressure is directly obtained from the method by simulating the transfer of particles between a liquid and a gas phase. The Gibbs Ensemble Monte Carlo (GEMC) method has been applied to


Page 4 of 38

Page 5 of 38

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


binary solvent mixtures, 45,46 but to our knowledge, not to electrolyte solutions. A fourth method directly calculates the osmotic pressure as a surface force is applied on an aqueous electrolyte solution separated from pure water by an ideal semi-permeable membrane. 47 The applied force on the membrane directly leads to the osmotic pressure of the solution, but also causes a lot of noise and the simulation requires the calculation of the force on an immobile object (the membrane). The new method presented here considers osmotic equilibrium of the solvent between the solution and its adjacent vapor phase. No membrane is needed. The calculation of the activity coefficient is possible because the vapor phase allows the measurement of the chemical potential of the solvent. It is the simulation version of the isopiestic measuresments of activity coefficients conducted by Rard et al. 10–12 It can also be applied on mixtures via the thermodynamic integration of the Gibbs-Duhem relation. It is also possible to describe non-polar solvent. One of the major advantages of the methods presented here is that it can be applied on any given system for any simulation code as long as a liquid-vapor equilibrium can be simulated. Current MD studies of liquid-vapor interfaces cover structural effects and interfacial density profiles, 48–51 thermodynamic properties, 52,53 and the statistics of evaporation for Lennard-Jones (LJ) fluids. 54–56 In this article the method is developed for concentrated aqueous lanthanoid salt solutions. The hydration properties of lanthanoid cations (Ln3+ ) have been part of various experimental 10,57–61 and theoretical studies 62–68 within the last decades. For Dy3+ and other rare earth cations, the behavior in aqueous − solutions containing chloride (Cl− ), nitrate (NO− 3 ) and perchlorate (ClO4 ) anions has been

well investigated. 67,68 However, to the authors best knowledge, MD simulations of liquidvapor equilibria have never been conducted to directly access the osmotic and activity coefficients of aqueous electrolyte solutions at different salt concentrations. Molecular dynamics simulations of the bulk liquid phase and the liquid-vapor interfaces of aqueous dysprosium nitrate salt (Dy(NO3 )3 ) solutions have been conducted at different concentrations. The new method based on MD simulations with explicit polarization and taking advantage of osmotic



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

(a) McMillan-Mayer: aw asalt

VMM Bulk Liquid

(b) Kirkwood-Buff:

g(r)

aw asalt

Bulk Liquid

(c) Grand Canonical Monte-Carlo: aw asalt

Nsalt(μi), Nwater(μi) Bulk Liquid

(d) Osmotic Membrane:

𝐅

aw asalt

FSurface

Membrane

(e) Liquid-Vapor Equilibrium: 𝜌w ∗ 𝜌w

Nvapor

aw asalt Vapor

Liquid - Vapor Interface

Figure 1: Schematic representations of various simulation methods used for accessing water and salt activities (a) McMillan-Mayer, (b) Kirkwood-Buff, (c) Monte-Carlo, (d) Osmotic membrane, and (e) Liquid-Vapor Equilibrium .

equilibria is introduced and the statistics of evaporation, time-averaged density profiles and hydration properties of the involved species are especially taken into account in order to quantify the accuracy of the method.


Page 6 of 38

Page 7 of 38

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 osmotic coefficient of a solvent φS (bE ) at a certain molality of the solute bE is related to the activity coefficient of the solute γE (bE ) through the Gibbs-Duhem relation. 2 Equation 1 provides a relation between both coefficients by Z

b=bE

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

φS (b) − 1 db , b

(1)

based on an integration of the molality b up to a certain solute concentration bE . Both expressions ln γE (bE ) and φS (b) are generally fitted from an extended form of the DebyeH¨ uckel equation. 18 The osmotic coefficient of the solvent is obtained by √ p p |z+ z− |Ab k × bE 1 √ φS (bE ) = 1 − (1 + Bb k × bE ) − 2 ln(1 + Bb k × bE ) − (Bb )3 k × bE 1 + Bb k × bE j X i + βi × (bE )i , i + 1 i=1 (2) where Ab and Bb are DH parameters on a molality scale, |z+ z− | is the absolute product of ionic charges, and k × bE is an expression for the ionic strength. The sum in the last term adds fitting parameters βi up to the i-th order. The corresponding activity coefficient formula 18 gives

√ j |z+ z− |Ab k × bE X √ βi × (bE )i , ln γE (bE ) = − + 1 + Bb k × bE i=1

(3)

The first term in Equation 3 equals the extended Debye-H¨ uckel equation for the activity coefficient and covers the dilute regime of electrolyte solutions. The second term sums up additional fitting terms up to the i-th order, which are necessary to describe the activity coefficient at high solute concentrations bE . Equation 4 is used to obtain the parameter Ab



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 8 of 38

from the DH theory Ab =

e3

p

2πNA ρS (T ) 3/2

3/2

(4π0 )3/2 × kB × T 3/2 × r

,

(4)

where ρS (T ) is the density of the solvent at temperature T , NA is the Avogadro constant, e is the elemental charge, 0 is the vacuum permittivity, r is the relative permittivity of the solvent, and kB is the Boltzmann constant. The parameter Bb can be evaluated with Equation 5 s Bb = rIon,aq × Bb∗ = rIon,aq ×

2e2 NA ρS (T ) , kB T r 0

(5)

where rIon,aq is the effective hydrated radius of the ions. rIon,aq depends on the Debye-H¨ uckel parameters Ab and Bb and is accessible by fitting experimental and simulated data of the two coefficients.

Determination of osmotic coefficients from osmotic equilibria For an aqueous electrolyte solution in osmotic equilibrium with the solvent vapor, the osmotic coefficient φS (bE ) is defined from the activity of the solvent at a certain solute concentration aS (bE ) as shown in Equation 6

φS (bE ) =

− ln aS (bE ) , νMS bE

(6)

where ν is the amount of ions per formula unit of solute and MS is the molar mass of the solvent. The ratio of the partial vapor pressure at a certain concentration pS (bE ) over the standard vapor pressure p S provides the activity of the solvent aS (bE )

aS (bE ) =

pS (bE ) nS (bE ) ρS (bE ) = = , pS nS ρ S

(7)

where nS (bE ) and ρS (bE ) are the corresponding number and mass gas phase density at the solute concentration bE within a given gas phase volume, respectively. n S and ρS are the

number and mass density in the gas phase for the pure solvent, respectively. Vapor densities 8 ACS Paragon Plus Environment

Page 9 of 38

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


and vapor pressure provide equivalent expressions for the activity since the equation of state for the gas phase corresponding to our water model used in this work (POL3 model 69,70 ) follows the ideal gas law in good agreement (Figure S1 in the Supporting Information). The reference vapor density ρ S simulated thus needs to be precisely determined for every volatile compound taking into account the statistics of the simulations as described in the statistics subsection. For the POL3 model of water, the vapor pressures simulated exceed the experimental vapor pressures by a factor of at least two. However, the offset remains constant when increasing the solute concentration. This provides solvent activities aS (bE ) very similar to those observed in experiments. 11 The number density of solvent molecules present in the gas phase nS,i and the vapor mass density for the i-th run for a fixed concentration can be calculated from atomistic simulations when the liquid-vapor interface is simulated. The vapor phase density ρS,i is defined by

ρS,i =

nS,i MS hNS,i,gas iMS = , NA Vgas,i NA

(8)

where hNS,i,gas i is the average amount of solvent molecules in the gas phase and Vgas,i is the gas phase volume of the i-th production run. Equation 9 describes the mean amount of solvent molecules present in the vapor phase for the time span from t0 = 0 to t t KX hN S,i,gas (t)i = NS,i,gas (t) , n t0 =0

(9)

where K is a norm factor with respect to the interfacial areas and NS,i,gas (t) the amount of solvent molecules present in the gas phase at the n-th recorded step for the i-th run. The average amount of solvent molecules in the vapor phase hNS,i,gas i is obtained by calculating the arithmetic mean of NS,i,gas (t), while neglecting the first several nanoseconds due to the equilibration process. For an enhanced precision of the method, the outcome of j different trajectories for each solute concentration is used leading to j independent values of hNS,i,gas i



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 10 of 38

and ρS . Equation 10 provides the average vapor density j

ρ S

1X = ρ j i=1 S,i

j

and

1X ρS (bE ) = ρS,i (bE ) . j i=1

(10)

The amount of different trajectories j is usually higher for the reference state than for a given solute concentration. The resulting standard deviation σρS related to each density ρS (bE ) and ρ S is used for a selective criterion for all ρS (bE ) ± σρS (bE ) with respect to the standard vapor density ρ S |ρ S − ρS (bE )| > σρS (bE ) ,

(11)

which excludes all results that are too close to the reference vapor pressure of the solvent. This criterion is necessary to avoid unreliable results for activities and osmotic coefficients due to the lack of sensitivity of the method at low solute concentrations. Only results within a sufficient distance to the reference value are used for fitting the curves towards the simulated results. Equation 11 usually neglects the results in the diluted domain, since the changes in the vapor pressure are smaller than the sensitivity of the method. This shortcoming is covered by the advantageous properties in the first part of the extended Debye-H¨ uckel concept (Equations 2 and 3), which are used for fitting the remaining data. This constrained fitting approach allows neglecting the results simulated for dilute solutions.

Poisson model used for Statistics of Evaporation The method introduced herein and its precision can be understood thanks to a Poisson model. Statistical effects are observed during the analysis of molecular dynamics simulations of liquid-vapor equilibria. 54–56 For different initial geometries of the j production runs, slightly different mean amounts of solvent in the gas phase NS,i,gas (t) are obtained due different quantities of solvent molecules S evaporated and a varying presence time in the gas phase τx . This leads to a time-dependent error σNS,gas (t) for the amount of solvent in the gas phase NS,gas (t), since the quantity of gaseous solvent molecules in the i-th run NS,i,gas (t) out of j 10 ACS Paragon Plus Environment

Page 11 of 38

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


runs at a give time t differs. The corresponding relation is j

1X NS,gas (t) = NS,i,gas (t) . j i=1

(12)

Due to the use of periodic boundary conditions in MD simulations, each evaporated solvent molecule moves through the gas phase within a transition time τx until it recondenses in the same bulk liquid. A solvent molecule, which arrives at any lateral plane, is translated to the opposite plane. Thus, a statistical model is developed for predicting the time-dependent error of the method σN,gas (t) neglecting the very rare collisions of molecules in the gas phase. Neglecting crashes allows the description of the evaporation by a Poisson process and the velocity distribution of the molecules in x-direction by a Maxwell-Boltzmann velocity distribution. The Poisson probability distribution function for emitting k solvent molecules within the time of simulation t is P (k|λ, t) =

(λ × t)k −λ×t e , k!

(13)

where λ is the evaporation or phase-transition rate of the solvent molecules for a given interval ∆t with t = n × ∆t. This rate decreases with an increasing salt content in the aqueous phase (λ(bE ) ∝ aW (bE )). The corresponding number of transferred solvent molecules hk(t)i during time t is calculated through hk(t)i = λ × t ,

(14)

and the time-dependent standard deviation σhk(t)i is given by

σhk(t)i =

√

λ×t.


(15)


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 12 of 38

The velocity probability distribution p(vx ) for a movement of the particles in x-direction is independent from the duration of the simulation mS vx2 mS × vx exp − , p(vx ) = kB T 2kB T

(16)

but depends on the temperature T and the mass of the evaporated solvent mS . The mean time of flight in x-direction hτx i is obtained by integrating this probability distribution 1 hτx i = lx × h i = lx × vx

Z

vx =∞

vx =0

1 p(vx )dvx = lx × vx

r

πMS 2RT

(17)

where lx is the total length of the gas phase Vgas in x-direction. h1/vx i is the average of the inverse of the velocity and it should be distinguished from the inverse of the average of p the velocity 1/hvx i = (2MS )/πRT . After a simulation time t = n × ∆t it is considered that the total amount of evaporated solvent molecules hk(t)i during time t lies within the confidence interval hk(t)i −

σhk(t)i σhk(t)i < hk(t)i < hk(t)i + . 2 2

(18)

On average, a solvent molecule stays during a time hτx i in the gas phase. Consequently, the mean amount of solvent molecules in the gas phase hNS,gas i is

hNS,gas i =

hτx i hτx i × hk(t)i = × λ × t = hτx i × λ t t

(19)

Equation 19 shows that the expected mean amount of solvent molecules in the gas phase is constant, if each solvent molecule evaporated returns to the bulk after a transfer time τx . From the Equations 15 and 18 the confidence interval of hNS,gas (t)i is obtained √ √ hτx i × λ hτx i × λ √ √ hNS,gas (t)i − < hNS,gas (t)i < hNS,gas (t)i + . 2 t 2 t


(20)

Page 13 of 38

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


Thus, the standard deviation of hNS,gas i reads √ hτx i × λ √ σNS,gas (t) = , t

(21)

which approaches zero for large simulation times t. It should be noted that the standard deviation of the mean time of flight στx is neglected. This is valid for sufficient large simulation times t. Subsequently, the relative error of the method δNS,gas (t) is given by

δNS,gas (t) =

σNgas (t) 1 =√ , hNgas i t×λ

(22)

It is possible to simulate the behavior of NS,i,gas (t) for i virtual runs for a given evaporation rate λ for any molar mass MS with the Poisson model. A more detailed section about the statistics of this approach and the Poisson method is found in Supporting Information. The relative error of the model in Equation 22 is inversely proportional to the square root of the length of the trajectory t and of the emission rate λ. The emission rate is proportional to the interface area. Therefore, the precision can be increased either by performing longer simulations or by simulating an enlarged interfacial area. If the thickness of the liquid phase is constant, the area is proportional to N . All the MD algorithms are at the very best proportional to N . Thus, if the computational time is multiplied by a factor α, in that √ idealized case, λ is also multiplied by α and the relative error is divided by α. The same effect on the precision can be obtained by multiplying the simulation time by the same factor α. This later route to improve the accuracy of the method is in fact much better because most of the MD programs actually scale proportionally to t without any problem, contrary to the N proportionality, which cannot be practically obtained. Consequently, considering the difficulty to scale MD algorithm proportional to N , the best approach to improve the precision consists in increasing the simulation time and not the size of the interface. The thickness of the liquid phase and the interfacial area have to be as small as possible.



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

Simulation details Molecular dynamics Classical molecular dynamics simulations of aqueous dysprosium nitrate salts at different concentrations ranging from 0.5 to 6.5 mol kg−1 and in pure water have been carried out with SANDER14, a module of AMBER14, 71 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. 72 The systems were equilibrated at 298.15 K and 1 bar (0.1 MPa) for at least 2 ns, and production runs were afterwards conducted for at least 15 ns. All atomic coordinates were written to the trajectory file every ps. Here, the van der Waals energy is described by a 12-6 Lennard-Jones potential. The water molecules were described by the rigid POL3 model, 69,70 which takes into account the polarization. Lennard-Jones parameters for the Dy3+ cation were set to reproduce the experimental hydration properties 67 (Table 1). Atomic polarizability used for the Dy3+ ion has been calculated by Clavaguéra et al. 63 . This polarizable force field for describing the Dy3+ cation in water provides structural and dynamical properties in good agreement with experimental 57,59,60 and theoretical 62–68 results. The NO− 3 anion was described by the polarizable model defined in our previous study, 67 providing structural properties of NO− 3 in aqueous solutions in reasonable agreement with published experimental 73 and calculated data. 74 Within 3.8 ˚ A of NNO−3 around six water molecules are found in the first hydration shell at low salt concentration (0.5 mol kg−1 ). Table 2 lists the compositions and the dimensions of the simulation boxes created with the PACKMOL package. 75 The equilibration stage was conducted for 5.0 to 10.0 ns in the NPT canonical ensemble at 298.15 K and 1.0 bar using a τP of 0.1 ps, where τP is the relaxation time of the Berendsen pressure coupling. 76 Several production runs have been performed by changing either the initial coordinates of the atoms or the degree of parallelization. The amount of threads was


Page 14 of 38

Page 15 of 38

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 1: Parameters used for the MD simulations

a

atom ii a σii b qi c αd Dy3+ 0.330 3.118 +3.000 0.728 NNO−3 0.711 3.25 -0.632 0.530 ONO−3 0.879 3.067 +0.896 0.434 HH2 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

Table 2: Characteristics of the MD simulation boxes bE a 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 5.0 6.0 6.5

N H2 O b 3300 3552 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000

N Dy3+ c 27 54 81 108 135 162 189 216 270 324 351

N NO−3 d 81 162 243 324 405 486 567 648 810 972 1053

tEq. e 10 5 10 10 10 10 10 10 10 10 10 10 10

tProd. f 15 15 15 15 15 15 15 15 15 15 15 15 15

#Runs g 3 5 3 3 3 3 3 3 3 3 3 3 3

ρSim. h 0.996 0.996 1.132 1.253 1.364 1.466 1.560 1.648 1.722 1.795 1.916 2.011 2.056

ρExp. i 0.997 0.997 1.139 1.268 1.388 1.496 1.594 1.682 1.765 1.840 1.976 2.110 2.177

V Bulk j 50.3 × 52.5 × 50.4 × 50.7 × 50.8 × 51.6 × 53.3 × 54.0 × 54.1 × 54.2 × 55.6 × 57.0 × 57.7 ×

50.3 52.0 50.1 51.0 51.2 52.4 53.0 52.9 53.5 54.0 55.6 57.3 57.4

× × × × × × × × × × × × ×

50.4 53.0 49.4 50.4 50.9 52.2 52.1 53.4 54.1 54.0 55.4 57.2 57.6

Dy(NO3 )3 concentration on a molality scale bE in mol kg−1 . b Number of water molecules. c e Number of Dy3+ cations. d Number of NO− 3 anions. Equilibration time in ns. f Production time in ns. g Number of production runs. h Density in g cm−3 calculated from MD simulations in the NPT ensemble. i Experimental density in g cm−3 from Reference. 11 j Liquid bulk volume (lx × ly × lz ) in ˚ A3 .

a

varied for the five runs of the box containing 3552 water molecules. For all other simulation boxes, three different equilibrated coordinate files of an identical composition have been used for each salt concentration. Subsequently, the simulation boxes have been extended to a finite size of 495.0 ˚ A in x-direction to obtain the liquid-vapor interfaces leading to a centered equilibrated bulk liquid phase and two adjacent vacuum areas with a width of around 225 ˚ A each (Figure 2). Note that this vacuum space will hereafter be referred as the vapor phase. Similar partitioning between liquid and gaseous phase is found in the work of Taylor et al. 48



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 16 of 38

y

495 Å

x z ~ 50 Å

~ 50 Å : H2O : NO3Vapor

Liquid

Vapor

: Dy3+

Figure 2: Schematic representation of a simulation box representing an aqueous Dy(NO3 )3 salt solution (1.0 mol kg−1 ) in equilibrium with water vapor. The dimension of the yz-plane depends on the composition of the bulk and the length in x-direction is 495.0 ˚ A for every bulk composition.

and Yuet et al. 52 Those geometries are used to produce liquid-vapor equilibria in the NVT ensemble for 15.0 ns at 298.15 K. The cut-off radius for the all equilibration and production runs is 12.0 ˚ A.

Data analysis All the data analysis has been performed using the 41 production runs (Table 2). Equilibrium between the liquid and the vapor phase is obtained after around 5 ns for each production run, therefore only the results of the last 10.0 ns are taken into account for calculating the mean amount of evaporated water hNS,i,gas i for all concentrations. The determination of hNS,i,gas i requires the locations of the two liquid vapor interfacial yz-planes perpendicular to the x-axis. The two interfacial planes are located at 200 and 295 ˚ A in x-direction, which provides two gas phases of identical volume. All oxygen atoms of water molecules (OH2 O ) that are within the defined range from 0 to 200.0 and 295.0 to 495.0 ˚ A are counted for each of the 104 frames of a trajectory and the arithmetic mean is calculated. The slight 16 ACS Paragon Plus Environment

Page 17 of 38

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


variation of the size of the yz-plane caused by the different used initial atomic coordinate files is taken into account during the calculation of the vapor phase density ρS,i of each run i. Time averaged density profiles of all compounds involved of the selected production runs are obtained from the coordinate files by calculating the arithmetic mean of a given compound within an interval of 1.0 ˚ A.

Results and Discussion Liquid phase The comparison of the densities at different concentrations after the equilibration process with the experimental densities of aqueous Dy(NO3 )3 solutions 77 shows a good agreement between the experiments and the simulations for dilute solutions (Table 2). With an increasing salt content, the densities calculated from MD simulations are underestimated compared to the experimental values extrapolated in the concentrated domain. Table 3 provides the structural properties of the solutes like the first maxima of the radial distribution functions (RDFs) and coordination numbers (CNs) as a function of the salt concentration. Increasing the Dy(NO3 )3 concentration leads to a remarkable decrease of the ion-ion distances with respect to the first maximum peaks of the RDF rNN −NN and rDy−Dy , respectively. Less water molecules are available per ion pair and thus the ion pair density in the liquid phase increases. The simulated hydration properties of Dy3+ at salt concentrations below 1.5 mol kg−1 are in agreement with the experimental results, 57,59,60 since 7.9 water molecules and 0.1 nitrate anions are found in its first coordination shell (Table 3). Increasing the ionic strength and therefore decreasing the water-to-ion pair ratio changes the hydration prop−1 erties of Dy3+ . Almost three water molecules are replaced by three NO− 3 at 6.5 mol kg ,

whereas only a slight decrease of the first rDy−OW , rDy−ON , rDy−NN , and rNN −OW maximum peaks is observed. This behavior, i.e., the formation of Contact Ion Pair (CIP), has already been observed experimentally for lanthanoid salts. 10 On the contrary, the distance between 17 ACS Paragon Plus Environment


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 18 of 38

Table 3: Hydration properties of Dy(NO3 )3 in aqueous solution at 298.15 K calculated at different molalities of salt bE a 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 5.0 6.0 6.5

rDy−OW b 2.40 2.40 2.39 2.39 2.39 2.39 2.39 2.39 2.39 2.38 2.38

CNDy−H2 O c 7.9 7.9 7.9 7.8 7.8 7.6 7.6 7.3 6.7 5.6 5.2

rDy−NN d 3.64 3.63 3.63 3.62 3.62 3.62 3.62 3.61 3.61 3.61 3.61

rDy−ON e 2.44 2.43 2.43 2.43 2.43 2.42 2.42 2.42 2.42 2.41 2.41

CNDy−NO3 f 0.1 0.1 0.1 0.2 0.2 0.4 0.4 0.7 1.3 2.4 2.7

rNN −OW 3.65 3.64 3.64 3.63 3.63 3.63 3.63 3.63 3.64 3.67 3.67

g

rNN −NN h 5.91 5.74 5.66 4.13 4.10 4.07 4.05 4.05 4.04 4.05 4.07

rDy−Dy i 10.9 10.5 10.3 8.66 8.52 8.45 8.37 8.27 8.09 6.48 6.52

Dy(NO3 )3 concentration on a molality scale bE in mol kg−1 . b First maximum peak of the Dy-OW RDF in ˚ A. c Number of water molecules in the Dy3+ first coordination sphere. d First maximum peak of the Dy-NN RDF in ˚ A. e First maximum peak of the Dy-ON RDF in ˚ A. f Number of nitrate anions in the Dy3+ first coordination sphere. g First maximum peak of the NN -OW RDF in ˚ A. h First maximum peak of the NN -NN RDF in ˚ A. i First j maximum peak of the Dy-Dy RDF in ˚ A. First maximum peak of the OW -OW RDF in ˚ A. Subscripts N and W for nitrate and water, respectively.

a

two water molecules rrOW −OW increases with the solute concentration due to the lower amount of available solvent molecules. More Dy(NO3 )3 salt dissolved leads to a transition from an aqueous solution with properly hydrated ions towards a more salt-like structure, where the solvent (H2 O) is more and more replaced by NO− 3 ions in the cation hydration shell. Plotted RDFs and CNs are part of the Supplementary Information (Figures S2 to S9). The influence of the water-to-ion pair ratio on the structure of the aqueous phase and the role of a liquid-vapor interface is shown in a selection of time-averaged mass density profiles along the x- and y-axes of the simulation boxes (Figures 3, 4, and 5). Detailed time-averaged mass density profiles for all the concentrations investigated are part of the Supplementary Information (Figures S10 to S20). The density profiles in y-direction correspond to the behavior of a bulk liquid due to the periodicity of the simulation boxes, whereas in x-direction the density profiles are strongly influenced by the presence of two liquid-vapor interfaces. At


rOW −OW j 2.74 2.74 2.74 2.75 2.76 2.77 2.80 2.82 2.84 2.85 2.86 2.86

Page 19 of 38

1.0

Ρ g cm-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


0.5 mol kg-1

HaL

HbL

0.8

Dy3+ OH2 O NNO3

0.6

-

0.4 0.2 0.0 225

dx Å 245

265

5

15

dy Å 25

35

45

Figure 3: Time-averaged density profiles along the (a) x-axis and (b) y-axis for Dy3+ (green), OH2 O (blue), and NNO−3 (red) in a 0.5 mol kg−1 aqueous solution of Dy(NO3 )3 .

low concentrations (until 3.5 mol kg−1 ) only water molecules are located at the interface and hydrated ions are found at a given distance with respect to the interface (Figures 3 and 4). At 0.5 mol kg−1 no preferred position of the ions in the water is observed, since enough water molecules is available for the hydration of the ion pairs (Figure 3). A thick layer of several ˚ A is found at the interface. As the amount of available water per ion pair decreases, a well defined structure is observed for all species. For 3.5 mol kg−1 a thinner water layer is still found at the interface, but the profile in x-direction exhibits a pronounced structuring compared to the rather flat density profile in y-direction (Figure 4).Beyond the interfacial water layer, a well-defined pattern of preferred positions with respect to the interface is observed for all species. Such a structured system is the product of a reduced amount of solvent molecules per ion pair and allows sharing the solvent molecules between different ions. This structure is energetically more favored than having ions directly at the interface. In the y-direction no explicit pattern is observed when no liquid-vapor interface is present. Finally, at a very high concentration Dy(NO3 )3 , e.g. 6.5 mol kg−1 , the distinct structure around the charged species and the interfacial water layer disappears, as shown in Figure 5. The presence of an interface in x-direction does not change the shape of the resulting density profile compared to the results in the y-direction. Only 2.1 water molecules are available for a single Dy(NO3 )3 salt, implying thus that ionic interactions now dominate. The closer the ions are 19 ACS Paragon Plus Environment


Ρ g cm-3

1.0

3.5 mol kg-1

HaL

HbL

0.8 0.6 0.4

Dy3+ OH2 O NNO3

0.2

-

0.0 225

dx Å 245

265

5

15

dy Å 25

35

45

Figure 4: Time-averaged density profiles along the (a) x-axis and (b) y-axis for Dy3+ (green), OH2 O (blue), and NNO−3 (red) for a 3.5 mol kg−1 aqueous solution of Dy(NO3 )3 .

1.0

Ρ g cm-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

Page 20 of 38

6.5 mol kg-1

HaL

HbL

0.8 0.6 0.4

Dy3+ OH2 O NNO3

0.2

-

0.0 225

dx Å 245

265

5

15

dy Å 25

35

45

Figure 5: Time-averaged density profiles along the (a) x-axis and (b) y-axis for Dy3+ (green), OH2 O (blue), and NNO−3 (red) in a 6.5 mol kg−1 aqueous solution of Dy(NO3 )3 .

packed in the bulk liquid, the more water molecules are available for the hydration shell. Missing water molecules in the Dy3+ hydration first shells are now replaced by oppositely charged ions causing a more and more salt like structure disturbed by the remaining water molecules. The snapshots of the liquid phases for 0.5, 3.5 and 6.5 mol kg−1 in x-direction in contact with two vapor phases are shown in Figure 6. The initial thick water layer at the interface thins for an increasing salt content. A transition from an aqueous electrolyte solution (0.5 mol kg−1 ) with a constant distribution of the ions to a more structured system (3.5 mol kg−1 ) with preferred ion positions towards a more chaotic and more salt-like bulk


Page 21 of 38

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


: H2O

:

(a)

(b)

(c)

NO33+

: Dy : Liquid : Vapor y

z x

3.5 mol kg-1

0.5 mol kg-1

6.5 mol kg-1

Figure 6: Parts of snapshots of the simulation boxes representing the liquid phase (blue background) in contact with the two vapor phases (green background) for (a) 0.5 mol kg−1 , (b) 3.5 mol kg−1 , and (c) 6.5 mol kg−1 aqueous Dy(NO3 )3 solutions.

phase (6.5 mol kg−1 ) is observed. Density profiles and hydration properties for systems with an increased x-dimension of the liquid phase are part of the Supporting Information (Tables S1 to S3 and Figures S21 to S24).

Water in the vapor phase The quantity of water transferred from the liquid to the gaseous state at equilibrium conditions depends on the salt concentration in the liquid phase. However, the transition time from the liquid phase through the gas phase back to the liquid phase only depends on the temperature and the size of the vapor phase. In order to understand the statistics of this osmotic equilibrium, data sets (8) from the Poisson model algorithm were used to model the results obtained from MD simulations of individual production runs (8) for pure water (Figures 7 and 8). This showed that each run converges towards a slightly different mean amount of water in the gas phase due to statistical effects related to the evaporation process. The input paramters for the Poisson model, the mean residence time of a solvent molecule in vapor hτx i and the evaporation rate λ, are derived from the results of the molecular dynamics simulations for pure water. The mean time of flight hτx i is in both approaches 135.2 ps. The mean amount of water present in the gas phase hNS,gas i from the MD simulations 21 ACS Paragon Plus Environment


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

Figure 7: The mean amount of water in the vapor phase NS,i,gas (t) obtained from Molecular Dynamics (MD) and the Poisson model and the corresponding confidence intervals (colored areas limited by dashed lines).

consequently allows to calculate the emission rate λ with Equation 19. For a simulated hNS,gas i = 2.29 , an emission rate λ = 0.017 ps−1 is obtained. Using the results for λ and hτx i from MD simulation in the Poisson model consequently provides an average of 2.32 water molecules present in the gas phase. Figure 7 shows the comparison of the mean amount of water and its variation σNS,gas (t) at a given time t for the results obtained from both molecular dynamics and the Poisson model. The similarity of the blue- and red-colored curves shows that the mean amount of water evaporated in the vapor is described by the Poisson model as introduced in statistics section. Thus, the results obtained by MD results are reproduced by (i) a combination of a Poisson process for describing the phase transfer and (ii) a Maxwell Boltzmann distribution for determining the velocities of the solvent molecules in the gas phase. Both curves exhibit a variance due to the statistical nature of this process. Figure 8 shows a decrease for the standard deviation of the number of water molecules in the gas phase σNS,gas (t) for MD simulation and the Poisson model. Equation 21 with λ = 0.017 ps−1 and hτx i = 135.2 ps provides a good approximation compared to the results of MD simulation and Poisson model. The error for both, molecular dynamics and the Poisson model rapidly decreases within the first 5 ns, but remains then almost constant for the last 10 ns due to the mathematical nature


Page 22 of 38

Page 23 of 38

1.50 1.25

ΣNS,gas HtL

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


Molecular Dynamics Poisson model Equation 21

1.00 0.75 0.50 0.25 0.00

0

2

4

t ns

6

8

10

12

14

Figure 8: Standard deviations σNS,gas (t) from Molecular Dynamics and the Poisson Model (solid line) and the corresponding approximation (dashed line) based on the Poisson model (Equation 21) for λ = 0.017 ps−1 and hτx i = 135.2 ps.

of the Equations 21 and 22, which allows the calculation of the precision of the method for a known emission rate. The emission rate λ = 0.017 ps−1 derived from simulations leads to a precision of 90 % after 5.9 ns. For 15 ns, an accuracy of 94 % is obtained, and 590 ns of production time are necessary to minimize the relative error with a precision of 99%. The corresponding relative error of the method (Equation 21) can be decreased either by drastically increasing the production time t or the phase transition rate λ. An increase of λ for a given salt concentration requires a larger liquid vapor interface, which consequently requires bigger simulation boxes. As already mentioned, increasing the total simulation time is more efficient to minimize the error. An increased precision of this method comes along with a drastic increase of computational cost.

Water activities, osmotic and activity coefficients Molecular dynamics simulations provide the vapor phase mass densities ρS,i,gas , which is used to calculate the activity of water in the liquid phase aS,i,liq (bE ) and the corresponding osmotic coefficients φS,i,liq for each of the i-th production runs. The resulting arithmetic means and standard deviations of those values are given in Table 4. The resulting vapor phase mass densities at concentrations up to 2.0 mol kg−1 are very similar compared to 23 ACS Paragon Plus Environment


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 4: Gas phase properties of aqueous Dy(NO3 )3 simulated at 298.15 K. bE a ρS,gas ± σρS,gas b aS,liq ± σaS,liq c φS,liq ± σφS,liq d 0.0 0.066 ± 0.007 0.5 0.068 ± 0.007 1.027 ± 0.110 -0.640 ± 2.909 1.0 0.066 ± 0.006 0.994 ± 0.088 0.126 ± 1.227 1.5 0.070 ± 0.004 1.057 ± 0.064 -0.505 ± 0.570 2.0 0.053 ± 0.007 0.805 ± 0.107 1.549 ± 0.946 2.5 0.053 ± 0.006 0.805 ± 0.094 1.231 ± 0.662 3.0 0.049 ± 0.006 0.741 ± 0.088 1.407 ± 0.565 3.5 0.044 ± 0.003 0.663 ± 0.049 1.636 ± 0.294 4.0 0.040 ± 0.004 0.600 ± 0.063 1.792 ± 0.375 5.0 0.034 ± 0.002 0.517 ± 0.029 1.835 ± 0.154 6.0 0.024 ± 0.008 0.370 ± 0.114 2.368 ± 0.670 6.5 0.022 ± 0.004 0.332 ± 0.066 2.382 ± 0.409 a Dy(NO3 )3 concentration on a molality scale bE in mol kg−1 . b Vapor phase mass density ρS,gas and the standard deviation σρS,gas in g m−3 . c Water activity in the liquid phase aS,liq and the standard deviation σaS,liq . d Osmotic coeffcient of water in the liquid phase φS,liq and the standard deviation σφS,liq . Subscript S for the solvent water. the reference mass densities ρ S,gas and thus very sensitive with respect to the precision of the method in general. Equation 11 was introduced as a selective criterion to sort out values similar to the reference density. If this criterion is applied, then the simulated results from 0.5 to 1.5 mol kg−1 are not part of the fitting process of the osmotic coefficient. The resulting values for 0.5 and 1.5 mol kg−1 provide water activities higher than 1.0 and negative osmotic coefficients. All results for φS,i,liq at higher concentrations are subsequently used for fitting the osmotic coefficients thanks to Equation 2 (Table 4). The osmotic coefficients for Dy(NO3 )3 aqueous solutions calculated from MD simulations are compared with the experimental ones. 11 Both data sets have been fitted with the parameters tabulated in Table 5 for two different dielectric constants of water, where r = 78.5 is the experimental value 78 at 298.15 K and r = 126.0 is for the rigid POL3 water model. 33,79 The variation of r influences the Debye-H¨ uckel parameters Ab and Bb derived from Equations 4 and 5, respectively, and requires an adjustment of the third fit parameter, the radius of dissolved ionic species rIon,aq . For a dielectric constant of 78.5, the fitted radius of the dissolved ionic species corresponds exactly to the hydrated cation radius of Dy3+ at infinite dilution. 80 The parameters Ab , 24 ACS Paragon Plus Environment

Page 24 of 38

Page 25 of 38

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 5: Fitting parameters for the Equations 2 and 3 at 298.15 K Ab a B∗b b rIon,aq c β1 d β2 e χ2f r (H2 O) = 78.5g sim. 1.17272 3.28231 ×109 4.75 0.49065 -0.00834 0.990 h 9 exp. 1.17272 3.28231 ×10 4.75 0.53846 -0.01096 0.999 r (H2 O) = 126.0i sim. 0.57620 2.59077 ×109 2.60 0.60492 -0.01836 0.989 h 9 exp. 0.57620 2.59077 ×10 2.60 0.65273 -0.02096 0.999 a 1/2 −1/2 b Debye-H¨ uckel parameter Ab in kg mol . Debye-H¨ uckel parameter B∗b in kg1/2 mol−1/2 m−1 . c Fitted ionic radius of the dissolved ion rIon,aq in ˚ A. d First-order fitting −1 e parameter β1 in kg mol . Second-order fitting parameter β2 in kg2 mol−2 . f Chi-square χ2 of the goodness of the fits. g Dielectric constant of water r (H2 O) at 298.15 from Fernández et al. 78 . h Experimental values from Reference. 11 i Dielectric constant of POL3 water r (H2 O) at 298.15 from Reference. 79

Bb , and rIon,aq determine the fit for dilute aqueous electrolyte solutions. The concentrated domain is described by the higher-order terms with the fit-parameters β1 and β2 , which are fitted towards experimental and simulated data. Additional higher-order fitting terms rather complicate the fit of the simulated data. The corresponding curves are depicted for the two dielectric constants, respectively (Figure 9). The osmotic coefficients simulated for salt concentrations lower than 1.5 mol kg−1 exhibit an offset from the experimental data, whereas all values at concentrations above 1.5 mol kg−1 show a good agreement with the data published by Rard and Spedding 11 (Figure 9). Due to the nature of Equation 2 and the identical parameters Ab , Bb and rIon,aq , the curves fitted for experimental and simulated data are very similar for low concentrations. At higher concentrations, the curves obtained from the simulated data provide slightly higher osmotic coefficients for both dielectric constants. Figure 10 compares the experimental data 11 with the simulated results based on Equation 3 for the activity coefficient of the electrolyte γE . Since the simulations provide higher osmotic coefficients compared to the experiments, the activity coefficients obtained from molecular dynamics are slightly higher for r = 78.5 and 126.0. Both sets of fit parameters allow expressing the activity coefficients in good agreement with the experimental data (Table 5). Nevertheless, using r = 78.5 allows reproducing experimental activity coefficients in a very



HaL Εr = 78.5

ΦS,aq HbEL

3 2 1 0 -1 -2 -3

ΦS,aq HbEL

0.0

3 2 1 0 -1 -2 -3

0.5

1.0

1.5

2.0

2.5

HbL Εr = 126

Exp. data MD data Fit HMDL 0.0

0.5

mol

1.0

bE

12

1.5 12

kg

2.0

2.5

-12

Figure 9: Experimental 11 and simulated (MD) osmotic coefficients compared with the DH fit of the the MD results for (a) r = 78.5 and (b) r = 126.0 (Table 5). Simulation results, which do not fulfill Equation 11 are colored gray.

ΓE HbEL

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 26 of 38

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0.0

Exp. data Εr = 78.5 Εr = 126

0.5

mol

1.0

bE

12

1.5 12

kg

2.0

2.5

-12

Figure 10: Experimental activity coefficients 11 γE (open circle) and results simulated using Equation 3 (solid lines) based on the fit parameters (Table 5) for the dielectric constants r 78.5 (blue) and 126 (red).

accurate way. Detailed tables for a comparison of osmotic coefficients, solvent activities and activity coefficients of Dy(NO3 )3 from both experiments and MD simulations are given in the Supporting Information (Tables S4 and S5).


Page 27 of 38

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


Conclusion Liquid-vapor equilibria of Dy(NO3 )3 aqueous solutions at different concentrations have been performed by means of molecular dynamics simulations using explicit polarization. Thermodynamic properties such as osmotic coefficients of water and activity coefficients of the salts dissolved have been calculated and are accessed in very good agreement with the experimental values 11 up to high concentrations with an extended form of the Debye-H¨ uckel equation. 18 The polarizable force field for water and the ions dissloved allow reproducing the vapor phase densities providing the activity of the solvent with respect to the vapor phase density of the pure solvent. In addition, the densities of the bulk liquid phases simulated at different concentrations are in very good agreement with the experiments. 77 Since the evaporation and recondensation processes of water molecules follow a Poisson process and a Maxwell-Boltzmann velocity distribution, respectively, it is of crucial interest understanding the origin of the noise occurring during the simulations for the amount of solvent molecules found in the gas phase in order to control the precision of this approach. Therefore, a statistic model has been introduced to estimate the error of this method and to understand the computational limitations of the model. Significantly increased precision comes along with a drastic increase of computational cost due to either larger liquid phases or to longer production runs. As a result of those constraints, the method provides results not as accurate as expected in the diluted domain (typically below 2.0 mol kg−1 ). However, activities and osmotic coefficients similar to the experiment are obtained when increasing the concentration of salts in the solution above 1.5 mol kg−1 . The evolving replacement of water molecules in the first Dy3+ hydration sphere by nitrate anions is observed at higher salt concentrations indicating a transition from hydrated towards more salt-like systems. However, the most structured density profiles with respect to the liquid-vapor interface are found for medium Dy(NO3 )3 concentrations of typically around 3.5 mol kg−1 . Up to a concentration of 3.5 mol kg−1 enough water is available to shield the ionic species leading to a certain pattern next to the interface. This pattern of peaks and minimum disappears at high concentrations due 27 ACS Paragon Plus Environment


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 28 of 38

to the lack of water for hydrating correctly the ion pairs and the transition to a structure mostly determined by ionic interactions. The method herein proposed, which is valid for concentrated solutions because of the inevitable noise in the dilute regime, is especially interesting because it is complementary to the McMillan Mayer approach, 28,29 which is valid only in the dilute domain. This fundamental description of an approach based on osmotic equilibrium between a vapor and a gas phase and the possibility to access thermodynamic properties of concentrated and even saturated aqueous electrolyte solutions will provide helpful tools for interpreting interfacial phenomena and for improving hydrometalurgical separation systems. This study paves the way for future investigations dealing with more complex and hardly accessible systems such as solvent mixtures, mixture of electrolyte solutions, and organic solvent phases as involved in the separation chemistry.

Supporting Information Available The equation of state of gaseous POL3 water at 298.15 K (Figure S1), a more detailed description of the statistics of evaporation and the herein introduced Poisson model, the radial distribution functions and coordination numbers for the simulation boxes (Figures S2 to S9), time averaged density profiles (Figures S10 to S20), finite size effects with respect to the thickness of the liquid phase (Tables S1 to S3, Figures S21 to 24), input settings for the AMBER MD simulation package, and a comparison of simulated and experimental results (Table S4 and S5) 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. 28 ACS Paragon Plus Environment

Page 29 of 38

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) Kunz, W.; Barthel, J. M. G.; Krienke, H. Physical Chemistry of Electrolyte Solutions: Modern Aspects (Topics in Physical Chemistry, Vol. 5, ed. by Deutsche Bunsengesellschaft); Steinkoppf: Darmstadt, 1998. (2) Stokes, R. H.; Robinson, R. A. Electrolyte Solutions - Second Revised Edition; Dover Publications: Mineola, New York, 2002; p 571. (3) Binnemans, K.; Tom, P.; Blanpain, B.; Gerven, T. V.; Yang, Y.; Walton, A.; Buchert, M. Recycling of Rare Earths : A Critical Review. J. Clean. Prod. 2013, 51, 1–22. (4) Massari, S.; Ruberti, M. Rare Earth Elements as Critical Raw Materials : Focus on International Markets and Future Strategies. Resour. Policy 2013, 38, 36–43. (5) Rydberg, J.; Cox, M.; Musikas, C.; Choppin, G. R. Solvent Extraction: Principles and Practice; Marcel Dekker, Inc.: New York, 2004. (6) Bonin, B.; Bouquin, B.; Dozol, M.; Lecomte, M.; Forestier, A. In Treatment and Recycling of spent Nuclear Fuel: Actinide Partitioning Application to Waste Management; Parisot, J.-F., Ed.; CEA Saclay and Groupe Moniteur: Paris, 2008. (7) Meridiano, Y.; Berthon, L.; Crozes, X.; Sorel, C.; Dannus, P.; Antonio, M. R.; Chiarizia, R.; Zemb, T. Aggregation in Organic Solutions of Malonamides: Consequences for Water Extraction. Solvent Extr. Ion Exch. 2009, 27, 607–637. (8) Dufrêche, J. F.; Zemb, T. Effect of Long-Range Interactions on Ion Equilibria in LiquidLiquid Extraction. Chem. Phys. Lett. 2015, 622, 45–49. (9) Bley, M.; Siboulet, B.; Karmakar, A.; Zemb, T.; Dufrêche, J.-F. A Predictive Model of Reverse Micelles solubilizing Water for Solvent Extraction. J. Colloid Interface Sci. 2016, 479, 106–114. 29 ACS Paragon Plus Environment


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

(10) Rard, J. A.; Shiers, L. E.; Heiser, D. J.; Spedding, F. H. Isopiestic Determination of the Activity Coefficients of Some Aqueous Rare Earth Electrolyte Solutions at 25 C. 3. The Rare Earth Nitrates. J. Chem. Eng. Data 1977, 22, 337–347. (11) Rard, J. A.; Spedding, F. H. Isopiestic Determination of the Activity Coefficients of Some Aqueous Rare-Earth Electrolyte Solutions at 25 ◦ C. 5. Dy(NO3 )3 , Ho(NO3 )3 , and Lu(NO3 )3 . J. Chem. Eng. Data 1981, 26, 391–395. (12) 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. (13) Robinson, R. A.; Lim, C. K. The Osmotic and Activity Coefficients of Uranyl Nitrate, Chloride, and Perchlorate at 25 C. J. Chem. Soc. 1951, 1840–1843. (14) Rush, R. M.; Johnson, J. S. lsopiestic Measurements of the Osmotic and Activity Coefficients for the Systems HCIO4 + OU2(ClO4)2 + H2O and NaClO4 + UO2(ClO4)2 + H2O at 25 ◦ C. J. Chem. Thermodyn. 1971, 3, 779–793. (15) Apelblat, A.; Korin, E. The Vapour Pressures of Saturated Aqueous Solutions of Sodium Chloride, Sodium Bromide, Sodium Nitrate, Sodium Nitrite, Potassium Iodate, and Rubidium Chloride at Temperatures from 227 K to 323 K. J. Chem. Thermodyn. 1998, 30, 459–471. (16) Debye, P.; H¨ uckel, E. Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. Phys. Z. 1923, 24, 185–206. (17) H¨ uckel, E. The Theory of Concentrated, Aqueous Solutions of Strong Electrolytes. Phys. Z. 1925, 26, 93–147. (18) 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. 30 ACS Paragon Plus Environment

Page 30 of 38

Page 31 of 38

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


(19) Bronsted, J. N. Studies on Solubility. IV. The Principle of the Specific Interaction of Ions. J. Am. Chem. Soc. 1922, 44, 877–898. (20) Guggenheim, E. A.; Turgeon, J. C. Specific Interaction of Ions. Trans. Faraday Soc 1955, 51, 747–761. (21) Pitzer, K. S. Thermodynamics of Electrolytes. I . Theoretical Basis and General Equation. J. Phys. Chem. 1973, 77, 268–277. (22) Pitzer, K. S.; Kim, J. J. Thermodynamics of Electrolytes. IV. Activity and Osmotic Coefficients for Mixed Electrolytes. J. Am. Chem. Soc. 1974, 96, 5701–5707. (23) Pitzer, K. S. In Act. coeffcients Electrolyte Solut., 2nd ed.; Pitzer, K., Ed.; CRC Press: Boca Raton, 1991; pp 75–153. (24) Waisman, E.; Lebowitz, J. L. Mean Spherical Model Integral Equation for Charged Hard Sphere. II. Results. J. Chem. Phys. 1972, 56, 3093–3099. (25) Waisman, E.; Lebowitz, J. L. Mean Spherical Model Integral Equation for Charged Hard Spheres I. Method of Solution. J. Chem. Phys. 1972, 56, 3086–3093. (26) Rasaiah, J. C.; Friedman, H. L. Integral Equation Methods in the Computation of Equilibrium Properties of Ionic Solutions. J. Chem. Phys. 1968, 48, 2742–2752. (27) Kunz, W.; Belloni, L.; Bernard, O.; Ninham, B. W. Osmotic Coefficients and Surface Tensions of Aqueous Electrolyte Solutions: Role of Dispersion Forces. J. Phys. Chem. B 2004, 108, 2398–2404. (28) McMillan, W. G.; Mayer, J. E. The Statistical Thermodynamics of Multicomponent Systems. J. Chem. Phys. 1945, 13, 276–305. (29) Simonin, J.-P. Study of Experimental-to-McMillan-Mayer Conversion of Thermodynamic Excess Functions. J. Chem. Soc., Faraday Trans. 1996, 92, 3519–3523. 31 ACS Paragon Plus Environment


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

(30) Lyubartsev, A. P.; Laaksonen, A. Osmotic and Activity Coefficients from Effective Potentials for Hydrated Ions. Phys. Rev. E. 1997, 55, 5689–5696. (31) 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. (32) 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. (33) 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. (34) 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. (35) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Solutions. I. J. Chem. Phys. 1951, 19, 774–777. (36) Hall, D. G. Kirkwood-Buff Theory of Solutions. Trans. Faraday Soc 1971, 67, 2516– 2524. (37) Ben-Naim, A. Inversion of Kirkwood-Buff Theory of Solutions - Application to WaterEthanol System. J. Chem. Phys. 1977, 67, 4884–4890. (38) Newman, K. E. Kirkwood-Buff Solution Theory: Derivation and Applications. Chem. Soc. Rev. 1994, 23, 31–40. (39) Gee, M. B.; Cox, N. R.; Jiao, Y.; Bentenitis, N.; Weerasinghe, S.; Smith, P. E. A Kirkwood-Buff derived Force Field for Aqueous Alkali Halides. J. Chem. Theory Comput. 2011, 7, 1369–1380. 32 ACS Paragon Plus Environment

Page 32 of 38

Page 33 of 38

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


(40) Weerasinghe, S.; Smith, P. E. A Kirkwood-Buff derived Force Field for Sodium Chloride in Water. J. Chem. Phys. 2003, 119, 11342–11349. (41) 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. (42) 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. 2015, 407, 76–83. (43) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813–826. (44) Panagiotopoulos, A. Z. Monte Carlo Methods for Phase Equilibria of Fluids. J. Phys. Condens. Matter 2000, 12, R25–R52. (45) L´ısal, M.; Smith, W. R.; Nezbeda, I. Accurate Computer Simulation of Phase Equilibrium for Complex Fluid Mixtures. Application to Binaries Involving Isobutene, Methanol, Methyl tert -Butyl Ether, and n -Butane. J. Phys. Chem. B 1999, 103, 10496–10505. (46) L´ısal, M.; Smith, W. R.; Nezbeda, I. Accurate Vapour-Liquid Equilibrium Calculations for Complex Systems using the Reaction Gibbs Ensemble Monte Carlo Simulation Method. Fluid Phase Equilib. 2001, 181, 127–146. (47) Luo, Y.; Roux, B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett 2010, 1, 183–189. (48) Taylor, R. S.; Dang, L. X.; Garrett, B. C. Molecular Dynamics Simulations of the Liquid/Vapor Interface of SPC/E Water. J. Phys. Chem. 1996, 100, 11720–11725.



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

(49) Jungwirth, P.; Tobias, D. J. Specific Ion Effects at the Air / Water Interface. Chem. Rev. 2006, 106, 1259–1281. (50) Wick, C. D.; Kuo, I. F. W.; Mundy, C. J.; Dang, L. X. The Effect of Polarizability for understanding the Molecular Structure of Aqueous Interfaces. J. Chem. Theory Comput. 2007, 3, 2002–2010. (51) 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. (52) Yuet, P. K.; Blankschtein, D. Molecular Dynamics Simulation Study of Water Surfaces: Comparison of Flexible Water Models. J. Phys. Chem. B 2010, 114, 13786–95. (53) Réal, F.; Vallet, V.; Flament, J. P.; Masella, M. Revisiting a Many-Body Model for Water based on a Single Polarizable Site: From Gas Phase Clusters to Liquid and Air/Liquid Water Systems. J. Chem. Phys. 2013, 139, 114502. (54) Holyst, R.; Litniewski, M. Evaporation into Vacuum: Mass Flux from Momentum Flux and the Hertz-Knudsen Relation revisited. J. Chem. Phys. 2009, 130, 074707. (55) Cheng, S.; Lechman, J. B.; Plimpton, S. J.; Grest, G. S. Evaporation of Lennard-Jones fluids. J. Chem. Phys. 2011, 134, 224704. (56) Holyst, R.; Litniewski, M.; Jakubczyk, D. Molecular Dynamics Test of the HertzKnudsen Equation for Evaporating Liquids. Soft Matter 2015, 11, 7201–7206. (57) Habenschuss, A.; Spedding, F. H. The Coordination (Hydration) of Rare Earth Ions in Aqueous Chloride Solutions from X-Ray Diffraction. I. TbCl3 ,DyCl3 , ErCl3 ,TmCl3 ,and LuCl3 . J. Chem. Phys. 1979, 70, 2797. (58) Helm, L.; Cossy, C.; Merbach, A. E. Oxygen- 17 Nuclear Magnetic Resonance Kinetic


Page 34 of 38

Page 35 of 38

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


Study of Water Exchange on the Lanthanide(III) Aqua Ions. Inorg. Chem. 1988, 27, 1973–1979. (59) Helm, L.; Foglia, F.; Kowall, T.; Merbach, A. E. Structure and Dynamics of Lanthanide Ions and Lanthanide Complexes in Solution. J. Phys. Condens. Matter 1994, 6, A137– A140. (60) 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. (61) Näslund, J.; Lindqvist-Reis, P.; Persson, I.; Sandström, M. Steric Effects Control the Structure of the Solvated Lanthanum(III) Ion in Aqueous, Dimethyl Sulfoxide, and N,N-Dimethylpropyleneurea Solution. An EXAFS and Large-Angle X-ray Scattering Study. Inorg. Chem. 2000, 39, 4006–4011. (62) Kowall, T.; Foglia, F.; Helm, L.; Merbach, A. E. Molecular Dynamics Simulation Study of Lanthanide Ions Ln3+ in Aqueous Solution including Water Polarization. Change in Cordination Number from 9 to 8 along the Series. J. Am. Chem. Soc. 1995, 117, 3790–3799. (63) 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. (64) Duvail, M.; Souaille, M.; Spezia, R.; Cartailler, T.; Vitorge, P. Pair Interaction Potentials with Explicit Polarization for Molecular Dynamics Simulations of La3+ in Bulk Water. J. Chem. Phys. 2007, 127, 034503. (65) Duvail, M.; Spezia, R.; Vitorge, P. A Dynamic Model to explain Hydration Behaviour along the Lanthanide Series. ChemPhysChem 2008, 9, 693–696.



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

(66) 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. (67) 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. (68) 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. (69) 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. (70) 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. (71) 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.; Goetz, A. W.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Paesani, F.; Roe, D. R.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker,; Wang, J.; Wolf, R. M.; Wu, X.; Kollman, P. A. Amber 14. University of California, San Francisco, 2014. (72) 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. (73) Caminiti, R.; Licheri, G.; Piccaluga, G.; Pinna, G. On NO− 3 –H2 O Interactions in Aqueous Solutions. J. Chem. Phys. 1978, 68, 1967–1970. 36 ACS Paragon Plus Environment

Page 36 of 38

Page 37 of 38

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


(74) Dang, L. X.; Chang, T.-M.; Roeselova, M.; Garrett, B. C.; Tobias, D. J. On NO− 3 –H2 O Interactions in Aqueous Solutions and at Interfaces. J. Chem. Phys. 2006, 124, 066101. (75) 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. (76) 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. (77) Spedding, F. H.; Shiers, L. E.; Brown, M. A.; Baker, J. L.; Guitierrez, L.; McDowell, L. S.; Habenschuss, A. Densities and Apparent Molal Volumes of Some Aqueous Rare Earth Solutions at 25 C. III. Rare Earth Nitrates. J. Phys. Chem. 1975, 79, 1087–1096. (78) Fernández, D. P.; Mulev, Y.; Goodwin, A. R. H.; Levelt Sengers, J. M. H. A Database for the Static Dielectric Constant of Water and Steam. J. Phys. Chem. Ref. Data 1995, 24, 33–70. (79) Wang, J.; Cieplak, P.; Cai, Q.; Hsieh, M.-J.; Wang, J.; Duan, Y.; Luo, R. Development of Polarizable Models for Molecular Mechanical Calculations. 3. Polarizable Water Models Conforming to Thole Polarization Screening Schemes. J. Phys. Chem. B 2012, 116, 7999–8008. (80) Ruas, A.; Guilbaud, P.; Auwer, C. D.; Moulin, C.; Simonin, J.-P.; Turq, P.; Moisy, P. Experimental and Molecular Dynamics Studies of Dysprosium(III) Salt Solutions for a Better Representation of the Microscopic Features Used within the Binding Mean Spherical Approximation Theory. J. Phys. Chem. A 2006, 110, 11770–11779.



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

Graphical TOC Entry

Graphical TOC Entry


Page 38 of 38

Simulating Osmotic Equilibria: A New Tool for ... - ACS Publications

Recommend Documents