Ionic Liquid Structures from Large Density Functional Theory

Sep 10, 2010 - Dhumal , N. R.; Kim , H. J.; Kiefer , J. J. Phys. Chem. A 2009, 113 ...... SECIL - Software Engine for Configuration of Ionic Liquids. ...
0 downloads 0 Views 196KB Size
Download PDF
J. Phys. Chem. C 2010, 114, 20577–20582

20577

Ionic Liquid Structures from Large Density Functional Theory Calculations Using Mindless Configurations† Knut Angenendt and Patrik Johansson* Department of Applied Physics, Chalmers UniVersity of Technology, SE-41296, Go¨teborg, Sweden ReceiVed: May 31, 2010; ReVised Manuscript ReceiVed: August 20, 2010

Three different popular imidazolium based ionic liquids (ILs); EMI-BF4, EMI-PF6, and EMI-TFSI, have been modeled by DFT calculations (B3LYP/6-311+G*) using large, up to 130 atom cluster models, for a better understanding of the structure and ion-ion interactions in these ILs and ILs in general. Particular emphasis has been put on the role of appropriate starting structures and how the present large models differ from the ion-pair models of ILs generally used. The system size normalized ion-ion interaction energies are shown to converge rapidly, and conformational equilibria and higher order properties like IR spectra are shown to be valuable as quality criteria. The explicit inclusion of an IL environment by the large cluster approach is also compared to using an implicit, continuum, strategy via SCRF C-PCM calculations. Introduction Ionic liquids (ILs), i.e., materials liquid below 100 °C despite being composed only of ions, are currently a very active field of science.1,2 This is primarily due to these materials having a combination of physical and chemical properties interesting for extremely varied applications such as solvents for synthesis with special templating properties,3 as solvents/electrolytes for lithium batteries and supercapacitors,2,4,5 as drugs,6 etc. However, scientific interest is also spurred by the lack of understanding of ILs at the most fundamental molecular level, e.g., the ion-ion interactions, ultimately responsible for the observed properties.7–9 One way to address this level of IL structure is by molecular simulations, where indeed methods like classical molecular dynamics (MD) has been applied for many systems,8–17 but also quantum mechanical (QM) calculations, predominantly ab initio or density functional theory (DFT), have been used to a very large extent.7,8,16–42 Methods like classic MD are powerful for dynamics, statistics, and ensemble properties, but the creation of accurate enough force fields for ILs is far from straightforward,10 and the combination of QM-MD is quite computationally demanding; therefore, studies are often limited to a single IL and ILs composed of small anions like Cl-.16,43,44 On the other hand, QM methods are often applied to quite reduced models of the size of only tens of atoms. This can still be useful to directly assist the interpretation of specific experimental observations for liquid and disordered systems with predominant short-range order, as demonstrated for similar systems to ILs (ionic glasses).45 We note, however, that quite recently also the use of QM methods like FMO/EFP, which are inherently aimed toward larger systems, have been demonstrated useful for ILs.46 Despite the simplicity, even gas-phase models of isolated ions have been successfully applied to connect with and explain experimental data for ILs.35–39,41 More common, however, is the usage of 1:1 cation-anion ion pairs as IL models, with the vast amount of work targeting specific interactions and properties.16–26 There are very few ion-pair ab initio/DFT studies directed at explaining structure-property connections for ILs in general. From a general point of view, it can even be argued that gas†

Part of the “Mark A. Ratner Festschrift”. * Corresponding author. E-mail: [email protected].

phase 1:1 cation-anion ion pairs are the worst possible molecular model choice to address bulk IL properties, even worse than isolated ions. The reason is that ion pairs introduce highly local and anisotropic interactions, while isolated ions do not and thus may better represent the reduced cation-anion interaction strength and directionality in ILs because the longrange Coulombic screening in ILs arguably is strong. Nevertheless, if progress in our understanding of ILs is to be made, the cation-anion interaction must be addressed in some way, but perhaps not by ion pairs. Two alternatives are to include the IL ionic environment either implicitly by, e.g., continuum methods,29 or explicitly, by the use of a “super-molecular” (cluster) approach (although this is a badly defined concept for pure ILs as there is no solvent/ solute difference). Looking at the cluster approach, hitherto, we know of only a few such studies.16,30–33 Most prominently, Ludwig et al. have been using IL clusters with up to 16 ion pairs (N ) 16).30–33 Additionally, some authors have used IL dimers (N ) 2),28,29 which seem better than simple 1:1 ion-pair models for experimental data interpretation. A second factor, to be added to the prime factor of atomic system size/number of IL units, is the compromise between choosing these ion pairs, dimers, or larger constructs as models, and the accuracy of the QM computational chemistry algorithm possible to apply in practice. Only a few groups have focused on the absolute best choice of method/algorithm for ILs and then primarily argued the use of ion pairs and an energetic point of view. There are those who advocate for MPx or higher ab initio methods,16,22,23 and others in support of DFT with dispersion included.20 Still, however, most studies use either simple Hartree-Fock (HF) or regular DFT methods to model ILs. The trade-off is sometimes very noticeable, as for the large clusters by Ludwig et al. where their studies were performed at HF level using very small basis sets with no polarization functionality.30–33 However, we here stress in addition a third factor, one that we feel is sometimes neglected: the role of the starting guess(es). As ab initio/DFT calculations almost exclusively are performed as deterministic searches for the (local) minimum energy structure(s), it is very important to emphasize and realize that the starting guess(es) have a direct influence on the results

10.1021/jp104961r  2010 American Chemical Society Published on Web 09/10/2010

20578

J. Phys. Chem. C, Vol. 114, No. 48, 2010

obtained. For ILs this notion is especially important for several reasons: First, many IL cations and anions can have stable conformations7,19–22,36,38,47 and, assembled together, configurations19–24,29–31,33 (this is indeed partly the reason for ILs being liquid); the ions are difficult to pack efficiently, and the importance of the complex potential energy surfaces (PES) has been proven.8,9 Therefore, even for ion pairs there might be many stable configurations that should be properly taken into account by the starting guesses. Of course, for larger IL models this becomes even more important as the complexity of the PES increases. Second, in the IL modeling literature, strikingly often the method or rationale of constructing the starting guesses is omitted or nonexistent even for the simple systems of separate ions, the conformational space,18,20,25,26 or the case of larger structures, also the configurational space.18,20,26,30,31 Even if well described, it is usually a limited and most often a directed effort by the investigator,7,19,21–24,28,32,33,36,39 and also it very often lacks experimental validation. To circumvent these issues, we here use a new code, SECIL,48 especially designed for automated and “mindless”49 construction of starting configurations of ionic systems and creation of input files for transfer to QM codes. The main algorithm in SECIL rotates and translates each ion randomly, with the only criteria being a fit based on a single minimum atom-atom distance between ions of opposite charge. The only input required is the internal coordinates for the separate IL ions and the number of ions. While the configurations produced by SECIL are far from any minimum energy structures, and thus the subsequent computational effort increases, the large gain is the automated and “mindless” sampling. With all three factors above in mind, SECIL generated configurations are here created for three popular imidazolium based ionic liquids. The scope is to present the strategy and to outline future qualitative analysis possibilities for ILs in general via selected examples. The largest clusters have no less than 130 atoms, which still makes our chosen DFT method, B3LYP,50–52 with a sufficiently large basis set to both accommodate polarization and diffuse electron densities, 6-311+G*, a feasible alternative with the advantages of reasonable accuracy, data to compare with the literature, and the possibility to reliably generate IR and Raman data to compare with experiments. The last is important for experimental validation of the IL local structures.

Angenendt and Johansson

Figure 1. Normalized binding energies, ∆Ebind/N, as a function of the configuration system size N for the three ILs at the two computational levels: (a) HF/6-31G*; (b) B3LYP/6-311+G*. Only the most stable configuration for each IL and value of N is included.

amounts to between 230 and 1100 core hours (Xeon E5430, 2.66 GHz).54 In total, 75 DFT optimized structures for N > 1 were obtained using >100 000 core hours. Due to the inherent randomness of the SECIL configurations only for ca. two-thirds it is possible to obtain initial SCF convergence; a ratio only slightly lower for the larger values of N. However, once an initial SCF convergence is achieved, very few geometry optimizations fail, and all minima obtained at the HF level produce minima also at the B3LYP level. The influence of the surrounding for 1:1 ion pairs has also been modeled implicitly by applying a continuum method added to B3LYP/6-311+G*, SCRF C-PCM,55,56 with parameter sets for both dichloroethane (DCE), ε ) 10.4, and water ε ) 78.4, respectively, to obtain the free energy of ion pair solvation (∆G) vs the sum of solvating the ions separately. SECIL was used to generate 10 ion pair configurations per IL. Similarly, the binding energy, ∆Ebind, is defined as the electronic energy difference with respect to the sum of the energies of the separate ions in the gas phase, and appropriately normalized by N, ∆Ebind/N. For a few selected structures, of different N and ILs, and for C-PCM ion pairs, IR and Raman spectroscopy data were computed via the Hessian by the second and third derivatives of the energies at the B3LYP/6-311+G* level in the harmonic approximation. No large imaginary frequencies were obtained, which is reassuring because it is not possible to calculate the Hessian for all obtained structures due to the computational demands.

Computational Methods SECIL48 was used to create starting configurations for three different ILs; EMI-BF4, EMI-PF6, and EMI-TFSI, where EMI is 1-ethyl-3-methyl imidazolium and TFSI is bis(trifluoromethanesulfonyl)imide. For the former two, using the nomenclature (cation:anion)N, Nmax was 5, while for the latter Nmax was limited to 3 for computational reasons. This results in configurations of maximum 120, 130, and 102 atoms, respectively. As a starting point, only the nonplanar (np) conformer was used for the EMI cation38 and the C2 symmetry (C2 or trans) conformer for the TFSI anion.38,47 All ion starting geometries were structures optimized at B3LYP/6-311+G*. After the automated generation of starting configurations and input files directly from SECIL, all configurations were transferred to Gaussian0353 calculations at supercomputers for optimization at two levels; ab initio HF/ 6-31G* and subsequently DFT B3LYP/6-311+G*. No BSSE correction was applied for the present purpose of only qualitative comparisons and due to the large Coulombic interactions present. As a measure of the computational effort needed, the B3LYP step for a configuration of EMI-BF4 with N ) 3 (72 atoms)

Results and Discussion The main aim is as outlined above to show the importance of including an appropriate procedure for generating starting structures to finally obtain configurations that mimic the ILs. Therefore, not only do we show here how the convergence with system size is affected by the configuration sampling, including how large errors can be expected when a single ion pair or even a dimer configuration is used as the measure, but also, to be a product of the computational level, we show how an average property like density changes, as well as detailed properties like specific distances and conformational equilibrium, and finally we probe how the presented configurations can explain experimental data from IR and Raman spectroscopy (features not possible to mimic using ion pair models and also highly sensitive to the configuration of the IL model). 1. Configuration and ∆Ebind Convergence and Size. 1.1. HF/6-31G*. As we will compare differently sized clusters to find convergence with system size, a normalization procedure is needed. In Figure 1 normalized maximum binding energies,

Large Ionic Liquid DFT Models ∆Ebind, as functions of cluster system size are represented for all different configurations. Let us first comment briefly on the fact that the normalized binding energy increases with N —should it not decrease if the ion pair models affect the ions too strongly? The reason for the opposite is due to the fact that for larger N the ions can interact with several other ions, thus the number of interactions increases more rapidly than the system size as measured by N, and thus the total binding energy increases, while per interaction it indeed is lowered. As a simple example, for N ) 2 we divide the obtained maximum binding energy by 2, but each anion will interact with both cations and vice versa. Thus the number of ion-ion interactions is more likely to be 4. Due to the complexity of the ILs, however, we have chosen not to try to determine whether there are interactions between the ions, as this becomes very difficult when N is larger, but the change from a highly anisotropic to a more isotropic interaction scheme should be beneficial from a more realistic IL structure point of view. Thus it is important to distinguish the convergence of system size from ion-ion interactions. From Figure 1a it is quite clear how much the ion-ion interaction can be overestimated by the use of ion pairs as models (N ) 1), while N ) 5 seems to give system size convergence within a few kJ mol-1 with respect to N ) 4-5, for both EMI-BF4 and EMI-PF6, converging to -450 and to -416 kJ mol-1, respectively. Again, note the difference between increasing binding energy as a function of system size and ion-ion interactions most likely decreasing. A few kJ mol-1 is the approximate difference between the most stable and the second most stable configurations, as a measure of uncertainty, and thus a further increase in N, with huge amounts of computational effort accompanying, can hardly be motivated on these grounds. The order of coordination strength between the different ILs is constant and furthermore also the percentage difference in ∆Ebind, the EMI-BF4 IL having the strongest total ion-ion interaction and EMI-PF6 (91.5-93.5%) and EMI-TFSI (88.5-89.5%) of this value, respectively. As it is highly beneficial, if possible, to rely on lower computational levels, we will next correlate the quantitative results from the HF/631G* calculations with our final DFT B3LYP/6-311+G* calculations. 1.2. B3LYP/6-311+G*. As mentioned in the Computational Methods, for all the configurations where the HF/6-31G* level produced a minimum energy structure, there was subsequently successfully a B3LYP/6-311+G* optimized structure obtained. In Figure 1b the obtained maximum binding energies, ∆Ebind, as functions of cluster system size are represented for all different configurations. There are some very noticeable differences toward Figure 1a; the binding energies decrease in absolute value and the convergence with N from 4 to 5 is better. There are also some differences toward Figure 1a that makes this next, and more accurate, level, if not needed, at least recommended; the spread in binding energies for the different values of N is slightly reduced, now being 92.9-94.5% for EMIPF6 and 91.2-92.4% for EMI-TFSI. The more important feature, however, is that the difference between the different configurations for the same choice of IL and N decreases. The importance of using several starting guesses is still emphasized by the vastly different ∆Ebind that can be obtained. For the DFT calculations the largest differences are obtained for the N ) 2 configurations, with the least stable configurations being 13, 33, and 19 kJ mol-1, the same order of ILs as above, higher in energy than the minimum, a difference on the order of 3-9%. This difference, however, quickly becomes smaller for increas-

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20579 ing N. Thus, this should serve as a warning for using dimer models for ILs without first making a proper scan of the configurational space for the system. Another set of comparisons is the evolution of the binding energies with N for each IL. For EMI-BF4 the value is down to -416 kJ mol-1 at N ) 5 and convergence, via N ) 3, is at -408 kJ mol-1 (98%) compared with the value for the ion pair model at N ) 1, -351 kJ mol-1 (84%). For EMI-PF6 the value is down to -393 kJ mol-1 at N ) 5 and convergence, via N ) 3, is at -385 kJ mol-1 (98%) compared with the value for the ion pair model at N ) 1, -326 kJ mol-1 (83%). For EMI-TFSI the value at N ) 1, -320 kJ mol-1, is 85% of the value at N ) 3, -377 kJ mol-1; a comparison gives 86 and 85%, respectively, for the two former ILs. Thus the method does not bias any particular chemistry out of these three ILs for faster or slower convergence. On the opposite, if we only are interested in comparing the relative total ion-ion interaction strengths, the ion pair model is now verified to be a good alternative, of course still with the demand of a proper starting configuration sampling. As it has been pointed out previously that B3LYP is not the best functional to obtain the best absolute binding energies,22 we again stress that we are here looking at trends for N and between only three different ILs. Nevertheless, with B3LYP the same order of anion dependence as from the MP2 calculations in ref 22 is obtained and also the anion dependence percentages with the BF4- anion as baseline, using the most stable C1mpyr values of ref 20, are not far apart (94 vs 93% and 92 vs 91%). Other studies using the ion pair approach obtain binding energies for EMI-BF4 in the range 350-375 kJ mol-1, HF and MP2 with various basis sets, and about 92% of that for EMI-PF6.19 However, the exact values for all studies using ion pairs should depend on the starting configuration scheme. Further comparing with the literature, there are very few studies on convergence, as there are few studies on larger models than ion pairs. Ludwig et al. mainly focus on the correlation of binding energies with either vibration frequencies30 or thermodynamic data,31 but also show slow convergence for the IL MMI-SCN, with convergence first for N ) 12 and maximum binding energies earlier, for N ) 10.31 This might be a consequence of the HF/3-21G method, with its very small basis set. They do, however, correlate all clusters toward two different choices of ion pairs and their energies, and not to the individual ions. The latter concept is applied in ref 32 where charged protic IL clusters are calculated and convergence within a few kJ mol-1 is obtained for N ) 9 and higher, but with a plateau already for N ) 4 and 5, thus in the region of N where we find convergence. Thus, as our calculations are not extended further than this plateau, we cannot be absolutely sure if our approach is converged in this sense, but as neutral systems arguably behave differently than charged ones and aprotic ILs are different from protic ILs, no strong conclusions based on this should be made. A better comparison can be made with the study on the MMICl IL where approximate binding energy convergence within 1-2 kJ mol-1 vs ions is found for N ) 6 and up (Table XI in ref 34), although a continuing increase is seen for even higher values of N. However, the total increase in binding energy with N is much smaller than for the here presented ILs, only 20 kJ mol-1 from N ) 1 to the “infinite” sized, probably due to the monatomic Cl- anion of the IL. To conclude on energy studies, it is difficult to find straightforward literature comparisons, but the convergence seems faster than or equal to that in previous studies, and in addition, the relative total ion-ion interaction strength using simple ion pair models has been verified by the larger models.

20580

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Angenendt and Johansson

Figure 2. One configuration of EMI-PF6, N ) 3, as obtained: (a) directly from SECIL; (b) as first optimized by HF/6-31G*; (c) as subsequently optimized (B3LYP/6-311+G*). The IL density increase was ∼5% and ∼3%, respectively.

2. Configuration Geometries and IL Structure. One purpose of the present study is to briefly discuss the role of the configuration generation for the final overall IL structures. As an example, in Figure 2a-c the changes of a single configuration of EMI-PF6, with N ) 3 and 78 atoms, between three different computational stages are visualized: directly as obtained from SECIL, as optimized using HF/6-31G*, and finally as the B3LYP/6-311+G* optimized structure. From this example, there are a few details to be extracted that pertain to the entire study. First, there is an increase in density from Figure 2a to Figure 2b, by ca. 5% as evaluated by the 12 largest xyz values with respect to origo, and thus the generated structures from SECIL can be thought of as more gas-phase like or at a higher temperature. There is also a quite noticeable rearrangement of cations vs anions. This is highly beneficial to enable the SECIL generated configurations to locate well-separated minima, although time-consuming. Second, from Figure 2b to Figure 2c there is also a density increase, by 3%, but less of an overall structural change. Thus, from this point of view, the last step is a refinement. Third, the B3LYP step might nevertheless be crucial for a correct description of the physical properties of the IL. A structural and physical properties example, an analysis to be more elaborately treated in the future,57 is the very often discussed role of the most acidic proton of the EMI cation, situated at the C2 carbon, (C2)H, for the ion-ion interactions, especially whether there is hydrogen bonding or not.14,15,18,19,21,25,27,29,30,32,34,36,42 For the connection of the (C2)H atoms and the F atoms of PF6- there is a small, but significant, difference between the HF and B3LYP levels of calculation. At the HF level there are two short (C2)H · · · F distances, 2.12 and 2.28 Å, notably belonging to the same PF6- anion, and the subsequent B3LYP optimization renders these even shorter, 2.07 and 2.10 Å, respectively. This type of bridging by an anion can of course not be obtained from studies of 1:1 ion pairs (where we do find one ∼2.1X Å distance) but may have large implications on the IL properties. As a second example, which also exemplifies the importance of the starting structures, both the EMI cation and the TFSI anion are well-known to have two stable conformers38,47 (p and np and C1 and C2, respectively) at room-temperature in the EMITFSI IL, with the C2 and np conformations of TFSI and EMI dominating, 75% and 87%, respectively.38 For EMI-TFSI and N ) 2, conformational data for 2 × 2 × 10 ions are available for analysis from the B3LYP/6-311+G* level. As the SECIL generation used only the two most stable conformations, C2 and np, any new conformation for either TFSI or EMI is a sign of a less rigid path toward the minimum energy structure than usually encountered. For TFSI, there is a distinct distribution of conformers obtained with 12 C2 and 7 C1, and only one TFSI anion in state between (the C-S-N-S dihedral angles being +96° and -143°). The C2 conformers have C-S-N-S dihedral

angles within 85-130°, and many in a state close to what is obtained for an isolated anion, i.e., both ca. 90°. The C1 conformers all have one C-S-N-S dihedral larger than 163° and the other ∼100°, as expected from the vacuum PES.47 Thus, even starting with only the C2 conformer, a slightly lower amount of C2 vs C1 (60/35) than experiment (75/25) is obtained.38 For EMI there is a worse agreement with experimental data: no conversion from the np to the p (planar) conformer is obtained; all EMI keep very rigid. Thus, for modeling of EMI-based ILs with the correct conformational state of EMI, a direct input to SECIL of ca. 1/10 of the p conformer is necessary. The present results can be a sign of a higher barrier for conversion between np and p conformers than for C2 and C1, as the relative energy differences are very similar.38 3. Explicit vs Implicit Models. To model the local ion environment within an IL there are two main strategies: implicitly and explicitly. For average properties the implicit strategy applied to, e.g., ion pairs may serve excellently for many purposes. However, for properties like the free energies of solvation, ∆G, as obtained from the C-PCM ion pair calculations (not shown), about the same maximum value (∼25 kJ mol-1) is obtained regardless of the choice of IL and solvent parameters. If any IL ion pair shows a stronger ion-ion interaction, it is the EMI-PF6 in water, but only with a few kJ mol-1. Thus, this is in contrast to the reliably predicted order for all explicit models: EMI-BF4 > EMI-PF6 > EMI-TFSI. However, as we now also include parameters for the solvents any direct energy comparisons cannot be made, as they are on different scales. Nevertheless, implicit models are arguably an excellent way to in an efficient way model an isotropic environment, something to be expected in ILs at large. A property highly sensitive to anisotropy is the IR and Raman selection rules for ions with high symmetry; thus this can be used as a probe of ion-ion interactions in ILs. Recently, in the IR spectrum of EMI-PF6, the fully symmetric stretching ν1 of the octahedral PF6- anion was observed.58 This is interesting as the A1 g mode (at ∼740 cm-1 experimentally) according to the selection rules should be IR-inactive. Similarly, the δ(BCN) bending mode of the tetrahedral B(CN)4- anion has been observed in the Raman spectrum of the EMI-B(CN)4 IL, although this mode should be Raman-inactive for a strict Td symmetry.37 Thus, both these anions must be at least slightly perturbed in the ILs, but there are no noticeable splits or accompanying shifts, and thus the interactions of, e.g., the PF6- anion with the EMI cations are weak enough to preserve at least a quasi-octahedral symmetry and the wavenumbers of a “free” anion. Can these fundamental observations be corroborated by either implicit or explicit models and, conversely, be used as a test of the modeling approaches? In Figure 3 artificial IR spectra generated from the calculations are shown for the region for the ν1 PF6- mode. The computed frequency for the PF6- anion in a vacuum is ∼678 cm-1 and

Large Ionic Liquid DFT Models

Figure 3. Comparison of artificial IR spectra for EMI-PF6 in the region of the PF6- A1g mode obtained from (a) N ) 1, (b) N ) 2, (c) N ) 3, and (d) N ) 1 with C-PCM parameters for DCE and (e) N ) 1 with C-PCM parameters for water. The calculated intensities are convoluted by a Gaussian of FWHM ) 2 cm-1. The vertical solid and dotted lines depict the zero-intensity position for a PF6- anion in vacuum and using C-PCM (within 0.5 wavenumbers for both DCE and water), respectively. The asterisks indicate an internal EMI mode. The spectra have been arbitrarily shifted along the y-axis for visuability.

for both C-PCM parameter sets ∼686-7 cm-1, with no scaling factors applied to the computed data. Experimentally this mode is superimposed on a broader absorption due to the EMI cation,58 but in the artificial spectra this EMI band appears distinct to the left. From Figure 3a it is clear that a simple gas-phase 1:1 ion pair model indeed produces an IR intensity and thus makes the ν1 band observable, but it also produces a shift from ∼678 cm-1 for the free PF6- anion to 671 cm-1, a shift that should have been noticeable in the IR experiment, even taking into account the EMI band. Thus the ion-pair model creates too much anisotropy. Using implicit solvation, both the C-PCM artificial spectra are excellent in the sense that the shift of the ion pair is severely reduced, down to