Models of Ion Solvation Thermodynamics in ... - ACS Publications

Aug 21, 2015 - ... in water, and the other based on high-level quantum calculations on ion–solvent dimers and fitting to a Buckingham-type potential...
0 downloads 0 Views 977KB Size
Article pubs.acs.org/JPCB

Models of Ion Solvation Thermodynamics in Ethylene Carbonate and Propylene Carbonate Ayse Arslanargin,† August Powers,‡ Thomas L. Beck,*,†,‡ and Steven W. Rick§ †

Department of Physics, University of Cincinnati, Cincinnati, Ohio 45221, United States Department of Chemistry, University of Cincinnati, Cincinnati, Ohio 45221, United States § Department of Chemistry, University of New Orleans, New Orleans, Louisiana 70148, United States ‡

ABSTRACT: Ethylene carbonate (EC) and propylene carbonate (PC) are organic solvents used extensively in energy storage applications such as lithium-ion batteries and supercapacitors. Using statistical mechanical theory and computer simulations, this paper compares and contrasts the thermodynamics of ion solvation in EC and PC with the behavior observed in water. The EC and PC solvents are modeled with the AMBER (GAFF) force field. Ion−solvent interactions are treated with two pointcharge models: one using an existing Lennard-Jones ion parameter set optimized for solvation in water, and the other based on high-level quantum calculations on ion− solvent dimers and fitting to a Buckingham-type potential form. The second model produces a coordination number for the Li+ ion in closer agreement with experiment. Neither model yields consistently accurate solvation thermodynamic quantities (free energies, enthalpies, and entropies), however. The simulations and thermodynamic analysis illustrate key physical aspects of the solvation; the studies also point to necessary modifications of these simple models. In particular, the calculations show that polarization and associated dispersion forces are important and that well-optimized polarizable or quantum models are likely required to accurately reproduce condensed-phase properties of ions in these technologically important solvents.



INTRODUCTION The global demand for renewable energy sources is rapidly increasing. Since many of these sources exhibit intermittency, improvements in energy storage are a crucial aspect of ongoing and future renewable energy developments. Lithium-ion batteries (LIBs) and supercapacitors are expected to play important roles in renewable energy generation and in hybrid and electric vehicles as electrochemical energy storage systems.1−4 Organic solvents such as ethylene carbonate (EC), propylene carbonate (PC), and dimethyl carbonate (DMC) are widely used in the liquid electrolytes of these storage systems.5 Several studies have been conducted on the electrodes of LIBs, whereas less effort has been directed at the investigation of the liquid electrolytes.6 The electrolyte structure affects the efficiency of the ion transport both in the bulk and near interfaces, however, and improvements to the solvent structure are essential for battery performance enhancements. Despite a number of studies of these electrolytes, the solvation structure of the Li+ ion in alkyl carbonate solutions is still debated,7 with a consensus building that the coordination number in EC and PC is in the range 4−5 (see below).8 EC and PC (see Figure 1), the subjects of this study, are polar aprotic solvents that are utilized in the electrolytes of LIBs and supercapacitors due to their capability to dissolve large amounts of inorganic salts. These solvent molecules are notable for their large dipole moments (4.61 D for EC and 4.81 D for PC), large dielectric constants (89.8 for EC and 64.9 for PC), and high boiling points (close to 515 K for both). EC is not a © XXXX American Chemical Society

Figure 1. Images of the EC (a) and PC (b) molecules.

liquid at room temperature, however; its melting point is 36.4 °C, whereas PC has a melting point of −49 °C. It is known that mixing them can enhance the battery performance and cyclability,9,10 but PC has a negative impact on the graphite electrode in commercial LIB systems.7 Neither EC nor PC form hydrogen bonds, and both are expected to solvate cations Special Issue: Bruce C. Garrett Festschrift Received: July 16, 2015 Revised: August 18, 2015

A

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

Article

The Journal of Physical Chemistry B

Symmetry-Adapted Perturbation Theory. Symmetryadapted perturbation theory (SAPT) provides methods for calculating noncovalent interaction energies, wherein the perturbative approach creates a natural decomposition into physically meaningful components.45 The Hamiltonian for the dimer is given as

more strongly than anions due to the high concentration of negative partial charge on the carbonyl oxygen and diffuse positive charge distribution on the alkyl end of the molecule. We observe direct evidence of this asymmetry below. A large number of experimental11−19 and theoretical20−30 studies have been performed to investigate the solvation of ions in these solvents. Most of these studies have analyzed the coordination structure of ions in the EC/PC solutions (and in related cosolvents31) and conductivity and viscosity measurements. There have been several thermodynamic studies of ion-EC/PC systems.11−13,32−34 Free energy and enthalpy evaluations are important for understanding a wide range of physical and chemical phenomena.35,36 These thermodynamic quantities also provide a challenging test of the accuracy of force field models. Various theoretical and numerical methods have been developed for evaluating free energy changes.36 Recently the local molecular field theory (LMFT) of Weeks and co-workers,37,38 developed to study the interplay of local and long-ranged forces in determining liquid structure, has been extended to free energy calculations by Beck.39 The theory allows for a step-by-step physical partitioning of the hydration free energy into cavity formation, dispersion, local electrostatic, and far-field electrostatic parts; this partitioning has proven helpful in several studies of ion hydration.39−42 In the present study, we use this method to calculate the ion solvation free energies in EC and PC. From their temperature dependence we obtain the enthalpies and entropies. Here we focus on the thermodynamic properties of potassium halide salts in order to compare with two recent extensive experimental studies of their solvation in EC and PC.11,12 We first present MP2-level quantum chemical results for the binding energy curves of single ions with EC and PC. In addition, the various energetic contributions to the binding were obtained using symmetry-adapted perturbation theory (SAPT). The results give insights into the relative roles of exchange, electrostatics, induction, and dispersion for each ion−solvent pair interaction. The binding energy curves were fit to a modified Buckingham potential form for subsequent simulations and free energy studies. For the simulations, we first explored the basic structural properties of LiCl solutions in order to compare them with previous classical and quantum studies and experiments. The single-ion free energies of solvation were then obtained for the potassium salts using the LMFT approach discussed above. Significant qualitative differences were observed in relation to ion solvation in water. The results are compared with the experimental data in refs 11 and 12. The results indicate important roles for polarization and dispersion in the ion solvation. We discuss the solvation behavior in terms of specific ion effects as evidenced by a “volcano” plot similar to the data of Collins, Neilson, and Enderby for ions in water,43 and we close by exhibiting the asymmetry between cation and anion solvation through electrostatic linear-response fluctuation analysis. The asymmetry is the opposite of that observed in water,44 with the cation more strongly solvated than the anion (for similar ion sizes).

H = FA + FB + V + WA + WB

(1)

where F is the Fock operator of a monomer, V is the interaction between the monomers, and W is the fluctuation potential of a monomer. Using this perturbation, the interaction energy is given as46 ∞

E int =





(nkl) (nkl) + Eexch ) ∑ ∑ ∑ (Epol

(2)

n=1 k=0 l=0

where n is the order in V, k is the order in WA, and l is the order in WB. The various levels of SAPT are denoted by the order of truncation in the monomer fluctuation potentials (i.e., k + l). The typically attractive polarization term is balanced by the repulsive exchange term that arises from the wave function antisymmetry.46 Common levels of SAPT include SAPT0 and SAPT2. The lowest level, SAPT0, can efficiently handle large systems;46 however, its performance relies on a fortuitous large error cancellation in the dispersion terms that is caused by using a small basis set.47 SAPT2 has an accuracy similar to MP2 (see below), and is useful for examining electrostatic-dominated systems.46 Due to its favorable performance, as well as its availability in the Psi448 and SAPT201249 packages, SAPT2 has become a widely used tool for analyzing noncovalent interactions. The components of the interaction energy are grouped into electrostatic (elst), exchange (exch), induction (ind, including some degree of charge transfer), and dispersion (disp) terms. More detailed theoretical descriptions of the SAPT method are given in refs 45−47, including discussion of the relation of the expansion of eq 2 to the physical terms listed above. Solvation Free Energy. The solvation free energy of an ion interacting with the surrounding medium can be expressed with the potential distribution theorem (PDT)50 βμex = −ln⟨exp( −β ΔU )⟩0

(3)

where β = 1/kT and ΔU is the interaction energy of the ion with the solvent molecules. The zero subscript indicates that the solute and the solvent molecules are sampled independently. The inverse form of the PDT can be written as in ref 39 βμex = ln⟨exp(β ΔU )⟩

(4)

in which the solute is fully coupled to the solvent molecules in the sampling. Equation 4 can be rewritten by inclusion of an arbitrary model potential M: βμex = βμMex + ln⟨exp(β(ΔU − M )⟩

(5)

For example, the model potential may be chosen as the uncharged ion van der Waals (vdW) potential or a purely repulsive (but continuous) cavity potential. The LMFT partitioning approach divides the electrostatic interactions between the solute and the solvent into local and far-field parts as in the Ewald prescription.39 The solvation free energy including electrostatic partitioning is then expressed as



THEORY In this section we briefly discuss the SAPT approach for analyzing molecular interactions and the theory behind our free energy calculations. B

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

Article

The Journal of Physical Chemistry B

mined from the equilibrium points for the binding energy curves. Molecular Dynamics Simulations. Simulation Details. All of the simulations were performed with the GROMACS molecular dynamics package.58 The simulations were conducted for 215 solvent molecules (EC or PC) and 1 solute (ion) over a range of temperatures (313−453 K for EC and 298−438 K for PC in 20 K intervals). The solutes in these systems are the Li+, K+, F−, Cl−, and Br− ions. The temperature was maintained constant using the V-rescale thermostat. A time step of 1 fs was used. Each system was equilibrated for at least 500 ps before production runs of 2 ns or longer. All systems were modeled in a cubic box with periodic boundary conditions. The particle-mesh Ewald method was used for the electrostatic interactions with a real-space damping length η−1 of 4 Å for cations and 5.5 Å for anions. These length scales have been varied to optimize the overlap of the resulting interaction energy distributions as in ref 39. The Generalized AMBER force field (GAFF)59 is employed to describe bonded and nonbonded interactions excluding the ion−solvent vdW interactions. The GROMACS compatible topologies were generated using the ACPYPE script.60 Intermolecular interactions between the solvent molecules are described with the Lennard-Jones (LJ) potential and electrostatic interactions generated by partial charges. Intramolecular interactions are included for atoms joined by three covalent bonds (1−4 interactions). Solvent and Ion−Solvent Interactions. The AMBER (GAFF) force field has been employed previously in simulation studies of the EC and PC solvents themselves.22−24,29,61 It has been successful in reproducing several key properties of EC and PC such as the density, dielectric constant, and surface tension.24,61 For example, in ref 24 the calculations of the thermal expansion coefficient and the dielectric constant for liquid PC agree well with the values from experiment,62,63 and in ref 61, several other thermodynamic properties of EC were well-reproduced. It has been shown that the PC/graphite contact angle is consistent with experiment after a reduction on the graphite-C atom and PC LJ interactions.24 An exception is the enthalpy of vaporization of the two solvents for which the AMBER (GAFF) force field predicts values over 20% too large in magnitude (see ref 61; the PC value was obtained from our own calculations). We note that a 20% reduction in the partial charges on EC and PC (resulting in a corresponding reduction in the GAFF dipole values of 6.2 D for EC and 6.4 D for PC to values of 4.95 and 5.20 D, not much larger than the experimental gas phase values of 4.61 and 4.81 D) leads to vaporization enthalpies of 59.5 kJ/mol (EC) and 62.6 kJ/mol (PC), close to the experimental values of 61.0 and 61.5 kJ/mol. These results suggest the effective dipole moments on the EC and PC molecules may be too large in the GAFF point-charge model. Alteration of the charges will likely degrade the good agreement with some other thermodynamic quantities, however. In addition, due to the large mean polarizabilities of EC (6.8 Å3) and PC (8.7 Å3) compared with water (1.47 Å3), complex induction and dispersion effects can be expected11 that are difficult to handle with these simple point-charge models (this is also true for the ion−solvent interactions). Our initial goal was to use the AMBER (GAFF) without modification of the ion−solvent vdW interactions, and test the ability of the force field by calculation of the ion−solvent

ex βμex = βμvdW + ln⟨exp(β ΔUloc,es)⟩ΔUloc

+ ln⟨exp(β ΔUfar,es)⟩

(6)

where the first term is the vdW contribution, the second term is the local electrostatic part of the free energy, and the last term is the far-field contribution. The ΔUloc sampling in the second term involves the solute locally coupled (electrostatically) to the solvent and the vdW component. The vdW part of the free energy can be further partitioned into repulsive (cavity formation) and attractive dispersion parts. The last term involves electrostatic interactions at long range (outside the first shell). Further discussion of the partitioning approach and the details of its implementation can be found in ref 39. The real single-ion solvation free energy change for moving an ion from vacuum into the bulk includes a contribution from any potential shift across the water interface and can be written in two formats.42 In periodic boundary conditions, the real free energy is ex μex = μint + qϕsp

(7)

μex int

where is the intrinsic free energy (in the absence of a bulk interface), while ϕsp is the surface potential across the liquid/ vapor interface. Alternatively, it can be defined as μex = μ bex + qϕnp

(8)

μex b

where is the bulk ion solvation free energy and ϕnp is the net potential at the center of a neutral ion. The net potential is the sum of a local contribution due to an inhomogeneous solvent distribution near an uncharged vdW particle and the distant surface potential: ϕnp = ϕlp + ϕsp. The bulk value μex b can be obtained when the local potential contribution at the center of an uncharged solute (due to nearby solvent molecules) is removed from μex int (the free energies computed in periodic ex boundaries), μex b = μint − qϕlp. All surface potential effects cancel out when considering ion pair quantities.



COMPUTATIONAL DETAILS Quantum Chemistry Calculations. Binding Energy Curves. Quantum chemical binding energy data were calculated at the MP2/aug-cc-phwCVDZ level, using the Gaussian 09 package.51 The aug-cc-phwCVDZ basis set is shown in Table 1. Table 1. Composition of the aug-cc-phwCVDZ Basis Set atom

basis set

H Li C, N, O, F, Cl, Br K

aug-cc-pVDZ52 cc-pCVDZ53 aug-cc-pwCVDZ54 Feller’s CVDZ55,56

A rigid (unrelaxed) potential energy scan was performed along the ion−carbonyl oxygen vector for each ion−solvent dimer structure, starting from a structure optimized at the same level of theory. The ion−oxygen separation was adjusted from req − 0.4 Å to req + 4.0 Å, in increments of 0.1 Å (similar to the procedure in ref 57); at each point, a counterpoise-corrected binding energy was calculated. SAPT Interactions. Interaction energies for the ion−solvent dimers were calculated at the SAPT2/aug-cc-phwCVDZ level (using the Psi4 package) so as to provide an energetic decomposition as discussed above. The structures used were the same MP2/aug-cc-phwCVDZ optimized structures deterC

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

Article

The Journal of Physical Chemistry B

Figure 2. Li+−EC (a) and Cl−−EC (b) binding energy curves with Lennard-Jones and modified Buckingham potentials, and MP2 calculations.

binding energy curves, ion coordination numbers, and some thermodynamic properties. The AMBER (GAFF) force field does not have its own set of ion parameters, however; the AMBER ion parameter set can be included, but in a comparison to computed ion−solvent binding energy curves, these parameter sets gave poorer agreement with MP2-level calculations than the LJ model chosen below (which is already in significant disagreement with the quantum chemical results). Due to the relatively consistent success in modeling ions in water observed in ref 64, the initial ion LJ parameters were taken from that work. Using this model, our calculations for the coordination numbers of the Li+ ion in EC and PC showed a deviation from the experimental values, and the free energies of potassium salts in EC and PC were not in agreement with the experiments (see below for further details and results). The binding energies for ion−solvent pairs were of significantly smaller magnitude than the binding energies computed at the quantum level (see Figure 2). The LJ potential has been commonly used in classical simulations because of its computational efficiency. For carbonate-based solvents, however, it has been shown that the repulsive part of the vdW interactions is better described with the exp-6 form rather than (1/r12).25,65 As a first step toward improving the AMBER (GAFF) force field for the ion− solvent interactions, we have thus employed a modified exp-6 Buckingham potential65 for the vdW interactions between the ion and the solvent molecules UvdW = Aij exp( −Bij rij) −

Cij r6

⎛ 12 ⎞12 ⎟⎟ + D⎜⎜ ⎝ Bij rij ⎠

A = 6ε(exp λ)/(λ − 6) B = λ/r* C = ελ(r *)6 /(λ − 6)

where ε is the potential energy at the interatomic separation at the minimum r* and λ is the steepness parameter. In order to decrease the number of parameters, we have specified the λ value as 12.74 (it is usually taken within the interval 12−1666). The σ and ϵ parameters of the solvent atoms (for the ion− solvent interactions) were taken from the GAFF set, and used to reproduce the solvent A, B, and C parameters with the specified λ value (as shown above). The ion vdW parameters (Table 2) were obtained by fitting the ion−solvent dimer Table 2. Buckingham vdW Parameters for the EC/PC−Ion Interactionsa

(9)

Bii3Bjj3

CiiCjj

B (Å−1)

C (kcal/mol Å−6)

C3 OS C O H1 HC Li+ K+ F−(EC) F−(PC) Cl− Br−

33 161.6 51 535.6 26 068.5 63 659.9 4760.2 4760.2 20 651.2 39 809.9 35 407.5 30 721.7 45 755.2 40 548.1

3.750 4.250 3.750 4.303 5.160 4.808 6.703 3.475 4.100 4.510 3.039 2.773

319.5 234.3 251.1 267.0 6.8 10.3 6.1 601.7 198.3 96.8 1544.7 2625.4

binding energies from the molecular mechanics calculations to the binding energies calculated at the MP2/aug-cc-pwCVDZ level as previously discussed. As seen from Table 2 we have obtained the same ion−solvent parameters for both EC and PC systems excluding the F− ion. A Monte Carlo simulated annealing (MCSA) method was used to obtain the nonlinear fits to the dimer binding energies from quantum chemistry calculations. In the MCSA approach, the parameters were initially set to zero and then randomly sampled. To obtain a good fit, the root-mean square error (rmse) between the quantum

Bij6

⎛ ⎞1/6 2 ⎜ ⎟ Bij = ⎜ −6 −6 ⎟ ⎝ Bij + Bij ⎠ Cij =

A (kcal/mol)

The EC/PC values were derived by comparison with the AMBER (GAFF) force field. The ion parameters were fit to the MP2-level quantum chemical energy curves using the Monte Carlo Simulated Annealing algorithm discussed in the text.

The term D(12/Bijrij) , with D = 0.0005 kcal/mol for all pairs, is essentially zero for distances larger than typical interatomic separations, but becomes the dominant repulsive term at rij < 1 Å, where the Buckingham form becomes unphysically negative. The Aij, Bij, and Cij parameters are given by Waldman−Hagler (WH) combining rules Aii A jj

atom type

a

12

Aij =

(11)

(10)

The A, B, and C parameters can be expressed as D

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

Article

The Journal of Physical Chemistry B mechanics and molecular dynamics binding energies was minimized by variation of the parameters: f=

To provide a further test of our results, we have calculated the solvation enthalpies from the differences in potential energies between the solution and pure solvent systems:

∑i wi(Ei(fit) − Ei)2 ∑i wi

hex = ⟨UN + S⟩N + S − ⟨UN⟩N (12)

= ⟨ΔU ⟩ + ⟨UN⟩N + S − ⟨UN⟩N

Here i labels the points from the quantum data, Ei is the reference binding energy, and Ei(fit) is the resulting classical binding energy for the fitted parameter set. The weight wi is a Boltzmann weighting function (with the MP2-level energies in the exponent, thus optimizing the fit near the potential minimum) that has been applied to all points for more accurate fitting. The selection of the parameters is based on the Metropolis algorithm. In this algorithm, as the initial temperature T is decreased (annealing), random changes in the energy are accepted with a Boltzmann probability condition. The fitted ion−solvent vdW parameters are listed in Table 2. Figure 2 shows the comparison of the dimer binding energies for Li+−EC (Figure 2a) and Cl−−EC (Figure 2b) obtained from force fields and MP2 quantum chemistry calculations. The modified exp-6 Buckingham potential with fitted repulsion− dispersion parameters yields relatively good agreement with the binding energies computed at the MP2 level. We do note the nontrivial deviation of the fitted potential from the quantum result at larger distances for anions. The fitting for those cases could be improved by adding terms to the potential (such as a charge-induced dipole term); it is clear from our free energy results below, however, that the free energies and enthalpies are already somewhat too negative with models obtained as discussed above, so we do not pursue these issues further here. vdW Free Energies. The vdW contribution to the free energy includes two parts: cavity formation and the attractive dispersion term. The cavity formation free energy was calculated employing a repulsive Fermi-type potential40 λb U (r , λ ) = 1 + exp[a(r − λrc)]

= ⟨ΔU ⟩ + USR

Here USR is “the solvent reorganization energy” due to the addition of ion S. Using eq 15, the entropy is then s ex = − k ln⟨exp[β(ΔU − ⟨ΔU ⟩)]⟩ +

∂βμex ∂β

N

U=

N

N

∑ ∑ qiqjφEW (rij) + ∑ i=1 j=1

(17)

i=1

qi2ξ 2L

(18)

where rij = rj − ri + n is a vector in the unit cell with the realspace lattice vector n and L is the lattice size. The final term is the self-energy correction with ξ = −2.837 297. The position dependent electrostatic potential φEW(r) is defined using Ewald’s method35

(13)

erfc(η|r + nL|) + φEW (r) = ∑ |r + nL | n π − 3 2 Lη

2

∑ k≠0

2

4π e−k /4η i k . r e L3 k 2 (19)



where k = L n is the Fourier-space lattice vector. In ref 44 the electrostatic potential, ϕ, at the ion site (r = 0) is defined as N

ϕ=

M

∑ ∑ qi φEW (ri ) i=1 α=1

α

α

(20)

where the double sum extends over all atom sites in EC or PC. In this study ϕ was calculated for both ion−EC and ion−PC systems. The finite size corrected mean mc and the fluctuation fc of the potential at the ion site are described as44

(14)

Plots of β μ versus β display linear behavior over a relatively large temperature range, and the slope produces the enthalpy. Once the solvation free energy and enthalpy changes are known, the entropy of solvation can be computed directly: ex

hex − μex s = T

USR T

which is the sum of an energetic fluctuation term (the first nonzero cumulant is at second order) and a purely mechanical term that is typically positive and thus counters the negative fluctuation term. It is worth noting from these exact statistical mechanical formulas that the solvent reorganization energy is physically present in both the enthalpy and entropy, but cancels when constructing the free energy. Thus, reliable evaluation of the free energy requires accurate physical sampling of the system and accurate solute−solvent interaction energies, but the computed value does not depend directly on solvent− solvent interaction energies. It is interesting that, at least for the case of ions in water, the near cancellation of the fluctuation and solvent reorganization terms leads to small magnitude entropies, often of the same magnitude as those for rare gases.41 Electrostatic Potential. The Coulomb energy of a system with charges qj at positions rj (j = 1, ..., N) in a cubic volume V = L3 can be expressed as

and applying thermodynamic integration (TI) for the transition λ = 0 → 1. As the λ parameter increases the potential grows in a repulsive cavity-forming particle. The rc parameter was chosen as a value that yields the Fermi potential slightly inside the repulsive vdW potential wall. The a and b parameters were chosen as 10 Å−1 and 20 kcal/mol, respectively. A final perturbation step was employed to compute the free energy change to go from U(r, λ = 1) to the vdW repulsive wall. The electrostatic portions of the free energy are computed using the same methods discussed in ref 39. Solvation Enthalpies and Entropies. The decomposition of the solvation free energy into entropic and enthalpic contributions can provide additional insights into the solvation process. We have calculated the enthalpy of solvation using the van’t Hoff equation as follows

hex =

(16)

mc = e⟨ϕ⟩ + qeξ fc = βe 2⟨(ϕ − ⟨ϕ⟩)2 ⟩ − e 2ξ

ex

(21)

and calculated from molecular dynamics simulation trajectories for a range of charges from −e to +e in 0.25e increments.

(15) E

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

Article

The Journal of Physical Chemistry B



where R is the minimum after the first peak in the rdf gij(r). We have obtained a coordination number of 4 EC or PC molecules in the first shell around the Li+ ion using the modified Buckingham potential, while it is 6 for the Lennard-Jones potential. As seen from Figure 3a, using a fitted exponential repulsive term in the vdW potential yields a shift to shorter distances in the rdf of Li+ which appears to provide a more realistic representation of the repulsive wall compared with the quantum chemical data. The closer solvent approach, and strong binding to the Li+ ion, reduces the available space for additional solvent molecule insertion into the inner shell. This reduces the coordination number from 6 to 4. We note that ab initio quantum simulations27,28 have indicated that the first peak in the rdf is centered very close to 2 Å, which is situated between the curves for the Lennard-Jones and Buckingham potentials shown in Figure 3a. There have been some controversies on the coordination number of Li+ ion in EC or PC. Our result for a fourcoordinated Li+ ion in the solution (approximately tetrahedrally coordinated, see below) is in agreement with Raman intensity experiments,14,17 several ab initio studies,26−28 and classical molecular dynamics calculations.20,21,30 However, Hyodo and Okabayshi16 have reported the solvation number of Li+ ions in the first solvation shell to be 4−4.9 depending on the concentration of the solution. They also claimed that if there is only one solvation species, the solvation number of Li+ ions is 4 in LiClO4/EC solution. Kameda et al. have found a coordination number of 4.5 in the neutron diffraction studies of LiPF6 in PC,19 the same value as inferred from X-ray absorption spectroscopy and simulation in ref 7. The first-shell solvent occupation number for Cl− is about 19 which is in agreement with previous force field predictions.20,21 In our fitted potential, the interactions between the Cl− ion and the solvent molecules have also been enhanced relative to the Lennard-Jones binding curves, but we do not see a major change in the solvent structure around the anion (which shows a relatively weak solvation in both EC and PC). We also note that the computed coordination numbers for the K+ ion are 7 for the Lennard-Jones model and 6.25 for the fitted Buckingham model (and the rdfs show a trend similar to that in Figure 3a). Thermodynamics of Ion Solvation in EC and PC. The computed results for the single-ion free energies and enthalpies are presented in Table 4. We have included the free energies and enthalpies obtained by employing both LJ and Buckingham potentials. As discussed above, the ion free energies have been

RESULTS AND DISCUSSION SAPT Interactions. Energetic decompositions of the ion− solvent dimer interactions from the SAPT2 method are presented in Table 3. The decompositions are only for the Table 3. SAPT2/aug-cc-phwCVDZ Interaction Energies for Ion−Solvent Dimersa

a

ion

Eelst

Eexch

Li+ K+ F− Cl− Br−

−41.63 −28.04 −40.83 −28.23 −26.43

13.55 10.08 30.97 19.02 18.05

Li+ K+ F− Cl− Br−

−43.50 −29.01 −39.81 −26.90 −25.08

13.81 10.41 31.70 19.27 18.31

Eind EC Solvent −21.52 −8.43 −17.05 −8.21 −7.01 PC Solvent −21.09 −8.78 −19.03 −9.20 −7.85

Edisp

ESAPT

EMP2

−0.57 −1.73 −5.25 −5.43 −5.57

−50.17 −28.12 −32.17 −22.85 −20.97

−49.57 −27.42 −30.62 −22.12 −20.24

−0.50 −1.77 −5.79 −6.02 −6.21

−51.27 −29.14 −32.93 −22.85 −20.84

−51.20 −28.45 −30.84 −21.82 −19.88

All energies are given in kcal/mol.

ion−solvent pairs, so no many-body effects due to the collective behavior in the solvent are present. Still, some interesting points can be noted in terms of the ion specificity. First, each of the electrostatics, exchange, and induction terms display significant ion specificity. On the other hand, while the cation and anion dispersion terms are quite different (larger for anions as expected), the variation between the anions is not large. The magnitude of the exchange term for the fluoride ion is anomalously large. It is clear that induction is an important contributor to the interaction energy for both cations and anions, while dispersion is more important for anions. The total interaction energies constructed from the various parts of the SAPT2 energy eq 2 are quite close to the total MP2 binding energies. Coordination Numbers. The radial distribution functions (rdfs) of the Li+ and Cl− ions with the carbonyl oxygens of EC are displayed in Figure 3. The rdfs with PC solvent molecules show similar behaviors. The coordination numbers of the solvent molecules surrounding the solvation shell of the ions are computed from the following integral nij(R ) = 4πρj

∫0

R

gij(r )r 2 dr

(22)

Figure 3. (a) Li+−carbonyl oxygen rdfs and coordination numbers with Lennard-Jones and modified Buckingham potentials for ion-EC vdW interactions. (b) Cl−−carbonyl oxygen rdf and coordination number with modified Lennard-Jones and Buckingham potentials. F

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

Article

The Journal of Physical Chemistry B

is not surprising given the deeper well (especially for anions) in the model based on the quantum calculations (Figure 2). The ion enthalpies listed in Table 4 were computed using van’t Hoff analysis. Alternatively, we have calculated (for the Buckingham potential model) the solvation enthalpies using eq 16. The calculated ion enthalpies in EC are −102.5, −90.5, −74.6, −71.0 kcal/mol for the K+, F−, Cl−, and Br− ions, respectively; in PC the enthalpies are −102.3, −90.1, −72.0, −67.0 kcal/mol for the K+, F−, Cl−, and Br− ions, respectively. This approach provides qualitatively consistent results (up to 2−4 kcal/mol difference) with the values obtained by using van’t Hoff analysis. The relatively good agreement between the two different approaches is encouraging. In the remainder of the paper we use the van’t Hoff analysis results since these are expected to display less statistical uncertainty than the results from eq 16. For the fitted Buckingham model we display the computed entropies, which are consistently negative and of substantial magnitude. As we will see below, these computed results deviate seriously from experiment for solvation in EC. Since the average solute−solvent binding energies are available, eq 16 also produces the solvent reorganization energies (shown in Table 4). The strong ion−solvent interactions can lead to decreased solvent−solvent attractive interactions along with increased repulsive interactions, so the USR values are positive. The variations of the magnitudes of USR also roughly track the variations of the solvation free energies and enthalpies. To our knowledge there are no direct experimental studies of the single-ion free energies and enthalpies of transfer from water to EC. Marcus69 has reported the single-ion free energies and enthalpies of the transfer from water to PC. We have extracted the single-ion free energies and enthalpies using the transfer energies and the values for ion hydration free energies and enthalpies in bulk water from ref 70 with the net potential contribution of ϕnp = ±11.6 kcal/mol.42 However, the resulting single-ion free energies are more positive for K+ (−73.2 kcal/ mol), while they are more negative for halide ions (−102.6, −76.7, −72.6 kcal/mol for F−, Cl−, and Br−, respectively) than in our calculations. In a later study,71 Marcus discussed the possible presence of water impurities found in other solvents,

Table 4. Single-Ion Free Energies and Enthalpies Calculated for the Two Models at T = 313 K for EC and PCa ion

μex LJ

hex LJ

K+ F− Cl− Br−

−84.0 −67.0 −57.8 −50.5

−93.7 −77.7 −71.2 −65.8

K+ F− Cl− Br−

−87.1 −64.0 −51.3 −47.1

−100.9 −77.1 −65.5 −63.6

μex Buck EC Solvent −87.8 −82.4 −68.2 −61.1 PC Solvent −87.9 −83.5 −63.7 −59.5

hex Buck

sex Buck

USR

−97.6 −90.3 −75.4 −70.3

−31.3 −25.2 −23.0 −29.4

60.8 45.2 41.9 38.8

−99.1 −91.9 −71.6 −70.9

−35.8 −26.8 −25.2 −36.4

58.0 45.1 41.0 39.6

a

All free energies and enthalpies are in kcal/mol. The last two columns show the entropies (cal/(mol K)) and solvent reorganization energies USR (kcal/mol) for the Buckingham model. The estimated one standard deviation errors for the free energies are 0.2 kcal/mol, and those for the enthalpies are 1.2 kcal/mol.

shifted by a ϕlp value computed at the center of an uncharged vdW particle (around ±3 kcal/mol for ions in EC and ±2 kcal/ mol for ions in PC, depending on the charge, and ϕlp is negative). The shift ϕlp due to simulating periodic boundaries is necessary for comparison with bulk solvation free energies as discussed in ref 42. It is clear from these data that the cation is more strongly solvated than the anions. The first indication of major differences from ion solvation in water can be seen for the K+/F− pair (which are of similar size); the F− anion is much more strongly hydrated than the K+ cation67 due to the strong hydrogen bonding between water and the anion. Preferential solvation of the cation in EC/PC is due to the direct interaction with the “harsh” carbonyl oxygen. The asymmetry in EC/PC is not as strong as in water, however. We note that, by analyzing the bulk quantities, shifts due to surface potential effects are eliminated; those shifts may have an impact on interfacial behavior of the ions but do not alter ion behavior in the bulk solution.68 It is also clear that there is typically a downward shift in the free energies and enthalpies going from the LJ model to the fitted Buckingham model. This

Table 5. Free Energy, Enthalpy, and Entropy Changes of Solvation Results for Potassium Salts Computed at T = 313 K for EC and PCa salt

μex LJ

hex LJ

sex LJ

μex Buck

KF

−151.8(19.8)

−172.1(−3.6)

−64.9(−74.8)

EC Solvent −170.4(1.2)

KCl

−143.2(5.2)

−166.3(−16.6)

−73.8(−69.6)

−156.4(−8.0)

KBr

−135.9(7.1)

−160.9(−11.0)

−79.9(−57.9)

−150.1(−7.1)

KF

−151.7

−178.3

−85.0

PC Solvent −171.4

KCl

−140.4(3.6)

−168.7(−3.7)

−90.4(−24.6)

−152.4(−8.4)

KBr

−136.0(4.6)

−166.3(−5.4)

−96.8(−31.9)

−148.0(−7.4)

hex Buck

sex Buck

μex exp

hex exp

−187.9(−19.4) −179.4(−10.9) −173.3(−23.6) −164.8(−15.1) −169.1(−19.2) −160.6(−10.7)

−55.9(−65.8) −28.8(−38.7) −54.0(−49.8) −26.8(−22.6) −60.7(−38.7) −33.5(−11.5)

−171.6

−168.5

9.9

−148.4

−149.7

−4.2

−143.0

−149.9

−22.0

−190.8 −182.3 −173.2(−8.2) −164.7(0.3) −171.1(−10.2) −162.6(−1.7)

−62.0 −34.8 −66.5(0.7) −39.3(26.5) −73.8(−8.9) −46.7(18.2)

−144.4

−165.0

−65.8

−140.6

−160.9

−64.9

sex exp

a

Free energies and enthalpies are in kcal/mol, while entropies are in cal/(mol K). The values in parentheses are the deviations from the experimental values. The second row for each salt includes an 8.5 kcal/mol shift in the enthalpy (for the Buckingham model) to provide a rough estimate of the change in the solvent reorganization energy with the inclusion of induced dipole-induced dipole interactions between the solvent molecules. G

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

Article

The Journal of Physical Chemistry B

some collective enhancement of the molecular dipoles occurs in the condensed phase, the AMBER (GAFF) charges may be increased beyond the actual values in order to mimic other thermodynamic quantities besides the vaporization enthalpies. On the other hand, the molecular polarizabilities of EC and PC are large, and thus the strong field of nearby ions can be expected to induce significant dipoles in both EC and PC, in addition to the existing bulk values. In our Buckingham model, the pair interaction curve is well-reproduced, and thus the impact of the induced dipole is effectively present in the ion− solvent pair interaction. The model does not include the resulting interactions between the solvent induced dipoles in the neighborhood of the ion, however (which is expected to give a positive energetic contribution). We have examined the relevant angular distributions to determine the degree of tetrahedral order in the solvation shell around the Li+ ion as an example, and have found that both EC and PC display a significant degree of tetrahedral order (indicated by narrow peaks around the tetrahedral angles, and less for PC). If the induced solvent dipoles are roughly pointing toward the ion, a net repulsive contribution to USR should result. Our results for ion pair solvation in EC fully support this notion; the computed entropies are far too negative compared with those from experiment. We made a rough estimate of this effect by placing inducible point dipoles on the carbonyl carbon with the experimental gas phase polarizabilities (the sites respond to the fixed-charge external fields during the simulation), and we estimated a shift of +8.5 kcal/mol in the enthalpy due to induced dipole−induced dipole interactions. With the dipole located at the carbonyl carbon site, the impact of the anions is small, but this is likely an artifact of our model, since the polarization is likely distributed over the whole molecule. While this correction does not completely offset the strong deviation of the entropies for solvation in EC (which may be remedied by the inclusion of additional polarizable sites), it does move the results in the correct direction. Likely, the inclusion of induced polarization for the solvent molecules will alter the molecular configurations near the ion too, which may lead to better agreement with the ab initio rdfs discussed above. It is interesting to note that when a similar correction is applied to the solvation enthalpies in PC, the predicted entropies are then too positive. It is clear from the experimental data11,12 that the solvation entropies in EC are much closer to zero than those in PC. This is likely due to the fact that liquid EC is significantly more ordered than liquid PC (due to the presence of the asymmetric methyl group on PC).11,12 For example, the major decrease in the melting point for PC relative to EC would appear to be due to this effect. Thus, the introduction of ions into EC involves both a strong destructuring of the solvent and an increased (but changed) order around the ions. It appears that these effects largely cancel. In PC, however, the liquid is already significantly disordered, so the ordering of the PC molecules around the ions leads to negative solvation entropies. The local solvation shell still may display slightly more disorder relative to EC, thus reducing the impact of induced dipoles on the computed enthalpies (through the solvent reorganization term). This discussion has focused on induction (induced dipole) interactions; a proper treatment of many-body dispersion may also be required for higher accuracy results.11 As a crude test of the impact of the existing permanent dipoles on EC and PC, we calculated the solvation free energy

which may cause unreliable transfer free energies. The experiments in refs 11 and 12 probe the solubilities of the salts directly into the purified EC and PC solvents and thus do not suffer from the inclusion of water around the ions in the organic phase. Thus, below we will compare our simulated data for ion pairs to the measured results without attempting to separate the quantities into single-ion values. The solvation free energies, enthalpies, and entropies of potassium salts in EC and PC are summarized in Table 5. Experimental values are obtained from ref 11 (EC data) and ref 12 (PC data). The experimental free energy of solvation values in EC and PC are not available in refs 11 and 12. We have used the following procedure to get those values for comparison. The standard Gibbs free energy of solvation, ΔGsolv, can be calculated through the cycle in Figure 4 as ΔGsolv = −(ΔGsoln + ΔG latt)

(23)

Figure 4. Thermodynamic cycle shows the relationship between the free energies of solution and solvation and the lattice free energy.

where ΔGsoln is the Gibbs free energy of solution and ΔGlatt is the lattice free energy of the solid. Experimental ΔGsoln values were obtained from refs 11 and 12, and a critical selection of ΔGlatt of salts has been given in refs 13 and 72. As discussed above, the LJ model does not reproduce the local solvation structure observed in experiments. The computed solvation free energies are generally somewhat too positive, while the enthalpies are too negative, producing large errors in the negative direction for the solvation entropies (for both EC and PC). The fitted Buckingham model, on the other hand, yields a more physically realistic local solvation structure. The computed solvation free energies move in the correct direction, but overshoot the experimental values, while the errors in the enthalpies are even more negative than for the LJ model. The resulting entropies are far too negative for ion solvation in EC, but are of decent accuracy for solvation in PC. Thus, neither of the models provides a robust treatment of the ion solvation thermodynamics. We would suggest that a first default criterion for more accurate models of ion solvation is to reproduce the measured local solvation structure to decent accuracy. The fitted Buckingham model, based on quantum chemical calculations, is an improvement over the LJ model in this sense, although ab initio simulations27,28 indicate the actual first shell physical behavior in the rdf lies between the two models for the Li+ ion case. A clear focus for the development of future models involves an accurate treatment of polarization. The AMBER (GAFF) point-charge model, while it performs well for several thermodynamic properties of the pure solvents, may contain overenhanced effective molecular dipoles as evidenced by the overestimation of the vaporization enthalpies for both pure solvents. (Other point-charge models have produced more accurate vaporization enthalpies.73) While it is expected that H

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

Article

The Journal of Physical Chemistry B of the K+ ion in EC with the charges reduced by 20%. The ion− solvent pair potential was again fit to the quantum chemical data, and the free energies were recomputed. The free energy actually becomes more negative relative to the default EC parameter set (by 11 kcal/mol); this appears to occur due to an even more severe contraction of the inner solvation shell. Thus, there does not appear to be a simple fix in terms of altering the EC and PC force field charges. To close the discussion of ion solvation thermodynamics, we present results related to specific ion effects. Collins et al.43 have shown how cations and anions behave in water depending on their water affinity. In the so-called volcano plot, kosmotrope−kosmotrope (strongly hydrated) and chaotrope−chaotrope (weakly hydrated) ion pairs tend to stay together in water, while kosmotrope−chaotrope pairs dissociate extensively. Figure 5 shows the relationship between ΔHsoln and

water. This suggests that cations are more solvated in EC and PC, whereas anions are more hydrated in water. Despite these differences, the plot suggests that specific ion behavior can occur in non-hydrogen-bonding organic solvents as well as in water, as discussed in ref 11. We note that only one side of the volcano appears in our limited plot; more diverse behavior would be expected with a wider range of cations. Asymmetric Cation and Anion Solvation. Finally, we examine the cation/anion solvation asymmetry following the same analysis presented in ref 44. That analysis for ion solvation in water examined the electrostatic potential at the center of a sphere of charge q as a function of q, and also presented the potential fluctuations as a function of q. In linearresponse theory (equivalent to a Gaussian model), the slope of the potential versus q curve is proportional to the fluctuations. In addition, the solvation free energy is roughly half of the energy fluctuations divided by kT. Figure 6a shows the dependence of the average potential mc on the charge of a vdW particle with the size of the Cl− ion (the Buckingham parameters are given in Table 2). The plot displays indications of near-linear behavior, with some evident deviation (the magnitude of the slope increasing with more positive q). For q = 0, ⟨ϕ⟩ is −3.4 kcal/(mol e) for EC, and −1.9 kcal/(mol e) for PC. These values are simply the local potentials ϕlp of the solvent molecules for the Cl−-sized particle. We have also calculated the surface potential of both EC and PC from cluster simulations without any solute. The far-field part of the electrostatic potential at a chosen arbitrary point near the center of mass of the cluster gives the surface potential.42 The ϕsp value for EC is −0.8 kcal/(mol e), and for PC it is −3.9 kcal/(mol e), yielding net potential values of −4.2 kcal/(mol e) and −5.8 kcal/(mol e) for EC and PC, respectively. The value for the SPC/E water model is −4.3 kcal/(mol e), while a more robust estimate for water that does not rely on direct computation of the potential is −9.9 kcal/ (mol e).68 Figure 6b shows the fluctuations fc as a function of q/e. The larger fluctuations in the positive charge range imply that the structures with the carbonyl oxygen pointing toward the ion involve more intense electrostatic interactions. For negative charges, the ion is increasingly solvated by the opposite end of the molecule on which the positive partial charges are more distributed, and thus the fluctuations are smaller. This indicates that positive ions are more strongly solvated compared to

Figure 5. Relationship between the enthalpy of solution and the difference between the single solvation enthalpies of K+ and halide ions (F−, Cl−, and Br−) for both EC (○) and PC (△). ΔHsoln values are obtained from refs 11 and 12. ΔHsoln for KF is estimated from a thermodynamic cycle similar to that in Figure 4: ΔHsoln = −(ΔHsolv + ΔHlatt), where ΔHlatt is the lattice enthalpy used from ref 13, and ΔHsolv is our calculated value.

the difference in the solvation enthalpies of K+ and halide ions as in ref 43. This plot requires single-ion data on the x-axis; we use our computed enthalpy values for the differences. KF is the least soluble ion pair in EC and PC, while it is more soluble in water. KBr and KCl also show opposite behavior in contrast to

Figure 6. (a) Average electrostatic potential mc at the position of a chloride-sized vdW particle in EC and PC as a function of its charge. (b) Fluctuation of the electrostatic potential fc at the position of a chloride-sized vdW particle in EC and PC as a function of its charge. I

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

Article

The Journal of Physical Chemistry B

realistic solvation structure near the ions, but neither model produced accurate predictions of the solvation thermodynamic quantities in conjunction with the AMBER (GAFF) model for the EC and PC solvents. The results suggest that one important contributor to the deviations from experiment is the lack of realistic polarization in the present models. The EC and PC molecules have relatively large polarizabilities; in the strong fields of nearby ions, these molecules can thus be expected to exhibit large induced dipoles. Such polarization effects are expected to shift the computed enthalpies to more positive values (through the solvent reorganization energy term) and as a result also shift the entropies to more positive values. This effect is expected to be particularly important for solvation in EC, a more structured liquid than PC. In addition, it is likely that a better treatment of many-body dispersion effects would lead to more accurate predictions.11 Major progress has been made in developing more accurate polarizable models for organic solvents such as EC and PC25 and ions in them,65 but those models have not been tested against ion solvation thermodynamic measurements in dilute solutions. (The models have been successfully tested versus various experimental measurements in studies of ion solvation in acetonitrile75,76 and in high-concentration electrolytes relevant to lithium-ion batteries.77) A future goal should also be ab initio computation27,28 of the thermodynamic quantities to provide a more rigorous benchmark test of the approximate models. In addition, quantum models could yield insights into average effective dipole moments for the solvent molecules in the condensed phase (and other molecular properties) that may be helpful in developing more accurate classical models that could be used for simulations of large systems. The present work is an initial step aimed at motivating and laying the foundations for such future studies.

negative ions, which is the opposite of the behavior in water. This interpretation fits in with the discussion in ref 18. Figure 7 shows the average electrostatic potential at the center of an uncharged purely repulsive cavity as a function of

Figure 7. Electrostatic potential at the center of a vdW particle cavity in EC and PC as a function of the cavity radius. The dashed line shows the size of the particle (Cl− size).

radius in periodic boundaries (the cavity potential is taken as the WCA LJ-based form). The potential is negative, increases with increasing radius, and then decreases, indicating the solvent structure near the cavity evolves with increasing cavity size. This indicates that carbonyl oxygens of solvent molecules point toward the cavity to some extent, again supporting the notion that cation solvation is favored. A potential of the opposite sign is observed in classical models of water (such as SPC and SPC/E),74 where the water hydrogens point toward the cavity center to some degree.



CONCLUSIONS The solvation of ions in the organic solvents EC and PC has been examined with a combination of statistical mechanical theory and simulation. The goals were to gain a better understanding of solvation in these technologically important liquids, and to compare and contrast the results with ion solvation in water. The theoretical results were compared with recent experiments on alkali halide ions in EC and PC.11,12 It was observed that the solvation behavior in the organic solvents displays several physical differences from solvation in water; for example, there is an asymmetry of solvation that leads to stronger solvation of cations relative to anions, reversed from the behavior in water. There exist specific ion effects in the solvation, but they assume a different form than in water due to the different asymmetry.43 Ab initio quantum calculations that employ the SAPT partitioning of interaction energies showed ion specificity in all of the interactions, including exchange, electrostatics, induction, and dispersion. While the ion interactions with the EC and PC molecules include small dispersion contributions for the cations, the dispersion energies are larger for the anions. The different anions do not exhibit substantial ion specificity for these dispersion interactions, however, implying the other portions of the energy produce greater ion specificity than dispersion proper. Two classical ion−solvent force field models were employed in the simulations used to predict solvation thermodynamic quantities: an LJ model taken from recent studies of ions in water,64 and a Buckingham model fit to the high-level quantum chemical results. The second model led to a more physically



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +1 513 556 4886. Fax: +1 513 556 9239. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Lawrence Pratt, Pierandrea Lo Nostro, and Barry Ninham for many helpful discussions. We gratefully acknowledge financial support of this research by the National Science Foundation (T.L.B., CHE-1011746 and CHE-1266105; S.W.R., CHE-1301072), and a large grant of computing time at the Ohio Supercomputer Center where most of these calculations were performed. We dedicate this paper to honoring the career of Bruce Garrett, who has contributed heavily in many ways to the development of theoretical and computational chemistry and energy science.



REFERENCES

(1) Goodenough, J. B. Electrochemical Energy Storage in a Sustainable Modern Society. Energy Environ. Sci. 2014, 7, 14−18. (2) Etacheri, V.; Marom, R.; Elazari, R.; Salitra, G.; Aurbach, D. Challenges in the Development of Advanced Li-ion Batteries: A Review. Energy Environ. Sci. 2011, 4, 3243−3262. (3) Goodenough, J. B.; Park, K.-S. The Li-ion Rechargeable Battery: A Perspective. J. Am. Chem. Soc. 2013, 135, 1167−1176. J

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

Article

The Journal of Physical Chemistry B (4) Sharma, P.; Bhatti, T. A Review on Electrochemical Double-layer Capacitors. Energy Convers. Manage. 2010, 51, 2901−2912. (5) Winter, M.; Brodd, R. J. What are Batteries, Fuel Cells, and Supercapacitors? Chem. Rev. 2004, 104, 4245−4270. (6) Xu, K. Nonaqueous Liquid Electrolytes for Lithium-Based Rechargeable Batteries. Chem. Rev. 2004, 104, 4303−4418. (7) Smith, J. W.; Lam, R. K.; Sheardy, A. T.; Shih, O.; Rizzuto, A. M.; Borodin, O.; Harris, S. J.; Prendergast, D.; Saykally, R. J. X-Ray Absorption Spectroscopy of LiBF4 in Propylene Carbonate: A Model Lithium Ion Battery Electrolyte. Phys. Chem. Chem. Phys. 2014, 16, 23568−23575. (8) Allen, J. L.; Borodin, O.; Seo, D. M.; Henderson, W. A. Combined quantum chemical/Raman spectroscopic analyses of Li+ cation solvation: Cyclic carbonate solventsdEthylene carbonate and propylene carbonate. J. Power Sources 2014, 267, 821−830. (9) Aurbach, D.; Gofer, Y.; Ben-Zion, M.; Aped, P. The Behaviour of Lithium Electrodes in Propylene and Ethylene Carbonate: The Major Factors that Influence Li Cycling Efficiency. J. Electroanal. Chem. 1992, 339, 451−471. (10) Zhang, S.; Xu, K.; Allen, J.; Jow, T. Effect of Propylene Carbonate on the Low Temperature Performance of Li-ion Cells. J. Power Sources 2002, 110, 216−221. (11) Peruzzi, N.; Ninham, B. W.; Lo Nostro, P.; Baglioni, P. Hofmeister Phenomena in Nonaqueous Media: The Solubility of Electrolytes in Ethylene Carbonate. J. Phys. Chem. B 2012, 116, 14398−14405. (12) Peruzzi, N.; Lo Nostro, P.; Ninham, B.; Baglioni, P. The Solvation of Anions in Propylene Carbonate. J. Solution Chem. 2015, 44, 1224−1239. (13) Salomon, M. Thermodynamics of Ion Solvation in Water and Propylene Carbonate. J. Phys. Chem. 1970, 74, 2519−2524. (14) Morita, M.; Asai, Y.; Yoshimoto, N.; Ishikawa, M. A Raman Spectroscopic Study of Organic Electrolyte Solutions Based on Binary Solvent Systems of Ethylene Carbonate with Low Viscosity Solvents which Dissolve Different Lithium Salts. J. Chem. Soc., Faraday Trans. 1998, 94, 3451−3456. (15) Kameda, Y.; Umebayashi, Y.; Takeuchi, M.; Wahab, M. A.; Fukuda, S.; Ishiguro, S.-i.; Sasaki, M.; Amo, Y.; Usuki, T. Solvation Structure of Li+ in Concentrated LiPF6− Propylene Carbonate Solutions. J. Phys. Chem. B 2007, 111, 6104−6109. (16) Hyodo, S.-A.; Okabayashi, K. Raman Intensity Study of Local Structure in Non-aqueous Electrolyte Solutions-I. Cation-solvent Interaction in LiClO4/Ethylene Carbonate. Electrochim. Acta 1989, 34, 1551−1556. (17) Hyodo, S.-A.; Okabayashi, K. Raman Intensity Study of Local Structure in Non-aqueous Electrolyte Solutions-II. Cation-solvent Interaction in Mixed Solvent Systems and Selective Solvation. Electrochim. Acta 1989, 34, 1557−1561. (18) Wood, R.; Craft, Q. Freezing Points, Osmotic Coefficients, and Activity Coefficients of Some Salts in Ethylene Carbonate: A High Dielectric Constant Solvent without Hydrogen Bonding. J. Solution Chem. 1978, 7, 799−812. (19) Kameda, Y.; Umebayashi, Y.; Takeuchi, M.; Wahab, M. A.; Fukuda, S.; Ishiguro, S.-i.; Sasaki, M.; Amo, Y.; Usuki, T. Solvation Structure of Li+ in Concentrated LiPF6-Propylene Carbonate Solutions. J. Phys. Chem. B 2007, 111, 6104−6109. (20) Li, T.; Balbuena, P. B. Theoretical Studies of Lithium Perchlorate in Ethylene Carbonate, Propylene Carbonate, and Their Mixtures. J. Electrochem. Soc. 1999, 146, 3613−3622. (21) Soetens, J.-C.; Millot, C.; Maigret, B. Molecular Dynamics Simulation of Li+BF−4 in Ethylene Carbonate, Propylene Carbonate, and Dimethyl Carbonate Solvents. J. Phys. Chem. A 1998, 102, 1055− 1061. (22) Yang, L.; Fishbine, B. H.; Migliori, A.; Pratt, L. R. Molecular Simulation of Electric Double-Layer Capacitors Based on Carbon Nanotube Forests. J. Am. Chem. Soc. 2009, 131, 12373−12376. (23) Yang, L.; Fishbine, B. H.; Migliori, A.; Pratt, L. R. Dielectric Saturation of Liquid Propylene Carbonate in Electrical Energy Storage Applications. J. Chem. Phys. 2010, 132, 044701.

(24) You, X.; Chaudhari, M. I.; Pratt, L. R.; Pesika, N.; Aritakula, K. M.; Rick, S. W. Interfaces of Propylene Carbonate. J. Chem. Phys. 2013, 138, 114708. (25) Borodin, O.; Smith, G. D. Development of Many-Body Polarizable Force Fields for Li-Battery Components: 1. Ether, Alkane, and Carbonate-Based Solvents. J. Phys. Chem. B 2006, 110, 6279− 6292. (26) Leung, K.; Budzien, J. L. Ab initio Molecular Dynamics Simulations of the Initial Stages of Solid-electrolyte Interphase Formation on Lithium Ion Battery Graphitic Anodes. Phys. Chem. Chem. Phys. 2010, 12, 6583−6586. (27) Ganesh, P.; Jiang, D.-E.; Kent, P. R. C. Accurate Static and Dynamic Properties of Liquid Electrolytes for Li-Ion Batteries from Ab initio Molecular Dynamics. J. Phys. Chem. B 2011, 115, 3085−3090. (28) Bhatt, M. D.; Cho, M.; Cho, K. Density Functional Theory Calculations and Ab initio Molecular Dynamics Simulations for Diffusion of Li+ within Liquid Ethylene Carbonate. Modell. Simul. Mater. Sci. Eng. 2012, 20, 065004. (29) Hamad, I. A.; Novotny, M. A.; Wipf, D. O.; Rikvold, P. A. A New Battery-charging Method Suggested by Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2010, 12, 2740−2743. (30) Masia, M.; Probst, M.; Rey, R. Ethylene Carbonate-Li+: A Theoretical Study of Structural and Vibrational Properties in Gas and Liquid Phases. J. Phys. Chem. B 2004, 108, 2016−2027. (31) Afroz, T.; Seo, D. M.; Han, S.-D.; Boyle, P. D.; Henderson, W. A. Structural Interactions within Lithium Salt Solvates: Acyclic Carbonates and Esters. J. Phys. Chem. C 2015, 119, 7022−7027. (32) Muhuri, P. K.; Ghosh, S. K.; Hazra, D. K. Solubilities of Some Alkali-metal Salts, Tetraphenylarsonium Chloride, and Tetraphenylphosphonium Bromide in Propylene Carbonate at 25 °C Using the Ion-selective Electrode Technique. J. Chem. Eng. Data 1993, 38, 242− 244. (33) Jones, J.; Anouti, M.; Caillon-Caravanier, M.; Willmann, P.; Lemordant, D. Thermodynamic of LiF Dissolution in Alkylcarbonates and Some of Their Mixtures with Water. Fluid Phase Equilib. 2009, 285, 62−68. (34) Jones, J.; Anouti, M.; Caillon-Caravanier, M.; Willmann, P.; Lemordant, D. Lithium Fluoride Dissolution Equilibria in Cyclic Alkylcarbonates and Water. J. Mol. Liq. 2010, 153, 146−152. (35) Beck, T. L.; Paulaitis, M. E.; Pratt, L. R. The Potential Distribution Theorem and Models of Molecular Solutions;. Cambridge University Press: New York, 2006. (36) Chipot, C.; Pohorille, A. Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer-Verlag: Berlin, 2007. (37) Rodgers, J.; Weeks, J. Local Molecular Field Theory for the Treatment of Electrostatics. J. Phys.: Condens. Matter 2008, 20, 494206. (38) Rodgers, J.; Weeks, J. Interplay of Local Hydrogen-bonding and Long-ranged Dipolar Forces in Simulations of Confined Water. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 19136−19141. (39) Beck, T. L. Hydration Free Energies by Energetic Partitioning of the Potential Distribution Theorem. J. Stat. Phys. 2011, 145, 335−354. (40) Arslanargin, A.; Beck, T. L. Free Energy Partitioning Analysis of the Driving Forces that Determine Ion Density Profiles Near the Water Liquid-vapor Interface. J. Chem. Phys. 2012, 136, 104503. (41) Beck, T. L. A Local Entropic Signature of Specific Ion Hydration. J. Phys. Chem. B 2011, 115, 9776−9781. (42) Beck, T. L. The Influence of Water Interfacial Potentials on Ion Hydration in Bulk Water and Near Interfaces. Chem. Phys. Lett. 2013, 561−562, 1−13. (43) Collins, K. D.; Neilson, G. W.; Enderby, J. E. Ions in Water: Characterizing the Forces that Control Chemical Processes and Biological Structure. Biophys. Chem. 2007, 128, 95−104. (44) Hummer, G.; Pratt, L. R.; Garcìa, A. E. Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100, 1206−1215. (45) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (46) Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. Levels of Symmetry Adapted Perturbation Theory (SAPT). I. K

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

Article

The Journal of Physical Chemistry B Efficiency and Performance for Interaction Energies. J. Chem. Phys. 2014, 140, 094106. (47) Hohenstein, E. G.; Sherrill, C. D. Wavefunction Methods for Noncovalent Interactions. Wiley Interdisc. Rev.: Computat. Mol. Sci. 2012, 2, 304−326. (48) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; et al. Psi4: An Open-source Ab initio Electronic Structure Program. Wiley Interdisc. Revs: Comput. Mol. Sci. 2012, 2, 556−565. (49) Bukowski, R.; Cencek, W.; Jankowski, P.; Jeziorski, B.; Jeziorska, M.; Korona, T.; Kucharski, S. A.; Lotrich, V. F.; Misquitta, A. J.; Moszynski, R. et al. SAPT2012: An Ab Initio Program for Many-Body Symmetry-Adapted Perturbation Theory Calculations of Intermolecular Interaction Energies; 2015; http://physics.udel.edu/szalewic/SAPT. (50) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (51) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2013. (52) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (53) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. V. Core-valence Basis Sets for Boron through Neon. J. Chem. Phys. 1995, 103, 4572−4585. (54) Peterson, K. A.; Dunning, T. H. Accurate Correlation Consistent Basis Sets for Molecular Core-valence Correlation Effects: The Second Row Atoms Al-Ar, and the First Row Atoms B-Ne Revisited. J. Chem. Phys. 2002, 117, 10548−10560. (55) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (56) Feller, D.; Glendening, E. D.; Woon, D. E.; Feyereisen, M. W. An Extended Basis Set Ab initio Study of Alkali Metal Cation-water Clusters. J. Chem. Phys. 1995, 103, 3526−3542. (57) Allinger, N. L.; Durkin, K. A. Van der Waals Effects between Hydrogen and First-row Atoms in Molecular Mechanics (MM3/ MM4). J. Comput. Chem. 2000, 21, 1229−1242. (58) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (59) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General AMBER Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (60) da Silva, A. W. S.; Vranken, W. F. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. (61) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8, 61−74. (62) Piekarski, H.; Kubalczyk, K.; Wasiak, M. Volumes, Heat Capacities, and Compressibilities of the Mixtures of Acetonitrile with N,N-Dimethylacetamide and Propylene Carbonate. J. Chem. Eng. Data 2010, 55, 5435−5440. (63) Côté, J.-F.; Brouillette, D.; Desnoyers, J.; Rouleau, J.-F.; StArnaud, J.-M.; Perron, G. Dielectric Constants of Acetonitrile, γbutyrolactone, Propylene carbonate, and 1,2-dimethoxyethane as a Function of Pressure and Temperature. J. Solution Chem. 1996, 25, 1163−1173. (64) Horinek, D.; Herz, A.; Vrbka, L.; Sedlmeier, F.; Mamatkulov, S. I.; Netz, R. R. Specific Ion Adsorption at the Air/water Interface: The Role of Hydrophobic Solvation. Chem. Phys. Lett. 2009, 479, 173−183. (65) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463−11478.

(66) Aseyev, G. G. Electrolytes: Supramolecular Interactions and NonEquilibrium Phenomena in Concentrated Solutions; CRC Press: Boca Raton, FL, 2014. (67) Marcus, Y. Ion Solvation; John Wiley, New York, 1985. (68) Pollard, T.; Beck, T. L. The Thermodynamics of Proton Hydration and the Electrochemical Surface Potential of Water. J. Chem. Phys. 2014, 141, 18C512. (69) Marcus, Y.; Kamlet, M. J.; Taft, R. W. Linear Solvation Energy Relationships: Standard Molar Gibbs Free Energies and Enthalpies of Transfer of Ions from Water into Nonaqueous Solvents. J. Phys. Chem. 1988, 92, 3613−3622. (70) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R. The Proton’s Absolute Aqueous Enthalpy and Gibbs Free Energy of Solvation from Cluster-Ion Solvation Data. J. Phys. Chem. A 1998, 102, 7787−7794. (71) Marcus, Y. Some Thermodynamic Aspects of Ion Transfer. Electrochim. Acta 1998, 44, 91−98. (72) Criss, C. M.; Salomon, M. In Physical Chemistry of Organic Solvent Systems; Covington, A., Dickinson, T., Eds.; Springer: New York, 1973; pp 253−329. (73) Silva, L. B.; Freitas, L. Structural and Thermodynamic Properties of Liquid Ethylene Carbonate and Propylene Carbonate by Monte Carlo Simulations. J. Mol. Struct.: THEOCHEM 2007, 806, 23−34. (74) Ashbaugh, H. S. Convergence of Molecular and Macroscopic Continuum Descriptions of Ion Hydration. J. Phys. Chem. B 2000, 104, 7235−7238. (75) Seo, D. M.; Borodin, O.; Han, S.; Boyle, P. D.; Henderson, W. A. Electrolyte Solvation and Ionic Association II. Acetonitrile-Lithium Salt Mixtures: Highly Dissociated Salts. J. Electrochem. Soc. 2012, 159, A1489−A1500. (76) Seo, D. M.; Borodin, O.; Han, S.; Ly, Q.; Boyle, P. D.; Henderson, W. A. Electrolyte Solvation and Ionic Association I. Acetonitrile-Lithium Salt Mixtures: Intermediate and Highly Associated Salts. J. Electrochem. Soc. 2012, 159, A553−A565. (77) McOwen, D. W.; Seo, D. M.; Borodin, O.; Vatamanu, J.; Boyle, P. D.; Henderson, W. A. Concentrated electrolytes: decrypting electrolyte properties and reassessing Al corrosion mechanisms. Energy Environ. Sci. 2014, 7, 416−426.

L

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