Density Functional Computations and Molecular Dynamics

Subscriber access provided by FLORIDA STATE UNIV

Article

Density Functional Computations and Molecular Dynamics Simulations of the Triethylammonium Triflate Protic Ionic Liquid Juan F. Mora Cardozo, Tatsiana Burankova, Jan Peter Embs, Antonio Benedetto, and Pietro A. Ballone J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10373 • Publication Date (Web): 29 Nov 2017 Downloaded from http://pubs.acs.org on December 3, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 34 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

The Journal of Physical Chemistry

Density Functional Computations and Molecular Dynamics Simulations of the Triethylammonium Triflate Protic Ionic Liquid Juan F. Mora Cardozo(1) , T. Burankova(1) , J. P. Embs(1) , A. Benedetto(2,3) , and P. Ballone(4),∗ (1) Laboratory for Neutron Scattering and Imaging, Paul Scherrer Institute, Villigen PSI, Villigen 5232, Switzerland (2) School of Chemistry, University College Dublin, Belfield, Dublin 4, Ireland (3) School of Physics, University College Dublin, Belfield, Dublin 4, Ireland and (4) Italian Institute of Technology, Via Morego 30, 16163 Genova, Italy

Abstract Systematic molecular dynamics simulations based on an empirical force field have been carried out for samples of triethylammonium trifluoromethanesulfonate (triethylammonium triflate, [TEA][Tf]), covering a wide temperature range 200 K ≤ T ≤ 400 K and analysing a broad set of properties, from self-diffusion and electrical conductivity to rotational relaxation and hydrogen-bond dynamics. The study is motivated by recent quasielastic neutron scattering and differential scanning calorimetry measurements on the same system, revealing two successive first order transitions at T ∼ 230K and T ∼ 310K (on heating), as well as an intriguing and partly unexplained variety of sub-diffusive motions of the acidic proton. Simulations show a weakly discontinuous transition at T = 310K, and highlights an anomaly at T = 260 K in the rotational relaxation of ions, that we identify with the simulation analogue of the experimental transition at T = 230K. Thus, simulations help identifying the nature of the experimental transitions, confirming that the highest temperature one corresponds to melting, while the one taking place at lower T is a transition from the crystal, stable at T ≤ 260 K, to a plastic phase (260 ≤ T ≤ 310 K), in which molecules are able to rotate without diffusing. Rotations, in particular, account for the sub-diffusive motion seen at intermediate T both in the experiments and in the simulation. The structure, distribution and strength of hydrogen bonds are investigated by molecular dynamics and by density functional computations. Clustering of ions of the same sign, and the effect of contamination by water at 1 % wgt concentration are discussed as well.

∗

Corresponding Author: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

I.

Page 2 of 34

INTRODUCTION

The reaction of a Brønsted acid with a Brønsted base normally results in the transfer of a proton from the acid to the base, leaving behind a cation-anion pair. Compounds of this type whose melting temperature is below 100 ◦ C are rightful members of the vast family of room temperature ionic liquids (RTILs, 1), of which they represent the protic ionic liquid (PIL) sub-class.2,3 The relevant impact of hydrogen bonding (HB) on structure, diffusion, viscosity and electrical conductivity4 often prevents the usage of PILs in the most popular applications of RTILs as solvents and lubricants. On the other hand, the sharable proton of PILs underlies many other applications in catalysis and especially in electrochemistry, as exemplified by their usage as anhydrous ionic conductors for fuel cells,5–7 in super-capacitors,8–10 in lithium-ion batteries11 and in solar cells.12 In this vast application context, PILs generally enjoy the advantage of being simple to prepare, and cheaper than several other electrolytes. To fulfill all PIL’s promises, however, much development is still needed. Most of it will consist primarily of very applied research and careful optimization, but it is apparent that also genuine basic science has an important role to play. For instance, the microscopic details of hydrogen bonding and proton conductivity are still partly unknown, and a broad scope investigation of these properties could greatly ease further progress in PILs applications. These considerations are illustrated and exemplified by the recent characterisation of dynamical properties of the triethylammonium-trifluoromethanesulfonate (triethylammonium-triflate, [TEA][Tf]) by elastic and quasi-elastic neutron scattering, complemented by the analysis of thermal properties by differential scanning calorimetry (DSC).13 Both natural (i.e., hydrogenated) and partially deuterated samples were used, greatly enhancing the ability of neutron scattering to focus on specific dynamical modes. The experimental study revealed two successive transitions, at temperatures that are made somewhat uncertain by sizeable hysteresis. On heating from the powder, the first transition takes place at T = 230 K, the second one at T ∼ 310 K. On cooling, the system transforms at T = 260 K and at T = 210 K. Diffusion coefficients have been estimated by neutron scattering, and compared to those from pulsed field gradient NMR,14,15 highlighting expected systematic differences due to the different time scale of the PFG-NMR and neutron scattering measurements. The quantitative difference extends to the qualitative description of self-diffusion as an activated Arrhenius process (neutron scattering13 ) or as a Vogel-Fulcher-Tamman process (PGE-NMR16 ). Last but not least, the analysis of scattering intensities revealed a variety of wide amplitude, sub-diffusive motions, whose precise nature

2


Page 3 of 34 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


cannot be unambiguously determined from the available experimental data. The present study has been designed to investigate the same system from a different point of view, using computational approaches. Molecular dynamics based on empirical force fields, in particular, has been applied to analyse a variety of properties ranging from thermodynamic functions to diffusion coefficients and electrical conductivity, to the structure, distribution and dynamics of HBs throughout the system. The investigation has covered a wide temperature range (200 K ≤ T ≤ 400 K), largely overlapping with the experimental one (200 K ≤ T ≤ 440 K). Since, to the best of our knowledge, the crystal structure of [TEA][Tf] is not available, our simulations concern only the cooling branch of the thermal cycle (heating from powder and cooling from the liquid) considered in the experimental paper. Plots of thermodynamic functions such as the average potential energy ad volume per ion pair reveal a phase change at T = 310 K, that manifests itself as a weakly first order transition, accompanied by the freezing of diffusion and electrical conductivity. At variance from experiments, no other transition is apparent in the thermodynamic functions computed by simulation. At T ∼ 260 K, however, simulations display a marked anomaly in the time correlation function of the orientation of both ion types. Taking into account the large difference of characteristic time scales of simulations and of experiments, we identify the continuous transition at T = 260K with the (discontinuous) experimental one occurring at T = 230 K on heating, and at T = 210 on cooling. Simulations, in addition, provide a wealth of information on the thermodynamic, structural and dynamical properties of the system over a wide temperature range. These analyses allow us to attribute the first transition transition to melting at T ∼ 310K, while T = 260 K marks the transformation between the solid (at T ≤ 260 K) and a plastic phase (at T > 260 K) in which [TEA]+ and [Tf]− ions rotate without diffusing. Molecular rotations could also account for the sub-diffusive motions seen at intermediate temperatures both in experimental and in simulation results. The computation of linear dynamical coefficients shows that in our model diffusion is an activated Arrhenius process, with an activation free energy of the order of 10 kJ/mol at low temperature, and 40 kJ/mol at high temperature. Comparison of the computed electrical conductivity with the prediction of the NernstEinstein relation reveals a significant pairing of cations and anions, that only at the highest simulation temperature (T = 400 K) approaches a predominantly dissociated state. The analysis of simulation trajectories shows that hydrogen bonding is not dramatically affected by temperature in the 200 K ≤ T ≤ 400 K range. Hydrogen bonding, on the other hand, gives origin to intriguing complexes made of up to three cations and a single anion, whose association is likely to affect the ionic and molecular mobility.17–19 3



Page 4 of 34

This unusual structural feature, together with many other observations made upon analysing MD trajectories, raised a number of questions that we address by ab-initio (density functional) computations. To complete our analysis, we present simulations for samples contaminated by water at low concentration (1 % in weight). The results show that this addition does not affect significantly thermodynamic properties, ruling out the possibility that the two transitions seen in experiments correspond to the same transformation in dry and in hydrated domains of the nominally pure [TEA][Tf] sample. On the other hand, even the low amount of water molecules added to the simulated system significantly affects the system dynamics, increasing dynamical coefficients and decreasing relaxation times. In recent years, [TEA][Tf] and closely related ammonium-based protic ionic liquids have been simulated by other groups, see, for instance, Ref. 20,21, using approaches similar to ours. To a large extent, the result from the different studies are compatible with each other. Our study, however, is the first one to analyse in detail the temperature evolution of a broad range of properties, giving us the opportunity to reproduce and identify the phase transitions seen in experiments.13

II.

METHOD

The bulk of our computations consisted of molecular dynamics simulations based on an empirical force field for the [TEA][Tf] ion pair, whose schematic structure is shown in Figure 1 The force field is of the all-atom OPLS/Amber type.22,23 The starting point for its parametrisation was provided by the potential of Ref. 24 for [Tf]− (see also Ref. 25), while [TEA]+ can be described using the extensive tabulation of force field parameters given in Ref. 26. Fine tuning of the potential, and of the [TEA]+ , [Tf]− cross interaction has been carried out using the results of our density functional computations (see Sec. IV) for the ground state geometry and harmonic dynamics of the [TEA][Tf] molecule. Both vibrational eigenvalues and eigenvectors have been used to this effect. Atomic charges have also been verified and re-tuned using the electrostatic potential (ESP)27,28 values provided by our density functional computations for the isolated ions and for the neutral [TEA][Tf] pair. The full list of parameters is reported in Table S1 to Table S4 of S.I. Simple chemical considerations suggest that no atom in [TEA][Tf] has the simultaneous H-donor and H-acceptor ability needed for the onset of a Grotthus proton conduction channel. This observation is the major justification for our usage of a fixed-topology force field model, avoiding the complications of a reactive force field approach.29

4


Page 5 of 34 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: Schematic structure of [TEA]+ and [Tf]− .

All classical simulations with the empirical force field have been carried out in the isothermalisobaric (NPT) ensemble using the Gromacs package.30 Constant temperature was enforced using the Nose’-Hoover thermostat with a characteristic relaxation time of 2 ps, while the constant pressure of 1 atm was enforced using the Parrinello-Rahman extended Lagrangian formulation, again with a characteristic time of 2 ps. All samples are periodic in space, and the long range Coulomb potential and forces are computed using particle-particle-particle-mesh Ewald sums.31 Bond distances, bending and, of course, dihedral angles are all flexible. Simulations have been carried out for dry [TEA][Tf], and for samples slightly contaminated by water. The TIP4P model was used for water, whose cross interactions with [TEA]+ and [Tf]− has been estimated using Berthelot’s rule. Our force field and molecular dynamics set up are equivalent to those of Ref. 20,21, despite a few different choices. The most important one is the rescaling of ion charges, summing up to some 80% of the formal value, adopted both in Ref. 20 and in Ref. 21, but not used here. Rescaling has been shown to improve the value of linear dynamical coefficient compared to experiments,32 but, despite some intuitive justification, it is still an ad-hoc adjustment. Moreover, we did not rescale charges because the equilibrium density of the simulated samples already slightly underestimates the experimental value (see following section). Remarkably, the underestimation of transport coefficients by our model does not seem to affect the location and quality of the transitions that represent the real interest of our study. Our level of ab-initio corresponds to density functional theory, with the Perdew-Becke-Ernzerhof (PBE) approximation for the exchange and correlation energy.33 Computations have been carried out using primarily the Gaussian 09 package,34 which uses Gaussian basis functions, and, as such, is suitable for all-electron computations. In all Gaussian 09 computations, a 6-311++G(2d,2p) basis has been used. Gaussian 09 has been used also to check total energies and especially energy differences among the different species using the B3LYP approximation35 for the exchange correlation function. Lacking a precise experimental thermo-chemical benchmark, the results of BPE and B3LYP turn out to be broadly equivalent. Other computations have been carried out using a different computer package,

5



Page 6 of 34

i.e., CPMD,36 based on plane waves and soft pseudopotentials.37 CPMD is especially optimised for molecular dynamics simulations, that have been extensively used in our computations. Ground state geometries were optimised by the so called eigenvector following method38 in Gaussian 09, and by simulated annealing39 in CPMD. This last approach, more intuitive than the first one, is, however, computationally time consuming. Dispersion interactions have been included at the empirical level40 in our CPMD computations, and neglected in the Gaussian 09 computations, upon verifying by CPMD that their exclusion did not change significantly the equilibrium geometry nor the cohesive properties. By and large, and despite the slight differences in their approach, the results of CPMD and Gaussian 09 are mutually compatible whenever a comparison is possible.

III. A.

CLASSICAL MD SIMULATIONS: RESULTS Dry [TEA][Tf ] samples

Simulations have been carried out for samples of 125 ion pairs, enclosed in a cubic box periodically repeated in space. The sample size, although relatively small, is sufficient to investigate properties of homogeneous systems. At the same time, the limited sample size allows us to extend simulations to the long time (∼ 100 ns) required to investigate diffusion and electrical conductivity in a disordered system that exhibits long relaxation times up to fairly high temperature. In our study we first equilibrate our sample at T = 400K and P = 1 atm during 15 ns. The equilibrated sample is apparently liquid-like, and it is used to start a cooling stage down to T = 200K, with a cascade of equilibration runs of 5 ns at T = 380K, 360K, 340K, 320K, 310K, 300K, 290K, 280K, 270K, 260K, 240K, 220K, 200K. Each of the samples obtained in this way is propagated in time in 5ns blocks repeated until equilibration has been achieved, and sufficient statistics has been accumulated. The change from equilibration to production runs is decided by monitoring the drift of diagnostic properties such as volume and potential energy. Needless to say, equilibration is longer at low T than at high T , and longest at the intermediate temperatures around T = 300 K, where thermodynamic and dynamical anomalies manifest themselves. All simulations are carried out with a time step of 1 fs, and configurations are stored on the trajectory file every 5 ps. Shorter runs (from 300 ps to 900 ps) with more frequent sampling (one configuration every 0.1 ps) have been performed to monitor fast processes such as the re-orientation of H bonds. In all cases, statistics runs lasted at least 40 ns. At each temperature below T = 300K,

6


Page 7 of 34 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


the sum of equilibration and production stages lasted more than 100ns. We emphasise that all the results presented below concern the cooling of the sample from the liquid phase. No heating simulations have been attempted because no crystal configuration at low T is available. Ab-initio searches for the crystalline ground state have been carried out (see Sec. III of S.I.) but the result is not sufficiently validated to justify extensive simulations of its melting. The average potential energy per ion pair hEi as a function of T is shown in panel(a) of Figure 2. The hEi(T ) relation is remarkably linear both at high and at low T , with however a weak anomaly at Ta = 310 ± 5 K (the a subscript stays for anomaly). The anomaly represents a quantitatively small but nevertheless clear drop in the average potential energy below its linear extrapolation (from high temperature), taking place at at 310K, accompanied by an equally slight change in the constantpressure specific heat Cp (T ), that on cooling across Ta decreases from Cp = 0.41 kJ/mol/K to Cp = 0.39 kJ/mol/K. To highlight the tiny anomaly in the average potential energy, in the inset of Figure 2 we plot the difference δE(T ) between hEi(T ) and the high-temperature linear interpolation α + βT of the average potential energy, showing that δE is virtually zero above T ∼ 310 K, and drops to about −0.8 kcal/mol below the transition temperature. The slight discontinuity in the average potential energy at Ta identifies a weakly first order transition. We anticipate that the analysis of dynamical transport coefficients shows that Ta also corresponds to the quenching of diffusion and electrical conductivity, thus suggesting that the transition at Ta = Tm = 310 K is in fact melting, joining a liquid and a solid-like phase. We emphasise that in this context solid-like is not equivalent to crystal-like. We also anticipate that the structure of the solid just below Tm is highly disordered, partly because of the sluggish dynamics and slow relaxation due to the molecular nature of the ions, but, even more, because the equilibrium phase resulting from cooling below Tm lacks the orientational ordering of a genuine crystal phase. The observation of the anomaly in hEi(T ) is remarkable. First order transitions are rarely observed on cooling molecular liquids by computer simulations, whose likely result is to turn the sample into a glass,41 even for systems less complex than [TEA][Tf] (see, for instance, the cooling branch in Figure 2, Ref. 42 for SiO2 ). A completely analogous behaviour is displayed by the average volume hV i(T ), shown in panel (b) of Figure 2. We note in passing that the density of the simulated sample is ρ(T ) = 1.34 g/cm3 at T = 200, and ρ(T ) = 1.20 g/cm3 at T = 400K, compared to the experimental value in the liquid phase at 403K of ρ403K = 1.33 g/cm3 (Ref. 14) or ρ400K = 1.16 g/cm3 (Ref. 43). The underestimation of the density by our model is one major reason why we do not re-scale the charges to a value below the nominal one. 7



360

0

0.0

350 -0.5

-1.0 200

250

300

350

400

340

[K]

〈V〉

T

[A 3 ]

〈E〉 −α−βT

[kJ/mol]

30

[kJ/mol]

0.5

〈E〉(T)

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 8 of 34

330

-30 320

(a)

-60 200

250

300

T

350

(b) 400

310 200

250

300

T

[K]

350

400

[K]

Figure 2: (a) Average potential energy per ion pair as a function of temperature. The statistical error bar is of the order of 0.5 kJ/mol, smaller than the size of the dots. Green line: linear interpolation to the four lowest-T points. Red line: linear interpolation to the four points of highest T . Inset: deviation of hEi(T ) from the high-temperature linear interpolation hEihigh (T ) = α + βT . (b) Average volume per ion pair as a function of temperature. The error bar on the volume is ∼ 0.1 ˚ A3 , smaller than the size of the dots. Green and red lines defined as in panel (a). No clear sign of a second transition a T < 310 K is apparent in either hEi(T ) or hV i(T ). The microscopic structure has been characterised primarily through the radial distribution functions. Both ions are fairly large and somewhat flexible, and we decided to represent cations by the coordinates of their nitrogen atom, and to look at correlations among N-O, N-N and O-O. The results are shown in Figure 3. By and large, the set of three distribution function are recognisable as the radial distribution functions of an ionic system, with the N-O (cation-anion) component displaying the strongest correlation (highest peak), in opposition of phase with the N-N and O-O functions, reaching their maxima where gN O (r) has its minima. The molecular nature of the ions, however, gives origin to a variety of sharp details that are not found in the radial distribution functions of simpler ionic systems. The N-O radial distribution function, for instance, displays a remarkable peak at short range (r ∼ 2.72 ˚ A ) corresponding to the N-O distance in their tight hydrogen bonding. The running coordination number nN O (r) is defined as: Z r nN O (r) = 4πρk r2 gN O (r)dr

(1)

0

When computed with ρk ≡ ρO it gives the average number of oxygen atoms around each N atoms. When computed with ρk ≡ ρN = ρO /3 it gives the average number of nitrogen atoms around each O 8


Page 9 of 34 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


atoms. The running coordination nN O (r) has a plateau for 3 ≤ r ≤ 3.6 (inset of Figure 3), corresponding to the deep minimum in gN O (r) following the H-bonding peak. When computed with ρk = ρO , the plateau value measures the average number of NH terminations engaged in a hydrogen bond. At T = 300 K this value is nearly exactly nN O = 1.00, and its dependence on T is mild. The nearly full saturation of all H-bonding sites and the weak dependence of the nN O plateau on T emphasise the relevance of the H-bonding mechanism. The secondary peak on the side of the main N-O peak at r = 5.22 ˚ A is a geometric consequence of the presence of three oxygen atoms on the -SO3 side of [Tf]− . On a slightly longer length scale and at low temperature (T ≤ 310), an inflection point in nN O (r) at r = 5.4 ˚ A marks a second weaker shell closing (see Figure S2 in SI), corresponding to the minimum in gN O (r) following the main broad peak at r ∼ 4.7 ˚ A . The nN O coordination of nearly 8 at the √ shell closing, together with the next-nearest neighbour peak at r = 2 × 4.7 ∼ 7 ˚ A suggests an NaCltype geometry for the distribution of cations and anions in space. This structural hint is supported by visual inspection of simulation snapshots and by the density functional computations reported in Sec. IV. Of the three radial functions displayed in Figure 3, gN N (r) displays the widest correlation hole, partly because of the like-charge repulsion among the N-H termination of [TEA]+ , and, more importantly, because of the N-H location at the centre of a tightly packed crown of ethyl chains. The O-O radial distribution function displays a prominent intra-molecular peak at r ∼ 2.63 ˚ A , that has been removed from Figure 3 because of its partial overlap with the short range N-O peak attributed to H-bonding. All these major structural features display a relatively weak temperature dependence from T = 400 K down to T = 200 K, see Figure S3 in SI, although below T = 310 K sharper feature start to develop especially in gN O (r), pointing to the propagation of partial ordering to the medium length scales covered by simulation.44 Nevertheless, all g(r)’s look closer to those of glassy systems than to those of crystalline samples. Weaker features in the radial distribution functions are more sensitive to temperature. For instance, down to T ∼ 310K, the N-F radial distribution function does not display any clear peak at short range that could be attributed to a N-F hydrogen bond, however weak. Below T = 300 K such a sharp peak starts to develop at rN F = 4.7 ˚ A , becoming apparent and prominent at T = 200K (see Figure S4 in S.I.). This feature, probably due to the low negative charge on the F atoms, could point to a weak H-bond joining N-H—F atoms, and the analysis of triplet configurations of this type indeed finds a 9



Page 10 of 34

few candidates within a cut-off N-F distance of 5.2 ˚ A and N − H − − − F angle < 30◦ . The long bond distance, and the sensitivity to temperature, however, suggest that this bonding mechanism has little effect on structural and dynamical properties of the system. Analysis of all other gXY (r) shows no other sharp feature at H-bonding distances in the radial distribution functions, confirming that the N-H–O case, and the very weak N-H–F case at low T , are the only H-bonding possibilities in [TEA][Tf]. The structural counterpart of the radial distribution function closer to experimental data is the set of partial structure factors Sij (k), where {i, j} are atom-type indices. Of the multitude of possible {i, j} combinations in [TEA][Tf], we selected again the N-O, N-N and O-O pairs to characterise the arrangement of ions in reciprocal space. To identify the role of packing and of Coulomb interactions, we display in Figure S5 and Figure S6 of S.I. the density-density and the charge-charge combination of the three Bathia-Thornton structure factors:45,46 Sρρ (k) =

X

Sαβ (k)

(2)

Zα Zβ Sαβ (k)

(3)

α,β

SZZ (k) =

X α,β

where {α, β ≡ N, O}, Sαβ (k) are partial structure factors, with the fictitious charges Z = 1, and Z = −1/3 attributed to N and O, respectively, to identify them as the cation and the anion of a globally neutral system. Inspection of Sρρ (k) and SZZ (k) clearly shows the prominent role of Coulomb forces in determining the structure of the system. The vanishing of SZZ (k) at low k is due to the perfect screening of charge in this ionic system.47 The low value of Sρρ (k) in the same k → limit is due to the low compressibility of this ionic system not far from its triple point.45 More intriguing is the sizeable pre-peak displayed by Sρρ (k) at k ∼ 0.7 ˚ A

−1 ,

that highlights the formation of ion clusters, dynamically modulating the

density and temporarily breaking the system homogeneity. The third function SρZ , describing density-charge cross correlations, is less prominent than the first two functions, pointing to a rather weak coupling of density and charge fluctuations that is found in many molten salts systems. The computation of linear dynamical coefficients confirms the presence of slow relaxation processes at all simulated temperatures, and the glassy translational dynamics of the system below T = 310K. The diffusion constant of both ions has been computed from the time dependence of their mean 10


Page 11 of 34

3 6

nNO(r)

4

gNO (r), gNN (r), gOO (r)

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

2

0 2

3

4

5

r [A]

1

T=300K 0 0

4

8

r

12

16

[A]

Figure 3: Radial distribution function of N-O (blue line), N-N (green line) and O-O (red line) pairs in [TEA][Tf] at T = 300K, P = 1 atm. Inset: short range portion of the running coordination number nN O (r), obtained by integration of the radial distribution function (see text). square displacement, which is defined as: M SDα (t) =

1 X h|ri (t + t0 ) − ri (t0 )|2 it0 Nα

(4)

i∈α

where α labels the ion type. In our computation, in particular, α = 1 labels the coordinates of the N atom in [TEA]+ , while α = 2 identifies the coordinates of the three oxygen atoms in [Tf]− . The average over the reference t0 spans the entire time interval covered by the simulation (i.e., at least 40 ns minus the time t itself). Typical results for the {M SDα (t), α = 1, 2 } at T = 360 K, T = 300K and T = 240 K are shown in Figure S7 of S.I. The diffusion coefficient is computed from the Einstein relation: 1 M SDα (t) t→∞ 6t

Dα = lim

(5)

where the limiting slope is in fact approximated by the linear coefficient of the linear interpolation over the 10 ≤ t ≤ 12 ns range. The choice of the time interval is crucial, and in our computation it has been selected upon inspecting the M SDα (t) curves at all T (See again examples in Figure S7 in S.I.). Beyond 12 ns, long-time fluctuations in M SDα (t) appear at the lowest temperatures. At 11



Page 12 of 34

intermediate temperatures 260 K ≤ T ≤ 310 K, M SDα (t) displays a non-negligible curvature up to ∼ 8 ns. Remarkably, the onset of the linear behaviour is faster at T ≤ 260 K than at intermediate temperatures (see the T = 240 K M SDα (t) again in Figure S7 of S.I.). These observation will become relevant in the discussion of molecular rotations reported below. Because of the downward curvature of M SDα (t) on the ns time scale, moving the range of the fit to lower times increases, even significantly, the computational estimate of the diffusion constant. In a similar way, the diffusion coefficients estimated by neutron scattering are systematically higher than those measured by PFGNMR, since the time probed by the former approach is much shorter than for the latter. To be precise, the neutron scattering estimates should be referred to as short-time diffusion coefficients. Given this dependence of the diffusion coefficients on the measurement time scale, and because of the systematic underestimation of Dα by unpolarisable, formal-charge force field used in our simulations, we will not comment in detail on the quantitative correspondence of experimental and computational results. For our purposes, the analysis of diffusion scaling with temperature is more significant. The results for the diffusion constants D+ and D− of cations and anions are shown in Figure 4 in a form immediately amenable to an analysis of activation (free) energy estimated by an Arrhenius relation. To give an idea of the diffusion rate of ions in liquid [TEA][Tf], we observe that, despite the similarity of the melting temperatures, at T = 320 K the diffusion coefficient of cations is three orders of magnitude slower than that of H2 O in liquid water at the same temperature, while each ion is roughly six times heavier than H2 O. In agreement with NMR measurements,14 we also observe that the cation diffuses somewhat faster than the anion. This might be due to the fact that the anion is slightly heavier than the cation, but the ratio of the square root of the masses accounts for only 20% of the difference, implying that molecular interactions play a relevant role in slowing down [Tf]− . In Figure 4, once again, we find two different linear ranges at low and at high temperature, with a non-strictly-monotonic behaviour in a temperature interval centred on the transition temperature Ta = 310 K. For cations, the high temperature branch corresponds to an activation free energy of 41.5 kJ/mol, that decreases to 16 kJ/mol at low T . For anions, the high temperature branch corresponds to an activation free energy of 45 kJ/mol, that decreases to 25 kJ/mol at low T . It might be useful to remark that the anomaly in D+ , D− on cooling through Ta does not unambiguously correspond to the drop of diffusion rate expected at solidification. Also the decrease of the activation energy on cooling from the liquid to the solid phase are somewhat counter-intuitive, and we think that the precise determination of D+ , D− around and below Ta is affected by transient effects and by the short time scale of simulations. Nevertheless, the low value of the diffusion constants at 12


Page 13 of 34

-6 10 -7 10 -8

log[D- ]

10 -9 10 -10 10 -11 10 -7 10 -8 10

log[D+ ]

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


-8 10 -9 10 -10 10 2.5

3.0

3.5

4.0

4.5

5.0

-1

1000 / T [K ] Figure 4: Arrhenius plot for the diffusion constant. D+ and D− in [cm2 /s]. Tm still identifies the low temperature phase as a solid. A second crucial linear dynamical coefficient is the electrical conductivity, computed from the simulation trajectories according to: 1 1 1 1 lim Π(t) = lim κ= kB T V t→∞ 6t kB T V t→∞ 6t

2 + * X qi [ri (t + t0 ) − ri (t0 )] i

(6)

t0

where V is the system volume and kB is the Boltzmann constant. Because of the partial cancellation of contributions from anions and cations, the time dependent average in Eq. 6 is more noisy than the mean square displacement in Eq. 5 (see examples in Figure S8 in S.I.), and, at equal simulation length, the linear range of Π(t) is restricted to a narrower time 13



Page 14 of 34

interval than in the M SD(t) case. Nevertheless, to avoid biasing the result, we fit the time dependent average in Eq. 6 on the same 10 ns ≤ t ≤ 12 ns already used to estimate diffusion constants. Once again the result is reported in the form of an Arrhenius plot, see Figure 5, displaying the same two linear ranges of the diffusion constants. The activation free energy for the electrical conductivity turns out to be 18 kJ/mol at low T , and 41 kJ/mol at high T . An important aspect of linear transport is represented by the well known relation of diffusion constants and electrical conductivity. The so-called Nernst-Einstein (NE) relation, in particular, provides a simple estimate of the electrical conductivity based on the diffusion coefficients: κN E =

F2 ρ 2 2 ν+ z+ D+ + ν− z− D− RT Mw

(7)

where F is the Faraday constant, R is the gas constant, ρ is the mass density and Mw the molecular weight, ν+ , ν− are the number of cations and anions in a formula unit, z+ and z− are the corresponding ions’ charge, Such a relation assumes that ions diffuse independently, performing uncorrelated random walks (see Sect.10.5 in Ref. 45). Of course this is a very idealised approximation, and the deviation of the actual conductivity with the NE estimate is a measure of correlation. Following the time honored practice, this correlation, or ”pairing” is quantified by the parameter ∆, defined as: κ = κN E (1 − ∆)

(8)

The results are shown in Figure 6. From the picture it is apparent that cation-anion pairing is always a relevant aspect in the dynamics of [TEA][Tf] over a broad range of temperatures, at least up to T = 400. AT T = 310 K the pairing ∆ has a sudden change, paralleling the behaviour observed in most other observables. Below T ∼ 310 K pairing reaches up to 95 % of the ions. At these low temperatures, however, both diffusion and conductivity are already reduced to very low rates by the sluggish dynamics of [TEA][Tf]. A preliminary and partial summary of the results presented up to this point could be stated as follows. The analysis of thermodynamic functions, of diffusion and electrical conductivity point to Ta = 310 K as the temperature at which our simulated systems turn solid on cooling from the liquid state. Above this temperature ions diffuse relatively freely, and carry electric current, despite nonnegligible ions’ pairing. Below this temperature diffusion and electrical conductivity are both very low, and the sudden change of activation free energy tells that the molecular mechanism underlying ions’ translations has changed across the transition. Besides translations, the rotational relaxation of ions represents another basic aspect of the system dynamics. Because of the ions flexibility, rotations are not uniquely defined. Our choice is to 14


Page 15 of 34

0 10

Nernst-Einstein

-1 10

Simulation -2 -2 10

κ

[S/m]

10

-3 10 -4 10 2.5

3.0

3.5

4.0

4.5

5.0

1000 / T [1/K]

Figure 5: Arrhenius plot for the electrical conductivity at pressure P = 1 atm.

1.0

0.8

∆(T)

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


0.6

0.4

κ = κ NE (1 - ∆) 0.2 200

250

300

T

350

400

[K]

Figure 6: Parameter ∆ measuring the pairing of [TEA]+ and [Tf]− ions. Estimated error bar of the order of 5% at T ≤ 310K, slightly lower above T = 310K. characterise this property by analysing the orientation of the N-H termination of [TEA]+ and of the (i)

C-S covalent bond of [Tf]− . To this aim, we associate a unit vector (o1 ) to the N-H covalent bond

15



Page 16 of 34

(j)

in [TEA]+ and another one (o2 ) to the S-C bond at the centre of [Tf]− . In the definition above, (i) and (j) refer to individual ions in the sample. The rotational relaxation time is estimated by first defining the time correlation functions: Θα (t) =

1 X hoi (t + t0 ) · oi (t0 )it0 Nα

(9)

i∈α

where Nα is the (equal) number of cations or anions in the system. Typical results are displayed in Figure 7 A better insight into these properties might be obtained by looking at the logarithm of the autocorrelation function, that becomes a straight line in the case of a simple exponential relaxation. An example of such plot is given in Figure S9 of S. I. Beyond a relatively short time of the order of 0.2 ns, the logarithm of Θα (t) is indeed a straight line, pointing to a relatively long-time (up to ∼ 10 ns at T = 310 K) exponential relaxation process. The 10−1 ns width of the short-time feature is covered by vibrational modes affecting dihedral angles in the covalent backbone of the ions (see Sect. IV). For this reason, we think that this fast relaxation channel corresponds to rotations of the single N-H or C-S bond within each ion, although a completely unambiguous identification is not available, and perhaps cannot be found. Then, it is tempting to identify the longer time relaxation with the rotation of the defining N-H and C-S bonds as part of the whole ions. The analysis of this feature is completed by performing a fit of the linear portion of log Θα (t). The inverse of the linear coefficient represents the opposite of a relaxation time τ that we report again on an Arrhenius-type plot in Figure 8. The result is similar to all the other plots obtained in our study up to now, with distinct low-T and high-T ranges. The crossing between the two regimes, however, takes place at a temperature ∼ 50 K lower than the melting of the simulated samples. In other words, both ions are still able to rotate well below the temperature at which the centre of mass of cations and anions is no longer able to diffuse significantly. The picture provided by the simulation, therefore, is that of a plastic phase arising from the solid some 50 ◦ C below the melting transition, being replaced at Tm = 310K by the homogeneous liquid phase. A similar sequence of transformations is fairly common in molecular systems, and it has been studied by computational approaches since the early days of computer simulation.48 Plastic crystal (or, equivalently, rotor phases) are well known in ionic liquids,49 and in protic ionic systems50 in particular, and [TEA][Tf] thus joins the list of such systems. Ion association and rotational relaxation are likely to be greatly affected by the distribution and properties of hydrogen bonds in the [TEA][Tf] system. In our analysis of simulation trajectories, hydrogen bonding (HB) is defined purely in terms of geometry. We aim at characterising the N-H—O 16


Page 17 of 34

1.0 1.0

Θ + (t)

Θ - (t)

T=280 K

Θ + (t);

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


T=320 K

0.5

T=360 K

0.5 T=400 K 0.0 0

2

4

6

t [ns]

T=360 0.0 0

2

4

6

t [ns] Figure 7: Time auto-correlation function of the N-H direction in [TEA]+ (red lines), and of the C-S bond in [Tf]− (blue line).

hydrogen bonds joining cations and anions, and the gN O radial distribution function already provides an estimate (hrN O i ∼ 2.75 ˚ A at T = 300 K) of the N-O separation in a typical H-bond. In our definition, a triplet of atoms consisting of N-H on the cation side, and an oxygen on the anion side, is a HB if the distance N-O is less than 3.2 ˚ A , and the angle between the separation vectors N − > H and H − −− > O is less than 30◦ . Analysis of trajectories shows that at very low T (T = 200 K) nearly all N-H terminations are engaged in an HB and the number of hydrogen bonds decreases only slowly with increasing temperature, consistent with an estimated strength of about 40 kJ/mol (see Sec. IV), significantly higher than kB T = 2.5 kJ/mol at T = 300 K. Based on the complementary structure of cations and anions, and on intuition, it is tempting to think of the system as made primarily of cation-anion unique pairs. The analysis of simulation snapshots, however, identifies anions with more than one hydrogen bond stretching to different cation partners.17–19 To quantify this qualitative observation, we compute the fraction of anions having 0, 1, 2 and 3 HBs. The cumulative number of HB per 100 [TEA][Tf] ion pairs and its partition in multiplebond contributions is displayed in Figure 9. The figures shows that with increasing temperature the

17



3 10

τ+

2

[ns]

10

τ−

2 10

τ

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 18 of 34

1 10

0 10 2.5

3.0

3.5

4.0

4.5

5.0

-1

1000 / T [K ] Figure 8: Time correlation function of cation and anion orientations. Dots and squares represent simulation results. Lines are a guide to the eye.

largest fractional change is due to fraction of anions forming three simultaneous HBs. An example of a triple-bonded cluster is represented in Figure 10 In our force field model, N in [TEA]+ , with its low positive charge qN = 0.025 e, is unlikely to act as the acceptor in an additional H-bond besides the one donated to oxygen, as discussed in the previous paragraphs. In fact, no close association of N and H is found in the simulation trajectories besides the pair connected by a primary covalent bond. Because of this observation, it is difficult to imagine proton transport paths going beyond vehicular transport, and approaching the Grotthus mechanism familiar from proton mobility in water. This aspect will be discussed in some more detail in the density functional section. Besides the static properties discussed up to now, molecular dynamics gives us access to the kinetics of H-bond formation and dissociation. The probability distribution of the breaking time for individual H-bonds has been computed by identifying all N-H—O hydrogen bonds at a given reference time t0 , and monitoring the persistence of the bond at time t + t0 , until it breaks. Similarly to what has been done with the the other time-dependent properties, the result is averaged over the initial time t0 . Also in this case, it is possible to distinguish a fast relaxation channel in the 0.2 ns time range, and a slower relaxation channel stretching over a few ns. More in detail, the analysis of snapshots shows that up to T = 400 K ions connected by H-bonds stay close to each other for a fairly long 18


Page 19 of 34

100

3-HB 80

2-HB (100)

60

nHB

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


40

1-HB

20

0 200

240

280

320

360

400

T [K]

Figure 9: Cumulative number of HB per 100 [TEA][Tf] ions pairs as a function of T , accounting for the contribution from anions accepting one HB, corresponding to the area below the green line; two HB’s, represented by the area between the green and blue line; or three HBs, corresponding to the strip between the blue and the red lines. Thus, the red line represents the total number of HB per 100 ion pairs.

Figure 10: Cluster consisting of an anion accepting three HBs from neighbouring cations. A thin dash line identifies hydrogen bonds.

19



Page 20 of 34

time. The H-bond itself, however, is broken and reformed at a much faster rate as a result of intraion vibrations.51,52 To decrease the effect of this flickering noise on our analysis, we disregarded all breaking event of H-bonds that reform within 0.4 ps. The result of this mild filtering is shown in Figure S10 of S.I. Only the slowest relaxation process is retained, whose time evolution appears to be a simple exponential. The time constant of the exponential is identified with the genuine lifetime of hydrogen bonds τHB , whose temperature dependence is displayed in Figure 11. Remarkably, the crossover between the low-T and the high-T regime in τHB takes place at T = 260K, consistently with the fact that the freezing of rotations, and the stabilisation of HBs are two aspects of the same dynamical transition. Our protocol to compute τHB is similar to what is currently used in the analysis of H-bonding by simulation, extensively discussed especially in the case of water. Also our results are qualitatively similar to those obtained and discussed in classical papers on this subject53,54 , with only a somewhat longer time scale due to the larger size and mass of [TEA][Tf] ions with respect to water. To monitor possible gradual changes, or even sudden transformations of the geometry of ions as a function of T we computed the principal momenta of inertia of cations and anions at all temperatures. The results confirm basic information already apparent from visual inspection of snapshots. At all temperatures [TEA]+ is globular, and [Tf]− is only somewhat more elongated, justifying the classification of the intermediate phase of [TEA][Tf] as a plastic solid instead of a liquid crystal. More in detail, [Tf]− is a prolate ellipsoid, with I1 = 200 ± 15, I2 ∼ I3 = 360 ± 15 where masses have been expressed in atomic mass units and lengths in ˚ A . [TEA]+ is somewhat less symmetric than [Tf]− , and strictly speaking it should be classified as a spheroid. However, there is a clear separation between the two lowest momenta of inertia and the highest one, and for this reason we we still consider it to be an oblate ellipsoid, with I1 ∼ I2 = 230 ± 3, and I3 = 380 ± 8. The probability distributions for the three momenta of inertia of [TEA]+ and [Tf]− are shown in Figure S11 of S.I. This information could provide a basis for a coarse grained modelling of the system.

B.

The effect of water at low concentration on the structure and dynamics of [TEA][Tf ]

Most if not all RTILs are hygroscopic, including those that are listed as hydrophobic when dissolved in water.55 It is crucial,therefore, to investigate the effect of water at relatively low concentration on a wide range of [TEA][Tf] properties. To this aim, we added 17 water molecules to our samples, corresponding to 1 % weight concentration, placing them at random positions within the simulation box, only taking care of avoiding close atomic contacts.

20


Page 21 of 34

2 10

τHB [ps]

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 10

0 10 2.5

3.0

3.5

1000/T

4.0

4.5

5.0

[K-1]

Figure 11: Lifetime of N-H–O hydrogen bonds as a function of inverse temperature. The straight lines are a guide to the eye. τHB is measured in ps. These new samples have been equilibrated and simulated following precisely the same protocol used for the dry samples, practically doubling the number of simulations and the CPU time requested by our computations. The results largely conform to our expectations, but they are nevertheless useful to fully characterise the properties of [TEA][Tf]. First of all, the average potential energy and volume per ion pair follow closely the behavior seen for the dry samples, apart from a nearly parallel shift due to the extra cohesion and volume due to the added water molecules. Because of this similarity, the results are shown in Figure S13 and Figure S14 of S.I. Then, the OW-OW radial distribution function (OW being the water oxygen) show that at 340 ≤ T ≤ 400K water molecules are distributed homogeneously throughout the system, without a clear tendency to segregate. To be precise, up to T = 400 K gOW OW (r) shows a sharp peak at r ∼ 2.72 ˚ A , due to water-water hydrogen bonding (see Figure S15 in S.I.). The computation of the average coordination number, however, shows that H-bonding among water molecules is quantitatively modest, and only one or two H-Bonded water dimers are typically found at any given time in our samples of 125 [TEA][Tf] and 17 water molecules. Below T ∼ 310 K a liquid-like modulation appears at intermediate range (6 ≤ r ≤ 16 ˚ A ) in the gOW OW (r) radial distribution function, while the sharp peak due to HBs is somewhat reduced. Below 21



Page 22 of 34

T = 300K, however, a sudden transition takes place, changing the shape of the radial distribution function and completely suppressing the water-water HB peak. More detailed analysis of simulation trajectories shows that water molecules form stable associations with pairs of anions, in which water is the HB donor (see Figure S16 in S.I.). A minority of water molecules forms even more complex structures made of two [Tf]− anions and one [TEA]+ cation acting as proton donor, as shown in Figure 12. The structural features described in the previous paragraphs are reflected in (or, equivalently, reflect) the dynamical properties of the simulated samples. First of all, at T ≥ 340 K,water diffuses significantly faster than [TEA]+ or [Tf]− , as can be seen in Table 1. At T ≤ 300 K, however, the tight association of water and [Tf]− anions implies that the two species diffuse at the same low rate. At all temperatures, however, cations and anions diffuse faster in water contaminated samples than in dry samples (see again Table 1), and the enhancement of diffusion increases with increasing temperature. As a result of the increased diffusion rate, also electrical conductivity is enhanced by the addition of water. Remarkably, the increase of conductivity is even larger than the increase of diffusion, as can be seen in Table 2. Moreover, we find that in wet samples above T = 310K, the conductivity computed by MD exceeds the Nernst-Einstein estimate. This result, that does not violate any strict condition, is increasingly found by simulations.56 It could be rationalised by thinking that at sufficiently hight density and fluidity (low viscosity), the motion of one ion along a given direction forces the motion of neighbouring counter-ions in the opposite direction, enhancing electric current beyond diffusion. A similar effect in quantum systems is known as backflow.57 We point out, however, that this discussion is elementary and based on intuition. Rigorous statistical mechanical definitions and concepts might affect both the result and its interpretation. All together, the results on diffusion and conductivity show that water addition changes the properties of the hydrogen-bonding pattern, affecting the degree of ion pairing in their diffusion at medium and high temperature.

IV.

DENSITY FUNCTIONAL COMPUTATIONS AND AB-INITIO SIMULATION

The analysis of the simulation trajectories raises many questions about details of the structure, relative strength of the hydrogen bonding and of Coulomb forces, stability of multiply charged clusters, etc., that we address here through the usage of density functional computations and a limited set of ab-initio (Car-Parrinello) simulations.

22


Page 23 of 34 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


TABLE 1: Diffusion constant of [TEA], [Tf]− and water molecules as a function of temperature in a system made of 125 [TEA][Tf] ion pairs and 17 water molecules (wet samples). Values for dry samples made of 125 ion pairs are reported for a comparison. All values in [cm2 /s] Wet samples Type

T = 240 K

T = 280 K

T = 300 K

T = 320 K

T = 360 K

[TEA]+ 2.1 ± 0.3 × 10−9 5.1 ± 0.2 × 10−9 7.1 ± 0.2 × 10−9 2.2 ± 0.1 × 10−8 1.2 ± 0.1 × 10−7 [Tf]−

5.1 ± 0.5 × 10−10 2.3 ± 0.2 × 10−9 3.8 ± 0.2 × 10−9 1.0 ± 0.1 × 10−8 6.3 ± 0.1 × 10−8

H2 O

5.2 ± 0.3 × 10−10 2.5 ± 0.2 × 10−8 4.3 ± 0.1 × 10−8 2.1 ± 0.1 × 10−7 8.7 ± 0.1 × 10−7 Dry samples

Type +

[TEA]

[Tf]−

T = 240 K

T = 280 K −9

2.5 ± 0.4 × 10

4.4 ± 0.2 × 10

T = 300 K −9

4.9 ± 0.2 × 10

T = 320 K −9

T = 360 K −8

1.9 ± 0.1 × 10

9.7 ± 0.1 × 10−8

5.1 ± 0.5 × 10−10 2.3 ± 0.2 × 10−9 3.8 ± 0.2 × 10−9 1.0 ± 0.1 × 10−8 6.3 ± 0.1 × 10−8

TABLE 2: Electrical conductivity as a function of temperature computed by simulation (Sim., Eq. 6) and estimated by the Nernst-Einstein approximation (NE, Eq. 7) in a system made of 125 [TEA][Tf] ion pairs and 17 water molecules (wet samples). Values for dry samples made of 125 ion pairs are reported for a comparison. All values in Siemens/m. Wet samples Method Sim. NE

T = 240 K

T = 280 K

T = 300 K

T = 320 K

1.1 ± 0.4 × 10−3 4.5 ± 0.4 × 10−3 5.3 ± 0.2 × 10−3 −3

6.3 ± 0.5 × 10

1.5 ± 0.1 × 10

−2

2.05 ± 0.1 × 10

−2

T = 360 K

0.155 ± 0.01 5.5 ± 0.1 × 10

−2

0.35 ± 0.01 0.271 ± 0.05

Dry samples Method

A.

T = 240 K

T = 280 K

T = 300 K

T = 320 K

T = 360 K

Sim.

7.3 ± 0.5 × 10−4 1.6 ± 0.3 × 10−3 1.75 ± 0.2 × 10−3 9.7 ± 0.4 × 10−3 0.103 ± 0.05

NE

1.1 ± 0.4 × 10−2 2.3 ± 0.2 × 10−2 2.6 ± 0.2 × 10−2 6.6 ± 0.3 × 10−2 0.198 ± 0.01

The molecule, the ions, and the energetics of ionisation and ion dissociation

Computations have been carried out to determine, first of all, the ground state geometry of the different species; the energetics of the neutral versus ionic state of TEA and TF and their combination in the [TEA][Tf] molecule; the geometry and strength of the hydrogen bond contributing to the cation-anion binding energy (together with Coulomb and dispersion forces). A picture of [TEA][Tf] in its gas-phase equilibrium configuration as determined by Gaussian 09 is shown in Figure 13. The closest contact between the two ions corresponds to a short and presumably 23



Page 24 of 34

Figure 12: Typical configuration of a water molecule donating H-bonds to two [Tf]− anions, and accepting a third one from a [TEA]+ cation. A thin dash line identifies hydrogen bonds. strong hydrogen bond connecting the nitrogen atom on [TEA]+ to one of the oxygen atoms on [Tf]− . The corresponding N-O distance is 2.63 ˚ A according to Gaussian 09, and 2.59 ˚ A according to CPMD. The N-H–O angle is 172◦ or 168◦ according to Gaussian 09 and CPMD, respectively. Such a short bond length suggests that the potential energy surface for the proton in the HB has a single minimum, in contrast to the double well potential of longer and looser hydrogen bonds.58 In experiments, such a feature could be probed by vibrational spectroscopy, and it could also affect the proton conductivity at low and intermediate temperature. The computational verification of the single-minimum property of the proton energy is illustrated below. The analysis of vibrational frequencies carried out by CPMD shows that the frequency of the N-H stretching goes from 3357 cm−1 in the [TEA]+ ion to 2206 cm−1 in the [TEA][Tf] molecule. Together with the short N-O distance in [TEA][Tf], this very large downwards shift emphasises the strength of the hydrogen bond. This is further confirmed by the fact that at 1550 cm−1 the N-H–O mode is the stiffest of all bendings. For the sake of completeness, we mention that 15 C-H stretching modes give origin to a narrow band from 2950 cm−1 and 3100 cm−1 , bending modes involving C-H pairs are found between 1150 and 1450 cm−1 . Modes localised on [Tf]− tend to have lower frequency than those on [TEA]+ . C-F stretching modes, for instance occur at 1050 - 1140 cm−1 , while S-N stretching is found at 820 - 900 cm−1 . Below 550 cm−1 the spectrum is made of a variety of modes affecting primarily dihedral angles. While the hydrogen bond is the most apparent structural feature, the cohesive energy of the

24


Page 25 of 34 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 13: Ground state geometry of [TEA][Tf] from density functional computations. PBE exchange correlation approximation, and Gaussian 09 computation (see Sec. II). The dash line represents a hydrogen bond connecting the cation to the anion.

ionic [TEA][Tf] molecule seems to be dominated, as expected, by Coulomb forces. This statement is supported by the following analysis, representing the formation of the ionic [TEA][Tf] pair out of the neutral TEA and Tf molecules in terms of a few more elementary steps. According to G09, the energy of the reaction: TEA + Tf

−→

[TEA][Tf]

(gas phase)

is ∆E = 1.19 eV, or 27.4 kcal/mol. In our convention, positive ∆E means that the reaction runs spontaneously from left to right. This relatively modest binding energy results from balancing the energy required to move one proton from the gas phase neutral Tf molecule to the gas phase neutral TEA molecule, and the larger gain of Coulomb energy in bringing the two gas phase ions together, making one globally neutral but ionic [TEA][Tf] molecule. The proton transfer reaction, in particular, could be seen as the combination of two simpler steps: Tf

−→

[Tf]− +p

(gas phase)

and: TEA + p

[TEA]+

−→

(gas phase)

where p is a proton, and the combination of these two steps requires ∆E = −2.88 eV, or −66.4 kcal/mol.

25



Page 26 of 34

Then, the energy gain in joining [TEA]+ to [TF]− according to the reaction: [TEA]+ + [TF]−

−→

[TEA][TF]

(gas phase)

is ∆E = 4.07 eV, or 93.9 kcal/mol. This last and large ∆E, in particular, can be seen as the Coulomb interaction of the two ions, and apparently determines the stability of the ionic [TEA][Tf] gas-phase molecule. The computed stability of the two neutral and isolated TEA and Tf molecules with respect to the gas phase ions [TEA]+ and [Tf]− is expected, since the ionised state of even the most ionic compounds (NaCl, for instance) is always due to the mutual interaction of the two ions. Additional details on the stability of the ionic versus the neutral state of TEA and Tf units are given below. All the energies reported in this sub-section have been computed by G09. Computations carried out by CPMD give results fully compatible with those of G09, but somewhat more uncertain for some residual difficulty in treating globally charged systems by plane waves. Moreover, all cohesive energies do not contain the zero-point contribution from vibrations. However, all vibrational frequencies have been computed and we verified that addition of the zero-point energy does not change significantly the results for the cohesive energy. Consistently with the room temperature ionic liquid character of [TEA][Tf], the molecular cohesive energy ∆E = 1.19 eV is modest on the scale of covalent or ionic systems, and larger than typical cohesive energies of H-bonds, however strong. On the other hand, this energy is still large on the scale of thermal energies at about room temperature. The analysis of the charge distribution provided by CPMD for all neutral and anionic species confirms that N tends to have a slightly positive charge, and supports our previous statement saying that N is unlikely to act as a proton acceptor in a H-bond. This, in turn, makes it difficult to imagine proton transfer paths in [TEA][Tf] reminiscent of the Grotthus mechanism in water. The picture is different in closely related compounds based on di-amines, like those discussed in Ref. 14,15 An alternative way for the proton to escape strict vehicular motion on board of [TEA]+ would require its transfer back to a [Tf]− ions, producing a neutral TEA and Tf pair, with the proton moving together with Tf before being transferred to a different neutral TEA unit. To determine whether the proton transfer is energetically possible and the associated barrier not overwhelming in the temperature range of interest (200 ≤ T ≤ 400 K), the molecular energy of [TEA][Tf] has been computed upon constraining the distance of the hydrogen bound to N in [TEA]+ , progressively moving it towards the closest oxygen on [Tf]− . All other degrees of freedom are relaxed to their minimum energy. Whenever the two ions [TEA]+ and [Tf]− are in close contact, E(H − −O) has a unique minimum, corresponding to the equilibrium distance in the [TEA][Tf] ground state (see Figure 14). 26


Page 27 of 34

[kcal/mol]

15

10

5

∆E

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


0 0.8

1.2

1.6

dO--H

2.0

[A]

Figure 14: Energy variation along the constrained coordinate measuring the dH−−O distance. Decreasing dH−−O means moving the proton towards the nearest [Tf]− oxygen, giving origin to a neutral TEA and Tf pair. Blue dots: computations for the [TEA][Tf] molecule. Red filled squares: computation for the periodic structure shown in Figure S18 of S.I. The full and dash black curves are a guide to the eye.

In other terms, the presence of TEA at short distance from Tf, destabilises the neutral pair, which is instead metastable when the two ions are well separated (not shown in the figures). The single minimum in the E versus dH−−O distance (see Figure 14) is one further point supporting the strict vehicular mechanism as the only relevant proton conductivity mechanism in [TEA][Tf].

V.

SUMMARY AND CONCLUSIONS

Molecular dynamics simulations based on a simple, unpolarisable force field model, covering the 200 K ≤ T ≤ 400K temperature range, and lasting ∼ 100 ns provide crucial insight to identify the transitions in [TEA][Tf] seen by calorimetry and in neutron scattering measurements.13 The experimental results, in particular, reveal two transitions on cooling, taking place on heating at T = 230K and T = 310K on heating, and at T = 260 K and T = 210K on cooling. Calorimetry shows that both transitions are first order, and connect the liquid phase at T ≥ 310 K to the crystal phase at T ≤ 210 K, with an intermediate solid phase of uncertain identification in between. The 27



Page 28 of 34

macroscopic time scale of experiments makes the two transitions broadly reversible, despite sizeable hysteresis. Simulations of [TEA][Tf] on cooling from T = 400 down to T = 200K show that the homogeneous liquid undergoes a weakly first order transition at T = 310 K. At the transition temperature the diffusion coefficient of both ions and the electric conductivity are simultaneously reduced to very low rates, while the rotation of the ions still represents an effective relaxation channel. These properties identify T = 310 K as the melting temperature (Tm = 310 K) of the simulated system, and the phase at 260 K ≤ T ≤ 310K as a plastic solid. The freezing of the system into a solid configuration is confirmed by quantitatively modest but qualitatively revealing changes in the structure as described by radial distribution functions. The observation of solidification by simulation is remarkable, especially for a molecular system of the [TEA][Tf] complexity, that usually become glassy on cooling at the high rates of molecular dynamics. No other discontinuous transition is apparent in the thermodynamic functions of the simulated sample down to T = 200K. Simulation of the heating cycle is prevented by the lack of an experimentally determined crystal structure. In the simulated samples, however, T = 260 K marks the change in the rotational diffusion of ions, that virtually stops at this temperature. The significance of this observation is emphasised by the change of the activation energy for rotational diffusion, pointing to a qualitative change in the rotational dynamics of the ions. We conclude that these effects are the simulation counterpart of the first order transformation seen in the experiment ∼ 50 K below melting, since the change from discontinuous to continuous transitions is a rather common consequence of the high cooling rates of molecular dynamics. The atomistic resolution of the force field model then provides a wealth of additional microscopic information on the thermodynamic functions, the structure, and the linear dynamical coefficients of [TEA][Tf] samples. For instance, comparison of electrical conductivity computed during the simulations with the estimate based on the Nernst-Einstein model allows us to quantify the effect of correlation in the diffusion of ions, that we express in terms of ion pairing. An important aspect in the system properties is represented by the structure and kinetics of the population of hydrogen bonds in [TEA][Tf] that greatly affect measurable properties of the system. Also in this case, molecular dynamics trajectories allow to quantify a variety of parameters, that could be compared with experiments, and also with previous computations on the same system.20 Most of the quantities that we analyse depend on activated processes, that can be easily identified by their characteristic dependence on temperature. Arrhenius-type plots on the log-(1/T) plane, in particular, show that these properties 28


Page 29 of 34 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


display piecewise-linear ranges, whose linear coefficient allows us to estimate activation free energies for a variety of processes in the three phases displayed by the simulated samples of [TEA][Tf]. Detailed analysis of simulation snapshots reveals the formation of unexpected complexes made by up to three cations connected by H-bonding to a single anion. In the condensed phase, overall neutrality is enforced by screening from the neighboring ions, and locally charged clusters are thermodynamically stable. Density functional computations for gas-phase species show that only the smallest multiply charged clusters are stable or at least metastable, but also highlight sizeable stabilising effect arising from ions’ correlation, and from the formation of hydrogen bonds. In this respect, and despite quantitative differences, our results are reminiscent of those from previous studies, devoted to similar charged complexes in other protic room-temperature ionic liquids.17,18 Vibrational properties of protic ionic liquids closely related to [TEA][Tf] have been measured experimentally.59–61 The quantitative comparison of computational and experimental data on vibrational properties would be highly valuable, but would also require a focused effort to quantitatively analyse the relative infrared or Raman weight of each mode, and to include anharmonic effects, that are all important for the low-frequency modes that carry information on inter-ion interactions. Such a detailed analysis has not been carried out yet, but remains as an appealing direction of quantitative research for the near future. The role of water contamination on the properties of [TEA][Tf] sample has been investigated by simulations, revealing no major effect on the thermodynamics of the system, but providing an intriguing view of characteristic changes in transport coefficients and in the structure due to the water ability to form strong H-bonds donated to the [Tf]− anion, and, to a lesser extent, accepted from the [TEA]+ cations. It is found, for instance, that water enhances the ion mobility down to the solidification point, increasing the diffusion coefficient of ions, and, even more, the electrical conductivity. Analysis of simulation snapshots reveals intriguing local clusters made of [TEA]+ and, again, multiple [Tf]− anions, joined by a bridging water molecule. Besides these positive aspects of our study, it is important also to mention the limitations. It is well known that unpolarisable force field models are certainly only a first approximation to the potential energy surface of the system. Quantitative inaccuracies and even qualitative failures due, for instance, to the lack of polarisability in the model are apparent, and are not easy to overcome. The fixed bonding topology adopted by the force field limits the simulation ability to analyse proton transport mechanisms besides vehicular. This limitation might not be crucial for [TEA][Tf], since ab-initio analysis shows that a change of bond topology from [TEA][Tf] is indeed very unlikely even for the acidic N-H termination, but, in general it is a strong drawback of the rigid ion model. Approaches to 29



Page 30 of 34

overcome the problem range from reactive force fields,29 to ab-initio simulations. These last, however, are still expensive and limited to time scales orders of magnitude shorter that needed for a system such as [TEA][Tf] at its triple point. Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: xxxxxxx. Definition of the force field; miscellaneous figures; search for a low temperature periodic structure of TEA][Tf]. (PDF)

REFERENCES 1

Hallett, J. P.; Welton, T. Room-temperature ionic liquids: Solvents for synthesis and catalysis. 2. Chem. Rev. 2010, 111, 3508-3576.

2

Greaves, T. L.; Drummond, C. J. Protic ionic liquids: properties and applications Chem. Rev. 2008, 108, 206-237.

3

Greaves, T. L.; Drummond, C. J. Protic ionic liquids: Evolving structure-property relationship and expanding applications. Chem. Rev. 2008, 108, 206-237.

4

Ghandi, K. A Review of ionic liquids, their limits and applications. Green and Sustainable Chemistry 2014, 4, 44-53.

5

Matsuoka, H.; Nakamoto, H.; Susan, M. A. B. H.; Watanabe, M. Brønsted acid-base and -polybase complexes as electrolytes for fuel cells under non-humidifying conditions. Electrochimica Acta 2005, 50, 4015-4021.

6

Nakamoto, H.; Watanabe, M. Brønsted acid-base ionic liquids for fuel cell electrolytes Chem. Commun. 2007, 24, 2539-2541.

7

MacFarlane, D. R.; Tachikawa, N.; Forsyth,M.; Pringle J. M.; Howlett, P. C.; Elliott, G. D.; Davis, J. H.; Watanabe, M.; Simon, P. Angell, C. A. Energy applications of ionic liquids. Energy Environ. Sci. 2014, 7, 232-250.

8

Mayrand-Provencher, L.; Lin, S.; Lazzerini, D.; Rochefort, D. Pyridinium-based protic ionic liquids as electrolytes for RuO2 electrochemical capacitors. Journal of Power Sources 2010, 195, 5114-5121.

9

Timperman, L.; Skowron, P.; Boisset, A.; Galiano, H.; Lemordant, D.; Frackowiak, E.; Beguin, F.; Anouti, M. Triethylammonium bis(tetrafluoromethylsulfonyl)amide protic ionic liquid as an electrolyte for electrical double-layer capacitors Phys. Chem. Chem. Phys. 2012, 14, 8199-8207.

10

Timperman, L.; Galiano, H.; Lemordant, D.; Anouti M. Phosphonium-based protic ionic liquid as electrolyte for carbon-based supercapacitors. Electrochem. Commun. 2011, 13, 1112-1115.

11

Menne, S.; Pires, J.; Anouti, M.; Balducci, A. Protic ionic liquids as electrolytes for lithium-ion batteries. Electrochem. Commun. 2013, 31, 39-41.

30


Page 31 of 34 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


12

Abate, A.; Hollman, D. J.; Teuscher, J.; Pathak, S.; Avolio, R.; D’Errico, G.; Vitiello, G.; Fantacci, S. Snaith, H. J. Protic ionic liquids as pdopant for organic hole transporting materials and their application in high efficiency hybrid solar cells. J. Am. Chem. Soc. 2013, 135, 13538-13548.

13

Burankova, T.; Hempelmann, R.; Fossog, V.; Ollivier, J.; Seydel, T.; Embs, J. Proton diffusivity in the protic ionic liquid triethylammonium triflate probed by quasi-elastic neutron scattering J. Phys. Chem. B 2015, 119, 10643-10651.

14

Iojoiu, C.; Hana, M.; Molmeret, Y.; Martinez M.; Cointeaux, L.; El Kissi, N.; Teles, J.; Leprˆetre, J.-C.; Judeinstrein, P.; Sanchez, J.-Y. Ionic liquids and their hosting by polymers for HT-PEMFC membranes. Fuel Cells 2010, 10, 778-789.

15

Iojoiu, C.; Martinez, M.; Hanna, M.; Molmeret, Y.; Cointeaux, L.; Lepretre, J.-C.; El Kissi, N.; Guindet, J.; Judeinstein, P.; Sanchez.J.-Y. PILs-based Nafion membranes: a route to high-temperature PEFMCs dedicated to electric and hybrid vehicles. Polymers for Advanced Technologies 2008, 19 1406-1414.

16

Martinez, M.; Molmeret, Y.; Cointeaux, L.; Iojoiu, C.; Lepretre, J.-C.; Kissi, N.E.; Judeinstein, P.; Sanchez, J.-Y. Proton conducting ionic liquid-based proton exchange membrane fuel cell membranes: the key role of ionomer-ionic liquid interaction. J. Power Sources 2010, 195, 5829-5839.

17

Knorr, A.; Ludwig, R. Cation-cation clusters in ionic liquids: cooperative hydrogen bonding overcomes like-charge repulsion. 2015, Sci. Rep., 5, 17505.

18

Strate, A.; Niemann, T.; Ludwig, R. Controlling the kinetics and thermodynamic stability of cationic clusters by the addition of molecules or counterions. 2017 Phys. Chem. Chem. Phys., 19, 18854-18862.

19

Knorr, A.; Fumino, K.; Bonsa, A.-M.; Ludwig, R. Spectroscopic evidence of jumping and pecking of cholinium and H-bond enhances cation-cation interaction in ionic liquids. Phys. Chem. Chem. Phys. 2015, 17, 3097830982.

20

Sunda, A. P.; Mondal, A.; Balasubramanian, S. Atomistic simulations of ammonium-based protic ionic liquids: steric effects on structure, low frequency vibrational modes and electrical conductivity. Phys. Chem. Chem. Phys. 2015, 17, 4625-4633.

21

Nasrabadi, A. T.; Gelb, L. D. Structural and transport properties of tertiary ammonium triflate ionic liquids: A molecular dynamics study. J. Phys. Chem. B 2017, 121, 1908-1921.

22

Jensen, K. P.; Jorgensen, W. L. Halide, ammonium, and alkali metal ion parameters for modeling aqueous solutions. J. Chem. Theory Comput. 2006, 2, 1499-1509.

23

Jorgensen W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236.

24

Lopes, J. N. C.; Padua, A. A. H. Molecular force field for ionic liquids composed of triflate or bistriflylimide anions. J. Phys. Chem. B, 2004, 108, 16893-16898.

25

Lopes, J. N. C.; Padua, A. A. H. CL&P: Generic and systematic force field model for ionic liquids modeling Theor. Chem. Acc. 2012, 131, 1129.

26

Cornell, W.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A 2nd generation force-field for the simulation of proteins, nucleic acids,

31



Page 32 of 34

and organic-molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. 27

Breneman, C.M.; Wiberg, K.B. Determining atom-centered monopoles from molecular electrostatic potentials. J. Comput. Chem. 1990; 11, 361-373.

28

Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129-145.

29

Hartke, B.; Grimme, S. Reactive force fields made simple. Phys. Chem. Chem. Phys. 2015, 17, 16715-16718.

30

Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel moleculardynamics implementation. Comp. Phys. Comm. 1995, 91, 43-56

31

Eastwood, J. W.; Hockney, R. W.; Lawrence, D. N. P3M3DP - The 3-dimensional periodic particle-particleparticle-mesh program. Computer Physics Communications 1980, 19, 215-261.

32

Mondal, A.; Balasubramanian, S. Quantitative prediction of physical properties of imidazolium based room temperature ionic liquids through determination of condensed phase site charges: A refined force field J. Phys. Chem. B 2014, 118, 3409-3422.

33

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865-3868.

34

Gaussian-09, Revision E.01, Frisch, M. J.; et al., Gaussian, Inc., Wallingford CT, 2009.

35

Becke, A. D. Densityfunctional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98 5648-5652.

36

CPMD, http://www.cpmd.org/ (last accessed on 11/28/2017), Copyright IBM Corp 1990-2015, Copyright MPI f¨ ur Festk¨orperforschung Stuttgart 1997-2001.

37

Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43, 1993-2006.

38

Wales, D. Insight into reaction coordinates and dynamics from the potential energy landscape. J. Chem. Phys. 2015, 142, 130901.

39

Kirkpatrick, S.; Gelatt Jr, C.D.; Vecchi, M. P. Optimization by simulated annealing. Science 1983, 220, 671-680.

40

Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comp. Chem. 2006, 27, 1787-1789.

41

Forero-Martinez N. C.; Cortes-Huerto, R.; Ballone, P. The glass transition and the distribution of voids in room-temperature ionic liquids: a molecular dynamics study. J. Chem. Phys. 2012, 136, 204510.

42

Cabriolu, R.; Ballone, P. Thermodynamic properties and atomistic structure of the dry amorphous silica surface from a reactive force field model, Phys. Rev. B, 2010, 81, 155432.

43

Belieres, J.-P.; Angell, C. A. Protic ionic liquids: preparation, characterization, and proton free energy level representation. J. Phys. Chem. B 2007, 111, 4926-4937.

44

Mamontov,

E.;

Luo,

H.;

Dai,

S.

Proton

dynamics

in

N,N,N,N-tetramethylguanidinium

bis(perfluoroethylsulfonyl)imide protic ionic liquid probed by quasielastic neutron scattering. J. Phys. Chem. B 2009, 113, 159-169.

32


Page 33 of 34 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


45

Hansen, J.-P. and McDonald, I. R. Theory of Simple Liquids, Elsevier, Amsterdam (2006).

46

Bathia, A.B.; Thornton, D. E. Structural aspects of the electrical resistivity of binary alloys, Phys. Rev. B 1970, 2, 3004-3012.

47

Stillinger, F. H.; Lovett, R. General restriction on distribution of ions in electrolytes. J. Chem. Phys. 1968, 49, 1991-1994.

48

Ciccotti, G.; Ferrario,M.; Memeo, E.; Meyer, M. Structural transition on cooling of plastic adamantane: A molecular-dynamics study Phys. Rev.Lett. 1987, 59, 2574-2577.

49

Luo,J.-S.; Conrad, O.; Vankelecom, I. F. J. Imidazolium methanesulfonate as a high temperature proton conductor. J. Materials Chem. A 2013 1

50

Luo, J.-S.; et al. 1,2,4-Triazolium perfluorobutanesulfonate as an archetypal pure protic organic ionic plastic crystal electrolyte for all-solid-state fuel cells. Energy & Environ. Sci. 2015 8, 1276-1291.

51

Hammerich,A.D.; Buch, V. An alternative near-neighbor definition of hydrogen bonding in water. J. Chem. Phys. 2008, 128, 111101.

52

Voloshin, V.P.; Naberukhin, Yu. I. Hydrogen bond lifetime distribution in computer simulated water. Journal of Structural Chemistry 2009, 50, 78-89.

53

Luzar, A.; Chandler, D. Effect of environment on hydrogen bond dynamics in liquid water. Phys. Rev. Lett. 1996, 76, 928-931.

54

Luzar A. Resolving the hydrogen bond dynamics conundrum. J. Chem. Phys. 2000, 113, 10663-10675.

55

Cammarata, L.; Kazarian, S. G.; Satter, P. A.; Welton, T. Molecular states of water in room temperature ionic liquids. Phys. Chem. Chem. Phys. 2001, 3, 5192-5200.

56

Armstrong, J. A.; Ballone, P. Computational verification of two universal relations for simple ionic liquids. Kinetic properties of a model 2:1 molten salt. J. Phys. Chem. B 2011 115, 4927-4238.

57

Pines, D.; Nozieres, P. Theory of Quantum Liquids: Normal Fermi liquids, Addison-Wesley, Boston (1989).

58

Kreuer, K. D. Proton conductivity: materials and applications. Chem. Mater. 1996, 8, 610-641.

59

Fumino, K.; Reichert, E.; Wittler, K.; Hempelmann, R.; Ludwig, R. Low-frequency vibrational modes of protic molten salts and ionic liquids: detecting and quantifying hydrogen bonds. Angew. Chem., Int. Ed. 2012, 51, 6236-6240.

60

Fumino, K.; Fossog, V.; Wittler, K.; Hempelmann,R.; Ludwig, R. Dissecting anion-cation interaction energies in protic ionic liquids. Angew. Chem., Int. Ed. 2013, 52, 2368-2372.

61

Fumino, K.; Reimann, S.; Ludwig, R. Probing molecular interaction in ionic liquids by low frequency spectroscopy: Coulomb energy, hydrogen bonding and dispersion forces. Phys. Chem. Chem. Phys. 2014, 16, 21903-21929.

33


E〉(T) 〈

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 34 of 34

[kJ/mol]


T

TOC graphics

34


[K]

Density Functional Computations and Molecular Dynamics

Recommend Documents