Relative Complexation Energies for Li+ Ion in Solution: Molecular

Dec 23, 2009 - Relative complexation energies for the lithium cation in acetonitrile and diethyl ether have been studied. Quantum-chemical calculation...
0 downloads 11 Views 738KB Size
Download PDF
J. Phys. Chem. A 2010, 114, 973–979

973

Relative Complexation Energies for Li+ Ion in Solution: Molecular Level Solvation Versus Polarizable Continuum Model Study Andrzej Eilmes* and Piotr Kubisiak Faculty of Chemistry, Jagiellonian UniVersity, Ingardena 3, 30-060, Krako´w, Poland ReceiVed: July 31, 2009; ReVised Manuscript ReceiVed: NoVember 17, 2009

Relative complexation energies for the lithium cation in acetonitrile and diethyl ether have been studied. Quantum-chemical calculations explicitly describing the solvation of Li+ have been performed based on structures obtained from molecular dynamics simulations. The effect of an increasing number of solvent molecules beyond the first solvation shell has been found to consist in reduction of the differences in complexation energies for different coordination numbers. Explicit-solvation data have served as a benchmark to the results of polarizable continuum model (PCM) calculations. It has been demonstrated that the PCM approach can yield relative complexation energies comparable to the predictions based on molecular-level solvation, but at significantly lower computational cost. The best agreement between the explicit-solvation and the PCM results has been obtained when the van der Waals surface was adopted to build the molecular cavity. 1. Introduction Solvation of inorganic ions in organic solvents plays a great role in chemistry and biochemistry since many chemical reactions occur in solution. Modeling of the ion-solvent interactions is therefore of great importance. Quantum-chemical calculations can provide information about the structure of ion-solvent complexes and about their binding energies, helping to determine the relative stabilities of the complexes with different coordination of the ion. However, typical quantum-chemical calculations represent the situation in the vacuum and, strictly speaking, pertain to the complexes in the gas phase rather than in the liquid. In solution the properties of the complex are affected by the solvent molecules located outside the first solvation shell, and it is reasonable to expect that the presence of bulk solvent should lead to changes in the relative complex stabilities. Effects of this kind are highly relevant for ion-polymer interactions in polymer-based solid electrolytes. In this case the polymer matrix surrounding the ion acts as the solvent, influencing the interactions in the complex. In quantum-chemical studies of such systems relatively small molecules are used to model the polymer, for example, in the investigations of ion-polymer interactions in electrolytes based on poly(ethylene oxide), aliphatic esters or oligoglyme molecules commonly serve this purpose.1-3 Similar calculations were reported for gel electrolytes based on poly(methyl metacrylate);4 ab initio calculations were also employed to study the solvation structures of Li+ in nonaqueous solvents used in battery electrolytes.5,6 As recognized recently, solvent effects play an important role in liquid and polymer electrolytes used for lithium batteries, modifying the strength of cation-anion interactions7,8 and the relative complexation energies of the cation9 or inducing the shifts of selected vibrational frequencies of the polymer.9,10 Information on distribution functions, solvation shell structure, relative abundance of structures characterized by different coordination numbers (CN) or on ion movements is commonly * To whom correspondence should be addressed. E-mail address: eilmes@ chemia.uj.edu.pl, Fax: +48 12 6340515.

obtained from molecular dynamics (MD) simulations. Vast literature is available on this subject, a series of papers on poly(ethylene oxide)11-13 being a prime example. On the other hand, this approach is not very well suited to learn about the relative stability of particular ion-polymer complexes, which, although not directly available from experiment, is important to understanding the details of ion-solvent interactions. Such data are more easily obtainable from typical quantum-chemical calculations; the disadvantage is that in that case they refer to the gas phase, not to the polymer environment. The relatively inexpensive expedient to acquire information about the relative energies of ion-polymer complexes with different CN not in vacuum, but in a dielectric medium, is to apply some kind of effective-solvent approach. We have done it recently, presenting a polarizable continuum model (PCM) study9 of the ion-oligoglyme interactions in a model polymer electrolyte. The results confirmed the intuitive expectation that the energy differences between complexes with different coordination numbers are significantly suppressed in solution (polymer matrix) compared to their vacuum counterparts. However, the validity of the detailed conclusions is limited by the inherent shortcomings of the applied methodology. First, it is not obvious in any individual case to what extent a set of discrete solvent molecules may be approximated by a polarizable continuum; second, the results suffer from some arbitrariness because of many possible parametrizations of the continuum model. The objective of the present paper is to verify the validity of our previous conclusions against a more fundamental approach where the polymer chains are explicitly represented by diethyl ether molecules. Admittedly, comparison between complexation energies calculated within the explicit-solvent model and the PCM approach is relevant not merely to the situation in polymer electrolytes, but samples a more general problem of modeling ion interactions with organic solvents, addressed in numerous quantum-chemical studies on the interactions of mono- or divalent ions with acetonitrile,14 acetone, tetrahydrofuran, and diethyl ether,15 acetone and N-methylacetamide,16 methanol,17

10.1021/jp907359a  2010 American Chemical Society Published on Web 12/23/2009

974

J. Phys. Chem. A, Vol. 114, No. 2, 2010

short chain alcohols,18 and aliphatic19 or aromatic amines.20 However, our primary goal is to attain conclusions specifically for the solvent mimicking the polymer matrix in a solid electrolyte; as we have no claim to generality, systematic investigation of a large variety of possible organic solvents is beyond our present scope. Nevertheless, in order to check the degree of (expected) sensitivity of the results to the dielectric constant of the medium, we decided to extend our study to acetonitrile, viewed as a representative highly polar solvent. In our present study we adopt a strategy similar to other works combining explicit/implicit solvent approach.21-23 We apply MD to obtain sets of complexes differing in the ion coordination number (Section 2.1), and for these structures we subsequently calculate the complexation energies by quantum chemical methods (Section 2.2). These results serve as reference data. The observed trends are then compared to the results of PCM calculations with different parametrization of the molecular cavity (Section 2.3). 2. Calculations 2.1. Molecular Dynamics Simulations. For our study of the energetics of Li+ complexation, we choose two commonly used organic solvents: one with a high dielectric constantsacetonitrile (ACN), and the other much less polarsdiethyl ether (DEE). Diethyl ether resembles the model ether or oligoglyme molecules used to study the ion-polymer interactions in solid electrolytes based on poly(ethylene oxide); consequently, the conclusions drawn for the DEE case are presumably transferable to polymer solid electrolytes. Acetonitrile is a typical polar nonaqueous solvent, and its interactions with ions have been extensively modeled both quantum-chemically14,24-27 and by MD28-30 simulations. We include it in our studies to obtain some estimates of the effect of solvent polarity. As a first step we prepared a series of structures of the complexed Li+ ion. Tinker v. 4.231 was used to perform the MD simulations of the lithium ion solvated in ACN or DEE. For both solvents, parametrization of bonded interactions (bonds, angles, and dihedral angles) and Buckingham-type van der Waals interactions were taken from the MM3(2000) force field.32 The electrostatic part of the potential energy was restricted to charge-charge interactions in vacuum (ε ) 1). Partial charges on the acetonitrile molecule were obtained from the Merz-Kollman fit to the electrostatic potential calculated using the B3LYP functional and the aug-cc-pVDZ basis set. The charges resulting from this procedure were: -0.497e on N, and 0.448e and -0.468e on C(sp) and C(sp3), respectively, and 0.1723e on H atoms. These values are close to the partial charges used in much more advanced potentials applied to model microsolvation of Na+ and I- in acetonitrile clusters14 or to perform MD simulations for monovalent ions in ACN.30 For diethyl ether we originally tried charge parametrization derived for poly(ethylene oxide)-based electrolyte11 (-0.28e on O atom, with other charges slightly modified to ensure charge neutrality of the DEE molecule), since the main reason of choosing DEE in this work was to model a solvent resembling the polymer electrolyte matrix. However, preliminary simulations showed that the negative charge on the oxygen atom is apparently not large enough to significantly reorient the ether molecules surrounding the Li+ ion, which leads to weak coordination of the ion by oxygen atoms. Therefore, we increased the charge on the oxygen atom to -0.48e (close to the partial-charge fit to the electrostatic potential in B3LYP/ aug-cc-pVDZ calculations) and modified accordingly other charges: -0.03e on all C atoms and 0.06e on hydrogens.

Eilmes and Kubisiak

Figure 1. Average number of N or O atoms found within distance r to the Li+ ion.

Our model systems consisted of one Li+ ion and 1000 ACN or 600 DEE molecules. The simulations were performed with periodic boundary conditions in the NVT ensemble at 300 K; the density was fixed at 0.786 g/cm3 for ACN and 0.713 g/cm3 for DEE, which yields the size of the periodic box as 44.25 Å for ACN and 46.95 Å for DEE. The cutoff of 9 Å was used for the nonbonded interactions including the electrostatic ones. The time step of the simulation was 1 fs. After at least 100 ps of equilibration snapshots of the system were collected every 0.1 ps for about 200 ps. Figure 1 shows the average number of nitrogen or oxygen atoms found within a given distance to the lithium ion. The dependence for ACN agrees well with the results obtained in ref 30 from MD simulations with nonpolarizable force field. The Li+ coordination number (CN) reaches 5 for Li-N distance about 2.8 Å, then rises slowly and for distances up to 4.5 Å the CN is close to 6; finally it starts to increase rapidly above Li-N separation of 5 Å. In the case of DEE, the number of oxygen atoms coordinating the ion reaches 4 at about 3 Å and remains practically constant until it begins to increase above the Li-O distance of 5 Å. Smaller CN values for DEE result apparently from the larger size of DEE and steric repulsion between C2H5 groups. A CN value close to 4 obtained in our simulations is comparable to CN ) 4-5 measured for Li+ in a model poly(ethylene oxide) electrolyte.33 Atomistic structures of the Li+ complexes with ACN or DEE obtained from the MD simulations were used in the next step as input geometries in quantum-chemical calculations of the binding energies. 2.2. Quantum-Chemical Study of Li+ Solvation at the Molecular Level. In this section we want to investigate how the relative energies of Li+ complexes with different CN evolve when the complex is embedded in increasing number of solvent molecules. To this end, we need a definition of the Li+ coordination numberswhich is straightforward in vacuum calculations (where CN usually equals the number of solvent molecules used to build the complex), but is not so easily defined in the solvent where there is a distribution of solvent molecules near the ion. Nevertheless, one may try to count the coordinating atoms (N or O) within a fixed radius from the cation and use the resulting “apparent” CN to distinguish between the structures with more or less efficiently coordinated Li+ ion. For our purposes the radius used to calculate the Li+ CN should allow for broad range of possible CN values, therefore it should be smaller than the distance at which the plateaus visible in Figure 1 appear. For ACN we used 2.60 Å, which gives the range of CN from 2 to 6 (with some exceptional cases of CN ) 1), and the average CN about 4.24. A similar radius

Complexation Energies for Li+ Ion in Solution

J. Phys. Chem. A, Vol. 114, No. 2, 2010 975

of 2.65 Å for DEE results in CN varying between 2 and 4 with the mean value 3.69. Admittedly, this arbitrary choice leads to somewhat artificial CN values, as the real difference between two structures with Li+ 4 or 6 may consist, for instance, in a mere tiny shift of two ACN molecules; nevertheless, the differences in the CN should be noticeable when the averages over series of structures are compared. On the basis of this argument, we selected randomly from the MD results sets of structures corresponding to different values of the “effective” Li+ coordination number. For each CN we used 20 structures with ACN and 10 structures with DEE. Then for each structure we performed a series of energy calculations for Li+ solvated in CN + M solvent molecules; with M ) 0, 2, 5, 10, 15, and 20 for ACN and M ) 0, 5, 10, and 15 for DEE. The case M ) 0 corresponds to typical vacuum calculations; with increasing M the original complex is solvated by an increasing number of explicit solvent molecules. All quantum-chemical results reported in this paper were obtained using the Gaussian 0334 package. In the calculations of complexation energies we applied the DFT methodology with the B3LYP functional and the 6-31+G** basis set. There are many possible ways of defining the complexation energy ∆Ec for the complex of Li+ with n solvent molecules. Here we apply a less common definition (used, e.g., to calculate the binding energies of Li+-diglyme complexes in ref 3): c ∆Ec ) Ecompl - (ELi + En(solv) )

(1)

where Ecompl and ELi are the energies of the complex and of the c isolated lithium ion, respectively and En(solv) is the energy of n solvent molecules at the geometry of the complex. The advantage of this definition is that the fluctuations of the interactions between solvent molecules partially cancel between c ; in effect, the scatter in the complexation Ecompl and En(solv) energies calculated according to this formula is suppressed. Likewise, the problem with fluctuations of the distortion energy of solvent molecules, making the trend in ∆Ec less pronounced for increasing M, is eliminated. It should be noted in passing that the definition of eq 1 yields larger values of the complexation energy (especially for higher CN values) compared to alternative definitions where the energy of the isolated molecule at equilibrium geometry is used as a reference, because the destabilizing effects of the distortion in the complex and in the solvent molecules cancel out. This fact is of no relevance in the present context, since we are interested not in particular values of the binding energies for given CN values (which have been already calculated elsewhere1-3,24-27), but rather in the changes of the complexation energy when more solvent molecules are added to the system. Therefore, it seems reasonable to choose the present definition of ∆Ec, as it makes the trend most easily noticeable. To reduce the computational effort, we refrained from routinely correcting the complexation energies for the basis set superposition error. Instead, we performed test calculations using the counterpoise method for M ) 0, 2, 5 and for selected structures with the largest M values. The errors turn out to result mainly from the change of the Li+ energy upon addition of the solvent molecules basis sets; consequently, they should saturate when the number of molecules in the aggregate exceeds some threshold value, expected to be quite low. Indeed, the BSSE errors generally amount to less than ∼1.5 kcal/mol, and for M larger than 2 they cease to depend on the CN. Accordingly, inclusion of the correction would change the relative complexation energies by no more than 1 kcal/mol.

Figure 2. Average complexation energies ∆Ec for Li+ complexes with CN +M ACN (top) or DEE (bottom) molecules calculated at the B3LYP/6-31+G** level for increasing M. Insets: nonaveraged data for M ) 0 and for the largest M used in calculations; subsets corresponding to different CN values have been horizontally shifted for clarity.

In Figure 2 we present the changes in the complexation energy ∆Ec upon increasing the number M of ACN or DEE molecules beyond the first solvation shell. The values displayed for each CN are averaged over the whole series. The table containing the numerical values and standard deviations is available in the Supporting Information. The stabilization energy of the complex as a whole increases with M, as expected for Li+ solvation in an increasing amount of the solvent. It is readily noticeable that the differences between the average complexation energies for different CN are the largest for M ) 0 (complex in vacuum) and decrease when more solvent molecules are added. Although some deviations resulting from poor statistics for relatively short series are discernible (especially for ACN), the energy differences seem to saturate for M > 5. To give some information about the scatter of complexation energies for given CN, nonaveraged data for M ) 0 and for the largest M are displayed in the insets. For M ) 0 (standard “vacuum” calculations) the Li+-(ACN)2 complex is on the average about 75 kcal/mol less stable than the complex with the highest CN: Li+-(ACN)6. In the solvent this difference is reduced to about 15 kcal/mol. Likewise, the energy differences between Li+-(DEE)2 and Li+-(DEE)4 decrease in the solvent to ∼15 kcal/mol from 40 kcal/mol in vacuum. One may define the reduction factor of relative energies as the ratio of the spread of complexation energies in vacuum to the spread obtained for the largest value of M: R ) ∆Emax(M ) 0)/∆Emax(M ) Mmax), where ∆Emax is the difference between the average complexation energies for the minimum (CN ) 2) and the maximum coordination number (CN ) 6 for ACN and CN ) 4 for DEE). This factor is significantly larger for the more polar ACN solvent (about 5 for ACN compared to ∼3 for DEE). Although the differences between the mean values of the complexation energy decrease, the scatter of individual energies for given CN increases when more solvent molecules are added

976

J. Phys. Chem. A, Vol. 114, No. 2, 2010

to the system (c.f. insets in Figure 2 and Table S1 in the Supporting Information). This is no matter for concern and was to be expected. In fact, different arrangements of the solvent molecules beyond the first solvation shell can either stabilize or destabilize the complex; the larger their number, the larger the overall amplitude of the ensuing fluctuations which add to those within the first solvation shell. In principle, these two contributions could be separated by a judicious modification of the computational procedure: it would be possible to freeze the nearest CN solvent molecules in MD simulations while allowing others to move, then generate a series of MD trajectory snapshots, and average the complexation energy to obtain its mean value for a particular structure of the complex. These mean values, including the effect of the fluctuations of the solvent, would allow for straightforward comparison of the complexation energies not only for different CN, but also for different arrangements of the CN solvent molecules closest to Li+, giving additional insight into the microscopic physics of the problem. However, the advantage would be mostly illusory anyway, since the corresponding highly specific data are not accessible to direct experimental probing. Moreover, such a procedure would be computationally very demanding. Therefore, we restrict ourselves to the simpler approach presented above, which automatically accounts for the inevitable fluctuations both outside and inside the first solvation shell and without the necessity of referring to individual configurations directly produces the values of ∆Ec averaged for a given CN. It is worthwhile to note that calculations of this kind are feasible not only for cations but also for cation-anion complexes or even for larger ion-counterion aggregates. Ion pairs also exhibit solvation shell structure, with the difference that the perturbation induced by the counterion makes it more difficult for the solvent to enter the first solvation shell, reducing the CN values of individual ions (as found, e.g., for NaI in ACN14). The technical difference is that in order to properly reproduce the solvation of the whole complex a larger number of molecules has to be used in the calculations. The reduction factor obtained for cation-anion pairs is expected to be smaller, since the effect is mainly electrostatic and is bound to be smaller for an electrically neutral complex. By and large, quantum-chemical calculations for an ion complexed by a small number of nearest solvent molecules (M ) 0) are relatively fast (even with geometry optimization of the whole system) and can easily provide optimized structures and complexation energiessin fact, they are routine tasks in theoretical studies of ion complexation. Nevertheless, one should be careful with the relative binding energiesstheir values obtained from vacuum calculations are appropriate for gas-phase complexes but greatly overestimate the energy differences in a real liquid, while just those are usually of primary interest. The results of this section show that the relative binding energies may be reduced several times, depending on the dielectric constant of the solvent (e.g., about 5 times for ACN). The vacuum results are still of great value, but for proper description of ion energetics in solution they need to be rescaled. On this view, it would be tempting to take advantage of the speed of quantum-chemical calculations for an ion complexed to a small number of explicit solvent molecules and find a reasonable estimate of the scaling factor by means of an inexpensive implicit-solvent model. In the next section we probe this possibility. 2.3. Polarizable Continuum Model Calculations. Preliminary Test: SolWation Energy in PCM. To include the effect of implicit solvent we applied the Gaussian 03 implementation of

Eilmes and Kubisiak TABLE 1: Values of the Solvation Energy (In kcal/mol) for Li+ in Acetonitrile and Diethyl Ether Calculated within the PCM Approach at the B3LYP/6-31+G** Level for Different Models (Atomic Radii/Surface) Used to Build the Molecular Cavity surface radii

SES

UA0 UFF UAHF UAKS Pauling Bondi

-129.3 -129.3 -91.6 -91.6 -115.0 -84.4

UA0 UFF UAHF UAKS Pauling Bondi

-101.9 -101.9 -77.8 -77.8 -90.5 -66.2

VdW

SAS

-129.3 -129.3 -129.3 -129.3 -115.0 -84.4

-35.9 -35.9 -35.9 -35.9 -33.0 -24.8

-101.9 -101.9 -101.9 -101.9 -90.5 -66.2

-17.8 -17.8 -17.8 -17.8 -15.8 -9.8

ACN

DEE

the PCM methodology, using the default parameters for acetonitrile (ε ) 36.64) and diethyl ether (ε ) 4.335). There are many possible options to build the molecular cavity, which may lead to different solvation energies. Therefore, it is desirable to check the effect of different cavity parametrizations in order to pinpoint the parameter sets that are unlikely to yield correct results for relative complexation energies and to exclude them from further use in this context. A suitable test may consist in comparing the values of solvation energy for the lithium cation (i.e., energetic effect of transferring the ion from vacuum into the solvent) obtained in the PCM approach with those resulting from the explicit solvent model. We began this test with the PCM calculations of Li+ solvation energy for different cavity parametrizations. We tested all available combinations of the atomic radii with the type of the surface representing the solute-solvent boundary, leading to the total of 18 different specifications of the cavity. Among the atomic radii under consideration were three parametrizations with explicit hydrogen atoms, namely the Pauling, Bondi, or universal force field (UFF) radii, and three united atom topological models (where the hydrogen atoms are incorporated into the sphere of the heavy atom to which they are bonded), based on the UFF radii (UA0 - Gaussian 03 default), on the radii optimized for HF calculations (UAHF) and on the radii optimized for DFT (UAKS). The explored surface types were: solvent excluding surface (SES - G03 default), van der Waals surface (VdW) and solvent accessible surface (SAS). Table 1 shows the calculated values of the PCM corrections to the Li+ energy in the solvent, including all contributions (electrostatic and nonelectrostatic terms: cavitation, repulsion, and dispersion energies). As readily seen, the values of solvation energy strongly depend on the cavity parameters. For ACN the weakest stabilization amounts to about -25 kcal/mol (SAS surface, Bondi radii), whereas the strongest reaches -130 kcal/ mol. The corresponding values for DEE are -10 and -102 kcal/ mol, respectively. For both solvents the smallest solvation energies were obtained for SAS, because the molecular cavities built within this model are the largest. The strongest stabilization was obtained for SES or VdW surfaces combined with UFFbased atomic radii (UA0 or UFF). Our next step was to estimate the solvation energy of the Li+ ion from molecular level solvation calculations similar to those described in the preceding section. To this end we

Complexation Energies for Li+ Ion in Solution calculated quantum-chemically the total energies for the Li+ ion embedded in a cluster of 25 ACN molecules or (alternatively) of 18 DEE molecules, each for a set of 50 geometries randomly generated from MD trajectories. In the same way we calculated the total energies of the aggregates consisting of 25 ACN or of 18 DEE molecules without the lithium ion (each for a set of 50 geometries selected from supplementary MD simulations of the pure solvents). The solvation energy of the lithium ion was then calculated as the difference between the mean values of the energies of the Li+-solvent and the puresolvent systems. For ACN we obtained -118 kcal/mol; for DEE the solvation effect is smaller and amounts to -88 kcal/mol. As high accuracy was of no consequence in this context and we were merely seeking a reasonable estimate, we did not perform the time-consuming extrapolation to infinite system size and took the values obtained for Li+ dissolved in 25 ACN or 18 DEE molecules as the lower bounds for the actual solvation energies. The best agreement between our explicit-solvation estimates and the solvation energies calculated from the PCM model (Table 1) is obtained for the SES or VdW surface combined with Pauling atomic radii. However, this may well be an artifact resulting from the finite size of our systems; one should expect that if the infinite solvent bulk were properly accounted for, the solvation energy values would increase toward the results obtained for the SES (or VdW) surface with UA0 (UFF) radii or for the VdW surface with UAKS radii. Evidently, the values resulting from PCM calculations with the SAS surface are seriously underestimated, so that this parametrization is unlikely to perform well in the calculations of complexation energies. Based on the above test for solvation energies, we eliminated the SAS surfaces from further analysis, focusing on the most promising surface/atomic radii combinations, which are the following: SES/UA0 (the Gaussian 03 default), SES/UAKS, SES/Pauling, VdW/UFF (new default in Gaussian 09), VdW/ UAKS, and VdW/Pauling. Complexation Energy in PCM. With this prior knowledge, we turned again to our original task, which was the application of PCM methodology to reduce the computational effort involved in complexation energy calculations. We performed the PCM calculations for the same series of structures as used in Section 2.2 to calculate the energies for Li+-(solvent)CN complexes with M ) 0 additional solvent molecules. The complexation energies were calculated according to eq 1 with all PCM corrections (electrostatic and total nonelectrostatic) added to all energies. The binding energy values defined in this way are smaller than in vacuum, because in the solvent the stabilization of a free Li+ ion is larger than the stabilization of the complex. At this point, we deliberately make no extra allowance for the basis set superposition errors, assuming that they would not exceed the values calculated in the preceding section for molecular solvation. This assumption is based on the results existing for other systems, where BSSE calculated within the PCM approach do not differ significantly from the vacuum values.35 Usually in the study of ion complexation in a given solvent the relatiVe energies of the complexes with different CN are more relevant rather than their absolute values, therefore in Figure 3 we display the complexation energies (averaged over the whole series) relative to the most stable structure, that is, shifted in such a way that the average complexation energy for the complex with the highest CN is set equal to zero (the nonshifted value for the most stable complex is displayed below each data set, and the complete set of original results is available

J. Phys. Chem. A, Vol. 114, No. 2, 2010 977

Figure 3. Average complexation energies for Li+-(solvent)CN complexes obtained from the explicit-solvation approach and from the PCM calculations. PCM parameters are indicated as surface/radii combinations. Energies are relative to the complex with the largest CN, the absolute values of the mean complexation energy for CNmax (in kcal/ mol) are displayed below each data set. Symbol codes as in Figure 2.

in Supporting Information). These results are to be compared with the corresponding data from vacuum calculations (M ) 0) and with the results of molecular solvation study for the largest used M, which serve as a reference for complexation energies in solution. It is apparent that the PCM calculations reduce the spread in relative complexation energies compared to standard vacuum calculations. The reduction factor depends on the parameters of the model used to build the molecular cavity. For both solvents the reduction of the energy differences is smaller for the SES surface, and the smallest is obtained for the SES/UAKS surface/radii combination. Therefore, the SES calculations, while improving the results compared to vacuum calculations (M ) 0), still yield significantly overestimated relative complexation energies in the solvent. Much better agreement between the explicit-solvation and the PCM results is found for the VdW surface, especially for VdW/Pauling in ACN; for DEE the energy differences are still slightly too large. Figure 3 emphasizes mainly the gap between the smallest and the largest binding energy. To inspect the pattern of energy levels corresponding to different CN values it is convenient to rescale the energy values. We applied to the data of Figure 3 a scaling factor f ) ∆Emax(molec. solvat.)/∆Emax, with ∆Emax defined in Section 2.2 Therefore, by construction, for all scaled sets the energy gap between the complexes with CNmin and CNmax is the same as for the results of explicit solvation calculations (for which f ) 1 by definition). The rescaled data are displayed in Figure 4, and the scaling factor is given below each data set (the values of the scaling factor lower than 1 indicate that the energy differences are larger than in the reference molecular-solvation data). Although the pattern of energy values are presumably affected by the fluctuations that are not satisfactorily averaged out in the relatively short sampling series, it seems that in the SES

978

J. Phys. Chem. A, Vol. 114, No. 2, 2010

Eilmes and Kubisiak the differences in the discrepancies for different CNs subsiding for larger M values. In the PCM approach the differences between test MP2 results and B3LYP are even smaller. Apparently, application of the much more time-consuming MP2 methodology changes the relative complexation energies by no more than 2 kcal/mol. A similar effect may be expected for dispersion-corrected DFT functionals. Evidently, the accuracy offered by the less sophisticated approximations used in this paper is perfectly sufficient for our present purposes. 3. Conclusions

Figure 4. Relative average complexation energies rescaled to reproduce the energy differences obtained from the explicit-solvation approach (CN + Mmax data, for details see text). The appropriate scaling factor is given below each data set. Symbol codes as in Figure 2.

model the energies for the Li+-ACN complexes with small CN are shifted upward, the effect being most pronounced for SES/UA0. A similar shift is observed in the SES/UA0 data for DEE. Noticeably, the scaled energy differences obtained for the VdW surface (for all types of atomic radii) agree quite well with the reference pattern of explicit-solvation results. This concurs with the fact that the scaling factors for these three parameter combinations are the closest to unity, suggesting that in the PCM calculations the VdW surface with Pauling or UFF radii are the most promising for obtaining reasonable estimates of relative binding energies is solution. A remarkable difference between the best PCM results and the explicit-solvation approach is the scatter of complexation energies. As already mentioned in Section 2.2, complexation energies for a given CN, calculated in the explicit-solvent model, vary in quite broad range. The dispersion of the PCM results is smaller, especially for ACN, where it is reduced 3-5 times (the numerical values and the plot for VdW PCM calculations are available in Supporting Information). This is a consequence of the fact that in order to obtain the mean value of the complexation energy, the explicit-solvation results have to be averaged over individual configurations (each of which may yield a significantly different energy), whereas in the PCM approach the averaging is already included in the effective solvent. All quantum-chemical calculations in this work employed density functional theory (DFT) with the B3LYP functional. This provokes the question how the results would change if other computational methods were used, especially those taking dispersion into account. As a test we performed MP2 explicitsolvation calculations for selected ACN (M ) 0, 2, 10) and DEE (M ) 0, 2, 5) complexes. For ACN the stabilization energies calculated within the MP2 method are 1 - 2 kcal/mol smaller than the B3LYP values, with no major differences in the discrepancies for different coordination numbers (for all tested M). In the case of DEE, MP2 predicts the stabilization about 1-3 kcal/mol stronger than the DFT calculations, with

We have presented a study of relative complexation energies of Li+ ion in two commonly used solvents, ACN and DEE. Combining MD simulations of the structure of ion-solvent complexes with quantum-chemical calculations we have obtained average complexation energies for different CN values. Molecular-level representation of solvation, embodied in explicit solvent molecules surrounding the complex, has allowed us to determine the trends in relative binding energies. Our results indicate that with increasing number of included explicit solvent molecules the differences in complexation energy between the complexes of different CN decrease. The changes in relative binding energies seem to saturate above 10 explicit solvent molecules used in the calculations. Overall, the surrounding solvent reduces the complexation energies, as intuitively expected. The reduction factor depends on the dielectric constant of the solvent and is larger for more polar solvents. For ACN, it reaches 5, whereas for less polar DEE it is about 3. This scaling must inevitably be taken into account when the results of quantum-chemical calculations in vacuum are to be used for drawing conclusions about ion complexation in solution. The explicit-solvation approach allows one to examine thoroughly the relation between the structure of the ion-solvent complex and the changes in the binding energies that result from increasing number of molecules and from their varying arrangement around the ion. Such calculations require initial MD simulations to obtain information about the structure of the solvated adduct, and subsequently quantum-chemical calculations for the systems which, by these standards, are relatively large (at least 10-15 solvent molecules). Consequently, solvation calculations at molecular level are very time-consuming, and for routine studies of ion complexation they seem less feasible than the standard vacuum calculations limited to an ion complexed to a small number of nearest solvent molecules. Seeking a methodology that could provide information about the relative binding energies in solution at a cost comparable to quantum-chemical computation in vacuum, we have tested the PCM description. This approach, combining explicit solvation with implicit-solvent model, had been previously applied to study various properties, for exmple, the solvent effect on NMR shielding constants,21 solvent-induced shifts in electronic22 or vibrational23 spectra, or solvent effect on reactivity indices.36 Our present work investigates the performance of the method for relative complexation energies. In our PCM study we have used various combinations of atomic radii and molecular surface modeling the molecular cavity in the solvent. For all parametrizations we have observed a decrease of the complexation energy differences compared to the results of vacuum calculations, albeit less pronounced in most cases than predicted by the explicit-solvation model. The best reproduction of the energy scaling factor and the structure of energy levels for different coordination numbers has been found for van der Waals molecular surface.

Complexation Energies for Li+ Ion in Solution The observed effects result mainly from the electrostatic interactions of the ion with the surrounding molecules (which is the cause of their dependence on the dielectric constant of the solvent). Accordingly, electrostatic PCM corrections to the energy are the most important contribution, whereas the nonelectrostatic terms are much smaller and partially cancel. This indicates that also other implicit-solvent models like COSMO37 (which is essentially electrostatic) should lead to similar conclusions. The electrostatic nature of the effect implies that similar results could be obtained for watersthe most common inorganic solventsas the electrostatic interactions are known to be predominant for ion hydration.38 The high dielectric constant of water, exceeding the ACN value more than two times, should produce an even larger reduction factor. However, indiscriminate application of this analogy may be misleading, since, for example, hydrogen bonds may influence the balance. In fact, for increasing number of water molecules explicitly hydrating the ion, the hydrogen-bonding interactions compete with the ion-water interaction.39 On the other hand, for larger ions the contribution from dispersion interactions increases,38 introducing a new relevant factor. Effectively, modeling of ion hydration is more challenging than the case of organic solvents. One would need to perform MD simulations with the force field properly describing the water molecules and their interactions with the ions, and then to apply high-level quantum-chemical methodology,40 making allowance in the implicit-solvent parametrization for hydrogen bonding and dispersion interactions in water. Our results, exemplified on two organic solvents, demonstrate that PCM calculations with van der Waals surface and Pauling or UFF atomic radii are capable of providing low-cost estimates of the relative complexation energies in solution that are comparable to the results of the time-consuming explicitsolvation approach. In fact, the gain in the computational time is substantial: PCM calculations took on the average about 1.5 times longer than the vacuum calculations, but were at least 60-200 times faster than molecular-solvation calculations (for the largest used M). In this context, PCM computations emerge as an extremely useful theoretical tool, substantially enhancing insight into the solvent effect on energetics of ion complexation, and ultimately contributing to a better understanding of ionic processes in liquid phase. Acknowledgment. The authors express their gratitude to Professor P. Petelenz for valuable comments on the manuscript. The Gaussian 03 calculations were performed in the ACK “Cyfronet” Computing Centre (computational grants MNiSW/ SGI3700/UJ/070/2009, MNiSW/SGI4700/UJ/070/2009 and MNiSW/IBM_BC_HS21/UJ/070/2009). Supporting Information Available: Tables of average complexation energies (with standard deviations) for Li+ in ACN and DEE. Plot of individual (nonaveraged) energies for selected data sets. This material is available free of charge via the Internet at http://pubs.acs.org/. References and Notes (1) Johansson, P.; Gejji, S. P.; Tegenfeldt, J.; Lindgren, J. Solid State Ion. 1996, 86-88, 297. (2) Sutjianto, A.; Curtiss, L. A. J. Phys. Chem. A 1998, 102, 968.

J. Phys. Chem. A, Vol. 114, No. 2, 2010 979 (3) Johansson, P.; Tegenfeldt, J.; Lindgren, J. Polymer 1999, 40, 4399. (4) Johansson, P.; Edvardsson, M.; Adebahr, J.; Jacobsson, P. J. Phys. Chem. B 2003, 107, 12622. (5) Blint, R. J. J. Electrochem. Soc. 1997, 144, 787. (6) Johansson, P.; Jacobsson, P. J. Power. Sources 2006, 153, 336. (7) Johansson, P.; Jacobsson, P. Solid State Ion. 2006, 177, 2691. (8) Johansson, P. Phys. Chem. Chem. Phys. 2007, 9, 1493. (9) Eilmes, A.; Kubisiak, P. J. Phys. Chem. A 2008, 112, 8849. (10) Eilmes, A. Solid State Ion. 2008, 179, 458. (11) Borodin, O.; Smith, G. D. J. Phys. Chem. B 2003, 107, 6801. (12) Borodin, O.; Douglas, R.; Smith, G. D.; Trouw, F.; Petrucci, S. J. Phys. Chem. B 2003, 107, 6813. (13) Borodin, O.; Smith, G. D.; Douglas, R. J. Phys. Chem. B 2003, 107, 6824. (14) Nguyen, T.-N. V.; Hughes, S. R.; Peslherbe, G. H. J. Phys. Chem. B 2008, 112, 621. (15) Jarek, R. L.; Miles, T. D.; Trester, M. L.; Denson, S. C.; Shin, S. K. J. Phys. Chem. A 2000, 104, 2230. (16) Peschke, M.; Blades, A. T.; Kebarle, P. J. Am. Chem. Soc. 2000, 122, 10440. (17) Megyes, T.; Gro´sz, T.; Radnai, T.; Bako´, I.; Pa´linka´s, G. J. Phys. Chem. A 2004, 108, 7261. (18) Rodgers, M. T.; Armentrout, P. B. J. Phys. Chem. A 1999, 103, 4955. (19) Salter, T. E.; Ellis, A. M. J. Chem. Phys. 2007, 127, 144314. (20) Rao, J. S.; Sastry, G. N. J. Phys. Chem. A 2009, 113, 5446. (21) Aidas, K.; Møgelhøj, A.; Kjær, H.; Nielsen, Ch. B.; Mikkelsen, K. V.; Ruud, K.; Christiansen, O.; Kongsted, J. J. Phys. Chem. A 2007, 111, 4199. (22) Meng, S.; Ma, J. J. Phys. Chem. B 2008, 112, 4313. (23) Mennucci, B.; da Silva, C. O. J. Phys. Chem. B 2008, 112, 6803. (24) Cabaleiro-Lago, E. M.; Rı´os, M. A. Chem. Phys. 2000, 254, 11. (25) Gre´goire, G.; Brenner, V.; Millie´, P. J. Phys. Chem. A 2000, 104, 5204. (26) Pejov, L. Int. J. Quantum Chem. 2002, 86, 356. ˇ ernusˇa´k, I.; Aranyosiova´, M.; Volla´rova´, O.; Velicˇ, D.; Kirdajova´, (27) C O.; Benko, J. Int. J. Quantum Chem. 2009, 109, 2365. (28) Cabaleiro-Lago, E. M.; Rı´os, M. A. Chem. Phys. 1998, 236, 235. (29) Gua`rdia, E.; Pinzo´n, R. J. Mol. Liq. 2000, 85, 33. (30) Spångberg, D.; Hermansson, K. Chem. Phys. 2004, 300, 165. (31) Tinker Molecular Modeling Package, v. 4.2; http://dasher.wustl.edu/ tinker/. (32) The set of parameters distributed with the Tinker program was used. The Tinker MM3(2000) parameters were provided by Allinger N. L. , University of Georgia. (33) Mao, G.; Saboungi, M.-L.; Price, D. L.; Badyal, Y. S.; Fischer, H. E. Europhys. Lett. 2001, 54, 347. (34) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, ReV. E01; Gaussian, Inc.: Wallingford CT, 2004. (35) Cammi, R.; Olivares del Valle, F. J.; Tomasi, J. Chem. Phys. 1988, 122, 63. (36) Jaramillo, P.; Pe´rez, P.; Fuentealba, P.; Canuto, S.; Coutinho, K. J. Phys. Chem. B 2009, 113, 4314. (37) Klamt, A.; Schu¨u¨rmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 799. (38) Lee, H. M.; Tarakeshwar, P.; Park, J.; Kołaski, M. R.; Yoon, Y. J.; Yi, H.-B.; Kim, W. Y.; Kim, K. S. J. Phys. Chem. A 2004, 108, 2949. (39) Miller, D. J.; Lisy, J. M. J. Am. Chem. Soc. 2008, 130, 15381. (40) Rao, J. S.; Dinadayalane, T. C.; Leszczynski, J.; Sastry, G. N. J. Phys. Chem. A 2008, 112, 12944.

JP907359A