Atomistic Description of Binary Lanthanoid Salt Solutions: A Coarse

ARTICLE pubs.acs.org/JPCB

Atomistic Description of Binary Lanthanoid Salt Solutions: A Coarse-Graining Approach John Jairo Molina,*,†,§ Magali Duvail,‡ Jean-Franc-ois Dufrêche,*,§ and Philippe Guilbaud*,‡ †

Physicochimie des Electrolytes, Colloides et Sciences Analytiques, UMR 7195, UPMC Universite Paris 06, F-75005 Paris, France RadioChemistry and Processes Department, Nuclear Energy Division, Commissariat a l’Energie Atomique, F-30207 Bagnols sur Ceze, France § Institut de Chimie Separative de Marcoule, UMR 5257, CEA-CNRS, Universite Montpellier 2, Site de Marcoule Bâtiment 426, BP 17171, 30207 Bagnols-sur-Ceze Cedex, France ‡

bS Supporting Information ABSTRACT: The experimental difficulties inherent to the solution chemistry of actinoids and lanthanoids have led to the use of a wide variety of models, from the microscopic to the macroscopic scale, in an attempt to represent their solution properties. Molecular dynamics (MD) simulations, with explicit solvents, have been successfully used to describe the structural characteristics, but the limits on the accessible length and time scales do not allow for an equivalent description of the macroscopic properties. In this study, we propose a multiscale approach, based on MD simulation results, to study the thermodynamic and structural properties of a series of lanthanoidchloride aqueous solutions. An inversion procedure, based on the approximate hypernetted chain (HNC) closure and the StillingerLovett sum rules for ionic liquids, is used to obtain the effective ionion potentials from MD-generated radial distribution functions (RDF). Implicit solvent Monte Carlo (MC) simulations are then performed to compute the osmotic coefficients of the salt solutions. This coarse-grained strategy provides accurate effective pair potentials for the lanthanoid salts, derived from an atomic model. The method presented here is an attempt to bridge the gap between MD and the thermodynamic properties of solutions that are experimentally measured.

’ INTRODUCTION An accurate description of the fundamental structural and thermodynamic properties of actinoid (An) and lanthanoid (Ln) salt solutions is indispensable for many important industrial applications. The reprocessing of spent nuclear fuel provides such an example, in which case the actinoid (An3þ) and lanthanoid (Ln3þ) activity coefficients must be known up to high concentrations, in order to properly control the liquid liquid extraction procedure used in the separation process. However, this data is relatively rare, especially for actinoid salt solutions at high concentrations, due to the complex chemical nature of the solutions (e.g., the different number of oxidation states and structural forms) and the experimental difficulty of working with radioactive elements. Therefore, any predictive theory capable of providing a suitable physical interpretation to the available experimental data is of great interest. MD simulations provide a convenient method for studying the physical and chemical properties of An3þ 13 and Ln3þ 49 aqueous salt solutions; but the computational limits on the r 2011 American Chemical Society

system sizes that can be simulated do not allow for a direct measurement of all thermodynamic properties. Spatial correlations (i.e., radial distribution functions) are easily obtained, but free energies and osmotic coefficients are prohibitively costly to measure. To overcome these difficulties, it is usual to employ simpler, implicit solvent, models. Due to the high charge of the cation (þ3), these models must explicitly allow for ion association. As such, the primitive hard-sphere representations commonly used for simple 11 electrolytes are inadequate. Among the theories that have been successfully used to study these systems, the binding mean spherical approximation10,11 (BiMSA) is perhaps the most popular. Derived using the Wertheim formalism, BiMSA12,13 provides an explicit mathematical solution for the equilibrium thermodynamic properties of an arbitrary mixture of hard spheres, with Received: June 21, 2010 Revised: February 9, 2011 Published: March 28, 2011 4329

dx.doi.org/10.1021/jp1110168 | J. Phys. Chem. B 2011, 115, 4329–4340

The Journal of Physical Chemistry B an arbitrary number of bonds, in a dielectric medium. Previous studies have shown that this analytical theory is capable of reproducing the osmotic coefficients for Ln3þ salts with a restricted number of adjustable (fitting) parameters, some of which can be related to microscopic properties of the system.7,14,15 In order to test the predictive power of such a theory, it is necessary to establish the physical (microscopic) meaning of the undetermined free parameters as well as provide a systematic method for determining them nonempirically (i.e., without fitting the model). Of the available parameters, the hydrated cation diameter at infinite dilution (σ(0) þ ) would seem to be a perfect candidate, since it can be directly related to experimental data and simulation results. However, the relation between the two, as well as the precise value that should be used, is not easily determined. The σ(0) þ diameter corresponds to a hydrated cation whose size exceeds the first solvation shell, making it unreachable with experiments; but including the first 7 two solvation shells leads to a σ(0) þ parameter that is too large. An external cutoff criteria (energetic, geometric, or dynamic) is required in order to determine which water molecules are to be included within this second hydration sphere. Additionally, the evolution of the σ(0) þ diameter along the lanthanoid(III) series is still not firmly established. David et al.1619 have assumed, based on previous calculations on the alkali metals, that this hydrated cation diameter should increase as its ionic radius decreases. However, for the lanthanoid(III) and actinoid(III) series, this hypothesis is questionable, since the high charge on the cation will induce a strong ordering of the water molecules in the first two coordination spheres along the series,8,20 in contrast with the case of alkali metals. The purpose of this work is to develop an alternative numerical solution, which uses the “exact” results provided by MD simulations, to construct an implicit solvent description that can be used to efficiently compute the thermodynamic properties of interest. The resulting coarse-grained model removes the unwanted solvent degrees of freedom, treating the ions as charged solvated particles moving in a dielectric medium. This representation allows for a much simpler calculation of the ion thermodynamic properties and could eventually be used to unambiguously determine the free parameters of an analytic model such as BiMSA, as has been done for the case of 11 electrolytes.21 We will focus specifically on lanthanoid solutions, since these systems are the most likely candidates for predictive applications in actinoid chemistry. In fact, Ln3þ are often regarded as An3þ salt analogues, and their structural20,2231 and thermodynamic3238 properties have been the subject of extensive studies.

’ METHODS Microscopic Model. MD Simulations. We have performed MD simulations of five representative lanthanoid chloride Ln3þC1aqueous solutions (Ln = La, Nd, Eu, Dy, Lu) at two different concentrations, c = 0.5 and 1.0 mol kg1, using explicit polarization in the NPT ensemble. The simulations were carried out with SANDER 10, a module of AMBER 10.39 Long-range interactions were calculated using the particle-mesh Ewald method.40 The temperature and pressure were controlled using a Berendsen41 thermostat and barostat, with a target temperature and pressure of T = 300 K and P = 1 atm, and a coupling constant of 1 ps. The equations of motion were numerically integrated using a 1 fs time step. All the systems were equilibrated during at least 100 ps, with the lanthanoid salts initialized as first shell

ARTICLE

Table 1. Ln3þ Parameters Used for the MD Simulationsa Ln3þ

εii

σii

Rb

La3þ

0.272

3.564

1.134

Nd3þ

0.297

3.430

0.955

Eu3þ

0.322

3.271

0.823

Dy3þ

0.331

3.118

0.728

Lu3þ

0.356

2.887

0.623

Energies are in kJ mol1, distances in Å, and atomic polarizabilities in Å3. b Atomic polarizabilities have been calculated by Clavaguera et al.45 a

dissociated,8 and production runs of 2 ns were subsequently performed. The van der Waals energy was described by a 126 Lennard-Jones (LJ) potential. The LJ parameters for the Dy3þ and Nd3þ cations were taken from ref 8. The remaining Ln3þ LJ parameters (Table 1) were adjusted using the same method, in order to reproduce the experimental values for the structural properties of hydration24,25 (coordination number and firstcoordination sphere distance), with the POL3 polarizable water model. The polarizable model defined by Smith et al.42 was used for the Cl ion, and the POL343,44 model was used to describe the water molecules. The atomic polarizabilities for the Ln3þ ions are those given by Clavaguera et al.45 The LJ parameters for the Ln3þCl interactions were obtained using standard LorentzBerthelot mixing rules. For systems at 0.5 (1.0) mol kg1, we have used a periodically replicated cubic simulation box containing 887 (1498) water molecules along with 8 (27) Ln3þ and 24 (81) Cl ions. The average lateral size of the simulation box was 31 and 37 Å, for c = 0.5 and c = 1.0 mol kg1, respectively. As is customary when studying the equilibrium properties of liquids, we have focused our attention on the RDF (g(r)) computed from the MD simulations. These spatial correlation functions give the probability of observing a particle at a distance r from another particle (located at the origin), irrespective of the positions of all other particles. It is not surprising then that they can be used to characterize the particle interactions. In fact, there is a uniqueness theorem which states that for a measured g(r) there exists a unique effective pair potential46 (i.e., a system of particles interacting with this pair potential, at the same density, will reproduce the g(r)). At infinite dilution, the interaction potential for a pair of particles is directly given by the g(r) UðrÞ

¼

FN s f0

kB T ln gðrÞ

ð1Þ

where FN = N/V, with N the total number of particles. PMF Calculations. At finite concentrations, a relation analogous to eq 1, serves as the definition of the so-called potential of mean force (PMF) Ψ(r) ΨðrÞ ¼ kB T ln gðrÞ

ð2Þ

which, as its name suggests, gives the force felt by two particle at a distance r from each other, averaged over the positions of all the other particles. If we could obtain accurate estimates for g(r) at infinite dilution, the corresponding PMF would give us the effective ionion potential averaged over the water configurations. Since these potentials characterize the average solvation effects on the ions, they can be used to perform implicit solvent simulations to measure thermodynamic quantities not directly accessible from explicit solvent simulations. Performing such calculations directly is not recommended, but there exist several ways of reducing the computational complexity involved. 4330

dx.doi.org/10.1021/jp1110168 |J. Phys. Chem. B 2011, 115, 4329–4340

The Journal of Physical Chemistry B

ARTICLE

We have used two different inversion schemes to obtain the solvent averaged ionion potentials: (1) a PMF calculation at infinite dilution, using an umbrella sampling technique, which is “exact”; and (2) a hypernetted chain (HNC) integral equation approximation using g(r)’s computed at finite concentrations. The reduced statistical uncertainties obtained from finite concentration simulations ultimately allow for a more accurate estimate of the interaction potentials. This proves to be important at large distances, where we assume the potential reduces to a Coulombic interaction in order to set the zero of our potential energy. Strictly speaking, the HNC potentials obtained from finite concentration calculations are density dependent and do not represent the “real” pair potential; however, a comparison between both sets of potentials (PMF vs HNC) shows that the density dependence is not significant (at least at c = 0.5 mol kg1). The PMFs for Ln3þCl, Ln3þLn3þ, and ClCl in water at infinite dilution were computed using additional umbrella sampling simulations,47 with the weighted histogram analysis method (WHAM).4850 This procedure computes the PMFs as a free energy difference, giving the effective pair potentials up to an additive constant. The interparticle distance (the reaction coordinate) was varied from 2.0 to 12.0 Å, in steps of 0.5 Å; it was constrained using harmonic biasing potentials of 2, 10, and 5 kcal mol1 Å2 for Ln3þCl, Ln3þLn3þ, and ClCl, respectively. In the case of Ln3þCl, it was found that additional calculations (with a stronger bias) were required to accurately represent the potential barrier separating the contact ion pair and the solvent-separated ion pair. Specifically, simulations for interparticle distances of 3.25, 3.5, and 3.75 Å for Nd3þCl and 3 and 3.25 Å for Dy3þCl were performed, with a force constant of 10 kcal mol1 Å2. All the umbrella-sampling simulation protocols were optimized to ensure a good overlap of the constraint windows, and therefore a good representation of the reaction pathways. These simulations were carried out using SANDER1039 in the NVT ensemble at 300 K. The interaction potentials were the same ones used in the unconstrained MD simulations. All the systems (or windows) were previously equilibrated at room temperature and pressure for at least 50 ps, in order to obtain the correct density, and production runs were subsequently run for 1 ns. These constrained simulations were performed for a single ion pair (Ln3þLn3þ, Ln3þCl, and ClCl), in a periodically replicated cubic simulation box containing approximately 1000 water molecules. A lateral box size length of 31 Å was used, corresponding to the average size obtained from the NPT equilibration. Continuous Solvent Model. Inverting the Effective Pair Potential. Using the PMF procedure detailed above, we were able to obtain the effective (average) solutesolute pair potentials up to an additive constant. In principle, this is not a problem, since we expect the effective interactions to be Coulombic at large distances: the limiting value of the effective potential should be that of two charged particles in a dielectric medium. However, the poor statistical accuracy obtained for Ψ(r) at long distances, together with artifacts arising from the use of periodic boundary conditions, make it impossible to fix the reference potential energy unambiguously. To overcome this issue, we use the approximate HNC closure51 eff βUνμ ðrÞ ¼ hνμ ðrÞ ceff νμ ðrÞ log½hνμ ðrÞ þ 1

to compute the effective potential Ueff νμ(r) from the corresponding hνμ(r) and ceff νμ(r) obtained from the MD simulations at c = 0.5 mol kg1; in eq 3, hνμ(r) = gνμ(r) 1 is the pair correlation function and ceff νμ(r) the effective direct correlation function. Only hνμ(r) = gνμ(r) 1 is directly measured during the simulations, since ceff νμ(r) can be obtained from hνμ(r) by means of the OrnsteinZernike (OZ) integral equation51 Z 0 0 dr0 ceff hνμ ðrÞ ¼ ceff ðrÞ þ ð4aÞ νμ νγ ðr ÞFγ hγμ ðjr r jÞ

∑γ

which in Fourier space reduces to hνμ ðkÞ ¼ ceff νμ ðkÞ þ

∑γ ceffνγðkÞFγ hγμ ðkÞ

ð4bÞ

where Fγ = Nγ/V is the number density of species γ. If hνμ(k) is known for all k, cνμ(k) follows directly from eq 4b; but computation of hνμ(k) requires knowledge of hνμ(r) up to values of r where hνμ(r) has converged to its asymptotic value of 0. Since we can only sample values of r up to L/2 (L is the length of the cubic simulation box), it is necessary to extend the measured hνμ(r) to distances r > L/2 in a manner consistent with the long-range Coulombic nature of the solutesolute interactions. We follow an inversion procedure similar to that presented in ref 52, adapted to two-component systems with long-range Coulomb interactions. For this purpose, we use the well-known StillingerLovett sum rules for ionic fluids51 Z dr xμ zμ gνμ ðrÞ ¼ zν FN ð5aÞ

∑μ

FN

∑ν

xν zν

∑μ

Z r 2 dr xμ zμ gνμ ðrÞ ¼ 6ΛD 2

∑ν xν zν2

ð5bÞ

where xν = Fν/FN is the density fraction, zν = qν/e the valence (qν is the charge and e the elementary unit of charge), and ΛD the Debye screening length ΛD 2 ¼

βe2 FN ε0 εr

∑ν xν zν 2

ð6Þ

We assume that the dielectric constant of the solution εr is given by the bulk value of the water model used in the MD simulations (106 for POL3). A concentration-dependent dielectric constant could be used to quantitatively correct for the many-body contributions to the effective pair potential,53 but this has not been considered here. Since we are interested in the long-range behavior of hνμ(r), it is again more comfortable to work in Fourier space and consider the short k, long-wavelength, limit of hνμ(k) instead. A straightforward calculation shows that eq 5a and 5b can be expressed as FN FN

∑μ xμ zμ hνμ ðk ¼ 0Þ ¼

zν

∑ν xν zν ∑μ xμ zμ ðrk 2hνμðkÞÞjk ¼ 0 ¼ 6ΛD 2 ∑ν xν zν2

ð7aÞ ð7bÞ

If hνμ(k) is a smooth, well-behaved function, it allows a Taylor expansion in even powers of k at small k ð0Þ ð2Þ 2 ð4Þ 4 ð6Þ 6 ð8Þ 8 10 hmod νμ ðkÞ ¼ aνμ þ aνμ k þ aνμ k þ aνμ k þ aνμ k þ O ðk Þ

ð3Þ

ð8Þ 4331


The Journal of Physical Chemistry B where ðnÞ ¼ 4π aνμ

ARTICLE

Z ¥ ð1Þn dr hνμ ðrÞr 2n þ 1 ð2n þ 1Þ! 0

and it is assumed that all a(n) νμ (r) are finite (i.e., hνμ(r) decay sufficiently rapidly at large r). The following extrapolation procedure was used: we fitted the small k values of hνμ(k) to the model form given in eq 8, subject to the constraints imposed by the StillingerLovett conditions, eqs 7a and 7b. To avoid discontinuities when computing the Fourier transforms, we truncated the MD-generated hmd νμ (r) at the node rc closest to L/2, such that hmd νμ (rc) = 0, and defined the following pair correlation functions 8 < hmod ðkÞ k e k∼ νμ ð9Þ hνμ ðkÞ ¼ ∼ : hdat k>k νμ ðkÞ where

( hdat νμ ðrÞ

¼

hmd νμ ðrÞ 0

r < rc r g rc

ð10Þ

dat ~ and hmod νμ (k) is obtained by fitting hνμ (k), with k chosen to ensure continuity. The best results were obtained by fitting the functions in the following order: h12(k), h22(k) ,h11(k), where 1 and 2 refer to the cation and anion, respectively. Since the zeroth-order (0) coefficient a(0) 12 for the fit of h12(k) uniquely determines a11 (k) (k) (eq 7a), the second fit has only four free parameters and a(0) 22 (2) (2) ({a(2i) 22 ,i=1,4}). Finally, the second-order coefficients a12 and a22 (2) are used to determine a11 (eq 7b), so that the fit of h11(k) only requires three free parameters ({a(2i) 11 ,i=2,4}). The advantage of this procedure can be fully appreciated by inspecting the BhatiaThornton (BT) structure factors51,54

∑ν ∑μ SνμðkÞ

ð11aÞ

∑ν ∑μ zμ SνμðkÞ

ð11bÞ

∑ν ∑μ zν zμ SνμðkÞ

ð11cÞ

SNN ðkÞ ¼ SNZ ðkÞ ¼ SZZ ðkÞ ¼

where Sνμ(k) = xνδνμ þ FNxνxμh(k) is the partial structure factor. The BT structure factors, which give information on the density and concentration correlations within the fluid, possess well-defined long-wavelength limits, which are independent of the specific short-range potentials (i.e., they are a consequence of the long-range Coulomb interactions). Specifically, at small k, we have SNN(k) ∼ k0, SNZ(k) ∼ k2, and SZZ(k) ∼ k2, with SZZ(k) > 0 for all values of k. In Figure 3 (solid lines) we give the partial and BT structure factors obtained directly from hmd νμ (r) for a solution of Dy3þCl at 0.5 mol kg1. It is immediately clear from the BT structure factors that the long-range correlations do not correspond to those of an ionic fluid, since neither SNZ nor SZZ goes to zero at small k. A comparison between the structure dat factors corresponding to hmd νμ (r), hνμ (r), and hνμ(r), the original, truncated, and fitted pair correlation functions, respectively, is given in Figure 3 for the case of Dy3þCl at 0.5 mol kg1. We see that the difference between the original and the truncated structure factors is only appreciable for k j 0.5 Å1, and between truncated and fitted for k j 0.25 Å1. The latter, though

Figure 1. Radial distribution functions calculated from molecular dynamics simulations of Ln3þCl salt solutions, at 0.5 mol kg1 (top) and 1.0 mol kg1(bottom).

considerable at small k (especially for like ion interactions), is barely noticeably when comparing the functions in real space (not shown). This is to be expected, since the fitting procedure is just designed to extrapolate the functions at long-range, in such a way as to satisfy the exact small k limits. In a previous work,55 we obtained the short-range ionion potential Usr(r) by directly fitting the asymptotic limit of Ueff(r) to a Coulombic potential, but this has the disadvantage of not providing a consistent representation of the electrostatic interactions (e.g., a different dielectric constant was obtained for the three pair potentials). Implicit Solvent Monte Carlo Simulations. The effective ionion pair potentials Ueff νμ(r) obtained through eq 3, from the full-atom RDFs measured during the MD simulations at c = 0.5 mol kg1, were used in implicit solvent NVT Monte Carlo (MC) simulations to study the solute thermodynamic and structural properties. Specifically, we have focused on the osmotic coefficient Φ and the radial distribution function g(r). The potentials Ueff νμ(r) are naturally separated into a long-range Coulombic interaction and a short-range solvent averaged interaction Usrνμ(r) zν zμ e2 ð12Þ 4πε0 εr r The long-range Coulomb interactions were computed using the Ewald summation technique, while the short-range interactions sr eff ðrÞ ¼ Uνμ ðrÞ Uνμ

4332



ARTICLE

Figure 2. Potential of mean force profiles for Ln3þLn3þ, Ln3þCl, and ClCl obtained from umbrella sampling simulations for Ln3þ = Nd3þ and Dy3þ.

are computed using a simple truncation scheme (rc = 13 Å), and standard periodic boundary conditions were used. Simulations were performed over a concentration range from 0.01 to 1.44 mol L1, with a dielectric constant εr = 106 corresponding to the POL3 water model used in the original MD simulations. Additional simulations were performed with εr = 78.4 in order to evaluate the water model dependence. All systems were composed of 54 salt molecules, for 216 total particles. Equilibration runs of 5 104 MC steps were followed by production runs of 104 blocks, of 102 MC steps each, to compute the desired block averages. During every MC step, an attempt is made to perform a random single-particle displacement for all particles in the system, attempted moves are accepted/rejected using the Metropolis scheme.56 Osmotic coefficients were calculated using the pressure energy relation for ionic systems57 Φ¼

βP βPsr ÆU lr æ ¼ 1þ þ FN 3N FN

ð13Þ

where Æ 3 3 3 æ denotes an ensemble average. The short-range contribution to the pressure Psr is computed using the virial formula56 βPsr 1 DU sr Φsr ¼ ¼ r ð14Þ i 3N i Dri 3 FN

∑

while the long-range contribution is given directly by the average electrostatic energy ÆUlræ. It should be noted that Psr refers to the pressure contributions arising from the non-Coulombic interactions only; they should not be confused with the real-space (shortrange) part of the electrostatic potential used in the Ewald summation. Convergence for the osmotic coefficients was obtained after roughly 7 105 MC steps, which corresponds to a few days of computing time on a single desktop CPU. To achieve the same level of accuracy from atomic simulations would have required weeks of computing time.

’ RESULTS AND DISCUSSION Binary Solutions of Ln3þ Chloride Salts. We present the MD

results for the five lanthanoid chloride salt solutions at 0.5 (1.0) mol kg1. As mentioned in our previous experimental and

Figure 3. Partial (top) and BhatiaThornton (bottom) structure factors obtained during the extrapolation procedure for the case of Dy3þCl at 0.5. mol kg1 (solid) original RDF, (dash) truncated RDF, (circle) small k fit (subject to the StillingerLovett conditions) of the hμν(k).

computational studies,7,8 no inner-sphere complexes are observed between Ln3þ cations and Cl anions at 0.5 mol kg1. In the Ln3þCl radial distribution functions, the distance of the first peak corresponds to the second coordination shell (SCS) (Figure 1). A decrease in this Ln3þCl distance was observed through the lanthanoid series, from about 5.14 Å for La3þ to 4.82 Å for Lu3þ; which follows a corresponding decrease in the Ln3þOH2 distance (Supporting Information, Figures S1 and S2).5860 This decrease in the Ln3þCl distance occurs alongside a decrease in the average number of anions within the Ln3þ SCS, from 1.15 to 0.94 Å for La3þ and Lu3þ, respectively (Supporting Information, Table S1). As already described in a previous study, two distinct positions for the Cl anion in the Ln3þ SCS are observed, corresponding to two different H-bonds between anions and water molecules located in the Ln3þ first coordination shell.8 The same observations (decrease of the Ln3þCl distance, accompanied by a decrease in the number of Cl anions within the Ln3þ SCS) are valid at 1 mol kg1, except for La3þ, for which some inner-sphere association was observed (Figure 1). However, this first-shell association 4333



Figure 4. Effective short-range HNC pair potentials derived from RDFs at 0.5 mol kg1; the average potentials for each of the pair interactions are also shown.

between Cl and La3þ (average distance dLaCl = 2.73 Å) is weak, with only one La3þCl inner-sphere pair observed during 300 ps, despite the 81 Cl anions present in the solution. No characteristic structure is observed in the Ln3þLn3þ RDFs, whatever the cation or the concentration (Figure 1). The strong electrostatic repulsion between cations favors configurations in

ARTICLE

which they are located as far as possible from each other. On the other hand, the ClCl RDFs show significant anionanion structuring in the solution, whatever the concentration or the Ln3þ cation (Figure 1). A well-defined first peak in the RDF is observed around 5 Å in all cases, which corresponds to an outersphere association in the Ln3þ second coordination shell. It is clear that the results depend on the atomic potentials used and that different parametrizations will give rise to different values. In our case, we have parametrized the LJ potentials using the variation of the ionic radii along the series,61 assuming a coordination number of 9 for La3þ and Nd3þ, and 8 for Eu3þ, Dy3þ, and Lu3þ, as explained in a previous study.62 Since there is still no clear consensus on the coordination number for the intermediate elements in the series, it is necessary to adopt such a convention. These “simple” potentials are enough to reproduce the experimental values for first hydration sphere structure, namely the Ln3þO distance and the coordination number (see Supporting Information). Although more sophisticated potentials exist, such as those recently developed by Beuchat et al.63 using quantum chemical calculations, they are not necessarily more accurate for studying macrosopic thermodynamic properties (such as the activity coefficients), since they are generally parametrized to reproduce hydration properties on small water clusters. The computation of the thermodynamic properties of the solution is more time consuming due to the long-range forces which require, in our case, an implicit solvent description based on the McMillanMayer theory. It is difficult to predict how a change in the atomic potentials will affect the final results, unless the corresponding effective potentials (obtained from MD simulations using the different atomic force fields) are compared with each other. The potentials we use have two main advantages: (i) they are simple enough to allow for the calculation of the thermodynamic properties at a reasonable time cost (typically several weeks for one electrolyte at a given concentration) and (ii) they are consistent with available experimental data, as described below. Effective Potentials. The infinite dilution PMF’s obtained from umbrella sampling MD simulations for Nd3þNd3þ, Dy3þDy3þ, and ClCl are shown in Figure 2. The potentials have all been shifted in such a way that Ψvμ(r=12Å) = 0. The cationcation potentials both present the same monotonic decrease, although the Nd3þNd3þ interaction is clearly more repulsive. The cationanion potentials exhibit the same overall behavior, with the maxima (minima) of the Nd3þCl interactions shifted to larger r values, which also indicates a more repulsive potential. Also worth mentioning is the large potential barrier J15kBT separating the ion pair and the solvent-separated ion pair (the outer-sphere pair), which explains why no innersphere complexes were observed. The HNC potentials obtained from the MD-generated g(r) are given in Figure 4. The cationcation potentials present a hard-core repulsion which reproduces the expected ordering in the series (La < Nd < Eu < Dy < Lu). The long-range, small energy, fluctuations (which can become attractive for Eu and Lu) are a consequence of the uncertainties present in the original RDFs. The cationanion potentials also reproduce the expected ordering, in both the height (depth) and position of the first maxima (minima), with the exception of the Nd3þCl potential which seems to completely underestimate the ion association (this can be already be seen in the MD g(r)). The Cl Cl potential is essentially the same for all the solutions considered, and the average of the five estimates was the potential used in the 4334



ARTICLE

Figure 5. Comparison between the infinite dilutions PMF, the HNC effective potential, and the asymptotic Coulomb potential for Nd3þCl and Dy3þCl.

Figure 7. Variation of the cationanion RDFs obtained from implicit solvent MC simulations as a function of salt concentration for solutions of Dy3þCl.

Figure 6. Bjerrum association constants obtained from the effective cationanion HNC potentials (eq 15). The open symbols show the values of K0(r) as a function of the cutoff distance r, and the solid lines eff show the value of the integrand (4π2 exp[βULn3þC1 ]) as a function of r.

Table 2. Association Constants K0 Obtained from the CationAnion Effective HNC Potentials and from BiMSA Fits to the Experimental Osmotic Coefficients (in L mol1)

a

Ln3þ

K0 (HNC)

K0 (BIMSA)a

La3þ Nd3þ

6.1 4.2

3.05 2.59

Eu3þ

5.6

2.55

Dy3þ

4.7

2.39

Lu3þ

4.2

2.08

Association constants obtained by Ruas et al.15

MC simulations. The average Ln3þCl and Ln3þLn3þ potentials are also given; these have been used in separate calculations in order to evaluate the sensitivity of the measured thermodynamic data to variations in the potentials.

Figure 8. Ln3þCl RDF obtained from MC simulations at two different concentrations for the series of lanthanoids.

A comparison between the effective HNC potentials and the exact PMF calculations is given in Figure 5, with the PMFs shifted to align with the first peak in the HNC potential. Very good agreement is obtained, especially at short distances, which shows that the concentration dependence of the effective potential is still very weak at c = 0.5 mol kg1. The poor asymptotic behavior of the PMFs, with respect to the Coulomb potential, is 4335



ARTICLE

Figure 9. Comparison between RDFs obtained from implicit solvent (MC) and explicit solvent (MD) simulations for solutions of Dy3þCl at 0.5 mol kg1 (left) and 1.0 mol kg1 (right).

due to the uncertainty in the WHAM calculations and to artifacts arising from the use of periodic boundary conditions. The former is due to the limited sampling that can be achieved for large ionion separations, while the latter results from the interaction with ions from all the periodic replicas of the simulation cell. The advantage of the approximate HNC inversion is clear, since it has allowed us to obtain the correct long-range behavior. The main drawback to this procedure is the fact that all information of correlations beyond the cutoff rc is lost. We expect that this effect will only be important for the cationcation interactions, since the cationanion and anionanion RDFs have already converged to their asymptotic value around rc=10 Å. A possible way to correct for this would be to use independent measurements of hνμ(k) and hνμ(r) in the extrapolation procedure.52 Additionally, the cationanion effective potentials can be used to quantify the amount of ion pairing in the different solutions, by

means of the Bjerrum association constant K0, where K0 is defined as the configuration integral for a pair of “associated” ions Z d eff K0 ðdÞ ¼ dr 4πr 2 exp½βULn ð15Þ 3þ Cl ðrÞ 0

where d is a characteristic length, which distinguishes between free and paired ions. There is no unique way of defining d, and we have chosen to use the value d = d* corresponding to the inflection point in K0(d). Figure 6 shows K0(d), as well as the Boltzmann factors appearing in the integral, as a function of the cutoff distance d. Our estimated values for K0 = K0(d*) are given in Table 2, along with those obtained by Ruas et al.15 after fitting the experimental osmotic coefficients to a BiMSA theory, in which the cation diameters, the association constant, and the (concentration dependent) dielectric permittivity were used as 4336


The Journal of Physical Chemistry B fitting parameters. Although we have used a vastly different approach, we are able, without any fitting parameters, to reproduce the same variation across the series of lanthanoids, obtaining values for K0 which differ only by a factor of 2. Implicit Solvent Solutions. We present the results obtained from implicit solvent MC simulations using the effective HNC ionion potentials described above. Figure 7 shows the variation of the cationcation RDF as a function of the concentration for the case of Dy3þCl. As a reference, we have also included the infinite dilution limit for the g(r) (eq 1). It is worth mentioning that, even at 0.01 mol L1, the solution is far from what can be considered as dilute. As the concentration is increased, the peaks in the g(r) are lowered, until a seemingly stable state is reached at c = 0.5 mol L1. The same concentration dependence is observed for all the other lanthanoid cations considered here. The variation of the RDFs as a function of the cation is the same as that observed in the original MD-generated RDFs, and is shown in Figure 8 for two concentrations, c = 0.25 and 1.06 mol L1. A comparison between MD- and MC-generated RDFs is given in Figure 9 for the case of Dy3þCl. At 0.5 mol kg1 (=0.51 mol L1) the agreement is seen to be excellent for the cationanion and anionanion RDFs. The cationcation comparison is less favorable, specifically at long distances, but this was to be expected, since we have truncated the cationcation RDFs at a node rc ∼ 11 Å in order to perform the inversion. We obtain the same results at c = 1.0 mol kg1 (=1.06 mol L1), although there is a slight underestimation of the cationanion pairing. This difference is most noticeable in the second maxima of the first peak. These results validate the coarse-graining strategy, although they seem to indicate that n-body effects start to become important at molar concentrations. The osmotic coefficients obtained from the MC simulations are compared to the experimental values32 in Figure 10 and Table 3, with the appropriate LewisRandall to McMillanMayer reference scale conversion.64 It is evident that the results exhibit an overestimation of the repulsion (underestimation of the attraction); but the comparison is a bit misleading, since the dielectric constant of our model (106) does not correspond to that of real water (78.4). A more appropriate test would be to compare the MC osmotic coefficients to the corresponding MD values, which can be directly obtained from the MD generated RDFs by using the KirkwoodBuff theory of solutions.65,66 However, this approach requires that the RDFs be known with a level of accuracy that is not easily obtained in simulations, hence the importance of developing coarse-grained potentials for such systems. Considering the fact that the thermodynamics of our system is governed primarily by the electrostatic interactions, we can expect the difference between both sets of data to be due, in large part, to the difference in dielectric constant (εr = 106 and 78.4). In fact, at infinite dilution, where the DebyeH€uckel limiting law becomes exact, the slope of the osmotic curve is inversely proportional to the dielectric constant: φ ∼ εr3/2. Since the cationanion interactions provide the dominant contributions, a screened Coulomb potential will lead to a decrease (increase) in the ion association (repulsion), and thus to an increase in the osmotic coefficient. Additional MC simulations were performed for εr = 78.4, with the same effective short-range potentials as before, in order to verify this hypothesis. These results are also presented in Figure 10, where it is clearly shown that a reduction in the dielectric constant produces a shift in the osmotic coefficient which is salt independent. The agreement with the experimental

ARTICLE

Figure 10. Comparison between experimental and MC osmotic coefficients ΦMM in the McMillanMayer reference frame: (left) experiments, (center) MC at εr = 106, and (right) MC at εr = 78.4.

results is improved, particularly at low concentrations (as expected), even though the effective potentials were derived for a solvent with a different dielectric constant. So far, no error analysis has been provided for our results, making it difficult to assess the reliability of the measured thermodynamic data. Although the MC results can be considered exact for a given set of pair potentials, with the relative error of the osmotic coefficients being less than 1%, the errors inherent to the inversion procedure must be taken into account. With this in mind, we have used two additional sets of effective potentials; for the first, the (reference) salt-specific Ln3þCl potentials are replaced by the average Ln3þCl potential (all other potentials being left unchanged), while for the second, the same substitution is performed on the Ln3þLn3þ potentials. The average potentials are shown in Figure 4. The osmotic coefficients obtained from MC simulations at c = 0.51 mol L1 (εr = 106), using these modified potentials, are compared with the reference values in Figure 11. The results obtained using an average Ln3þCl potential (i.e., varying only the Ln3þLn3þ potential) are seen to be essentially equivalent for all the lanthanoids, with a relative error of =5% with respect to the mean; even though the energy differences among the Ln3þLn3þ potentials are of the order of 1 kBT, and the onset of the short-range hard-core repulsion can vary by up to 1 Å. Since the energy fluctuations present in the short-range Ln3þLn3þ potentials are introduced during the inversion procedure, due to the noise in the original RDFs (see Figure 1), this variation in φ can be considered as an estimate of the error due to the uncertainty in these interactions. Although small, this change in φ is nonetheless surprising, since the shortrange energy fluctuations are roughly an order of magnitude weaker than the Coulombic repulsion which dominates the cationcation interactions. The reason for this is simple: variations in the cationcation potentials become increasingly important as the concentration is increased, due to packing effects within the fluid. A straightforward calculation will help illustrate this situation: consider a hard-sphere representation for the system, with the hard-core diameters σνμ given by the energy barriers in the total effective potentials βUeff νμ. Taking σClCl = 4.5 Å (Ueff ClCl(σClCl) = 1 kBT) and σLnCl = 6 Å (corresponding to the solvent-separated ion pair), we obtain σLnLn = 7.5 Å for the 4337



ARTICLE

Table 3. Osmotic Coefficients in the MM Reference Frame, as a Function of Salt Concentration (in Molarity Scale), Obtained from MC Simulationsa c (mol/L)

La3þ

Nd3þ

Eu3þ

Dy3þ

Lu3þ

0.01

0.8978/0.839 (0.841)

0.906/0.853 (0.836)

0.898/0.840 (0.836)

0.902/0.845 (0.844)

0.908/0.853 (0.841)

0.04 0.09

0.849/0.773 (0.788) 0.834/0.749 (0.777)

0.872/0.799 (0.785) 0.873/0.792 (0.777)

0.851/0.775 (0.788) 0.836/0.751 (0.781)

0.861/0.786 (0.795) 0.855/0.770 (0.787)

0.875/0.801 (0.791) 0.880/0.795 (0.784)

0.16

0.839/0.750 (0.787)

0.901/0.812 (0.789)

0.841/0.750 (0.796)

0.873/0.781 (0.800)

0.908/0.816 (0.800)

0.25

0.863/0.767 (0.813)

0.953/0.856 (0.816)

0.861/0.763 (0.825)

0.911/0.813 (0.830)

0.960/0.859 (0.831)

0.36

0.903/0.802 (0.853)

1.030/0.926 (0.855)

0.894/0.791 (0.867)

0.972/0.866 (0.875)

1.034/0.924 (0.876)

0.51

0.973/0.864 (0.915)

1.152/1.040 (0.916)

0.951/0.839 (0.931)

1.069/0.955 (0.943)

1.145/1.026 (0.944)

0.64

1.044/0.929 (0.975)

1.272/1.153 (0.974)

1.009/0.889 (0.993)

1.168/1.045 (1.009)

1.253/1.126 (1.009)

0.81

1.152/1.030 (1.060)

1.447/1.318 (1.057)

1.097/0.969 (1.080)

1.311/1.179 (1.102)

1.405/1.268 (1.101)

1.06 1.21

1.344/1.208 (1.199) 1.477/1.336 (1.289)

1.738/1.599 (1.191) 1.933/1.786 (1.279)

1.254/1.111 (1.221) 1.368/1.215 (1.311)

1.554/1.409 (1.253) 1.721/1.567 (1.350)

1.656/1.504 (1.250) 1.824/1.661 (1.348)

1.44

1.712/ (1.436)

2.259/ (1.420)

1.572/ (1.459)

2.004/ (1.509)

2.106/ (1.507)

The first and second values correspond to the results obtained for a dielectric constant of 106 (POL3) and 78 (experimental), respectively. The experimental values are given in parentheses. a

Figure 11. Variation of the osmotic coefficient at c = 0.51 mol kg1 due to changes in the cationcation and cationanion effective potentials: (circles) salt-specific reference potentials, (squares) reference Ln3þCl and average Ln3þLn3þ potentials, and (triangles) average Ln3þCl and reference Ln3þLn3þ potentials.

diameter of the hydrated cation; this corresponds to a packing fraction of =22% at 1 mol L1. The results for φ obtained using the average Ln3þLn3þ potentials are considerably closer to the reference values, reproducing the same ordering across the series, providing further evidence that the dominant contributions come from the cationanion interaction. A similar analysis was performed to study the sensitivity of φ to variations in the depth of the first minimum of the Ln3þCl potential (corresponding to the contact ion pair (CIP)). For small variations (j0.5kBT), the corresponding change in φ was always less than the error given by the uncertainty in the Ln3þLn3þ interaction potentials. Considering the fact that the Ln3þCl potentials are known with greater accuracy than the Ln3þLn3þ potentials, the estimated errors of ∼5% in φ due to the latter can be used as a conservative estimate for the error of the final osmotic coefficients (Figure 10). Given the good agreement between the RDFs obtained using the coarse-grained HNC potentials and those obtained from MD simulations, as well as the

favorable comparisons between these potentials and the exact PMF calculations, we assume that the errors due to the approximate nature of the HNC closure are negligible compared to the other sources of error. It is worth noting that the osmotic coefficient curves (Figure 10) roughly follow the same ordering as the association constants K0 (Table 2), with La3þ the most associated and Lu3þ the least. The inversion of Eu3þ and La3þ is understood to be a consequence of the error in the Ln3þLn3þ potentials; as expected, it does not appear in the series shown in Figure 11 which uses the same average Ln3þLn3þ potential for all the cations. However, the variation in the effective cation size would seem to favor an inverse ordering in K0, as evidenced by the shift in the first peak of the Ln3þCl RDF (Figures 1 and 8), since La3þ presents a larger size than Lu3þ (and thus a less attractive cationanion interaction). But, one must also consider the change in the strength of the effective cationanion interaction (Figure 4), which is also seen to decrease along the series.

’ CONCLUSIONS We have used numerical simulations to study the thermodynamic and structural properties of a series of aqueous lanthanoid chloride solutions. To achieve this, we have employed a coarse-graining strategy to derive an implicit solvent model from a molecular description. Finite concentration molecular dynamics simulations were performed, using a polarizable model of lanthanoid solutions, to compute the ionion radial distribution functions. These radial distribution functions were then inverted, using an approximate HNC closure, to obtain effective ionion interaction potentials. Although in principle they are density dependent, the difference with the exact (infinite dilution) potentials was shown to be relatively small. These solvent-averaged potentials were then used to perform implicit solvent Monte Carlo simulations. This allowed us to compute the osmotic coefficients, as a function of the concentration, for a representative series of lanthanoids (La, Nd, Eu, Dy, Lu), a task that would have been virtually impossible using a molecular model for both solute and solvent particles. The correspondence between the explicit and implicit solvent models (molecular dynamics vs Monte Carlo) was found to be very good, validating 4338


The Journal of Physical Chemistry B our inversion procedure, although a slight concentration dependence was observed when comparing the results obtained at 1.0. mol kg1. Directly from these inverted potentials, we were able to compute ion association constants, which are in agreement with those given by Ruas et al.,15 who consider them as fitting parameters within a BiMSA10 model. We arrive at essentially the same results starting from a molecular description of the system. The osmotic coefficients computed using our implicit solvent model show good agreement with the available experimental data,32 especially considering how sensitive this quantity is to the interaction potentials. The increase in the activity (salt activity and osmotic coefficient) across the series, as a function of the concentration, is globally reproduced. A decrease in the effective (solvated) cation “diameter” throughout the series is observed in both the Ln3þCl radial distribution functions and effective potentials, accompanied by a decrease in the strength of the cationanion solvent-averaged interaction. These two competing effects combine to give a net decrease in the ion association. This can also be interpreted in terms of the solvation of the contact ion pair, with a strong solvation corresponding to a high association constant. Although the solvent averaged (shortrange) potentials characterizing the strength of this solvation are small (roughly an order of magnitude smaller than the Coulomb potential in the CIP region), the osmotic coefficient is very sensitive to variations in these potentials. In order to distinguish between different salts on the basis of their activities, it is necessary to finely tune both the cationanion and the cationcation potentials. Even though the cationanion interactions provide the dominant contribution, seemingly small variations in the strongly repulsive cationcation interaction can lead to nonnegligible changes in the osmotic coefficient due to packing effects within the fluid.

’ ASSOCIATED CONTENT

bS

Supporting Information. Table S1 gives the hydration properties of Ln3þ cations obtained from MD simulations at 0.5 and 1.0 mol kg1. Figures S1 and S2 show the ionoxygen RDFs obtained at 0.5 and 1.0 mol kg1, respectively. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected] (J.J.M.); jean-francois.dufreche@ univ-montp2.fr (J.-F.D.); [email protected] (P.G.).

’ ACKNOWLEDGMENT This work was supported by the “Agence Nationale pour la Recherche” AMPLI (Approche Multi-echelle du Partage de solutes entre phases LIquides). The authors thank Olivier Bernard for useful discussions and the referees for helping us correct our PMF results and bringing to our attention several references on experimental and simulation studies of actinoid and lanthanoid solutions. ’ REFERENCES (1) Galbis, E.; Hernandez-Cobos, J.; Den Auwer, C.; Le Naour, C.; Guillaumont, D.; Simoni, E.; Pappalardo, R. R.; Marcos, E. S. Angew. Chem., Int. Ed. 2010, 49, 3811–3815. (2) Yang, T.; Bursten, B. E. Inorg. Chem. 2006, 45, 5291–5301.

ARTICLE

(3) Hagberg, D.; Bednarz, E.; Edelstein, N. M.; Gagliardi, L. J. Am. Chem. Soc. 2007, 129, 14136–14137. (4) Benay, G.; Schurhammer, R.; Wipff, G. Phys. Chem. Chem. Phys. 2010, 12, 11089–102. (5) Chaumont, A.; Wipff, G. Inorg. Chem. 2009, 48, 4277–4289. (6) D’Angelo, P.; Zitolo, A.; Migliorati, V.; Persson, I. Chem.—Eur. J. 2010, 16, 684–692. (7) Ruas, A.; Guilbaud, P.; Den Auwer, C.; Moulin, C.; Simonin, J.P.; Turq, P.; Moisy, P. J. Phys. Chem. A 2006, 110, 11770–11779. (8) Duvail, M.; Ruas, A.; Venault, L.; Moisy, P.; Guilbaud, P. Inorg. Chem. 2010, 49, 519–530. (9) Petit, L.; Vuilleumier, R.; Maldivi, P.; Adamo, C. J. Phys. Chem. B 2008, 112, 10603–10607. (10) Simonin, J.-P.; Bernard, O.; Blum, L. J. Phys. Chem. B 1998, 102, 4411–4417. (11) Vilarino, T.; Bernard, O.; Simonin, J.-P. J. Phys. Chem. B 2004, 108, 5763–5770. (12) Blum, L.; Bernard, O. J. Stat. Phys. 1995, 79, 569–583. (13) Bernard, O.; Blum, L. J. Chem. Phys. 1996, 104, 4746–4754. (14) Ruas, A.; Simonin, J.-P.; Turq, P.; Moisy, P. J. Phys. Chem. B 2005, 109, 23043–23050. (15) Ruas, A.; Moisy, P.; Simonin, J.-P.; Bernard, O.; Dufrêche, J.; Turq, P. J. Phys. Chem. B 2005, 109, 5243–5248. (16) David, F. H.; Vokhmin, V. New J. Chem. 2003, 27, 1627–1632. (17) David, F. H.; Vokhmin, V. J. Phys. Chem. A 2001, 105, 9704–9709. (18) David, F. H.; Vokhmin, V.; Ionova, G. J. Mol. Liq. 2001, 90, 45–62. (19) David, F. H.; Fourest, B. New J. Chem. 1997, 21, 167–176. (20) Duvail, M.; Vitorge, P.; Spezia, R. J. Chem. Phys. 2009, 130, 104501–13. (21) Molina, J.; Dufrêche, J.; Salanne, M.; Bernard, O.; Jardat, M.; Turq, P. Phys. Rev. E 2009, 80, 065103. (22) Cossy, C.; Helm, L.; Merbach, A. E. Inorg. Chem. 1988, 27, 1973–1979. (23) N€aslund, J.; Lindqvist-Reis, P.; Persson, I.; Sandstr€om, M. Inorg. Chem. 2000, 39, 4006–4011. (24) Allen, P.; Bucher, J. J.; Shuh, D. K.; Edelstein, N. M.; Craig, I. Inorg. Chem. 2000, 39, 595–601. (25) D’Angelo, P.; De Pan_lis, S.; Filipponi, A.; Persson, I. Chem.— Eur. J. 2008, 14, 3045–3055. (26) Rao, L.; Tian, G. Inorg. Chem. 2009, 48, 964–970. (27) Kowall, T.; Foglia, F.; Helm, L.; Merbach, A. E. J. Am. Chem. Soc. 1995, 117, 3790–3799. (28) Floris, F. M.; Tani, A. J. Chem. Phys. 2001, 115, 4750–4765. (29) Clavaguera, C.; Pollet, R.; Soudan, J. M.; Brenner, V.; Dognon, J. P. J. Phys. Chem. B 2005, 109, 7614–7616. (30) Villa, A.; Hess, B.; Saint-Martin, H. J. Phys. Chem. B 2009, 113, 7270–7281. (31) Brendebach, B.; Banik, N. L.; Marquardt, C. M.; Rothe, J.; Denecke, M.; Geckeis, H. Radiochim. Acta 2009, 97, 701–708. (32) Spedding, F. H.; Weber, H. O.; Saeger, V. W.; Petheram, H. H.; Rard, J. A.; Habenschuss, A. J. Chem. Eng. Data 1976, 21, 341–360. (33) Rard, J. A.; Weber, H. O.; Spedding, F. H. J. Chem. Eng. Data 1977, 22, 187–201. (34) Rard, J. A.; Shiers, L. E.; Heiser, D. J.; Spedding, F. H. J. Chem. Eng. Data 1977, 22, 337–347. (35) Rard, J. A.; Miller, D. G.; Spedding, F. H. J. Chem. Eng. Data 1979, 24, 348–354. (36) Rard, J. A.; Spedding, F. H. J. Chem. Eng. Data 1981, 26, 391–395. (37) Rard, J. A.; Spedding, F. H. J. Chem. Eng. Data 1982, 27, 454–461. (38) Rard, J. A. J. Chem. Eng. Data 1987, 32, 92–98. (39) Case, D. A. AMBER 10; University of California: San Francisco, 2008. (40) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. 4339



ARTICLE

(41) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (42) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757–3766. (43) Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. 1995, 99, 6208–6219. (44) Meng, E.; Kollman, P. J. Phys. Chem. 1996, 100, 11460–11470. (45) Clavaguera, C.; Dognon, J. Chem. Phys. 2005, 311, 169–176. (46) Chayes, J. T.; Chayes, L. J. Stat. Phys. 1984, 26, 471–488. (47) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. (48) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1992, 13, 1011–1021. (49) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1339–1350. (50) Souaille, M.; Roux, B. Comput. Phys. Commun. 2001, 135, 40–57. (51) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1986. (52) Molina, J.; Pierleoni, C.; Capone, B.; Hansen, J.-P.; Santos de Oliveira, I. Mol. Phys. 2009, 107, 535–548. (53) Hess, B.; Holm, C.; van der Vegt, N. J. Chem. Phys. 2006, 124, 164509–8. (54) Bhatia, A. B.; Thornton, D. E. Phys. Rev. B 1970, 2, 3004–3013. (55) Molina, J.; Duvail, M.; Guilbaud, P.; Dufrêche, J. J. Mol. Liq. 2010, 153, 107–111. (56) Frenkel, D.; Smit, B. Understanding Molecular Simulations; Academic Press: New York, 1996. (57) Hummer, G.; Gronbech-Jensen, N.; Neumann, M. J. Chem. Phys. 1998, 109, 2791–2997. (58) Persson, I.; D’Angelo, P.; De Pan_lis, S.; Sandstr€om, M.; Eriksson, L. Chem.—Eur. J. 2008, 14, 3056–3066. (59) Cossy, C.; Helm, L.; Powell, D. H.; Merbach, A. E. New J. Chem. 1995, 19, 27–35. (60) Duvail, M.; Spezia, R.; Vitorge, P. Chem. Phys. Chem. 2008, 9, 693–696. (61) Shannon, R. D. Acta Crystallogr. A 1976, 32, 751–767. (62) Duvail, M.; D’Angelo, P.; Gaigeot, M.-P.; Vitorge, P.; Spezia, R. Radiochim. Acta 2009, 97, 339–346. (63) Beuchat, C.; Hagbert, D.; Spezia, R.; Gagliardi, L. J. Phys. Chem. B 2010, 114, 15590–15597. (64) Lee, L. L. J. Mol. Liq. 2000, 87, 129–147. (65) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774–777. (66) Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1987, 86, 5110–5116.

4340


Atomistic Description of Binary Lanthanoid Salt Solutions: A Coarse

Recommend Documents