Article pubs.acs.org/JCTC
Comparison of the Effective Fragment Potential Method with Symmetry-Adapted Perturbation Theory in the Calculation of Intermolecular Energies for Ionic Liquids Samuel Y. S. Tan* and Ekaterina I. Izgorodina* School of Chemistry, Monash University, 17 Rainforest Walk, Clayton, Victoria 3800, Australia S Supporting Information *
ABSTRACT: The effective fragment potential (EFP) method that decomposes the interaction energy as a sum of the five fundamental forceselectrostatic, exchange-repulsion, polarization, dispersion, and charge transferwas applied to a large test set of ionic liquid ion pairs and compared against the stateof-the-art method, Symmetry-Adapted Perturbation Theory (SAPT). The ion pairs include imidazolium and pyrrolidinium cations combined with anions that are routinely used in the field of ionic liquids. The aug-cc-pVDZ, aug-cc-pVTZ, and 6311++G(d,p) basis sets were used for EFP, while SAPT2+3/aug-cc-pVDZ provided the benchmark energies. Differences between the two methods were found to be large, and strongly dependent on the anion type. For the aug-cc-pVTZ basis set, which produced the least errors, average relative errors were between 2.3% and 18.4% for pyrrolidinium ion pairs and between 2.1% and 27.7% for imidazolium ion pairs for each individual energetic component (excluding charge transfer), as well as the total interaction energy. Charge transfer gave the largest relative errors: 56% and 63% on average for pyrrolidinium- and imidazolium-based ion pairs, respectively. Scaling of the EFP components against SAPT2+3 showed improvement for polarization (induction) and dispersion terms, thus indicating potential for the development of cost-effective alternatives for intermolecular induction and dispersion potentials for ionic liquids. constructed, such as an AMOEBA (a multipolar force field) force field for imidazolium-based ionic liquids,11 and the force field developed by Choi et al. specifically for 1-butyl-3methylimidazolium tetrafluoroborates.12 The latter parametrized intermolecular interaction terms to reproduce SAPT energies. This resulted in a better force field that performed well in the prediction of thermodynamic and transport properties of the ionic liquid when compared to experiment. The Choi force field is expected to be transferable due to the parametrization from first principles. Nevertheless, while classical molecular simulations are suitable for large systems, methods based on the first-principles of theory provide the benefit of a stronger understanding and the ability to systematically improve errors. To this end, this study will consider two robust methods that decompose the interaction energy into fundamental componentssymmetryadapted perturbation theory and the effective fragment potential methodwith the goal of seeing how they compare for ionic liquids, and to gain insight into how the effective fragment potential method performs for charged species. Symmetry-adapted perturbation theory (SAPT)13−15 is the state-of-the-art method for calculating the five fundamental components of intermolecular interactions: electrostatics, exchange-repulsion, induction, dispersion, and charge transfer.
1. INTRODUCTION Intermolecular interactions have a significant effect on the physical and chemical properties of semi-Coulombic condensed systems such as ionic liquids, because of the importance of noncovalent interactions between molecules. The calculation of the interaction energy requires accounting for nonspecific interactions such as the Coulomb interaction that dictates much of the intermolecular dynamics of ions, as well as other specific interactions such as hydrogen bonding, π−π stacking, and van der Waals interactions. 1−3 The noncovalent interactions in ionic liquids are dominated not only by electrostatics (Coulomb) but also by dispersion and induction (also known as polarization).4,5 The complex interplay of all these interactions affects their thermodynamic and transport properties in a nonlinear fashion5 and therefore, characterization of the intermolecular dynamics of ionic liquids is a challenging task.3 This challenge of predicting properties for ionic liquids has been addressed by several groups previously using classical force fields.6,7 The advantage of classical methods is that simulations can be run on large systems with relatively little computational effort. Induction was shown to be very important in the prediction of physical properties of ionic liquids, with polarizable force fields that explicitly include induction effects performing better for these complex systems.2,8−10 Furthermore, methods that offer an energy decomposition of the intermolecular interaction have been © XXXX American Chemical Society
Received: February 7, 2016
A
DOI: 10.1021/acs.jctc.6b00141 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
surfaces and interaction energy components. The average rootmean-square deviation (RMSD) between EFP and CCSD(T) was 0.49 kcal/mol with a range of 0.31−0.66 kcal/mol.48,49 For hydrogen-bonded and stacked DNA base pairs, Smith et al. found fair to excellent agreement between EFP and CCSD(T) results,50 and further emphasized the importance of a method that accurately treats the dispersion and polarization terms. While the electrostatic term dominates in hydrogen-bonded systems, its contribution to the binding energy is largely offset by exchange-repulsion, meaning that dispersion and polarization interactions play an important part in the stabilization of these systems. EFP is able to account for these interactions at a much lower cost, compared to MP2, CCSD(T), and SAPT, while being in reasonable agreement with these methods.50 Hands and Slipchenko have also pioneered a study validating the EFP method for polar heterogenic systems51 such as hydrogen-bonded tert-butanol−water mixtures. While the EFP method has been shown to compare well with other higher level correlated methods of ab initio theory for small clusters, much of its potential lies in the possibility of applying it to large clusters, which are typically treated with force fields and molecular dynamics. They found that while EFP accurately reproduced the experimental structure of pure water and water−tert-butanol mixtures, hydrogen-bonding patterns are very sensitive to the model potential used. Furthermore, the significance of the polarization energy was emphasized in the bulk compared to dimer systems where electrostatics and repulsion terms dominated. Flick et al. undertook a systematic study on the performance of EFP for the prediction of interaction energies in noncovalently bound complexes.52 The results were compared against a raft of semiempirical and correlated wave functionbased methods as well as methods of density functional theory. In their study, the two well-known databases, S22 and S66 described in the work of Hobza et al.,53,54 were used. These benchmark databases consist of a well-balanced series of intermolecular complexesfor example, complexes between small molecules such as water and ammonia, DNA base pairs, and amino acid pairs. These results showed that EFP is a reliable method to give a balanced description for different types of interactions, performing well for dispersion-dominated complexes and complexes with a mixture of interaction types. The accuracy of EFP, when compared with the benchmark CCSD(T)/CBS energies, approached that of second-order Møller−Plesset (MP2) perturbation theory method, which is widely used for studying intermolecular complexes. The EFP method reproduced the interaction energies well, with a mean unsigned error of 2.5 kJ mol−1 for the S66 set and 3.8 kJ mol−1 for the S22 set. This corresponds to a relative error of 11%− 12% in interaction energy. They found that the main sources of error in EFP came from the underestimation of the polarization and Coulomb components, which they attributed to insufficient treatment of short-range charge-penetration effects. This is partially offset by the underestimation in the exchangerepulsion term. For the S22 set, Flick et al. also found that the strongly interacting hydrogen bonding complexes had below average accuracy, with an underestimation of 10%−15% in the electrostatic energy.52 This finding was attributed to incomplete treatment of charge penetration. The electrostatic term for complexes with mixed interaction types was better treated since the magnitude of the interaction was smaller, and chargepenetration effects did not dominate. This is not observed in
The analysis of these components provides important insight into how they affect the structure and physicochemical properties of the chemical system in consideration.16,17 For example, SAPT was used to study dispersion18 and induction forces in small organic molecules and organic crystals,13,14,19 “weak, medium, and strong charge-transfer complexes” such as C2H4 for electron donors and F2 as acceptors, done at the CCSD(T)/CBS limit,20 and hydrogen bonding in water clusters,21 as well as the helium dimer potential,22 and π−π interactions in benzene.23,24 However, while accurate, this method is very expensive computationally, scaling to at least N5 when truncated at second order (where N is the number of basis functions). Within the second order, electrostatics and exchange-repulsion are taken into account in the first order of the interaction potential, whereas the dispersion and induction terms incorporate second-order terms. SAPT can also be made to account for intramolecular electron correlation effects within each interacting molecule. The general effective fragment potential (EFP2, referenced hereafter as simply EFP) method was developed by Gordon et al.25−29 as a computationally inexpensive method to model interaction energies and energetic components of intermolecular interactions. This method was originally created to model solvents (e.g., water),30−32 but has been later generalized for noncovalently bound complexes.33,34 It belongs to a class of fragmentation methods and is a purely ab initio-based method, without any empirical parameters, with each term being developed independently of the rest. Each term in the EFP method represents an individual fundamental component of interaction energyelectrostatics, exchange, polarization (i.e., induction), dispersion, and charge transfer, with the interaction energy being calculated as a sum of these terms. The EFP interaction energy is directly comparable with correlated wave function-based methods such as coupled cluster and perturbation theory. Due to its cost-effective formulation, EFP has found numerous applications in chemistry, such as studying solvent effects for which it was originally formulated.30−32,35−38 Other applications include solvent-induced shifts in the electronic spectra of uracil,39 spectroscopy of enzyme active sites,40,41 noncovalent π−π and hydrogen-bonding interactions in DNA strands using (nucleobase oligomers),34 hydrogen bonding,42 and other noncovalently bound systems.27,43,44 EFP has been applied to reliably study a broad range of intermolecular complexes. For example, EFP relative energies of hydrogenbonded complexes in (MeOH/H2O)n clusters, for n = 2, ..., 8 were in good agreement with MP2.32 In the case of many configurations of styrene dimers, exhibiting hydrogen bonding and π−π stacking interactions, EFP generally reproduced MP2 geometries and binding energies.45 In the work of Slipchenko et al.,46 the EFP method was found to perform very well for the π−π stacked benzene dimers in three different configurations. For the most challenging configuration, the parallel-displaced benzene dimer, the difference between EFP and CCSD(T)/ aug-cc-pVQZ was only 0.3 kcal/mol, which is a smaller error than typical MP2 at a much lower computational cost. In addition, the EFP method was shown to recover ∼70% of the charge penetration energy, arising due to the decrease of the classical electrostatic energy, as a result of electron density overlap.46 In further work examining π−π stacking,47 studying dimers of benzene, pyridine, and benzene−pyridine, the EFP method was found to be in good agreement with CCSD(T), MP2, and SAPT2 results for the calculation of potential energy B
DOI: 10.1021/acs.jctc.6b00141 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
A brief overview of the theory is given in the next section. The computational methodology is outlined in Section 3, followed by discussion of the results in Section 4.
the S66 dataset as the hydrogen-bonding complexes included in the set have weaker interactions, compared to those in S22. Ionic liquids generally have strong noncovalent interactions, and many ionic liquids have significant hydrogen bonding. In light of this, EFP might underestimate the strong interactions present in ionic liquids. Induction (polarization) was another term that was significantly underestimated (4−5 kcal/mol) for hydrogenbonded complexes.52 The relative energies for dispersiondominated complexes have even larger errors. Because of the small magnitudes of their polarization energies, the associated absolute errors are not large (largest 2 kcal/mol). This is attributed to the limitations of the multipole approximation for dispersion-dominated systems. The underestimation of the electrostatic and polarization terms discussed previously is offset in part by underestimation in the exchange-repulsion term, with typical relative errors of 10%−20%. This underestimation is hypothesized to arise from the neglect of correlation effects in EFP. Although charge transfer is negligible for dispersiondominated neutral complexes, it may play a more important role in strongly hydrogen-bonded charged complexes. In ionic liquids, net charge transfer has been demonstrated to occur regularly and is seen to be an important contribution to ion dynamics.3 The magnitude of its stabilization effect on the total interaction energy is challenging to estimate from the computational point of view. EFP is an attractive option for studying ionic liquids, since it decomposes the interaction energy into physically meaningful components. Energy decomposition provides further insight on the relationship between calculated energies and physicochemical properties.5 Understanding this relationship will tremendously aid in the prediction of properties and design of new ionic liquids. Furthermore, EFP is computationally costeffective, scaling as N2.52 To date, no systematic study has been done on the suitability of the EFP method for charged species such as ionic liquids. While this method was not originally designed for charged species, its computational efficiency renders it a good candidate. Based on the findings of Flick et al.,52 the three terms, Coulomb, exchange-repulsion, and polarization, are expected to show larger errors for ionic liquids. In addition, charge transfer might also be treated inaccurately due to the significant orbital overlap in ionic liquids, as is indicated by fractional charges that ionic liquid ions adopt in large-scale clusters.1,55−58 The chargetransfer term has been shrouded in controversy, and it is still unclear whether it should be included as a stabilizing energy term or whether an actual transfer of electronic charge is necessary for performing calculations involving accurate structure and dynamics.59−64 This is due to charge transfer being inherently tied with polarization, leading to varied formulations ranging from charge transfer being an artifact of the incompleteness of the basis set in SAPT to theories such as natural energy decomposition analysis,65,66 where it is a significant force in intermolecular attraction. Accurate treatment of these components is crucial in ionic liquids, where these are expected to constitute the bulk of the total interaction energy. This study identifies how well EFP performs for ionic liquids for the prediction of interaction energies in single ion pairs of ionic liquids. The test set consists of a large number of cation−anion combinations at various configurations. Imidazolium- and pyrrolidinium-based cations were coupled with routinely used anions, such as chloride and tetrafluoroborate.
2. THEORETICAL BACKGROUND 2.1. SAPT. SAPT was first used by London et al.67 to describe the intermolecular interaction operator as a multipole expansion. The theory has been further improved and refined and is the current benchmark for calculating the intermolecular interaction energy between two molecules. Jeziorski, Moszynski, and Szalewicz have presented a comprehensive description of the theory elsewhere.68 In the context of this work, “dimer” refers to an ion pair, while “monomer” refers to an individual ion. Note that, within the SAPT formulation, the intermolecular interaction energy, which is defined as the difference between the total energies of the dimer and constituent monomers, is calculated free of basis set superposition error. In order to obtain physically sound concepts of intermolecular forces such as the electrostatic, dispersion and induction, a nonsymmetric decomposition of the Hamiltonian is used. This means that electrons are no longer indistinguishable, and the corresponding zeroth-order wave function no longer obeys the Pauli exclusion principle. As a result, antisymmetrization is required, making the antisymmetrized wave function no longer an eigenfunction of the unperturbed sum of the Hamiltonian’s of constituent monomers, HA + HB, where A and B are monomers. To circumvent the issue, the symmetry-adapted perturbation procedure is applied to keep the HA + HB sum as the unperturbed operator while still utilizing the antisymmetrized wave function. The SAPT method has the Hamiltonian partitioned as H = FA + FB + WA + WB + V (1) where FA and FB are the Fock operators for monomers A and B, respectively. Similarly, WA and WB are the differences between the exact Coulomb operator and the Fock operator for each monomer, whereas V contains all the intermolecular terms. SAPT perturbs all of WA, WB, and V through various orders to calculate the individual energy terms. The different energy components are grouped to produce five fundamental forces:The superscripts in parentheses denote the perturbation order
of V and W = WA + WB, respectively. In this work, the 2+3 truncation was used in the SAPT expansion.17 In the equations above, blue singly underlined terms refer to second order, whereas red doubly underlined terms represent third order . The first order has only electrostatics and exchange terms, while induction and dispersion occur in the second order. Also present in the second order is quenching via exchangeC
DOI: 10.1021/acs.jctc.6b00141 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
water molecule has five occupied orbitals and 60 virtual orbitals in the 6-31++G(3df,2p) basis set. Thus, the calculation for the CT term is usually 20−30 times slower than that for the other terms.75 In order to trim the expense of CT calculations, quasiatomic minimal-basis-set orbitals77 are used that include the valence virtual orbitals. The latter ensures recovery of the most important CT interactions in the virtual space. Dispersion is treated by the sum of two terms,
repulsion in the intramolecular contributions to electrostatics and exchange. SAPT2+ further includes intramolecular electron correlation terms pertaining to dispersion. The third order, SAPT2+3, consists of additional terms for dispersion, as well as quenching of induction and dispersion by third-order exchange. In perturbation theory, the induction energy can be separated into two categories: those involving excitations from the occupied orbitals of a molecule to virtual orbitals of the same molecule, and excitations from the occupied orbitals of a molecule to the virtual orbitals of another molecule.69 The latter is known as the charge-transfer (CT) energy. To calculate this interaction, the SAPT induction energy from the monomer basis set, where no charge transfer is permitted, is subtracted from the SAPT induction energy from the dimer basis set. Furthermore, note that induction and dispersion include exchange components due to the quenching of forces, as a result of proximity of the interaction species and non-negligible orbital overlap.68,70−74 2.2. EFP. The effective fragment potential method is an ab initio-based potential method that models the intermolecular interactions of noncovalently bound systems using a costeffective formulation.26−28,33,34 In the EFP method, the system is broken into fragments. Typically each of the interacting molecules is a fragment. In this study, since only ion pairs are considered, individual ions are treated as single fragments. Each fragment is treated separately at the Hartree−Fock level of theory in order to generate potentials. The cation and anion potentials then are allowed to interact and the interaction energy is decomposed into individual components. The EFP method partitions the interaction energy in the following way: Etotal interaction = EElst + EPol + EDisp + ERepl + ECT
E Disp =
C6 6
R
+
C8 R8
(8)
The first term is the induced dipole−induced dipole interaction. In the EFP method, the coefficients for this C6 term are calculated through the interactions between pairs of localized molecular orbitals of each ion using the timedependent Hartree−Fock method, with the C8 coefficients being approximated as one-third of those of C6.78 This expression is corrected for short-range charge penetration effects through a distance-dependent damping function. Note that this dispersion term was formulated by comparison with SAPT dispersion as the benchmark. The exchange-repulsion is also calculated using a static localized molecular orbital basis by expanding the intermolecular overlap integral, with truncation at the quadratic term for exchange-repulsion. Eijexch = − 4
⎛ Sij 2 −2 ln|Sij| − 2Sij⎜⎜ ∑ FikASkj + R ij π ⎝k∈A
⎛ Z 1 − 2Sij 2⎜⎜ ∑ I + 2 ∑ + R R kj k∈A ⎝ I ∈ A Ij
(7)
Here, EElst is the electrostatic component, EPol the induction (polarization) component, EDisp the dispersion component, ERepl the exchange-repulsion component, and ECT the chargetransfer component. Coulomb, induction, and dispersion are considered long-range interactions. They decay as R−n, with n = 1 for Coulomb, n = 2−4 for induction, and n = 6 for dispersion. The short-range interactions that decay exponentially are exchange-repulsion and charge transfer. In the EFP formulation, the Coulomb interaction uses Stone’s distributed multipolar analysis16 truncated at the octopole term. To correct for charge penetration effects arising from orbital overlap, a damping function is employed. This overlap-based screening has exponential dependence on the separation.46 Induction, also known as polarization in the EFP method, is the effect of inducing a dipole moment in a molecule by the electric field of another. This term is treated with dipole polarizability tensors located at the centroids of localized bond and lone pair orbitals of the molecules.75 To prevent “polarization collapse” at short intermolecular distances, EFP employs Gaussian-type damping, because of its mathematical simplicity and independence of the choice of Coulomb damping.76 Polarization in EFP is analogous to induction in SAPT with one exception. The SAPT induction term also contains the charge-transfer energy, whereas it is calculated separately in the EFP method and is defined as the interaction between the occupied orbitals on one EFP fragment with the virtual orbitals of another fragment. For charge transfer, the EFP method uses a second-order perturbation at the HF level of theory.75 As these calculations involve the virtual orbitals, it becomes slower as the number of basis functions increases. For example, the
∑ J∈B
ZJ R iJ
⎞
∑ FjlBSil − 2Tij⎟⎟ ⎠
l∈B
+ 2∑ l∈B
⎞ 1 1 ⎟ − R il R ij ⎟⎠ (9)
where A,B are the effective fragments, i, j, k and l are the LMOs, and I, J are the nuclei. S refers to the intermolecular overlap integral, and T to the kinetic energy integral. The Fock matrix element is represented by F.34 It is expected that higher-order correlation effects are not well accounted for in second-order exchange-repulsion. While the computational costs for each component varies, depending on system size and complexity, the most expensive interactions to calculate by means of EFP are generally the exchange-repulsion and charge-transfer interactions. These two components might be more than five times as computationally demanding than the other three components, which are of roughly the same cost, relative to each other. Originally, the exchange-repulsion and charge transfer were designed with optimizations for neutral molecules; therefore, it is suggested that these terms might not perform as well for charged species such as ionic liquids. These interactions will be stronger, because of greater orbital overlap among ions. While stronger interaction energies in ionic liquid ion pairs might not result in higher relative errors, absolute errors would be expected to be larger. In comparing EFP and SAPT, Table 1 describes which terms from each method will be compared against each other. Note that, because of the difference in the definitions of charge transfer, in the EFP and SAPT approaches, the EFP polarization term was directly compared to the SAPT induction without inclusion of the CT energy. D
DOI: 10.1021/acs.jctc.6b00141 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation Table 1. Energetic Components from SAPT and EFP Compared with Each Other, and Abbreviation of These Components Used in the Text component
SAPT name
EFP name
abbreviation
electrostatics exchange-repulsion induction/polarization dispersion charge transfer
Eelectrostatics Eexchange Einduction Edispersion Echarge‑transfer
EElst ERepl EPol EDisp ECT
EElst EExch EInd EDisp ECT
3. METHODOLOGY 3.0.1. Chemical Systems Studied. The chemical systems studied were single ion pairs of ionic liquids. Routinely used anions such as tetrafluoroborate (BF4−), bromide (Br−), chloride (Cl−), dicyanamide (Dca−), mesylate (Mes−), tosylate (Tos−), hexafluorophosphate (PF−6 ), and bis{(trifluoromethyl)sulfonyl}amide (NTf−2 ) were used in this study. These anions were combined with N-alkyl-N′-pyrrolidinium (denoted here as Cnmpyr) and 1-methyl-3-alkyl-imidazolium (denoted here as Cnmim) cations with varying alkyl chain from methyl, ethyl, propyl to butyl. The names and abbreviations of the different cations and anions are tabulated in Table 2. Furthermore, in the text, halides are abbreviated as “Hal”, whereas the rest of the anions are referred to as typical ionic liquid anions (TILAs).
Figure 1. Different configurations of [C3mim][Br].
Table 2. List of Cations and Anions cation 1-methyl-3-methylimidazolium 1-methyl-3-ethylimidazolium 1-methyl-3-propylimidazolium 1-methyl-3-butylimidazolium N,N-dimethylpyrrolidinium N-ethyl-N′-methylpyrrolidinium N-propyl-N′-methylpyrrolidinium N-butyl-N′-methylpyrrolidinium
abbreviation
anion
abbreviation
C1mim+
tetrafluoroborate
BF−4
C2mim+
bromide
Br−
C3mim+
chloride
Cl−
C4mim+
dicyanamide
Dca−
C1mpyr+
mesylate
Mes−
C2mpyr+
NTf−2
C3mpyr+
bis{(trifluoromethyl) sulfonyl}amide hexafluorophosphate
C4mpyr+
tosylate
Tos−
PF−6
For each cation−anion combination, different configurations of these ion pairs were incorporated. These configurations differ by how the anion interacts with the cation. In the imidazolium-based cation, the anion can interact with the cation above and below the imidazolium ring; these configurations are referred to as p1 and p4, respectively. When the anion interacts with the cation in the C2−H bond plane, this configuration is referred to as in-plane interactions and is denoted hereafter in the text as p2 and p3. The different configurations for [C3mim][Br] are presented in Figure 1 as an example. For the case of pyrrolidinium-based ion pairs, the configurations studied are different as the anion has a tendency to interact with the nitrogen center of the cation from three energetically favorable positions denoted here as p1, p2, and p3.5 The NTf−2 anion has multiple interaction sites such as the central nitrogen and the oxygens on the sulfonyl groups, as shown in Figure 2. Therefore, there are more ion-pair configurations corresponding to the anion interacting through
Figure 2. Different configurations of [C2mpyr][NTf2].
different centers, compared to halides. In the case of the NTf−2 anion, there are six possible configurations (for more details, see Figure 2). The [Cnmim] X series of ion pairs, where X represents chloride or bromide, were optimized at the MP2/aug-cc-pVDZ E
DOI: 10.1021/acs.jctc.6b00141 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
interaction energy is 7.3 kJ mol−1, with a standard deviation of 1.2 kJ mol−1. Excluding charge transfer, the differences between the two basis sets for each energetic component (that is, electrostatics, exchange, induction, and dispersion) range between −7.0 kJ mol−1 and 4.8 kJ mol−1. Considering the low standard deviations (0.5−1.5 kJ mol−1), this indicates that these two basis sets differ consistently for each energetic component. Note that the aug-cc-pVTZ basis set requires tremendous amounts of CPU timein most cases, more than double the amount required for the aug-cc-pVDZ basis set. For example, for [C4mim][Cl], aug-cc-pVTZ required 329 CPU h and 22 GB of memory, compared to 26 h and 4 GB for aug-cc-pVDZ. Taking the computational expense into account, the aug-cc-pVDZ is the largest basis set possible for many of the bulky ionic liquid ion pairs studied, such as [C4mpyr][NTf2]. SAPT2+3 calculations on the intermolecular complexes in the S22 and S66 datasets53,54,87 were also performed using augcc-pVDZ.
level, while for the other TILAs, MP2/6-31+G(d,p) was used. For the [Cnmpyr] series of ion pairs, geometry optimization was performed at the B3LYP/6-31+G(d) level. The geometries of these configurations have previously been published by our group.5,79 3.0.2. Software. The PSI4 quantum chemistry package was used for the SAPT2+3 calculations.17 All SAPT calculations were performed using the aug-cc-pVDZ basis set, unless stated otherwise.5 The GAMESS-US software package was used to perform the EFP calculations.80,81 Three basis sets were used for EFP calculations: aug-cc-pVDZ, aug-cc-pVTZ, and the Pople basis set 6-311++G**. Note that, for the bromide anion, since it is not included in the 6-311++G** basis set, 6-311G** basis functions were used instead. Three basis sets were employed to determine how consistently the method performs for this series of basis sets, and observe the effect of the basis set on the EFP performance for ionic liquids. The Avogadro program was used to visualize molecules.82 For data analysis, the R statistical programming language was used,83 with the ggplot2 plotting package and the data.table package.84,85 3.0.3. Basis Sets. While three basis sets were used to give an indication of basis set dependency for the EFP method, only the aug-cc-pDVZ basis set was used for SAPT2+3. This is consistent with the results presented by Parker et al., where, for SAPT2+3, the basis set with the lowest mean absolute error was aug-cc-pVDZ. (See Tables III and IV in ref 15, pp 094106− 094108)15 This holds true across all the different subsets, which are categorized based on the dominating interactions, such as hydrogen bonding, dispersion, and mixed. To further validate this choice, many test calculations were performed using aug-cc-pVTZ, indicating that aug-cc-pVDZ gave satisfactory accuracy. The calculations using aug-cc-pVTZ were run on select representative systems, largely with halide anions, namely, [Cnmim]X, where X = Cl and Br. The rest of the ion pairs for which SAPT2+3/(aug-cc-pVTZ) calculations were run are [Cnmpyr]X, where n = 1, 3 and X = BF−4 and Cl−. Halidesin particular, chlorideswere selected due to smaller system sizes (monatomic anions). These are strongly bound to the cation, representing challenging systems from the theoretical point of view.86 The differences between these two basis sets are reported in the Supporting Information. The consistency between the two basis sets (E(aug-cc-pVTZ) − E(aug-cc-pVDZ)) is immediately apparent: with the triple-ζ basis set, electrostatics is always weaker compared to the double-ζ basis set, by 2.7 kJ mol−1 on average (3.0 kJ mol−1 for imidazolium systems and 1.9 kJ mol−1 for pyrrolidinium systems). Exchange follows the same trend, with the larger basis set giving a larger exchange energy of 4.9 kJ mol−1, on average. It is not surprising that the larger aug-ccpVTZ basis set leads to greater recovery of the dispersion energy, by 5.6 kJ mol−1 on average. There is no such pattern observed in the induction energy, with both basis sets giving excellent agreement to each other at
