Cohesive Properties of Ionic Liquids Calculated from First Principles

Subscriber access provided by Nottingham Trent University

Thermodynamics

Cohesive Properties of Ionic Liquids Calculated from First Principles Ctirad #ervinka, Martin Klajmon, and Vojt#ch Štejfa J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00625 • Publication Date (Web): 22 Aug 2019 Downloaded from pubs.acs.org on August 29, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Cohesive Properties of Ionic Liquids Calculated from First Principles Ctirad Červinka*, Martin Klajmon, Vojtěch Štejfa †

Department of Physical Chemistry, University of Chemistry and Technology Prague, Technická

5, CZ-166 28 Prague 6, Czech Republic

*Corresponding author: [email protected]

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 46

ABSTRACT. Low volatility of ionic liquids (ILs), being one of their most valuable properties, is also the principle factor making reliable measurements of vapor pressures and vaporization (or sublimation) enthalpies of ILs extremely difficult. Alternatively, vaporization enthalpies at the temperature of the triple point can be obtained from the enthalpies of sublimation and fusion. While the latter can be obtained calorimetrically with a fair accuracy, the former is in principle accessible through ab initio computations. This work assesses the performance of the first-principles calculations of sublimation properties of ILs. Namely, 3 compounds, coupling the 1-ethyl-3-methylimidazolium cation

[emIm]

with

either

tetrafluoroborate

[BF4],

hexafluorophosphate

[PF6],

or

bis(trifluoromethylsulfonyl)imide [NTf2] anions were selected for a case study. A computational methodology, originally developed for molecular crystals, is adopted for crystals of ILs. It exploits periodic density functional theory (DFT) calculations of the unit-cell geometries and quasiharmonic phonons and many-body expansion schemes for ab initio refinements of the lattice energies of crystalline ILs. The vapor phase is treated as the ideal gas whose properties are obtained combining the rigid rotor – harmonic oscillator model with corrections from the one-dimensional hindered rotors, and molecular-dynamics simulations capturing the contributions from the interionic interaction modes. Although the given computational approach enables to reach the chemical accuracy (4 kJ mol−1) of calculated sublimation enthalpies of simple molecular crystals, reaching the same level of accuracy for ionic liquids proves challenging as crystals of ionic liquids are bound appreciably stronger than common molecular crystals – the underlying cohesive energies of solid ionic liquids is up to one order of magnitude larger. Still, combination of the mentioned computational and experimental frameworks results in a novel promising scheme that is expected to generate reliable and accurate temperature-dependent data on sublimation (and vaporization) of ILs.

2


Page 3 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION Ionic liquids possess a vast potential for numerous technologies (gas capture or separation, smart electrolytes, polymer or nanoparticles exfoliation media, etc.) thanks to their unique properties such as low volatility, large electrochemical window or a boundless structural variability.1-2 Broader exploitation of their beneficial characteristics is impeded by their cost, limited availability of the physico-chemical data or an insufficient understanding of ILs-related phenomena. Concomitantly, literature on ILs is growing unprecedentedly, although the quality of the published data is sometimes at least questionable.3 Known for a century and considered absolutely nonvolatile for decades, ILs were proved to possess a nonzero saturated vapor pressure only a decade ago.4-5 Still, vapor pressures of typical ILs are orders of magnitude lower when compared to common molecular solvents. Substitution of conventional volatile solvents in large-scale processes with suitable ILs would decrease evaporation of solvents into the atmosphere. On the other hand, the extremely low volatility is a principal obstacle impeding reliable and reproducible measurements of their vapor pressures, sublimation enthalpy (ΔsubH) or vaporization enthalpy (ΔvapH).6 The most common computational tool to study thermodynamic properties of ILs are currently the molecular dynamics (MD) simulations.7-8 Accuracy of the MD-based results strongly depends on the appropriateness of the force field model and on the quality of the phase space sampling. The latter factor along with the well-known slow internal dynamics of ILs causes need for generating long trajectories (tens of nanoseconds) to obtain statistically sound thermodynamic properties for the liquid phase of ILs.9 Although the force-field parameters are commonly derived from firstprinciples calculations which significantly improves the predictive power of MD simulations,10 it still remains challenging to capture some many-body effects such as polarizability or charge transfer accurately within the MD.11 There are numerous literature studies focusing on calculations

3



of phase change enthalpies and related properties for ILs,12-17 with results often exhibiting inacceptable scatter and interlaboratory disagreement. A typical discussion in this filed concerns e. g. the existence of a trend shift or prevailing linearity of ΔvapH within the homologous series of alkylimidazolium based ILs.17-19 Overall, an uncertainty of such thermodynamic data for ILs exceeding 10 % should be expected.17 Published studies that use static quantum-chemical tools to characterize ILs focus on the optimum ionic geometry and spectral characteristics, refining atomic parameters for classical force-field models,10,

20-21

or ionic interactions.22-26 There were attempts aiming at modeling the cohesive

energy of liquid ILs with embedded fragment-based quantum-chemical calculations for ionic clusters of various size.27-31 Unlike for crystals, the finite-cluster approximations for liquids always face issues related to the limited cluster size or insufficient phase space sampling. To our knowledge, no study using fragment-based or periodic quantum-chemical approaches for crystalline ILs exists. These approaches underwent a rapid development and emerged as robust and fairly accurate tools in the field of molecular crystals. The experimental crystal structure is the only a priori requirement for such calculations. Note that the Cambridge Structural Database contains millions of entries on crystal structures, including numerous crystalline ILs.32 A rigorous methodology, capable of accurate determinations of the volatility of ILs, which would be based on more reliable approaches (combining experiments and first-principles calculations) than the direct measurements of vapor pressures of ILs or (crude) semi-empiric computational approaches would certainly found use in the future design of chemical processes and selection of suitable ILs. Especially if there are millions of relevant ILs systems, while the thermodynamic data exist for no more than hundreds of ILs.33 Low volatility of ILs arises due to the strong electrostatic cohesion of the ions.34 As most ILs are derived from organic or inorganic molecules, a condensed phase of an IL can be regarded as a system of charged molecules. In other words, solid ILs represent an intermediate step between 4


Page 4 of 46

Page 5 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simple inorganic ionic materials and conventional molecular materials35 which enables to model the ILs systems with similar computational methodologies as in case of molecular systems.10 A special attention must be paid to the treatment of the electrostatic interactions such as long-range charge/dipole interactions, induction effects due to the polarizability, or charge transfer between the ions. Still, dispersion interactions and site-specific hydrogen or halogen bonds play a considerable role in the subtle interplay of all molecular interactions in systems containing ILs.10, 36

The above-used term cohesive properties unites the cohesive/lattice energy of a crystal, its ΔsubH, and ΔvapH, which is related to ΔsubH through the enthalpy of fusion (ΔfusH) at the temperature of the triple point. While ΔfusH of ILs is measurable with a reasonable experimental uncertainty,37 ΔsubH of a molecular crystal can be calculated with a fair (chemical or better) accuracy from first principles from the cohesive energy and phonon-based properties.38-39 High level ab initio fragment-based combinations52,

techniques,40-47 57-59

periodic

quantum-chemical

calculations,48-56

or

their

along with the quasi-harmonic approximation60 capture these finite-

temperature properties satisfactorily for molecular crystals.50, 61-68 It has become a relatively straightforward task to calculate the thermodynamic properties of the ideal gas using quantum-chemical (mostly DFT) calculations of molecular parameters and statistical-thermodynamic models.69-70 Uncertainty of the resulting heat capacities and entropies amounts to a few units of percent for small and relatively rigid molecules.71-72 The complexity of this task, however, significantly increases with molecular flexibility and number of conformers of the given molecule. This approach has already been employed to calculate the thermodynamic properties of the ideal-gaseous phase of several imidazolium-based ILs, treating the mutual cationanion degrees of freedom as hindered rotations or harmonic oscillations.73 Although the authors report that their computed entropies are reasonably consistent with the experimental data, weaknesses of such an approach are the high number of conformations of the whole ion pair, and 5



possibly, an ambiguous definition of the interionic hindered rotors. An alternative to this fully statistical-thermodynamic approach consists in using the statistical-thermodynamic model only for isolated ions while their interaction and resulting contributions to the thermodynamic properties are extracted from molecular dynamics simulations. Three archetypal ILs were selected for investigation in this work, namely [emIm][BF4], [emIm][PF6], and [emIm][NTf2], since various independent experimental data on their thermophysical properties, relevant in the context of this work, exist in the literature. This work demonstrates that an approach exploiting high-accuracy first-principles calculations of ΔsubH with critically assessed experimental ΔfusH represents a novel promising scheme capable of generating reliable and accurate temperature-dependent data on liquid – vapor equilibrium of ILs. Although this work is exclusively computational, it manifests the novelty of the presented approach consisting in a unique symbiosis of computational chemistry and experiment, benefitting from the strong points of both − reliable calculations of ΔsubH at arbitrary finite temperatures and accurate calorimetric determination of ΔfusH. Furthermore, it enables to avoid the difficulties related to direct determination of ΔvapH for ILs − prohibitive cost of ab initio molecular dynamics for liquids, dependence on the quality of the force field for classical molecular dynamics, or difficult experimental setup and a substantial experimental uncertainty. Uncertainty of vaporization properties obtained by the presented approach is expected not to exceed the threshold valid for corresponding experimental determinations. Moreover, there is no temperature or pressure limitation for the range of applicability of such a combined approach, unlike for purely experimental determinations of vapor pressures of ILs that are feasible only at sufficiently high temperatures.

6


Page 6 of 46

Page 7 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2. COMPUTATIONAL METHODS 2.1

Crystalline phase

All computations for crystalline ILs were initiated from the experimental crystal structures culled from the CSD database.32 Corresponding structures are listed in Table 1 and illustrated in Figure 1. Electronic structure of the crystals was treated using the periodic PBE-D3(BJ) calculations74-75, the projector augmented wave (PAW) technique,76 coupled with standard PAW potentials,77 400 eV cut-off for the plane-wave energy, and with a Monkhorst-Pack type Γ-point-centered mesh for k-space sampling,78 as implemented in VASP, version 5.4.1.79 Number of k-points in directions of the inverse unit cell vectors were set to roughly 30/a k-points where a stands for the length of the unit cell vector in the corresponding direction.62 Using this computational setup, unit-cell geometries were optimized with only the space group of symmetry constrained. Following the quasi-harmonic approximation,60 dependence of the cohesive energy on unit-cell volume Ecoh(V) was extracted from performing additional 15 fixed-volume optimizations, with unit-cell vectors gradually scaled by factors ranging to 0.95−1.09. TABLE 1 Experimental crystal structures used as input for computational treatment of the crystal phases of ILs, along with their computationally optimized counterparts obtained minimizing either the static electronic cohesive energy (at 0 K) or the total quasi-harmonic (QHA) Helmholtz energy. Ionic liquid

CSD refcode

Space group

Za

Data set

a, Å

b, Å

[emIm][BF4]

LAZRIO01

P21/c

4

Experimentb

8.653

9.285

13.217 121.36

Staticc

8.68

9.18

13.43

120.9

QHAd

8.58

9.00

13.29

120.9

Experimente

8.627

9.035

13.469 101.92

Staticc

8.66

8.89

13.60

101.3

QHAd

8.58

8.74

13.46

101.0

[emIm][PF6]

HAYBUE02

P21/c

4

7


c, Å

β, deg.


[emIm][NTf2]

RENSEJ01

Pna21

8

Page 8 of 46

Experimentf

27.713 7.0343 15.769

90.00

Staticc

27.56

6.90

16.25

90.0

QHAd

27.39

6.85

16.11

90.0

Number of ionic pairs forming a unit cell. Structure determined at 100 K in ref. 80 c Static unit-cell parameters (based on PBE-D3(BJ)/PAW electronic energy only) at 0 K. d Quasi-harmonic unit-cell parameters (based on PBE-D3(BJ)/PAW geometries and phonons and MP2C-F12/avdz + pB3LYP refinement of cohesive energies) at the temperatures of the experimental structure determination. e Structure determined at 173 K in ref. 81 f Structure determined at 100 K in ref. 82

a

b

Phonon properties were obtained using the finite-displacement supercell approach,83 replicating the optimized unit cells to supercells exceeding 10 Å size in all three dimensions. Energies related to all symmetry-inequivalent perturbed supercells (with atom displaced from the optimum geometry) were computed to evaluate the force constants for all fundamental vibrational modes which were subsequently used to construct the density of phonon states, as implemented in the Phonopy code.84 Phonons and the vibrational Helmholtz energy Avib were computed within the harmonic approximation for five unit-cell volumes around the minimum of the Ecoh(V) curve. Avib data points, obtained in this way, were then fitted by a quadratic function for individual temperatures in the range 0 K to 500 K. At this point, it is necessary to cover the whole expected volume range so that the quadratic Avib(V) function is not extrapolated to larger volumes. When extrapolations are needed, a linear fit of Avib(V) represents a safer alternative.57,

62

The total

temperature and volume-dependent Helmholtz energy of the crystal Acr was obtained as: Acr(T, V) = Avib(T, V) + Ecoh(V),

(1)

which is finally fitted by the Murnaghan equation of state at the individual temperatures.85 To allow for a sufficient flexibility of the fitting functions, expansion and compression branches of the Acr(V), being smoothly joined at the minimum point, were fitted independently. All thermodynamic properties of interest were then derived from Acr from the essential thermodynamic relations. Since it holds p = − δA / δV for the total pressure, equilibrium molar volume of the crystals at zero 8


Page 9 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pressure is accessible through searching for the minima of the Acr(V) curves at the individual temperatures. Furthermore, the isobaric capacity can be derived as: H   2G Cp   G  TS   T 2 , T T T

(2)

which becomes Cp = − T δ2A / δT 2 assuming the pressure is zero. To obtain a deeper understanding, the Ecoh data points obtained for all unit-cell volumes were further refined by performing single-point electronic energy calculations of the unit cell with various higher levels of theory, using the PBE-D3(BJ)/PAW(400 eV)-optimized unit-cell geometries.57,

59

For this purpose, either plane-wave based calculations with the hard PAW

potentials at the PBE-D3(BJ)/PAW(1000 eV) and B3LYP-D3(BJ)/PAW(1000 eV) levels, or calculations exploiting the hybrid many-body interaction (HMBI) model, developed for molecular crystals,38, 40, 46 were performed. The HMBI model selects a reference molecule in the crystal lattice for both the cation and anion and determines the neighboring molecules along with the contact distances to the reference one, taking the symmetry of the lattice into account.86 In this work, interactions of proximate molecular (ionic) pairs within a cut-off distance 10 Å were treated by the counterpoise-corrected explicitly correlated87-88 dispersion-coupled89 second-order perturbative method MP2C-F12 with the aug-cc-pVDZ basis set and the fixed amplitude 3C ansatz, as implemented in Molpro, version 2012.1.90 This level of theory has recently been appointed as the bronze standard method for ab initio calculations of non-covalent interaction energies.91 Distant interactions and many-body effects were treated at the cheaper periodic Hartree-Fock (pHF) level or periodic B3LYP-D3(BJ,ATM) level using the CRYSTAL code.92 Atomic basis set pob-TZVP93 was used here for all calculations in CRYSTAL. This basis set is reported to be refined for reliable periodic DFT and pHF calculations of electronic structure of solids, including ionic materials. 2.2

Vapor phase

9



Vapor phase of the studied ILs was assumed to be an ideal gas of neutral ionic pairs which is an appropriate approximation due to the low vapor pressures of ILs. First, thermodynamic properties of ensembles of bare cations and anions were computed. DFT calculations of optimum molecular geometry and vibrational frequencies were performed using the PBE-D3(BJ) functional and the split-valence 6-311+G(d,p) basis set. This particular functional was selected for the sake of consistency with the computations for crystals. For rigid [BF4] and [PF6] ions, the rigid rotor – harmonic oscillator (RRHO) model was used to calculate the intramolecular (vibrational) g0 g0 contributions to the heat capacity of the ideal gas of ions ( Cintra,c and Cintra,a ). No frequency scaling

was applied here, again due to consistency with the phonons of crystals. Thermodynamic properties of the flexible ions [emIm] and [NTf2], both reported to form more conformations, were modeled using the conformation mixing model, assuming that the pure ion is an equilibrium mixture of individual conformers. Electronic energies of the individual conformers were calculated at the PBE-D3(BJ)/6-311+G(d,p) level of theory. CCSD(T)/CBS refinements of the conformer energies were performed for the heat capacity data set which was used for correlations of experimental data. Population of the conformers was derived from their temperature-dependent molar Gibbs energies (including the electronic, zero-point vibrational, entropic and thermal contributions). Strongly anharmonic internal degrees of freedom, corresponding to hindered rotations of the methyl groups in [emIm] and of the trifluoromethyl groups in [NTf2], were treated using the one-dimensional hindered rotor (1D-HR) model.70, 94 Barriers to these internal three-fold rotations were calculated at the PBE-D3(BJ)/6-311+G(d,p) level of theory and the corresponding reduced moments of inertia were obtained from the optimized ionic geometry and using the formula for an asymmetric uncoupled rotating top.95 g0 To capture the interionic contributions to the heat capacity Cinter , which is expected to be strongly

anharmonic, classical molecular dynamics simulations of the ionic pair and of both isolated ions 10


Page 10 of 46

Page 11 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


were performed in the LAMMPS software96 (version 31 Oct 2015) using an all atom nonpolarizable OPLS force-field developed by Canongia Lopes and Pádua (hereafter CL&P forcefield).10 Namely, force field parameters were adopted from ref.97 for [emIm] and [BF4], from ref.98 for [NTf2] and from ref.99 for [PF6]. Nosé-Hoover thermostat was used to control the temperature of the NVT ensembles. After a 100 ns equilibration run, a trajectory covering 1000 ns with a time step of 1 fs and a sampling frequency of 10,000 fs (10 ps) was generated for each temperature in the range 20−500 K. Temperature steps of 10, 20 and 25 K were used to cover the ranges 20−100, 100−200 and 200−500 K, respectively. For the ion pairs and the isolated ions, the average total energy Etot and its statistical variance σ2(Etot), representing fluctuations of Etot which give birth to the heat capacity, were evaluated at each temperature. Subsequently, the contribution to the heat capacity of the vapor due to the cation-anion interaction was calculated at the individual temperatures from σ2(Etot) of the ion pair, the cation, and the anion, using the relation: g0 inter

C



pair cat ani )   2 ( Etot )   2 ( Etot )  2 ( Etot

RT 2

.

(3)

g0 The obtained Cinter data in the considered temperature range were smoothed by polynomial

functions, as shown further in Section 4. Finally, the total isobaric heat capacity of the ionic pair C pg0,pair was obtained as: g0 g0 g0 C pg0,pair  Cintra,c  Cintra,a  Cinter  4R .

(4)

Note that the 4R term arises as a sum of transforming isochoric CVg0 to isobaric C pg0 (a R term) and the contributions of the overall translation and rotation (both 3/2 R terms) which need to be included only once for the whole ion pair. As this approach is not capable of capturing the entropy of the ion pair directly, the following scheme was adopted for calculations of entropy of gaseous ILs. The most stable conformation of the ionic pair was searched for using the simulated annealing procedure using the same CL&P 11



force filed. A random initial configuration of the ionic pair at 500 K was gradually cooled down to 0.001 K over a 1 μs trajectory, giving a cooling rate of 0.5 K ns−1. The obtained annealed ionic pair geometry at 0.001 K was subsequently re-optimized at the PBE-D3(BJ)/PAW level in VASP and its harmonic vibrational frequencies were computed using the finite difference method. The key assumptions here are that the ionic pair still dominantly exists in its most stable conformation (realizable as two mirror images) at a sufficiently low temperature, namely at 10 K, and that the contributions of the translation and rotation of the whole ionic pair are already converged to the 3/2 R values at this (sufficiently high) temperature. Then, the original RRHO model was used to g0 calculate the total entropy of the whole ionic pair S pair at 100 kPa, using the PBE-D3(BJ)/PAW

g0 frequencies. Finally, this S pair (10 K) value was temperature-adjusted using C pg0,pair obtained using

equation (4) to include the conformation flexibility also in the entropy. To calculate ΔsubH and psub, the electronic energy of the isolated ion pair in its PBE-D3(BJ) optimized geometry was further refined either at the MP2C-F12/aug-cc-pVDZ level or using various DFT/PAW computational setups to test the consistency of the electronic energies of vapor and crystal.

12


Page 12 of 46

Page 13 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FIGURE 1. Illustration of the crystal structures (upper row) and optimized annealed isolated ionic pair geometries (bottom row) for the studied ILs: left – [emIm][BF4], middle – [emIm][PF6], right – [emIm][NTf2].

3

Reference experimental data

3.1

Melting point temperatures and Enthalpies of Fusion

Table S1 summarizes literature studies on phase behavior of the discussed ILs. An adiabaticcalorimetric study of [emIm][NTf2] was performed by Paulechka et al.,100 who described four monotropic polymorphs all melting within a short temperature range of 15 K. Since the data from adiabatic calorimetry data are generally of lowest achievable uncertainty and the data in ref.100 is more detailed than the other calorimetric studies on [emIm][NTf2], we do not list and use the other works for benchmarking purposes. There are numerous studies reporting the temperatures of the melting point (Tfus) of [emIm][PF6] and [emIm][BF4], see Table S1. Although most of these use differential scanning calorimetry (DSC), very few information on ΔfusH is available. The data on melting points of ILs are generally more scattered than for non-ionic organic compounds, which is mainly attributed to the purity of the samples related to problematic purification and their hygroscopicity. Despite the known effect of water on the properties of ILs,101 many of the studies do not report the water contents and do not treat the compounds adequately. Some doubts remain about ΔfusH of [emIm][BF4], for which the literature values are in a large disagreement and the data on Tfus and water contents do not serve as a clue either. Although a high water content was reported by Van Valkenburg et al.,102 we give confidence to this more complete and thorough study. The weighted averages of Tfus and ΔfusH (for [emIm][NTf2] values only from ref.100) are used in the following calculations. 3.2

Condensed Phase Heat Capacities

13



Page 14 of 46

Table S2 lists literature studies of the isobaric heat capacities of the condensed phases of the studied ILs. Heat capacities of [emIm][NTf2] measured by Paulechka et al.100 using the adiabatic calorimetry are considered to be more reliable than any DSC measurements. Data by Ferreira et al. 103

extend the availability of C p to significantly higher temperatures, being in a good agreement

with the adiabatic data100 in the overlapping interval at the same time. For [emIm][BF4], only multiple DSC determinations of C p exist. The Tian-Calvet DSC results, generally possessing lower uncertainty than those from heat-flux DSC, are all in good mutual agreement and are accepted for the development of the reference data. The data by Waliszewski et al. 104 cover the longest temperature range, but seem to be offset by 1% compared to the other four Tian-Calvet DSC. Such systematic error is tolerable for a calorimeter of this type, but the data in ref.104 were adjusted for higher reliability of the subsequent correlation. Only one reliable source of C p was found for [emIm][PF6] in the work by Serra et al.,105 reporting C p of both liquid and crystalline phases. Existence of a second (crII) polymorph of [emIm][PF6]

(monotropically related to crI), which exothermically transforms to crI between 270 and 290 K, was reported together with it’s C p .105 3.3

Saturated vapor pressures

Among the three discussed ILs, [emIm][NTf2] is definitely the most frequently studied compound with more than 15 papers containing determination of its volatility listed in Table S3. Most of the available data can be, however, addressed to three research groups only: Verevkin et al. 112(University

106-

of Rostock), Santos et. al.12, 113-114 (University of Porto), and Licence et al.34, 115

(University of Nottingham). Still, available data can be sometimes burdened by the following deficiencies: i) the authors limit their results on ΔvapH, although vapor pressures could be determined and presented from the same setup; ii) the authors are so eager to obtain ΔvapH at 298 K that they sometimes omit to present ΔvapH relevant to the original temperature interval of the 14


Page 15 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


measurement. Furthermore, the ΔvapCp(T) = C pg − C pl term required for the correction to 298 K is hardly available for most ILs and the authors utilize temperature-independent ΔvapCp values obtained by too crude approximations; and iii) the authors tend to be overly optimistic considering uncertainty of their data. 3.4

SimCor treatment and resulting correlations

Selected literature vapor pressure data and vaporization enthalpies (Table S3, in bold), selected literature liquid phase heat capacities (Table S2, in bold), and the ideal-gas heat capacities calculated in this work (Table S4) were treated with the Simultaneous Correlation (SimCor) routine116 to obtain a consistent thermodynamic description. See SimCor description in Section S2 in the SI for details on the treatment of the experimental thermodynamic data. The Clarke and Glew equation117 with four parameters was used for the SimCor description, which corresponds to a linear temperature dependence of the  vapC p term. The Clarke and Glew equation has the following form:

 vapG 0   p  1 1 T R ln 0     vap H 0        vapC p0     1  ln  p   T   T 0 T  T      vapC p        2 ln     2  T     T

   

,

(5)

where p0 = 100 kPa is a reference pressure, θ = 298.15 K is a reference temperature, R is the molar gas constant, and Δ vap G 0 and Δ vap H 0 are the standard molar vaporization Gibbs energy and enthalpy, respectively.

15



Page 16 of 46

FIGURE 2. Deviations of literature vapor pressure data from the SimCor results (pcor) for [emIm][NTf2]. Left: full scale, right: zoomed. , Zaitsau et al.112, , Emel'yanenko et al.108, , Lovelock et al. Rocha et

al.113,

, , Heym et

al.118,

, Ahrenberg et

al.106,,

Santos et

al.114.

115,

,

Data sets given with empty

symbols were excluded from the SimCor correlations.

[emIm][NTf2] vapor pressure data from four sources108, 112-114 were found to be in an acceptable agreement (within 30 %) as can be seen in Figure 2.  vap H to be included in the SimCor were selected based on the following criteria: i)  vap H were reported in sources, where vapor pressures were not available, ii)  vap H were preferably obtained by methods with higher reliability (vaporization calorimetry, QCM, and Knudsen effusion), iii)  vap H are mutually consistent with accepted vapor pressure data and with ΔvapCp. The temperature dependence of ΔvapCp displayed in Figure 3 indicates consistency with  vap H decreasing in the temperature range 350 to 450 K, but remaining almost constant over the temperature range 550 to 650 K. Finally, we selected 5 consistent sources12, 106-107, 109-110 of Δ vap H out of 11 available values.

16


Page 17 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FIGURE 3. Thermal quantities resulting from SimCor for [emIm][NTf2]. Left: vaporization enthalpy Δ vap H ; , Santos et al.12, , Armstrong et al.119, , Tolstoguzov120, , Luo et al.121, , Wang, et al.122,

, Lovelock et al.115, , Verevkin et al.110, , Verevkin et al.109, , Ahrenberg et al.107, , Ahrenberg et al.106, , Dunaev et al.123. Right: difference in heat capacity ΔvapCp; calorimetry + ab initio: , C pl from Paulechka et al.100, , C pl from Ferreira et al.103, approximations: , Zaitsau et al.112, Emel'yanenko et al.108, and Esperança et al.124, , Armstrong et al.119, +, Santos et al.12, , Rocha et al.113, , Verevkin et al.111, , Ahrenberg et al.106. , SimCor result. Data sets given with empty symbols were excluded from the SimCor correlations. See discussion in Section S3 for origin of the literature ΔvapCp values.

Parameters of equation (5) developed by the SimCor method with validity down to the triple point temperature of the ILs of interest are presented in Table 2.  vapG 0 is not available for [emIm][BF4] and [emIm][PF6], since no vapor pressure data were correlated. Temperature dependence of  vap H and  sub H for all three ILs is discussed in more detail in Section S3 of the SI and illustrated in Figures S1-S3. TABLE 2. Parameters of the Clarke and Glew equation (5), derived using the SimCor method at the reference temperature θ = 298.15 K and pressure p0 = 100 kPa. (Tmin – Tmax)a Δ G 0 / vap / -1 kJ·mol K [emIm] [BF4] [emIm] [PF6]

283 – 432

n/a

334 – 435

n/a

Δ vap H 0 /

Δ vap C p0 /

kJ·mol-1

J·K-1·mol-1

129.803 ± 2.2 141.238 ± 4.0

−68.23 ± 2.8 −85.213 ± 22

17


Δ vapC p0 T

/

J·K-2·mol-1 0.1687 ± 0.13 0.0363 ± 0.25

σb / mPa

σrb

n/a

n/a

n/a

n/a


[emIm] [NTf2]

257 – 625

80.650 ± 0.3 c

126.459 ± 1.0

−103.43 ± 4.0

Page 18 of 46

0.3174 ± 0.05

34

7.8 %

The parameters of the Clarke and Glew equation (5), developed by SimCor method are valid over a combined temperature interval of the input thermodynamic data. a

b

σ is the standard deviation of the fit defined as   

  p  n

2

i 1

i

1/ 2

( n  m)  , where Δp is the difference between the 

experimental and the smoothed values, n is the number of experimental points used in the fit and m is the number of adjustable parameters of Clarke and Glew equation (5), σr is the relative standard deviation of the fit defined as 1/2

 r    i 1   ln p  i ( n  m)  .   n

2

c The

uncertainties quoted are estimated expanded uncertainties (k=2) of the thermodynamic quantities and are reported with more digits to avoid round-off errors in calculations based on the correlation.

4. RESULTS AND DISCUSSION 4.1

Crystal properties at the absolute zero

When the calculated vibrational frequencies of crystalline and gaseous IL are compared, appreciable shifts (up to dozens of cm−1) in the position of the peaks can be observed, see Figures S4-S6 for details for all three ILs. Sublimation of the ionic pair, leading to the loss of the crystal field and also to a significant conformation change, translates to changes of frequencies of both C−H/N−H stretching and lower-frequency deformation modes, affecting thus both the zero-point energy and thermal contributions to sublimation enthalpy. These sublimation-related frequency shifts do not exhibit any positive/negative bias so they compensate out to a certain extent. TABLE 3 Cohesive properties (in kJ∙mol−1, given per ionic pair) of the three studied ILs at 0 K. Calculations based on optimized unit-cell geometries, phonons, and isolated ion pair geometries obtained at the PBE-D3(BJ)/PAW(400 eV) level unless noted otherwise and unit-cell energy – volume profiles refined at various levels of theory. [emIm][BF4] Method

−Elata

−EIPb

−Ecohc

∆EZPEd

∆subH(0 K)

PBE-D3(BJ)/PAW(400 eV)

483.9

322.0

161.7

−5.1

156.6


532.8

370.9

160.3

−4.9

155.4

18


Page 19 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


B3LYP-D3(BJ)/PAW(1000 eV)

564.8

363.9

199.0

−7.0

192.0

HMBI (MP2C-F12+pHF)

583.3

367.6

213.3

−8.3

205.0

HMBI (MP2C-F12+pB3LYP)

556.6

367.6

186.6

−7.3

179.3

[emIm][PF6] Method

−Elata

−EIPb

−Ecohc

∆EZPEd

∆subH(0 K)


516.2

339.1

172.8

−6.0

166.8


517.2

316.2

172.5

−6.0

166.5


555.1

359.3

212.6

−7.4

205.2

HMBI (MP2C-F12+pHF)

583.0

347.7

218.5

−8.2

210.3


558.6

347.7

194.1

−6.9

187.2

[emIm][NTf2] Method

−Elata

−EIPb

−Ecohc

∆EZPEd

∆subH(0 K)


495.5

322.7

159.1

4.2

163.3


493.6

326.2

157.8

3.9

161.7


534.7

325.7

185.4

−3.7

181.7

HMBI (MP2C-F12+pHF)

575.0

344.1

199.1

−4.5

194.6


550.0

344.1

174.1

−2.9

171.2

Lattice energy related to formation of the crystal from individual ions (in their gas phase conformations). Geometry of the isolated ions optimized only at the PBE-D3(BJ)/6-311+G(d,p) level of theory in Gaussian. b Interaction energy of the ionic pair in the most stable gas phase conformation. Geometry of the isolated ions optimized only at the PBE-D3(BJ)/6-311+G(d,p) level of theory in Gaussian. c Cohesive energy related to formation of the crystal from individual ion pairs (in their gas phase conformations). d Difference of the zero point vibrational energies of the gas and crystal. a

Table 3 summarizes obtained lattice energy Elat (corresponding to formation of the crystals from individual ions), energy of formation of ionic pairs EIP, cohesive energy Ecoh (corresponding to formation of the crystals from individual ion pairs), and sublimation enthalpy ∆subH(0 K). Ab initio calculations confirmed the prevailing ionic character of crystalline ILs, whose Elat range from −470 19



to −585 kJ∙mol−1, being closer to lattice energies of inorganic salts of larger monovalent ions (e.g. CsI)125 than to nonionic molecular crystals. Since the vapor of ILs is formed by neutral ionic pairs,126-127 energy of formation of these pairs at 0 K from isolated ions EIP needs to be accounted for to obtain Ecoh. This contribution, facilitating the sublimation, compensates out roughly two thirds of Elat. Resulting values of Ecoh calculated at the individual levels of theory exhibit a relatively large scatter (up to 65 kJ∙mol−1 or 40 % of Ecoh). It is not that surprising when we consider that only the corresponding pair interaction energies EIP are scattered over some 45 kJ∙mol−1, see Table 3. Pair interionic interaction energies for all ionic dimers which are separated by less than 10 Å in the crystal lattice of [emIm][BF4], calculated using MP2C-F12/aug-cc-pVDZ, HF/pob-TZVP,93 and B3LYP-D3(BJ,ATM)/pob-TZVP levels of theory are depicted in Figure S7. MP2C-F12 calculations predict (up to 6 kJ∙mol−1) less repulsive interactions of the same ions and more attractive interactions of the counterions at the closest contacts below 4 Å than B3LYP. This illustrates that the fragment-based additive approaches should be in principle capable to correct periodic DFT calculations of Ecoh for ionic materials. At longer distances, individual methods predict closer values of pair interaction energies, which still, however, exceed 100 kJ∙mol−1 in absolute value at the interionic separation 10 Å. This proves that a straightforward adoption of the computational procedures developed for non-ionic molecular crystals (rather depreciating the role of long-range regime) can be inadequate for predictions of the cohesion of crystalline ILs for which adequate treatment of both the short-range and the long-range interactions plays an important role. To assess the importance of the dispersion corrections within the MP2C method (i. e. the MP2C – MP2 difference), we evaluated the MP2-F12 cohesive energies for the structures corresponding to the minima of the Ecoh(V) curves. The given dispersion corrections to the lattice energy amount to 3.5 kJ∙mol−1, 3.3 kJ∙mol−1, and 12.4 kJ∙mol−1 for [emIm][BF4], [emIm][PF6], and [emIm][NTf2], respectively. Such positive values indicate that the bare MP2 overestimates the dispersion 20


Page 20 of 46

Page 21 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interactions also for the ionic species when the ions are formed by delocalized π-electron systems, e. g. the imidazolium core. Higher value of the correction for [emIm][NTf2] can be attributed to abundant π bonds in its anion. For a detailed statistics on the dispersion corrections, see Table S6. Table 3 reveals that the arguably lowest level of theory used, the PBE-D3(BJ)/PAW method with the 400 eV plane-wave energy cutoff yields the least bound crystals (except [emIm][NTf2]) while the remaining methods (using either a larger plane-wave basis set, or a hybrid DFT functional, or including ab initio wavefunction treatment of proximate ionic clusters) give significantly larger cohesive energies in absolute value. Comparison with experimental data is given in section 4.3. There are several factors that could cause this behavior. Comparison of Ecoh calculated at the B3LYP-D3(BJ)/PAW and HMBI (MP2C-F12+pB3LYP) levels of theory in Table 3 shows that recalculating the interaction energies of proximate ionic pairs in the crystal lattice by MP2C-F12 leads to lower Ecoh by 10-20 kJ∙mol−1. Such trends of calculated Ecoh show that: i) correcting the proximate pair interaction energies improves accuracy of Ecoh, and ii) many-body effects (mostly polarizability of the ions) are captured sufficiently in the HMBI model even for ionic materials. To suppress the basis set issues in the HMBI calculations and to minimize the uncertainty related to using the double-zeta basis set, explicit electronic correlation via the MP2C-F12 method was used here. As the plane-wave based DFT calculations using the B3LYP and PBE functionals yield very different Ecoh values (differing by up to 40 kJ∙mol−1), the quality of the underlying optimized-unit cell geometries needs to be investigated further. Although it has become common and considered as fairly accurate for non-ionic molecular crystals to optimize their unit-cell geometry with a lowerlevel method and to subject these geometries to higher-level single-point electronic energy calculations,58-59 this approach could fail in case of ionic (yet molecular) crystals of ILs. Pair interaction energies of proximate ions in the crystal lattice often exceed ±300 kJ∙mol−1 (see Figure S7), being significantly higher than in the case of non-ionic molecular crystals for which the closest 21



molecular pair interactions ranges from units of kJ∙mol−1 for non-polar molecules to dozens of kJ∙mol−1 for hydrogen-bonded clusters.41 A slight variation of the unit-cell geometry (equivalent to using a geometry previously optimized at a lower level of theory to rationalize the computational cost) can thus translate into higher error in the cohesive energy of ionic crystals, determining the significantly higher sensitivity of the single-point energy calculations of the cohesive energy. To test this hypothesis, unit-cell geometries for crystalline [emIm][PF6] and geometries of isolated [emIm][PF6] pairs were further optimized at the PBE-D3(BJ)/PAW level using the hard PAW potentials and plane-wave energy cut-offs 700 eV and 1000 eV. Also, geometries of isolated [emIm][PF6] ionic pairs were optimized using the PBE-D3(BJ) functional with the 6-311+G(d,p) basis set. Table 4, showing Ecoh calculated for [emIm][PF6] using various computational setups, reveals that the given hypothesis does partially explain the observed behavior of calculated Ecoh. If both the unit-cell and the isolated ion pair geometries are optimized with any size of the PAW basis set, variation of the plane-wave energy cut-off for subsequent single-point Ecoh refinements yields Ecoh values differing by no more than 1-2 kJ∙mol−1. However, when the underlying vapor-phase ion pair geometry is optimized using the same PBE-D3(BJ) functional but with the 6-311+G(d,p) basis set instead, selection of the plane-wave energy cut-off for the single-point Ecoh refinement becomes the accuracy determining step (Ecoh scattered over an interval of 30 kJ∙mol−1). In this case, smaller PAW basis sets (400 eV and 700 eV cut-off) yield lower |Ecoh| results which are also in the closest agreement with experimental ∆subH (as is discussed in section 4.3). This can be attributed to the fact that the smallest PAW basis set used here is the closest in size to the split-valence Gaussian atomic orbital 6-311+G(d,p) basis set. TABLE 4 Cohesive energy Ecoh (in kJ∙mol−1) of [emIm][PF6] obtained by PBE-D3(BJ)/PAW calculations with various plane-wave energy cutoff values and underlying unit-cell geometries. Ion pair optimization Unit cell optimization

PBE-D3(BJ)/6-311+g(d,p) 400 eV 700 eV 1000 eV 22


PBE-D3(BJ)/PAW 400 eV 700 eV 1000 eV

Page 22 of 46

Page 23 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Single-point Ecoh calculations 400 eV

177.0

177.1

176.8

172.8

172.6

172.4

700 eV

169.5

169.8

169.9

172.5

172.7

172.8

1000 eV

200.9

201.4

201.6

171.7

172.1

172.3

−Ecoh

−Ecoh

When the ion-pair geometry of [emIm][PF6] was optimized at the B3LYP-D3(BJ)/6-311+G(d,p) and MP2/aug-cc-pVDZ levels of theory, mimicking the levels of the theory used for single-point calculations of Ecoh for crystals in Table 3, and the resulting ion-pair electronic energies were combined with the electronic unit-cell energies obtained with equivalent methods, Ecoh values amounting to −211.3 kJ∙mol−1 and −206.1 kJ∙mol−1 were obtained, respectively. In the case of B3LYP results, switching from single-point energy calculations to proper geometry optimizations for the ion pairs does not lead to a substantial improvement (1.3 kJ∙mol−1 change in Ecoh). This suggests that the B3LYP functional favors somewhat different unit-cell geometry than the PBE one, being in agreement with the conclusions of Izgorodina et al.128 about the importance of the interplay of dispersion, electrostatic, induction and exchange interactions and their treatment within particular DFT functionals for geometries and interactions energies of ionic pairs of ILs. Altogether, this translates to considerably outlying Ecoh obtained from B3LYP single-point Ecoh refinements on PBE geometries. Furthermore, cohesive energy values given in Table 3 indicate a systematic trend in the row of increasing rate of inclusion of the exact electronic exchange in the model (PBE – B3LYP – HF). Among these, PBE results are the lowest in absolute value which proves also to be in the closet agreement with experiment (see section 4.3). As inclusion of more exact exchange leads to larger errors here, this behavior cannot be directly attributed to the delocalization error which is typical for generalized-gradient approximation (GGA) functionals overbinding and artificially favoring ionic crystals over neutral co-crystal structures.129

23



To investigate the varying performance of individual methods further, interaction energies of the ionic pairs from the PBE-D3(BJ)/PAW optimized [emIm][BF4] crystal lattice were calculated at the PBE-D3(BJ), B3LYP-D3(BJ), HF, and CCSD(T)-F12 levels of theory, always with the augcc-pVDZ basis set. Figure 4 illustrates trends of the interaction energies of the [BF4]∙∙∙[BF4] and [BF4]∙∙∙[emIm] ionic pairs, compared to the reference CCSD(T)-F12 results. Calculated interaction energies indicate that models including higher rates of exact exchange exhibit worse performance, most importantly PBE interaction energies are in a better agreement with the CCSD(T)-F12 counterparts than the B3LYP results at interionic separation smaller than 4-5 Å. This is in contrast with the findings of the benchmark studies of interactions of ionic liquids by Izgorodina et al.,130131

reporting that GGA functionals exhibit appreciably poorer performance than the hybrid

functionals for calculations of interaction energies in ILs-related systems. This discrepancy can be, however, possibly explained by using the PBE-optimized geometries in this work of the ion pairs which can be less compatible with B3LYP interaction energy calculations.128 While interaction energies of the [BF4]∙∙∙[BF4] ionic pair calculated by all methods converge to the reference CCSD(T)-F12 results at longer interionic separation, results of all of DFT, HF, and MP2C calculations of [BF4]∙∙∙[emIm] interactions exhibit a linear drift from the CCSD(T)-F12 reference at longer interionic separations. This could be an imprint of the charge transfer occurring only between the cation – anion ionic pair which is possibly better captured by CCSD(T)-F12. The charge transfer varies the effective partial charges at the ions determining the long-range electrostatics. Analogous plots of the calculated interactions of the [emIm] and [NTf2] ions are given in Figure S8 in the case of which PBE and B3LYP exhibit a comparable performance.

24


Page 24 of 46

Page 25 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FIGURE 4. Percentage deviations of the interaction energies of the [BF4]∙∙∙[BF4] (left) and [BF4]∙∙∙[emIm] (right) ionic pairs calculated by various methods from the reference results of the CCSD(T)-F12/aug-cc-pVDZ calculations.

To conclude, the size of the plane wave basis set is not the key factor affecting the accuracy of the calculated Ecoh only if the geometries of both phases are optimized consistently using exactly the same level of theory, including the basis set type. Prohibitive cost of performing unit-cell optimizations with hybrid functionals or ab initio wavefunction techniques hinders us to investigate the effect of the quality of unit-cell geometries thoroughly. Nevertheless, our results indicate that using single-point Ecoh refinements can lead to misleading Ecoh results for ionic crystals. In contrast to non-ionic crystals, the chemical entities forming the vapor phase of ILs are not single covalentlybound molecules but ionic pairs. Consequently, the quality of the geometry of the vapor phase ionic pairs becomes an important factor strongly affecting the overall accuracy of Ecoh. 4.2 Ideal-gas properties Conformation analysis at the PBE-D3(BJ)/6-311+G(d,p) level revealed that both the [emIm] cation and the [NTf2] anion form two dynamically stable conformations if isolated – [emIm] with the ethyl group in planar or non-planar configuration; and [NTf2] with the trifluoromethyl groups being 25



Page 26 of 46

aligned in a cis or trans position. Conformational energies calculated in this work, listed in Table 5, are in a good agreement with previous conformational studies.73 TABLE 5 Calculated energies (kJ mol−1) of individual conformations of isolated [emIm] cation and [NTf2] anion. Data correspond to PBE-D3(BJ)/6-311+G(d,p) optimized geometries whose energies were further refined at the given levels of theory. [emIm] Conformation

[NTf2]

Non-planar

Planar

Trans

Cis

2

1

2

2

PBE-D3(BJ)/6-311+G(d,p)

0.00

1.72

0.00

3.40

B3LYP-D3(BJ)/6-311+G(d,p)

0.00

2.17

0.00

3.75

MP2-F12/aug-cc-pVDZ

0.00

2.25

0.00

3.80

CCSD(T)/CBS

0.00

1.81

0.00

3.59

Isomers

FIGURE 5. Contributions of the interaction cation – anion to the heat capacity of IL vapor calculated using molecular dynamics simulation with the CL&P all-atom force field. Solid lines are correlations of the data using polynomial functions. The dashed gray line shows the 4.5 R hightemperature limit.

26


Page 27 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Reliable molecular-dynamics calculations of properties highly sensitive to numerical noise, such as the heat capacity and its individual contributions, require generating extremely long (1 μs) trajectories. In case of isolated ions or an ion pair, this is still computationally feasible despite the relatively poor sampling of the phase space for such a small simulated ensemble. Figure 5 illustrates g0 the calculated contributions of the cation – anion interaction to the total heat capacity ( Cinter ) of the

IL vapor. Although the particular computed data points exhibit appreciable scatter, there are visible g0 trends in Cinter (T) behavior that can be further interpreted. Analyzing the numbers of individual

types of degrees of freedom within an ionic pair, there are six degrees corresponding to mutual movement of the cation and the anion. Three of these could be considered as low-frequency harmonic oscillations of the centers of mass of the bare ions in three Cartesian directions while the other three modes could be approximated as mutual hindered rotations of the ions. Following this simplification, the high-temperature limit of the heat capacity for such six degrees of freedom is 4.5 R (37.4 J K−1 mol−1). Distribution of the allowed energy levels of hindered rotations typically translates to temperature trend of the corresponding heat capacity contribution exhibiting a maximum at lower temperatures, and only after that gradually converging to the high-temperature limit for the heat capacity of a free rotor (1/2 R).132 Both features can be observed for the calculated g0 (T) curve for [emIm][PF6]. As the hexafluorophosphate anion is the most spherical structure Cinter

among the studied anions, the approximate model of three mutual vibrations of rigid ions and three hindered rotations of the anion with respect to the cation provides the best interpretation for g0 [emIm][PF6]. Similar behavior can be seen also for [emIm][BF4], only with a slower Cinter decay

at higher temperatures. Low sphericity of the [NTf2] anion, and possibly higher barriers to its g0 rotation with respect to the cation, can explain the different shape of the [emIm][NTf2] Cinter (T)

curve (convergence is expected at much higher temperatures than considered).

27



For [emIm][NTf2], data on C pg0,pair calculated in this work can be compared with the results of Paulechka et al.73 who also used DFT-based calculations of molecular parameters which were subsequently processed by the statistical-thermodynamic models to obtain the thermodynamic properties of ideal-gaseous ILs. Standard molar entropy S g0 of [emIm][NTf2], equal to 972.9 J K−1 mol−1 at 470 K and 100 kPa as calculated in ref.73 is 40 J K−1 mol−1 lower (equivalent to 4%) than the value calculated in this work (see Table S5 for S g0 at selected temperatures). This discrepancy arises from different approximations in the utilized computational models. In contrast to this work, Paulechka et al.73 used: i) a combination of experimental and scaled B3LYP/6-31+G(2df,p) frequencies although the development of the given scale factors for the lowest frequency modes suffered from limited availability of experimental data; ii) only a very limited number of conformations of the ionic pair assumed to be present in the vapor phase; iii) the one-dimensional hindered rotor model to account for only one of the degrees of freedom corresponding to the cation – anion interaction. Paulechka et al.73 report experimental values of S g0 at 470 K, based on thermodynamic data from refs.12,

34, 73, 112

in the range 997 – 1027 J K−1 mol−1, while the

corresponding S g0 value obtained by the SimCor method in this work amounts to 992.5 J K−1 mol−1. This means that the errors of S g0 values calculated in both studies possess opposite signs, however being comparable in their absolute values. The SimCor method enabled us to derive the experimental value S g0 (298.15 K) = 782.3 J K−1 mol−1, compared to which the value calculated in ref.73 is underestimated by −4.4 J K−1 mol−1 while the value obtained in this work is overestimated by 20 J K−1 mol−1. On the basis of the comparison of calculated and experimental S g0 , the estimated standard uncertainty of calculated C pg0,pair and S g0 amounts to 3%.

28


Page 28 of 46

Page 29 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.3 Finite-temperature crystal properties Literature experimental data for benchmarking the accuracy of calculated thermodynamic properties of crystalline ILs are scarce. Figure 6 shows comparison of plane-wave based and HMBI-calculated molar volumes at zero pressure with experimental data.80-82 While full PBED3(BJ)/PAW calculations overestimate the molar volumes for all three ILs (by 2-8%), both B3LYP-D3(BJ)/PAW

and

MP2C-F12/avdz+pB3LYP

HMBI

single-point

refinements

underestimate the volumes by 2-5% and 1-4%, respectively. This indicates that the HMBI model (using MP2C-F12/avdz+pB3LYP) works reliably to predict finite-temperature molar volumes of crystalline ILs. Figures S9-S11 reveal that the treatment of the long-range forces within the HMBI model with the pHF method underestimates the molar volumes more than in the case of pB3LYP. Calculated molar volumes of crystalline ILs at selected temperatures are listed in Table S7.

FIGURE 6. Comparison of experimental (diamonds, ♦ – ref.80, ♦ – ref.81, ♦ – ref.82) and quasiharmonic molar volumes calculated for the studied crystalline ILs at zero pressure using the PBED3(BJ)/PAW optimized geometries and phonons. Dotted, solid and dashed lines correspond to no single-point refinement, B3LYP-D3(BJ)/PAW refinement and HMBI MP2C-F12/aug-ccpVDZ+pB3LYP refinements of the Ecoh(V) curves, respectively.

29



Ordering of molar volumes of all studied ILs is correctly captured by theory, with [emIm][NTf2] exhibiting the loosest crystal packing due to lower symmetry and larger size of its anions. To our knowledge, no temperature-dependent experimental data on unit cell volumes or crystal densities exist, so the calculated thermal expansivity cannot be verified. Regardless of the computational method, calculations here predict that [emIm][BF4] and [emIm][PF6] exhibit the largest and the smallest thermal expansivity, respectively. This trend can be tracked to the slopes of the Ecoh(V) curves, out of which the one for [emIm][BF4] has the least steep expansion branch. Figure 7 compares the predicted compressibility of the three ILs at 312.8 K which is the only temperature at which any experimental data are available (for [emIm][PF6]133). PBE-D3(BJ)/PAW level predicts the steepest V(p) curve for all three ILs while B3LYP-D3(BJ)/PAW gives the least steep curves. The only two experimental points for [emIm][PF6] indicate that the compressibility is predicted reasonably, at least in the sub-GPa region.

FIGURE 7. Comparison of experimental (diamonds, ♦ – ref.133) and quasi-harmonic molar volumes calculated for the studied crystalline ILs at 312.8 K using the PBE-D3(BJ)/PAW optimized geometries and phonons. Dotted, solid and dashed lines correspond to no single-point refinement, B3LYP-D3(BJ)/PAW refinement and HMBI MP2C-F12/aug-cc-pVDZ+pB3LYP refinements of the Ecoh(V) curves, respectively.

30


Page 30 of 46

Page 31 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8 compares the computed and experimental isobaric heat capacities for crystalline and vapor phases. Calculated C pcr are appreciably underestimated when compared to experiment (11% or 32 J∙K−1∙mol−1 for [emIm][PF6]134 and 7% or 27 J∙K−1∙mol−1 for [emIm][NTf2]100 at 270 K). This usually is a result of two factors – overestimated DFT-calculated vibrational frequencies and underestimated thermal expansivity originating from too steep Ecoh(V) curves. Such an uncertainty in C pcr is estimated to translate to a 7 kJ∙mol−1 error in the thermal contribution to the sublimation enthalpy at 300 K, being already a non-negligible source of error. As all of the calculated C pcr data sets are based on the same phonons, their scatter (due to variations of the shape of the Ecoh(V) curves) is small, not reaching the deviations from the experiment, see Figures S12-S14. Calculated heat capacity data for crystalline ILs at selected temperatures are listed in Table S8.

FIGURE 8. Comparison of calculated (PBE-D3(BJ)/PAW, dashed lines) and experimental (diamonds, ♦ – data for both polymorphs from ref.134 , ♦ – data from ref.100) isobaric heat capacities for the crystalline ILs, and of isobaric heat capacities calculated for the vapor (PBE-D3(BJ)/6311+G(d,p)+MM, solid lines). Gas phase heat capacities obtained only within the rigid rotor – harmonic oscillator are given with dotted lines.

31



When the heat capacities of the crystalline and vapor phases are compared, C pg0 and C pcr of [emIm][BF4] and [emIm][PF6] exhibit similar behavior as for most non-ionic crystals, as visible in Figure 8. At very low temperatures, C pg0 is higher than C pcr due to facile excitation of translation and rotation modes in the vapor phase, whereas C pcr becomes higher than C pg0 only above some threshold temperature (154 K and 135 K for [emIm][BF4] and [emIm][PF6], respectively). This g0 trend is further amplified here due to the shapes of the Cinter (T) curves, see Figure 8. Heat capacities

of crystalline and ideal-gaseous [emIm][NTf2] do not follow this usual pattern with C pg0 being somewhat higher than C pcr at all temperatures of interest, which occurs mainly due to the higher g0 contribution to C pg0 of [emIm][NTf2]. Cinter

FIGURE 9. Calculated trends of ∆subH of three studied ILs (solid lines) along with the experimental ∆subH and its expanded uncertainty (diamonds), see Section 3 for a discussion of the reference experimental data. Optimized geometries and phonons obtained at the PBE-D3(BJ)/PAW level, vapor phase treated by the QM/MM approach described above.

This behavior of heat capacities propagates to temperature trends of the sublimation enthalpies which are illustrated in Figure 9. Only the temperature-dependent sublimation enthalpies of

32


Page 32 of 46

Page 33 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[emIm][BF4] and [emIm][PF6] exhibit a maximum. Nevertheless, it holds for all three studied ILs that the sublimation enthalpy at the temperature of the triple point is higher than at 0 K. Table 6 compares calculated and experimental sublimation enthalpies at the temperatures of the triple point. Note that Table 6 summarizes only results calculated consistently at the PBE-D3(BJ)/PAW level of theory as these sublimation enthalpies are in the closest agreement with experiment. In contrast, single-point B3LYP or MP2C energy refinements were shown above to outperform the PBE calculations in terms of accuracy of the molar volumes and heat capacities of crystalline ILs. This is probably due to the fact that the molar volumes and heat capacities do not depend on the value of the cohesive energy, but both are affected by the slope (first derivative) of the cohesive energy – volume Ecoh(V) curve for the crystal. Although the single-point energy refinements were observed to give misleading cohesive energies for the given ionic materials, the slopes of the Ecoh(V) curves could be affected negligibly or even beneficially at the same time, so that additional compensation of errors due to the electronic cohesive energy and due to the phonons occurs, leading to more accurate densities and heat capacities. TABLE 6 Cohesive properties (given per ionic pair) of the three studied ILs at the temperature of fusion. All energetic quantities given in kJ∙mol−1. [emIm][BF4], Tfus = 287.0 K (average from refs.102, 135-137) ∆Htherma

 sub H Tfus 

psub Tfus 

 vap H Tfus 

Calculationb

2.8

158.2

4.87∙10−10 Pa

146.9c

Experimentd

n/a

140.2 ± 1.3

n/a

130.6 ± 1.3

Method

[emIm][PF6], Tfus = 334.1 K (average from refs.134, 138-139) ∆Htherma


psub Tfus 


Calculationb

0.1

166.5

4.26∙10−9 Pa

145.4c

Experimentd

n/a

156.0 ± 2.3

n/a

138.2 ± 2.3

Method

33



Page 34 of 46

[emIm][NTf2], Tfus = 271.4 K100 ∆Htherma


psub Tfus 


Calculationb

3.9

165.6

7.68∙10−13 Pa

149.3e

Experimentd

n/a

151.2 ± 0.7

(4.6 ± 0.3)∙10−12 Pa

129.3 ± 0.5

Method

a

g0

cr

Thermal contribution to sublimation enthalpy at Tfus, given as an integral of the C p − C p term.

Data calculated using the PBE-D3(BJ)/PAW level of theory with the plane-wave energy cutoff 400 eV for optimizing the geometries and calculating phonons, and 1000 eV for refining the cohesive energies. c Estimated using the MD result on ∆ H from ref.17 fus d See Section 3 and Sections S1-S3 for details on sources of experimental data. e Estimated using the MD result on ∆ H from this work. fus b

Sublimation enthalpies at the temperature of the triple point calculated at the PBE-D3(BJ)/PAW level of theory are overestimated for all three ILs compared to experiment. Namely, the errors of calculated ΔsubH amount to 18, 11, and 14 kJ∙mol−1 (13%, 6.7%, and 9.5%) for [emIm][BF4], [emIm][PF6], and [emIm][NTf2], respectively. Such an uncertainty still exceeds the chemicalaccuracy level (4 kJ∙mol−1) targeted for high-level ab initio calculations of thermodynamic properties of crystals. However, the percentage errors range to values that are typical for the uncertainties of cohesive properties calculated for non-ionic molecular crystals.39, 53-55, 59 While there are no literature data on vapor pressures for [emIm][BF4] and [emIm][PF6], a relatively good agreement of the psub value obtained for [emIm][NTf2] from correlation of experimental data and from calculation at the PBE-D3(BJ)/PAW level of theory can be seen in Table 6. Calculated psub at the melting point temperature is underestimated only by a factor of 6, indicating the computational error being well within a single order of magnitude which fulfills the criterion for a successful first-principles prediction of psub.61, 140 Since a 14 kJ∙mol−1 error in calculated ΔsubH occurs for [emIm][NTf2], such a good agreement in terms of psub should be, however, attributed to a fortuitous error compensation between calculated ΔsubH and ΔsubS. In other words, the error of psub amounts to its minimum value in the vicinity of the triple-point temperature. Note that this

34


Page 35 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


error is strongly temperature-dependent due to the 14 kJ∙mol−1 error in calculated ΔsubH which predetermines the slope of psub. First-principles calculations of ΔfusH for molecular crystals remain unaffordably expensive. To provide a comprehensive comparison of the theoretical and experimental phase change enthalpies, molecular dynamics-based ΔfusH values were subtracted from the PBE-D3(BJ)/PAW ΔsubH values to obtain calculated ΔvapH at the temperature of the melting point, the ultimate goal of our research effort. Comparison with experiment in Table 6 reveals that the errors of such calculated ΔvapH amount to 16, 7, and 20 kJ∙mol−1 (12%, 5.2%, and 15%) for [emIm][BF4], [emIm][PF6], and [emIm][NTf2], respectively, meaning that inclusion of the MD-based ΔfusH somewhat alters the errors compared to the error values for ΔsubH, but without a significant augmenting or lessening trend. Expanded uncertainties of experimental ΔvapH given in Table 6 are up to one-order-of-magnitude smaller than the errors of calculated ΔvapH. Note that these experimental uncertainties are mainly based on the original uncertainty values reported by the authors (that are generally overly optimistic) and we used only thermodynamically consistent experimental data to derive these uncertainties within the simultaneous correlation procedure. Such a selective critical approach is strongly data-demanding and requires a thorough understanding of the underlying experimental techniques. Scatter of all published data on ΔvapH (see Figures 3, S1-S3), regardless of their mutual consistency, is fully comparable with the computational errors obtained in this work. 5. CONCLUSIONS We adapted the first-principles computational methodology, originally developed for molecular crystals, to be used for predictions of the cohesive energy and sublimation enthalpy and pressure for crystalline ionic liquids. The results calculated for three archetypal ionic liquids reveal that it is possible to target comparable accuracy levels (chemical or better) for the given ionic crystals of 35



molecular ions as is commonly done for non-ionic molecular crystals following the quasi-harmonic first-principles methodologies. For non-ionic molecular crystals, it is crucial to use highlycorrelated state-of-the-art ab initio wavefunction techniques or reliably parametrized semi-empiric dispersion corrections in DFT to accurately account for the interactions of proximate molecules in the crystal lattice, while the treatment of the long-range interactions usually plays a lesser role, yet being important to reach the sub-kJ mol−1 accuracy. Due to the prevailing Coulombic character and a large magnitude of the cohesive forces and their long-range character in crystals of ionic liquids, significance of the long-range interactions matches the importance of the short-range interactions, unlike for non-ionic crystals. An appropriate treatment of the long-range and many-body forces proves to be of an utmost importance to converge the cohesive properties of crystalline ionic liquids. Although the molar volumes and isobaric heat capacities of the crystals are captured with a fair accuracy using both periodic DFT-D or ab initio HMBI approaches, the computations of the cohesive properties exhibit a considerable sensitivity to any inconsistencies between the solid and gas phases, such as different levels of theory used for geometry optimizations or calculations of vibrational frequencies. Importantly, performing single-point energy refinements for geometries previously optimized with a cheaper method is shown to give misleading cohesive energies. Calculations of the thermodynamic properties of the ideal-gaseous phase of ionic liquids, combining the DFT calculations of molecular parameters of isolated ions, essential statisticalthermodynamic models assuming existence of several conformers of bare ions and classical molecular-dynamics simulation of the cation – anion interaction proved to be capable of providing reasonable results that are in a fair consistency with experimentally obtained data. The interionic contribution to the vapor phase heat capacity is the only part of our predictions of the sublimation enthalpy that is not covered exclusively from first principles. A future development of this methodology, exploiting ab initio molecular dynamics, is therefore expected.

36


Page 36 of 46

Page 37 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In context of the large magnitude of the cohesive forces binding the crystals of ILs, the observed computational uncertainty of finite-temperature ΔsubH data below 20 kJ∙mol−1 (or around 10%) can be regarded as a partial computational success. The resulting uncertainty is fully comparable with the typical spread of available experimental data on vaporization enthalpies of ionic liquids obtained by various research groups and different experimental techniques. Our results confirm that a combination of ΔsubH predicted from first principles and ΔfusH obtained experimentally can yield relatively accurate data on ΔvapH for ionic liquids at the temperatures of the melting point. Still, future development will be required to improve the ab initio cohesive properties of ionic (yet molecular) crystals further. In this field, highly accurate (but affordable for larger systems at the same time) periodic quantum-chemical calculations capable to treat the long-range, many-body and dispersion interactions at the same time with a high accuracy would be desirable above all.

ASSOCIATED CONTENT Supporting Information The supporting information is available free of charge at DOI: 1. Details on literature experimental data on thermodynamic properties of ILs. 2. Description of the SimCor procedure and underlying thermodynamic relationships. 3. Evaluation of the temperature dependence of vaporization enthalpies. 4. Calculated properties of ideal-gaseous ILs. 5. Calculated properties of crystalline ILs.

AUTHOR INFORMATION Corresponding Author * E-mail: [email protected]

37



ACKNOWLEDGMENT The authors acknowledge financial support from the Czech Science Foundation (GACR No. 1904150Y). The access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the program "Projects of Large Infrastructure for Research, Development, and Innovations" (LM2010005) and the CERIT-SC under the program Centre CERIT Scientific Cloud, part of the Operational Program Research and Development for Innovations, Reg. no. CZ.1.05/3.2.00/08.0144 is highly appreciated.

REFERENCES (1) MacFarlane, D. R.; Kar, M.; Pringle, J. M., Fundamentals of ionic liquids: from chemistry to applications. Wiley‐VCH: 2017. (2) Plechkova, N. V.; Seddon, K. R., Applications of ionic liquids in the chemical industry. Chem. Soc. Rev. 2008, 37, 123-150. (3) Deetlefs, M.; Fanselow, M.; Seddon, K. R., Ionic liquids: the view from Mount Improbable. RSC Advances 2016, 6, 4280-4288. (4) Paulechka, Y. U.; Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A., Vapor pressure and thermal stability of ionic liquid 1-butyl-3- methylimidazolium Bis(trifluoromethylsulfonyl)amide. Thermochim. Acta 2005, 439, 158-160. (5) Earle, M. J.; Esperanca, J. M. S. S.; Gilea, M. A.; Lopes, J. N. C.; Rebelo, L. P. N.; Magee, J. W.; Seddon, K. R.; Widegren, J. A., The distillation and volatility of ionic liquids. Nature 2006, 439, 831-834. (6) Rocha, M. A. A.; Santos, L. M. N. B. F., First volatility study of the 1-alkylpyridinium based ionic liquids by Knudsen effusion. Chem. Phys. Lett. 2013, 585, 59-62. (7) Kirchner, B.; Holloczki, O.; Lopes, J. N. C.; Padua, A. A. H., Multiresolution calculation of ionic liquids. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2015, 5, 202-214. (8) Maginn, E. J., Atomistic simulation of the thermodynamic and transport properties of ionic liquids. Acc. Chem. Res. 2007, 40, 1200-1207. (9) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J., Molecular modeling and experimental studies of the thermodynamic and transport properties of pyridinium-based ionic liquids. J. Phys. Chem. B 2006, 110, 2821-2832. (10) Canongia Lopes, J. N.; Padua, A. A. H., CL&P: A generic and systematic force field for ionic liquids modeling. Theor. Chem. Acc. 2012, 131, 1129. (11) Dequidt, A.; Devemy, J.; Padua, A. A. H., Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator. J. Chem. Inf. Model. 2016, 56, 260-268. (12) Santos, L. M. N. B. F.; Lopes, J. N. C.; Coutinho, J. A. P.; Esperanca, J. M. S. S.; Gomes, L. R.; Marrucho, I. M.; Rebelo, L. P. N., Ionic liquids: First direct determination of their cohesive energy. J. Am. Chem. Soc. 2007, 129, 284-285. (13) Rane, K. S.; Errington, J. R., Using Monte Carlo Simulation to Compute Liquid-Vapor Saturation Properties of Ionic Liquids. J. Phys. Chem. B 2013, 117, 8018-8030.

38


Page 38 of 46

Page 39 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(14) Rai, N.; Maginn, E. J., Critical behaviour and vapour-liquid coexistence of 1-alkyl-3methylimidazolium bis(trifluoromethylsulfonyl)amide ionic liquids via Monte Carlo simulations. Faraday Discuss. 2012, 154, 53-69. (15) Koeddermann, T.; Paschek, D.; Ludwig, R., Ionic liquids: Dissecting the enthalpies of vaporization. Chemphyschem 2008, 9, 549-555. (16) Kelkar, M. S.; Maginn, E. J., Calculating the enthalpy of vaporization for ionic liquid clusters. J Phys. Chem. B 2007, 111, 9424-9427. (17) Červinka, C.; Pádua, A. A. H.; Fulem, M., Thermodynamic Properties of Selected Homologous Series of Ionic Liquids Calculated Using Molecular Dynamics. J. Phys. Chem. B 2016, 120, 2362-2371. (18) Verevkin, S. P.; Zaitsau, D. H.; Emel'yanenko, V. N.; Yermalayeu, A. V.; Schick, C.; Liu, H.; Maginn, E. J.; Bulut, S.; Krossing, I.; Kalb, R., Making Sense of Enthalpy of Vaporization Trends for Ionic Liquids: New Experimental and Simulation Data Show a Simple Linear Relationship and Help Reconcile Previous Data. J. Phys. Chem. B 2013, 117, 6473-6486. (19) Vilas, M.; Rocha, M. A. A.; Fernandes, A. M.; Tojo, E.; Santos, L. M. N. B. F., Novel 2alkyl-1-ethylpyridinium ionic liquids: synthesis, dissociation energies and volatility. Phys. Chem. Chem. Phys. 2015, 17, 2560-2572. (20) Ishizuka, R.; Matubayasi, N., Effective Charges of Ionic Liquid Determined SelfConsistently through Combination of Molecular Dynamics Simulation and Density-Functional Theory. J. Comput. Chem. 2017, 38, 2559-2569. (21) Zhang, Y.; Maginn, E. J., A Simple AIMD Approach to Derive Atomic Charges for Condensed Phase Simulation of Ionic Liquids. J. Phys. Chem. B 2012, 116, 10036-10048. (22) Addicoat, M. A.; Stefanovic, R.; Webber, G. B.; Atkin, R.; Page, A. J., Assessment of the Density Functional Tight Binding Method for Protic Ionic Liquids. J. Chem. Theory Comput. 2014, 10, 4633-4643. (23) Grimme, S.; Hujo, W.; Kirchner, B., Performance of dispersion-corrected density functional theory for the interactions in ionic liquids. Phys. Chem. Chem. Phys. 2012, 14, 48754883. (24) Halat, P.; Seeger, Z. L.; Acevedo, S. B.; Izgorodina, E. I., Trends in Two- and ThreeBody Effects in Multiscale Clusters of Ionic Liquids. J. Phys. Chem. B 2017, 121, 577-588. (25) Hunt, P. A., Quantum Chemical Modeling of Hydrogen Bonding in Ionic Liquids. Top. Curr. Chem. 2017, 375, 59. (26) Malberg, F.; Brehm, M.; Holloczki, O.; Pensado, A. S.; Kirchner, B., Understanding the evaporation of ionic liquids using the example of 1-ethyl-3-methylimidazolium ethylsulfate. Phys. Chem. Chem. Phys. 2013, 15, 18424-18436. (27) Liu, J. F.; He, X., Accurate prediction of energetic properties of ionic liquid clusters using a fragment-based quantum mechanical method. Phys. Chem. Chem. Phys. 2017, 19, 2065720666. (28) Ludwig, R., Thermodynamic properties of ionic liquids - a cluster approach. Phys. Chem. Chem. Phys. 2008, 10, 4333-4339. (29) Strate, A.; Niemann, T.; Ludwig, R., Controlling the kinetic and thermodynamic stability of cationic clusters by the addition of molecules or counterions. Phys. Chem. Chem. Phys. 2017, 19, 18854-18862. (30) Seeger, Z. L.; Kobayashi, R.; Izgorodina, E. I., Cluster approach to the prediction of thermodynamic and transport properties of ionic liquids. J. Chem. Phys. 2018, 148, 193832. (31) Halat, P.; Seeger, Z. L.; Acevedo, S. B.; Izgorodina, E. I., Trends in Two- and ThreeBody Effects in Multiscale Clusters of Ionic Liquids. Journal of Physical Chemistry B 2017, 121, 577-588. 39



(32) Groom, C. R.; Bruno, I. J.; Lightfoot, M. P.; Ward, S. C., The Cambridge Structural Database. Acta Cryst. B 2016, 72, 171-179. (33) Kazakov, A.; Magee, J. W.; Chrico, R. D.; Paulechka, E.; Diky, V.; Muzny, C. D.; Kroenlein, K.; Frenkel, M. NIST Standard Reference Database 147: ILThermo. (34) Armstrong, J. P.; Hurst, C.; Jones, R. G.; Licence, P.; Lovelock, K. R. J.; Satterley, C. J.; Villar-Garcia, I. J., Vapourisation of ionic liquids. Phys. Chem. Chem. Phys. 2007, 9, 982-990. (35) Shi, R.; Wang, Y. T., Dual Ionic and Organic Nature of Ionic Liquids. Sci. Reports 2016, 6, 19644. (36) Izgorodina, E. I.; Golze, D.; Maganti, R.; Armel, V.; Taige, M.; Schubert, T. J. S.; MacFarlane, D. R., Importance of dispersion forces for prediction of thermodynamic and transport properties of some common ionic liquids. Phys. Chem. Chem. Phys. 2014, 16, 72097221. (37) Serra, P. B. P.; Ribeiro, F. M. S.; Rocha, M. A. A.; Fulem, M.; Ruzicka, K.; Coutinho, J. A. P.; Santos, L., Solid-liquid equilibrium and heat capacity trend in the alkylimidazolium PF6 series. J. Mol. Liq. 2017, 248, 678-687. (38) Beran, G. J. O.; Wen, S.; Nanda, K.; Huang, Y.; Heit, Y., Accurate and Robust Molecular Crystal Modeling Using Fragment-Based Electronic Structure Methods. Top. Curr. Chem. 2014, 345, 59-93. (39) Červinka, C.; Fulem, M., State-of-the-Art Calculations of Sublimation Enthalpies for Selected Molecular Crystals and Their Computational Uncertainty. J. Chem. Theory Comput. 2017, 13, 2840-2850. (40) Beran, G. J. O.; Nanda, K., Predicting Organic Crystal Lattice Energies with Chemical Accuracy. J. Phys. Chem. Lett. 2010, 1, 3480-3487. (41) Červinka, C.; Fulem, M.; Růžička, K., CCSD(T)/CBS Fragment-Based Calculations of Lattice Energy of Molecular Crystals. J. Chem. Phys. 2016, 144, 064505. (42) Nolan, S. J.; Bygrave, P. J.; Allan, N. L.; Manby, F. R., Comparison of the Incremental and Hierarchical Methods for Crystalline Neon. J. Phys. Condens. Matter 2010, 22, 074201. (43) Ringer, A. L.; Sherrill, C. D., First Principles Computation of Lattice Energies of Organic Solids: the Benzene Crystal. Chem. Eur. J. 2008, 14, 2542-2547. (44) Rosciszewski, K.; Paulus, B.; Fulde, P.; Stoll, H., Ab Initio Coupled-Cluster Calculations for the FCC and HCP Structures of Rare-Gas Solids. Phys. Rev. B 2000, 62, 5482-5488. (45) Sode, O.; Hirata, S., Second-Order Many-Body Perturbation Study of Solid Hydrogen Fluoride. J. Phys. Chem. A 2010, 114, 8873-8877. (46) Wen, S.; Beran, G. J. O., Accurate Molecular Crystal Lattice Energies from a Fragment QM/MM Approach with On-the-Fly Ab Initio Force Field Parametrization. J. Chem. Theory Comput. 2011, 7, 3733-3742. (47) Yang, J.; Hu, W. F.; Usvyat, D.; Matthews, D.; Schutz, M.; Chan, G. K. L., Ab Initio Determination of the Crystalline Benzene Lattice Energy to Sub-kilojoule/mole Accuracy. Science 2014, 345, 640-643. (48) Brandenburg, J. G.; Grimme, S., Organic Crystal Polymorphism: A Benchmark for Dispersion-Corrected Mean-Field Electronic Structure Methods. Acta Cryst. B 2016, 72, 502513. (49) Carter, D. J.; Rohl, A. L., Benchmarking Calculated Lattice Parameters and Energies of Molecular Crystals Using van der Waals Density Functionals. J. Chem. Theory Comput. 2014, 10, 3423-3437. (50) Deringer, V. L.; George, J.; Dronskowski, R.; Englert, U., Plane-Wave Density Functional Theory Meets Molecular Crystals: Thermal Ellipsoids and Intermolecular Interactions. Acc. Chem. Res. 2017, 50, 1231-1239. 40


Page 40 of 46

Page 41 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(51) Klimeš, J., Lattice Energies of Molecular Solids from the Random Phase Approximation with Singles Corrections. J. Chem. Phys. 2016, 145, 094506. (52) Maschio, L.; Usvyat, D.; Schutz, M.; Civalleri, B., Periodic Local Moller-Plesset Second Order Perturbation Theory Method Applied to Molecular Crystals: Study of Solid NH3 and CO2 Using Extended Basis Sets. J. Chem. Phys. 2010, 132, 134706. (53) Moellmann, J.; Grimme, S., DFT-D3 Study of Some Molecular Crystals. J. Phys. Chem. C 2014, 118, 7615-7621. (54) Otero-de-la-Roza, A.; Johnson, E. R., A Benchmark for Non-Covalent Interactions in Solids. J. Chem. Phys. 2012, 137, 054103. (55) Reilly, A. M.; Tkatchenko, A., Understanding the Role of Vibrations, Exact Exchange, and Many-body van der Waals Interactions in the Cohesive Properties of Molecular Crystals. J. Chem. Phys. 2013, 139, 024705. (56) Reilly, A. M.; Tkatchenko, A., Van der Waals Dispersion Interactions in Molecular Materials: Beyond Pairwise Additivity. Chem. Sci. 2015, 6, 3289-3301. (57) Červinka, C.; Beran, G. J. O., Ab initio thermodynamic properties and their uncertainties for crystalline alpha-methanol. Phys. Chem. Chem. Phys. 2017, 19, 29940-29953. (58) Červinka, C.; Beran, G. J. O., Ab initio prediction of the polymorph phase diagram for crystalline methanol. Chem. Sci. 2018, 9, 4622-4629. (59) McKinley, J. L.; Beran, G. J. O., Identifying pragmatic quasi-harmonic electronic structure approaches for modeling molecular crystal thermal expansion. Faraday Discuss. 2018, 211, 181-207. (60) Stoffel, R. P.; Wessel, C.; Lumey, M.-W.; Dronskowski, R., Ab Initio Thermochemistry of Solid-State Materials. Angew. Chem. Int. Edit. 2010, 49, 5242-5266. (61) Červinka, C.; Fulem, M., Probing the Accuracy of First-Principles Modeling of Molecular Crystals: Calculation of Sublimation Pressures. Cryst. Growth Des. 2019, 19, 808-820. (62) Červinka, C.; Fulem, M.; Stoffel, R. P.; Dronskowski, R., Thermodynamic Properties of Molecular Crystals Calculated within the Quasi-Harmonic Approximation. J. Phys. Chem. A 2016, 120, 2022-2034. (63) Heit, Y. N.; Beran, G. J. O., How Important Is Thermal Expansion for Predicting Molecular Crystal Structures and Thermochemistry at Finite Temperatures? Acta Cryst. B 2016, 72, 514-529. (64) Heit, Y. N.; Nanda, K. D.; Beran, G. J. O., Predicting Finite-Temperature Properties of Crystalline Carbon Dioxide from First Principles with Quantitative Accuracy. Chem. Sci. 2016, 7, 246-255. (65) Hirata, S.; Gilliard, K.; He, X.; Li, J.; Sode, O., Ab Initio Molecular Crystal Structures, Spectra, and Phase Diagrams. Acc. Chem. Res. 2014, 47, 2721-2730. (66) George, J.; Deringer, V. L.; Wang, A.; Muller, P.; Englert, U.; Dronskowski, R., Lattice Thermal Expansion and Anisotropic Displacements in Alpha-Sulfur from Diffraction Experiments and First-Principles Theory. J. Chem. Phys. 2016, 145, 234512. (67) Brandenburg, J. G.; Potticary, J.; Sparkes, H. A.; Price, S. L.; Hall, S. R., Thermal Expansion of Carbamazepine: Systematic Crystallographic Measurements Challenge Quantum Chemical Calculations. J. Phys. Chem. Lett. 2017, 8, 4319-4324. (68) Buchholz, H. K.; Hylton, R. K.; Brandenburg, J. G.; Seidel-Morgenstern, A.; Lorenz, H.; Stein, M.; Price, S. L., Thermochemistry of Racemic and Enantiopure Organic Crystals for Predicting Enantiomer Separation. Cryst. Growth Des. 2017, 17, 4676-4686. (69) Irikura, K. K.; Frurip, D. J., Computational Thermochemistry: Prediction and Estimation of Molecular Thermodynamics. American Chemical Society: Washington, DC, 1998. (70) Pfaendtner, J.; Yu, X.; Broadbelt, L. J., The 1-D Hindered Rotor Approximation. Theor. Chem. Acc. 2007, 118, 881-898. 41



(71) Červinka, C.; Fulem, M.; Růžička, K., Evaluation of Accuracy of Ideal-Gas Heat Capacity and Entropy Calculations by Density Functional Theory (DFT) for Rigid Molecules. J. Chem. Eng. Data 2012, 57, 227-232. (72) Červinka, C.; Fulem, M.; Růžička, K., Evaluation of Uncertainty of Ideal-Gas Entropy and Heat Capacity Calculations by Density Functional Theory (DFT) for Molecules Containing Symmetrical Internal Rotors. J. Chem. Eng. Data 2013, 58, 1382-1390. (73) Paulechka, Y. U.; Kabo, G. J.; Emel'yanenko, V. N., Structure, Conformations, Vibrations, and Ideal-Gas Properties of 1-Alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide Ionic Pairs and Constituent Ions. J. Phys. Chem. B 2008, 112, 15708-15717. (74) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (75) Grimme, S.; Ehrlich, S.; Goerigk, L., Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456-1465. (76) Blöchl, P. E., Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953-17979. (77) Kresse, G.; Joubert, D., From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method. Phys. Rev. B 1999, 59, 1758-1775. (78) Monkhorst, H. J.; Pack, J. D., Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188-5192. (79) Hafner, J.; Kresse, G.; Vogtenhuber, D.; Marsman, M. Vienna Ab-initio Simulation Package 5.4.1, 5.4.1; 2014. (80) Matsumoto, K.; Hagiwara, R.; Mazej, Z.; Benkic, P.; Zemva, B., Crystal structures of frozen room temperature ionic liquids, 1-ethyl-3-methylimidazolium tetrafluoroborate (EMImBF4), hexafluoroniobate (EMImNbF6) and hexafluorotantalate (EMImTaF6), determined by low-temperature X-ray diffraction. Solid State Sci. 2006, 8, 1250-1257. (81) Reichert, W. M.; Holbrey, J. D.; Swatloski, R. P.; Gutowski, K. E.; Visser, A. E.; Nieuwenhuyzen, M.; Seddon, K. R.; Rogers, R. D., Solid-state analysis of low-melting 1,3dialkylimidazolium hexafluorophosphate salts (ionic liquids) by combined x-ray crystallographic and computational analyses. Cryst. Growth Des. 2007, 7, 1106-1114. (82) Paulechka, Y. U.; Kabo, G. J.; Blokhin, A. V.; Shaplov, A. S.; Lozinskaya, E. I.; Golovanov, D. G.; Lyssenko, K. A.; Korlyukov, A. A.; Vygodskii, Y. S., IR and X-ray Study of Polymorphism in 1-Alkyl-3-methylimidazolium Bis(trifluoromethanesulfonyl)imides. J. Phys. Chem. B 2009, 113, 9538-9546. (83) Parlinski, K.; Li, Z. Q.; Kawazoe, Y., First-Principles Determination of the Soft Mode in Cubic ZrO2. Phys. Rev. Lett. 1997, 78, 4063-4066. (84) Togo, A.; Tanaka, I., First Principles Phonon Calculations in Materials Science. Scr. Mater. 2015, 108, 1-5. (85) Murnaghan, F. D., The Compressibility of Media under Extreme Pressures. Proc. Natl. Acad. Sci. U.S.A. 1944, 30, 244-247. (86) Heit, Y.; Beran, G. J. O., Exploiting Space-Group Symmetry in Fragment-Based Molecular Crystal Calculations. J. Comput. Chem. 2014, 35, 2205–2214. (87) Werner, H. J.; Manby, F. R., Explicitly correlated second-order perturbation theory using density fitting and local approximations. J. Chem. Phys. 2006, 124, 054114. (88) Werner, H.-J.; Adler, T. B.; Manby, F. R., General orbital invariant MP2-F12 theory. J. Chem. Phys. 2007, 126, 164102. (89) Pitoňák, M.; Heßelmann, A., Accurate Intermolecular Interaction Energies from a Combination of MP2 and TDDFT Response Theory. J. Chem. Theory Comput. 2010, 6, 168-178. 42


Page 42 of 46

Page 43 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(90) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M., Molpro: A GeneralPurpose Quantum Chemistry Program Package. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 242-253. (91) Burns, L. A.; Marshall, M. S.; Sherrill, C. D., Appointing silver and bronze standards for noncovalent interactions: A comparison of spin-component-scaled (SCS), explicitly correlated (F12), and specialized wavefunction approaches. J. Chem. Phys. 2014, 141, 234111. (92) Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C. M.; Civalleri, B.; Maschio, L.; Rérat, M.; Casassa, S.; Baima, J.; Salustro, S.; Kirtman, B., Quantum-mechanical condensed matter simulations with CRYSTAL. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2018, 8, e1360. (93) Peintinger, M. F.; Oliveira, D. V.; Bredow, T., Consistent gaussian basis sets of TripleZeta valence with polarization quality for solid-State Calculations. J. Comput. Chem. 2013, 34, 451-459. (94) East, A. L. L.; Radom, L., Ab Initio Statistical Thermodynamical Models for the Computation of Third-Law Entropies. J. Chem. Phys. 1997, 106, 6655-6674. (95) Pitzer, K. S., Energy Levels and Thermodynamic Functions for Molecules with Internal Rotation .2. Unsymmetrical Tops Attached to a Rigid Frame. J. Chem. Phys. 1946, 14, 239-243. (96) Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1-19. (97) Canongia Lopes, J. N.; Deschamps, J.; Padua, A. A. H., Modeling ionic liquids using a systematic all-atom force field. J Phys. Chem. B 2004, 108, 2038-2047. (98) Canongia Lopes, J. N.; Padua, A. A. H., Molecular force field for ionic liquids composed of triflate or bistriflylimide anions. J Phys. Chem. B 2004, 108, 16893-16898. (99) Kaminski, G. A.; Jorgensen, W. L., Host-guest chemistry of rotaxanes and catenanes: application of a polarizable all-atom force field to cyclobis(paraquat-p-phenylene) complexes with disubstituted benzenes and biphenyls. J. Chem. Soc., Perkin Trans. 2 1999, 11, 2365-2375. (100) Paulechka, Y. U.; Blokhin, A. V.; Kabo, G. J.; Strechan, A. A., Thermodynamic properties and polymorphism of 1-alkyl-3-methylimidazolium bis(triflamides). J. Chem. Thermodyn. 2007, 39, 866-877. (101) Marsh, K. N.; Deev, A.; Wu, A. C. T.; Tran, E.; Klamt, A., Room temperature ionic liquids as replacements for conventional solvents – A review. Korean J. Chem. Eng. 2002, 19, 357-362. (102) Valkenburg, M. E. V.; Vaughn, R. L.; Williams, M.; Wilkes, J. S., Thermochemistry of ionic liquid heat-transfer fluids. Thermochim. Acta 2005, 425, 181-188. (103) Ferreira, A. F.; Simões, P. N.; Ferreira, A. G. M., Quaternary phosphonium-based ionic liquids: Thermal stability and heat capacity of the liquid phase. J. Chem. Thermodyn. 2012, 45, 16-27. (104) Waliszewski, D.; Stępniak, I.; Piekarski, H.; Lewandowski, A., Heat capacities of ionic liquids and their heats of solution in molecular liquids. Thermochim. Acta 2005, 433, 149-152. (105) Serra, P. B. P.; Ribeiro, F. M. S.; Rocha, M. A. A.; Fulem, M.; Růžička, K.; Coutinho, J. A. P.; Santos, L. M. N. B. F., Solid-liquid equilibrium and heat capacity trend in the alkylimidazolium PF6 series. Journal of Molecular Liquids 2017, 248, 678-687. (106) Ahrenberg, M.; Beck, M.; Neise, C.; Kessler, O.; Kragl, U.; Verevkin, S. P.; Schick, C., Vapor pressure of ionic liquids at low temperatures from AC-chip-calorimetry. Phys. Chem. Chem. Phys. 2016, 18, 21381-21390. (107) Ahrenberg, M.; Brinckmann, M.; Schmelzer, J. W.; Beck, M.; Schmidt, C.; Kessler, O.; Kragl, U.; Verevkin, S. P.; Schick, C., Determination of volatility of ionic liquids at the nanoscale by means of ultra-fast scanning calorimetry. Phys. Chem. Chem. Phys. 2014, 16, 29712980. 43



(108) Emel'yanenko, V. N.; Verevkin, S. P.; Heintz, A., The gaseous enthalpy of formation of the ionic liquid 1-butyl-3-methylimidazolium dicyanamide from combustion calorimetry, vapor pressure measurements, and ab initio calculations. J. Am. Chem. Soc. 2007, 129, 3930-3937. (109) Verevkin, S. P.; Ralys, R. V.; Zaitsau, D. H.; Emel’yanenko, V. N.; Schick, C., Express thermo-gravimetric method for the vaporization enthalpies appraisal for very low volatile molecular and ionic compounds. Thermochim. Acta 2012, 538, 55-62. (110) Verevkin, S. P.; Zaitsau, D. H.; Emelyanenko, V. N.; Heintz, A., A new method for the determination of vaporization enthalpies of ionic liquids at low temperatures. J Phys. Chem. B 2011, 115, 12889-12895. (111) Verevkin, S. P.; Zaitsau, D. H.; Emel'yanenko, V. N.; Yermalayeu, A. V.; Schick, C.; Liu, H.; Maginn, E. J.; Bulut, S.; Krossing, I.; Kalb, R., Making sense of enthalpy of vaporization trends for ionic liquids: new experimental and simulation data show a simple linear relationship and help reconcile previous data. J Phys. Chem. B 2013, 117, 6473-6486. (112) Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A.; Paulechka, Y. U.; Tschersich, A.; Verevkin, S. P.; Heintz, A., Experimental vapor pressures of 1-alkyl-3-methylimidazolium Bis(trifluoromethylsulfonyl)imides and a correlation scheme for estimation of vaporization enthalpies of ionic liquids. J. Phys. Chem. A 2006, 110, 7303-7306. (113) Rocha, M. A.; Lima, C. F.; Gomes, L. R.; Schroder, B.; Coutinho, J. A.; Marrucho, I. M.; Esperanca, J. M.; Rebelo, L. P.; Shimizu, K.; Lopes, J. N.; Santos, L. M., High-accuracy vapor pressure data of the extended [C(n)C1im][Ntf2] ionic liquid series: trend changes and structural shifts. J Phys. Chem. B 2011, 115, 10919-10926. (114) Santos, L. M. N. B. F.; Ferreira, A. I. M. C. L.; Štejfa, V.; Rodrigues, A. S. M. C.; Rocha, M. A. A.; Torres, M. C.; Tavares, F. M. S.; Carpinteiro, F. S., Development of the Knudsen effusion methodology for vapour pressure measurements of low volatile liquids and solids based on a quartz crystal microbalance. J. Chem. Thermodyn. 2018, 126, 171-186. (115) Lovelock, K. R.; Deyko, A.; Licence, P.; Jones, R. G., Vaporisation of an ionic liquid near room temperature. Phys. Chem. Chem. Phys. 2010, 12, 8893-8901. (116) Růžička, K.; Majer, V., Simultaneous Treatment of Vapor Pressures and Related Thermal Data Between the Triple and Normal Boiling Temperatures for n-Alkanes C5-C20. J. Phys. Chem. Ref. Data 1994, 23, 1-39. (117) Clarke, E. C. W.; Glew, D. N., Evaluation of thermodynamic functions from equilibrium constants. Trans. Faraday Soc. 1966, 62, 539-547. (118) Heym, F.; Korth, W.; Etzold, B. J. M.; Kern, C.; Jess, A., Determination of vapor pressure and thermal decomposition using thermogravimetrical analysis. Thermochim. Acta 2015, 622, 917. (119) Armstrong, J. P.; Hurst, C.; Jones, R. G.; Licence, P.; Lovelock, K. R.; Satterley, C. J.; Villar-Garcia, I. J., Vapourisation of ionic liquids. Phys Chem Chem Phys 2007, 9, 982-90. (120) Tolstoguzov, A. B., Investigation of thermal evaporation of imidazolium-based ionic liquids. Mass-Spektrom. 2007, 4, 283-288. (121) Luo, H.; Baker, G. A.; Dai, S., Isothermogravimetric Determination of the Enthalpies of Vaporization of 1-Alkyl-3-methylimidazolium Ionic Liquids. J. Phys. Chem. B 2008, 112, 10077-10081. (122) Wang, C.; Luo, H.; Li, H.; Dai, S., Direct UV-spectroscopic measurement of selected ionic-liquid vapors. Phys. Chem. Chem. Phys. 2010, 12, 7246-7250. (123) Dunaev, A. M.; Motalov, V. B.; Kudin, L. S.; Butman, M. F., Molecular and ionic composition of saturated vapor over EMImNTf 2 ionic liquid. J. Mol. Liq. 2016, 219, 599-601. (124) M. S. S. Esperança, J.; Canongia Lopes, J. N.; Tariq, M.; Santos, L. M. N. B. F.; Magee, J. W.; Rebelo, L. P. N., Volatility of Aprotic Ionic Liquids — A Review. J. Chem. Eng. Data 2010, 55, 3-12. 44


Page 44 of 46

Page 45 of 46 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(125) Atkins, P. W., Shriver & Atkins' inorganic chemistry. Oxford University Press: Oxford; New York, 2010. (126) Dong, K.; Zhao, L. D.; Wang, Q.; Song, Y. T.; Zhang, S. J., Are ionic liquids pairwise in gas phase? A cluster approach and in situ IR study. Phys. Chem. Chem. Phys. 2013, 15, 60346040. (127) Kirchner, B.; di Dio, P. J.; Hutter, J., Real-World Predictions from Ab Initio Molecular Dynamics Simulations. In Multiscale Molecular Methods in Applied Chemistry, 2012; pp 109153. (128) Izgorodina, E. I.; Seeger, Z. L.; Scarborough, D. L. A.; Tan, S. Y. S., Quantum Chemical Methods for the Prediction of Energetic, Physical, and Spectroscopic Properties of Ionic Liquids. Chem. Rev. 2017, 117, 6696-6754. (129) LeBlanc, L. M.; Dale, S. G.; Taylor, C. R.; Becke, A. D.; Day, G. M.; Johnson, E. R., Pervasive Delocalisation Error Causes Spurious Proton Transfer in Organic Acid–Base Co‐Crystals. Angew. Chem. Int. Ed. 2018, 57, 14906-14910. (130) Izgorodina, E. I.; Bernard, U. L.; MacFarlane, D. R., Ion-Pair Binding Energies of Ionic Liquids: Can DFT Compete with Ab Initio-Based Methods? J. Phys. Chem. A 2009, 113, 70647072. (131) Zahn, S.; MacFarlane, D. R.; Izgorodina, E. I., Assessment of Kohn–Sham density functional theory and Møller–Plesset perturbation theory for ionic liquids. Phys. Chem. Chem. Phys. 2013, 15, 13664-13675. (132) Červinka, C.; Fulem, M.; Štejfa, V.; Růžička, K., Analysis of uncertainty in calculation of ideal-gas thermodynamic properties using one-dimensional hindered rotor (1-DHR) model. J. Chem. Eng. Data 2017, 62, 445-455. (133) Taguchi, R.; Machida, H.; Sato, Y.; Smith, R. L., High-Pressure Densities of 1-Alkyl-3methylimidazolium Hexafluorophosphates and 1-Alkyl-3-methylimidazolium Tetrafluoroborates at Temperatures from (313 to 473) K and at Pressures up to 200 MPa. J. Chem. Eng. Data 2009, 54, 22-27. (134) Serra, P. B. P.; Ribeiro, F. M. S.; Rocha, M. A. A.; Fulem, M.; Růžička, K.; Coutinho, J. A. P.; Santos, L. M. N. B. F., Solid-liquid equilibrium and heat capacity trend in the alkylimidazolium PF6 series. J. Mol. Liq. 2017, 248, 678-687. (135) Nishida, T.; Tashiro, Y.; Yamamoto, M., Physical and electrochemical properties of 1alkyl-3-methylimidazolium tetrafluoroborate for electrolyte. J. Fluorine Chem. 2003, 120, 135141. (136) Shirota, H.; Mandai, T.; Fukazawa, H.; Kato, T., Comparison between Dicationic and Monocationic Ionic Liquids: Liquid Density, Thermal Properties, Surface Tension, and Shear Viscosity. J. Chem. Eng. Data 2011, 56, 2453-2459. (137) Wachter, P.; Schreiner, C.; Schweiger, H.-G.; Gores, H. J., Determination of phase transition points of ionic liquids by combination of thermal analysis and conductivity measurements at very low heating and cooling rates. J. Chem. Thermodyn. 2010, 42, 900-903. (138) Domańska, U.; Marciniak, A., Solubility of 1-Alkyl-3-methylimidazolium Hexafluorophosphate in Hydrocarbons. J. Chem. Eng. Data 2003, 48, 451-456. (139) Sifaoui, H.; Ait-Kaci, A.; Modarressi, A.; Rogalski, M., Solid–liquid equilibria of three binary systems: 1-Ethyl-3-methylimidazolium hexafluorophosphate+2-phenylimidazole, or 4,5diphenylimidazole or 2,4,5-triphenylimidazole. Thermochim. Acta 2007, 456, 114-119. (140) Červinka, C.; Beran, G. J. O., Towards reliable ab initio sublimation pressures for organic molecular crystals – are we there yet? Phys. Chem. Chem. Phys. 2019, 21, 14799-14810.

45



Table of contents

46


Page 46 of 46

Cohesive Properties of Ionic Liquids Calculated from First Principles

Recommend Documents