Ion Pairing in Alkali Nitrate Electrolyte Solutions - The Journal of

Feb 22, 2016 - Ion hydration (free) energies and radii of ion first hydration shell are often used to parametrize force field of ions.(25, 28) We vali...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Ion Pairing in Alkali Nitrate Electrolyte Solutions Wen Jun Xie, Zhen Zhang, and Yi Qin Gao* Institute of Theoretical and Computational Chemistry, College of Chemistry and Molecular Engineering, Beijing National Laboratory of Molecular Sciences, and Biodynamic Optical Imaging Center, Peking University, Beijing 100871, China S Supporting Information *

ABSTRACT: In this study, we investigate the thermodynamics of alkali nitrate salt solutions, especially the formation of contact ion pairs between alkali cation and nitrate anion. The ion-pairing propensity shows an order of LiNO3 < NaNO3 < KNO3. Such results explain the salt activity coefficients and suggest that the empirical “law of matching water affinity” is followed by these alkali nitrate salt solutions. The spatial patterns of contact ion pairs are different in the three salt solutions studied here: Li+ forms the contact ion pair with only one oxygen of the nitrate while Na+ and K+ can also be shared by two oxygens of the nitrate. In reproducing the salt activity coefficient using Kirkwood−Buff theory, we find that it is essential to include electronic polarization for Li+ which has a high charge density. The electronic continuum correction for nonpolarizable force field significantly improves the agreement between the calculated activity coefficients and their experimental values. This approach also improves the performance of the force field on salt solubility. From these two aspects, this study suggests that electronic continuum correction can be a promising approach to force-field development for ions with high charge densities.

■

computationally,6 very few studies touch on ion pairs in alkali nitrate bulk solutions. One reason is that the force fields of halide ions are extensively and relatively well developed, while that for nitrate ion is not.7 Nevertheless, there is a trend to study ion pairing in complex ions including charged groups in biomolecules, which is of great importance to reveal the underlying molecule mechanism of the well-known “Hofmeister series”.4,8−10 In this study, we use MD simulation to investigate ion pairing in three alkali nitrate salt solutions at different moderate concentrations. The electronic continuum correction (ECC)11 that effectively takes electronic polarization effect into account in a mean-field way is used to reparameterize the nonpolarizable force field of both alkali cations and nitrate anion, the accuracy of which is assessed by comparing experimental activity coefficients with their calculated values from MD simulations using Kirkwood−Buff theory.12,13 The ECC approach has already been successfully used for different ions, especially those with high charge densities where electronic polarization is thought to be important.14−16 Even for salts of low charge density, this approach was shown to qualitatively reproduce the water diffusion trends which is misrepresented using nonpolarizable force field.17 It is also shown here that the nonpolarizable force field with ECC significantly improves its capability of reproducing experimental salt activity coefficients and solubilities. With such an accurate force field, the ion-pair-

INTRODUCTION Alkali nitrate salts including lithium nitrate (LiNO3), sodium nitrate (NaNO3), and potassium nitrate (KNO3) abound naturally on earth.1 These salts are involved in many life processes. Both alkali cations and nitrate anion can be assimilated by plants. In addition, dietary nitrate ion in the body can form nitrite ion and then be reduced to nitric oxide, which reduces hypertension.2 Since all of the above processes take place in an aqueous environment, studying the thermodynamic behavior of alkali nitrates in aqueous solution can be very important for further industrial and biological applications. One core concept in our understanding of salt solutions is ion pairing, which describes the formation of ion pairs between cation and anion in salt solutions.3,4 In liquids, thermal agitation tends to keep ions diffused away from each other. When the electrostatic attraction energy between the oppositely charged ions is larger than thermal fluctuation, they can form distinct entities such as ion pairs, which is thought to influence several properties of electrolyte solutions, including activity coefficients of some salts at moderate and high concentrations. 5 Quantifying the ability of ion-pair formation in different salt solutions is thus essential to understanding the microscopic structure and thermodynamics of electrolyte solution. Ion pairs can be detected by different experiments, such as conductometry and spectroscopic measurements.3 In addition, molecular dynamics (MD) simulations are extensively used in such studies, which are complementary to experiments and can provide a detailed molecular picture of the ion pairs. Contrary to alkali halide solutions which have been largely investigated © XXXX American Chemical Society

Received: November 3, 2015 Revised: February 9, 2016

A

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

using classical MD simulation. In the simulation system, a large water box with 39304 water molecules were first constructed. Different numbers of ion pairs (600, 900, and 1200) were then used to displace water molecules at random positions. The salt concentration is 0.87, 1.33, and 1.81 m (mol/kg), respectively. The solubility at ambient temperature (25 °C) for LiNO3, NaNO3, and KNO3 is 14.8, 10.7, and 3.8 m, respectively.20 All MD simulations were performed using the AMBER 12 suite of programs under periodic boundary conditions.21 After the initial equilibration, we carried out 60 ns NPT ensemble calculations for each system with a time step of 2 fs. The Berendsen weak-coupling barostat was used to control the pressure at 1 bar, and the pressure relaxation time was 2.0 ps.22 The temperature was controlled using Langevin dynamics at 300 K and a collision frequency of 5.0 ps−1. The neighbor list was updated whenever any atom had moved more than 1.0 Å since the last list update.21 The SHAKE algorithm was used to constrain all bonds involving hydrogen atoms.23 A cutoff of 1.0 nm was used for van der Waals interactions. A long-range dispersion correction based on an analytical integral assuming an isotropic, uniform bulk particle distribution beyond the cutoff was added to the van der Waals energy and pressure.21 The Particle Mesh Ewald was applied to treat long-range electrostatic interactions and the truncation between real space and reciprocal space was 1.0 nm. Lorentz/Berthelot mixing rules were used to evaluate Lennard-Jones parameters between different atom types. Water was described with the TIP3P model,24 which is extensively used, especially in the simulation of biomolecules where ion pairing can play important roles. Force-field parameters for alkali cation and nitrate anion were taken, respectively, from the work of Joung et al.25 and Baaden et al.,26 all of which are nonpolarizable. However, these force fields failed to reproduce the activity coefficient qualitatively (see more details later). Such a result suggests the weakness of the standard nonpolarizable force field in studying thermodynamics of these nitrate solutions. Recently, a rigorous and efficient mean-field way to include electronic polarization effects was proposed by Leontyev and Stuchebrukhov, which is called ECC.11 In this approach, the partial charges in nonpolarizable force field are scaled by a factor of 1/(ϵel)1/2, where ϵel is the relative permittivity of the solvent. The ϵel is due to the electronic (high-frequency) dielectric constant and for water it is 1.78.27 In doing so, the electronic polarization effect is effectively compensated. It should be mentioned that many water models like TIP3P used here have already taken the electronic polarization into account by fitting a variety of experimental results. We then performed MD simulations using the above nonpolarizable force fields with the ECC for ions, in which the partial charges of the ions are rescaled by a factor of 1/(1.78)1/2 (we call it ECC force field). The force field parameters used in this study are summarized in Tables S1 and S2. The ECC force field semiquantitatively or even quantitatively captures the characteristics of alkali nitrate activity coefficients as shown in the following. Ion hydration (free) energies and radii of ion first hydration shell are often used to parametrize force field of ions.25,28 We validated these properties for the parameters of nitrate anion in studying the hydration dynamics of nitrate anion.29 The results of alkali cations are presented in the SI and also coincide well with previous results. Based on these observations, we conclude that the ECC force field parameters of alkali cation and nitrate anion can be combined to study the thermodynamics of alkali

formation propensity and the molecule structure of ion pairs in different alkali nitrate solutions are investigated and discussed. Significant differences between halides and nitrates are noticed, which largely arise from the different molecular geometry of the anions.

■

EXPERIMENTAL ACTIVITY COEFFICIENTS OF ALKALI NITRATE SALT SOLUTIONS In electrolyte solution, chemical potentials of cations and anions cannot be measured separately due to the electroneutrality restriction. Thus, both ions are considered as one chemical species, a “salt molecule”.13,18 The chemical potential of salt molecule S is μ = μ * + RT ln aS = μ * + RT ln γ ρ (1) S

S

S

S S

where μS* is its chemical potential in the standard state, aS is the activity of the salt solute, ρS is the salt concentration, and γS is the activity coefficient. In infinite dilution, the salt obeys Henry’s law and is in the standard state, where the activity coefficient is defined as unity. With increasing salt concentrations, the electrolyte solution gradually deviates from ideality. The experimental activity coefficients of alkali nitrate solutions are shown in Figure 1. It is apparent that the

Figure 1. Activity coefficients of LiNO3, NaNO3, and KNO3 at different salt concentrations (data taken from ref 5).

behavior of LiNO3 in aqueous solution is quite different from the other two alkali nitrates. For LiNO3, the salt activity coefficient decreases first and then increases with the addition of salt. From the Debye−Hückel theory, the decrease of the activity coefficient is due to the long-range electrostatic forces. Its increase at even higher concentrations is due to the shortrange interactions between different entities in the solution.5 In contrast, the activity coefficients decrease at all salt concentrations for NaNO3 and KNO3. It is believed that ion pairs form in KNO3 salt solution based on the simple argument that ion pairing tends to lower the salt activity coefficient.19 However, there remain several questions. What is the molecular basis of such an argument? Are there also significant amounts of ion pairs in NaNO3 solution as in KNO3 solution? What does the ion pair look like in different alkali salt solutions? This work is aimed to address these questions.

■

METHODS Alkali nitrate salt solutions including LiNO3, NaNO3, and KNO3 with three different salt concentrations were investigated B

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B nitrate salt solutions although they are from two different sources. Kirkwood−Buff theory relates microscopic radial distribution function (RDF) between different components in the solution with macroscopic thermodynamic properties, e.g., chemical potential, partial molar volume, and compressibility.12,13 The computed RDFs in MD simulations can be used to yield the thermodynamic properties with Kirkwood−Buff theory. We use this strategy to quantify the force field and analyze thermodynamics information on nitrate solutions. The alkali cation and nitrate anion are together treated as the “salt molecule”. The derivative of the activity coefficient of salts with respect to the electrolyte molar concentration in this two component electrolyte solution systems (water and salt) can be expressed as ⎛ ∂ ln a ⎞ ⎛ ∂ ln γ ⎞ s s ⎟⎟ = 1 + ⎜⎜ ⎟⎟ as′ = ⎜⎜ ρ ρ ln ln ∂ ∂ ⎝ ⎝ s ⎠ p,T s ⎠ p,T =

1 1 + ΔNss − ΔNws

ΔNjs = ρs Gjs

Gjs = 4π

∫0

∞

(j = s or w)

(g js(r ) − 1)r 2dr

(2) (3)

(4)

where a′s is the derivative of the salt activity, ΔNjs is the excess coordination number of component s around component j, and Gjs is the Kirkwood−Buff integral (KBI) between j and s, which is defined as the spatial integral of the RDF gjs. A positive value of ΔNjs corresponds to an excess of component j around s, while a value of ΔNjs less than zero means a depletion of component j around s.

■

RESULTS AND DISCUSSION 1. Salt Activity Coefficients Using ECC Force Field. We first calculated the excess coordination number of salt and water around salt molecules in the simulations using the ECC force field, denoted as ΔNss and ΔNws, respectively. As illustrated in Figure 2, these two numbers for the three different alkali nitrate solutions at 1.33 m differ a lot even in terms of their relative order. For LiNO3 solution, the value of ΔNss is smaller than ΔNws. From eq 2, it is predicted that for this solution the salt activity derivative is larger than unity. In contrast, the salt activity derivative would then be less than unity in NaNO3 and KNO3 solutions. These numbers at 0.87 and 1.81 m salt concentration are presented in Table S4. The sign and the relative order of ΔNss and ΔNws in different salt concentrations are the same, and only the values differ. If the force fields are accurate enough, there should not be any precipitation in all the above simulations. However, precipitation occurs in the KNO3 solution at 1.81 m which is about half of the experimental salt solubility. The Kirkwood−Buff theory is not expected to work for such systems and the associated activity coefficients are not calculated. The discussions about salt solubility were presented in the last part of this section. One can certainly integrate over the activity derivative to obtain the salt activity. However, the activity derivative which is related to excess numbers of different components around salt has a more straightforward physical meaning. Thus, the activity derivative instead of the activity coefficient itself is used here to interpret the data. The predicted activity derivatives using

Figure 2. Excess coordination number of salt around salt (ΔNss) and water around salt (ΔNws) in the 1.33 m salt solutions using ECC force field: (A) LiNO3, (B )NaNO3, (C) KNO3.

Kirkwood−Buff theory are compared with experimental values in Figure 3. The predicted values coincide well with the experiments in the following aspects: the relative order of the three alkali salt solutions follows LiNO3 > NaNO3 > KNO3; the value in LiNO3 solution is larger than unity in contrast with the other two; with the increase of salt concentration, the activity derivative increases for LiNO3, while it decreases for NaNO3 and KNO3, especially for KNO3. Kirkwood−Buff theory is an exact theory derived from statistical mechanics in the μVT ensemble without any other assumptions. In principle, KBIs should be calculated in open systems, but in practice, they are usually obtained from closed boundary NPT or NVT simulations. The fluctuations at a length of several molecular diameters are important in the convergence of KBIs since they are spatial integrals of RDFs over infinite distances. The simulation system must be large enough in order to allow the fluctuations which can be unaffordable for extensive calculations. Usually, empirical finite C

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

ion-pairing propensity, is much higher in KNO3 solution than that in LiNO3 or NaNO3 solution. Meanwhile, Na+ is more prone to forming ion pairs than Li+. This order of formation of contact ion pair which decreases with the increasing cation charge density is in accordance with that between divalent metals and nitrate anion.31 In elucidating salt solution structures, the diffraction method has limitations since several interatomic correlations overlap in the range of the radial distance.32,33 In our case, the overlap of the distance correlations between water oxygen-nitrate nitrogen pair and cation-nitrate nitrogen pair causes the difficulty in quantifying contact ion pair properties. Compared with an earlier simulation study of NaNO3 salt solution,33 we obtained similar RDF peak positions of contact ion pairs. In our simulation, the peak of the RDF between Na+ and nitrate nitrogen (oxygen) positions at 3.45 (2.45) Å. The corresponding values in ref 33 are 3.48 and 2.43 Å, respectively. Nitrate is a complex anion with one central nitrogen atom surrounded by three identically bonded oxygen atoms in a trigonal planar arrangement. All three oxygen atoms can act as contact ion pairing site. In our calculation, a cation is considered to form contact ion pairs with nitrate oxygen provided it lies within the distance of the first minimum of the RDF between the cation and nitrate oxygen. We quantified the spatial distribution of cations around nitrate anion and show the results in Figure 5. The distribution of cations in the first

Figure 3. Values of activity derivatives from ECC force field simulation. Experimental values (numerical fitting to experiments) are plotted in lines. Predicted values from MD simulation data using Kirkwood−Buff theory are plotted: LiNO3 (square), NaNO3 (circle), KNO3 (triangle). The associated error bars obtained from separating the 60 ns simulation data of each system into three equal blocks are all within 0.02, which is very small and not shown in this figure.

size corrections are necessary when the simulation box is small.30 In this study, we chose to use simulation systems that are large enough and also a long enough simulation time to get converged excess numbers (also KBIs) without any corrections. The convergence of the calculated excess numbers can be seen in Figure 2. As seen from Figure 2, the values of ΔNws differ slightly in the three different electrolyte solutions, whereas ΔNss differs significantly. This comparison suggests that the interaction between ions, instead of the hydration of ions, determines the salt activity, although the hydration abilities of different alkali cations are quite different and lithium has the strongest hydration. From the electroneutrality condition for a 1:1 salt one can show that ρSG± = ρSGss + 1 = ΔNss + 1, where G± is the KBI between cation and anion. Thus, we can quantify the interaction between cation and nitrate anion directly through the value of ΔNss. ΔNss increases in the order of LiNO3, NaNO3 and KNO3, which means that the propensity of cation locating in the vicinity of nitrate anion increases in this order. 2. Thermodynamic Properties of Alkali Nitrate Salt Solutions in ECC Force Field Simulation. From the RDF between different alkali cations and nitrate anion (Figure 4), it can also be easily seen that K+ tends to form contact ion pairs with nitrate. The first peak in the RDF, which represents the

Figure 5. Illustration of alkali cations (represented by the blue isosurface) sitting around the nitrate anion in the corresponding 1.33 m alkali nitrate solutions: (A) Li+; (B) Na+; (C) K+. The isosurface of Li+ above illustrates the density of Li+ in this region is 3.7 times the bulk salt density, while this multiple for Na+ and K+ is 6.4.

hydration shell of nitrate can be decomposed into two classes: in the first one the cations are in direct contact with only one oxygen atom of the nitrate (monovalent coordination type), and in the second cations are shared by two oxygen atoms (bivalent coordination type). It is clearly that the spatial distribution of alkalis is quite different in the three nitrate solutions. In LiNO3 aqueous solution, lithium ion only forms contact ion pairs with one oxygen atom of nitrate anion. In contrast, both sodium and potassium also form contact ion pairs simultaneously with two oxygen atoms of nitrate anion. Cations that are in contact with all the three oxygen atoms of nitrate locate well above the molecular plane with a population of nearly zero and can be neglected in all three alkali nitrate solutions. The numbers of alkali cations residing around one nitrate anion in the different systems are summarized in Table 1. It can be seen that in all systems studied here there are fewer shared contact pairing cations than those in contact with only one oxygen. In the LiNO3 solution, the numbers of the two types of nitrate bound cations are both smaller than the corresponding values in the other two systems, but its

Figure 4. RDF between alkali cations and nitrate anion in the 1.33 m salt solutions. D

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 1. Numbers of Alkali Cations around One Nitrate in Different Salt Solutions concn (M)

coord type

Li+

Na+

K+

0.87

monovalent bivalent monovalent bivalent monovalent bivalent

0.06 0.01 0.09 0.01 0.13 0.01

0.24 0.11 0.36 0.17 0.49 0.24

0.31 0.26 0.46 0.40

1.33 1.81

monovalent/bivalent ratio is the largest. Also, the difference between NaNO3 and KNO3 systems mainly manifests in the number of shared contact pairing cations: the difference for monovalent cations is 0.07−0.10 while this number is 0.15− 0.23 for bivalent ones, with the dominance of the latter. When a cation enters into the two oxygen region, it loses more hydration water. The numbers of hydration water of cation, when it is monovalently or bivalently coordinated with one nitrate anion, are summarized in Table 2. Although different Table 2. Numbers of Hydration Water around Cation under Different Coordinating Conditions coord type

Li+

Na+

K+

fully hydrated monovalent bivalent

4.0 3.1 2.8

5.5 4.7 4.1

7.3 6.4 5.8

cations have different numbers of hydration water when they are fully hydrated, on average they all lose 0.8−0.9 hydration water molecules when they are monovalently coordinated to one nitrate anion. This number is 1.2−1.5 in the bivalent case. Since K+ has the lowest charge density among the three cations studied here, it is more easily to lose hydration water. The data above also suggests that the dehydration ability, other than the electrostatic interactions between cations and nitrate anion, dominates ion pairing probability. We also note here that all ion pairing values in 1.33 m are nearly the corresponding values averaged for 0.87 and 1.81 m. Such a relation also reflects that the ion pairing depends almost linearly on the salt concentration at least for the systems and the concentration range studied here. Until now, we focused on the analyses of single ion pairs. When different ion pairs share common ions, they are considered to belong to the same ion cluster. Increasing evidence shows that ion clusters play an important role in crystal nucleation.34 We then calculated the ion cluster size distribution for different alkali nitrate solutions (Figure 6). The cluster size is denoted by the sum of cation and anion numbers. In all cases, most of the ions are fully hydrated and do not involve in any contact ion pairs statistically. Although the ion clusters form and dissolve constantly, the ratio of fully hydrated ions also shows the contact ion pairing propensity. Consistent with the above analysis, for the solutions of the same concentration, the fraction of fully hydrated ions decreases in the order of LiNO3, NaNO3, and KNO3. In 0.87 m salt solution, the formation probability of the single ion pair which consists one cation and one anion (cluster size = 2) in NaNO3 is less than in KNO3. For a salt cluster of size 3, it can contain one cation and two anions or two cations and one anion. It is found that in all the systems the former is more likely to be

Figure 6. Cluster size distribution at three different salt concentrations: (A) 0.87 m; (B) 1.33 m; (C) 1.81 m.

observed (Table S7). In addition, about 80% salt clusters of size 4 contain two cations and two anions (Table S7). The ion pair separated by the thickness of only one water molecule is defined as solvent separated ion pair (SSIP). The fractions of SSIPs are summarized in Table 3. The majority of the fully hydrated ions mentioned above form SSIPs in all systems. The fractions of SSIPs also follow the order LiNO3 > NaNO3 > KNO3. The formation of SSIPs in LiNO3 increases along with the concentration from 61% to 77%. In contrast, the Table 3. Fractions of SSIPs in Different Salt Solutions

E

concn (M)

Li+

Na+

K+

0.87 1.33 1.81

0.61 0.72 0.77

0.49 0.47 0.41

0.38 0.32

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

In the simulation of LiNO3 solution using a nonpolarizable force field, ΔNss is shown to be larger than ΔNws. The predicted activity derivative is less than unity, which is in conflict with experimental value. The reason can be easily seen from the RDF between Li+ and nitrate anion (Figure 8). When

SSIPs in NaNO3 and KNO3 slightly decrease with the increasing of concentrations. Raman spectroscopy has been used to probe the solution structure of nitrate salt solutions.35−38 The in-phase symmetric stretching band v1 at ∼1049 cm−1 is sensitive to the cation type and salt concentration. Although the vibration frequencies of nitrate in different surrounding environments differ slightly, a component analysis of this band in the range of 1020−1080 cm−1 could yield the relative amount of different ionic species in NaNO3.37 The fraction of SSIPs from this spectroscopic component analysis in the concentration range studied in our simulation is around 50%. This coincides with our calculated results, which are 41%−49% (Table 3). In large ion clusters, nitrate anions in its center will vibrate with a frequency slighted different from that of free hydration nitrate, CIP or SSIP. The fraction of large ion clusters in NaNO3 below 1.2 m is negligible from Raman spectroscopy.37 In our simulation, a similar result was obtained. Ion clusters with a size above 6 are ∼1% in NaNO3 at 0.87 m, while this fraction is ∼5% at 1.33 m and ∼9% at 1.81 m. 3. Comparison between the Performance of ECC Force Field and Nonpolarizable Force Field. It was shown above that the ECC force field performs well in the calculations of thermodynamic properties of salt solutions. The results for nonpolarizable force field at different concentrations are presented here for comparison (Figure 7).As for the ECC

Figure 8. RDF between Li+ and nitrate anion in the 1.33 m salt solutions using ECC force field and nonpolarizable force field.

the nonpolarizable force field is used, the first peak of the RDF is too high which means a large population of contact ion pairs. This is also consistent with the above comparison of solution structures. From the calculated activity coefficient, we can conclude that the nonpolarizable force field for systems containing Li+ and nitrate anion is not accurate. For NaNO3 and KNO3 salt solutions, the KBIs calculated using ECC force field are similar to those calculated using the nonpolarizable force field. The RDF between Na+ and nitrate anion is shown in Figure 9 for the nonpolarizable force field

Figure 7. Values of activity derivatives from nonpolarizable force field simulation. All symbols are identical to those Figure 3. In three simulation systems, the salts precipitate, and thus, the salt activities are not calculated.

force field, KNO3 salt solution at 1.81 m precipitates when the nonpolarizable force field is used. In addition, precipitation also occurs in the simulation at KNO3 at 1.33 m and NaNO3 at 1.81 m in nonpolarizable force field simulations, although all the salt concentrations in our simulations are below the experimental solubility limit. In this sense, the ECC force field improves to a limited extent the performance of NaNO3 and KNO3 ion force fields. The properties of ion pairing and cluster formation using nonpolarizable force field are presented in Table S6 and Figure S1. The main difference is that ion pairing propensity of LiNO3 salt solution using nonpolarizable force field is two times that using ECC force field. The nonpolarizable force field used here fails to qualitatively reproduce the activity coefficients of LiNO3 salt solutions which will be elaborated in the following. The comparison between nonpolarizable and ECC force field could give clues in improving the accuracy of force field. We compared these two force fields at 1.33 m and more details are given below.

Figure 9. RDF between Na+ and nitrate anion in the 1.33 m salt solutions calculated using ECC force field and nonpolarizable force field.

simulation. In this simulation, the nitrate anion is surrounded by a stronger and more ordered cation shell than that seen from the simulation using ECC force field. Since the excess numbers mentioned above is the spatial integral of RDF, the error cancelation over the cation−anion distance makes the activity derivative nearly the same when ECC and nonpolarizable force fields are used for NaNO3 salt solution. Representative snapshots of ion distribution in KNO3 solution obtained using nonpolarizable and ECC force fields are presented in Figure 10. It can be clearly seen that ions described by the nonpolarizable force field precipitate and form F

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

simulation details are the same as those in the above simulations of salt solutions. In the ECC force field simulation, the salt solution contains 1743 water molecules and 941 ion pairs, which results in a concentration of about 30 m. In the nonpolarizable force field simulation, 3411 water molecules and 737 ion pairs were used, and the initial concentration is about 12 m. Both systems were shown to be supersaturated, and the ions gradually attach onto the surface of the crystal nuclei. The concentrations of the remaining salt solutions are presented in Figure 11. The [001] growth direction of the

Figure 10. Snapshot of ions in the 1.33 m KNO3 salt solution using different force fields: (A) salt precipitates in the simulation using nonpolarizable force field; (B) distribution of ions is homogeneous when ECC force field is used.

several large ion clusters (Figure 10A). Using the ECC force field, the KNO3 solution is homogeneous at 1.33 m, which is significantly below the experimental solubility limit (Figure 10B). The above comparison between the nonpolarizable and ECC force fields clearly shows that including the electronic polarization is important in reproducing the thermodynamics of salt solutions. Nonpolarizable force field results in a too high population of contact ion pairs for LiNO3 salt solution. Polarizable force field including the ECC force field is thought to stabilize the first hydration shell of ions and further inhibit ion clustering, especially for ions with high charge density.15 4. Salt Solubilities Using ECC Force Field and Nonpolarizable Force Field. In the development of most ion force fields, salt solubility is not set as a target property. A good performance of ion force field on salt solubility is essential for studies of crystal nucleation. The correctness of salt solubility requires an accurate description of the thermodynamics of both its liquid solution phase and solid crystal phase. The solubility of NaCl was studied recently, which was found to depend quite sensitively on the force field.39,40 Using some force fields of NaCl in water, the solubility is underestimated by a factor of 4.40 As shown above, the ECC force field of LiNO3 performs well in the thermodynamics of salt solution at all three concentrations studied here. One might expect this force field to reproduce the salt solubility. Next, we compared the salt solubility of LiNO3 using the ECC force field and the nonpolarizable force field. There are two commonly used methods of calculating salt solubilities in computer simulations. One method relies on the estimation of chemical potential of the simulated system in both solid crystal and salt solution.39,40 At the solubility limit, the chemical potential at two different phases should be equal. The other method (direct coexistence method) straightforwardly simulates the crystal growth in the presence of a large crystal nuclei.40 At supersaturation, salt ions tend to continually change into the solid phase around the crystal nuclei until the solubility limit. These two methods have been shown to be in reasonable agreement with each other for NaCl.40 In this study, we used the direct coexistence method to estimate the solubility of LiNO3. The solid phase which has a calcite rhombohedral structure41 and contains 600 ion pairs was immersed in a supersaturated salt solution. The solid phase could grow along its [001] direction. Semianisotropic pressure scaling was applied to control the pressure.21 All of the other

Figure 11. Molal concentration of LiNO3 versus simulation time for ECC force field and nonpolarizable force field. Dashed lines indicate the value at which the molality reaches a plateau. The last 100 ns simulation data were used to estimate the salt solubility.

LiNO3 crystal we chose has a net dipole which enables large initial growth rates. During equilibration, some ions have already grown into the solid crystal, which makes the concentrations at 0 ns in Figure 11 less than their initial values. After the initial decrease, the salt concentrations fluctuate, which indicates the salt solutions and crystal solids are in dynamic equilibrium. Simulation data (800 ns) were obtained for each system. The calculated solubility using nonpolarizable force field is 7.3 ± 0.2 m. In contrast, this value using ECC force field is 11.5 ± 0.4 m, within 25% of the experimental value. It can be concluded that the ECC force field for LiNO3 also significantly improves the description of its solubility. Considering the activity coefficients for NaNO3 and KNO3 are not exact although their qualitative behavior are right, we did not calculate the salt solubilities of NaNO3 and KNO3. Nevertheless, from the emergence of precipitation one can still conclude that the ECC improves the performance of force fields on salt solubilities. For NaNO3 and KNO3, the improvement of salt activity coefficients upon the usage of ECC force field is less obvious than the case of LiNO3. This is not surprising since the ECC was invented to treat the electronic polarization effect which is less severe for Na+ and K+. Such a complexity also reveals the role of van der Waals parameters which were not investigated in this work.42 It is likely that further parametrization can improve the force field’s performance on the activity coefficients and/or salt solubilities.

■

CONCLUSIONS In this study, we investigated carefully ion-pair formation in three alkali nitrate salt solutions at different salt concentrations using MD simulations and Kirkwood−Buff theory. It was found that the ion-pairing propensity decreases in the order of KNO3, G

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

■

NaNO3, and LiNO3, which provides an explanation of their differences in salt activity coefficients. The formation of ion pairs is the result of a competition between the electrostatic attraction of oppositely charged ions and the hydration of these ions, which can be roughly rationalized by the “Law of Matching Water Affinities”: cations and anions that have similar absolute hydration free energies tend to form contact ion pairs.43−45 Generally speaking, large and small anions tend to form ion pairs with large and small cations, respectively. This law does provide qualitative explanations to many observations, including the solubility of different salts. However, it was shown that this simple law fails to describe the interaction of alkali cations with acetate anion,8 which illustrates the complexity of ion interactions in salt solutions, especially complex ions like oxyanions. For nitrate anion, we found this law to work well: K+, which has the smallest charge density, tends to form contact ion pairs with the low charge density nitrate anion, while the ion-pair-formation tendency between Li+ and nitrate anion is much weaker. We also show that the polarization effect is essential for understanding the ion-pairing effect. ECC used here was used by Jungwirth et al. to study ion pairing in LiCl15 and CaCl246 salt solutions, where the charge density of cations are very high and can lead to severe ion clustering. In this study, we show that the ECC significantly improves the performance of force field in calculating LiNO3 salt activity coefficients and salt solubility. For NaNO3 and KNO3, the ECC has little effect on salt activity coefficients but improves the description of salt solubility although the van der Waals parameters were not optimized. Although the investigation of more ions with ECC is necessary for the general validation of this approach, such topics should be covered by a larger scale studies, which are beyond the scope of this study. The simple charge rescaling strategy without further adjustment of other force field parameters, including the van der Waals interaction, improves the result dramatically even for ion distribution between oil/ water interface.47 These results together also demonstrate that the correctness of partial charges is essential in ion force field parametrization. We also noticed that in several studies Kirkwood−Buff theory had been applied to parametrize the ion force field which in turn can be used to investigate the electrolyte solution.8,18 This approach takes the partial charges of ions as their valence, e.g., +1 for Na+ and −1 for Cl−. Such a treatment does not explicitly include the electronic screen of the electrostatic interactions between ions. This negligence of electronic polarizability is not easy to be compensated by adjusting only the van der Waals parameters, especially for salt solutions. The full charge ion force fields combined with commonly used water models which has already taken the electronic polarizability into account are, in principle, selfinconsistent, and it would be unsurprising to encounter conceptual difficulties in reproducing the experimental salt activity. The salt activity is a complex function of all the interaction parameters used in simulations. Both partial charges and other force field parameters should be parametrized carefully to mimic the experiments and reproduce the ion interactions at the microscopic level. This study shows that correcting the partial charges for ions can dramatically improve the capability of the simulations in reproducing salt activity coefficients.

Article

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b10755. More details about the force field as well as the results using different kinds of force fields (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Key Basic Research Foundation of China (2012CB917304) and NSFC (91027044, 21233002, and 21125311). We thank Alexei Stuchebrukhov for helpful discussions on the force field parameters.

■

REFERENCES

(1) Laue, W.; Thiemann, M.; Scheibler, E.; Wiegand, K. W. Nitrates and Nitrites. In Ullmann’s Encyclopedia of Industrial Chemistry; WileyVCH: Weinheim, 2000. (2) Nitrite and Nitrate in Human Health and Disease; Bryan, N. S., Loscalzo, J., Eds.; Humana Press: New York, 2011. (3) Marcus, Y.; Hefter, G. Ion Pairing. Chem. Rev. 2006, 106, 4585− 4621. (4) Jungwirth, P. Ion Pairing: From Water Clusters to the Aqueous Bulk. J. Phys. Chem. B 2014, 118, 10333−10334. (5) Robinson, R. A.; Stokes, R. H. Electrolyte Solutions; Dover Publications: New York, 2002. (6) Fennell, C. J.; Bizjak, A.; Vlachy, V.; Dill, K. A. Ion Pairing in Molecular Simulations of Aqueous Alkali Halide Solutions. J. Phys. Chem. B 2009, 113, 6782−6791. (7) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. Rational Design of Ion Force Fields Based on Thermodynamic Solvation Properties. J. Chem. Phys. 2009, 130, 124507. (8) Hess, B.; van der Vegt, N. F. A. Cation Specific Binding with Protein Surface Charges. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 13296−13300. (9) Xie, W. J.; Gao, Y. Q. A Simple Theory for the Hofmeister Series. J. Phys. Chem. Lett. 2013, 4, 4247−4252. (10) Xie, W. J.; Liu, C. W.; Yang, L. J.; Gao, Y. Q. On the Molecular Mechanism of Ion Specific Hofmeister Series. Sci. China: Chem. 2014, 57, 36−47. (11) Leontyev, I.; Stuchebrukhov, A. Accounting for Electronic Polarization in Non-Polarizable Force Fields. Phys. Chem. Chem. Phys. 2011, 13, 2613−2626. (12) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Solutions I. J. Chem. Phys. 1951, 19, 774−777. (13) Ben-Naim, A. Molecular Theory of Solutions; Oxford University Press: Oxford, 2006. (14) Pegado, L.; Marsalek, O.; Jungwirth, P.; Wernersson, E. Solvation and Ion-Pairing Properties of the Aqueous Sulfate Anion: Explicit Versus Effective Electronic Polarization. Phys. Chem. Chem. Phys. 2012, 14, 10248−10257. (15) Pluharova, E.; Mason, P. E.; Jungwirth, P. Ion Pairing in Aqueous Lithium Salt Solutions with Monovalent and Divalent Counter-Anions. J. Phys. Chem. A 2013, 117, 11766−11773. (16) Kohagen, M.; Pluhařová, E.; Mason, P. E.; Jungwirth, P. Exploring Ion−Ion Interactions in Aqueous Solutions by a Combination of Molecular Dynamics and Neutron Scattering. J. Phys. Chem. Lett. 2015, 6, 1563−1567. (17) Kann, Z. R.; Skinner, J. L. A Scaled-Ionic-Charge Simulation Model That Reproduces Enhanced and Suppressed Water Diffusion in Aqueous Salt Solutions. J. Chem. Phys. 2014, 141, 104507. H

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (18) Weerasinghe, S.; Smith, P. E. A Kirkwood-Buff Derived Force Field for Sodium Chloride in Water. J. Chem. Phys. 2003, 119, 11342− 11349. (19) Ninham, B. W.; Nostro, P. L. Molecular Forces and Self Assembly: in Colloid, Nano Sciences and Biology; Cambridge University Press: Cambridge, 2010. (20) CRC Handbook of Chemistry and Physics, 95th ed.; Haynes, W. M., Ed.; CRC Press: Boca Raton, FL, 2014. (21) AMBER 12; University of California, San Francisco, 2012. (22) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (23) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (24) Price, D. J.; Brooks, C. L. A Modified TIP3P Water Potential for Simulation with Ewald Summation. J. Chem. Phys. 2004, 121, 10096− 10103. (25) 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. (26) Baaden, M.; Burgard, M.; Wipff, G. TBP at the Water-Oil Interface: The Effect of TBP Concentration and Water Acidity Investigated by Molecular Dynamics Simulations. J. Phys. Chem. B 2001, 105, 11131−11141. (27) Leontyev, I. V.; Stuchebrukhov, A. A. Electronic Continuum Model for Molecular Dynamics Simulations of Biological Molecules. J. Chem. Theory Comput. 2010, 6, 1498−1508. (28) Dang, L. X.; Chang, T. M.; Roeselova, M.; Garrett, B. C.; Tobias, D. J. On NO3−-H2O Interactions in Aqueous Solutions and at Interfaces. J. Chem. Phys. 2006, 124, 066101. (29) Xie, W. J.; Yang, Y. I.; Gao, Y. Q. Dual Reorientation Relaxation Routes of Water Molecules in Oxyanion’s Hydration Shell: A Molecular Geometry Perspective. J. Chem. Phys. 2015, 143, 224504. (30) Ganguly, P.; van der Vegt, N. F. A. Convergence of Sampling Kirkwood-Buff Integrals of Aqueous Solutions with Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9, 1347−1355. (31) Xu, M.; Larentzos, J. P.; Roshdy, M.; Criscenti, L. J.; Allen, H. C. Aqueous Divalent Metal-Nitrate Interactions: Hydration Versus Ion Pairing. Phys. Chem. Chem. Phys. 2008, 10, 4793−4801. (32) Chialvo, A. A.; Vlcek, L. NO3− Coordination in Aqueous Solutions by 15N/14N and 18O/natO Isotopic Substitution: What Can We Learn from Molecular Simulation? J. Phys. Chem. B 2015, 119, 519−531. (33) Megyes, T.; Balint, S.; Peter, E.; Grosz, T.; Bako, I.; Krienke, H.; Bellissent-Funel, M. C. Solution Structure of NaNO3 in Water: Diffraction and Molecular Dynamics Simulation Study. J. Phys. Chem. B 2009, 113, 4054−4064. (34) Gebauer, D.; Kellermeier, M.; Gale, J. D.; Bergstrom, L.; Colfen, H. Pre-Nucleation Clusters as Solute Precursors in Crystallisation. Chem. Soc. Rev. 2014, 43, 2348−2371. (35) Mathieu, J. P.; Lounsbury, M. The Raman Spectra of Metallic Nitrates and the Structure of Concentrated Solutions of Electrolytes. Discuss. Faraday Soc. 1950, 9, 196−207. (36) Vollmar, P. M. Ionic Interactions in Aqueous Solution - A Raman Spectral Study. J. Chem. Phys. 1963, 39, 2236−2248. (37) Yu, J. Y.; Zhang, Y.; Tan, S. H.; Liu, Y.; Zhang, Y. H. Observation on the Ion Association Equilibria in NaNO3 Droplets Using Micro-Raman Spectroscopy. J. Phys. Chem. B 2012, 116, 12581−12589. (38) Irish, D. E.; Nelson, D. L.; Brooker, M. H. Quasilattice Features of Concentrated Aqueous LiNO3 Solutions. J. Chem. Phys. 1971, 54, 654−657. (39) Sanz, E.; Vega, C. Solubility of KF and NaCl in Water by Molecular Simulation. J. Chem. Phys. 2007, 126, 014507. (40) Aragones, J. L.; Sanz, E.; Vega, C. Solubility of NaCl in Water by Molecular Simulation Revisited. J. Chem. Phys. 2012, 136, 244508.

(41) Wu, X.; Fronczek, F. R.; Butler, L. G. Structure of LiNO3 PointCharge Model and Sign of the 7Li Quadrupole Coupling-Constant. Inorg. Chem. 1994, 33, 1363−1365. (42) Aragones, J. L.; Rovere, M.; Vega, C.; Gallo, P. Computer Simulation Study of the Structure of Licl Aqueous Solutions: Test of Non-Standard Mixing Rules in the Ion Interaction. J. Phys. Chem. B 2014, 118, 7680−7691. (43) Collins, K. D. Ions from the Hofmeister Series and Osmolytes: Effects on Proteins in Solution and in the Crystallization Process. Methods 2004, 34, 300−311. (44) Collins, K. D. Charge Density-Dependent Strength of Hydration and Biological Structure. Biophys. J. 1997, 72, 65−76. (45) Collins, K. D. Ion Hydration: Implications for Cellular Function, Polyelectrolytes, and Protein Crystallization. Biophys. Chem. 2006, 119, 271−281. (46) Kohagen, M.; Mason, P. E.; Jungwirth, P. Accurate Description of Calcium Solvation in Concentrated Aqueous Solutions. J. Phys. Chem. B 2014, 118, 7902−7909. (47) Vazdar, M.; Pluharova, E.; Mason, P. E.; Vacha, R.; Jungwirth, P. Ions at Hydrophobic Aqueous Interfaces: Molecular Dynamics with Effective Polarization. J. Phys. Chem. Lett. 2012, 3, 2087−2091.

I

DOI: 10.1021/acs.jpcb.5b10755 J. Phys. Chem. B XXXX, XXX, XXX−XXX