Article pubs.acs.org/JPCB
How Is Diffusion of Neutral and Charged Tracers Related to the Structure and Dynamics of a Room-Temperature Ionic Liquid? Large Deviations from Stokes−Einstein Behavior Explained Juan C. Araque,† Sharad K. Yadav,† Michael Shadeck,‡ Mark Maroncelli,‡ and Claudio J. Margulis*,† †
Department of Chemistry, University of Iowa, Iowa City, Iowa 52242, United States Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802, United States
‡
S Supporting Information *
ABSTRACT: The deviations from Stokes−Einstein hydrodynamics of small solutes are more pronounced in ionic liquids than in conventional solvents (J. Phys. Chem. B 2013 117 (39), 11697). Small neutral solutes diffuse much faster than expected, whereas small charged solutes diffuse much slower. This article attempts to establish a link between the local friction experienced by tracer solutes and the polar/apolar structure of ionic liquids. We find that small neutral solutes probe locally “stiff” (mostly charged, high electrostriction) regions and locally “soft” (mostly apolar, low electrostriction) regions. These regions of high and low friction are associated with cage and jump regimes. Enhanced neutral tracer mobility in the low friction regions associated with the cationic apolar component has an important bearing on the large positive deviations from Stokes−Einstein behavior. In contrast, diminished charged tracer mobility involves long caging dynamics separated by jump events often triggered by the loss and recovery of counterions.
1. INTRODUCTION A deeper mechanistic understanding of the diffusion process in the case of charged and neutral tracers in modern ionic liquids (ILs) is desirable because it is at the core of numerous technological areas such as charge transport for energy conversion or storage devices,1−3 diffusion of reactants and intermediates in homogeneously catalyzed reactions,4−7 or selective recovery of solutes in liquid−liquid extraction processes.8,9 When the size ratio between the solute and solvent is small, the detailed liquid landscape becomes relevant and deviations from Stokes−Einstein (SE) behavior apparent.10 Experimental evidence is available in common solvents,11 polymers,12,13 and supercooled liquids.14 This phenomenon becomes extreme in the case of ILs15−17 because the liquid landscape is often heterogeneous and friction, normally a property of the solvent, becomes dependent on the charge of the solute and its location. One example of such anomalous behavior occurs in bimolecular electron transfer reactions. In ILs, it has been experimentally found that if reactant species are small and neutral, rates are anomalously high,5,18−24 whereas if they are charged, reactions are anomalously slow.25−28 Here anomalous refers to a comparison with what would be expected from theories that assume friction is independent of the solute charge and the same everywhere in the solvent. The current work attempts to provide mechanistic insight into intriguing phenomena depicted in Figure 1. In Figure 1b, the blue and red lines are fits of the ratio between the observed and SE predicted friction (defined in eq 1) to an extensive set © XXXX American Chemical Society
of experimental tracer diffusion data, derived from a wide variety of experimental data collected in refs 15 and 16. Open circles correspond to values derived in the current work from molecular dynamics (MD) simulations and filled circles from NMR measurements of different solutes, all in 1-butyl-1methylpyrrolidinium bis(trifluoromethanesulfonyl)imide ([Pyrr4,1+][NTf2−]). SE behavior occurs for both neutral and charged solutes in the limit of large solute to solvent volume ratio Vu/Vv , but when this ratio is smaller than one, significant deviations occur. kBT /6πηobsR solute ζobs D = SE = ζSE Dobs Dobs
(1)
We will show that in the case of small neutral solutes diffusion is likely not dominated by solute−solvent interactions but instead by solvent−solvent interactions. Because of the polar/apolar nature of the solvent, small neutral tracers sample regions that are stiff, where solvent−solvent friction is high and regions that are soft where it is low. Neutral tracer mobility is therefore enhanced by diffusion through the soft regions of the liquid. Stiff and soft solvent locations naturally correlate with the polar or apolar components of the medium. Soft, mostly apolar solvation plays a defining role on the faster than expected diffusion of neutral tracers. Received: February 2, 2015 Revised: March 24, 2015
A
DOI: 10.1021/acs.jpcb.5b01093 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
different way, the tracer is at all times an intrinsic part of the IL charge network. On average, a small charged solute will experience significantly higher friction than would be expected by the macroscopic viscosity; therefore, diffusion is retarded.
2. METHODOLOGY In the present study we consider dilute solutions of spherical (single atom/ion) and quasi-spherical (tetrahedral symmetry) solutes. The rationale for this choice is to avoid undesirable rotational−translational coupling as well as to diminish to some extent the effect of dipolar and other multipolar interactions. As the solvent, we have selected 1-butyl-1-methylpyrrolidinium bis(trifluoromethanesulfonyl)imide [Pyrr4,1+][NTf2−] because substantial experimental and computational structural data are available31,32 which shows there is no first sharp diffraction peak in this liquid. The absence of a first sharp diffraction peak or prepeak is indicative of an absence of preformed apolar regions separated by charged regions.33,34 Work on the mechanisms of diffusion in systems where there is a prepeak indicative of polarity alternation will be the subject of a future publication. As we will show later, the absence of a polarity alternation prepeak does not preclude preferential diffusion in regions that are mostly polar or mostly apolar. The molecular structure of all tracers as well as the solvent is schematically depicted in Figure 1a. As can be seen from parts a and b of Figure 1, the chosen combination of solutes and solvent results in sampling of Vu/Vv in the full relevant range between 0.08 to 1.0. In our simulations both the solvent and solutes are represented at the all-atom level. A list of the sources of force field parameters is given in Table S1 of the Supporting Information (SI). Model parameters for bonded and nonbonded interactions were chosen according to the OPLS-AA35 prescription. As can be seen from Table S1 in the SI, the [Pyrr4,1+] potential is derived from the Canongia Lopes and Padua (CLP)36,37 parameters, whereas the [NTf2−] parameters are those of Köddermann, Paschek, and Ludwig (KPL).38 The latter parametrization has been shown to provide a much better description of transport properties without resorting to the typical method of charge scaling.39 Such improvement was achieved by fine-tuning both Lennard-Jones parameters (εij and σij) while leaving the CLP partial atomic charges unmodified.38 Data presented in Figure S1 in the SI show that this combined CLP(cation)−KPL(anion) formulation gives a 2-fold deviation of the calculated diffusion coefficients (D+s and D−s ) from experimental values at 303 K,40 whereas a 20-fold deviation is obtained using CLP with charges scaled by 80% and a 40-fold deviation using the original CLP with full charges. We show in Figures S2 and S3 in the SI that these parameters acceptably reproduce the liquid structure factor and the viscosity when compared with experiment. The study of dynamical properties of tracer molecules has always been problematic because simulations must be performed at low enough solute concentrations such that solvent properties are minimally altered but statistics are sufficient so that average results are meaningful. In the case of our studies the tracer concentration was set at ∼100 mM, which is only twice the experimental value used in the NMR refs 24 and 15. In all cases, 10 tracer solutes were simulated in a box of 343 solvent ion pairs. When simulating the positively charged solutes (seen in Figure 1a), an equal number of “extra” [NTf2−] anions were included to maintain charge neutrality.
Figure 1. (a) Schematic representation of: (top) molecular structure of the neat IL constituents 1-butyl-1-methylpyrrolidinium cation [Pyrr4,1+] and bis(trifluoromethanesulfonyl)imide anion [NTf2−] and (bottom) molecular structure of the solutes considered ordered according to their molecular volumes and classified as neutral or charged. The color coding of the atoms is as follows: carbon, cyan; hydrogen, white; nitrogen, blue; oxygen, red; sulfur, yellow; fluorine, green. (b) Translational friction ratio, ζobs/ζSE, for spherical and quasispherical solutes with either neutral or positive (+1) charge. Comparison is between MD simulation estimates (at 400 K) and experimental data from current work (at 298 K), except data for Li+ in [Pyrr4,1+][NTf2−] at 303 K, which is from ref 29. Blue and red lines are fits to diffusion NMR data in ILs from refs 15 and 16: (ζobs/ζSE)neutral = {1 + 1.93(Vu/Vv)−1.88}−1 for dilute neutral solutes and (ζobs/ζSE)±1 = (ζobs/ζSE)neutral+0.198(Vu/Vv)−0.679 for dilute charged (±1) species, respectively. All volumes are computed using the method described in ref 30 and are listed in Tables S3 and S4 in the Supporting Information. For clarity, not all experimental tracer data points are labeled; actual experimental and computational friction values are provided in the same Tables.
In contrast, the large positive deviations in the case of small charged solutes are correlated with strong caging that results in an enlarged hydrodynamic radius of the tracer. Put in a B
DOI: 10.1021/acs.jpcb.5b01093 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B [Pyrr4,1+][NTf2−] has an experimental viscosity of 74 cP at 298 K, whereas it is ∼4.5 cP at 400 K.24 Because higher ionic mobilities provide better statistics for all of our calculations, we decided to carry out all of our production runs at 400 K. Simulated viscosities for the neat liquid and solute systems are provided in Table S2 of the SI. Two types of production runs are needed for the mechanistic studies described here. Constant temperature and volume (NVT) calculations are required to properly determine the liquid viscosities needed in the accurate calculation of the SE friction (eq 1). Constant energy and volume (NVE) calculations are instead required to estimate diffusion constants that are not spuriously contaminated by the effect of thermostats and barostats. All of our simulations were performed using the GROMACS MD package,41,42 version 4.5.5, at the double-precision level. Equations of motions were integrated using a leapfrog algorithm43 with a time step of 1 fs during equilibration and with a 2 fs step during long production runs. The size of our cubic simulations boxes ranged between 5.7 and 5.9 nm depending on the system, and in all cases periodic boundary conditions were applied in all three dimensions. Nonbonded pair interactions were considered using the minimum image convention with van der Waals and Coulomb (real-space sum) cutoff radii of 1.50 nm during equilibration and switch/cutoff radii of 1.18/1.20 nm during production. Longer ranged electrostatic forces and energies beyond the specified cutoff were calculated with the particle mesh Ewald (PME) method44,45 with B-spline interpolation of order six and fourierspacing parameter set to 0.08(equilibration)/0.12(production) nm. The accuracy of the reciprocal space sum, both during equilibration and production, was kept by assigning a tolerance level of 1 × 10−5 to the realspace electrostatic interaction at the cutoff distance. Additionally, all covalent bonds involving hydrogen atoms were constrained to their equilibrium values using the linear constraint solver (LINCS) algorithm.46,47 A careful equilibration protocol was followed to ensure proper sampling of transport properties during production runs. We began our simulations from a lattice. Solute molecules were then interdispersed within the IL crystal to achieve the prescribed concentration. This initial setup was energyminimized using the steepest descent algorithm. The system was then simulated at 0, 10, and 100% charges in the NpT ensemble. During this last step a temperature-annealing protocol was implemented that raised the temperature to 500 K and lowered it back to 400 K, after which it was held constant for another 2 ns. For each tracer system, the final equilibrated positions and velocities were used as initial conditions for a production NVE run of 30 ns. Energy conservation was carefully verified. Typical energy drifts were found to be
