System-Size Dependence of Electrolyte Activity Coefficients in

Subscriber access provided by UNIV OF TASMANIA

Article

System-Size Dependence of Electrolyte Activity Coefficients in Molecular Simulations Jeffrey M. Young, and Athanassios Z. Panagiotopoulos J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b09861 • Publication Date (Web): 03 Jan 2018 Downloaded from http://pubs.acs.org on January 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a 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 31 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

System-Size Dependence of Electrolyte Activity Coefficients in Molecular Simulations Jeffrey M. Young and Athanassios Z. Panagiotopoulos∗ Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, USA E-mail: [email protected]

1

ACS Paragon Plus Environment

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

Abstract A systematic study of the dependence of electrolyte activity coefficients on simulation system size has been undertaken. Using implicit-solvent simulations for which calculations with low statistical uncertainty are feasible, it was found that the chemical potential for a NaCl model depends strongly on simulation system size at concentrations up to about 0.3 mol/L; system-size effects at higher concentrations are much smaller. Similar trends were confirmed in systems with an explicit solvent. System-size effects on the chemical potential, when uncorrected, can lead to systematic errors in the activity coefficient greater than 10%. The rigorous method to correct for such systemsize effects is to perform multiple simulations at each concentration and extrapolate to infinite system size. Unfortunately, this becomes impractical for explicit-solvent simulations at low concentrations, because of computational limitations that lead to large statistical uncertainties of the results. Somewhat counterintuitively, we find that lower systematic errors for the Henry’s law reference chemical potential are obtained by using simulations at higher concentrations, for which system-size effects are much smaller, to obtain estimates for the reference chemical potential. This is the case even though at these higher concentrations deviations from the Debye-H¨ uckel limiting law (or its empirical extensions) are greater than at lower concentrations.

Introduction Aqueous salt solutions are found in many systems relevant for biology, 1,2 carbon sequestration, 3,4 or geology. 5,6 These electrolyte solutions can be investigated using molecular simulations which allow the study of extreme temperatures and pressures and provide information about the solution’s behavior on a molecular scale. Good molecular models for electrolyte solutions should be able to predict both transport and thermodynamic properties. Among thermodynamic properties, perhaps the most important is the chemical potential (partial molar free energy of a component in a mixture) which determines chemical reactivity. Ac-

2


Page 2 of 31

Page 3 of 31 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


tivity coefficients are also measures of the chemical potential, expressing deviations from their ideal-solution values. Chemical potentials are often calculated for molecular models using Widom’s test particle insertion method 7 or grand canonical Monte Carlo. 8 These approaches, which require direct particle insertions into an existing system, fail for explicit-solvent models of electrolyte solutions because of their strong interactions and highly structured character. Other methods (such as gradual insertion of ions or measuring the pressure difference across a virtual semipermeable membrane) to reliably calculate solubility and activity coefficients of electrolyte solutions have only recently been established. 9–18 Interestingly, many commonly used models for electrolyte solutions have been found to be lacking in their prediction of chemical potentials and activity coefficients. 11,12,19 The calculation of activity coefficients for electrolyte solutions requires the simulation of low concentrations where electrostatic screening by ions is minimal and long-range electrostatic interactions dominate. Therefore, there are good reasons to believe that the size of the simulation box may influence the calculated activity coefficients, as suggested recently by Thompson and Sanchez. 20 Many studies have found a system-size dependence in other properties calculated from molecular simulations including critical parameters, 21 solid free energy, 22 surface tension, 23 diffusion coefficients, 24 nucleation, 25 and solubility. 15 The possible presence of finite-size effects is inherent in the nature of electrostatic interactions and cannot be circumvented by a clever choice of computational method, even if the exact magnitude of the effect will depend on the statistical ensemble used. As the creation of new electrolyte models is an active field of research, 26–31 being able to reliably calculate the activity coefficients of electrolytes will allow for the development of new models targeted to accurately predict solubility and activity coefficients over a range of temperatures. There have been a number of theoretical investigations into the effect of system size on the chemical potential of electrolytes 20,32–35 which we summarize here. Figueirido et. al. 32 and Hummer et. al. 33 created corrections for the effect of system size on the infinite-dilution free energy of one ion from simulations in a fixed volume (with side length L) and dielectric

3



Page 4 of 31

continuum. Figueirido et. al. assumed that the ion was a point charge with only Coulomb interactions, giving the first correction term in Equation 1 given below. 32 Hummer et. al. extended the analysis to account for the size of the ion, since the dielectric constant within the ion is not the same as in the system, giving an additional

R2 L3

dependence (where R is

the radius of the ion being charged). 33 Sakane et. al. 34 and Hunenberger and McCammon 35 derived a similar expression for inserting a single ion, as well as an expression for the insertion of an ion pair, for which case there is an additional term to account for the interactions with the periodic images of the complementary ion, the last term in Equation 1. 34,35 Recently, Thompson and Sanchez derived an expression for the system-size dependence of free energy at finite concentrations which, at infinite dilution, reduces to the same expression given in the other studies. 20 These studies have focused on ions in a dielectric continuum, and mostly at infinite dilution, so they may not be quantitatively applicable to calculations for the more pertinent explicit-solvent models. For our purposes, the important result from the prior theoretical studies is the expression for the infinite system-size chemical potential from the value at a finite system size when using Ewald sums, which add up the interactions of the ions in the simulation box with the interactions of all of the surrounding periodic images, to infinity. Assuming that the ions are point charges in a medium with a constant dielectric constant and the medium surrounding the macroscopic assembly of periodic boxes has an infinite dielectric constant (conducting boundary conditions) an approximation for the chemical potential for a infinite box from the calculated value at a finite box side length (L) is 34,35 2 2 ξ q+ + q− q+ q− µ(∞) = µ(L) − − Ψ(r+− ) 8π0 L

(1)

where ξ = −2.837297 is a constant that describes the Coulomb energy of a single charged particle in a cubic lattice, 0 is the permittivity of vacuum, is the relative permittivity of the system, q+ and q− are the charges of the ions being inserted into the system, and r+− is the distance between them. Ψ(r) is the effective potential between charges in a periodic

4


Page 5 of 31 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


system due to the Ewald summation, 34 " # X erfc (η|r + n|) X 4π 1 |k|2 π + Ψ(r) = exp − 2 + ik · r − 3 2 3 |k|2 4π0 n |r + n| L 4η Lη k6=0

(2)

where η is a convergence parameter and the sums are over all periodic cells with real space vector n and Fourier space vector k =

2πn . L2

Assuming the distance between two ions is either

constant (ions paired) or proportional to L (ions separated), this potential decays faster than 1/L as L increases. Therefore, the dominant system-size effect at large L is just the first correction term in Equation 1: 32–35 2 2 ξ q+ + q− µ(∞) = µ(L) − 8π0 L

(3)

These corrections are only valid at infinite dilution 33,34 and have been previously used to correct for system size effects of hydration free energy. 32–34,36,37 At higher concentrations, it is expected that the chemical potential will no longer scale as 1/L; 20 however, previous work has found that the chemical potential of ions in the canonical ensemble scales with 1/L at a concentration as high as 1.003 mol/L. 38 As shown later in this work, these low concentrations corrections scaling approximately with 1/L are the most important. In this paper, we investigate the effect of system size on the chemical potential and activity coefficients first using an implicit-solvent model, where the uncertainty is lower and the calculations are faster. We compare these simulations to a model with explicit water to ensure that the same trends hold true. Lastly, we propose a method to correct for this system-size effect for use in future simulations. The surprising outcome from the work is that it is preferable to use simulations at higher concentrations to obtain estimates of the Henry’s law reference chemical potentials, because system-size effects are much smaller at these conditions, and extrapolation of the chemical potential to infinite system size is possible.

5



Page 6 of 31

Methods The NaCl model developed by Joung and Cheatham 39 (JC) in conjunction with the SPC/E water model 40 were chosen for the present work. The reason for this choice is that this model combination has been extensively studied in prior literature. 11,12,41–45 Since only the shorter-ranged van der Waals interactions differ between most monovalent ion models, we expect that the general trends we observe are insensitive to the details of the electrolyte and water models used. The parameters for the models are given in Supporting Information. Cross interactions follow the Lorentz-Berthelot combining rules. Simulations were performed in both an implicit- and an explicit-solvent environment. The implicit water simulations consisted of Na+ and Cl− ions placed in a dielectric continuum with a dielectric constant of 73, equal to the dielectric constant of SPC/E water at 298 K. 12 The dielectric constant of SPC/E water was used instead of the experimental dielectric constant of water in order to allow for direct comparison between the implicit and explicit solvent models. The explicit solvent simulations used SPC/E water molecules for the solvent. The mean ionic activity coefficient γ of a 1:1 electrolyte is defined as

βµ = βµ† (P, T ) + 2 ln(mγ)

where β =

1 , kb T

(4)

m is the concentration in moles of salt per kg water, µ is the chemical poten-

tial at that concentration, and µ† is the Henry’s law reference chemical potential, a function of pressure and temperature. For the implicit-solvent simulations, since there is no water present, the concentration is instead represented as M (moles salt/liter solution). These two units of concentration are approximately equivalent since the mass of the ions is small at low concentrations and the solution density is around 1000 kg/m3 . When calculating the activity coefficient from Equation 4, the temperature and the concentration are input parameters in canonical-ensemble simulations, requiring only the calculation of the chemical potential at the given concentration and the Henry’s law reference chemical potential to solve for the

6


Page 7 of 31 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


activity coefficient. Alternatively, one could use grand canonical Monte Carlo (GCMC) simulations in which temperature and chemical potential are known and the concentration is a calculated (fluctuating) quantity, but one still requires an estimate of µ† to obtain the activity coefficients. In this work, both molecular dynamics (MD) and GCMC simulations were used. For the MD simulations, the chemical potential at a given concentration was calculated using the thermodynamic integration approach of Mester and Panagiotopoulos. 11,12 The free energy change of an ion pair insertion was determined using Bennet’s acceptance ratio. 46 The ion pair was inserted slowly, with the potential of the ion pair being ramped up stepwise. The interactions of the added ion pair (indices i, j) with the other molecules in the system can be expressed as " U =λ

# X

"

(Uik,Vdw + Ujk,Vdw ) + Uij,Vdw + φ

k6=i,j

# X

(Uik,Coul + Ujk,Coul ) + Uij,Coul

(5)

k6=i,j

where  Ulm,Vdw = 4lm 

!12

σlm 6 [0.5σlm (1

− λ) +

6 rlm ]

1/6

−

!6 

σlm 6 [0.5σlm (1

− λ) +

6 1/6 rlm ]



(6)

lm and σlm are the soft-core Lennard-Jones interactions between species l and m, and

Ulm,Coul =

ql qm 4π0 rlm

(7)

where ql and qm are the charges of the atoms and rlm is the distance between them. λ and φ are parameters in the range [0, 1]. As they change from 0 to 1, the inserted ion pair changes from not interacting with the system to fully interacting. Van der Waals forces for the inserted ion pair are turned on first, using a soft-core potential to avoid the discontinuity when they are first turned on. For the implicit-solvent simulations, 5 equally spaced values of λ are used (λ = [0, 0.25, . . . , 1]), while the explicit-solvent simulations used 11 equally spaced 7



Page 8 of 31

steps in λ. The Coulomb interactions are then turned on, using 5 equally spaced steps in φ for the implicit simulations and 21 steps for the explicit simulations. The free energy difference between each step in λ or φ is calculated using Bennet’s acceptance ratio, 46 which has been found to be one of the most efficient methods to determine free energy differences. 47–49 From the definition of chemical potential,

µ=

∆F ∂F ≈ ∂N ∆N

(8)

where F is the free energy, the chemical potential can be found using the finite difference approximation for the insertion of one ion pair. For this work, the concentration of the chemical potential is taken as the average of the concentrations before and after the insertion. The ideal chemical potential for the implicit-solvent calculations is

βµIG = ln

NNaCl Λ3Na+ NNaCl Λ3Cl− + ln V V

(9)

Here NNaCl is the number of ion pairs (also taken as the average of the insertion), Λ is the de Broglie wavelength, and V is the total volume of the system. The ideal chemical potential for the explicit-solvent simulations is the one used in Mester and Panagiotopoulos, 11,12

βµIG = µ0Na+ + µ0Cl− + 2 ln

NNaCl βP 0 hV i

(10)

where the standard chemical potentials µ0Na+ = 574.317 kJ/mol and µ0Cl− = −240.167 kJ/mol are from the NIST-JANAF tables 50 and the standard pressure P 0 is 1 bar. It should be noted that this choice of the standard chemical potentials only affects the reported values of the electrolyte chemical potential. The activity coefficients, representing a difference between the infinite-dilution and finite-concentration states, are unaffected by this choice. MD simulations were performed using GROMACS 5.1.3 51 in the canonical ensemble using a time step of 2 fs for 10 million steps, a cutoff of 2.0 nm (implicit solvent) or 0.9

8


Page 9 of 31 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


nm (explicit solvent) for Lennard-Jones and electrostatic interactions, and a Nose-Hoover thermostat with a coupling constant of 1 ps. The temperature of the simulations was 298.15 K. Explicit solvent simulations used a Parrinello-Rahman barostat set to a pressure of 1 bar with a 2 ps coupling constant. Standard long-range corrections were used for the LennardJones interactions, while particle mesh Ewald summations with a Fourier-space grid spacing of 0.1 nm, a 6th order interpolating function, and a relative tolerance of 1 × 10−6 were used for the long-range electrostatic interactions. For the GCMC simulations, which are only feasible with an implicit solvent, the chemical potential is an input, as stated earlier. To obtain the chemical potential at a specific concentration, the secant method was used to solve for the chemical potential that would give the desired concentration. Simulations were performed using Cassandra 1.0 52 which was modified to insert pairs of ions randomly into the system. The simulations were run for one million steps with data collection every thousand steps. The first 100,000 steps were discarded for equilibration. The temperature was set to 298.15 K. As for the MD simulations, a 2.0 nm cutoff was used for the interactions with standard long-range corrections to the Lennard-Jones interactions, and Ewald summations for the Coulomb interactions. Since µ† cannot be directly calculated, one means of obtaining it is by fitting a simulated low concentration chemical potential to the Davies equation, as in Mester and Panagiotopoulos. 12 The Davies equation 53 is:

ln γ = −A where, for a 1:1 electrolyte,

√

m √ − 0.2m 1+ m

√ qe3 2NA ρs A= 8π(0 kb T )3/2

(11)

(12)

where qe is the elementary charge, NA is Avogadro’s number, and ρs is the density of the solution. This is an empirical extension of the Debye-H¨ uckel limiting law that is valid to higher concentrations, predicting almost the exact experimental activity coefficient for NaCl

9



Page 10 of 31

at 0.01 mol/kg. 54 A is the same constant as in the Debye-H¨ uckel limiting law so that Equation 11 reduces to the Debye-H¨ uckel equation at infinite dilution. Combining equations 4 and 11, we can obtain µ† from

†

βµ = βµ(m) − 2 ln m + 2A

√

m √ − 0.2m 1+ m

(13)

where m is a suitably low concentration, and µ(m) is the simulation chemical potential at that concentration. Data for the chemical potentials at the lowest concentration studied, namely 0.006 mol/L, were used to obtain µ† for the implicit model. The use of the empirical Davies equation is not central to the purposes of the present study. The theoretical Debye-H¨ uckel limiting law could be used for extrapolating to infinite dilution, at the cost of an increased systematic error because of the deviations from the limiting law at lower concentrations. Statistical uncertainties were estimated by performing five independent runs of each state point. The uncertainty in the raw data was propagated to the parameters in the fits by assuming the data was normally distributed around each point. Random points drawn from these distributions were fit 1,000 times to estimate the uncertainty in the fit.

Results and Discussion Implicit-solvent simulations of NaCl have been performed to evaluate the effect of the system size on the chemical potential and activity coefficients in a simplified system. These implicitsolvent solutions allow us to study large system sizes and low concentrations without the significant computational cost of sampling solvent degrees of freedom. A concentration of 0.01 molal requires approximately 5,000 water molecules per ion pair for explicit-solvent simulations, whereas this concentration can be reached in an implicit solvent by simply having a mostly empty simulation box. Implicit-solvent simulations are faster and have lower statistical uncertainties relative to explicit-solvent simulations. In this work, we are not trying to replicate the activity coefficients in the explicit solvent as accurately as possible. 10


Page 11 of 31 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


Previous works have used distance or concentration dependent dielectric constants in implicit solvent models to obtain accurate activity coefficients. 55–59 Instead, we are investigating the effect of system size on a simple model with a constant dielectric constant. As stated earlier (Equation 3), the chemical potential of an ion pair in a periodic system is approximately proportional to the inverse of the box length, at large L and infinite dilution. At the lowest concentration studied, a 1/L dependence is still expected to dominate even though we are not quite at infinite dilution. Chemical potentials from the MD simulations were extrapolated to an infinite system size using a simple linear scheme,

µ(L) =

a + µ(∞) L

(14)

where a is the slope of the extrapolation. The side length of the box, L, was varied from 5 nm to 12.5 nm for most simulations, with some calculations with L ≈ 4 nm for the highest concentrations. The concentration was kept constant between system sizes so that each extrapolation could be performed at a given concentration. Figure 1a and Figure 1b show the extrapolations for concentrations of 0.0066 mol/L and 0.14 mol/L respectively, indicating that the assumption of linear scaling holds well for the MD results. Data for other concentrations are given in Supplementary Information. The value of the chemical potential at infinite system size (y-intercept of the lines) was constrained to be the same for all of the fits of data at a specific concentration. For Ewald summations, the dielectric constant surrounding the infinite system can be set to infinity (conducting boundary conditions) or to a specified value. Since it is possible that some of the system-size effect is due to the boundary conditions, we also performed simulations with the dielectric constant surrounding the system set to 73, the same as the dielectric constant of the system (∞ = 73 boundary conditions). Since the ∞ = 73 boundary conditions only affect the Coulomb interactions, it was assumed that the free energy change of adding Lennard-Jones interactions was the same as with conducting boundary conditions.

11


µ (kJ/mol)

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

µ (kJ/mol)


84.0 84.5 85.0 85.5 86.0 86.5 87.0 87.5 88.0

a) 0.0066 mol/L

72.85 72.90 b) 0.14 mol/L 72.95 73.00 73.05 73.10 MD Conducting BC MD ² = 73 BC 73.15 GCMC Conducting BC 73.20 0.00 0.05 0.10 0.15 0.20 0.25 1/L (1/nm) ∞

Figure 1: Extrapolation of chemical potential to infinite system size for implicit-solvent simulations at concentrations of a) 0.0066 mol/L and b) 0.14 mol/L. BC stands for boundary conditions. Error bars and shading indicate one standard deviation.

12


Page 12 of 31

Page 13 of 31 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


As seen in Figures 1a and b, ∞ = 73 boundary conditions lead to a higher chemical potential for all system sizes compared to conducting boundary conditions. However, since the Ewald summation boundary conditions do not matter for the chemical potential in an infinite system, both types of boundary conditions can be extrapolated to yield the same value at an infinite box side length. GCMC simulations were run up to a box side length of 17.5 nm to ensure that they converged to the same limit as the MD simulations. Since the previous expressions for the scaling of chemical potential with system size were derived at a constant N , 20,32–35 the chemical potential calculated from the GCMC simulations will not necessarily scale with 1/L. Instead, since properties are expected to differ by a factor of 1/N between statistical ensembles at a finite system size, the extrapolation for the GCMC simulations was done using µ(L) =

b a + 3 + µ(∞) L L

(15)

where a and b are fitted parameters. The 1/L3 term is really a 1/N dependence and can be significant for small N at low concentrations, unlike the 1/L3 terms neglected in Equation 3. Writing it as a 1/L3 dependence means that b is proportional to 1/ρ (where ρ is the density of the ions). For the lowest concentrations, the density is extremely small, causing b to be on the order of 300 kJ nm3 /mol whereas a is on the order of -1 kJ nm/mol. In Figures 1a and b, the system-size dependence decreases when the concentration is increased. Note that the scale on the y-axis changes drastically between Figures 1a and 1b. At the lowest concentration of 0.0066 mol/L, the system-size effect is very strong, about an order of magnitude larger than the uncertainty in the chemical potential. There could be a number of causes for this system-size effect. As explained earlier, there is an expected 1/L dependence of the chemical potential due to the periodicity introduced by the Ewald summation treatment for the long-range electrostatics. An additional cause for the systemsize effect could be the use of the finite difference method to estimate the chemical potential. At high concentrations, the number of ion pairs tends to be higher, meaning that the insertion 13



of one ion pair better approximates the derivative. However, at low concentrations, there are only a few ion pairs, so the insertion of one pair is no longer a good approximation for the derivative. This error can be thought of as an error in the concentration, instead of in the chemical potential, since, by the mean value theorem, the true concentration must be somewhere between the uncoupled and fully coupled concentrations. As previously noted, for this work, the concentration was taken as the midpoint of these two concentrations. There was no effect of the system size on the free energy of adding the van der Waals interactions, indicating that these results are not specific to our choice of model, but general for all monovalent ions. At low concentrations, the GCMC simulations deviate drastically from the MD simulations, possibly because the GCMC simulations, which have a variable concentration, frequently sample configurations with a very few ion pairs in the simulation box, leading to strong dipole correlations at the lowest concentrations. MD simulations, by contrast, have a fixed number of ion pairs and never sample these very dilute states. The limiting slope of the fitted curves (parameter a) is plotted for different concentrations in Figure 2. The slopes for the MD conducting and GCMC conducting simulations were fit with simple exponentials, while the limiting slope of the MD ∞ = 73 boundary condition simulations was fit with a sum of two exponentials. The magnitude of the limiting slope is large for low concentrations, but quickly decreases to near zero. For concentrations above 0.3 mol/L, there is essentially no system-size effect using conducting boundary conditions. This decrease in the system-size effect is most likely due to electrostatic screening from the additional ions in the system at higher concentrations. Both the GCMC and MD simulations with conducting boundary conditions show the same behavior of the limiting slope, different than the slopes from the simulations with ∞ = 73 boundary conditions. This agreement indicates that the 1/L system-size dependence is likely due to the Ewald summations. At higher concentrations, the extrapolation slopes with conducting boundary conditions appear to be systematically slightly negative; however it is hard to tell within the uncertainty. For low concentrations, ∞ = 73 boundary conditions give a better result than the conducting

14


Page 14 of 31

Page 15 of 31

boundary conditions. However, as the concentration is increased, the slope of the extrapolation for the ∞ = 73 boundary conditions decays to zero more slowly, although it does pass through a zero slope around 0.08 mol/L likely because of a cancellation of error. Vacuum boundary conditions where ∞ = 1 are another alternative, but initial calculations showed that the system-size effect would be even larger than with ∞ = 73 boundary conditions. The theoretical value (from Equation 3) for the slope at infinite dilution is slightly below the y-intercept of the exponential fits to the conducting boundary condition data. This difference is likely due to the higher order terms neglected in Equation 3. 2 1 a (kJnm/mol)

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 -1 -2 MD Conducting BC MD ∞ = 73 BC

-3 -4

GCMC Conducting BC Explict MD Conducting BC Slope from Eq. 3

-5 -6 0.0

0.5

√

1.0

1.5

2.0

2.5

Concentration

Figure 2: Parameter a in Equations 14 and 15 (limiting slope of the extrapolation) versus concentration. Concentration is in units of mol/L for implicit-solvent simulations and units of mol/kg for explicit-solvent simulations. Error bars and shading indicate one standard deviation. Simulations with explicit water were also performed to determine if the system-size effects seen for the implicit-solvent simulations are reproduced in this more realistic case. The systems studied had 500 to 5500 water molecules (L ≈ 2.5 - 5.5 nm). The range of L differed from the implicit-solvent simulations since larger simulations with an explicit solvent are prohibitively expensive. Because the explicit-solvent simulations are at a fixed pressure, there is a small uncertainty in the box size. This uncertainty is less than 0.001 nm for all of the simulations. These simulations were again extrapolated using Equation 14. As seen in

15



Figure 3, at a concentration of 0.056 mol/kg, the minimum possible with 500 water molecules, there is some effect of the box side length on the chemical potential of the salt. However, since the uncertainty is much larger in the explicit-solvent calculations, the extrapolation to infinite size is very uncertain. As noted by Mester and Panagiotopoulos, 12 the chemical potential of all of the system sizes agree within two standard deviations. However, this is due to the large uncertainty in the simulations, not the lack of a system-size effect. The slope of the extrapolation for three concentrations is plotted in Figure 2 along with the implicitsolvent data. The explicit extrapolations slopes are comparable to the extrapolation slopes of the implicit solvent, but are very uncertain. This indicates that it is difficult to perform an extrapolation to an infinite system size for an explicit-water model without introducing large amounts of error or requiring impractically long simulations.

405 406

µ (kJ/mol)

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 31

407 408 409 0.0

0.1

0.2 0.3 1/L (1/nm)

0.4

Figure 3: Extrapolation of chemical potential to infinite system size for explicit-solvent simulations at a concentration of 0.056 mol/kg. Error bars and shading indicate one standard deviation. Uncertainty in L due to volume fluctuations is smaller than symbol size. The strong system-size dependence of the chemical potential can lead to very inaccurate activity coefficients. From Equation 4, it is seen that the activity coefficients depend on the value of µ† , as well as the simulated chemical potential at a finite concentration. Using the method presented above, µ† was estimated from the chemical potentials obtained at the lowest concentration. Since the lowest concentration also has the largest system-size effect, 16


Page 17 of 31 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


activity coefficients uncorrected for system-size effects will be inaccurate for all concentrations. Other methods to determine µ† , such as fitting an activity coefficient model to all concentrations, will have similar problems since µ† is still sensitive to the lowest concentration simulations. Figure 4 shows the prediction of the activity coefficient with an implicit solvent using both the smallest box size in this work (L = 5 nm) and the chemical potentials extrapolated to infinite box length. Using a fixed finite box size gives unphysical activity coefficients that are substantially higher than those for the infinite system. At high concentrations the difference is about an order of magnitude larger than the uncertainty. Previous studies have used increasingly larger box sizes for lower concentrations, but only because a larger box was required to reach the desired concentration. 11,12 Using a similar method does give more accurate activity coefficients (as seen in Figure 4), but there are still systematic deviations from the true value at high concentrations. In addition, for explicit-solvent simulations, considerably smaller box sizes are typically used than those in Figure 4. Mester and Panagiotopoulos 11,12 used a box side length of about 5.3 nm for their largest system and around 2.5 nm for most of their calculations. Therefore, the system-size effects in these earlier studies could be even larger than those shown in Figure 4. Ideally, the chemical potential at every concentration could be extrapolated to infinite system-size. Unfortunately, this is not possible for low concentrations in explicit-solvent models which require 5,000 water molecules per ion pair to reach concentrations of around 0.01 mol/kg. An extrapolation to infinite system size will require simulations with multiple ion pairs and therefore a prohibitively large numbers of water molecules. Instead, the best way to minimize the system-size error is to use a higher concentration to obtain µ† from the Debye-H¨ uckel limiting law (or the Davies equation). A higher concentration has a smaller system-size effect and is possible to extrapolate to infinite system size. As an example, using our data at 0.056 mol/kg (lowest concentration possible with 500 water molecules), we can reanalyze the data in Mester and Panagiotopoulos, 11,12 as done in Figure 5. The extrapolated infinite-size chemical potential at 0.056 mol/kg, shown in Figure 3, is used with the Davies

17


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

lnγ


0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0

Page 18 of 31

Infinite system size L = 5 nm Varying system size

0.5

1.0

1.5

M 1/2

2.0

2.5

Figure 4: Activity coefficients for the implicit-solvent model. The varying system size method uses L = 10 nm for the lowest concentration, L = 7.5 for the second lowest, and L = 5 nm for all other concentrations. The infinite system size is fit using the expression in Mester and Panagiotopoulos. 11,12 M is in units of mol/L. Error bars indicate one standard deviation and are not displayed when smaller than the points. equation to obtain µ† . At a concentration of 0.05 mol/kg, the experimental value for ln γ is only 0.006 units different from that predicted by the Davies equation, and the activity coefficient predicted by the fit to the implicit model at infinite system size is 0.02 units different from the Davies equation. This error is considerably smaller than the 0.05-0.1 units of error between the fixed L = 5 and the infinite system-size activity coefficients, as shown in Figure 4, and is the same order of magnitude as the uncertainty of the implicit model. For the implicit model, fitting to the Davies equation at a similar concentration would give activity coefficients about halfway between those of the infinite system size and the varying system size in Figure 4. Figure 5 shows that the activity coefficients calculated using the extrapolated chemical potential at 0.056 mol/kg fall slightly below the data in Mester and Panagiotopoulos, as expected from Figure 4. The difference between the two calculations of ln γ is approximately 0.15 units at the high concentrations. The uncertainty in each data point was calculated by combining the uncertainty in the chemical potential at that concentration with the uncertainty in µ† using the standard propagation-of-error formula given as Equation 38 in

18


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


lnγ

Page 19 of 31

1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 0.0

Reanalyzed data Mester and Panagiotopoulos

0.5

1.0

1.5

2.0

2.5

m 1/2 Figure 5: Activity coefficients for the explicit-solvent model, as a function of molality, m (in mol/kg). Data points are from Refs. 11 and 12. The reanalyzed data have their concentration shifted to the midpoint of the ion pair insertion and are fitted to the functional form used by Mester and Panagiotopoulos. 11,12 Error bars indicate one standard deviation. Nezbeda, Mou˘cka, and Smith. 16 This gives values larger than those in Mester and Panagiotopoulos 11,12 due to the higher uncertainty in µ† . Note that the uncertainty given here for the data from Refs. 11 and 12 also includes the propagation of the uncertainty in both µ and µ† , giving values for the activity coefficient uncertainty different from those in the original articles. The higher uncertainty of the reanalyzed data means that most of the error represented in the reanalyzed data’s error bars is a systematic error that would affect all of the calculated activity coefficients. The reanalyzed activity coefficients at the two lowest concentrations studied by Mester and Panagiotopoulos are unrealistic because of the strong system-size effect at the low concentrations and are omitted from the figure. As seen in Figure 2, extrapolating to infinite system size is necessary up to a concentration of approximately 0.3 mol/kg. Above this concentration little error will be introduced by using the chemical potential of the finite system instead of that of the infinite system. Therefore, it is valid to plot the data from Mester and Panagiotopoulos above 0.3 mol/kg even though their simulations were all at a finite system size of 500 water molecules. The value for µ† obtained from this method is compared to that calculated in previous works in Table 1. The value for µ† calculated by Mester and Panagiotopoulos 11 was found

19



by fitting a concentration of 0.01 mol/kg to the Davies equation. That of Mou˘cka, L´ısal, and Smith 10 was calculated by a weighted regression simultaneously fitting the activity coefficients for both NaCl and H2 O to the whole range of concentrations. This likely reduces the system size effect since it considers all of the chemical potentials (and therefore gives a µ† less dependent on system size), taking into account the larger uncertainties of low concentration data. Table 1: Comparison of µ† for the JC-SPC/E Model source this work Ref. 10 Ref. 11 experiments 60

µ† (kJ/mol) -390.8(9) -391.278 -391.7(2) -393.133

Conclusions There is a large effect of the system size on the chemical potential of electrolyte solutions, due to the presence of long-range electrostatic interactions. MD simulations using Ewald summations with both conducting and ∞ = 73 boundary conditions predict the same chemical potential as GCMC simulations when extrapolated to an infinite system size. System-size effects can be seen most clearly in implicit-solvent simulations for which the uncertainty in the chemical potential is an order of magnitude smaller than for explicit-solvent simulations. The system-size dependence increases nearly exponentially as the concentration decreases, but reaches a finite value at the limit of infinite dilution. Since the lowest concentration is typically used in simulations to obtain µ† , not correcting for system-size effects leads to incorrect estimates for the activity coefficients at all concentrations. In implicit-solvent and explicit-solvent simulations, this system-size effect, if uncorrected, causes a systematic error of over 0.1 units in the logarithm of the activity coefficients (10% in the activity coefficient values). As the dielectric constant appears in the denominator of Equation 3, a decrease in

20


Page 20 of 31

Page 21 of 31 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 dielectric constant (due to using a different solvent or increasing the temperature) will cause the system size effect to increase. In explicit-solvent simulations it is not feasible to extrapolate to an infinite system size at very low concentrations, which require a large number of solvent (water) molecules to be present in the simulation box. Therefore, we recommend using a higher concentration to obtain µ† from the Davies equation (or a similar expression). This concentration should be as low as possible but still high enough to allow an extrapolation to infinite system size. For the specific example of NaCl in water studied here, we found it effective to extrapolate to infinite system-size chemical potential at a concentration of about 0.056 mol/kg (which corresponds to 1 ion pair per 1000 water molecules), and using this chemical potential with the Davies equation to estimate µ† . This concentration level does introduce some additional systematic error relative to lower concentrations caused by the Davies equation; however, this error is significantly smaller than the error due to the system-size effect. We find that calculations of activity coefficients need to rigorously extrapolate to infinite system-size chemical potentials up to approximately 0.3 mol/kg. Above this concentration level, systemsize effects are minimal for 1:1 electrolytes because of screening of electrostatic interactions; it is valid to use chemical potentials calculated at a finite system size without correction. For NaCl, experimental data for ln γ differ from the Davies equation by only 0.006 units at m = 0.05 mol/kg. For other salts, the approximate error from fitting to the Davies equation at a higher concentration can be estimated by comparing the Davies equation to the experimental activity coefficients. These deviations are expected to be of similar magnitude also for “reasonable” atomistic models of water and ions. These considerations are important for future evaluations of atomistic models with respect to their ability to properly represent activity coefficients of electrolyte solutions. As well as mean ionic activity coefficients, strong system size effects are also expected to be present for single ion activity coefficients. In addition, calculations of solubilities for sparingly soluble salts using the chemical potential route will also need to take into account these strong system-size effects

21



on calculated chemical potentials in the solution phase.

Acknowledgments Financial support for this work was provided by the Office of Basic Energy Sciences, U.S. Department of Energy, under Award DE-SC0002128. Computing resources were provided by the Princeton Institute for Computational Science and Engineering (PICSciE).

Supporting Information Included in the supporting information are the parameters for the models used in this work. In addition, the extrapolation of the chemical potential to infinite system size is shown for all of the concentrations studied, using both the implicit and explicit solvent. Also provided are tables of the chemical potential and activity coefficients for different system sizes and concentrations, as well as a graph comparing the simulated activity coefficients to experiments.

References (1) Dill, K. A. Dominant forces in protein folding. Biochemistry 1990, 29, 7133–7155. (2) Auffinger, P.; Cheatham, T. E.; Vaiana, A. C. Spontaneous formation of KCl aggregates in biomolecular simulations: A force field issue? J. Chem. Theory Comput. 2007, 3, 1851–1859. (3) Gale, J. Geological storage of CO2: What do we know, where are the gaps and what more needs to be done? Energy 2004, 29, 1329–1338. (4) De Silva, G. P. D.; Ranjith, P. G.; Perera, M. S. A. Geochemical aspects of CO2 sequestration in deep saline aquifers: A review. Fuel 2015, 155, 128–143. 22


Page 22 of 31

Page 23 of 31 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


(5) Sherman, D. M.; Collings, M. D. Ion association in concentrated NaCl brines from ambient to supercritical conditions: results from classical molecular dynamics simulations. Geochem. Trans. 2002, 3, 102–107. (6) Di Maio, C. Exposure of bentonite to salt solution: osmotic and mechanical effects. Géotechnique 1996, 46, 695–707. (7) Widom, B. Some topics in the theory of fluids. J. Chem. Phys. 1963, 39, 2808–2812. (8) Norman, G. E.; Filinov, V. S. Investigation of phase transitions by a Monte-Carlo method. High Temp. 1969, 7, 216. ˇ (9) Mouˇcka, F.; L´ısal, M.; Skvor, J.; Jirsák, J.; Nezbeda, I.; Smith, W. R. Molecular simulation of aqueous electrolyte solubility. 2. Osmotic ensemble Monte Carlo methodology for free energy and solubility calculations and application to NaCl. J. Phys. Chem. B 2011, 115, 7849–7861. (10) Mouˇcka, F.; Nezbeda, I.; Smith, W. R. Molecular simulation of aqueous electrolytes: Water chemical potential results and Gibbs-Duhem equation consistency tests. J. Chem. Phys. 2013, 139, 124505. (11) Mester, Z.; Panagiotopoulos, A. Z. Temperature-dependent solubilities and mean ionic activity coefficients of alkali halides in water from molecular dynamics simulations. J. Chem. Phys. 2015, 143, 044505. (12) Mester, Z.; Panagiotopoulos, A. Z. Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations. J. Chem. Phys. 2015, 142, 044507. (13) Manzanilla-Granados, H. M.; Saint-Mart´ın, H.; Fuentes-Azcatl, R.; Alejandre, J. Direct coexistence methods to determine the solubility of salts in water from numerical simulations. Test case NaCl. J. Phys. Chem. B 2015, 119, 8389–8396.

23



(14) Benavides, A. L.; Aragones, J. L.; Vega, C. Consensus on the solubility of NaCl in water from computer simulations using the chemical potetial route. J. Chem. Phys. 2016, 144, 124504. (15) Espinosa, J. R.; Young, J. M.; Jiang, H.; Gupta, D.; Vega, C.; Sanz, E.; Debenedetti, P. G.; Panagiotopoulos, A. Z. On the calculation of solubilities via direct coexistence simulations: Investigation of NaCl aqueous solutions and Lennard-Jones binary mixtures. J. Chem. Phys. 2016, 145, 154111. (16) Nezbeda, I.; Mouˇcka, F.; Smith, W. R. Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility. Mol. Phys. 2016, 114, 1665–1690. (17) Kohns, M.; Reiser, S.; Horsch, M.; Hasse, H. Solvent activity in electrolyte solutions from molecular simulation of the osmotic pressure. J. Chem. Phys. 2016, 144, 084112. (18) Kohns, M.; Schappals, M.; Horsch, M.; Hasse, H. Activities in aqueous solutions of the alkali halide salts from molecular simulation. J. Chem. Eng. Data 2016, 61, 4068–4076. (19) Mouˇcka, F.; Nezbeda, I.; Smith, W. R. Molecular force fields for aqueous electrolytes: SPC/E-compatible charged LJ sphere models and their limitations. J. Chem. Phys. 2013, 138, 154102. (20) Thompson, J. P.; Sanchez, I. C. System-size effects in ionic fluids under periodic boundary conditions. J. Chem. Phys. 2016, 145, 214103. (21) Panagiotopoulos, A. Z. Molecular simulation of phase coexistence: Finite-size effects and determination of critical parameters for two- and three-dimensional Lennard-Jones fluids. Int. J. Thermophys. 1994, 15, 1057–1072. (22) Polson, J. M.; Trizac, E.; Pronk, S.; Frenkel, D. Finite-size corrections to the free energies of crystalline solids. J. Chem. Phys. 2000, 112, 5339–5342. 24


Page 24 of 31

Page 25 of 31 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


(23) Errington, J. R. Evaluating surface tension using grand-canonical transition-matrix Monte Carlo simulation and finite-size scaling. Phys. Rev. E 2003, 67, 012102. (24) Yeh, I.-C.; Hummer, G. System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions. J. Phys. Chem. B 2004, 108, 15873–15879. (25) Wedekind, J.; Reguera, D.; Strey, R. Finite-size effects in simulations of nucleation. J. Chem. Phys. 2006, 125, 214505. (26) Kolafa, J. Solubility of NaCl in water and its melting point by molecular dynamics in the slab geometry and a new BK3-compatible force field. J. Chem. Phys. 2016, 145, 204509. (27) Benavides, A. L.; Portillo, M. A.; Chamorro, V. C.; Espinosa, J. R.; Abascal, J. L. F.; Vega, C. A potential model for sodium chloride solutions based on the TIP4P/2005 water model. J. Chem. Phys. 2017, 147, 104501. (28) Kiss, P. T.; Baranyai, A. A new polarizable force field for alkali and halide ions. J. Chem. Phys. 2014, 141, 114501. (29) Reiser, S.; Deublein, S.; Vrabec, J.; Hasse, H. Molecular dispersion energy parameters for alkali and halide ions in aqueous solution. J. Chem. Phys. 2014, 140, 044504. (30) Fuentes-Azcatl, R.; Barbosa, M. C. Sodium chloride, NaCl/: New force field. J. Phys. Chem. B 2016, 120, 2460–2470. (31) Perez Sirkin, Y. A.; Factorovich, M. H.; Molinero, V.; Scherlis, D. A. Vapor pressure of aqueous solutions of electrolytes reproduced with coarse-grained models without electrostatics. J. Chem. Theory Comput. 2016, 12, 2942–2949. (32) Figueirido, F.; Del Buono, G. S.; Levy, R. M. On finite-size corrections to the free energy of ionic hydration. J. Phys. Chem. B 1997, 101, 5622–5623. 25



(33) Hummer, G.; Pratt, L. R.; Garc´ıa, A. E. Ion sizes and finite-size corrections for ionicsolvation free energies. J. Chem. Phys. 1997, 107, 9275–9277. (34) Sakane, S.; Ashbaugh, H. S.; Wood, R. H. Continuum corrections to the polarization and thermodynamic properties of ewald sum simulations for ions and ion pairs at infinite dilution. J. Phys. Chem. B 1998, 102, 5673–5682. (35) H¨ unenberger, P. H.; McCammon, J. A. Ewald artifacts in computer simulations of ionic solvation and ion-ion interaction: A continuum electrostatics study. J. Chem. Phys. 1999, 110, 1856–1872. (36) Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; H¨ unenberger, P. H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing latticesum methods: An accurate correction scheme for electrostatic finite-size effects. J. Chem. Phys. 2013, 139, 184103. (37) Misin, M.; Fedorov, M. V.; Palmer, D. S. Hydration free energies of molecular ions from theory and simulation. J. Phys. Chem. B 2016, 120, 975–983. (38) Malasics, A.; Boda, D. An efficient iterative grand canonical Monte Carlo algorithm to determine individual ionic chemical potentials in electrolytes. J. Chem. Phys. 2010, 132, 244103. (39) Joung, I. S.; Cheatham, T. E. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B 2008, 112, 9020–9041. (40) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269–6271. (41) Joung, I. S.; Cheatham, T. E. Molecular dynamics simulations of the dynamic and

26


Page 26 of 31

Page 27 of 31 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


energetic properties of alkali and halide ions using water-model-specific ion parameters. J. Phys. Chem. B 2009, 113, 13279–13290. (42) Mouˇcka, F.; Nezbeda, I.; Smith, W. R. Chemical potentials, activity coefficients, and solubility in aqueous NaCl solutions: Prediction by polarizable force fields. J. Chem. Theory Comput. 2015, 11, 1756–1764. (43) Aragones, J. L.; Sanz, E.; Vega, C. Solubility of NaCl in water by molecular simulation revisited. J. Chem. Phys. 2012, 136, 244508. (44) Smith, W. R.; Mouˇcka, F.; Nezbeda, I. Osmotic pressure of aqueous electrolyte solutions via molecular simulations of chemical potentials: Application to NaCl. Fluid Phase Equilib. 2016, 407, 76–83. (45) Kobayashi, K.; Liang, Y.; Sakka, T.; Matsuoka, T. Molecular dynamics study of saltsolution interface: Solubility and surface charge of salt in water. J. Chem. Phys. 2014, 140, 144705. (46) Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 1976, 22, 245–268. (47) Shirts, M. R.; Bair, E.; Hooker, G.; Pande, V. S. Equilibrium free energies from nonequilibrium measurements using maximum-likelihood methods. Phys. Rev. Lett. 2003, 91, 140601. (48) Shirts, M. R.; Pande, V. S. Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration. J. Chem. Phys. 2005, 122, 144107. (49) Paliwal, H.; Shirts, M. R. A benchmark test set for alchemical free energy transformations and its use to quantify error in common free energy methods. J. Chem. Theory Comput. 2011, 7, 4115–4134. 27



(50) Chase, M. W. NIST-JANAF Thermochemical Tables. J. Phys. Chem. Ref. Data, Monograph No. 9 1998, American Chemical Society and American Institute of Physics. (51) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19–25. (52) Shah, J. K.; Maginn, E. J. A general and efficient Monte Carlo method for sampling intramolecular degrees of freedom of branched and cyclic molecules. J. Chem. Phys. 2011, 135, 134121. (53) Davies, C. W. The extent of dissociation of salts in water. Part VIII. An equation for the mean ionic activity coefficient of an electrolyte in water, and a revision of the dissociation constants of some sulphates. J. Chem. Soc. 1938, 2093–2098. (54) Hamer, W. J.; Wu, Y.-C. Osmotic coefficients and mean activity coefficients of uniunivalent electrolytes in water at 25C. J. Phys. Chem. Ref. Data 1972, 1, 1047–1100. (55) Lenart, P. J.; Jusufi, A.; Panagiotopoulos, A. Z. Effective potentials for 1:1 electrolyte solutions incorporating dielectric saturation and repulsive hydration. J. Chem. Phys. 2007, 126 . (56) Vincze, J.; Valiskó, M.; Boda, D. The nonmonotonic concentration dependence of the mean activity coefficient of electrolytes is a result of a balance between solvation and ion-ion correlations. J. Chem. Phys. 2010, 133, 154507. (57) Valiskó, M.; Boda, D. The effect of concentration- and temperature-dependent dielectric constant on the activity coefficient of NaCl electrolyte solutions. J. Chem. Phys. 2014, 140, 234508. (58) Valiskó, M.; Boda, D. Unraveling the behavior of the individual ionic activity coeffi-

28


Page 28 of 31

Page 29 of 31 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


cients on the basis of the balance of ion-ion and ion-water interactions. J. Phys. Chem. B 2015, 119, 1546–1557. (59) Valiskó, M.; Boda, D. Activity coefficients of individual ions in LaCl3 from the II+IW theory. Mol. Phys. 2017, 115, 1245–1252. (60) Wagman, D. D.; Evans, W. H.; Parker, V. B.; Schumm, R. H.; Halow, I.; Bailey, S. M.; Churney, K. L.; Nuttall, R. L. The NBS tables of chemical thermodynamic properties. J. Phys. Chem. Ref. Data 1982, 11, Suppl. 2.

29



TOC Graphic

lnγ

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.1 Infinite system size 0.0 Fixed system size 0.1 Varying system size 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.5 1.0 1.5 2.0 2.5

M 1/2

30


Page 30 of 31

lnγ

0.1ofJournal PageThe 31 31 of Physical Chemistry Infinite system size 0.0 Fixed system size 0.1 Varying system size 0.2 1 0.3 2 0.4 3 4 0.5 5 0.6 6 0.7 ACS Paragon Plus Environment 0.0 0.5 1.0 1.5 2.0 7 1/2 8 9

M

2.5

System-Size Dependence of Electrolyte Activity Coefficients in

Recommend Documents