Role of Water Model on Ion Dissociation at Ambient Conditions
RRole of Water Model on Ion Dissociation at Ambient Conditions
Alec Wills
1, 2, a) and Marivi Fernández-Serra
1, 2, b) Physics and Astronomy Department, Stony Brook University. Stony Brook, New York 11794-3800,United States Institute for Advanced Computational Science, Stony Brook, New York 11794-3800, United States (Dated: 29 January 2021)
We study ion pair dissociation in water at ambient conditions using a combination of classical and ab initio approaches.The goal of this study is to disentangle the sources of discrepancy observed in computed potentials of mean force. Inparticular we aim to understand why some models favor the stability of solvent-separated ion pairs versus contact ionpairs. We found that some observed differences can be explained by non-converged simulation parameters. However,we also unveil that for some models, small changes in the solution density can have significant effects on modifyingthe equilibrium balance between the two configurations. We conclude that the thermodynamic stability of contact andsolvent-separated ion pairs is very sensitive to the dielectric properties of the underlying simulation model. In general,classical models are very robust on providing a similar estimation of the contact ion pair stability, while this is muchmore variable in density functional theory-based models. The barrier to transition from solvent-separated to contaction pair is fundamentally dependent on the balance between electrostatic potential energy and entropy. This reflectsthe importance of water intra- and inter-molecular polarizability in obtaining an accurate description of the screenedion-ion interactions.
I. INTRODUCTION
Electrolyte solutions are ubiquitous in nature, and are ofgreat importance in many biological and industrial processes.As such, it is important to understand the nature of ion sol-vation and the effects that the ions have on the solvent. Inso doing, more effective and efficient applications can makeuse of accurate solvation knowledge to improve a variety ofalready-existing technologies, including drug delivery mech-anisms and crystallization methods.Despite their prevalence, it is not completely understoodthe role that ions have in regards to structural changes to sol-vent particles, in particular water. X-ray spectroscopic re-sults show that the alkali cations are generally chaotropic(structure-breaking), having effects on water similar to thatof increasing the solution temperature. Neutron diffractionstudies similarly show chaotropic effects, where the ions actby collapsing the second hydration shell inwards as pres-surizing the solution would do.
However, vibrational re-laxation studies show increased mode lifetimes for Li + andNa + cations, indicating some level of cosmotropic (structure-making) effects for smaller alkali cations. A better under-standing of the explicit solvation process will be beneficial inexplaining the apparent discrepancies between experimentalinterpretations, and yield insight into what exactly structure“making" and “breaking" entails.As the solvent is a critical component in any electrolyte so-lution, the choice of model for water alone is also a factor thatwill influence any simulated dynamics. Classical, empiricallyderived potentials for water molecules have been shown to beinadequate in reproducing experimental results for nonpolar,charged solute molecules, and the choice of charge localiza-tion in the model results in affects orientational dynamics in a) Electronic mail: [email protected] b) Electronic mail: [email protected] the solvation process. X-ray and neutron diffraction studiesdisagree on the Na − O distance of hydrated sodium ions by0.1 Å, and ab initio simulations have a wide range of predic-tions for this value that typically predict values an additional0.1 Å longer than experimental results. Coordination num-ber has an effect on these predicted distances, from whichionic radii are estimated. Charge localization in classical wa-ter models will affect stable configuration geometry, therebyaffecting coordination number. The inter-connectivity of allthese model dependent processes lead to an ever more com-plex web of model choices and parameters that prevents re-producibility in simulations and can obscure physical results.As it is the de facto method for examining system char-acteristics along a variable reaction coordinate, potentials ofmean force (PMFs) are widely used. In particular, NaCl isfrequently used as the model system due to the simplicity(monovalent constituent ions) and prevalence of the salt andits components in biological, industrial, and civilian usage.For NaCl PMFs in particular, there are three main featuresthat characterize the solvation process, shown in Fig. 1. Thecontact ion pair (CIP, shown in Fig. 1a), is where the twoions are Coulombically bound and are surrounded by solventmolecules. The solvent-separated ion pair (SSIP, shown inFig. 1b), is a stable state of solvation wherein the ions share acoordination shell of solvent between them. A transition bar-rier exists between these two states, typically attributed to thework required to separate (bring together) the ions and insert(remove) the intermediate solvent molecules. The height ofthis barrier allows us to characterize the relative stability ofthe separate solvation states, as we do in Tab. I.There are a variety of studies that have been done to re-construct the NaCl PMF with the reaction coordinate as theinter-ionic distance in different classical force fields.
Abinitio treatments of systems vary, ranging from Car-Parrinellomolecular dynamics (CPMD) treatments of the nuclearmotions, to QM/EFP dynamics where fixed geometry frag-ments are assigned electrostatic parameters based on ab ini-tio simulations, to a mix of both BOMD and CPMD for dif- a r X i v : . [ phy s i c s . c h e m - ph ] J a n (a)(b)(c) FIG. 1. (a)
The contact-ion pair (CIP) solvation state. Here, the ionshave no intermediate buffer molecules between them. (b)
Thesolvent-separated ion pair (SSIP) state, where the ions share asolvent hydration shell. (c)
The dissociated ion state, where the ionseach have their own hydration shells. ferent XC functionals in the same study. Two-dimensionalPMFs have also been examined, using coordination num-ber as the second reaction coordinate.
Classical force-field generated PMFs typically predict more stable CIP statesthan the SSIP state that precedes dissociation (shown inFig. 1c), although some conditions show the SSIP to be morefavorable. The results are not as clear for ab initio simula- tions, with results depending on the choice of XC functional,amongst other parameters.
Beyond the energetics, structural information is importantto understand the physical processes that occur during the sol-vation process. Initial neutron diffraction analysis led to re-sults indicating minimal affects beyond the first ion hydrationshell, while later results showed that ion effects extend be-yond the first solvation shell. Even further results suggestedthat cosmo/chaotropic effects have minimal bearing on molec-ular interactions between the solute and solvent, despite thechanges in solvent structure. X-ray absorption studies haveindicated cationic destabilization of the hydrogen bond net-work indicative of extended cationic effects on the solvent,while anionic effects from Cl − indicated insignificant struc-tural effects. However, GGA-level ab initio simulationshave suggested that Na + affects the solvent only locally (notbeyond the first solvation shell), while Cl − has longer-rangeeffects, directly contrasting with interpretations of experi-mental results. Moreover, nuclear quantum effects have beenshown to be an important inclusion to yield more accurateCl – solvation structures, and ionic polarizability also playsa significant role in whether solvation shells have the correctamount level of structure. In our study, we examine the effects of solution densityon the solvation state stability, to better understand how thesolvent thermodynamic properties, which are strongly depen-dent on the underlying solvent model, affect the solvationprocess. We address this question in a two-step approach.Firstly, we investigate the PMFs obtained by different clas-sical models for water. We evaluate separately the differencesthat can be associated to convergence and methodological ap-proaches from those intrinsic to the model used to describethe water-water and water-ion interactions. Then, using con-vergence criteria established within the first half of the work,we study how a first principles treatment of the interactionsmodifies what can be learned from purely classical models.This is done within the formalism of density functional theory(DFT), using two different approximations to the exchangeand correlation (XC) functional. Our choices are motivatedby the need to compare to previously published results, aim-ing to make a fair assessment of similar, but not identicalmethodological approaches. However our choices are funda-mentally driven because these two different ab initio modelspredict a rather distinct structural and dynamical descriptionof liquid water. Hence, we can evaluate whether a modelwhich predicts a low density liquid-type order of liquid wa-ter at normal conditions like PBE is more or less sensitiveto ionic chaotropic/cosmotropic effects than a model whichfavors high density liquid-order type, like the van der Waals(vdW) XC functional vdW-cx. Finally in the last part of the manuscript we discuss thefundamental differences observed between classical and DFTmodels. In the case of classical models, we analyze the en-tropic contributions to the free energy of solvation. We showthat the contact ion pair potential depth is hardly dependent onthe classical model for water or in small changes on the sol-vent density when using the same ion parameters. Howeverthe solvent-separated ion pair is, mostly due to small entropicdifferences between the models.This is not the case for DFT simulations with different XCpotentials. Here both CIP and SSIP potential depths dependon the model and solvent conditions. We argue that the rea-sons are both enthalpic and entropic and fundamentally drivenby changes in the polarizability and dielectric properties of thesolvent water molecules.
II. METHODOLOGY d Na Cl ( Å ) V ( e V ) NPTNVTNVE
FIG. 2. Testing the convergence of the potential of mean force indifferent simulation ensembles with TIP4P water and the defaultOPLS-AA ion parameters in
GROMACS . The agreement between thedifferent ensembles is almost exact, with a CIP depth of ∼
125 meV.Error estimates come from the bootstrapping techniquesimplemented in g-wham . Temperatures in the Boltzmann factors forthe NVE trajectories are the average temperatures of the simulation.
A. Potential of Mean Force
The potential of mean force (PMF) between two particles U ( r ) is defined as: U ( r ) = − kT ln g ( r ) , (1)where g ( r ) is the inter-particle radial distribution function. Insimulations, the evaluation of g ( r ) for charged particles in animplicit polar solvent like water is the computational bottle-neck for calculating U ( r ) , because infinitely long trajecto-ries are needed to completely explore the radial configura-tion space. In particular, electrostatic interactions make thehydration environments of some ionic species long-lasting, requiring biasing methods to force the exploration. One pop-ular method in literature is to use umbrella sampling and re-construct the PMF by weighted-histogram analysis. Anotherpopular alternative to this procedure is to directly evaluate thePMF by integrating constraint forces between the two par-ticles, along a discrete grid of interparticle distances. Both d Na Cl ( Å ) V ( e V ) WHAMDirect Constraints
FIG. 3. Testing the PMF convergence through different samplingmethods, here shown with umbrella sampling and directlyintegrating holonomic constraints. The exact agreement allows formore controlled sampling of the configuration space in our ab initio simulations. The solution used the TIP4P water model and defaultOPLS-AA ion parameters from
GROMACS trajectories. methods map the infinitely long trajectory problem into a fi-nite set of finite time molecular dynamic simulations of con-strained particles.When sampling with constraints, they must be accountedfor in the derivation of the free energy of separation betweenthe particles. Drawing from the statistical mechanics identityrelating the Helmholtz free energy F to the canonical partitionfunction Z F = − kT ln Z , we can suppose a different free energy W is related to a con-strained probability distribution P by W ( r o ) = − kT ln P r o ( r ) . (2)Here, P r o ( r ) is the conditional statistical average given by P r o ( r ) = h δ ( | r | − r o ) i = Z Z d N r e − β Φ ( { r i } ) δ ( | r − r | − r o ) . (3)This new definition is necessary because in order to use dif-ferent methods of PMF calculation, we must force the systeminto thermodynamically unfavorable configurations through P . We note that the definition in Eq. (3) implies that P r o ( r ) isa constrained partition function, measuring the likelihood ofthe interionic distance r taking a value of r o . Constraints areassumed to only depend on particle positions, not momenta,hence the integration over only position-space.In following through with relating the constrained free en-ergy W ( r ) to the PMF U ( r ) , the constraints bring about a cor-rection factor 2 kT / r that must be used to amend the averageforces, yielding dU ( r ) dr = − f ( r ) + kTr , (4)which is integrated to yield the PMF. This correction factorcan be interpreted as an energy correction added back on tothe free energy W ( r ) , arising from the constraints reducing theoverall entropy of the available phase space of the system. Of further note is that this correction factor need only be ap-plied to our ab initio simulations, where we sample the phys-ical forces on the ions (i.e. - ∂ r Φ ( r ) ). Since either alterna-tive PMF construction methods are used in the classical sim-ulations (e.g. the WHAM method), or holonomic constraintforces are used to integrate for the PMF, the correction fac-tor is not needed, as these methods directly yield the freeenergy. III. CONVERGENCE TESTING
When initially determining production simulation param-eters, the trade-off between efficiency and accuracy is key.In particular, it is not currently feasible to simulate ab ini-tio systems at the same speed (and therefore, total trajectorytime) as classical force fields would allow for the same sys-tem. For this reason, we examined the simulation ensemble,sampling method, distance discretization, time step, and sim-ulation length to try to find an appropriate balance betweencomputationally efficient ab initio simulations and accuratepotential of mean force data.Additionally, there is ambiguity in choice of simulation sys-tem physical parameters, most notably the number of watersand size of the simulation box. In pure water simulations ofconstant-volume, one typically chooses a box length from aconstant-pressure equilibration that reproduces the bulk char-acteristics of the water system, for example the density of wa-ter at which temperature you equilibrate. For an infinite bulksystem, introduction of ions into the solvent will have vanish-ing effect on this choice as a large bulk is overwhelmingly in-dependent of any perturbing interactions from the solute. Thisis not so for smaller systems.As ab initio molecular dynamics are much more compu-tationally intensive than classical force field dynamics, andsince the reconstruction of a PMF often increases the expenseby an at least an order of magnitude, we are forced to limitour system sizes to reasonable levels at the expense of bulkstatistics. Since we are limited in size, it is not clear whatphysical characteristic we should choose in determining boxsize. Moreover, the density of liquid water and ice is knownto be dependent on the concentration of dissolved salts, sochoice of simulation parameters might then affect extrapola-tion to bulk characteristics. For this study we chose to solvatea single NaCl molecule in 96 H O molecules. We use a boxsize of 14.373 Å for the ab initio simulations, chosen such thatthe number density of the system was close to that of liquidwater, n = .
033 Å − , corresponding to a solution concentra-tion of 0.56 M. With this choice, we can investigate the localsolvent changes near the ion interface while having the rela-tively limited bulk still approximate the infinitely dilute limit.Additionally, to investigate the effects of solution density onsolvation state stability, we choose another cubic box length of 14.725 Å yielding a number density of 0 . − , slightlygreater than that of ice, corresponding to a concentration of0.52 M. A. Ensemble
For ensemble testing, a single NaCl pair was solvated with96 TIP4P H O molecules using the
GROMACS simulationsuite. With a timestep of 0.5 femtoseconds, each NVE andNVT trajectory was underwent a 25 picosecond NVT equili-bration, while the NPT trajectory equilibrated with NPT con-ditions. Each trajectory then underwent a production run inthe corresponding ensemble for 250 picoseconds. For anygiven ionic distance window, there were eight initial randomseeds given to different equilibrations, resulting in 2 nanosec-onds of total trajectory time for each window.Since we are less limited by computational expense whenusing classical force fields, we first use umbrella sampling toreconstruct the PMF. Each restraining potential was centeredon a reference ion distance starting at 2 Å and ending at 6 Å, insteps of 0.1 Å. With our harmonically constrained windows,the WHAM equations are used to reconstruct the free en-ergy of solvation as the ions are separated. In particular, boththe the g-wham method implemented in GROMACS and theexternal wham program were used and yielded consistent re-sults. Error estimates for WHAM generated PMFs come fromthe bootstrapping method built into g-wham , and are shown aswindows around the average value lines in Fig. 2.We note the almost exact agreement of all the results. Withthe agreement of the ensembles to within error, to optimizesimulation efficiency in the ab initio simulations we choose touse the NVE ensemble in our ab initio simulations. This hasfurther benefit by allowing us to sample the true dynamics ofthe systems we simulate. B. Sampling Method
The WHAM method requires convergence of each sampledumbrella window histogram to properly converge the PMF.This means that windows centered around inter-ion distancesnearing the completely solvated range will require much moretime to accurately converge, as the configuration states inthese windows are more or less equally likely. Rather thanrisk this computational bottleneck, we instead choose to di-rectly constraint the distances in our ab initio simulations. Theresults comparing the WHAM and direct integration methodsare plotted in Fig. 3. Ion distances are constrained with theSHAKE algorithm, and the holonomic constraint forces aredirectly integrated to yield the same result as generated withumbrella sampling. C. Distance Discretization
For direct constraints PMF simulations, me must choose adiscretization of the mesh of inter ionic distances on which the d Na Cl ( Å ) V ( e V ) x = 0.1 Åx = 0.2 Åx = 0.3 Åx = 0.4 Å FIG. 4. PMF convergence testing parameterized by changing stepsize in the discretization of distance between the ions, using afixed-distance restraint in a TIP4P solution with default OPLS-AAion parameters in
GROMACS . Increments of 0.3 Å and 0.4 Å (greenand red, respectively) do not yield converged results, whileincrements of 0.2 Å and 0.1 Å (yellow and blue, respectively)converge to consistency.
PMF is integrated. A choice must be made between compu-tational cost and accuracy, which can significantly affect theintegrated potential results. Fig. 4 shows results for the di-rectly constrained NVE PMF from 2 to 6 Å as parameterizedby the distance mesh separating the two points.The results exactly agree whether using distance steps of0.1 Å or 0.2 Å . Some slight disagreement is seen when thestep is increased to 0.3 Å and increasing to 0.4 Å completelyruins any convergence. For the purposes of our ab initio sim-ulations, we chose a distance step of 0.2 Å to have a moreregular grid of sample points while keeping the number ofnecessary simulations at a reasonable level.
D. Simulation Length
With above parameters determined, we compared the totallength of necessary trajectory sampling needed to convergeto a consistent potential. For the given trajectories, config-urations were randomly sampled (without replacement) withvarying sizes to reconstruct the PMF. The results are shown inFig. 5. We find that randomly sampling only 5-10% of the fullsample size yields a well-converged and consistent PMF. Witha full sample size of 40,000 snapshots at each reference dis-tance, a subset of 4,000 configurations is sufficient to exactlyagree with the full trajectory result. Moreover, 2,000 config-urations agrees with minimal deviation. With less sampling,the results diverge more starkly and no longer agree within areasonable tolerance. For this reason we choose our ab initio sample size to be 20,000 snapshots at each reference distance,which, with a timestep of 0.5 femtoseconds, yields 10 picosec-ond trajectories across the 21 distances at which we constrainthe ions. d Na Cl ( Å ) V ( e V ) N = 200 N = 400 N = 1000 N = 2000 N = 4000N = 40,000 FIG. 5. Examining the convergence of the directly integrated NVEPMF through varying the total sample size of snapshots fromtrajectories of TIP4P water with OPLS-AA ion parameters in
GROMACS . The PMFs generated using 10% and 5% (purple and red,respectively) are consistent with the full trajectory’s PMF.Decreasing the subset of configurations used beyond these valuesshows deviations in the result. E. Ab Initio
Simulation Convergence
In Born-Oppenheimer DFT ab initio molecular dynamics,the time step of the simulation is determined by the need to ac-curately sample the shortest vibrational modes of the system.As classical models are usually rigid, classical molecular dy-namic simulations can afford a longer time step than quantummechanical simulations. Having a time step too large has beenshown to affect the vibrational structure of the bulk system,thus changing the effective hydrogen bond strength betweenwater molecules.
In particular, increasing the simulationtimestep has been shown to redshift the O − H stretching and,to a lesser extent, bending modes.In the Supplemental Information, Fig. S1 shows the vibra-tional spectra of an ab initio simulation of the same box of96 water molecules using the vdW-DF-cx exchange and cor-relation (XC) functional with timesteps of 0.5 and 1.5 fs, forthe same total amount of steps, as well as the peak resonancesof intermittent time steps from 1.0 to 1.5 fs. The redshift ofthe vibrational modes is clearly evident as the timestep in-creases, most notably for the O − H stretching modes. Red-shifted stretching mode peaks correspond to water with aweaker hydrogen bond network, and thus less structure be-tween water molecules. This is shown in Fig. S2, where weplot the O − O radial distribution functions for pure water forboth 0.5 and 1.5 fs simulation timesteps against experimentalresults. A weaker bulk system could have significant effects on anyreconstructed PMF in that system. In this case, as the ionsare pulled apart from each other, there is less strength in thefirst hydration shell to resist the separation, which would ef-fectively decrease the stability of the contact ion pair state.For this reason, we choose an ab initio time step of 0.5 fem-toseconds, to accurately simulate the hydrogen bond networkof the solvent. We believe that longer timestep sizes can sig-nificantly compromise the comparison of results obtainedwith the same model.We further note that there are many convergence propertiesof ab initio simulations which are basis-set dependent. In par-ticular, the ab initio framework utilized in this work uses nu-merical atomic orbitals (NAOs), which, while more computa-tionally efficient, lack clear systematic ways of improvementfor a fixed basis size that plane-wave methods may employ.Such convergence tests have been explored in great detail, inparticular for the liquid water bulk in our current system ofinterest. We do not discuss such properties in this work.
F. Single-/Multi-dimensional Reaction Coordinates
Previous works have shown that multidimensional re-action axes are necessary to accurately probe the tran-sition dynamics and dynamic correction factors to rateconstants. In particular, in addition to the ionic distance re-action coordinate, a coordinate involving the solvent structurenearby the ions might also be imposed. One way to achievean accurate sampling of solvent configurations is to apply aharmonic constraint to a coordination number around a givenion – in our case, around the Na + ion. To compare effects mul-tidimensional sampling have on our calculated PMF, we usedthe coordNum reaction coordinate from the COLVAR pack-age implemented in LAMMPS. This reaction coordinate isdefined as C ( G , G ) = ∑ i ∈ G ∑ j ∈ G − (cid:16) r ij d (cid:17) n − (cid:16) r ij d (cid:17) m (5)where G and G are atoms groups composed of the Na atomand O atoms in SPC/E water, respectively. Ion parameterswere chosen to match that of the solvent model. The cutoffdistance d was chosen to be 3.0 Å, just before the end of thefirst solvation shell of the Na − O radial distribution functionas seen in Fig. S9b. The exponents n and m maintained theirdefault values of 6 and 12, respectively.For the two dimensional FES reconstructed with umbrellasampling, the distance windows had centers from 2.0 to 6.0Å, in steps of 0.2 Å. The coordination number reaction coor-dinate had window centers from 3.0 to 8.0 in increments of0.25. The full two-dimensional free-energy heatmap is shownin Fig. S3. Of note is the much deeper CIP region than theSSIP, in agreement with our one-dimensional results. TheSSIP shows to be less stable than the CIP by a value on theorder of 60 meV.Moreover, integrating out the extra degree of freedom tocompare one-dimensional and two-dimensional PMFs showsno significant difference in the free energy change during sol-vation, as seen in Fig. 6. Indeed, while forcing the solventconfiguration sampling to occur through use of a harmonicpotential might speed up the process, a sufficient number ofsnapshots at a given reference distance should sample the sec-ond reaction coordinate space reasonably well, and without d Na Cl ( Å ) V ( e V ) CV: r CV: r , CN FIG. 6. A comparison between the two-dimensional PMF afterintegrating out the coordination number reaction coordinate (dotted)and the one-dimensional PMF using just the interionic distance. Thecoordination number degree of freedom is integrated out usingBoltzmann factor weighting across the sampled coordinationnumbers at a given distance. introducing potentially unphyisical configurations by forcingtoo many waters to be close to a given ion. For this reasoning,we choose to focus on the one-dimensional PMFs generatedthrough the interionic distance coordinate, to maximize effi-ciency in simulation expense while maintaining accuracy inour results.
G. Parameter Convergence Summary
So far we have evaluated how the PMF depends on primarysimulation parameters, independently of the chosen water-water and water-ion interaction model. The previous classi-cal and ab initio convergence tests give us insights as to whatparameters are appropriate for our studied system. We usesymplectic integrators for both classical and ab initio simu-lations, with no temperature or pressure couplings, so that wecan accurately simulate the true dynamics of the systems in anNVE ensemble. We choose to directly constraint the ab initio interionic distances from 2 to 6 Å in steps of 0.2 Å, to havebetter control over the mean force sampling convergence. Weuse a 0.5 fs timestep for our ab initio simulations and a to-tal time of 10 ps at each distance to be able to efficiently andaccurately reconstruct the PMF.
IV. RESULTS: MODEL AND DENSITY-DEPENDENTPOTENTIALS OF MEAN FORCE
The goal of our study is to understand the physical ori-gin for the observed disparity of relative stabilities of CIPand SSIP obtained with different models.
In particular wewant to disentangle the effects from water-water and ion-waterinteractions. To do so we choose to evaluate separately clas- d Na Cl ( Å ) V ( e V ) FIG. 7. The PMFs from the trajectories generated from directlyconstraining the ion distances in
GROMACS . There is seeminglyminimal effect of molarity on the classical CIP, as evidenced by the0.56 M (red) and 0.52 M (blue) having the same overall structure,with a minima of comparable depths. The SSIP region is slightlyaffected, however, as the higher density predicts a marginally morestable state than the lower density simulation. Error is estimated byrandomly partitioning the constraint forces and finding the standarddeviation of the partition means. sical and ab initio models using two different molarities, bothwithin the highly diluted limit. A summary table with PMFproperties for all our simulations is presented in Table I. Forboth classical and ab initio setups, we use two cubic boxes ofNaCl+96 H O, one with a side length of L = .
373 Å and onewith L = .
725 Å, to reproduce the number densities previ-ously mentioned. These lengths yield solution concentrationsof 0.56 M and 0.52 M, respectively. For the two different clas-sical models we also obtain PMFs with two different sets ofLennard-Jones parameters to describe the ion-water interac-tions.
A. Classical PMFs
At the classical level, we use OPLS-AA force-field pa-rameters included in GROMACS , with the ions solvated inTIP4P water.With the previously mentioned system compositions, we re-construct the NVE ensemble PMF of NaCl as shown in Fig. 7.We note the relatively equal depths in both conditions for theCIP state , on the order of 125 meV, both of which are muchmore stable than the SSIP state as we would expect from theclassical force field. We further note that the SSIP state in the0.56 M solution is slightly more stable than that of the 0.52 Msolution. As noted by a previous study, increased pressurehas the effect of stabilizing the SSIP through entropic contri-butions. Therefore it is to be expected that the smaller systemhas a more stabilized SSIP.Fig. 8 shows the PMFs for the 0.56 M simulation box butwith a combination of different water models and ion param- d Na Cl ( Å ) V ( e V ) FIG. 8. PMFs generated with different parameter choices in 0.56 Msolution. Note the enhanced stability of the SPC/E (blue) SSIPstates when compared to those in Fig. 7. Moreover, changing fromthe OPLS-AA parameters in the SPC/E solvent to the JC parameters(green) further enhances the stability of the SSIP, and affects theCIP location as well. eters. The OPLS-AA force field parameters are used to com-pare the effects of changing the water model from TIP4P toSPC/E. As shown in Table I, while for both models the CIPis the favored minimum, the SSIP is stabilized by 20 meV inSPC/E. This stabilization must be due to the changes in thewater model. We also evaluate how the change of ion parame-ters modify the stabilities. For the SPC/E simulations two dif-ferent sets of ion parameters are used – the OPLS-AA parame-ters and parameters optimized by Joung and Cheatham (JC). The JC parameters further stabilize the SSIP, while also shift-ing the location of the CIP. Both changes are expected becausethe JC parameters strengthen the water-ion interactions andincrease the size of the ions. These results for the classicaltransition barriers are summarized in Table I. We also simu-lated a system with the SPC water model and OPLS-AA ionparameters, shown in Fig. S4. The lack of drastic differencebetween the SPC and SPC/E simulations with the OPLS-AAion parameters can be attributed to the minimal difference be-tween the water models, as opposed to the more substantialchange in going from TIP4P to either of the SPC models. B. Ab initio
PMFs
The DFT-based PMFs are computed using the
SIESTA code, which uses numerical atomic orbitals as basis sets forvalence electrons and norm-conserving pseudopotentials. For all atoms, we used double- ζ +P basis sets. Input fileswith all used parameters are provided in the SI. In all ofthe ab initio simulations, we do 10 picosend NVE-ensembleruns with a timestep of 0.5 femtoseconds, seeded with coordi-nates from the corresponding box’s classical simulations. Foreach box size, we do a set of simulations using the PBE exchange-correlation functional in the generalized-gradient TABLE I. State transition barriers for classical and ab initio
PMFs.Model Concentration (M) V TS − V CIP (meV) V TS − V SSIP (meV) V SSIP − V CIP (meV) log K a log K e TIP4P+OPLS-AA 0.56 130 . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± − . ±
19 3.47 0.900.52 77 . ± . ± − ±
17 3.51 2.64vdW-DF-cx 0.56 116 . ± . . ± − ±
20 3.21 1.400.52 106 . ± . ± − ±
17 3.36 0.95PBE a . ± . − . ± . b . ± . − . ± . a NPT, 400 K, 1 atm b NPT, 1000 K, 11 GPa d Na Cl ( Å ) V ( e V ) (a) PBE d Na Cl ( Å ) V ( e V ) (b) vdW-DF-cx FIG. 9. (a)
The potentials of mean force from the fixed distance constraint method implemented in SIESTA on the PBE trajectories. In thiscase, the lower concentration PMF (blue) is the only one that predicts a markedly more stable SSIP than CIP. (b)
The potentials of mean forcefrom the fixed distance constraint method implemented in SIESTA on the vdW-DF-cx trajectories. The inclusion of van der Waals correctionsmakes the simulation less sensitive to density perturbations. Error is estimated by randomly partitioning the constraint forces and finding thestandard deviation of the partition means. approximation, and another set including van der Waals inter-actions through the vdW-DF-cx XC parameterization. Wemake this choice because we want to evaluate the differencesin PMF using two models which clearly produce two very dif-ferent types of liquid water structure. PBE is well known toproduce an overstructured, low-density-type liquid, andits equilibrium density is ≈ . g / cm . vdW-DF-cx on theother hand produces a high density liquid, with an equilib-rium density of ≈ . g / cm .The distance between the Na-Cl ions was constrained from2 to 6 Å, incremented in steps of 0.2 Å, and the PMF wasreconstructed by integrating the mean force projection ontothe vector separating the ions. Due to the fact that we samplethe physical forces on the ab initio ions as opposed to holo-nomic constraint forces, the forces integrated in these SIESTA simulations include a Jacobian correction factor. After inte-grating Eq. 4, the Jacobian factor in the case of a two-atom distance constraint takes the form 2 k B T ln r , with r the dis-tance between ions, and can be interpreted as an entropiccorrection. To seed the ab initio simulations, coordinates from the clas-sical trajectory were chosen as initial configurations for thedynamics at each reference distance. Only the last 15,000timesteps are used in the analysis, to allow the systems timeto equilibrate. The resulting data from the normal and lownumber density simulations are shown in Figs. 9a and 9b,respectively. The transition barrier energies are summarizedalongside the classical results in Table I.Contrary to the classical simulations, the DFT PMFs sta-bilities show a much stronger dependence on the density ofliquid water. In particular, the 0.52 M PBE solution predictsa much more stable SSIP state than CIP, whereas the higherconcentration PBE solution has equally favorable state stabil-ities. The 0.52 M PBE solution is the only result here to have d O O ( Å ) g OO ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M)Soper (2013) (a) Pure Water Solutions d O O ( Å ) g OO ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M)Soper (2013) (b) Ionic Solutions FIG. 10. (a)
The O − O radial distributions functions of pure waterusing TIP4P, vdW-DF-cx, and PBE (green, red, and blue,respectively) with the box size that corresponds to the 0.56 Msolutions with the ions removed. (b)
The therymodynamicallyaveraged O − O radial distributions functions of the ion solutions forthe TIP4P, vdW-DF-cx, and PBE simulations with concentrations of0.56 M. Both figures show experimental results at ambientconditions for comparison. The inclusion of the ions serves toincrease the structure in the ab initio cases. the favored SSIP, as predicted by others.
However, onestudy uses a timestep of 1.5 fs in their simulations, whichwe have shown weakens the hydrogen-bond network, therebyenhancing SSIP stability. All other DFT simulation resultshere predict relatively equal stability for the solvation states,with the vdW-DF-cx simulations predicting a slightly morestable SSIP. Moreover, the vdW inclusive functional yieldsPMFs that seem to be more resistant than PBE to changes inmolarity. It is, however, striking how much more sensitivethe DFT PMFs are to small changes in the simulation densitythan classical results. This can originate from the very differ-ent compressibilities of the DFT pure liquids. To explicitly show the network behavior in the differentdensities and XC functionals, we plot the O − O radial distribu-tion functions for pure water in Fig. 10a. Average thermody-namic values for the systems plotted are given in Table SII for the ionic solutions, and in Table SIII for the pure solutions.We immediately see that the inclusion of the ions serves tostructure both sets of system. However, the effect is similarfor both XC functionals. This indicates that the reason be-hind the high sensitivity of the PBE PMF to small changes indensity of water must be rooted in energetic and/or entropicreasons. We will explore this in the following section.Finally, under the approximation that we are near theinfinitely-dilute regime, we can reconstruct the ion-ion radialdistribution function by inverting Eq. 1 for g ( r ) , which wewill denote as g ∞ ( r ) to emphasize the approximation. Withthis, we can find the ion-ion association constant and solvationstate equilibrium constant. We adapt the referenced work’sequations for the radial cutoffs given in our simulations, yield-ing K a = π Z R D g ∞ ( r ) · r drK e = R R D R TS g ∞ ( r ) · r dr R R TS g ∞ ( r ) · r dr , where K a and K e are the association and equilibrium con-stants, respectively. R TS and R D are the transition state bar-rier radius and maximum simulation constraint radius, respec-tively. Formally, the lower integration limit of 0 to cover theCIP state must be replaced by an artificial cutoff of 2.0 Å,which is the smallest constraint distance we utilize. In oursimulations R D is 6.0 Å, while R TS varies depending on theXC functional and density. The results for the constants gen-erated from our simulations are shown in Table I. In all cases,the K e values predict a higher population of SSIP pairs, evenwhen the CIP is more energetically favored as is the case inthe simulations with OPLS-AA ion parameters. This is likelya result of the extended region over which the SSIP is stable.Through the larger volumetric probability, the ions are morelikely to be found in the SSIP regime.The K a values we calculate differ significantly from em-pirical fits to experimental values at lower solvent densitiesand higher temperatures, but behave consistently in that theassociation constant is expected to increase as the tempera-ture decreases and the ions can more readily overcome ther-mal fluctuations. We further note that our ab initio K a valuesare less than those predicted by a QM/EFP study of the samesystem, but again are consistent in that the QM/EFP studythermostatted their system to 300 K, while ours equilibratedthemselves to a higher temperature, thereby causing increasedthermal fluctuation strength. We note that the K e value of 2.64for the low density PBE simulation predicts an approximateratio of 1:400 of CIP to SSIP solvation state population, inagreement with the steep decrease in CIP stability.As structural changes induced by the ions are of consid-erable important in understanding the solvation process, wedicsuss the relevant results here but relegate a large amount offigures to the SI. Notably, we see in Fig. S9 the effects the ionshave on the neighboring solution through the thermodynami-cally averaged radial distribution functions of the ion-oxygenpairs. The presence of a clear second shell in the Na − O dis-tributions (Fig. S9b) is indicative of slightly extended pertur-0 d Na Cl ( Å ) F ( e V ) TIP4P+OPLS-AASPC/E+OPLS-AA (a) d Na Cl ( Å ) H ( e V ) TIP4P+OPLS-AASPC/E+OPLS-AA (b) d Na Cl ( Å ) T S ( e V ) TIP4P+OPLS-AASPC/E+OPLS-AA (c)
FIG. 11. (a)
The WHAM generated NVT PMFs for OPLS-AA ions solvated in TIP4P (red) and SPC/E (blue) water models. (b)
Thehistogram-binned total potential energies along the reaction coordinate, giving the enthalpic contribution to the free energy. (c)
The entropiccontribution to the free energy, taken as − T ∆ S = F − H . bations on the surrounding solvent molecules, what one mightsee referred to as a “cosmotropic” effect in that it introducesstructure into the surrounding solvent. An analogous effect ispresent in the Cl − O distributions (Fig. S9c), although not asstrongly as is the case with the Na − O pair.Additionally, indications of local cosmotropic effects arepresent in Table SIV, where the first hydration shell radius( r ) for the Na − O RDF is notably pulled in to lower dis-tances. Moreover, the PBE Na − O second coordination shellis markedly more contracted inwards than the vdW-DF-cxsimulations, for both molarities. This indicates that the PBEsolution is more prone to nonlocal cosmotropic perturbationsinduced by the Na + cation. Table SIV further shows that theCl − O RDFs show extended perturbation in solvation struc-ture, in agreement with previous results. The solvation radiifor both the first and second solvation shells ( r and r , re-spectively) are notably shifted to higher distances, indicativeof a weakening of the hydrogen bond network in the vicinityof the Cl – ion – a “chaotropic” response. V. DISCUSSION
From all the previous results two main observations can bemade: (i)
Classical PMFs using the same ion parameters reliablypredict the same potential depth for the contact ion pair. Thesolvent-separated ion pair minimum and depth, on the otherhand, is more dependent on the water model. Small changesin the density of the solution do not alter these results. Thesetwo outcomes are not totally unexpected. The CIP poten-tial depth and position is mostly dependent on the ion sizesand their charges. For rigid, non-polarizable ion models onewould expect them to not change, independently of the waterforce field. The change of the SSIP when the water modelis changed is more complex. The reasons could be entropic,enthalpic or a combination of the two. (ii)
DFT results show that both CIP and SSIP are stronglydependent both on the exchange and correlation approxima-tion and the density of the solution. The latter effect is partic-ularly large for PBE. The fact that the CIP is so sensitive to the XC model indicates that the ion size and polarizability arestrongly sensitive to the XC approximation, a conclusion sup-ported by previous studies. Furthermore, the solvent dielec-tric properties will also be heavily dependent on the XC ap-proximations employed, as well as the density and other ther-modynamic conditions of the solution.
The dependencewith the density indicates that there is an additional contri-bution to this polarizability on the density of the surroundingwater molecules. Hence these results point to the importantrole that the dielectric screening plays in shaping the PMF.Although in this study we do not aim to make an exhaustivestudy of how the dielectric screening of charges influences thePMF, we can gain some information of this by partitioning thefree energy into entropic and enthalpic contributions.In order to do this we can sample the potential energy of theconfigurations that produced the PMF. Because most of thesystem is composed of liquid water, large fluctuations dom-inate the potential energy and it is necessary to sample longtrajectories to observe the ion-ion potential energy contribu-tion. This limits our potential energy study to classical trajec-tories. Results are shown in Fig. 11b for TIP4P and SPC/E,both using the same ion parameters. As the simulations aredone at constant volume, the enthalpy is obtained using theaverage potential energy U . The free energy F is nothing elsebut the potential of mean force, from which we can extractthe entropic component through F − U = − T ∆ S , in a mannersimilar to previous studies. Thermostatting is required forthis decomposition, otherwise we don’t get realistic energeticpartitioning as the ions are separated, thus the partitioned sim-ulations were done in the NVT ensemble.Two interesting observations can me made. Firstly, theSSIP is more stable than the CIP when only the total potentialenergy (enthalpic) contribution is accounted for. This indi-cates the CIP stabilization at room temperature is dominatedby the entropic contribution. Moreover, the entropy is rela-tively constant along larger interionic separation distances, upto the transition barrier region between the SSIP and the CIPstates. Here, the two systems’ entropic contributions increasewith similar slope. Note we plot − T ∆ S in the figure, thusthe entropy is increasing in the regime between transition bar-rier peak inwards to the CIP. This entropic stabilization can1be associated to the entropy increase of the liquid by the re-lease of the water molecule(s) involved in the SSIP sharedsolvation shell. These results, obtained with models that onlydiffer on their electrostatic description of water, indicate thatthe model-dependent dielectric screening of the ions signifi-cantly contributes to the PMF’s differences. We note that thedielectric constant of TIP4P water is ∼
50 while that of SPC/Ewater is ∼ Larger screening effects will soften the en-thalpic contribution to the SSIP. At the same time, a largerdielectric constant will also result in a larger entropy increase,as seen in Fig. 11c, which shows that the SPC/E entropy islarger than that of TIP4P. Hence, this decomposition unveilsthat the stability of the two minima is a balance of oppositeinteractions, and small changes in the model can result in bigdifferences in the PMF. This might be even more significantin polarizable models, explaining why DFT-based simulationsare much more sensitive to small changes in simulation con-ditions and XC functionals.It is not feasible to obtain an equivalent figure for the DFTsimulations, given that one would need nanosecond lengthsimulations to obtain a potential energy plot discernible fromits fluctuations. It is possible to attempt to obtain an estimateof the entropic contribution using local connectivity fluctua-tions through Voronoi tesselations. We present such a studyin the Supplementary Information, but notions of entropy insuch characterizations are not readily connected to the phys-ical entropies we produce above. Nevertheless, the behaviorof local ionic volumes reveal markedly different behavior inclassical and ab initio systems, as shown in Fig. S5. Overall,these results indicate that the entropic effects are very similarfor all models. The leading conclusion for this points to theimportance of polarizability and overall dielectric propertiesof water as the source of the large variations of PMF energiesobserved across the models.
VI. CONCLUSION
We have reconstructed potentials of mean force for NaClsolutions of different molarities using solvent number densi-ties corresponding to normal and almost ice-like water, in boththe classical and ab initio regime. The goal of this study wasto understand the differences observed in PMFs computed us-ing different models, sampling methods, and simulation pa-rameters. Overall we find that when simulation-based param-eters are converged, classical, non-polarizable models pro-duce similar PMFs when using the same ion parameters. Inall classical simulations, the CIP solvation state is more stablethan the SSIP state under the conditions chosen. Moreover,the depth of the CIP free energy minimum is independent ofthe specific water model. The choice of water model does,however, have a greater effect on the SSIP stability. By de-composing the free energy into entropic and enthalpic contri-butions we show that differences between models are likelyrooted in the dielectric screening of the ions, which for non-polarizable models mostly affects the SSIP.On the other hand, ab initio simulations have much morevariability in the solvation state stabilities. In particular, the solvent density has a much more profound impact on the en-ergetically preferred state. These results point to a much largersensitivity of the polarizability and dielectric response in thesemodels to the thermodynamic conditions.A detailed study of how the dielectric screening influencesthe ion pair association and stability will be the subject of afuture research study.
ACKNOWLEDGMENTS
This material is based upon work supported by the Depart-ment of Energy under award number DE-SC0001137. In ad-dition, this research used resources of the National EnergyResearch Scientific Computing Center, a DOE Office of Sci-ence User Facility supported by the Office of Science of theU.S. Department of Energy under Contract No. DE-AC02-05CH11231. We also thank Stony Brook Research Comput-ing and Cyberinfrastructure, and the Institute for AdvancedComputational Science at Stony Brook University for accessto the high-performance SeaWulf computing system, whichwas made possible by a National Science Foundation grantNo. 1531492. I. Waluyo, D. Nordlund, U. Bergmann, D. Schlesinger, L. G. M. Petters-son, and A. Nilsson, “A different view of structure-making and structure-breaking in alkali halide aqueous solutions through x-ray absorption spec-troscopy,” The Journal of Chemical Physics , 244506 (2014). A. K. Soper and K. Weckström, “Ion solvation and water structure inpotassium halide aqueous solutions,” Biophysical Chemistry , 180–191(2006). R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, “Per-turbation of water structure due to monovalent ions in solution,” PhysicalChemistry Chemical Physics , 2959–2967 (2007). M. F. Kropman and H. J. Bakker, “Effect of ions on the vibrational re-laxation of liquid water,” Journal of the American Chemical Society ,9135–9141 (2004). R. C. Remsing, T. T. Duignan, M. D. Baer, G. K. Schenter, C. J. Mundy,and J. D. Weeks, “Water Lone Pair Delocalization in Classical and QuantumDescriptions of the Hydration of Model Ions,” Journal of Physical Chem-istry B , 3519–3527 (2018). J. J. Fifen and N. Agmon, “Ionic radii of hydrated sodium cation fromQTAIM,” The Journal of Chemical Physics , 034304 (2019). I. V. Khavrutskii, J. Dzubiella, and J. A. McCammon, “Computing accu-rate potentials of mean force in electrolyte solutions with the generalizedgradient-augmented harmonic Fourier beads method,” Journal of ChemicalPhysics (2008), 10.1063/1.2825620. J. Timko, D. Bucher, and S. Kuyucak, “Dissociation of NaCl in water fromab initio molecular dynamics simulations,” Journal of Chemical Physics (2010), 10.1063/1.3360310. M. K. Ghosh, S. Re, M. Feig, Y. Sugita, and C. H. Choi, “Interionic Hydra-tion Structures of NaCl in Aqueous Solution: A Combined Study of Quan-tum Mechanical Cluster Calculations and QM/EFP-MD Simulations,” TheJournal of Physical Chemistry B , 289–295 (2013). Y. Yao and Y. Kanai, “Free Energy Profile of NaCl in Water:First-Principles Molecular Dynamics with SCAN and ω B97X-V Ex-change–Correlation Functionals,” Journal of Chemical Theory and Com-putation , 884–893 (2018). S. Roy, M. D. Baer, C. J. Mundy, and G. K. Schenter, “Marcus Theory ofIon-Pairing,” Journal of Chemical Theory and Computation , 3470–3477(2017). R. Car and M. Parrinello, “Unified Approach for Molecular Dynamicsand Density-Functional Theory,” Physical Review Letters , 2471–2474(1985). C. Zhang, F. Giberti, E. Sevgen, J. J. de Pablo, F. Gygi, and G. Galli,“Dissociation of salts in water under pressure,” Nature Communications , 3037 (2020). R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, “Per-turbation of water structure due to monovalent ions in solution,” PhysicalChemistry Chemical Physics , 2959–2967 (2007). R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, “Hydra-tion of sodium, potassium, and chloride ions in solution and the conceptof structure maker/breaker,” Journal of Physical Chemistry B , 13570–13577 (2007). I. Waluyo, C. Huang, D. Nordlund, U. Bergmann, T. M. Weiss, L. G. M.Pettersson, and A. Nilsson, “The structure of water in the hydration shell ofcations from x-ray Raman and small angle x-ray scattering measurements,”The Journal of Chemical Physics , 064513 (2011). A. P. Gaiduk and G. Galli, “Local and Global Effects of Dissolved SodiumChloride on the Structure of Water,” The Journal of Physical ChemistryLetters , 1496–1502 (2017). J. Xu, Z. Sun, C. Zhang, M. DelloStritto, D. Lu, M. L. Klein, and X. Wu,“Importance of nuclear quantum effects on the hydration of chloride ion,”Physical Review Materials , L012801 (2021), arXiv:2009.07304. M. DelloStritto, J. Xu, X. Wu, and M. L. Klein, “Aqueous solvation of thechloride ion revisited with density functional theory: impact of correlationand exchange approximations,” Physical Chemistry Chemical Physics ,10666–10675 (2020). J. Perdew, K. Burke, M. Ernzerhof, “Generalized gradient approximationmade simple,” Phys. Rev. Lett. , 3865–3868 (1996). K. Berland, P. Hyldgaard, “Exchange functional that tests the robustnessof the plasmon description of the van der Waals density functional,” Phys.Rev. B , 1–8 (2014). M. Fritz, M. Fernández-Serra, and J. M. Soler, “Optimization of anexchange-correlation density functional for water,” Journal of ChemicalPhysics (2016), 10.1063/1.4953081, arXiv:1603.07302. J. G. Kirkwood, “Statistical mechanics of fluid mixtures,” The Journal ofChemical Physics , 300–313 (1935). M. Li, Z. Duan, Z. Zhang, C. Zhang, and J. Weare, “The structure, dy-namics and solvation mechanisms of ions in water from long time molecu-lar dynamics simulations: a case study of CaCl 2 (aq) aqueous solutions,”Molecular Physics , 2685–2697 (2008). S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Koll-man, “THE weighted histogram analysis method for free-energy calcula-tions on biomolecules. I. The method,” Journal of Computational Chemistry , 1011–1021 (1992). J.-L. Li, R. Car, C. Tang, and N. S. Wingreen, “Hydrophobic interactionand hydrogen-bond network for a methane pair in liquid water,” Proceed-ings of the National Academy of Sciences , 2626–2630 (2007). G. Ciccotti and M. Ferrario, “Holonomic Constraints: A Case for StatisticalMechanics of Non-Hamiltonian Systems,” Computation , 11 (2018). R. C. Dougherty, “Density of Salt Solutions: Effect of Ions on the ApparentDensity of Water,” The Journal of Physical Chemistry B , 4514–4519(2001). P. Novotný and O. Söhnel, “Densities of Binary Aqueous Solutions of 306Inorganic Substances,” Journal of Chemical and Engineering Data , 49–55 (1988). W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L.Klein, “Comparison of simple potential functions for simulating liquid wa-ter,” The Journal of Chemical Physics , 926–935 (1983). M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, andE. Lindahl, “GROMACS: High performance molecular simulations throughmulti-level parallelism from laptops to supercomputers,” SoftwareX ,19–25 (2015). J. S. Hub, B. L. de Groot, and D. van der Spoel, “g_wham—A FreeWeighted Histogram Analysis Implementation Including Robust Error andAutocorrelation Estimates,” Journal of Chemical Theory and Computation , 3713–3720 (2010). A. Grossfield, “Wham: an implementation of the weighted histogram anal-ysis method,” http://membrane.urmc.rochester.edu/content/wham/, version2.0.9. J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, “Numerical integrationof the cartesian equations of motion of a system with constraints: molecu-lar dynamics of n-alkanes,” Journal of Computational Physics , 327–341 (1977). M. Allesch, E. Schwegler, F. Gygi, and G. Galli, “A first principles simu-lation of rigid water,” Journal of Chemical Physics , 5192–5198 (2004),arXiv:0401267 [cond-mat]. M. Praprotnik and D. Janežiˇc, “Molecular dynamics integration and molec-ular vibrational theory. III. The infrared spectrum of water,” The Journal ofChemical Physics , 174103 (2005). K. Berland and P. Hyldgaard, “Exchange functional that tests the robustnessof the plasmon description of the van der Waals density functional,” Phys-ical Review B - Condensed Matter and Materials Physics , 1–8 (2014),arXiv:arXiv:1309.1756v1. A. K. Soper, “The Radial Distribution Functions of Water as Derived fromRadiation Total Scattering Experiments: Is There Anything We Can Say forSure?” ISRN Physical Chemistry , 1–67 (2013). F. Corsetti, M.-V. Fernández-Serra, J. M. Soler, and E. Artacho, “Optimalfinite-range atomic basis sets for liquid water and ice,” Journal of Physics:Condensed Matter , 435504 (2013). P. L. Geissler, C. Dellago, and D. Chandler, “Kinetic Pathways of Ion PairDissociation in Water,” The Journal of Physical Chemistry B , 3706–3710 (1999). R. G. Mullen, J.-E. Shea, and B. Peters, “Transmission Coefficients,Committors, and Solvent Coordinates in Ion-Pair Dissociation,” Journal ofChemical Theory and Computation , 659–667 (2014). G. Fiorin, M. L. Klein, and J. Hénin, “Using collective variables todrive molecular dynamics simulations,” Molecular Physics , 3345–3362(2013). S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynam-ics,” Journal of Computational Physics , 1–19 (1995), arXiv:nag.2347[10.1002]. I. S. Joung and T. E. Cheatham, “Determination of Alkali and Halide Mono-valent Ion Parameters for Use in Explicitly Solvated Biomolecular Simula-tions,” The Journal of Physical Chemistry B , 9020–9041 (2008). W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, “Development andTesting of the OLPS All-Atom Force Field on Conformational Energeticsand Properties of Organic Liquids,” J. Am. Chem. Soc. , 11225–11236(1996), arXiv:_barata Materials and Techniques of polychrome woodensculpture. H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term ineffective pair potentials,” The Journal of Physical Chemistry , 6269–6271(1987). M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, andS.-P. Daniel, “The SIESTA method for ab initio order-N materials,” J. Phys.Cond. Mat. , 2745–2779 (2002). F. Corsetti, M. Fernández-Serra, J. M. Soler, and E. Artacho, “Optimalfinite-range atomic basis sets for liquid water and ice,” Journal of Physics:Condensed Matter , 435504 (2013). J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approx-imation made simple,” Physical Review Letters , 3865–3868 (1996),arXiv:0927-0256(96)00008 [10.1016]. F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre, and M.-V. Fernández-Serra, “Room temperature compressibility and diffusivity of liquid waterfrom first principles,” The Journal of Chemical Physics , 194502 (2013),arXiv:1307.1645. J. Wang, G. Román-Pérez, J. M. Soler, E. Artacho, and M.-V. Fernández-Serra, “Density, structure, and dynamics of water: The effect of van derwaals interactions,” The Journal of chemical physics , 024516 (2011). J. Wang, G. Román-Pérez, J. M. Soler, E. Artacho, and M.-V. Fernández-Serra, “Density, structure, and dynamics of water: The effect of van derWaals interactions,” The Journal of Chemical Physics , 024516 (2011). A. A. Chialvo, P. T. Cummings, H. D. Cochran, J. M. Simonson, and R. E.Mesmer, “Na + Cl – ion pair association in supercritical water,” The Journalof Chemical Physics , 9379–9387 (1995). G. H. Zimmerman, H. Arcis, and P. R. Tremaine, “Limiting conductiv-ities and ion association constants of aqueous NaCl under hydrothermalconditions: Experimental data and correlations,” Journal of Chemical andEngineering Data , 2415–2429 (2012). D. Pan, L. Spanu, B. Harrison, D. A. Sverjensky, and G. Galli, “Dielectricproperties of water under extreme conditions and transport of carbonatesin the deep Earth,” Proceedings of the National Academy of Sciences ,6646–6650 (2013). D. Lu, F. Gygi, and G. Galli, “Dielectric Properties of Ice and Liquid Waterfrom First-Principles Calculations,” Physical Review Letters , 147601(2008). D. C. Elton and M. V. Fernández-Serra, “Polar nanoregions in water:A study of the dielectric properties of TIP4P/2005, TIP4P/2005f andTTM3F,” Journal of Chemical Physics (2014), 10.1063/1.4869110,arXiv:1401.5090. upplemental: Role of Water Model on Ion Dissociation at AmbientConditions
Alec Wills
1, 2, a) and Marivi Fernández-Serra
1, 2, b) Physics and Astronomy Department, Stony Brook University. Stony Brook, New York 11794-3800,United States Institute for Advanced Computational Science, Stony Brook, New York 11794-3800, United States (Dated: 29 January 2021)
I. COMPLEMENT TO THE DISCUSSION SECTION
In addition to information contained in the PMFs, the dif-ferent thermodynamic quantities tell us how the ions affect thebulk as a whole. As the different models have been parame-terized to yield different behaviors, pressure changes betweensystems at different densities and between pure solvent andsolution can tell us about the equilibrium states.Moreover, when examining the discrepancies betweenmodels, it is important to investigate differences on a smallerscale than bulk thermodynamic values. Thus, micropscopicsolvent characteristics can be investigated by examining theradial distribution functions between the ions and neighboroxygen atoms. Local and long-range effects are apparent inthe solvation shells that are defined by these distributions, andcan tell us how the ions affect the solvent water.In calculating system averages for the ion solutions, we usethe calculated PMF values to thermodynamically average themean system values across the ion separation range. That is,for a given simulation average h X ( r i ) i at ion constraint dis-tance r i , we find the average across simulations as h X i = ∑ Ni e − β U ( r i ) h X ( r i ) i ∑ Ni e − β U ( r i ) , (1)where U ( r i ) is the PMF value at constraint distance r i . Thisway, presented averages take into account the configurationalpreference towards the favored state. A. Pressure
The mean pressure of the simulated systems gives key in-formation as to the level of thermodynamic equilibrium we es-tablish with our chosen parameters. The thermodynamicallyaveraged pressures for the ion solution and pure water sys-tems, P ions and P pure , respectively, are given in Table II. Addi-tionally, the differences in pressure between the same systemswith and without the ion pair is included.We see that there is minimal, though nonzero, effect in thechange of TIP4P system pressure with the inclusion of theions, whereby their presence would induce an expansion ofthe box size through the small positive pressure. However, the a) Electronic mail: [email protected] b) Electronic mail: [email protected] change in simulation density affects the pressure as we wouldexpect. With an equilibrium mass density of 0 .
994 g cm − , decreasing the volume (and thus density) bring the system fur-ther from equilibrium, explaining the change in magnitude ofpressure.Analagous to the classical force field simulations, we seeslight differences in adding the ion pair to ab initio watersystems. However, the equilibrium densities of PBE andthe dispersion-inclusive liquid water systems are very differ-ent. The mass density of PBE water at mechanical equilib-rium (zero pressure) is around 0 .
87 g cm − , whereas function-als that include dispersion corrections typically equilibrate tohigher densities, on the order of 1 . − . In this man-ner, inclusion of the ions again would increase seek the sys-tem volume for the PBE solutions to varying degrees, whichour results confirm. Addition of the ions to both PBE systemsincreases the density of the box, bringing it further away fromequilibrium, resulting in higher pressures for the ionic solu-tions. The 0.52 M system, already close to equilibrium, doesnot see a large pressure difference when the ions are included.This directly contrasts with the 0.56 M solution, which is fur-ther away from equilibrium and results in a pressure changean order of magnitude greater than the 0.52 M difference. Thisbehavior is mirrored in the vdW-DF-cx liquid, where the 0.56M solution is the closer to equilibrium and so reacts less to theion addition/removal.
B. Voronoi Volume
Table I shows the thermodynamically averaged volumes forthe Voronoi cells of the ions and water molecules. The wateraverages confirm what is shown in Fig. 5 – namely, that themodels all predict the same occupied per H O molecule andits value is only dependent on the density of the system. Onthe other hand, the volumes of the Na and Cl ions are, withinerror, largely independent of density effects in both classicaland ab initio simulations, with the ab initio volumes havingslightly larger averages that agree between the different XCparameterization and molarity combinations. There is a smallincrease in the ab initio ion Voronoi volume with decreas-ing solution volume, but between XC functionals the averagesagree and are larger than the corresponding classical density. a r X i v : . [ phy s i c s . c h e m - ph ] J a n C. Voronoi Entropy
Tessellations like those generated in the Voronoi diagramsand their connectivity have a measure of information content.In particular, the partitions have an associated information en-tropy that characterize their disorder.We may use this fact to construct the information entropyof the Voronoi diagram at a given simulation step: S Vor = − ∑ i P i ln P i , (2)where P i is the fraction of polyhedra with i faces. In two di-mensions, the probability is chosen with respect to i edges ofthe Voronoi polygon, the only unambiguous connection be-tween any two polygons in the tessellation. Analogously, inthree dimensions, the unambiguous connectors between thepolyhedra are their faces. It is for this reason we construct theprobabilities in Eq. 2 with the face number of a polyhedron, asit is the direct measure of connectivity in the tessellation. Thesum over i is taken from i = i faces for a given tessellation. In a completely ordered system(all polyhedra the same), S Vor = , thus we expect increasingdisorder with increasing value of the information entropy. It isimportant to note that this entropy is not the thermodynamicentropy, however they have been related in some models. Fig. 7 shows the calculated average Voronoi entropies fromthe sampled 0.56 M configurations from each ab initio simu-lation, shifted by the average entropy from a trajectory of thecorresponding system without the ions. A similar plot includ-ing the 0.52 M trajectories is shown in Fig. 7. We see thatthe 0.52 M simulations tend to have, on average, more disor-der in the respective Voronoi diagrams than the correspondinghigh density solutions for the ab initio simulations. Moreover,for each density there tends to be marginally more disorder inthe PBE simulations than in the vdW-DF-cx. Inclusion of theions in the classical trajectories does not induce much vari-ability in the average Voronoi entropies. Furthermore, we findthe effect of increasing the box size has a noticeably greatereffect on the average Voronoi entropies of the PBE system,varying by up to ∼ .
07, while density changes have a lesspronounced effect on the vdW-DF-cx system. This corrobo-rates the behavior of the drastic behavior change in the PBEPMF.Table II contains the mean Voronoi entropies for the differ-ent pure water systems, S Vor , pure . The results indicate minimaldifference within the same XC for different densities. How-ever, we see that, despite the slightly more structured O − Oradial distribution functions for PBE shown in Fig. 8a whencompared to the vdW-DF-cx systems, the PBE Voronoi en-tropy is higher. This is surprising, given the ice-like nature ofPBE water. This indicates more variability in the polyhedragenerated by the Voronoi tessellation in PBE systems.
D. Radial Distribution Functions and Coordination Numbers
Insight into solute ion effects on the solvent can be gainedfrom examining the ion-oxygen and oxygen-oxygen radial distribution functions (RDFs). Initial neutron diffraction anal-ysis led to results indicating minimal affects beyond the firstion hydration shell, while later results showed that ion effectsextend beyond the first solvation shell. Even further resultssuggested that cosmo/chaotropic effects have minimal bear-ing on molecular interactions between the solute and solvent,despite the changes in solvent structure. X-ray absorptionstudies have indicated cationic destabilization of the hydro-gen bond network indicative of extended cationic effects onthe solvent, while anionic effects from Cl − indicated insignif-icant structural effects. However, GGA-level ab initio simulations have suggestedthat Na + affects the solvent only locally (not beyond thefirst solvation shell), while Cl − has longer-range effects, di-rectly contrasting with interpretations of experimental results.Specifically, it was shown that the Cl − ion weakens hydro-gen bonding up to the third solvation shell. We summarizeour simulation results in Table IV, where we list the coordi-nation numbers and distances of the first two radial minimain the ion-oxygen and oxygen-oxygen RDFs at different ionseparation distances and the thermodynamic averages.Notably, we see in Fig. 9 the effects the ions have onthe neighboring solution through the thermodynamically av-eraged radial distribution functions of the ion-oxygen pairs.The presence of a clear second shell in the Na − O distributions(Fig. 9b) is indicative of slightly extended perturbations on thesurrounding solvent molecules, what one might see referredto as a “cosmotropic” effect in that it introduces structure intothe surrounding solvent. An analogous effect is present in theCl − O distributions (Fig. 9c), although not as strongly as is thecase with the Na − O pair.Additionally, indications of local cosmotropic effects arepresent in Table IV, where the first hydration shell radius( r ) for the Na − O RDF is notably pulled in to lower dis-tances. Moreover, the PBE Na − O second coordination shellis markedly more contracted inwards than the vdW-DF-cxsimulations, for both molarities. This indicates that the PBEsolution is more prone to nonlocal cosmotropic perturbationsinduced by the Na + cation. Table IV further shows that theCl − O RDFs show extended perturbation in solvation struc-ture, in agreement with previous results. The solvation radiifor both the first and second solvation shells ( r and r , re-spectively) are notably shifted to higher distances, indicativeof a weakening of the hydrogen bond network in the vicinityof the Cl – ion – a “chaotropic” response.Fig. 9a also shows the O − O pair distributions of the ionicsolutions compared to experimental data. In all cases besidesthe TIP4P water model, the second shell is overstructuredcompared to experiment, further evidence for induced struc-tural effects in ab initio simulations that might bring aboutconflicting results in the PMF energetic ordering. Bulk struc-tural effects will influence the strength with which the solventinteracts with the solute ions, affecting the stability of eithergiven solvation state. C. Vega, C. McBride, E. Sanz, and J. L. F. Abascal, “Radial distributionfunctions and densities for the SPC/E, TIP4P and TIP5P models for liquidwater and ices Ih, Ic, II, III, IV, V, VI, VII, VIII, IX, XI and XII,” PhysicalChemistry Chemical Physics , 1450 (2005). M. J. Gillan, D. Alfè, and A. Michaelides, “Perspective: How good is DFTfor water?” Journal of Chemical Physics (2016), 10.1063/1.4944633. F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre, and M.-V. Fernández-Serra, “Room temperature compressibility and diffusivity of liquid waterfrom first principles,” The Journal of Chemical Physics , 194502 (2013),arXiv:1307.1645. M. Fritz, M. Fernández-Serra, and J. M. Soler, “Optimization of anexchange-correlation density functional for water,” Journal of ChemicalPhysics (2016), 10.1063/1.4953081, arXiv:1603.07302. E. Bormashenko, M. Frenkel, A. Vilk, I. Legchenkova, A. A. Fedorets,N. E. Aktaev, L. A. Dombrovsky, and M. Nosonovsky, “Characterizationof self-assembled 2D patterns with voronoi entropy,” Entropy (2018),10.3390/e20120956. V. Senthil Kumar and V. Kumaran, “Voronoi cell volume distribution andconfigurational entropy of hard-spheres,” The Journal of Chemical Physics , 114501 (2005). A. K. Soper and K. Weckström, “Ion solvation and water structure inpotassium halide aqueous solutions,” Biophysical Chemistry , 180–191(2006). R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, “Per- turbation of water structure due to monovalent ions in solution,” PhysicalChemistry Chemical Physics , 2959–2967 (2007). R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, “Hydra-tion of sodium, potassium, and chloride ions in solution and the conceptof structure maker/breaker,” Journal of Physical Chemistry B , 13570–13577 (2007). I. Waluyo, C. Huang, D. Nordlund, U. Bergmann, T. M. Weiss, L. G. M.Pettersson, and A. Nilsson, “The structure of water in the hydration shell ofcations from x-ray Raman and small angle x-ray scattering measurements,”The Journal of Chemical Physics , 064513 (2011). I. Waluyo, D. Nordlund, U. Bergmann, D. Schlesinger, L. G. M. Petters-son, and A. Nilsson, “A different view of structure-making and structure-breaking in alkali halide aqueous solutions through x-ray absorption spec-troscopy,” The Journal of Chemical Physics , 244506 (2014). A. P. Gaiduk and G. Galli, “Local and Global Effects of Dissolved SodiumChloride on the Structure of Water,” The Journal of Physical ChemistryLetters , 1496–1502 (2017). A. K. Soper, “The Radial Distribution Functions of Water as Derived fromRadiation Total Scattering Experiments: Is There Anything We Can Say forSure?” ISRN Physical Chemistry , 1–67 (2013). S UPPLEMENTARY T ABLE
I. Thermodynamically averaged Voronoi volumes for the atomic ions and water molecules.Model Molarity (M) h V Vor , Na i (Å ) h V Vor , Cl i (Å ) h V Vor , water i (Å )TIP4P 0.56 17 . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . UPPLEMENTARY T ABLE
II. Thermodynamic averages of the Voronoi entropies, pressures, and temperatures of different systems.Model Molarity (M) S Vor , ions S Vor , pure ∆ S Vor P ions (kBar) P pure (kBar) ∆ P (kBar) T (K)TIP4P 0.56 3 . ± .
02 3 . ± .
02 0 ± . − . ± . − . ± .
15 0 . ± . . ± . . ± .
02 3 . ± .
02 0 ± . − . ± . − . ± . − . ± . . ± . . ± .
03 2 . ± .
03 0 . ± .
04 2 . ± . . ± .
16 0 . ± . . ± . . ± .
03 2 . ± .
03 0 . ± .
03 0 . ± . . ± .
14 0 . ± . . ± . . ± .
03 2 . ± .
03 0 ± . − . ± . − . ± .
99 0 . ± . . ± . . ± .
03 2 . ± .
03 0 . ± . − . ± . − . ± . − . ± . . ± . Ion solution weighted averages are computed using Boltzmann factors weighted by the PMF values at a given distance. Individual trajectory averages forpressure and temperature values are calculated with time averaging 100 equally spaced subtrajectories to estimate error. The parameters for 0.56 M and 0.52 Msimulations when used in the pure water simulations (no ions) correspond to 0.97 g · cm − and 0.90 g · cm − , respectively. S UPPLEMENTARY T ABLE
III. Average temperatures for the pure water systems.Model Concentration ∗ (M) ρ (g · cm − ) T (K)TIP4P 0.56 0 .
97 297 . ± . .
90 300 . ± . .
97 399 . ± . .
90 423 . ± . .
97 379 . ± . .
90 382 . ± . ∗ Concentrations are given to clarify that the box parameters used for the pure water simulations are the same as those of the given molarity when the ions areincluded. S UPPLEMENTARY T ABLE
IV. Characteristics of the Cl − O, Na − O, O − O radial distribution functions at ion separations of 2.8, 4.8, and 6.0 Åand the thermodynamic averages of the values across all separations. N is the coordination number in the first solvation shell, while r and r are the minima (in Å) defining the first and second shell, respectively. Cl − O r r N r r N r r N r r N TIP4P 0.56 3.63 5.87 5.36 3.85 5.98 7.30 3.88 6.03 7.30 3.79 5.95 6.620.52 3.68 5.92 5.46 3.88 6.01 7.21 3.88 6.12 7.08 3.85 5.98 6.58PBE 0.56 3.21 5.50 2.21 4.07 6.15 7.74 3.99 6.03 7.04 3.15 5.92 2.020.52 3.21 5.36 1.89 3.77 5.89 5.14 3.57 5.87 4.43 3.79 5.89 5.62vdW-DF-cx 0.56 3.18 5.53 1.85 3.77 6.03 6.18 3.99 6.15 6.93 3.71 5.73 5.810.52 3.77 5.84 4.82 3.71 6.06 5.07 3.77 5.95 5.17 3.77 5.84 5.51 Na − O r r N r r N r r N r r N TIP4P 0.56 3.32 5.53 4.67 3.23 5.53 5.94 3.29 5.53 5.97 3.26 5.53 5.130.52 3.23 5.53 4.46 3.35 5.56 5.89 3.23 5.59 5.84 3.23 5.53 4.86PBE 0.56 3.01 5.19 3.59 3.23 5.25 4.94 3.18 5.22 5.04 3.21 5.31 4.640.52 3.40 5.64 3.75 3.15 5.56 4.12 3.26 5.31 4.80 3.21 5.56 4.44vdW-DF-cx 0.56 3.15 5.28 3.80 3.07 5.31 5.25 3.09 5.39 4.57 3.18 5.39 4.930.52 3.09 5.28 3.55 3.40 5.42 4.85 3.15 5.42 4.43 3.29 5.39 4.63 O − O r r N r r N r r N r r N TIP4P 0.56 3.40 5.61 4.91 3.35 5.59 4.68 3.35 5.61 4.68 3.40 5.61 4.910.52 3.40 5.70 4.60 3.40 5.67 4.60 3.40 5.67 4.59 3.40 5.70 4.60PBE 0.56 3.43 5.64 4.87 3.32 5.59 4.33 3.26 5.61 4.18 3.32 5.59 4.330.52 3.29 5.53 4.03 3.35 5.59 4.15 3.26 5.50 3.90 3.35 5.59 4.15vdW-DF-cx 0.56 3.29 5.47 4.27 3.35 5.59 4.45 3.32 5.56 4.27 3.32 5.56 4.290.52 3.26 5.53 3.87 3.32 5.47 3.95 3.29 5.56 3.91 3.29 5.53 3.96Experiment – – – – – – – – – – 3.39 5.58 4.724 Timestep redshifting. (cm ) I n t e n s i t y ( a r b . u . ) t = 0.5 fs t = 1.5 fs (a) t (fs) ( c m ) (b) t (fs) ( c m ) (c) S UPPLEMENTARY F IGURE (a) : The vibrational spectra of water,simulated with a 0.5 fs timestep (blue) and a 1.5 fs timestep (red). (b) : The resonance peak of the ν (HOH angle bending) mode astime a function of simulation timestep. (c) : The resonance peak ofthe ν + (stretching) modes as a function of the timestep. A largeredshift is evident as the timestep increases. r ( Å ) g OO ( r ) t=0.5 fst=1.5 fsSoper (2013) S UPPLEMENTARY F IGURE
2. The O − O radial distributionfunctions for pure water simulated with different timesteps at thesame mass density, plotted alongside experimental results. The firstshell of the 1.5 fs simulation is clearly destructured compared toexperiment, while the 0.5 fs shell accurately captures the first shellpeak. d Na Cl ( Å ) C oo r d . N u m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V ( m e V ) S UPPLEMENTARY F IGURE
3. The two-dimensional free energysurface resulting from the umbrella sampling as mentioned in themain text. d Na Cl ( Å ) V ( e V ) S UPPLEMENTARY F IGURE
4. PMFs generated with differentparameter choices in 0.56 M solution. There is no siginificantdifference in the predicted behavior between SPC/E and SPCsolvents with the same ion parameters. Note the enhanced stabilityof the SPC/E (blue) and SPC (red) SSIP states when compared tothose in main body TIP4P simulations. Moreover, changing fromthe OPLS-AA parameters in the SPC/E solvent to the JC parameters(green) further enhances the stability of the SSIP, and affects theCIP location as well. d Na Cl ( Å ) V V o r ( Å ) TIP4P
Na (0.56 M) Cl (0.56 M) (a) d Na Cl ( Å ) V V o r ( Å ) PBE
Na (0.56 M) Cl (0.56 M) (b) d Na Cl ( Å ) V V o r ( Å ) vdW-DF-cx Na (0.56 M) Cl (0.56 M) (c) S UPPLEMENTARY F IGURE (a-c) Voronoi volumes as a function of ionic distance for the 0.56 M TIP4P, PBE, vdW-DF-cx simulations.The red and blue lines correspond to Voronoi volume averages for Na and Cl, respectively. Average H O volumes are constant across thesimulations and are not plotted, but are presented in Table I. d Na Cl ( Å ) V V o r ( Å ) TIP4P
Na (0.52 M) Cl (0.52 M) (a) d Na Cl ( Å ) V V o r ( Å ) PBE
Na (0.52 M) Cl (0.52 M) (b) d Na Cl ( Å ) V V o r ( Å ) vdW-DF-cx Na (0.52 M) Cl (0.52 M) (c) S UPPLEMENTARY F IGURE (a-c) Voronoi volumes as a function of ionic distance for the 0.52 M TIP4P, PBE, vdW-DF-cx simulations.The red and blue lines correspond to Voronoi volume averages for Na and Cl, respectively. Average H O volumes are constant across thesimulations and are not plotted, but are presented in Table I. d Na Cl ( Å ) S V o r TIP4P (0.56 M)TIP4P (0.52 M) (a) d Na Cl ( Å ) S V o r PBE (0.56 M)PBE (0.52 M) (b) d Na Cl ( Å ) S V o r vdW-DF-cx (0.56 M)vdW-DF-cx (0.52 M) (c) S UPPLEMENTARY F IGURE (a-c) Voronoi entropies as a function of ionic distance for the 0.56 M TIP4P, PBE, vdW-DF-cx simulations.The solid and dashed lines correspond to Voronoi entropy averages for 0.56 M and 0.52 M solutions, respectively. d O O ( Å ) g OO ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M)Soper (2013) (a) Pure Water Solutions d O O ( Å ) g OO ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M)Soper (2013) (b) Ionic Solutions S UPPLEMENTARY F IGURE (a) The O − O radial distributionsfunctions of pure water using TIP4P, vdW-DF-cx, and PBE (green,red, and blue, respectively) with the box size that corresponds to the0.56 M solutions with the ions removed. (b)
Thetherymodynamically averaged O − O radial distributions functions ofthe ion solutions for the TIP4P, vdW-DF-cx, and PBE simulationswith concentrations of 0.56 M. Both figures show experimentalresults at ambient conditions for comparison. The inclusion of theions serves to increase the structure in the ab initio cases. d O O ( Å ) g OO ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M)Soper (2013) (a) O − O radial distribution functions. d Na O ( Å ) g N a O ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M) (b) Na − O radial distribution functions. d Cl O ( Å ) g C l O ( r ) vdW-DF-cx (0.56 M)PBE (0.56 M)TIP4P (0.56 M) (c) Cl − O radial distribution functions. S UPPLEMENTARY F IGURE
9. Thermodynamically averaged X − Oradial distribution functions for the 0.56 M solutions. In each, thevdW-DF-cx, PBE, and TIP4P solutions are red, blue, and green,respectively. (a)
Plotted with the simulated O − O radial distributionfunctions is experimental data. (b) The Na − O radial distributionfunctions. Notable is the clear presence of a second shell, indicativeof cosmotropic effects beyond the first solvation shell. (c)
TheCl − O radial distribution functions, indicating the similar structuraleffects as seen in the Na −−