THE JOURNAL OF CHEMICAL PHYSICS 142, 014504 (2015)
Sensing polarization effects through the analysis of the effective C6 dispersion coefficients in NaCl solutions Alfredo Guevara-García,a) Joel Ireta, and Marcelo Galván Departamento de Química, División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, A.P. 55-534, México D.F. 09340, México
(Received 24 October 2014; accepted 16 December 2014; published online 7 January 2015) Density functional theory based ab initio molecular dynamics is used to obtain microscopic details of the interactions in sodium chloride solutions. By following the changes in the atomic C6 coefficients under the Tkatchenko-Scheffler’s scheme, we were able to identify two different coordination situations for the Cl− ion with significant different capabilities to perform dispersion interactions. This capability is enhanced when the ion-ion distance corresponds to the contact ion-pair situation. Also, the oxygen and hydrogen atoms of the water molecules change their aptitudes to interact through van der Waals like terms when they are close to the cation region of the ion-pair. These results have interesting implications on the design of force fields to model electrolyte solutions. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905125] I. INTRODUCTION
Sodium chloride aqueous solutions are one of the classic textbook examples for strong electrolytes, i.e., the ions are considered to be completely dissociated and solvated by water molecules. Although we might expect this behavior to be true over a wide range of concentrations, there is experimental and theoretical evidence of ion-pairs and clusters formation at concentrations well below the saturation limit.1,2 The formation of ion-pairs and clusters is the origin of the differences observed between real and ideal solutions. On the other hand, and because of its importance in different branches of knowledge such as electrochemistry, biochemistry, and atmospheric sciences, there has been a considerable effort to develop accurate force fields for using them in molecular dynamics (MD) simulations to describe the behavior of inorganic salts in aqueous solution.3–6 Despite this effort, up to now none of the force fields at hand are able to describe, among other properties of electrolytes, the correct solubility at all concentrations. In this context, knowing the detailed microscopic behavior of ions and water molecules in solution is of crucial importance in order to understand ion pairing and clusters formation and obtain new guidelines to improve classical force fields; ab initio MD simulations based on Density Functional Theory (DFT) calculations are a good choice to achieve this goal. In spite of the fact that van der Waals (vdW) interactions are tiny, as compared to their charge-charge, ion-dipole, and dipole-dipole counterparts, they are important in the accurate description of ionic crystals.7 So, as vdW interactions might play an important role in the correct description of electrolyte aqueous solutions, it is important to include such effects to achieve their precise description. In particular, for the purposes of this work, the Tkatchenko-Scheffler’s (TS) method8 was used to bypass the known inability of commonly used density functional approximations to describe vdW interactions.9 This method includes a)Electronic mail: [email protected]
0021-9606/2015/142(1)/014504/5/$30.00
the vdW interactions by adding a pairwise interatomic C6 R−6 term to the DFT energy 1 eff −6 f damp(RAB,Reff ETotal = EDFT − (1) A ,RB )C6AB RAB, 2 A, B where RAB is the distance between atoms A and B, and f damp is a damping function that is described below. The C6AB coefficients are computed as C6AB =
eff eff 2C6A C6B eff
αB
eff
αA
αA
αB
eff eff C6A +
,
(2)
eff eff C6B
eff where α eff i and C6i are the effective static polarizabilities and atomic C6 coefficients inside a specific chemical environment.
The advantage of this specific method is that the C6AB coefficients do capture the effects of the environment surrounding each one of the involved atoms. This is done by rescaling the free-atom reference values of the static polarizability (Eq. (3)) Vieff * = free + α free i , Vi and the C6 atomic coefficients (Eq. (4)) α eff i
(3)
2
V eff free = * ifree + C6i (4) V , i by a factor that depends on the ratio of the effective volume that an atom occupies in the actual environment (Vieff ) and the volume of the isolated one (Vifree). The effective volume of an atom in the system is computed as proposed by Johnson and Becke10 eff Vi = r 3ωi (r)ρ(r)d 3r, (5) eff C6i
where ρ(r) is the total density of the system, r is the distance from the nucleus of atom i, and ωi is the Hirshfeld atomic
142, 014504-1
© 2015 AIP Publishing LLC
014504-2
Guevara-García, Ireta, and Galván
J. Chem. Phys. 142, 014504 (2015)
partitioning weight11 of atom i,
TABLE I. Parameters used in the TS method. free
ρfree (r) ωi (r) = i free . ρ j (r)
(6)
The damping function in Eq. (1) has the form eff f damp(RAB,Reff A ,RB ) =
1
) ,
(
1 + exp −d
R AB eff
eff
s R (R A +R B )
(7)
H O Na Cl
free
free
C 6i (a.u.)
Ri (a.u.)
αi (a.u.)
Reference state
6.500 15.600 1.670 351.000
3.099 3.194 2.494 3.855
4.500 5.400 1.060 36.997
0 0 1 −1
−1
eff where Reff A and RB are also rescaled values of the vdW radii of the free atoms 1/3
V eff (8) = * ifree + Rifree. , Vi In this article, we present the results of ab initio MD simulations for a sodium chloride solution using the DFTTS method. One of the interesting features of this method is that the effective C6 coefficients are feelers for the changes in the coordination environment of each atomic site. Johnson has shown,12 using the exchange-hole dipole moment model, that the dispersion coefficients are quite sensitive to the atomic environment. Along these lines, our analysis focuses on the induced changes in the Ceff 6 coefficients by the environments embodied in the ionic solution. Rieff
II. COMPUTATIONAL DETAILS
The TS method predicts too small lattice constants and too large cohesive energies when treating ionic systems.13 Buˇcko et al.14 analyzed this problem and realized that this failure was mainly produced by the arbitrary choice of the neutral atomic reference state for the Hirshfeld partition scheme, used in the original implementation of the TS method. They observed that the computed properties of ionic solids improved just by using the noninteracting ions as references in the original Hirshfeld partition involved in the TS method. As they mentioned the applicability of that approach is not general, so they proposed to use the Hirshfeld Iterative scheme15 instead of the standard Hirshfeld partitioning, in order to compute the relative volumes, obtaining a significant improvement in the calculated properties and keeping a moderate computational cost. In order to avoid this increase in the computational cost, we decided to use the noninteracting ions as references (for Na and Cl) in the original TS method knowing that the atoms in our system under study will not suffer large oxidation state changes. We used the atomic polarizabilities and the atomic C6 coefficients from Refs. 16 and 17. The Rfree for Na+ was i 8 computed using the protocol proposed in the TS method, i.e., the Rfree corresponds to the radius, where the electron i density value of Na+ equals the electron density value at the neon’s vdW radius. The Rfree for Cl− was then adjusted to i reproduce the NaCl lattice constant (5.63 Å)18 using the TS method in conjunction with the parameters already mentioned. The TS parameters used in this work are listed in Table I. The cohesive energy and bulk modulus for NaCl using this new set of parameters are 3.29 eV/atom and 30.59 GPa, respectively (experimental values: 3.31 eV/atom19 and 28 GPa20). This is
a significant improvement with respect to the values obtained with the standard parameters13 (3.97 eV/atom and 84 GPa). All the calculations were performed using the PerdewBurke-Ernzerhof (PBE)21 exchange-correlation functional in combination with the TS method (PBE-TS) as implemented in VASP 5.3.3 code13,22,23 and using the projector-augmented wave method.24 PBE functional plus dispersion corrections have been previously used to study liquid water25 and NaCl solutions.26 In these works,25,26 it has been observed that the inclusion of vdW corrections improves the description of these systems. Our model system contains 30 water molecules and a NaCl pair in a cubic cell of a = 9.9459 Å with periodic boundary conditions. This high concentration (1.69M) was used to promote ion-pairing events along the MD simulation. Only one k-point was used (Gamma point) and the energy cutoff was 700 eV. This cutoff guarantees convergence of the energy within 0.01 eV between different configurations. The last statement is important because the changes in the total vdW contributions are of this order of magnitude. The time step for the ab initio MD was 0.25 fs at 300 K. The initial structure for the molecular dynamics run was obtained from an equilibrated classical MD simulation. We relaxed the final configuration and performed a 20 ps MD using Grimme’s PBE-D2 method27 over this system. The final coordinates and velocities were used as the initial input for the PBE-TS MD using the cation and the anion reference densities for Na and Cl, respectively. The production trajectory thus obtained was 32.5 ps long.
III. RESULTS AND DISCUSSION
The data extracted from the analysis of our PBE-TS MD are in good agreement with neutron diffraction experiments28 of NaCl solutions. For example, the position of the first peak of the NaO and ClH radial distribution functions g(r), calculated using the whole production trajectory, are 2.32 and 2.18 Å, respectively. These values compare pretty well to the experimental ones, 2.34 and 2.19 Å, observed in a solution with a solute:solvent ratio of 1:83. A recent Car-Parrinello MD study26 using the Grimme-D2 method, that does not take into account the environment when computing the vdW contribution to the total energy, predicts peaks at 2.43 and 2.15 Å. Besides the possibility of including dispersion corrections to the standard GGA functionals, it has been pointed out that the effective atomic C6 coefficients (Ceff 6 ) obtained within the TS method can be used to monitor the impact of chemical environments on atomic species.29 This feature is particularly useful in this work, as one can see in Figure 1, where we
014504-3
Guevara-García, Ireta, and Galván
eff
FIG. 1. C6 (Cl) values in atomic units along a 32.5 ps molecular dynamics trajectory. The distance between the chloride and the closest sodium ion is also shown as measured within the minimum image convention. The SSIP label indicates the representative interval in which the ion pair is separated by around 5 Å, allowing solvent molecules reside between the ions. The CIP label indicates the situations in which the ions remain at a distance around 3 Å.
eff show the values of the Ceff 6 for the chloride ion, C6 (Cl), and the distance between the chloride ion and the nearest sodium ion along the production trajectory (130 000 configurations), obtained from the PBE-TS MD simulation. It can be noticed that the Ceff 6 of the chloride ion take different values depending on the distance to the sodium ion. We analyzed two intervals of the MD trajectory (see Figure 1), in order to obtain average values that clearly represent different configurations: the solvent separated ion pair (SSIP) structure having a Ceff 6 (Cl) average value of 207.60 a.u. and a mean distance between Na+ and Cl− of 4.97 Å, and the contact-ion pair(CIP) conformation, with average values of 227.64 a.u. and 2.91 Å respectively. This contrast in the behavior of the Ceff 6 (Cl) coefficient is clearly related to the polarization induced by the Na+ ion. One may conclude, from Eq. (1) and Figure 1, that different chemical environments for Cl− imply, through the value of Ceff 6 (Cl), a different contribution of the dispersion correction terms coming from the surroundings of this ion. If one tries to recast this finding in the classical MD language, one may say that the interaction pair potential between the chloride ion and its neighbors in the CIP configuration is different from the one in the SSIP situation. To verify that this finding was not an artifact of the simulation box size, we built a larger system containing 3 NaCl pairs and 123 water molecules in a cubic cell of 15.52 Å and ran a classical MD. We picked a structure that contained ions in the CIP and SSIP configurations simultaneously and optimized it using the PBE-TS method. The resulting Ceff 6 for the CIP and SSIP chlorides follow the same trend observed in the smaller system. It is worth mentioning that changes of the same magnitude in classical MD parameters can produce very different results. eff For the case of sodium ion’s Ceff 6 , C6 (Na), the difference between the CIP and SSIP configurations is smaller than for the chloride ion; the average value in the SSIP region is 1.83 a.u. and 1.79 a.u. in the CIP region. This observation is in agreement with what we could expect from the difference
J. Chem. Phys. 142, 014504 (2015)
FIG. 2. Average values of the effective C6 coefficients for the oxygen atom in eff atomic units, C6 (O)(a.u.), as a function of its distance to the closest ion. The average is displayed for the SSIP and for the CIP regimes (see Figure 1). This average is evaluated by considering all oxygen atoms at the given distance regardless if they are closer to the Na+ or Cl− ions; so, at a given distance the average may include oxygen atoms close to the cation and close to the anion. The radial pair distribution functions, g(r), for Na-O and Cl-O are also shown. The C6 reference value for oxygen atom used in the Tkatchenko-Scheffler method is 15.6 a.u.
in polarizabilities between the cation and the anion and their behavior in solution.30 As can be seen in Figure 2, the values for the oxygen’s Ceff 6 , Ceff (O), decay as a function of the distance to the closest ion 6 (Na+ or Cl−) up to what can be considered as the “bulk water value” of around 11.9 a.u. This limit value, as compared to the reference of 15.6 a.u. of the isolated oxygen atom, represents a decrease of the order of 24%, indicating an important change in the Ceff 6 (O) coefficient induced by the bulk water environment. The general behavior of Ceff 6 (O) as a function of the distance to the closest ion is similar for both regimes of separation of the ions, SSIP and CIP: when the water molecule is approaching the ion pair, its oxygen atom keeps the value of the Ceff 6 (O) around the one of bulk water up to a distance of 3 Å. For closer distances, the value increases between 13% and 20% with respect to the bulk like value. These results point towards the idea that the changes in Ceff 6 (O) depend mainly on the distance to any of the two ions rather than the distance between them; a similar behavior in the dipole moment of hydration water molecules was observed by Timko et al.31 Comparing the NaO and ClO radial pair distribution functions with the + average values for Ceff 6 (O) (Figure 2), the impact of the Na eff on C6 (O) is clearly illustrated for distances between 2 and 3 Å. The comparison also makes evident that the changes in Ceff 6 (O) for distances between 2.8 and 3.6 Å are connected to the interaction with Cl−. The small differences between the SSIP and CIP curves are related to changes induced in the ion pair environment by the interaction of the ions: when the two ions are close, there is a slightly increase in the Ceff 6 value of the oxygen atoms surrounding the chloride ion and a small decrease in the same quantity for the oxygen atoms in the sodium ion environment. These effects can be rationalized as the impact of one of the ions in the way the other charged species polarizes the water molecules in its neighborhood.
014504-4
Guevara-García, Ireta, and Galván
J. Chem. Phys. 142, 014504 (2015) TABLE II. Coordination numbers for Na+ and Cl−. Computational details (method, temperature (T in K), and concentration (C in M)) are also shown. Experimental uncertainties are reported in parenthesis.
Method
C (M)
T (K)
NaO
ClO
ClH
PBE033 PBE33 PBE26 PBE-D226 PBE-TSa PBE-TS SSIP PBE-TS CIP Expt.28
1.03 0.87 1.3 1.3 1.69 1.69 1.69 0.66
380 400 315 315 300 300 300 298
5.3 5.3 4.9 5.8 4.3 4.9 3.8 5.3(0.8)
6.4 5.9 5.9 6.2 5.3 5.4 4.6 6.9(1.0)
5.2 5.1 4.7 5.3 4.0 6.0(1.1)
a Considering
FIG. 3. Average effective C6 coefficient for the hydrogen atoms in atomic eff units, C6 (H)(a.u.), as a function of the distance between the H atom and the closest ion. The average is displayed for the SSIP and for the CIP regimes (see Figure 1). This average is evaluated by considering all hydrogen atoms at the given distance regardless if they are closer to the Na+ or Cl− ions; so, at a given distance the average may include hydrogen atoms close to the cation and close to the anion. The radial pair distribution functions, g(r), for Na-H and Cl-H are also shown. The C6 reference value for the hydrogen atom used in the Tkatchenko-Scheffler method is 6.5 a.u.
Figure 3 shows the average values for the hydrogen’s Ceff 6 , (H), together with the NaH and ClH radial pair distributions functions. There, it is evidenced that the hydrogen’s Ceff 6 values depend only on how close the hydrogen atoms are to the ions. The greatest changes in the hydrogen’s Ceff 6 occur for those atoms that belong to the sodium’s first water coordination shell. Again, the chloride ion has a smaller influence on the Ceff 6 (H) than the sodium ion. For this case, the change between the isolated reference value and the bulk like one is of 46%; and the largest change as the water molecule approach the ion pair is of 13%. Very recently, DiStasio et al.32 have shown that inclusion of exact exchange effects and vdW interactions are important to describe liquid water. We evaluated the Ceff 6 for the selected snapshots of our MD trajectory using the PBE0TS functional and observed that the behavior of the Ceff 6 follows the same trend as in the case of the PBE-TS method. With the aim of looking for differences between the SSIP and CIP regimes, the g(r) for both regions of the simulation were analyzed for the ClH and NaO pairs. The major difference between the SSIP and the CIP is the height of the first peak. The larger height in the SSIP reflects the greater contact surface area to interact with water molecules. There is not an appreciable change on the position of the peaks for all the above cases. Another common quantity computed to characterize NaCl solutions is the coordination number of the ions. Table II shows the coordination numbers computed in this work, the values of other selected works and the experimental ones. Several important features can be observed from Table II: in general, coordination numbers are very sensitive to the functional and temperature used during the simulation.26,33 It has been observed experimentally28 and in simulations4 that the coordination numbers decrease when the solute concentration increases. Our computed coordination numbers follow this tendency and fall in the experimental error window.28 As expected, the coordination numbers evaluated in the CIP region Ceff 6
the whole trajectory.
are smaller than the ones evaluated in the SSIP region. It has to be highlighted that the number of CIP and SSIP configurations observed in a specific solution depends on the concentration but the main polarization effects observed in this work should be less sensitive to this variable.
IV. CONCLUSIONS
Exploring and testing the idea of using Ceff 6 as a tool for analyzing the impact of the environment in the atomic regions of a typical ion pair give rise to several points to notice. On − the one hand, the significant change of the Ceff 6 of Cl when the cation is close to it indicates a clear influence of the cation on the capabilities of the anion to interact with its neighbors through dispersion like interactions. On the other hand, there are changes in Ceff 6 of the hydrogen and oxygen atomic regions in the water molecules of the first solvation sphere, mainly induced by the cation region of the ion pair environment. It is important to notice that the main changes induced by the cation environment in the C6 coefficients imply an increase of the corresponding dispersion attraction terms. Our findings suggest that, in the case of classical potentials, a polarizable water force field might be essential to describe the ion pair solvation. In addition, the modification induced by the cation on the anion in terms of the capabilities to perform dispersion like interactions should be taken into account; both observations could have impact in the development of force fields to address electrolyte systems. Whether these effects are a possible explanation of why classical force fields, that use a fixed set of parameters, in some cases favor hydration over ion association and vice versa,3,4 or if they are related to the high sensitivity of the NaCl nucleation mechanism to small changes in the force field, as has been observed within the context of classical MD,3 has to be tested but is beyond the objectives of this work. ACKNOWLEDGMENTS
We thank José Alejandre for useful comments and discussion and for providing us with the classical molecular dynamics equilibrated systems. This work was partially supported by CONACYT through the Grant Nos. 155698 and 154784. A.G. thanks CONACYT for a postdoctoral fellow-
014504-5
Guevara-García, Ireta, and Galván
ship. The authors thank the Laboratorio de Supercómputo y Visualización en Paralelo at UAM-I for computational resources. 1Y.
Marcus and G. Hefter, Chem. Rev. 106, 4585 (2006). Degrève and F. L. B. da Silva, J. Mol. Liq. 87, 217 (2000). 3J. Alejandre and J.-P. Hansen, Phys. Rev. E 76, 061505 (2007). 4J. Alejandre, G. A. Chapela, F. Bresme, and J.-P. Hansen, J. Chem. Phys. 130, 174505 (2009). 5F. N. Mendoza and J. Alejandre, J. Mol. Liq. 185, 50 (2013). 6F. Mouˇ cka, I. Nezbeda, and W. R. Smith, J. Chem. Theory Comput. 9, 5076 (2013). 7G.-X. Zhang, A. Tkatchenko, J. Paier, H. Appel, and M. Scheffler, Phys. Rev. Lett. 107, 245501 (2011). 8A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 9J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012). 10E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006). 11F. L. Hirshfeld, Theor. Chim. Acta 44, 129 (1977). 12E. R. Johnson, J. Chem. Phys. 135, 234109 (2011). 13T. Buˇ cko, S. Lebègue, J. Hafner, and J. G. Ángyán, Phys. Rev. B 87, 064110 (2013). 14T. Buˇ cko, S. Lebègue, J. Hafner, and J. G. Ángyán, J. Chem. Theory Comput. 9, 4293 (2013). 15P. Bultinck, C. Van Alsenoy, P. W. Ayers, and R. Carbó-Dorca, J. Chem. Phys. 126, 144111 (2007). 16G. D. Mahan, J. Chem. Phys. 76, 493 (1982). 2L.
J. Chem. Phys. 142, 014504 (2015) 17C.
Hättig and B.A. Heß, J. Chem. Phys. 108, 3863 (1998). W. C. Boswell, Proc. Phys. Soc. Sect. A 64, 465 (1951). 19M. W. Chase, Jr., NIST-JANAF Thermochemical Tables (AIP, New York, 1998). 20H. Spetzler, C. Sammis, and R. O’Connell, J. Phys. Chem. Solids 33, 1727 (1972). 21J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 22G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996). 23G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 24P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). 25J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter, and C. J. Mundy, J. Phys. Chem. B 113, 11959 (2009). 26A. Bankura, V. Carnevale, and M. L. Klein, J. Chem. Phys. 138, 014501 (2013). 27S. Grimme, J. Comput. Chem. 27, 1787 (2006). 28R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, J. Phys. Chem. B 111, 13570 (2007). 29M. Rossi, A. Tkatchenko, S. B. R. Rempe, and S. V. Varma, Proc. Natl. Acad. Sci. U. S. A. 110, 12978 (2013). 30J. J. Molina, S. Lectez, S. Tazi, M. Salanne, J.-F. Dufrêche, J. Roques, E. Simoni, P. A. Madden, and P. Turq, J. Chem. Phys. 134, 014511 (2011). 31J. Timko, D. Bucher, and S. Kuyucak, J. Chem. Phys. 132, 114510 (2010). 32R. A. DiStasio, B. Santra, Z. Li, X. Wu, and R. Car, J. Chem. Phys. 141, 084502 (2014). 33A. P. Gaiduk, C. Zhang, F. Gygi, and G. Galli, Chem. Phys. Lett. 604, 89 (2014).
18F.
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
Sensing polarization effects through the analysis of the effective C6 dispersion coefficients in NaCl solutions Alfredo Guevara-García,a) Joel Ireta, and Marcelo Galván Departamento de Química, División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, A.P. 55-534, México D.F. 09340, México
(Received 24 October 2014; accepted 16 December 2014; published online 7 January 2015) Density functional theory based ab initio molecular dynamics is used to obtain microscopic details of the interactions in sodium chloride solutions. By following the changes in the atomic C6 coefficients under the Tkatchenko-Scheffler’s scheme, we were able to identify two different coordination situations for the Cl− ion with significant different capabilities to perform dispersion interactions. This capability is enhanced when the ion-ion distance corresponds to the contact ion-pair situation. Also, the oxygen and hydrogen atoms of the water molecules change their aptitudes to interact through van der Waals like terms when they are close to the cation region of the ion-pair. These results have interesting implications on the design of force fields to model electrolyte solutions. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905125] I. INTRODUCTION
Sodium chloride aqueous solutions are one of the classic textbook examples for strong electrolytes, i.e., the ions are considered to be completely dissociated and solvated by water molecules. Although we might expect this behavior to be true over a wide range of concentrations, there is experimental and theoretical evidence of ion-pairs and clusters formation at concentrations well below the saturation limit.1,2 The formation of ion-pairs and clusters is the origin of the differences observed between real and ideal solutions. On the other hand, and because of its importance in different branches of knowledge such as electrochemistry, biochemistry, and atmospheric sciences, there has been a considerable effort to develop accurate force fields for using them in molecular dynamics (MD) simulations to describe the behavior of inorganic salts in aqueous solution.3–6 Despite this effort, up to now none of the force fields at hand are able to describe, among other properties of electrolytes, the correct solubility at all concentrations. In this context, knowing the detailed microscopic behavior of ions and water molecules in solution is of crucial importance in order to understand ion pairing and clusters formation and obtain new guidelines to improve classical force fields; ab initio MD simulations based on Density Functional Theory (DFT) calculations are a good choice to achieve this goal. In spite of the fact that van der Waals (vdW) interactions are tiny, as compared to their charge-charge, ion-dipole, and dipole-dipole counterparts, they are important in the accurate description of ionic crystals.7 So, as vdW interactions might play an important role in the correct description of electrolyte aqueous solutions, it is important to include such effects to achieve their precise description. In particular, for the purposes of this work, the Tkatchenko-Scheffler’s (TS) method8 was used to bypass the known inability of commonly used density functional approximations to describe vdW interactions.9 This method includes a)Electronic mail: [email protected]
0021-9606/2015/142(1)/014504/5/$30.00
the vdW interactions by adding a pairwise interatomic C6 R−6 term to the DFT energy 1 eff −6 f damp(RAB,Reff ETotal = EDFT − (1) A ,RB )C6AB RAB, 2 A, B where RAB is the distance between atoms A and B, and f damp is a damping function that is described below. The C6AB coefficients are computed as C6AB =
eff eff 2C6A C6B eff
αB
eff
αA
αA
αB
eff eff C6A +
,
(2)
eff eff C6B
eff where α eff i and C6i are the effective static polarizabilities and atomic C6 coefficients inside a specific chemical environment.
The advantage of this specific method is that the C6AB coefficients do capture the effects of the environment surrounding each one of the involved atoms. This is done by rescaling the free-atom reference values of the static polarizability (Eq. (3)) Vieff * = free + α free i , Vi and the C6 atomic coefficients (Eq. (4)) α eff i
(3)
2
V eff free = * ifree + C6i (4) V , i by a factor that depends on the ratio of the effective volume that an atom occupies in the actual environment (Vieff ) and the volume of the isolated one (Vifree). The effective volume of an atom in the system is computed as proposed by Johnson and Becke10 eff Vi = r 3ωi (r)ρ(r)d 3r, (5) eff C6i
where ρ(r) is the total density of the system, r is the distance from the nucleus of atom i, and ωi is the Hirshfeld atomic
142, 014504-1
© 2015 AIP Publishing LLC
014504-2
Guevara-García, Ireta, and Galván
J. Chem. Phys. 142, 014504 (2015)
partitioning weight11 of atom i,
TABLE I. Parameters used in the TS method. free
ρfree (r) ωi (r) = i free . ρ j (r)
(6)
The damping function in Eq. (1) has the form eff f damp(RAB,Reff A ,RB ) =
1
) ,
(
1 + exp −d
R AB eff
eff
s R (R A +R B )
(7)
H O Na Cl
free
free
C 6i (a.u.)
Ri (a.u.)
αi (a.u.)
Reference state
6.500 15.600 1.670 351.000
3.099 3.194 2.494 3.855
4.500 5.400 1.060 36.997
0 0 1 −1
−1
eff where Reff A and RB are also rescaled values of the vdW radii of the free atoms 1/3
V eff (8) = * ifree + Rifree. , Vi In this article, we present the results of ab initio MD simulations for a sodium chloride solution using the DFTTS method. One of the interesting features of this method is that the effective C6 coefficients are feelers for the changes in the coordination environment of each atomic site. Johnson has shown,12 using the exchange-hole dipole moment model, that the dispersion coefficients are quite sensitive to the atomic environment. Along these lines, our analysis focuses on the induced changes in the Ceff 6 coefficients by the environments embodied in the ionic solution. Rieff
II. COMPUTATIONAL DETAILS
The TS method predicts too small lattice constants and too large cohesive energies when treating ionic systems.13 Buˇcko et al.14 analyzed this problem and realized that this failure was mainly produced by the arbitrary choice of the neutral atomic reference state for the Hirshfeld partition scheme, used in the original implementation of the TS method. They observed that the computed properties of ionic solids improved just by using the noninteracting ions as references in the original Hirshfeld partition involved in the TS method. As they mentioned the applicability of that approach is not general, so they proposed to use the Hirshfeld Iterative scheme15 instead of the standard Hirshfeld partitioning, in order to compute the relative volumes, obtaining a significant improvement in the calculated properties and keeping a moderate computational cost. In order to avoid this increase in the computational cost, we decided to use the noninteracting ions as references (for Na and Cl) in the original TS method knowing that the atoms in our system under study will not suffer large oxidation state changes. We used the atomic polarizabilities and the atomic C6 coefficients from Refs. 16 and 17. The Rfree for Na+ was i 8 computed using the protocol proposed in the TS method, i.e., the Rfree corresponds to the radius, where the electron i density value of Na+ equals the electron density value at the neon’s vdW radius. The Rfree for Cl− was then adjusted to i reproduce the NaCl lattice constant (5.63 Å)18 using the TS method in conjunction with the parameters already mentioned. The TS parameters used in this work are listed in Table I. The cohesive energy and bulk modulus for NaCl using this new set of parameters are 3.29 eV/atom and 30.59 GPa, respectively (experimental values: 3.31 eV/atom19 and 28 GPa20). This is
a significant improvement with respect to the values obtained with the standard parameters13 (3.97 eV/atom and 84 GPa). All the calculations were performed using the PerdewBurke-Ernzerhof (PBE)21 exchange-correlation functional in combination with the TS method (PBE-TS) as implemented in VASP 5.3.3 code13,22,23 and using the projector-augmented wave method.24 PBE functional plus dispersion corrections have been previously used to study liquid water25 and NaCl solutions.26 In these works,25,26 it has been observed that the inclusion of vdW corrections improves the description of these systems. Our model system contains 30 water molecules and a NaCl pair in a cubic cell of a = 9.9459 Å with periodic boundary conditions. This high concentration (1.69M) was used to promote ion-pairing events along the MD simulation. Only one k-point was used (Gamma point) and the energy cutoff was 700 eV. This cutoff guarantees convergence of the energy within 0.01 eV between different configurations. The last statement is important because the changes in the total vdW contributions are of this order of magnitude. The time step for the ab initio MD was 0.25 fs at 300 K. The initial structure for the molecular dynamics run was obtained from an equilibrated classical MD simulation. We relaxed the final configuration and performed a 20 ps MD using Grimme’s PBE-D2 method27 over this system. The final coordinates and velocities were used as the initial input for the PBE-TS MD using the cation and the anion reference densities for Na and Cl, respectively. The production trajectory thus obtained was 32.5 ps long.
III. RESULTS AND DISCUSSION
The data extracted from the analysis of our PBE-TS MD are in good agreement with neutron diffraction experiments28 of NaCl solutions. For example, the position of the first peak of the NaO and ClH radial distribution functions g(r), calculated using the whole production trajectory, are 2.32 and 2.18 Å, respectively. These values compare pretty well to the experimental ones, 2.34 and 2.19 Å, observed in a solution with a solute:solvent ratio of 1:83. A recent Car-Parrinello MD study26 using the Grimme-D2 method, that does not take into account the environment when computing the vdW contribution to the total energy, predicts peaks at 2.43 and 2.15 Å. Besides the possibility of including dispersion corrections to the standard GGA functionals, it has been pointed out that the effective atomic C6 coefficients (Ceff 6 ) obtained within the TS method can be used to monitor the impact of chemical environments on atomic species.29 This feature is particularly useful in this work, as one can see in Figure 1, where we
014504-3
Guevara-García, Ireta, and Galván
eff
FIG. 1. C6 (Cl) values in atomic units along a 32.5 ps molecular dynamics trajectory. The distance between the chloride and the closest sodium ion is also shown as measured within the minimum image convention. The SSIP label indicates the representative interval in which the ion pair is separated by around 5 Å, allowing solvent molecules reside between the ions. The CIP label indicates the situations in which the ions remain at a distance around 3 Å.
eff show the values of the Ceff 6 for the chloride ion, C6 (Cl), and the distance between the chloride ion and the nearest sodium ion along the production trajectory (130 000 configurations), obtained from the PBE-TS MD simulation. It can be noticed that the Ceff 6 of the chloride ion take different values depending on the distance to the sodium ion. We analyzed two intervals of the MD trajectory (see Figure 1), in order to obtain average values that clearly represent different configurations: the solvent separated ion pair (SSIP) structure having a Ceff 6 (Cl) average value of 207.60 a.u. and a mean distance between Na+ and Cl− of 4.97 Å, and the contact-ion pair(CIP) conformation, with average values of 227.64 a.u. and 2.91 Å respectively. This contrast in the behavior of the Ceff 6 (Cl) coefficient is clearly related to the polarization induced by the Na+ ion. One may conclude, from Eq. (1) and Figure 1, that different chemical environments for Cl− imply, through the value of Ceff 6 (Cl), a different contribution of the dispersion correction terms coming from the surroundings of this ion. If one tries to recast this finding in the classical MD language, one may say that the interaction pair potential between the chloride ion and its neighbors in the CIP configuration is different from the one in the SSIP situation. To verify that this finding was not an artifact of the simulation box size, we built a larger system containing 3 NaCl pairs and 123 water molecules in a cubic cell of 15.52 Å and ran a classical MD. We picked a structure that contained ions in the CIP and SSIP configurations simultaneously and optimized it using the PBE-TS method. The resulting Ceff 6 for the CIP and SSIP chlorides follow the same trend observed in the smaller system. It is worth mentioning that changes of the same magnitude in classical MD parameters can produce very different results. eff For the case of sodium ion’s Ceff 6 , C6 (Na), the difference between the CIP and SSIP configurations is smaller than for the chloride ion; the average value in the SSIP region is 1.83 a.u. and 1.79 a.u. in the CIP region. This observation is in agreement with what we could expect from the difference
J. Chem. Phys. 142, 014504 (2015)
FIG. 2. Average values of the effective C6 coefficients for the oxygen atom in eff atomic units, C6 (O)(a.u.), as a function of its distance to the closest ion. The average is displayed for the SSIP and for the CIP regimes (see Figure 1). This average is evaluated by considering all oxygen atoms at the given distance regardless if they are closer to the Na+ or Cl− ions; so, at a given distance the average may include oxygen atoms close to the cation and close to the anion. The radial pair distribution functions, g(r), for Na-O and Cl-O are also shown. The C6 reference value for oxygen atom used in the Tkatchenko-Scheffler method is 15.6 a.u.
in polarizabilities between the cation and the anion and their behavior in solution.30 As can be seen in Figure 2, the values for the oxygen’s Ceff 6 , Ceff (O), decay as a function of the distance to the closest ion 6 (Na+ or Cl−) up to what can be considered as the “bulk water value” of around 11.9 a.u. This limit value, as compared to the reference of 15.6 a.u. of the isolated oxygen atom, represents a decrease of the order of 24%, indicating an important change in the Ceff 6 (O) coefficient induced by the bulk water environment. The general behavior of Ceff 6 (O) as a function of the distance to the closest ion is similar for both regimes of separation of the ions, SSIP and CIP: when the water molecule is approaching the ion pair, its oxygen atom keeps the value of the Ceff 6 (O) around the one of bulk water up to a distance of 3 Å. For closer distances, the value increases between 13% and 20% with respect to the bulk like value. These results point towards the idea that the changes in Ceff 6 (O) depend mainly on the distance to any of the two ions rather than the distance between them; a similar behavior in the dipole moment of hydration water molecules was observed by Timko et al.31 Comparing the NaO and ClO radial pair distribution functions with the + average values for Ceff 6 (O) (Figure 2), the impact of the Na eff on C6 (O) is clearly illustrated for distances between 2 and 3 Å. The comparison also makes evident that the changes in Ceff 6 (O) for distances between 2.8 and 3.6 Å are connected to the interaction with Cl−. The small differences between the SSIP and CIP curves are related to changes induced in the ion pair environment by the interaction of the ions: when the two ions are close, there is a slightly increase in the Ceff 6 value of the oxygen atoms surrounding the chloride ion and a small decrease in the same quantity for the oxygen atoms in the sodium ion environment. These effects can be rationalized as the impact of one of the ions in the way the other charged species polarizes the water molecules in its neighborhood.
014504-4
Guevara-García, Ireta, and Galván
J. Chem. Phys. 142, 014504 (2015) TABLE II. Coordination numbers for Na+ and Cl−. Computational details (method, temperature (T in K), and concentration (C in M)) are also shown. Experimental uncertainties are reported in parenthesis.
Method
C (M)
T (K)
NaO
ClO
ClH
PBE033 PBE33 PBE26 PBE-D226 PBE-TSa PBE-TS SSIP PBE-TS CIP Expt.28
1.03 0.87 1.3 1.3 1.69 1.69 1.69 0.66
380 400 315 315 300 300 300 298
5.3 5.3 4.9 5.8 4.3 4.9 3.8 5.3(0.8)
6.4 5.9 5.9 6.2 5.3 5.4 4.6 6.9(1.0)
5.2 5.1 4.7 5.3 4.0 6.0(1.1)
a Considering
FIG. 3. Average effective C6 coefficient for the hydrogen atoms in atomic eff units, C6 (H)(a.u.), as a function of the distance between the H atom and the closest ion. The average is displayed for the SSIP and for the CIP regimes (see Figure 1). This average is evaluated by considering all hydrogen atoms at the given distance regardless if they are closer to the Na+ or Cl− ions; so, at a given distance the average may include hydrogen atoms close to the cation and close to the anion. The radial pair distribution functions, g(r), for Na-H and Cl-H are also shown. The C6 reference value for the hydrogen atom used in the Tkatchenko-Scheffler method is 6.5 a.u.
Figure 3 shows the average values for the hydrogen’s Ceff 6 , (H), together with the NaH and ClH radial pair distributions functions. There, it is evidenced that the hydrogen’s Ceff 6 values depend only on how close the hydrogen atoms are to the ions. The greatest changes in the hydrogen’s Ceff 6 occur for those atoms that belong to the sodium’s first water coordination shell. Again, the chloride ion has a smaller influence on the Ceff 6 (H) than the sodium ion. For this case, the change between the isolated reference value and the bulk like one is of 46%; and the largest change as the water molecule approach the ion pair is of 13%. Very recently, DiStasio et al.32 have shown that inclusion of exact exchange effects and vdW interactions are important to describe liquid water. We evaluated the Ceff 6 for the selected snapshots of our MD trajectory using the PBE0TS functional and observed that the behavior of the Ceff 6 follows the same trend as in the case of the PBE-TS method. With the aim of looking for differences between the SSIP and CIP regimes, the g(r) for both regions of the simulation were analyzed for the ClH and NaO pairs. The major difference between the SSIP and the CIP is the height of the first peak. The larger height in the SSIP reflects the greater contact surface area to interact with water molecules. There is not an appreciable change on the position of the peaks for all the above cases. Another common quantity computed to characterize NaCl solutions is the coordination number of the ions. Table II shows the coordination numbers computed in this work, the values of other selected works and the experimental ones. Several important features can be observed from Table II: in general, coordination numbers are very sensitive to the functional and temperature used during the simulation.26,33 It has been observed experimentally28 and in simulations4 that the coordination numbers decrease when the solute concentration increases. Our computed coordination numbers follow this tendency and fall in the experimental error window.28 As expected, the coordination numbers evaluated in the CIP region Ceff 6
the whole trajectory.
are smaller than the ones evaluated in the SSIP region. It has to be highlighted that the number of CIP and SSIP configurations observed in a specific solution depends on the concentration but the main polarization effects observed in this work should be less sensitive to this variable.
IV. CONCLUSIONS
Exploring and testing the idea of using Ceff 6 as a tool for analyzing the impact of the environment in the atomic regions of a typical ion pair give rise to several points to notice. On − the one hand, the significant change of the Ceff 6 of Cl when the cation is close to it indicates a clear influence of the cation on the capabilities of the anion to interact with its neighbors through dispersion like interactions. On the other hand, there are changes in Ceff 6 of the hydrogen and oxygen atomic regions in the water molecules of the first solvation sphere, mainly induced by the cation region of the ion pair environment. It is important to notice that the main changes induced by the cation environment in the C6 coefficients imply an increase of the corresponding dispersion attraction terms. Our findings suggest that, in the case of classical potentials, a polarizable water force field might be essential to describe the ion pair solvation. In addition, the modification induced by the cation on the anion in terms of the capabilities to perform dispersion like interactions should be taken into account; both observations could have impact in the development of force fields to address electrolyte systems. Whether these effects are a possible explanation of why classical force fields, that use a fixed set of parameters, in some cases favor hydration over ion association and vice versa,3,4 or if they are related to the high sensitivity of the NaCl nucleation mechanism to small changes in the force field, as has been observed within the context of classical MD,3 has to be tested but is beyond the objectives of this work. ACKNOWLEDGMENTS
We thank José Alejandre for useful comments and discussion and for providing us with the classical molecular dynamics equilibrated systems. This work was partially supported by CONACYT through the Grant Nos. 155698 and 154784. A.G. thanks CONACYT for a postdoctoral fellow-
014504-5
Guevara-García, Ireta, and Galván
ship. The authors thank the Laboratorio de Supercómputo y Visualización en Paralelo at UAM-I for computational resources. 1Y.
Marcus and G. Hefter, Chem. Rev. 106, 4585 (2006). Degrève and F. L. B. da Silva, J. Mol. Liq. 87, 217 (2000). 3J. Alejandre and J.-P. Hansen, Phys. Rev. E 76, 061505 (2007). 4J. Alejandre, G. A. Chapela, F. Bresme, and J.-P. Hansen, J. Chem. Phys. 130, 174505 (2009). 5F. N. Mendoza and J. Alejandre, J. Mol. Liq. 185, 50 (2013). 6F. Mouˇ cka, I. Nezbeda, and W. R. Smith, J. Chem. Theory Comput. 9, 5076 (2013). 7G.-X. Zhang, A. Tkatchenko, J. Paier, H. Appel, and M. Scheffler, Phys. Rev. Lett. 107, 245501 (2011). 8A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 9J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012). 10E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006). 11F. L. Hirshfeld, Theor. Chim. Acta 44, 129 (1977). 12E. R. Johnson, J. Chem. Phys. 135, 234109 (2011). 13T. Buˇ cko, S. Lebègue, J. Hafner, and J. G. Ángyán, Phys. Rev. B 87, 064110 (2013). 14T. Buˇ cko, S. Lebègue, J. Hafner, and J. G. Ángyán, J. Chem. Theory Comput. 9, 4293 (2013). 15P. Bultinck, C. Van Alsenoy, P. W. Ayers, and R. Carbó-Dorca, J. Chem. Phys. 126, 144111 (2007). 16G. D. Mahan, J. Chem. Phys. 76, 493 (1982). 2L.
J. Chem. Phys. 142, 014504 (2015) 17C.
Hättig and B.A. Heß, J. Chem. Phys. 108, 3863 (1998). W. C. Boswell, Proc. Phys. Soc. Sect. A 64, 465 (1951). 19M. W. Chase, Jr., NIST-JANAF Thermochemical Tables (AIP, New York, 1998). 20H. Spetzler, C. Sammis, and R. O’Connell, J. Phys. Chem. Solids 33, 1727 (1972). 21J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 22G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996). 23G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 24P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). 25J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter, and C. J. Mundy, J. Phys. Chem. B 113, 11959 (2009). 26A. Bankura, V. Carnevale, and M. L. Klein, J. Chem. Phys. 138, 014501 (2013). 27S. Grimme, J. Comput. Chem. 27, 1787 (2006). 28R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, J. Phys. Chem. B 111, 13570 (2007). 29M. Rossi, A. Tkatchenko, S. B. R. Rempe, and S. V. Varma, Proc. Natl. Acad. Sci. U. S. A. 110, 12978 (2013). 30J. J. Molina, S. Lectez, S. Tazi, M. Salanne, J.-F. Dufrêche, J. Roques, E. Simoni, P. A. Madden, and P. Turq, J. Chem. Phys. 134, 014511 (2011). 31J. Timko, D. Bucher, and S. Kuyucak, J. Chem. Phys. 132, 114510 (2010). 32R. A. DiStasio, B. Santra, Z. Li, X. Wu, and R. Car, J. Chem. Phys. 141, 084502 (2014). 33A. P. Gaiduk, C. Zhang, F. Gygi, and G. Galli, Chem. Phys. Lett. 604, 89 (2014).
18F.
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
