Equilibrium Structure, Hydrogen Bonding, and Proton Conductivity in

Subscriber access provided by UNIV OF LOUISIANA

B: Fluid Interfaces, Colloids, Polymers, Soft Matter, Surfactants, and Glassy Materials

Equilibrium Structure, Hydrogen Bonding and Proton Conductivity in Half-Neutralised Diamine Ionic Liquids Juan F. Mora Cardozo, Jan Peter Embs, Antonio Benedetto, and Pietro Ballone J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b00890 • Publication Date (Web): 15 Mar 2019 Downloaded from http://pubs.acs.org on March 17, 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 41 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

Equilibrium Structure, Hydrogen Bonding and Proton Conductivity in Half-neutralised Diamine Ionic Liquids Juan F. Mora Cardoso(1) , J. P. Embs(1) , A. Benedetto(1,2,3,4) , †, and P. Ballone(5) † (1) Laboratory for Neutron Scattering and Imaging, Paul Scherrer Institute, Villigen PSI, Villigen 5232, Switzerland (2) School of Chemistry, UCD, Dublin, Ireland (3) School of Physics, UCD, Dublin, Ireland (4) Department of Sciences, University of Roma Tre, Via della Vasca Navale 84, 00146 Rome, Italy and (5) Italian Institute of Technology, Via Morego 30, 16163 Genova, Italy

Abstract Recent experiments on proton conducting ionic liquids point to half-neutralised diamine-triflate salts as promising candidates for applications in power generation and energy conversion electrochemical devices. Structural and dynamical properties of the simplest among these compounds are investigated by a combination of density-functional theory (DFT) and molecular dynamics (MD) simulation based on an empirical force field. Three different cations have been considered, consisting of a pair of amine-ammonium terminations joined by a short aliphatic segment -(CH2 )n - with n = 2, 3 and 4. First, the ground state structure, vibrational eigenstates and hydrogen-bonding properties of single ions, neutral ion pairs, small neutral aggregates of up to eight ions and molecularly-thin hydrogen bonded wires have been investigated by DFT computations. Second, structural and dynamical properties of homogeneous liquid and amorphous phases are investigated by MD simulations over the temperature range 200 K ≤ T ≤ 440 K. Structure factors, radial distribution functions, diffusion coefficient and electrical conductivity are computed and discussed, highlighting the inherent structural heterogeneity of these compounds. The core investigation, however, is the characterisation of connected paths consisting of cation chains that could support proton transport via a Grotthuss-type mechanism. Since simulations are carried out using a force field of fixed bonding topology, this analysis is based on the equilibrium structure only, using geometrical criteria to identify potential paths for proton conduction. Paths of connected cations can reach a length of 80 cations and 30 ˚ A , provided bridging oxygen atoms from triflate anions are taken into account. The effects of water contamination at the 1 % weight concentration on structure, dynamics and paths for proton transport are discussed.

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 41

INTRODUCTION

The increasing need for high power, high capacity electrochemical energy conversion devices such as batteries, supercapacitors, fuel and photovoltaic cells is driving the quest for innovative electrolytes.1 In this context, room temperature ionic liquids (RTILs) are posed to play an important role because of their fluidity at ambient conditions and at anhydrous composition, and because of general properties of RTILs such as low vapour pressure, high thermal stability, wide electrochemical window. Protic ionic liquids2 appear to be particularly promising, since protons represent a general charge and energy token, their low mass favouring high mobility, while the moderate cohesion, the selectivity and directionality of hydrogen bonding (HB) provide a range of opportunities for tailoring structural and dynamical properties of these materials. One further advantage of protic RTILs is the prospect of proton conduction through the analogue of Grotthuss’ mechanism in water, which could overcome the limitation in the mobility of whole ions, and persist at low temperatures because of quantum mechanical assistance. Protic RTILs, moreover, are easily synthesised by neutralising a Brønsted acid with a Lewis base. High density of protons and hydrogen bonds, however, is a necessary but not a sufficient condition for achieving Grotthuss conductivity, since the topology of the hydrogen bonding network is also relevant. To be able to contribute to the collective proton displacement, in particular, each unit involved in Grotthuss transport needs a double functionality as a donor and as an acceptor. Moreover, for conductivity, i. e., long range transport, the sequence of elementary links has to be connected into an extended chain or network. These requirements limit the number of potential candidates able to exceed the vehicular conductivity, which in RTILs is hampered by their high viscosity. A comprehensive experimental study3 recently identified aliphatic diamine species as promising candidates for proton conductivity at anhydrous condition. Compounds of this kind are obtained by neutralising the trifluoromethanesulphonic (triflic) acid by the Lewis base represented by the aliphatic diamine. The result of protonating only one of the two functions is a compound of low melting temperature, criss-crossed by a multitude of intra- and inter-ion hydrogen bonds. Measurement of diffusion coefficients by NMR revealed an excess mobility of acidic protons with respect to cations, pointing to a proton transport mechanism beyond vehicular. In our study we investigate structural and dynamical properties of the shortest aliphatic diamines, focusing, in particular, on the lightest member of this class, i.e., ethylenediamine, DAEt, in which the carbon chain is represented by -(CH2 )2 −. The corresponding salt is [DAEt][Tf], where [Tf]− is CF3 SO3 − . Less extensive computations for longer-chain diamines [-(CH2 )n -, n = 3, 4] are also

2


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


reported, aiming at a preliminary assessment of the role of the chain length. First, by density functional computation we identify basic structural properties, including the possibility of intra-ion and inter-ion hydrogen bonding. Again by DFT, we compute the potential energy profile upon moving protons along each type of hydrogen bond. Whenever this profile has a double-minimum shape, the barrier in between has been determined. Then, extensive MD simulations have been carried out, using a rigid ion, fixed topology force field model and molecular dynamics in the NPT ensemble to compute structural and dynamical properties on an extended temperature range covering liquid and glassy conditions. Already in the smallest [DAEt][Tf] aggregates, the results of the ab-initio stage highlight the expected predominance of Coulomb forces, but also emphasise the structural and cohesive role of hydrogen bonding. Besides strong intra-cation HBs, several HBs link cations to anions, but also cations to cations. In the majority of cases, the potential energy profile for protons moving along these HBs has a single minimum, but it is also very soft, and protons experience a fair degree of delocalisation, that in reality is certainly greatly enhanced by quantum effects. Both the intra-cation and the cation-cation bond present two energy minima, separated by a moderate barrier, compatible with proton transfer at T ∼ 300 K, that, once again, could be eased by quantum effects. An important issue addressed but not fully solved by our investigation is the precise mechanism that could underly proton transport along cation chains. Comparison of energy barriers for single-proton vs collective proton displacement computed by DFT on model 1D [DAEt][Tf] wires, however, suggests that proton conductivity could occur by propagation of a soliton-like excitation along a chain of hydrogen bonded cations. The temperature dependence of the thermodynamic properties of samples simulated by classical MD shows a broad change centred at T ∼ 310 K, about 40 K below the experimental melting point of these compounds. In our computations this change appears to be predominantly a glass transition with superimposed a weak first order change, whose structural origin is not apparent in our results. The transition is reflected in dynamical properties such as the diffusion coefficient of cations and anions, electrical conductivity, rotational dynamics of both ions. At all temperatures the equilibrium structure, as characterised by the structure factor and by the radial distribution function, displays a clear signature of nanostructuring, already found in a number of RTILs of comparable molecular weight and complexity.4 With the force field model underlying our MD simulations, we cannot follow in real-time the transfer of protons along HB and among ions, hence our explicit determination of proton mobility covers the vehicular mechanism only. However, we search and identify paths in the equilibrium structure suitable 3



Page 4 of 41

for macroscopic proton transport, although we are unable to quantify their impact on the transport rate. The results show that chains, or, better, clusters of up to 80 hydrogen bonded cations are represented in the equilibrium samples. The concentration of chains longer than 40 cations increases monotonically with decreasing temperature. Since the number ions in the sample is fixed, the trend towards large connected clusters is compensated by a slight relative decrease of the (large) population of clusters with less than 40 cations. The size probability distribution of all these aggregates is computed and discussed in comparison with the prediction of equilibrium polymerisation.5 Similarly to what has been done in experiments,3 we also investigate the effect of water at 1 % weight concentration on various properties of these compounds. As widely known from experiments and previous computations, contamination by water enhances the mobility of RTIL ions, and thus increases diffusion, electrical conductivity and also proton conductivity, both vehicular and, whenever active, by Grotthuss mechanism. The topics covered in the present paper are closely related to active research fields widely represented in the literature. Density functional computations, in particular, have been carried out to investigate proton transport in ionic-liquid doped membranes,6 and in an organic compound closely related to RTILs such as neutral imidazole chains.7,8 Molecular dynamics simulations based on an empirical force field have been used to study structural properties and ion transport in protic ionic liquids consisting of triflate anions and tertiary ammonium cations with alkyl chains of the [Nknm]+ type, where knm labels alkyl tails k-, n- and m-carbon long, such as [N113], [N122], [N133], [N223].9 The effect of water contamination on RTILs, and the effect of RTILs on water properties in water/RTIL mixtures has also been investigates many times, so much that we can cite only a few examples.10–13 Despite many general features found in all these systems, the full exploration of this broad subject is far from complete. Besides and beyond RTILs, a broader context for our study is the decades-long interest of the experimental and computational communities for proton conducting materials.14 For condensed systems, computational investigations in this field are the realm of density functional computations, and ab-initio simulation in particular. Again DFT and, less often, wave function-based quantum chemistry methods are used in the case of molecular models. Methods such as the empirical valence bond (EVB) model have also been applied with success to proton transport in extended systems.8 Narrowing the view to the activity of our group, the present study is related to a previous one devoted to triethylammonium triflate ([TEA][Tf]),15 whose results provide a useful comparison. Both studies are part of a broader investigation of proton conductivity in room temperature ionic liquids using both experiments and computational approaches.16 4


Page 5 of 41 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


II.

METHODS

The bulk of our computational study consists of molecular dynamics simulations based on an atomistic force field, whose functional form corresponds to the OPLS/Amber model,17,18 with fixed molecular topology, fixed atomic charges and no explicit atomic polarisability. Most parameters for [DAEt]+ and [DABu]+ were taken from Ref. 19, while for [Tf]− the parametrisation from Ref. 20 has been used. Atomic charges have been computed from our density functional computations, using the electrostatic potential model, ESP, and its refinement, RESP.21,22 For the sake of definiteness we give all parameters as Supporting Information. In the present study the rigid ion approximation is the major source of uncertainty in the computation of general chemical physics properties, while the the ability of the model to cover proton transport beyond the vehicular mechanism is limited by the fixed topology, although we endeavour to provide indirect information on the basis of geometrical criteria. Force-field-based MD simulations have been carried out using the Gromacs package23 with a time step of 1 fs. NPT conditions have been enforced introducing a Nosé-Hoover thermostat and ParrinelloRahman barostat, both with a relaxation time constant of 2 ps. Long-range Coulomb interactions are dealt with by the particle-particle-particle-mesh Ewald approach.24 All degrees of freedom were kept unconstrained, including stretching and bending. To simulate samples contaminated by water, the force field has been supplemented by the TIP4P potential.25 The exploration of structural features and hydrogen bonding in ion pairs and in very small clusters at the density functional level has been carried out by resorting to the pseudopotential-plane wave approach implemented in CPMD.26,27 In our computation we used the generalised gradient approximation of Perdew-Burke-Ernzerhof28 (PBE), norm-conserving pseudopotentials of the Troullier-Martins type29 and a kinetic energy cut-off of 120 Ry for the expansion of Kohn-Sham orbitals. Periodic boundary conditions have been applied in all computations concerning neutral species, while open boundary conditions30 have been used to treat systems with a net electrostatic charge. Dispersion interactions have been included using the empirical approach of Ref. 31. Molecular dynamics on the ab-initio potential energy surface has been used to optimise structures, mainly by short (∼ 10 ps) simulated annealing runs. Vibrational eigenvalues and eigenvectors have been determined by diagonalisazion of the dynamical matrix (i.e., the Hessian with mass factors), computed by the simplest finite-difference approximation to second derivatives. Vibrational frequencies, in particular, have been used to compute the zero-point energy of all species, whose contribution changes reaction energies by a few hundredth of an eV.

5



(a)

Page 6 of 41

(b) (c)

FIG. 1: Ground state geometry of the gas phase cations. (a) [DAEt]+ ; (b) [DAPr]+ ; (c) [DABu]+ .

III. A.

RESULTS Density functional computation of geometry and hydrogen bonding in neutral ion pairs

and small aggregates.

To quantify energy and structural properties of half-neutralised diamine salts we carried out density functional computations for single ions, neutral ion pairs, small neutral aggregates with up to eight ions, and molecularly thin hydrogen bonded wires. Three different cations have been considered, consisting of short aliphatic segments terminated by a pair of complementary amine-ammonium groups [H2 N-(CH2 )n -NH+ 3 , n = 2, 3, 4] in a bola configuration. To fix notation, we refer to the ethyl (n = 2) case as [DAEt]+ , to the propyl case (n = 3) as [DAPr]+ , to the butyl case (n = 4), as [DABu]+ . Each of these ions has been investigated alone and in combination with the trifluoromethanesulfonate, i.e., triflate anion, CF3 SO3 − , or [Tf]− . The ground state geometry of the three gas-phase cations optimised at the DFT-PBE level by a short simulated annealing is reported in Fig. 1. Basic geometrical parameters are listed in Tab. 1. The three structures show a tight pairing of the amine-ammonium group, with a bridging proton in between, forming a kind of ring configuration. The [DAEt]+ structure appears to be fairly strained, and the N-H—N combination would not be considered a proper hydrogen bond according to standard criteria used to analyse MD trajectories, since the N-H—N angle is 129◦ only. The strain is much less in [DAPr]+ and [DABu]+ , whose N-H—N angle rapidly approaches linearity, making the N-H— N combination a genuine hydrogen bond. The increasing strength of the HB with increasing n is accompanied (or, one could say, emphasised) by the increasing length of the covalent N-H bond, and decreasing length of the non-covalent H—N bond. The n-dependence of the N · · · N distance is

6


Page 7 of 41 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: Geometric parameters and N-H stretching frequency for the intra-ion H-bond in the gas-phase cations. Distances are in ˚ A , the angle is in degrees and the stretching frequency is in cm−1 . Compound

d(N-H)

d(N...N)

d(H—N)

ˆ − − N) (N − H −

ω(N − H)

[DAEt]+

1.10

2.56

1.71

129

2465

[DAPr]+

1.15

2.61

1.53

153

1850

[DABu]+

1.18

2.63

1.48

163

1721

(b)

(a)

FIG. 2: Lowest energy configuration of: (a) [DAEt][Tf]; (b) [DABu][Tf]. Dash lines represent hydrogen bonds.

somewhat counter-intuitive, since it grows with increasing n, while the HB appears less strained and thus more stable. The longer N· · ·N distance, however, is simply the geometric consequence of the increasing linearity of the N-H—N bond. The strengthening of the intra-cation HB with increasing carbon content is decisively confirmed by the estimate of the N-H stretching frequency, that decreases monotonically and significantly from [DAEt]+ to [DABu]+ . The frequency of this N-H stretching mode 7



Page 8 of 41

is low (at 2465 cm−1 ) already in [DAEt]+ . In [DAPr]+ and [DABu]+ the N-H stretching is no longer distinct from H-N-H bendings, but the highest frequency mode with significant stretching character is computed at 1850 cm−1 in [DAPr]+ and at 1721 cm−1 in [DABu]+ (see Tab. 1). Pairing each of these cations with the [Tf]− anion gives origin to very stable neutral ion-pairs, whose structure is illustrated by the ground state geometry of [DAEt][Tf] shown in Fig. 2 (a) and [DABu][Tf] in Fig. 2 (b). The ground state structure of [DAPr][Tf], shown in Fig. S3 in SI, is intermediate between these two. In all these pairs, Coulomb forces apparently play the major role, as emphasised by the close contact between the ionic moiety of [Tf]− , represented by -SO− 3 and the corresponding ionic moiety of the cation, represented by the -NH+ 3 termination, and, to a lesser extent, by -NH2 . Hydrogen bonds play a secondary but still important role. The intra-cation HB, in particular, persists also in the ion pairs, while the H2 N–NH+ 3 terminations act like a bidentate functional group donating two protons to [Tf]− , forming two HBs of unequal strength, the strongest originating from -NH+ 3 , the weakest from -NH2 . The geometric parameters defining the strongest of these two inter-ion HBs are collected in Tab. 2. The primary role of Coulomb forces and the complementary role of H-bonding is confirmed by similar computations for larger aggregates, up to four neutral ion pairs. In all cases, at T = 0 K ions arrange themselves on a somewhat distorted but easily recognisable NaCl-type lattice. Also in the case of clusters, cations extend at least two HBs to different neighboring anions, giving origin to a connected network of HBs. This is exemplified by the ground state structure of ([DAEt][Tf])4 illustrated in Fig. S4 of SI. The energetics of ionic neutral ion pairs is summarised by the reaction energies that follow, computed by DFT in the PBE gradient corrected approximation. HTf → [TF]− +p

∆E = −13.57 eV

DAEt+ p+ → [DAEt]+

∆E = 10.19 eV

[DAEt]+ + [TF]− → [DAEt][TF]

DAPr+ p+ → [DAPr]+

∆E = 5.08

∆E = 10.62 eV

[DAPr]+ + [TF]− → [DAPr][TF]

DABu+ p+ → [DABu]+

∆E = 5.09 eV

∆E = 10.76 eV

[DABu]+ + [TF]− → [DABu][TF]

∆E = 4.90 eV

where positive energies mean that the reaction runs spontaneously from left to right. 8


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


Although the triflic acid HTf is one of the strongest acids (a super-acid), its deprotonation to give the triflate anion [Tf]− requires a sizeable 13.57 eV, only partially offset by the energy gain (about 10 eV) in protonating the diamine to form the cation. A positive energy balance is restored only by pairing the two ions into neutral aggregates, releasing about 5 eV. The overall stability of the ionic pairs is measured by the reaction energy (T = 0K including zero-point energies, zpe): HTf + DAEt → [DAEt][TF]

∆E = 1.47 eV

HTf + DAPr → [DAPr][TF]

∆E = 1.76 eV

HTf + DABu → [DABu][TF]

∆E = 1.80 eV

Zero-point energy contributions are not negligible, their effect being to shift downwards ∆E, decreasing it from the estimate without zpe: ∆E = 1.70 eV for [DAEt][Tf]; ∆E = 2.14 eV for [DAPr][Tf]; ∆E = 2.09 eV for [DABu][Tf]. The slow increase of stability of the ion pairs with increasing length of the carbon chain, therefore, results from the compensation between the increasing energy gain through protonation of the diamine and the decreasing gain in joining cation and anion. Both these competing effects can be seen as primarily due to hydrogen bonding. All inter-ion HB are relatively unstrained (as shown by the N-H—O angle, see Tab. 2), but their loss of cohesion from [DAEt][Tf] to [DABu][Tf] is apparent, and likely to result from the competition of the intra-cation bond, whose antagonistic role is highlighted by the comparison of the geometric data in Tab. 1 and Tab. 3. Tab. 1, in particular, reports parameters of the intra-cation HB computed for the gas-phase cation, while Tab. 3 reports the same parameters computed for the neutral ion-pair. In Tab. 3 the non-monotonicity in the N · · · N distance as a function of n is likely to be due to the interplay with the secondary HB joining the cation and the anion in the neutral ion-pairs. Comparison of the ESP charges computed for all the ions and aggregates described until now and for different isomers of selected species, reveals a variety of effects not always unambiguously due to electrostatic polarisation. A detailed analysis, in fact, points to a sizeable charge redistribution within each ion due to changes in the geometry of bonds, or to the proximity of chemical groups, giving origin to a chemical polarisation effect in addition and beyond electrostatic polarisation. For the systems investigated in the present study, the effect is particularly relevant for the cations, due to the flexibility of their structure. Quantitative examples of these effects are reported in Fig. S6 of SI. The fixed charges of the rigid ion model do not account for these effects. Moreover, to avoid spurious 9



TABLE 2: Geometric parameters for the primary H-bond joining cation and anion. A weaker HB is also present. Distances are in ˚ A and the angle is in degrees. Compound

d(N-H)

d(N...O)

d(H—O)

ˆ − − O) (N − H −

[DAEt][Tf]

1.12

2.58

1.47

171

[DAPr][Tf]

1.11

2.58

1.48

170

[DABu][Tf]

1.09

2.66

1.58

170

TABLE 3: Geometric parameters for the intra-cation H-bond measured on the ground state geometry of the neutral ion pair. Distances are in ˚ A and the angle is in degrees. Compound

d(N-H)

d(N...N)

d(H—N)

ˆ − − N) (N − H −

[DAEt]+

1.05

2.66

1.94

122

[DAPr]+

1.05

2.74

1.91

133

[DABu]+

1.10

2.69

1.65

157

5

4

∆ U [kcal/mol]

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 10 of 41

3

2

1

0 -0.8

-0.4

0.0

0.4

0.8

x [A]

FIG. 3: Variation of the potential energy on moving the proton out of the ground state position along the direction joining the two N-atoms in the same [DAEt]+ cation. Green line: computed on the gas-phase cation. Blue line: computed on the cation joined to a [Tf]− anion in a neutral ion pair. The reaction coordinate x is the difference of the proton distances from the two N atoms.

distinctions among atoms formally equivalent such as the three protons on -NH3 or the two protons on -NH2 , the assigned charges are somewhat averaged over equivalent atoms and over low energy isomers.

10


Page 11 of 41 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


We realise that these choices and approximations might affect at least quantitatively features that are discussed in the force-field MD section, including the absolute and relative stability of different types of HBs. The basic geometric and bonding properties characterising the ground state structure and isomers of small aggregates are reflected in the energy variations upon moving the protons in HBs from their covalent to the non-covalent partner along each hydrogen bonds. The results are collected in Fig. 3 and Fig. 4. Each of the curves on these two figures has been computed by moving the proton (p) towards its non-covalently bonded partner (X) in steps of 0.04 ˚ A . At each step, the p-X distance is kept fixed by a constraint, and the energy minimised with respect to all other degrees of freedom. By symmetry, the potential energy along the reaction coordinate for the intra-cation HB always has a double minimum. The barrier between the two minima is below ∼ 1 kcal/mol for all gas-phase cations, and decreases slowly along the [DAEt][Tf], [DAPr][Tf], [DABu][Tf] sequence. The barrier is somewhat higher, but remains below 5 kcal/mol when measured along the same path in the neutral ion pairs (see again Fig. 3), partly because the N · · · N distance within the cation increases with respect to the isolated cation, and partly because the proton transfer has to be accompanied by a substantial re-arrangement of the relative orientation of the two ions. Nevertheless, both in the case of gas phase cations and of the ion pairs, the low barrier and proximity of the two minima along the potential energy curve suggest that the amine and ammonium terminations could be more symmetric than implied by their structural formula. Taking quantum effects into account, it is likely that one proton is equally shared between the two nitrogens on the same cation. This picture is compatible with the experimental 1 H NMR spectra3 which are unable to detect the presence of the free amine. The corresponding energy variations in moving the proton along the cation-anion primary hydrogen bond are collected in Fig. 4 (a). Also in this case, quantum mechanical effects are not included, although they are expected to be non-negligible. In this case, the reaction coordinate is the H — O distance, which is progressively reduced to plot the curves in panel (a) of Fig. 4. For all choices of cation, the potential energy curve has one minimum only, consistently with the short N · · · O distance.14 However, the energy needed to move the proton to a covalent bond distance (∼ 1.1 ˚ A) from the oxygen is low, and a double minimum could form whenever cation and anion move further from each other by thermal fluctuation. In particular, the energy to move the proton to about 1.1 ˚ A from the nearest oxygen in a N-H–O HB is ∼ 4 kcal/mol in [DAEt][Tf] and in [DAPr][Tf], growing slightly in [DABu][Tf] because the N · · · O separation is somewhat longer in this last compound. The computation of chemical reaction barriers by simple exchange-correlation approximations such as PBE might be affected by significant and non-systematic errors.32 To assess the uncertainty of the 11



8

O 6

∆ U [kcal/mol]

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 12 of 41

+ N H x

H

x

O S

H H N+

O

N H

H

4

2

(b)

(a) 0 1.6

1.4

1.2

1.6

1.4

1.2

x

x

FIG. 4: Energy variation on moving the proton out of the equilibrium position along the direction leading the proton to the non-covalently bonded partner in inter-ion HBs. Panel (a): primary cation-anion HB. The reaction coordinate x is the distance between the proton on: [DAEt]+ (blue line); [DAPr]+ (green line); [DABu]+ (red line) and the oxygen on [Tf]− . Panel (b): cation-cation HB in the model [DAEt][Tf] wire (see text). The coordinate x is the distance between the proton from -NH+ 3 and the nitrogen on the neighbouring NH2 -. Red line: single proton displacement. Blue line: in-phase displacement of all equivalent protons. In the latter case, the energy per moving proton is displayed.

PBE barriers computed in our study, we repeated the barrier determination using the hybrid B3LYP functional.33 Although this method might not be the best one among the many hybrid functionals that are available,34 its Hartree-Fock component represents a significant difference with respect PBE, and the discrepancy between the two barriers is expected to give a fair estimate of the uncertainty on the results we reported. In the case of the intra-cation proton transfer in neutral [DAEt][Tf] discussed above and illustrated in Fig. 3 (green line), the B3LYP barrier turns out to be ∆E = 5.2 kcal/mol,to be compared to the ∆E = 4.46 kcal/mol given by PBE. Although significant, this difference of less than 1 kcal/mol does not alter the general picture of hydrogen bonding and proton transfer dynamics given by PBE. This assessment is reinforced by considering that the zero-point energy of the proton in between the two nitrogens of [DAEt][Tf] (ω = 2565 cm−1 ) is already nearly 4 kcal/mol. The full energy profile for the proton transfer computed for this case using PBE and B3LYP is shown in Fig. S5.

12


Page 13 of 41 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


B.

Proton conductivity in model [DAEt][Tf ] wires

The positive outlook for proton conduction resulting from the low barriers towards proton displacement along single hydrogen bonds is tempered by the observation that Grotthuss transport across the system eventually requires proton jumps from cation to cation. To make this mechanism relevant, however, cations need to form chains like those schematically depicted in part (a) of Fig. 5, neutralised by a corresponding number of anions, as already suggested in Ref. 3. The type of hydrogen bond stabilising these chains is apparently less stable than those linking cations and anions, since the former join functional groups of comparable positive charge (assuming that the two terminations are nearly equivalent), while the latter join groups of opposite charge, thus enjoying a sizeable stabilisation by electrostatic forces. This qualitative picture is confirmed by the results of our simulations of nano-clusters, in which we observe a number of very stable cationanion HB, but no cation-cation in the ground state structure. Cation-cation HBs, however, occur in the ab-initio simulation of the [DAET][Tf] tetramer at room temperature, suggesting a stabilisation contribution from entropy. An example of such a bond is shown in Fig. S7 in SI. Even the intracation HB appears to be significantly more stable than the cation-cation link, although it joins the same type of terminations, probably because of the forced proximity of the two groups. However, in condensed phases, inter-cation hydrogen bonds might be stabilised by close correlation with anions, hence, estimating the stability of these bonds and quantifying energy barriers for proton transfer becomes a relevant task. To this aim, an artificial 1D system has been prepared, consisting of four [DAEt][Tf] ion pairs on a linear arrangement, periodically repeated in their longitudinal direction to form an extended wire (see Fig. 6). Because of computational convenience, the system is replicated also in the two transversal directions, but in this case the periodicity is such to prevent sizable interaction among replicas. Because of the topology of cations, able to accept a single HB at its -NH2 termination, chains are the most likely type of cation-cation aggregation in these protic RTILs, analogous but geometrically distinct from those discussed, for instance, in Ref. 35. The simultaneous acceptance / donation of HBs by each cation provides the crucial degree of cooperativity. Branching of chains in principle is possible, because of the ability of -NH+ 3 to donate more than one proton, but in that case cooperativity might act in reverse, weakening multiple HBs from the same -NH+ 3 group. The system prepared in this way turns out to be at least metastable, with an optimal N · · · N distance of 2.72 ˚ A . However, to prevent large deformations, or even the collapse of the wire into a compact cluster under the perturbation driving the proton transfer, a weak parabolic restrain has been

13



Page 14 of 41

(a) N

(2)

N N

(2)

(3)

N (3)

(b) N (3) N

N (2)

(2)

N (3)

(c) N (2)

N

(3)

N N (2)

(3)

N (3) N (2)

FIG. 5: Chain-like arrangement of [DAEt]+ cations (a), (b), (c), with: (b) bridging by [Tf]− and (c) bridging by [Tf]− and by water. Possible paths for proton jumps are indicated by dash lines and arrows. Blue dots: N; red dots: O; yellow dots: S. For clarity, the -CF3 moiety of [Tf]− is omitted.

FIG. 6: Linear chain built to investigate the stability and energy profile of hydrogen bonds joining cations.

14


Page 15 of 41 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


applied to the carbon atoms in the chain, keeping them close to the position found by the starting local optimisation of the structure. Remarkably, the lowest energy wire configuration we found has all the cations (all the anions), on the same side (on the opposite side) of an ideal line marking the longitudinal direction of the wire. A different geometry, with ions arranged in such a way to cancel the 2D dipole moment in the transversal directions is found to have higher potential energy. From this starting point, we carried out the same exploration of the potential energy along the reaction coordinate as it has been done for the other types of HB. In a first stage, a single proton has been displaced along the reaction coordinate. Once again, the results show that the potential energy curve has one minimum only (see the red line in Fig. 4, panel (b)), but a moderate amount of energy (∼ 5 kcal/mol) is sufficient to move the proton to within 1.1 ˚ A from the non-covalently bonded N atom. The absence of a second minimum on the potential energy curve (despite the equivalence of the two N-terminations) is not surprising, since moving a single proton creates a neutral DAEt molecule and a doubly charged [DAEt]++ cation, whose energy would be fairly high even if screened by a corresponding number of neighboring [Tf]− . A double energy minimum is recovered if all the (equivalent) bridging protons in the chain are moved together. The result of this computation is shown by the blue line in Fig. 4, panel (b). The curve is not symmetric with respect to the midpoint of the line joining the two nitrogen atoms because the local optimisation routine used to minimise the energy is unable to reach the minimum at the same energy of the original one, which, by symmetry, has to exist. To be precise, the plot reports the variation per moving proton of the potential energy along the reaction coordinate. Hence, the barrier opposing the shift of n protons quickly overcomes the thermal energy with increasing n, even though quantum effects decrease significantly the energy per moving proton. Matching the two pictures, we expect the propagation of protons by a soliton-like excitation, as the one postulated in Ref. 36 for the bi-squaric acid. These results can be summarised as follows. Density functional computations suggest mechanisms for the coherent transport of protons along chains of cations, compatible with an accessible energy scale. The precise mechanism and quantitative barrier could be determined by density functional computations at least for rather idealised model geometries, but such a complex undertaking is beyond the scope of our study. The major task left to molecular dynamics simulation, therefore, is to verify the presence of these chains, and to estimate their length and the frequency of their occurrence as a function of temperature.

15



C.

Page 16 of 41

Molecular dynamics simulations

Simulations based on the classical force field of Sec. II have been carried out for samples of [DAEt][Tf] and [DABu][Tf] consisting of 216 neutral ion pairs (4536 atoms for [DAEt][Tf], 5832 atoms for [DABu][Tf]) in a cubic box of ∼ 4 nm side with periodic boundary conditions applied. Most of the effort has been devoted to [DAEt][Tf] while a few simulations on [DABu][Tf] have been carried out for a comparison, aiming to assess the effect of the length of the aliphatic chain joining the two nitrogen terminations. All simulations have been performed at ambient pressure. The experimental melting temperature of [DAEt][Tf] is Tm = 351K, and its thermal stability extends up to 500 K, as experimentally verified by the absence of weight loss during thermogravimetry.3 In the [DAEt][Tf] case, simulations started in the liquid range at T = 440 K, then the temperature was progressively decreased in steps of 10 K down to T = 200K, well below the experimental melting temperature. At each temperature the [DAEt][Tf] sample has been equilibrated during 6 ns, and statistics accumulated over 200 ns. Thus, over the entire sequence of [DAEt][Tf] simulations from T = 440 K down to T = 200 K, the accumulated statistics covers 5 µs, for an average annealing rate of 48 × 106 K/s, not counting the equilibration runs of 6 ns at each change of temperature. In the case of [DABu][Tf], only four temperatures have been considered, i.e., T = 440 K, 400 K, 350 K and 300 K. To optimise equilibration, also in this case we moved from the highest to the lowest temperature, each time collecting 200 ns of statistics. Equilibration, however, has been longer (16 ns) because of the larger temperature variations at each stage. The average density hρ(T )i of [DAEt][Tf] at T = 400K turns out to be ρ = 1.44 g/cm3 to be compared with the experimental value ρexp = 1.6 g/cm3 . The 10 % difference is the same as the one found in Ref. 15 for [TEA][Tf], and corresponds to a 3 % overestimate of the average side of the cubic simulation box. For molecular species, and for RTILs in particular, an error of this size is not unusual even in DFT computations.37 Moreover, the discrepancy is comparable to the spread of experimental values in closely related ionic liquid systems. The density of [TEA][Tf] at T = 400 K, for instance, is quoted as ρ = 1.33 g/cm3 in Ref. 3 and as ρ = 1.16 g/cm3 in Ref. 38. The density of [DABu][Tf] computed at T = 400 K is ρ = 1.31 g/cm3 . To the best of our knowledge, no experimental value is available for a comparison. Despite the hint of a curvature, the hρi(T ) of [DAEt][Tf] is recognizably linear both at low (200 ≤ T ≤ 280 K) and at high T (330 ≤ T ≤ 440 K) (see Fig. 7(a)). Over the 280 ≤ T ≤ 330 K interval, hρi(T ) crosses over these two linear ranges of slightly different slope, with only a weak anomaly that is highlighted in the inset of Fig. 7(a) by subtracting from hρi(T ) the linear fit to the T ≥ 380 K

16


Page 17 of 41 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


data. The transformation has a weak first order component representing a largely incomplete ordering process, superimposed to what appears to be a predominantly glass transition. No long range order in the ion configuration can be detected by visual inspection of snapshots. This picture is likely to be determined by genuine properties of the material such as the high viscosity, but certainly it is also affected by the small system size and high cooling rate. The temperature dependence of the density is reflected in an analogous behaviour of the average potential energy, shown in Fig. 7(b), with the broad dip in the inset of panel (a) becoming an equally broad peak in the inset of panel (b). More detailed knowledge on the structure and cohesion is revealed by the structure factor Sαβ (k), defined as: Sαβ (k) =

1 hρα (k)ρβ (−k)i N

(1)

where α, β label atomic species, N is the total number of ions (twice the number of molecules), and: ρµ (k) =

X

e−ik·rj

(2)

j∈µ

rj being the position of ion j belonging to species µ. Since we aim at highlighting the ionic character of the system, we define the structure factor of a distribution of cations and anions, whose location is identified by the position of the ammonium nitrogen, and by the sulphur atom of [Tf]− . The partial structure factors S++ (k), S+− (k) and S−− (k) (see the definition in Ref. 39) are combined into a density-density Snn (k) and a charge-charge SZZ (k) structure factor (see again Ref. 39) as: Snn (k) =

1 [S++ (k) + 2S+− (k) + S−− (k)] 2

(3)

SZZ (k) =

1 [S++ (k) − 2S+− (k) + S−− (k)] 2

(4)

The results, displayed in Fig. 8 (a) and (b) for T = 350K, show a prominent peak in the charge-charge structure factor at k = 1.05 ˚ A−1 , emphasising the crucial role of Coulomb interaction driving the strict alternation of cations and anions. In simple ionic liquids, the basic periodicity of charge is roughly twice as long as the periodicity of mass, since charge repeats itself at the next nearest neighbour distance (in a NaCl-type structure). Hence, the major peak in the density-density structure factor is ˚−1 . This is indeed the position of the highest peak in Snn (k), but Fig. 8 shows expected at k = 2.1 A ˚−1 , k = 1 A ˚−1 , as well as a sizeable shoulder at k ∼ 0.6 A ˚−1 . The clear pre-peaks at k = 1.45 A last two pre-peaks, in particular, point to regularities in the same-ion position of periodicity 6 and 10 17



-20

〈ρ〉 −γ−δT [A-3]

〈E〉 −α−βT

0.003

1.6

-2.0 200 260 320 380 T

[K]

0.0 200 260 320 380 440 T

1.5

[K]

〈ρ〉

[kJ/mol]

-60

0.006

0.0 -0.5 -1.0 -1.5

[A -3 ]

[kJ/mol]

1.6

-40

〈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 18 of 41

-80

1.5 -100

(a) -120 200

260

320

T

(b) 380

440

1.4 200

260

[K]

320

T

380

440

[K]

FIG. 7: Average potential energy hEi and average density hρi as a function of temperature from MD simulations of [DAEt][Tf] at NPT conditions. In both panels the inset reports the plotted quantity minus its linear interpolation for T ≥ 330 K. Data collected on cooling.

˚ A , beyond the nearest neighbour distance, that might be interpreted as due to the formation of a mesophase. This feature is not unusual in RTILs,40–44 and could be interpreted in terms of the interplay of charged and neutral moieties in both ions. In other terms, to optimise both Coulomb and dispersion interactions, charged moieties cluster in tight configurations surrounded by the non-polar moieties. Depending on the relative size of these charged and neutral moieties, different length scales arise, giving origin to the rich structure of these systems, as discussed, for instance, in Ref. 45. In protic ionic liquids such as the systems investigated here, hydrogen bonding contributes to the nanostructuring of the liquid compound, as apparent from our results, and as already discussed in the literature.35,46 Nanostructuring in [DAEt][Tf] might also be enhanced by the low polarity of the fluorinated group in [Tf]− ,47 although in this case fluorination is reduced to its bare -CF3 minimum. The nanostructuring picture is confirmed and reinforced by the comparison with the results for [DABu][Tf], shown in Fig. 8 (b). The longer saturated chain of [DABu]+ reduces significantly the first peak of SZZ (k), emphasising the decreasing role of Coulomb interactions. At the same time, the pre-peaks in the density-density Snn (k) grow in size, and move to slightly lower k, because of the general expansion of the system caused by the insertion of the additional -CH2 -CH2 - segment. Apart from minor details, the structure factor of [DAEt][Tf] and [DABu][Tf] displays little dependence on T over the interval covered by our simulations, showing only the expected but slow increase in the amplitude of the Snn (k), SZZ (k) oscillations from T = 440 K down to T = 200K, while the position of peaks and dips remains the same. The complexity of the relative arrangement of ions at short and medium range is emphasised also by the radial distribution functions, reported in Fig. 9, showing even more clearly than the structure

18


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


Snn(k), SZZ(k)

Page 19 of 41

SZZ(k)

(a)

2

SZZ(k)

(b)

Snn(k)

Snn(k)

1

T=350 K

T=350 K

0 0

1

2

3

0

1

k [A-1]

2

3

4

k [A-1]

FIG. 8: Density-density and charge-charge structure factors at T = 350 K of: (a) [DAEt][Tf], and (b) [DABu][Tf]. The experimental melting point of [DAEt][Tf] is Tm = 351 K.3

factor the interplay of different length scales in the local and medium-range structure of the system. Also in this case, we compute and report the radial distribution of cations and anions, identified with − the position of nitrogen in -NH+ 3 and sulphur in [Tf] . Apparent features are the very prominent first

peak of g+− (r) at ∼ 3 ˚ A , the relatively simple behaviour of g−− (r), and the broad, nearly structure-less shape of g++ (r), due to the complex geometry of the cation. The shape of the second peak of g+− (r) at 6 ˚ A is the most apparent signature of nanostructuring in the radial distribution function. In simple molten salts such as NaCl close to their triple point, the second peak of g+− (r) reaches a value well above one.48 In [DAEt][Tf], instead, the second peak remains below one (g+− (r = 6) = 0.96), reflecting the depletion of the second cation-anion coordination shell, contradicting the system homogeneity at all length scales. These observations suggests a 12 ˚ A diameter for the average blob made of charged moieties, consistently with the pre-peak in S+− (k) at ∼ 0.6 ˚ A

−1 .

In [DABu][Tf], such a depletion

of the second cation-anion shell is more pronounced, with (g+− (r = 6) = 0.62), again consistent with the enhanced role of the mesophase in this compound, as revealed by the structure factor. The first peak of g+− (r) around a central -NH+ 3 integrates to a coordination number by S atoms of nearly 4, corresponding to a compact coordination of nearly 12 by oxygen atoms. Needless to say, the alkane chain, the neutral -NH2 termination of [DAEt]+ , and the (equally neutral) -CF3 termination of [Tf]− are virtually excluded from this ionic core. Analysis of trajectories also show that in the [DAEt]+ case, up to T = 440 K, the N-C-C-N backbone remains in the cis configuration, hence the two terminal N atoms of each cation remain within a HB distance (rc = 3.2 ˚ A ) from each other. The strain of the [DAEt]+ structure, however, prevents the alignment of the N-H—N combination to within ∼ 30◦ as required by most empirical 19



4

g++(r), g+-(r), g--(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

Page 20 of 41

g--(r)

3

g--(r)

g+-(r)

T=350 K 2

g+-(r)

T=350 K

g++(r)

g++(r)

1

[DAEt][Tf]

[DABu][Tf]

0 2

4

6

8

10

2

4

r [A]

6

8

10

12

r [A]

FIG. 9: Radial distribution functions of ions in [DAEt][Tf] (left panel) and in [DABu][Tf] (right panel).

criteria to identify a HB. Hence, no formal intra-cation HB is found along the trajectory, although these atoms are clearly bound by a sizeable attractive force, as confirmed by the DFT computations for the gas-phase cation. In what follows, the interaction among these atoms on the same cation is loosely referred to as a HB. As expected, the picture is different in [DABu][Tf], since in that case the gain of entropy made available by the longer chain overcomes the enhanced stability of the HB, and a majority of intra-cation HB opens up, being replaced by further HB to neighbouring ions (see Fig. S9 in SI). Thus, these results from simulation help addressing the question raised in Ref. 3 concerning the intraor inter-ion chelation of protons, which cannot be unambiguously answered by NMR measurements. Simulation confirms the guess proposed in that same reference, stating that the balance of intra- and inter-ion chelation depends on the length of the spacer chain. In our results, chelation is exclusively intra-cation in [DAEt][Tf], and predominantly but not exclusively intra-cation in [DABu][Tf]. The major reason of interest in these materials and in our simulations is the structure, energetics and topology of the network of HBs. Each cation forms, on average, more than one HB to neighboring anions, as can be seen in Fig. 10. Moreover, and perhaps more importantly, the multiple HBs centred on each cation reach out to different anions. Hence, the analysis of the HB network shows that, almost without exceptions, each sample represents a uniquely HB-connected set. This high connectivity of the HB network is certainly one of the reasons of the relatively slow diffusion of ions and high viscosity of these compounds.46 Here, as before, H-bonding is defined in terms of geometrical parameters only, requesting a distance between the electronegative pairs (N · · · N or N · · · O) of less than rc = 3.2 ˚ A , and a N-H—N angle deviating less that 30◦ from linearity. As expected, the average number of cation-anion (N-H—O) HBs increases with decreasing temperature, and the relative change over

20


Page 21 of 41 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 200 ≤ T ≤ 440 K range is higher than in the [TEA][TF] case.15 Together with strong Coulomb interactions, the network of HB is one of the reasons of the clustering of ionic moieties in the [DAEt][Tf] and (even more) [DABu][Tf] equilibrium liquid structure, giving origin to the pre-peaks in Snn (k) already discussed. As discussed in the ab-initio sub-section, hydrogen bonding could also join amine and ammonium nitrogens on neighbouring cations. According to our simulations, and using the geometric criteria listed above, the number of these HB in [DAEt][Tf] is low, although not completely negligible. This observation affects the number of paths open for proton conductivity through a Grotthuss-like mechanism. It is important to remark that these results from simulation might also be affected by limitations in the force field model, which treats the amine and ammonium terminations as distinct, while ab-initio computations suggest they are nearly equivalent. Moreover, the polarisability of -NH2 , not included in the force field, could also stabilise its H-bonding to -NH+ 3 . The quantitative analysis of paths open for proton conduction reported in the last part of this section will take these considerations into account. The long simulation trajectories allow a precise determination of diffusion properties down to low temperature. The diffusion coefficient of cations and anions is computed using Einstein’s relation, i.e., from the slope of the mean-square displacement per ion over the interval 60 ≤ t ≤ 100 ns, where the mean square displacement is apparently linear in time. As in [TEA][Tf] (see Ref. 15) and in agreement with experiments, cations diffuse faster than anions, but the difference is not large, and p might simply reflect (the square root of) the ratio of cations and anions ( 61/149, masses given in amu). As for many RTIL systems, the simulation value for the diffusion coefficient underestimates significantly the experimental value. At T = 400K, for instance, our computations give D+ = 3.13 10−7 cm2 /s, D− = 2.74 10−7 cm2 /s for cations and anions, respectively, while at the same temperature the experimental values3 are D+ = 1.36 10−6 cm2 /s, D− = 1.18 10−6 cm2 /s. Remarkably, the experimental value for the ammonium proton mobility is higher than the diffusion coefficient of the whole cation, emphasising the presence of Grotthuss mobility. The systematic underestimation of diffusion coefficients is a well know limitations of rigid-ion force fields, and could be corrected by rescaling the charges of cations and anions, for instance, to a total of ±0.8e.49 This corrections is not adopted here because it would spoil other properties, such as the average density or the electrical conductivity. A more systematic improvement could be provided by including polarisable force fields, but their usage is computationally much more demanding. More importantly, as already pointed out in Sec. II, charge polarisation in these systems might have a chemical origin besides the electrostatic one, due to structural isomerisation or to the change of nearest neighbours. Effects of this type can be included in bond-order force fields more naturally than in traditional polarisable force fields. 21



2.5

2.0

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

Page 22 of 41

1.5

1.0 200

240

280

320

360

400

440

T [K]

FIG. 10: Average number of N-H—O HB per cation (or, equivalently, per anion) as a function of temperature. HBs are defined by N · · · O separation ≤ 3.2 ˚ A and N-H—O angle ≥ 150◦ . The line is a guide to the eye.

An Arrhenius plot of the diffusion coefficients as a function of the inverse temperature shows a cross-over between two different linear regimes taking place at T = 320K (see Fig. 11). At high (d)

temperature we estimate an activation energy for diffusion of Ea = 4500 K. At low temperature the (d)

activation energy turns out to be Ea

= 2100 K. We quote only one value for cations and anions

because the difference cannot be resolved. The approximate Stokes-Einstein relation, already used to analyse dynamical data of ionic liquids,50 provides a way to estimate viscosity. We use the relation in the form: η=

kB T 6πDr

(5)

where η is the viscosity coefficient, D = (D+ + D− )/2 is the diffusion coefficient averaged over species, 6π depends on the choice of boundary conditions (no slip) for the viscous flow, and r measures the max of the first peak hydrodynamic radius of diffusing particles. We estimate r from the position r+− max /2. In this way, for [DAEt][Tf] we estimate η = 48 cP (≡ mPa · s) at in g+− , setting r = r+−

T = 400 K, η = 224 cP at T = 350 K, η = 24.8 103 cP at T = 260 K. Since they are obtained from the approximate Stokes-Einstein relation, and, moreover, the definition of r is uncertain, these values represent only an order of magnitude estimate. Nevertheless, the value at T = 400 K is within a factor of two from the experimental value ηexp = 27 cP.3 More in general, we observe that the computed values are consistent with the generally high viscosity of these compounds, and, at low temperature, comparable to those of glassy molecular systems as obtained by MD simulation. Before discussing paths of HBs suitable for Grotthuss proton transport, we collect here our results 22


Page 23 of 41

-6 10

[cm2/s]

D-8

D+

10

D+, 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


-10 10

2.5

3.0

3.5

1000/T

4.0

4.5

5.0

[K-1]

FIG. 11: Arrhenius plot of the diffusion coefficient of cations and anions in [DAEt][Tf]. Full lines represent the linear interpolation to the low-and high-temperature data for [DAEt]+ .

for a few other dynamical properties, such as the electrical conductivity and the rotational relaxation. The duration of our simulated trajectories, although long, is not sufficient to provide an estimate of conductivity of quality comparable to that of diffusion. The relatively poor convergence of conductivity has simple statistical mechanics motivations, and the problem is well known in simulation, especially for molecular ionic liquids.51 Nevertheless we plot our best estimate of conductivity, computed using an Einstein-type relation, based on the computation of the time-correlation function: 2 + * 1 1 1 X 1 lim Π(t) = lim κ= qi [ri (t + t0 ) − ri (t0 )] kB T V t→∞ 6t kB T V t→∞ 6t i

(6)

t0

The increased statistical uncertainty forced us to extract the asymptotic linear part of Π(t) on the time interval 40 ≤ t ≤ 60 ns. The value computed for [DAEt][Tf] at T = 400 K is 5.25 10−4 S/m to be compared to the experimental value κ = 3.5 10−4 S/m. The good agreement of computed and measured electrical conductivities from a simulation reproducing only qualitatively the diffusion dynamics is likely to reflect the compensation of errors. An Arrhenius plot of κ (see Fig. 12) shows a nearly linear dependence of log (κ) versus 1/T , with a change of slope possibly corresponding to the one seen in the diffusion coefficient. A simple estimate of electrical conductivity is provided by the Nernst-Einstein equation, using the ion mobility expressed by their diffusion coefficient, and relying on the assumption that the motion of

23



-2 10 -3 10

κ [S/m]

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 24 of 41

-4 10 -5 10 -6 10 -7 10

2

3

4

1000 / T

5

[K-1]

FIG. 12: Arrhenius plot of the electrical conductivity. Full dots: simulation results from Eq. 6. Empty squares: Nernst estimate of conductivity, Eq. 7. The straight lines represent the linear interpolation of the high and low temperature portions of the conductivity data.

ions is uncorrelated: κN E =

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

(7)

In this relation, 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. Comparison of κN E and κ allows us to estimate the pairing among ions, measured by the parameter ∆ = 1 − κ/κN . In [DAEt][Tf] pairing reaches up to 0.95 at T = 200 K, indicating the expected tight pairing of cations and anions at low temperature. Up to T = 440 K, ∆ does not decrease below ∆ = 0.8. As a comparison, the same parameter ∆ computed for [TEA][Tf] shows a clear change from tight pairing (∆ ≥ 0.8) to looser pairing (∆ down to less than 0.5) around T = 320K. The ∆ value estimated in the same way from the experimental diffusion and electric conductivity coefficients of [DAEt][Tf] is about 0.7 at T = 400 K,3 in semi-quantitative agreement with the computed values. The diffusion coefficient of cation and anion is higher in [DABu][Tf] than in [DAEt][Tf] (see Tab. 4), despite [DABu]+ being slightly heavier than [DAEt]+ . This trend is expected, since the melting temperature Tm decreases from [DAEt][Tf] to longer-chain diamine salts. The melting temperature of [DABu][Tf] is not known, but, for instance, Tm decreases by 42 K in going from [DAEt][Tf] to [DABu][Tf], and by 30 K from [DAEt][Tf] to [H2 N-(CH2 )10 -NH3 ][Tf]. The enhanced fluidity appar-

24


Page 25 of 41 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


ently is the reason also for the increase in electrical conductivity from [DAEt][Tf] to [DABu][Tf], despite the dilution of the charged moieties in this last compound. The lifetime of the HBs that join cations to anions in [DAEt][Tf] has been estimated by first identifying all HB in the system at time t0 , then following their persistence in time, and averaging over the initial time t0 . As before, the definition of an HB is purely geometrical, requesting a N · · · O distance less than 3.2 ˚ A and a N-H—O angle deviating less than 30◦ from linearity. Short-time fluctuations in the distance or angle could compromise the integrity of a HB. Similarly to what ha been done in Ref. 15, the broken bond is still counted as intact provided it reforms within 0.4 ps. On the other hand, bonds that break for a time longer than 0.4 ns are no longer counted at later time, even if they do reform. The same bond, however, might contribute again to the time correlation function starting at any later time t00 . In this way, for [DAEt][Tf] at T ≥ 300 K we compute lifetimes of the order of a few ps, decreasing exponentially with increasing temperature. Slightly longer (by ∼ 20 %) lifetimes are computed in [DABu][Tf] at the same temperature. In many ways, the results are equivalent to those reported in Ref. 15 for [TAE][Tf], and certainly similar to those found in many other simulation studies52,53 of non-RTIL hydrogen bonded systems, using classical force field but also ab-initio methods. In the case of RTILs, in particular, comparable lifetimes have been computed both for protic54 ionic liquids and for water in aprotic ionic liquids.12 For these reasons, we do not elaborate any further on this property. Since the two ions in [DAEt][Tf] are anisotropic, rotational diffusion represents an additional type of motion. The orientation of the anion is relatively easy to define, and more difficult for the cation, which is less regular an also less rigid. In what follows, the orientation of the cation is defined by the vector joining the nitrogen atom in -NH2 to the complementary nitrogen in -NH3 + , while the orientation of the anion is identified by the vector joining the carbon atom to sulphur in [Tf]− , representing a kind of ternary (almost cylindrical) axis for [Tf]− . The rotational diffusion is quantified by the average: Θα (t) =

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

(8)

i∈α

where Nα is the number of ions of α type, and oi is the vector defining the orientation of ion i, identified as explained in the previous sentences. At all temperatures covered by our simulation, this time correlation function is well represented by a stretched exponential: h i Θα (t) = A exp − (t/τ )β

(9)

as shown in Fig. 13, reporting both the simulation data and their fit by the expression in Eq. 9. To be precise, the interpolation does not reproduce a very narrow feature at times of the order of a few 25



0 10

Θ - (t) Θ + (t), Θ - (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 26 of 41

Θ + (t) -1 10

τ m = 2.75ns τ p = 0.375ns -2 10 0

5

10

t

15

20

[ns]

FIG. 13: Time autocorrelation function of ion orientation in [DAEt][Tf] at T = 350 K. The full lines give the interpolation of the simulation data with a stretched exponential.

ps, that might be due to intra-ion deformations. The exponent β turns out to be almost the same for cations and anions and close to 1/2, being β = 0.54 ± 0.01 for [DAEt]+ , and β = 0.57 ± 0.03 for [Tf]− . The quoted values are the average from all the fit at different temperatures, since there is no apparent trend in β over the entire temperature range. Although there is no reason to treat the time τ as the relaxation time of a single exponential relation, we nevertheless report the temperature dependence of τ on an Arrhenius plot (see Fig. 14). From the dimensional point of view, the slope of log τ vs 1/T is an activation energy, and there is a cross-over at Tc ∼ 324 K (more precisely, Tc = 327K for τp and Tc = 320 K for τm ). The process taking place with decreasing T at Tc = 324 K is the analog of a glass transition in the rotational motion, since ions rotate fairly easily at T > Tc , while rotations are progressively frozen below this temperature, but there is no orientationally ordered phase at low T to mark a genuine plastic transition.

D.

Paths for proton diffusion according to a Grotthuss-type mechanism

The double donor-acceptor functionality of each cation opens the way to proton conductivity through a Grotthuss-type mechanism. In our classical MD simulations based on a force field of fixed topology (i.e., non-reactive) we are unable to follow the process in real time. Ab-initio (DFT) simulation, which is inherently reactive, are too expensive for unbiased simulations, and current DFT approximations, while certainly reliable for systems evolving close to their low energy thermodynamic

26


Page 27 of 41

4 10

τm

[ns]

2 10

τp

τ

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 10

-2 10 2

3

4

1000/T [K-1]

FIG. 14: Relaxation time τ for rotations, estimated from the fit of simulation data by a stretched exponential in Eq. 9. The straight lines represent the linear interpolation of the low- and high-temperature values of the relaxation times of cations and anions.

basin, are less accurate and less tested for systems undergoing chemical reactions that involve the breaking and forming of covalent bonds. For this reason, our study is limited to the investigation of correlation in atomic positions that favour the opening of path for proton transfer along chains of cations. Hence, our analysis accounts only for the system geometry, and does not say anything about the kinetics of the proton jump process, even in an approximate way. We focus on [DAEt][Tf] since experiments point to this compound as the most promising candidate for fast proton conductivity. To account for the limitations in the force field, we select fairly loose criteria to define connectivity. First, since the DFT computations show that the ammonium and amine groups in [DAEt]+ are virtually equivalent to each other, in our analysis no distinction is made with respect to the two terminations. Moreover, connectivity is defined through the N · · · N distance of nitrogens on neighbouring cations only, disregarding the orientation of all N-H bonds, assuming that hydrogens can reorient themselves relatively quickly. Then, we define the separation of two cations as the shortest distance between a nitrogen on the first cation and a nitrogen on the second cation, irrespective of their belonging to the -NH2 or NH+ 3 termination. In other terms, the separation of two cations is the shortest of four N · · · N distances. Then, chains are identified as a sequence of cations related by a nearest-neighbour condition. A schematic view of the chains of interest is shown in Fig. 5 (a). Chains are terminated when the nearest-neighbour distance is longer than 3.6 ˚ A. At each temperature, we consider 1000 configurations extracted every 0.2 ns out of the 200 ns

27



10 10

[µm-3]

[µm-3]

7

7 10

c(n)

8 10

4 10

10 1 10 2

c(n)

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 28 of 41

4

6

8

10

n

6 10

T=440 K 5 10

T=200 K 0

20

40

60

80

n

FIG. 15: Logarithm of the concentration c(n) of chains length n, expressed in number of chains per inverse cubic micron. The main panel reports the results considering the possibility of oxygen atoms on [Tf]− shuttling protons between two [DAEt]+ cations, see text. Inset: results considering only the HBs linking cations directly, see text.

production trajectories. For each configuration, we construct chains of cations, starting in turn from each [DAEt]+ , and moving to the nearest cation neighbour, with distance measured by the shortest separation of N atoms on the two cations. In this progression, the growing chain is self-avoiding, meaning that we exclude re-crossings of the chain identified up to that point. Moreover, chains are mutually non-intersecting, meaning that no cation can belong to more than one chain. The chain is terminated when no cation is within 3.6 ˚ A from the last identified one, or when the following cation is the first of the chain, giving a closed loop. The results of this first restricted analysis is summarised in the inset of Fig. 15, reporting on a logarithmic scale the n dependence of the concentration c(n) of chains of length n, expressed in µm−3 . The figure concerns only a few chain sizes (up to n = 10), since longer chains are not found in the simulated samples. Once normalised, the curves in Fig. 15 can be interpreted as the size probability distribution: c(n) p(n) = P∞ n=1 c(n)

(10)

of an equilibrium polymerisation process.5 In this case, the probability distribution is a decreasing exponential of the size n. Both the concentration of chains of length n and the maximum length of chains in the sample increase slightly with decreasing T , reflecting more a general effect of thermal 28


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


expansion than a change in the system topology. Although conceptually interesting, the result is not very relevant in view of proton conductivity, since the possibility of moving by 10 − 15 ˚ A every time a Grotthuss event occurs will not enhance proton diffusion by much, taking into account that initiating the Grotthuss process requires the formation of a coordination defect, whose concentration is not large even at relatively high T . The picture changes if we include the possibility of an oxygen from [Tf]− shuttling a proton between two successive cations along the chain. Also in this case, the chain is built by a minimum-separation requirement, but this time it includes also protons jumping from N to O and back to a nitrogen on another cation. This bridging mechanism is sketched in Fig. 5 (b). While constructing the chain, the bridging anion is retained only if the distances N · · · O and O · · · N’ on another cation are both shorter than the separation to the nearest cation, and, in all cases, distances need to be less than 3.6 ˚ A. The result of this analysis is displayed in the main panel of Fig. 15. To emphasise the difference, only the results at the two extreme temperatures T = 200 K and T = 440 K are plotted in both panels of Fig. 15. The results show that chains of up to n ∼ 80 cations (plus a comparable number of bridging anions) can be found, although at relatively low concentration for n ≥ 40. For all the chains contributing to Fig. 15, as already stated, the longest separation among nitrogen atoms is less than 3.6 ˚ A . Notice that crossing such a gap would require the proton to travel by up to 1.5 ˚ A only, and the jump could also be assisted by quantum mechanical effects. An example of a chain in [DAEt][Tf] at T = 350 K, including 12 cations and 10 bridging anions is shown in Fig. 16. Because of the selfavoiding constraint, the subdivision of the cations into chains depends somewhat on the selection and ordering of the starting points. However, the difference in the size distribution function due to this dependence is limited to quantitative details. As already stated, the temperature dependence of the restricted analysis (inset of Fig. 15) is uniform on the entire n range, pointing to a stabilisation of chains with decreasing T . The temperature dependence of the extended analysis is more complex. At n > 40, the concentration of chains increases markedly with decreasing temperature, reproducing the trend of the restricted analysis. However, since the number of ions in the sample is fixed, the trend towards larger clusters at low T has to be compensated by the relative decrease in the concentration of short chains. With the extended definition of connectivity, including the shuttling by [Tf]− oxygens, the probability distribution of chain sizes p(n) is no longer a simple exponential, but it is well reproduced by a stretched exponential (see the fit in Fig.S12 in SI). This differs from the results of equilibrium

29



Page 30 of 41

FIG. 16: Chain of 12 cations and 10 bridging anions (see text) forming a connected set in which the longest link is less than 3.6 ˚ A . Two further bridging anions have been omitted for clarity.

polymerisation theory,5 predicting a Zimm-Schultz distribution: p(n) ∝ nγ−1 exp [−γ

n ] hni

(11)

where hni is the average length of chains. The difference might also be due to the small samples of our simulations, perhaps distorting the probability distribution at n > 40. On average, one fifth of the chains form closed loops. Moreover, apart from relatively short linear segments, chains tend to be coiled, and perhaps it would be more appropriate to talk about connected clusters instead of chains. The minimum distance condition used to define chains prevents the identification of branched structures, that, however, are possible and likely. Reducing the cut-off distance from 3.6 ˚ A to 3.2 ˚ A reduces the concentration of chains by a factor of two for 2 ≤ n ≤ 40. Beyond n = 40 the relative decrease in concentration is progressively larger, limiting the range of chain sizes to n ≤ 60. All the results in Fig. 15 refer to self-avoiding, mutually non-intersecting chains. The estimate including all possible chains irrespective of their mutual intersection might also be relevant, since it gives information on the connectivity and multiplicity of the proton transport network. If we count all paths, irrespective of their crossings, the chain concentration curve extends to larger n’s, up to the size of the entire sample. For any given size n, the chain concentration c(n) growth with respect to the self-avoiding case. At n = 40, the increase is by one order of magnitude. Once again, we emphasise that the relevance of each path has to be assessed a posteriori using a kinetic model. This difficult task is left for further investigations. Nevertheless, a sizeable concentration of paths is an important feature, since the Grotthuss mechanism, in general, requires an 30


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


initiation step. In water, for instance, the initiation step is represented by the auto-ionization of the H2 O molecule, which is a rare event. The initiation step in [DAEt][Tf] could be represented by the neutralisation reaction: [DAEt][Tf] → DAEt+HTf, which, in the gas-phase requires 1.47 eV (see Sec. III A). This energy could be reduced by correlations and screening in the condensed phase. For instance, the energy of the same charge defect in the ([DAEt][Tf])2 (see Fig. S8) is already reduced to 1.27 eV. Moreover, and more importantly, some other mechanisms of lower free energy could act as the initiation step. For instance, the close approach of an anion to a cation termination could shift the HB balance at the other termination, releasing of accepting a proton. Nevertheless, it is likely that the initiation step still remains a rare event. Hence, every initiation step that naturally occurs by fluctuation has to open a long path for the proton transfer, otherwise the Grotthuss channel is likely to be quantitatively unimportant. Unfortunately, no information on the initiation step is provided by the MD simulation based on a non-reactive force field. Further ab-initio computations on systems larger than nano-clusters are needed to clarify this important point.

E.

Water effect on the structure and dynamics

A preliminary exploration of the effect of water on the structure and dynamics of [DAEt][Tf] has been carried out, paralleling the experimental investigation reported in Ref. 3, and following a closely related computational analysis we reported in Ref. 15. In the present case, we added 25 water molecules to the sample of 216 [DAEt][Tf] neutral ion pairs, corresponding to a 1% weight water contamination. The sample preparation follows the same procedure used in Ref. 15, distributing molecules at random over the system volume, quenching to remove short contacts, and equilibrating at the target temperature. Simulations have been carried out at T = 440 K, T = 400 K, T = 350 K. Each production run lasted 200 ns, following 40 ns of equilibration. The longer equilibration time was meant to allow for water diffusion throughout the sample. Additional simulations, following a different schedule, have been carried out at T = 300 K and T = 250 K, starting from T = 350 K, decreasing T in steps of δT = 10 K, each time equilibrating during 60 ns, pausing at T = 300 K and T = 250 K for a 200 ns equilibration and 200 ns statistics. Comparison of the average density computed for dry and wet samples show that at T = 300 K the volume occupied by water in [DAEt][Tf] is v = 26.2 ± 0.4 ˚ A ˚ A

3

3

per molecule, to be compared to v = 30

in bulk water at the same conditions. The 10 % shrinking underlines the good affinity of water

for [DAEt][Tf], at least at low water concentrations. Over the 300 K ≤ T ≤ 440 K range, the radial distribution function for the water-oxygens (OW-OW) has a clear peak at r ∼ 2.72 ˚ A apparently due

31



Page 32 of 41

TABLE 4: Diffusion constants and electrical conductivity of dry [DAEt][Tf] and [DABu][Tf] samples. Data for [DAEt][Tf] at 1 % weight water contamination are given for a comparison. Diffusion constants in cm2 /s; electrical conductivity in S/m. [DAEt][Tf] Dry samples Type

T = 250 K

T = 300 K

T = 350 K

T = 400 K

T = 440 K

[DAEt]+ 3.74 ± 0.5 × 10−10 3.53 ± 0.3 × 10−9 6.57 ± 0.2 × 10−8 3.13 ± 0.05 × 10−7 8.13 ± 0.05 × 10−7 [Tf]−

2.94 ± 0.4 × 10−10 2.76 ± 0.2 × 10−9 4.42 ± 0.05 × 10−8 2.74 ± 0.05 × 10−7 6.64 ± 0.05 × 10−7 6.6 ± 1 × 10−7

κ

1.2 ± 0.1 × 10−5

8.9 ± 0.6 × 10−5

5.25 ± 0.4 × 10−4 14.5 ± 1.2 × 10−4

[DABu][Tf] Dry samples Type

T = 250 K +

T = 300 K

T = 350 K −9

1.39 ± 0.1 × 10−6

8.38 ± 0.9 × 10

[Tf]−

-

6.18 ± 0.7 × 10−9 9.24 ± 0.9 × 10−8 5.25 ± 0.5 × 10−7 1.27 ± 0.1 × 10−6

κ

-

3.7 ± 0.4 × 10−5

6.6 ± 0.5 × 10−4

4.96 ± 0.5 × 10

T = 440 K −7

-

[DAEt]

1.07 ± 0.2 × 10

T = 400 K −7

11.5 ± 1 × 10−4

17.3 ± 1 × 10−4

T = 400 K

T = 440 K

[DAEt][Tf] Wet samples Type

T = 250 K

T = 300 K

T = 350 K

[DAEt]+ 3.88 ± 0.4 × 10−10 7.00 ± 0.6 × 10−9 8.88 ± 0.7 × 10−8 4.89 ± 0.4 × 10−7 1.21 ± 0.1 × 10−6 [Tf]−

3.08 ± 0.4 × 10−10 5.43 ± 0.5 × 10−9 7.28 ± 0.6 × 10−8 3.85 ± 0.3 × 10−7 7.34 ± 0.6 × 10−7

H2 O

7.27 ± 0.8 × 10−9 2.56 ± 0.3 × 10−7 2.02 ± 0.2 × 10−6 5.51 ± 0.4 × 10−6 9.22 ± 0.7 × 10−6

κ

6.7 ± 0.3 × 10−7 1.2 ± 0.15 × 10−5 9.4 ± 0.7 × 10−5

7.3 ± 0.1 × 10−4

14.8 ± 1.1 × 10−4

to water-water hydrogen bonding (see Fig. S10 of SI). Because of the low density of water, however, the peak corresponds, on average, to just two hydrogen bonded water dimers, and the distribution of water − across the sample is homogeneous. The OW-nitrogen (in -NH+ 3 ) and OW-oxygen (on [Tf] ) radial

distributions also clearly show the signature of hydrogen bonding. In agreement with the experimental results of Ref. 11, the hydrogen bonding of water is particularly strong and geometrically well defined in the case of OW-[Tf]− . As expected, the association of water and [Tf]− is progressively enhanced with decreasing T . However, this association is never as strong as in [DAEt][Tf],15 in which [Tf]− ions practically sequester water at T < 300 K. Fig. S10 shows a selection of OW-OW, OW-cation and OW-anions radial distribution functions, documenting these structural features. Moreover, typical water-[Tf]− and [DAEt]+ -water configurations are illustrated in Fig. S11 of SI. The effect of water on general chemical physics properties of [DAEt][Tf] turns out to be similar to those found in samples of [TEA][Tf] protic ionic liquid also at 1 % water contamination.15 Both [DAEt]+ and [Tf]− diffuse faster in wet samples than in dry samples. Water diffuses faster than either [DAEt]+ and [Tf]− at T ≥ 300 K. Below this temperature the water advantage in mobility is

32


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


progressively reduced, because of its tight association with [Tf]− . Quantitative results can be found in Tab. 4. As in the [TEA][Tf] case,15 electrical conductivity increases with respect to the dry samples, but the relative variation is fairly small. At variance from the [TEA][Tf] case, the Nernst-Einstein conductivity remains above the direct estimate (Eq. 6) at all temperatures, showing that even in the presence of water, pairing of cations and anions remains tight. More interesting, in this case, are the consequences on the statistics of chains of HB that represent paths for the fast proton diffusion. First of all, we observe that the inclusion of water enhances the connectivity of the hydrogen-bonding network, even at temperatures such that water is tightly bound to [Tf]− . At all temperatures, in particular, every water molecules that donates more than one proton links to different [Tf]− anions. Conversely, the rare times a water molecule accepts two hydrogen bonds, these are donated by different cations. The identification of chains now including also the bridging by water molecules (see Fig. 5 (c)) shows a widening of the range of sizes present at significant concentration at all simulates temperatures. Moreover, again for n > 40, the density of chains c(n) is higher in the water-contaminated samples, as shown in Fig. S13 of SI. Both observations point to a potential increase of proton conductivity upon water absorption. A similar effect in other hydrogen bonded systems is well documented.55 Unfortunately, efficiency requirements on electrochemical devices such as fuel cells force them to operate at temperatures of the order of 400 K, i.e., above the boiling point of water.56 Nevertheless, computations could provide the way to extrapolate experimental diffusion and conductivity measurements for water contaminated specimens to the limit of ideally dry samples, relevant for power generation applications.

IV.

CONCLUSIONS

Half-neutralised diamine cations combined with the triflate anion have been investigated by density functional computations in the generalised gradient approximation (PBE) and by molecular dynamics based on an empirical force field. Three cations have been considered, consisting of a pair of amineammonium terminations, separated, or, better, joined by a short alkane segment -(CH2 )n - with n = 2, 3 and 4 carbon atoms, whose gas-phase, ground state structure is shown in Fig. 1. Density functional computations for the isolated cations, for neutral ion pairs and for small neutral clusters with up to eight ions identify low-energy geometrical motifs in these small aggregates, highlighting the primary role of Coulomb forces, as well as the secondary but still important role of hydrogen bonding. In all the aggregates we studied, each cation shows a very stable intra-ion hydrogen bond joining its amine and ammonium terminations. Despite strain in the smallest member of the sequence, this intra-cation

33



Page 34 of 41

HB is stable, short, and probably more symmetric than implied by the structural formula, especially taking quantum effects into account. Strong and stable HBs join each cation to multiple anions (when present), contributing to the determination of the ground state structure. Of particular interest, however, is a different class of HBs, connecting neighbouring cations. Because of electrostatics, these bonds are certainly less stable than cation-anion HBs, and, for this reason, they are not found in the ground state structure of very small aggregates, but they form during

∼ 10 ps long ab-initio simulations at room temperature.

Moreover, they could form at sizeable concentration in extended connected phases, stabilised by entropy and by correlation with anions, and in such a case they would support proton transport via the Grotthuss mechanism. To investigate proton displacement by this process, a chain-like geometry made of hydrogen bonded [DAEt]+ cations and neutralising [Tf]− anions has been prepared. The potential energy profile upon moving a single proton along the HB has a single minimum, although it is soft enough to allow large amplitude displacements at a modest energy, requiring ∼ 5 kcal/mol to bring the proton to 1.1 ˚ A from the non-covalently bonded nitrogen. The corresponding energy profile for the in-phase displacement of all protons along the HB chain has a double minimum shape. The barrier per proton separating the two minima is low, but, of course, it increases linearly with the number n of protons moving together, hence the barrier becomes high when moving many protons. Although we did not analyse this aspect in detail, it is likely that a semi-localised soliton-like mode combining single proton and collective displacement might provide a mechanism with a barrier compatible with the role of extra-vehicular proton transport measured in experiments on the same ionic liquid compound.3 Molecular dynamics simulations based on an atomistic empirical force field have been carried out to determine general chemical physics properties of [DAEt][Tf] and [DABu][Tf], and to explore possible paths for proton conduction through the Grotthuss mechanism. Simulations at ambient pressure for [DAEt][Tf] have been carried out for temperatures from 200 K to 440 K on a thin grid of δT = 10 K spacing. Long trajectories covering 200 ns production following at least 6 ns equilibration have been analysed to compute thermodynamic properties such as density and average potential energy, structural properties such as the structure factors, the radial distribution functions, the number and geometry of the most stable HBs, as well as dynamical properties such as the diffusion constant, the electrical conductivity, the relaxation time for ion rotations, the spatial distribution and the life-time of individual hydrogen bonds. In the [DAEt][Tf] case, both density and potential energy as a function of T show two linear ranges at high and low temperature, with a crossover interval centred at T = 310 K of ∼ 20 K half-width, identifying what appears to be a glass transition with super-imposed a weakly first order transition, 34


Page 35 of 41 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


possibly giving origin to a very incomplete degree of ordering that we have been unable to identify. The structure of both high- and low-temperature phases displays periodicity on a few different lengths, hinting to the separation of ionic and neutral domains on the 6 − 10 ˚ A scale. The temperature evolution of dynamical properties reflects the two phases and broad transition seen in the thermodynamic functions, as emphasised by the Arrhenius plot of diffusion coefficients, electrical conductivity and rotational relaxation time. The sizeable overestimation of ionic conductivity by the Stokes-Einstein relation, that persists up to the highest temperature in our simulation, points to the long-time association of cations and anions, consistent with the strength of Coulomb interactions, the directionality of hydrogen bonding, and the resulting high viscosity of these compounds. MD simulations of [DABu][Tf] samples at four temperatures (T = 440 K, 400 K, 350 K and 300 K) and ambient pressure give results similar to those obtained for [DAEt][Tf], with expected differences due to the relatively lower role of Coulomb forces and to the increased role of dispersion and packing effects with respect to [DAEt][Tf]. In particular, the tendency to a mesophase is enhanced in [DABu][Tf]. The intra-cation HB is somewhat destabilised by the entropy advantage of extending the -(CH2 )4 - segment and separating the amine-ammonium terminations, resulting in a dynamical equilibrium between the ring and chain configuration for the cations. Remarkably, ions diffuse (slightly) faster in [DABu][Tf], an its electrical conductivity is higher than for [DAEt][Tf]. The enhanced fluidity of [DABu][Tf] is not unexpected, since it is in line with the role of larger and less symmetric ions in decreasing the melting temperature of RTILs down to room temperature. Possibly because of the lower number density of cations, sequences of [DABu]+ ions separated by less than rc = 3.6 ˚ A are shorter (by size n) in [DABu][Tf] than in [DAEt][Tf], both including and excluding the possibility of bridging by oxygen atoms. The exploration of possible paths for proton conduction arguably is the most intriguing part of our study, since the double donor-acceptor capability of cations opens the way to proton transport via the Grotthuss mechanism, breaking the limit on vehicular conductivity set by the relatively high viscosity of these compounds. Unfortunately, size and time scales prevent the usage of ab-initio simulation at this stage, and our analysis has to rely on the equilibrium structure provided by molecular dynamics and the empirical force field. Starting from each cation in turn, we identify self-avoiding sequences moving from one cation to the next according to a condition of closest proximity, terminating whenever the distance to the nearest neighbour is longer than rc = 3.6 ˚ A . If we consider only cations, and we define their separation as the minimum distance of amine and ammonium terminations (no distinction between the two) on the two ions, only short chains of up to ∼ 10 cations can be identified. The abinitio investigation of the potential energy landscape for protons moving between a nitrogen atom 35



Page 36 of 41

on [DAEt][Tf] and an oxygen on [Tf]− , however, shows that long elongation proton displacements are possible within an energy threshold of a few kcal/mol. We take this information into account by including the possibility of oxygen atoms from [Tf]− shuttling protons across the gap between two cations, following the procedure described in Sec. III D. In this way, we find connected sequences of up to 80 cations, consisting of a few short linear segments and a majority of globular domains. Although our approach is silent on the dynamical aspects underlying these observation, we think that the results based on structure might provide a first preliminary explanation of the observed high rate of proton transport in [DAEt][Tf],3 exceeding the rate implied by the diffusion constant of cations. In a last stage of our study, we investigated the effect of water contamination at the 1 % in weight on the properties of [DAEt][Tf]. Once again, the results follow expectation but are nevertheless intriguing. At T ≥ 300 K water diffuses faster in [DAEt][Tf] that either ions, which, themselves, diffuse faster than in dry samples. At all temperatures, the water-anion interaction is stronger than water-cation. Below T = 300 K, this tight coupling to [Tf]− slows down water, reducing significantly its mobility advantage with respect to both ions. At all temperatures water enhances the electrical conductivity, both because it enhances diffusion and because it decreases (slightly) the strict correlation in the motion of cations and anions. By donating two HBs to two different [Tf]− anions, and accepting, on average, one HB from [DAEt]+ , water enhances the connectivity of the hydrogen bond network, and increases the length of cation chains, as defined in the previous section. In this way, water contamination tends to increase the mobility of protons both by enhancing cation diffusion and, possibly, by allowing longer paths to Grotthuss transfer. Although the cation species in our RTILs include only short -(CH2 )n - segments, our results on the effect of water on the HB network and on the relative geometry of water and anions is reminiscent of those presented in Ref. 57 for longer alkyl chain RTILs. In summary, our extensive computational study provides a comprehensive picture of halfneutralised diamine and triflate salts, outlining a number of quantitative and qualitative details. In many respects, it extends the recent computational analysis we carried out on [TEA][Tf], and is part of a combined experimental and computational investigation of protic ionic liquids focusing,in particular, on their proton conduction capability. The major limitation of our approach is represented by the rigid atomic charges, fixed-topology, classical mechanics picture underlying the force field model used for MD simulations. Full blown DFT simulations for a few hundred RTIL ions, covering ∼ 100 ns are still years away, but qualitative improvements of the force fields, and perhaps the introduction of valence bond (VB) models could greatly extend the scope of simulation in the investigation of proton conductivity in organic ionic systems. 36


Page 37 of 41 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


Acknowledgment - A.B. acknowledge support from (i) Science Foundation Ireland (grant n. 15-SIRG-3538), and (ii) The Italian Ministry of Education, University and Research (grant n. MIURDM080518-372). Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: xxxxxxx. It reports: force field parameters; ground state structure of the ([DAEt][Tf])n dimer and tetramer; miscellaneous figures.

References 1

Silvester, D. S.; Compton, R. G. Electrochemistry in room temperature ionic liquids: A review and some possible applications. Z. Phys. Chem. 2006, 220, 1247-1274. doi: 10.1524/zpch.2006.220.10.1247

2

Greaves, T. L.; Drummond, C. J. Protic ionic liquids: Properties and applications. Chem. Rev. 2008, 108, 206-237. doi: 10.1021/cr068040u

3

Iojoiu, C.; Hana, M.; Molmeret, Y.; Martinez, M.; Cointeaux, L.; El Kissi, N.; Teles, J.; Leprêtre, J.-C.; Judeinstein, P.; Sanchez, J.-Y. Ionic liquids and their hosting by polymers for HT-PEMFC membranes. Fuel Cells 2010, 10, 778-789. doi: 10.1002/fuce.201000026

4

Hayes, R.; Warr, G. G.; Atkin, R. Structure and nanostructure in ionic liquids. Chem. Rev. 2015, 115, 6357-6426. doi: 10.1021/cr500411q

5

Greer, S. C. Physical chemistry of equilibrium polymerization. J. Phys. Chem. B 1998, 102, 5413-5422. doi: 10.1021/jp981592z

6

Kumar, M.; Venkatnathan, A. Mechanism of proton transport in ionic-liquid-doped perfluorosulfonic acid membranes. J. Phys. Chem. B, 2013, 117, 14449-14456. doi: 10.1021/jp408352w

7

Kumar, M.; Venkatnathan, A. Quantum chemistry study of proton transport in imidazole chains. J. Phys. Chem. B 2015, 119, 3213-3222. doi: 10.1021/jp508994c

8

Chen, H.-N.; Yan, T.-Y.; Voth, G. A. A computer simulation model for proton transport in liquid imidazole. J. Phys. Chem. A 2009, 113, 4507-4517. doi: 10.1021/jp811156r

9

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. doi: 10.1021/acs.jpcb.6b12418

10

Seddon, K. R.; Stark, A.; Torres, M.-J. Influence of chloride, water, and organic solvents on the physical properties of ionic liquids. Pure Appl. Chem. 2000, 72, 2275-2287. doi: 10.1351/pac200072122275

11

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

37



12

Page 38 of 41

Koishi, T, Molecular dynamics study of the effect of water on hydrophilic and hydrophobic ionic liquids. J. Phys. Chem. B 2018, 122, 12342-12350. doi: 10.1021/acs.jpcb.8b07774

13

Kelkar, M. S.; Maginn, E. J. Effect of temperature and water content on the shear viscosity of the ionic liquid 1-Ethyl-3-methylimidazolium Bis(trifluoromethanesulfonyl)imide as studied by atomistic simulations. J. Phys. Chem. B 2007, 111, 4867-4876. doi: 10.1021/jp0686893

14

Kreuer, K.-D. Proton Conductivity: Materials and Applications. Chem. Mater. 1996, 8, 610-641. doi: 10.1021/cm950192a

15

Mora Cardozo, J. F.; Burankova, T.; Embs, J. P.; Benedetto, A.; Ballone, P. Density functional computations and molecular dynamics simulations of the triethylammonium triflate protic ionic liquids. J. Phys. Chem. B, 2017, 121, 11410-11423. doi: 10.1021/acs.jpcb.7b10373

16

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

17

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. doi: 10.1021/ja9621760

18

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, and organic-molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. doi: 10.1021/ja00124a002

19

Rizzo, R.C.; Jorgensen, W. L. OPLS all-atom model for amines: resolution of the amine hydration problem. J. Am. Chem. Soc. 1999, 121, 4827-4836. doi: 10.1021/ja984106u

20

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. doi: 10.1021/jp0476545

21

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

22

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

23

Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comp. Chem. 2005, 1701-1718. doi: 10.1002/jcc.20291

24

Eastwood, J. W.; Hockney, R. W.; Lawrence, D. N. P3M3DP - The 3-dimensional periodic particleparticle-particle-mesh program. Computer Physics Communications 1980, 19, 215-261. doi: 10.1016/00104655(80)90052-1

25

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. doi: 10.1063/1.445869

26

D. Marx and J. Hutter, in Modern Methods and Algorithms of Quantum Chemistry, Ed.: J. Grotendorst, John von Neumann Institute for Computing, J¨ ulich, NIC Series, Vol. 3, 329-477 (2000).

27

CPMD, http://www.cpmd.org/, Copyright IBM Corp 1990-2015, Copyright MPI f¨ ur Festkörperforschung

38


Page 39 of 41 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


Stuttgart 1997-2001. 28

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865-3868. doi: 10.1103/PhysRevLett.77.3865

29

Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane wave computations Phys. Rev. B 1991, 43, 1993-2006. doi: 10.1103/PhysRevB.43.1993

30

Hockney, R. W. The potential calculation and some applications. Methods Comput. Phys. 1970, 9, 136-211.

31

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

32

Andersson, S.; Gr¨ uning, M. Performance of density functionals for calculating barrier heights of chemical reactions relevant to astrophysics. J. Phys. Chem. A 2004, 108, 7621-7636. doi: 10.1021/jp040448c

33

Becke, A. D. Density-functional thermochemistry.III. The role of exact exchange. J. Chem. Phys. 1993, 5648-5652. doi: 10.1063/1.464913

34

Lynch, B.J.; Truhlar, D. G. How well can hybrid density functional methods predict transition state geometries and barrier heights? J. Phys. Chem. A 2001, 105, 2936-2941. doi: 10.1021/jp004262z

35

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

36

Takasu, I.; Sugawara, T.; Mochida, T. Dielectric response in bisquaric acid crystal: Possible generation of protonic soliton in a quasi-one-dimensional hydrogen-bonded system. J. Phys. Chem. B 2004, 108, 1849518499. doi: 10.1021/jp030540t

37

Del Pópolo, M. G.; Pinilla, C.; Ballone, P. Local and semilocal density functional computations for crystals of 1-alkyl-3-methyl-imidazolium salts. J.Chem. Phys. 2007, 126, 144705. doi: 10.1063/1.2715571

38

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. doi: 10.1021/jp067589u

39

Hansen, J.-P. and McDonald, I. R. Theory of Simple Liquids, 3rd ed.; Academic Press: Amsterdam, 2006.

40

Urahata, S. M.; Ribeiro, M. C. C. J. Chem. Phys. Structure of ionic liquids of 1-alkyl-3-methylimidazolium cations: A systematic computer simulation study. 2004, 120, 1855-1863. doi: 10.1063/1.1635356

41

Wang, Y.; Voth, G. A. Unique spatial heterogeneity in ionic liquids. J. Am. Chem. Soc. 2005, 127, 1219212193. doi: 10.1021/ja053796g

42

Canongia Lopes J. N. A.; Padua, A. A. H. Nanostructural organization in ionic liquids. J. Phys. Chem. B 2006 110, 3330-3335. doi: 10.1021/jp056006y

43

Wang, Y.; Voth, G.A. Tail Aggregation and Domain Diffusion in Ionic Liquids. J. Phys. Chem. B 2006, 110, 18601-18608. doi: 10.1021/jp063199w

44

Triolo, A.; Russina, O.; Bleif, H.-J.; Di Cola, E. Nanoscale segregation in room temperature ionic liquids. J. Phys. Chem. B 2007, 111, 4641-4644. doi: 10.1021/jp067705t

45

Araque, J. C.; Hettige, J. J.; Margulis, C. J. Modern room temperature ionic liquids, a simple guide to understanding their structure and how it may relate to dynamics. J. Phys. Chem. B 2015, 219, 1272712740. doi: 10.1021/acs.jpcb.5b05506

39



46

Page 40 of 41

Le Donne, A.; Adenusi, H.; Porcelli, F.; Bodo, E. Hydrogen bonding as a clustering agent in protic ionic liquids: like-charge vs opposite-charge dimer formation. ACS Omega 2018, 3, 10589-10600. doi: 10.1021/acsomega.8b01615

47

Ferreira, M. L.; Pastoriza-Gallego, M. J.; Ara´ ujo, J. M. M.; Lopes, J. N. C.; Rebelo, L. P. N.; Pi˜eiro, M. M.; Shimizu, K.; Pereiro, A. B. Influence of nanosegregation on the phase behaviour of fluorinated ionic liquids. J. Phys.Chem. C 2017, 121, 5415-5427. doi: 10.1021/acs.jpcc.7b000516

48

Ballone, P.; Pastore, G.; Tosi, M. Structure and thermodynamic properties of molten alkali chlorides. J. Chem. Phys. 1984, 81, 3174-3180. doi: 10.1063/1.448022

49

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. doi: 10.1021/jp500296x

50

Köddermann, T.; Ludwig, R.; Paschek, D. On the validity of Stokes-Einstein and Stokes-EinsteinDebye relations in ionic liquids and ionic-liquid mixtures. ChemPhysChem 2008, 9, 1851-1858. doi: 10.1002/cphc.200800102

51

Andreussi, O.; Marzari, N. Transport properties of room-temperature ionic liquids from classical molecular dynamics. J. Chem. Phys. 2012, 137, 044508. doi: 10.1063/1.4737388

52

Voloshin, V. P.; Naberukhin, Yu. I. Hydrogen bond lifetime distribution in computer simulated water. J. of Structural Chem. 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. doi: 10.1103/PhysRevLett.76.928

54

Zahn, S.; Thar,J.; Kirchner, B. Structure and dynamics of the protic ionic liquid monomethylammonium nitrate [CH3 NH3 ][NO3 ] from ab initio molecular dynamics simulations. J.Chem.Phys. 2010, 132, 124506. doi: 10.1063/1.3354108

55

Loerting, T.; Liedl, K. R. Water-mediated proton transfer: A mechanistic investigation on the example of the hydration of sulfur oxides. J. Phys. Chem. A 2001, 105, 5137-5145. doi: 10.1021/jp0038862

56

Alberti, G; Casciola, M. Composite membranes for medium-temperature PEM fuel cells. Annual Review of Materials Research 2003, 33, 129-154. doi: 10.1146/annurev.matsci.33.022702.154702

57

Bodo, E.; Mangialardo, S.; Capitani, F.; Gontrani, L.; Leonelli, F.; Postorino, P. Interaction of a long alkyl chain protic ionic liquid and water. J. Chem. Phys. 2014, 140, 2004503. doi: 10.1063/1.4876036

40


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


FIG. 17: TOC graphics

41


Equilibrium Structure, Hydrogen Bonding, and Proton Conductivity in

Recommend Documents