Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations

Subscriber access provided by University of Florida | Smathers Libraries

Article

Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations Brian Doherty, Xiang Zhong, Symon Gathiaka, Bin Li, and Orlando Acevedo J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00520 • Publication Date (Web): 07 Nov 2017 Downloaded from http://pubs.acs.org on November 9, 2017

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

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

Page 1 of 46

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

Journal of Chemical Theory and Computation

Revised version 10/12/17

Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations Brian Doherty, Xiang Zhong, Symon Gathiaka, Bin Li, and Orlando Acevedo* Department of Chemistry, University of Miami, Coral Gables, Florida 33146 E-mail: [email protected] Submitted May 18, 2017 Abstract: Our OPLS-2009IL force field parameters (J. Chem. Theory. Comput. 2009, 5, 10381050) were originally developed and tested on 68 unique ionic liquids featuring the 1-alkyl-3methylimidazolium [RMIM], N-alkylpyridinium [RPyr], and choline cations. Experimental validation was limited to densities and a few, largely conflicting heat of vaporization (DHvap) values reported in the literature at the time. Owing to the use of Monte Carlo as our sampling technique, it was also not possible to investigate the reproduction of dynamics. The [RMIM] OPLS-2009IL parameters have been revisited in this work and adapted for use in molecular dynamics (MD) simulations. In addition, new OPLS-AA parameters have been developed for multiple

anions,

i.e.,

AlCl4–,

BF4–,

Br–,

Cl–,

NO3–,

PF6–,

acetate,

benzoate

bis(pentafluoroethylsulfonyl)amide, bis(trifluoroethylsulfonyl)amide, dicyanamide, formate, methylsulfate,

perchlorate,

propionate,

thiocyanate,

tricyanomethanide,

and

trifluoromethanesulfonate. The computed solvent densities, heats of vaporization, viscosities, diffusion coefficients, heat capacities, surface tensions, and other relevant solvent data compared favorably with experiment. A charge scaling of ± 0.8 e was also investigated as a means to mimic polarization and charge transfer effects. The 0.8-scaling led to significant improvements for DHvap, surface tension, and self-diffusivity; however, a concern when scaling charges is the potential degradation of local intermolecular interactions at short ranges. Radial distribution

ACS Paragon Plus Environment


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

functions (RDFs) were used to examine cation-anion interactions when employing 0.8*OPLS2009IL and the scaled force field accurately reproduced RDF’s from ab initio MD simulations.

Introduction Ionic liquids are a remarkable class of solvents composed exclusively of ions, characteristically at room temperature.1,2 Cations, such as 1-alkyl-3-methylimidazolium [RMIM], can be fine-tuned via different functional groups, e.g. R = M (methyl), E (ethyl), B (butyl), etc., in order to enhance the degree of localized structuring with the anions in the liquid phase.3-6 Consequently, different permutations of a low symmetry organic cation and a weakly coordinating inorganic or organic anion can have dramatic consequences upon solvent properties, including the melting point, viscosity, density, electrical conductance, solvent polarity, and gas solubility.7-10 In addition to their massive customizability (with an estimated 1018 different combinations!),11-13 ionic liquids exhibit many attractive properties, including low vapor pressures and excellent thermal and chemical stabilities, that allow them to be practical alternatives to traditional organic solvents.14-17 The physical and chemical properties of ionic liquids are largely dependent upon solvent structure,3,18 including the organization resulting from standard ion-ion interactions and specific hydrogen-bonding between cations and anions. While not wholly understood, experimental19-21 and theoretical22-25 efforts have provided significant insight into the nature of hydrogen bonds within imidazolium-based ionic liquids. The C-H groups of the imidazolium ring, and in particular the H2-C2 group, have been shown to play a predominantly large role in dictating solvent properties26-28 and upon chemical reactivity when used as the reaction medium, e.g. the Diels-Alder reaction.29-31

2 ACS Paragon Plus Environment

Page 2 of 46

Page 3 of 46

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


To better understand the relationship between intermolecular interactions and macroscopic observations, extensive effort has been put forth by the ionic liquid simulation community to develop force field potentials for atomistic simulations beginning with LyndenBell in 2001.32 Early imidazolium-based ionic liquid force fields by the groups of Maginn,33,34 Berne,35 and Stassen36,37 helped shape initial thoughts as to the origin of their unique solvent properties. Considerable theoretical work has followed over the past decade; reviews are available.38-43 However, shortcomings in existing force fields, e.g., poor dynamics,44 have led to critical evaluations45,46 and a constant call for more accurate ionic liquid potentials.27,47,48 Alternatives to fixed charge force fields such as ab initio molecular dynamics,49-52 polarizable force fields,53-55 and first-principles force fields56-58 have varied advantages over their nonpolarizable counterparts, but generally share the same disadvantage, namely an increased cost in computational resources. When our own OPLS-based ionic liquid force field (i.e., OPLS-2009IL) was published in 2009 and tested upon 68 unique ion combinations, experimental validation was limited to densities and the few existing and largely conflicting heats of vaporization literature values.59 As a consequence of employing Monte Carlo (MC) as our sampling technique, it was not possible to investigate the reproduction of dynamics for the room temperature molten salts. In this work, the OPLS-2009IL parameters have been reinvestigated using molecular dynamics (MD) by adapting our parameters to the GROMACS software package60 and validated against a significantly larger and more reliable set of experimental data, including densities, heats of vaporization, viscosities, diffusion coefficients, heat capacities, and surface tensions. In addition, whereas OPLS-2009IL parameters for a fully flexible cation were previously reported, the anions were modeled as rigid



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

systems that only considered their nonbonded interactions. In this work, new bond, angle, and torsion parameters have been developed for the anions shown in Figure 1. Polarizable ionic liquid models, e.g. applying self-consistent inducible dipoles53,61 or Drude oscillators,54 have been shown to give major improvements to the description of dynamics.44,55,62 Still, the positive gains are counter-balanced by significant requirements in computer power. An alternative, cost-effective approach to account for polarization effects is the scaling of atomic charges to mimic the average charge screening caused by the polarization as well as the charge transfer effects.63,64 Morrow and Maginn were the first to report a partial charge assignment for the [BMIM][PF6] ion pair based on a non-integer molecular charge of ± 0.904 e.34 Subsequent work determined a molecular charge scaling of ± 0.8 e to be generally ideal for reproducing experimental ionic liquid-phase properties including enthalpy of vaporization, surface tension, and self-diffusion coefficients.65-68 However, Schröder noted some discrepancies at the local level between the scaled charges and polarizable models, e.g. radial distribution functions (RDFs) below distances of 8 Å for [EMIM][TfO].69 This brings into question the usefulness of down-scaling solvent charges, particularly, in mixed quantum and molecular mechanics (QM/MM) chemical reaction studies.70-72 In this work, a charge scaling of ± 0.8 e has also been investigated for the OPLS-2009IL force field and computed RDFs for intermolecular cation-anion and cation-cation interactions are compared to published RDFs derived from ab initio molecular dynamics and polarizable force fields.


Page 4 of 46

Page 5 of 46

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


Figure 1. Ionic liquid-forming ions. R = E (ethyl), B (butyl), or O (octyl).

Computational Methods OPLS-AA Force Field. In the non-polarizable OPLS-AA force field formalism,73 the total energy for each ionic liquid system are evaluated as a sum of individual energies for the harmonic bond stretching and angle bending terms, a cosine series for torsional energetics, and Coulomb and 12-6 Lennard-Jones terms for the nonbonded interactions (Eq. 1-4). The parameters are the force constants k, the ro and qo equilibrium bond and angle values, Fourier coefficients V, partial atomic charges, q, and Lennard-Jones radii and well-depths, s and e.



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

Page 6 of 46

Ebonds = å kb,i (ri - ro,i )2

(1)

Eangles = å kb,i (qi - qo,i )2

(2)

i

i

Etorsion = ∑[ 12 V1,i (1+ cos φi ) + 12 V2,i (1− cos 2φi ) + 12 V3,i (1+ cos3φi ) + 12 V4,i (1− cos 4φi )]

(3)

i

Enonbond = åå{ i

j >i

qi q j e2 rij

+ 4e ij [(

s ij rij

)12 - (

s ij rij

)6 ]}

(4)

Geometric combining rules for the Lennard-Jones coefficients were utilized, i.e., sij = (siisjj)1/2 and eij = (eiiejj)1/2. Nonbonded interactions were evaluated intermolecularly and for intramolecular atom pairs separated by three or more bonds. To apply the same parameters for both intra- and intermolecular interactions the 1,4-intramolecular interactions were reduced by a factor of 2. To maintain compatibility with OPLS-AA, published potentials were retained without change whenever appropriate. This includes assigning all standard bond stretching and angle bending force constants from OPLS-AA, which may also include some entries from the AMBER force field.74,75 All Lennard-Jones parameters also came from the OPLS-AA parameter set except when explicitly stated. Ab Initio Calculations. Individual anions were optimized using the Jaguar program76 at the Hartree-Fock (HF) theory level using the 6-31G(d) basis set, followed by single-point energy calculations using the local MØller-Plesset second-order perturbation (LMP2) method77,78 and the correlation-consistent polarized valence cc-pVTZ(-f) basis set.79 This LMP2/cc-pVTZ(-f)//HF/631G(d) method is common practice for OPLS-AA parameterization.80 Vibrational analytical frequency calculations at the HF/6-31G(d) were carried out to confirm all minima as stationary points. The ab initio derived ion geometries were used for the equilibrium bond and angle, ro and

qo, reference values in the force field. Partial charges were computed by fitting the molecular 6 ACS Paragon Plus Environment

Page 7 of 46

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


electrostatic potential (ESP) at the atomic centers. For a better description of the charge density, LMP2 dipole moments were also computed along with a coupled perturbed Hartree-Fock (CPHF) term. Charges were symmetrized for similar atoms and used for the Coulombic nonbonded force field partial charges. Molecular Dynamics. Molecular dynamic (MD) simulations were performed using the GROMACS 5.0.7 software package.60 Cubic boxes containing 500 ion pairs were constructed using the program Packmol.81 Periodic boundary conditions were applied to each box with longrange interactions handled with Particle-Mesh Ewald summations. Boxes were minimized using a steepest descent algorithm for 5000 steps. Equations of motion were integrated using the leapfrog algorithm with a time step of 1 fs. The temperature value of 298 K was kept constant using velocity rescaling with a stochastic term (v-rescale)82 and a constant pressure of 1.0 bar was maintained with the Berendsen coupling during an isothermal-isobaric ensemble (NPT) simulation for 5 ns of equilibration. [EMIM][PF6] was simulated at 333 K, as that ionic liquid is solid at room temperature.83 All covalent hydrogen bonds were constrained using the LINCS algorithm and a cutoff range for the short-range electrostatics was set to 13 Å. Production runs were set for an additional 40 ns. Radial distribution functions (RDF), combined distribution functions (CDF), and spatial distribution functions (SDF) were computed using the TRAVIS program.84

Results and Discussion Parameters. Non-polarizable force fields are routinely parameterized using quantum mechanically derived charges fit to molecular electrostatic potential (ESP) at the atom centers of isolated gas-phase ions or individual cation-anion pairs.75,85 The OPLS-2009IL partial charges



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

were derived in this fashion by using the charge density from LMP2/cc-pVTZ(-f)//HF/6-31G(d) calculations and LMP2 dipole moments computed along with a coupled perturbed Hartree-Fock (CPHF) term for each gas-phase ion.59 For the cations, multiple low-energy geometry configurations exist stemming primarily from alkyl side-chain dihedral rotations. Developing a single set of charges for [RMIM] that could accurately represent different alkyl lengths and orientations was challenging. ESP charges were initially computed for all existing [RMIM] energy-minimized stationary points, where R = Me, Et, and Bu. An average partial charge value for each atom in [RMIM] was developed by appropriately weighting the contribution of each ground-state structure to the overall conformational population, e.g., using a two-state model with a Boltzmann distribution for [BMIM]’s gauche and trans side-chain conformers. The use of averaged point cation charges from individually parameterized [MMIM], [EMIM], and [BMIM] cations was necessary to construct the fully transferable [RMIM] OPLS-2009IL force field (Figure 2). Charges and 12-6 Lennard-Jones parameters for [RMIM] are given in Table 1 and include the 0.8-scaled charges. The resultant geometries from the ab initio calculations on [RMIM] were used to obtain the equilibrium bond and angle reference values, ro and θo, respectively. Force constants were taken directly from published OPLS-AA parameters73 whenever possible. For the torsion parameters, the development of new Fourier coefficients involved direct adjustment of OPLSAA torsional parameters to reproduce the energy differences between conformational energy minima from LMP2/cc-pVTZ(-f) calculations. All [RMIM] bonded terms are given in the Supporting Information Tables S1 and S2. The [AlCl4], [BF4], [NO3], [PF6], and [TFO] anions originally published in the 2009 paper only listed nonbonded parameters, i.e., charges and Lennard-Jones terms, as the anions were modeled as rigid structures. In this work, new


Page 8 of 46

Page 9 of 46

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


intramolecular bond, angle, and torsion parameters were developed in a similar fashion to [RMIM] for the anions listed in Figure 1. The new anion parameters are available in Supporting Information Tables S3-S5. OPLS-AA anion parameters were taken directly from the literature when proper validation was reported. For example, Canongia Lopes and Pádua took great care in producing OPLS-AA charges and dihedral parameters that accurately reproduced MP2/ccpVTZ(-f) torsion energy scans for [NTf2].86 Any OPLS-AA parameters obtained from the literature are referenced within the Supporting Information.

Table 1. Charges (unscaled and scaled by 0.8) and Lennard-Jones Parameters for the 1-Alkyl-3methylimidazolium [RMIM] Cations.a Atom type

a

q (e)

s (Å)

e (kcal/mol)

2009IL

0.8*2009IL

CR

-0.09

-0.072

3.55

0.070

NA

0.22

0.176

3.25

0.170

CW

-0.24

-0.192

3.55

0.070

CM

-0.35

-0.280

3.50

0.066

CA

-0.17

-0.136

3.50

0.066

CS

-0.12

-0.096

3.50

0.066

CT

-0.24

-0.192

3.50

0.066

HR

0.21

0.168

2.42

0.030

HW

0.27

0.216

2.42

0.030

HM

0.18

0.144

2.50

0.030

HA

0.18

0.144

2.50

0.030

HS

0.06

0.048

2.50

0.030

HT

0.08

0.064

2.50

0.030

Atom types for [RMIM] given in Figure 2.



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 46

Figure 2. OPLS atom types for the 1-alkyl-3-methylimidazolium [RMIM] cations.

Density. Computed solvent densities for the ionic liquids at 298 K are given in Table 2 for the unscaled and 0.8-scaled OPLS-2009IL force field. An extensive literature search on imidazolium-based ionic liquids uncovered a large amount of experimental bulk phase data available for comparison. As much of the data had reported uncertainties, a weighted mean was computed per solvent by determining the weight for each data point from the inverse of its uncertainty. All individual experimental values, uncertainties, measurement methodologies, and references are given in the Supporting Information. The calculated deviation from experimental densities was relatively small for the unscaled OPLS-2009IL force field with an overall mean absolute error (MAE) of 2.7%. Individual percent errors compared to the experimental weighted mean are also reported for each ion combination. Whereas in some individual cases scaling the charges improved agreement with experiment, e.g., [BMIM][NO3], in other cases a larger deviation was observed (Table 2). Scaling the force field systematically reduced the density for each solvent pair and subsequently raised the MAE to 2.9%.


Page 11 of 46

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


Table 2. Calculated and Experimental Liquid Densities (g/cm3) at 298 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. 2009IL

% error

0.8*2009IL

% error

Exptl.b

[BMIM][ACE]

1.064

1.3

1.039

3.6

1.078

[BMIM][AlCl4]

1.201

3.0

1.153

6.9

1.238

[EMIM][BF4]

1.273

0.4

1.212

5.2

1.278

[BMIM][BF4]

1.196

1.1

1.150

4.9

1.209

[OMIM][BF4]

1.104

0.3

1.074

3.0

1.107

[BMIM][BNZ]

1.110

-

1.093

-

-

[BMIM][Br]

1.317

1.6

1.277

1.5

1.296

[BMIM][Cl]

1.080

0.0

1.048

3.0

1.080

[BMIM][ClO4]

1.281

2.3

1.245

0.6

1.253

[EMIM][DCA]

1.113

0.8

1.069

3.1

1.103

[BMIM][DCA]

1.077

1.6

1.039

2.0

1.060

[BMIM][HCOO]

1.084

1.0

1.050

2.1

1.073

[EMIM][MS]

1.268

1.1

1.229

4.1

1.282

[BMIM][MS]

1.208

0.3

1.177

2.9

1.212

[BMIM][NO3]

1.189

2.8

1.155

0.1

1.156

[EMIM][NPF2]

1.683

5.7

1.643

3.1

1.593

[BMIM][NPF2]

1.592

5.2

1.566

3.5

1.513

[EMIM][NTf2]

1.630

7.4

1.584

4.4

1.517

[BMIM][NTf2]

1.548

7.8

1.504

4.7

1.436

[EMIM][PF6]a

1.430

-

1.363

-

-

[BMIM][PF6]

1.356

0.5

1.315

3.5

1.363

[OMIM][PF6]

1.238

1.1

1.207

3.5

1.251

[BMIM][PROP]

1.051

1.5

1.026

0.9

1.035

[EMIM][SCN]

1.204

7.8

1.130

1.2

1.117

[BMIM][SCN]

1.107

3.3

1.068

0.4

1.072

[EMIM][TCM]

1.104

2.1

1.055

2.4

1.081

[BMIM][TCM]

1.068

2.1

1.028

1.8

1.046

[EMIM][TfO]

1.466

5.9

1.419

2.6

1.384

[BMIM][TfO]

1.368

5.3

1.333

2.6

1.299

Ionic Liquid

MAE (%) a

2.7

2.9

b

Simulation at 333 K. Weighted experimental average; all individual values, uncertainties, methodology, and references given in the Supporting Information.



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 46

Heats of Vaporization. Ionic liquids possess vaporization enthalpies that are generally an order of magnitude greater than that of molecular liquids due to the strong electrostatic interactions between the ions.9 Reproducing heats of vaporization, DHvap, is particularly important as it serves as a key representative of the average intermolecular interactions. When the OPLS-2009IL parameters were originally published,59 very few and frequently conflicting experimental DHvap values were available for comparison e.g., 154.887 and 191.688 kJ/mol for [BMIM][PF6]. Fortunately, recent advances in experimental techniques have led to a substantial number of newly reported ionic liquid DHvap measurements.89 Several of the more popular ion pair combinations now have multiple experimental DHvap values (see the Supporting Information). The DHvap values for the ionic liquids investigated in this work were computed using equation 5. ∆𝐻#$% = ∆𝐻'$( − ∆𝐻*+,-+. = 𝐸010$* 𝑔𝑎𝑠 − 𝐸010$* 𝑙𝑖𝑞𝑢𝑖𝑑 + 𝑅𝑇

(5)

Experimental evidence has shown that ionic liquids go into the vapor phase as an ion pair,9 therefore Etotal (gas) was computed using a single ion pair and Etotal (liquid) is the corresponding value in the condensed phase. The RT term is used in place of a PV-work term in the enthalpy. To calculate Etotal (gas), 10 ns simulations were performed using the same setup described for liquid simulations in the methods section. The computed vaporization enthalpies gave an overall MAE of 33.8% when using the unscaled charges (Table 3). Many other nonpolarizable force fields gave similar outsized DHvap predictions.41,47 The general trend seen among ionic liquid force fields are accurate densities, but dynamical properties that are too slow and vaporization enthalpies that are too large. As previously mentioned, a molecular charge scaling of ± 0.8 e improves predictions of ionic liquid-phase properties including enthalpy of vaporization, surface tension, and ion diffusion coefficients.65-68 Accordingly, the OPLS-2009IL 12 ACS Paragon Plus Environment

Page 13 of 46

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


force field benefited immensely from a 0.8-charge scaling and yielded a noticeable improvement in DHvap with an overall MAE of 11.2% (Table 3). A plot of calculated versus experimental DHvap is provided in Figure 3 as a visual representation of the extensive data. In the case of [BMIM] coupled to carboxylate containing anions, i.e., [HCOO], [ACE], and [PROP], the computed errors are disproportionately large relative to the rest of the predictions. For example, the 0.8-scaled force field yielded percent errors of 56.4%, 24.3%, and 29.5% for [BMIM] with [HCOO], [ACE], and [PROP], respectively, which are significantly higher than the more typical low single to double digit percent errors computed for the other ionic liquids (Table 3). Direct experimental determination of the molar enthalpy of vaporization for these solvents proved to be very difficult. As such, Xu et al. indirectly determined90 the DHvap using Kabo and Verekin’s empirical equation based on surface tension.87 These carboxylatecontaining anions yielded DHvap values that are lower, i.e., 95.4, 123.8, and 121.8 kJ/mol for [BMIM] with [HCOO], [ACE], and [PROP], respectively, than the ~135-165 kJ/mol values generally reported for other ionic liquid combinations. Accordingly, they have a disproportionally negative effect on the overall MAEs for both the unscaled and scaled OPLS2009IL predicted DHvap values. For example, removing the [BMIM] with [HCOO], [ACE], and [PROP] combinations lower the overall unscaled and 0.8-scaled DHvap MAEs from 33.8% and 11.2%, respectively, to 28.6% and 7.7%. Other solvent properties for ionic liquids containing the carboxylate anions were accurately reproduced with the [RMIM] parameters. For example, viscosity calculations for [BMIM][HCOO] yielded values of 39.3 and 38.0 cP for the unscaled and 0.8-scaled OPLS-2009IL, respectively, compared to the experimental value of 38 cP.91



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 14 of 46

Table 3. Calculated and Experimental Heat of Vaporization (kJ/mol) at 298 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. 2009IL

% error

0.8*2009IL

% error

Exptl.b

[BMIM][ACE]

192.5

55.5

153.9

24.3

123.8

[BMIM][AlCl4]

163.1

5.2

125.2

27.2

172.0

[EMIM][BF4]

183.7

32.1

135.1

2.8

139.0

[BMIM][BF4]

186.3

36.6

140.5

3.0

136.4

[OMIM][BF4]

201.3

24.2

159.0

1.9

162.0

[BMIM][BNZ]

211.7

-

177.4

-

-

[BMIM][Br]

211.3

39.1

157.7

3.8

151.9

[BMIM][Cl]

213.0

41.0

157.7

4.5

151.0

[BMIM][ClO4]

200.8

-

157.7

-

-

[EMIM][DCA]

182.4

16.6

137.1

12.3

156.4

[BMIM][DCA]

180.7

14.8

140.2

11.0

157.5

[BMIM][HCOO]

190.6

99.8

149.2

56.4

95.4

[EMIM][MS]

200.5

34.0

158.1

5.7

149.6

[BMIM][MS]

205.4

-

164.8

-

-

[BMIM][NO3]

202.7

24.9

160.0

1.4

162.3

[EMIM][NPf2]

169.7

24.9

143.9

6.0

135.8

[BMIM][NPf2]

175.7

-

151.5

-

-

[EMIM][NTf2]

186.1

38.9

146.9

9.6

134.0

[BMIM][NTf2]

190.4

41.3

154.8

14.9

134.7

[EMIM][PF6]a

188.3

34.5

138.9

0.8

140.0

[BMIM][PF6]

193.8

32.3

148.7

1.5

146.5

[OMIM][PF6]

209.2

35.1

166.9

7.8

154.9

[BMIM][PROP]

196.2

61.1

157.7

29.5

121.8

[EMIM][SCN]

178.3

16.6

125.9

17.7

152.9

[BMIM][SCN]

195.0

31.7

147.7

0.3

148.1

[EMIM][TCM]

158.7

14.5

121.6

12.2

138.6

[BMIM][TCM]

165.7

15.8

130.5

8.8

143.1

[EMIM][TfO]

190.0

37.7

148.1

7.4

137.9

[BMIM][TfO]

193.3

36.7

154.4

9.2

141.4

Ionic Liquid

MAE (%) a

33.8

11.2

b

Simulation at 333K. Weighted experimental average; all individual values, uncertainties, methodology, and references given in the Supporting Information.


Page 15 of 46

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


Figure 3. Calculated versus experimental heat of vaporization (DHvap) for 1-alkyl-3methylimidazolium [RMIM]-based ionic liquids at 298 K and 1 atm. Computed values for unscaled OPLS-2009IL shown as black squares and 0.8-charge scaling shown as red triangles.

Due to the low and often undetectable vapor pressures of ionic liquids at room temperature,92,93 most experimental data is determined at elevated temperatures, ca. 400–500 K, and adjusted to 298 K. Consequently, ionic liquids with moderate thermal stabilities, e.g., containing [BF4] and [PF6] anions, can readily decompose to yield spurious DHvap values. Recent work by Zaitsau et al. took great care in accurately measuring the DHvap for [RMIM][PF6] ionic liquids by employing a quartz-crystal microbalance (QCM) technique.94 Molecular dynamics simulations by the same authors provided further validation of the experimental work. As a further test of OPLS-2009IL, the Zaitsau et al. method for computing DHvap was performed here 15 ACS Paragon Plus Environment


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 16 of 46

and reported in Table 4. Briefly, MD was performed for [EMIM][PF6], [BMIM][PF6], and [OMIM][PF6] at the average temperature, Tav, of the elevated temperature range measured during the QCM experiments. To adjust the DHvap(Tav) values down to 298.15 K, eq. 6 was used, where '

'

@ * ∆* 𝐶%,? is the molar heat capacity difference between the gas and liquid phases, i.e., 𝐶%,? − 𝐶%,? . '

@ ∆𝐻#$% 298.15 𝐾 = ∆𝐻#$% 𝑇$# + ∆* 𝐶%,? (298.15 − 𝑇$# )

(6)

'

@ The experimental values of ∆* 𝐶%,? were taken directly from Zaitsau et al.94 The difference in %

error when computing the DHvap directly at 298 K (Table 3) or adjusting DHvap(Tav) back to 298 K from the elevated temperature simulation (Table 4) are minor for the [EMIM][PF6] and [BMIM][PF6] ionic liquids. However, a considerable improvement was found for [OMIM][PF6], with an error of 2.6% for the elevated temperature MD simulations compared to 7.8% for the direct 298 K simulation; the improved accuracy may be a consequence of enhanced sampling at 433.7 K of the longer alkyl side-chain.

Table 4. Heat of Vaporization (kJ/mol) Using the Elevated Temperature Method, including the @ Average Simulation Temperature (Tav) and the Heat Capacity Adjustment (∆*' 𝐶%,? ).

Ionic Liquid

Tav (K)

DHvap (Tav)

'

@ ∆* 𝐶%,? (J/mol.K)

DHvap (298.15 K)

Exptl.a

% Error

2009IL [EMIM][PF6]

435.2

178.3

-74

188.4

140.0 ± 2.8

34.6

[BMIM][PF6]

425.7

182.0

-74

191.4

146.5 ± 2.6

30.6

[OMIM][PF6]

433.7

193.7

-85

205.2

154.9 ± 2.8

32.5

0.8*2009IL [EMIM][PF6]

435.2

128.9

-74

139.1

140.0 ± 2.8

0.6

[BMIM][PF6]

425.7

134.9

-74

144.3

146.5 ± 2.6

1.5

[OMIM][PF6]

433.7

148.2

-85

159.0

154.9 ± 2.8

2.6

a

QCM method.94


Page 17 of 46

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


Heat Capacity. Experimental heat capacities for room-temperature ionic liquids have been critically evaluated and many are available for comparison.95 Heat capacity at constant pressure (Cp) can be calculated using classical statistical mechanics (eq. 7) by monitoring the enthalpy fluctuations throughout the simulation trajectory.96 𝐶%K*$(( =

𝜕𝐻M 𝑘O 𝑇 M

(7) P

In equation 7, H is the enthalpy, T is the temperature of the system, and kb is the Boltzmann’s constant. However, using a purely classical description, Cpclass, will substantially overestimate the value at room temperature, due to the vibrational contributions being treated classically as opposed to quantum mechanically.97 A two phase thermodynamic method can be used to correct for quantum mechanical effects by considering the distribution of normal modes, i.e., density of states (DoS), as a function of the normal-mode frequency n, which is calculated using the Fourier transform of the mass-weighted atomic velocity autocorrelation functions.97-100 The distribution is normalized to the total number of degrees of freedom present in the system and each normal mode is approximated as a quantum mechanical harmonic oscillator with frequency n. The quantum mechanical correction to the classical heat capacity can then be approximated by eq. 8. ,?

𝛿𝐶R

= 𝑘O

X @

𝑊 𝜈 − 1 𝐷𝑜𝑆 𝜈 𝑑𝜈

(8)

The W is a weighting function; a full derivation of the DoS method is available.100 The final expression for the quantum-corrected heat capacity is given by eq. 9. This treatment has provided accurate Cp results in previous reports, including for water,97 ice,101 and the [BMIM][Cl] ionic liquid.102 ,?

𝐶%K1YY = 𝐶%K*$(( + 𝛿𝐶R


(9)


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 46

The quantum mechanically corrected heat capacities, Cpcorr, were calculated here using a procedure similar to that described by Chen et al.102 To summarize, the Cpclass was computed from enthalpy fluctuations using a 50 ns long trajectory, with the initial 10 ns discarded for equilibration. Since the simulations employed bond constraints, and each harmonic bond potential contributes kb to the heat capacity, a term Nckb was added to Cpclass where Nc is the number of constraints. Normal mode distributions were then calculated from the autocorrelation functions using the g_dos utility program in GROMACS.100 The simulations used for the DoS calculations were extended to 100 ps and were run in the NVT ensemble with fully-flexible bonds. All atomic velocities were saved every 4 fs and used a basic time step of 0.5 fs. The individual Cpclass and dCVqm values computed in each ionic liquid simulation are given in the Supporting Information Table S6. The final quantum corrected heat capacities, Cpcorr, are reported within a temperature range of ca. 360-400 K (Table 5). The heat capacities computed using the OPLS-2009IL force field gave reasonable agreement with experimental measurements. The overall MAE for the OPLS-2009IL unscaled and 0.8-scaled versions were 11.7% and 12.5%, respectively. The MAE’s for the ionic liquids computed here compared well to previously simulated Cp values for [BMIM][Cl], where three different force fields gave errors within 10% of experiment for both the unscaled and charge-scaled parameters.102 It should be noted that several of the experimental data points in Table 5 were reported to have large uncertainties, e.g., [EMIM][PF6] listed a Cp of 320 ± 130 J/mol.K (see Supporting Information for all experimental uncertainty values). The computed errors in Cp are also comparable to widely used water models, e.g., Spregner et al. reported a 9.9% and 18.7% error in heat capacity predictions for TIP3P and SPC/E, respectively.75


Page 19 of 46

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


Table 5. Calculated (Quantum Corrected) and Experimental Heats Capacities (J/mol.K) for 1Alkyl-3-methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

T (K)a

2009IL

% error

0.8*2009IL

% error

Exptl.b

[BMIM][ACE]

403.2

463.2

0.5

405.2

12.1

461.0

[EMIM][BF4]

358.2

354.7

7.3

321.1

2.9

330.7

[BMIM][BF4]

370.0

436.6

9.0

427.0

6.6

400.5

[OMIM][BF4]

370.0

461.0

15.2

643.8

18.4

543.8

[BMIM][Br]

403.2

463.6

26.7

497.5

35.9

366.0

[BMIM][Cl]

403.2

337.8

0.8

337.3

0.9

340.4

[EMIM][DCA]

372.2

438.0

16.3

441.2

17.1

376.7

[BMIM][DCA]

370.0

423.5

4.6

415.3

2.5

405.0

[EMIM][NTf2]

403.2

474.5

13.6

479.5

12.7

549.0

[BMIM][NTf2]

413.2

490.5

22.1

463.1

26.5

630.0

[EMIM][PF6]

403.1

289.7

9.5

277.6

13.2

320.0

[BMIM][PF6]

400.1

400.7

14.9

395.4

16.0

470.9

[OMIM][PF6]

400.0

593.1

-

580.7

-

-

[EMIM][SCN]

372.2

422.1

4.3

411.3

1.6

404.7

[BMIM][SCN]

372.2

409.8

10.4

412.2

9.9

457.3

[EMIM][TCM]

372.2

471.9

17.1

458.3

13.7

403.1

[EMIM][TfO]

400.2

336.2

19.2

347.5

16.5

416.0

[BMIM][TfO]

403.2

435.2

7.2

438.2

6.5

468.9

MAE (%) a

11.7 b

12.5

Simulation and experimental temperature. Weighted experimental average; all individual values, uncertainties, methodology, and references given in the Supporting Information.



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 46

Viscosity. The shear viscosity of an organic solvent is often computed using an equilibrium MD simulation coupled to the Green-Kubo equation, i.e., the integral over time of the pressure tensor autocorrelation function.96 However, in the case of ionic liquids, the GreenKubo method has been shown to be unreliable due to numerical errors accumulated over the long simulation times required to equilibrate the molten salts.103-105 A time decomposition correction to the Green-Kubo approach has been reported that improves the reproducibility of shear viscosity values.106 However, an alternative approach is to use a non-equilibrium periodic perturbation method shown to properly handle the highly viscous nature of ionic liquids.107,108 An excellent review by Hess on how the periodic perturbation method is implemented is available.109 Briefly summarized, simulations are carried out in 3-D periodic cells with an external force in the x direction ax(z). According to the Navier-Stokes equation, this will consequently create a velocity field, u, of the liquid (eq. 10). 𝑝

[[0

+ 𝑝 𝑢 ⋅ ∇ u = pa − ∇p + η∇M 𝑢

(10)

Since the velocity field in the y and z direction will be zero because the force is in the x direction, ux can be simplified to equation 11. 𝑝

[-b (c) [0

= 𝑝𝑎d 𝑧 + 𝜂

[ g -b (c) [c g

(11)

The velocity field can be defined by a velocity profile, V, which is calculated during a simulation. As a result, the viscosity of the system can be obtained by equation 12. 𝜂=

h % R ig

(12)

Here, k = 2p/lz, where lz is the height of the box. Also, L is the acceleration amplitude of the external force, ax(z), see equation 13. 𝑎d 𝑧 = Λcos (kz)


(13)

Page 21 of 46

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


The viscosity of a system depends greatly on the acceleration amplitude given to the external force. In order to accurately predict the viscosity, the acceleration amplitude must be small enough so that the perturbation does not disturb the equilibrium of the system. However, it must be large enough to properly sample the system. Consequently, 4 to 8 acceleration amplitudes were selected varying from 0.01–0.24 nm/ps2. For each individual amplitude, a 10 ns simulation was performed as a continuation of the production run. Each acceleration amplitude would then provide a point viscosity. Plotting the viscosities versus the amplitudes then allows for the undisturbed system where L = 0 to be extrapolated, and the intercept is taken as the calculated viscosity. The calculated viscosities for the ionic liquids yielded reasonably good agreement with experiment producing MAE values of 9.5% and 11.0% for the unscaled and 0.8scaled OPLS-2009IL parameters, respectively (Table 6).



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 46

Table 6. Calculated and Experimental Viscosities (cP) at 298 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. 2009IL

% error

0.8*2009IL

% error

Exptl.b

[BMIM][ACE]

490.5

9.7

472.8

5.7

447.3

[BMIM][AlCl4]

27.4

1.5

14.4

46.7

27.0

[EMIM][BF4]

38.2

6.8

37.7

8.0

41.0

[BMIM][BF4]

98.5

10.3

97.8

10.9

109.8

[OMIM][BF4]

357.7

6.6

337.5

0.6

335.6

[BMIM][BNZ]

196.4

-

190.7

-

-

[BMIM][Br]

271.3

18.0

231.1

0.5

230.0

[BMIM][ClO4]

209.9

16.4

177.0

1.8

180.3

[EMIM][DCA]

16.5

5.6

17.9

14.5

15.6

[BMIM][DCA]

25.9

15.1

24.6

19.3

30.5

[BMIM][HCOO]

39.3

3.4

38.0

0.0

38.0

[EMIM][MS]

88.3

6.3

74.9

9.9

83.1

[BMIM][MS]

202.8

2.7

194.7

1.4

197.4

[BMIM][NO3]

169.5

5.6

157.7

1.7

160.5

[EMIM][NPf2]

54.2

-

61.4

-

-

[BMIM][NPf2]

129.3

11.5

129.7

11.9

115.9

[EMIM][NTf2]

28.6

15.7

24.9

26.5

33.9

[BMIM][NTf2]

58.5

14.3

49.9

2.5

51.2

296.0

-

286.3

-

-

[BMIM][PF6]

287.8

2.1

264.4

6.2

282.0

[OMIM][PF6]

754.8

4.5

734.3

1.6

722.6

[BMIM][PROP]

204.7

-

180.6

-

-

[EMIM][SCN]

24.1

1.7

18.9

23.0

24.5

[BMIM][SCN]

52.3

5.9

51.2

7.9

55.6

[EMIM][TCM]

21.0

47.7

16.3

14.9

14.2

[BMIM][TCM]

28.0

1.4

21.1

23.6

27.6

[EMIM][TfO]

41.3

5.9

40.1

8.7

43.9

[BMIM][TfO]

78.7

10.1

73.7

15.8

87.5

Ionic Liquid

[EMIM][PF6]

a

MAE (%) a

9.5

11.0

b

Simulation at 333 K. Weighted experimental average; all individual values, uncertainties, methodology, and references given in the Supporting Information.


Page 23 of 46

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


Self-diffusion Coefficient. Self-diffusivity, Ds, for both the cation and anion were calculated using the Einstein relation (eq. 14) through a mean-square displacement (MSD) of the center of mass position for each ion, ri.110,111 The bracketed term is the MSD over time, N is the number of ions, and the 1/6 factor accounts for the three-dimensional system. The absolute values of the ion diffusion constants, i.e., 10-11 m2/s, are generally two orders of magnitude smaller than water and other common solvents at ambient temperatures. Equivalently, Ds can also be given by the time integral of the velocity autocorrelation function, i.e., Green-Kubo expression; however, the values generated are typically higher than that of the Einstein relation due to the same numerical imprecision noted in the viscosity calculations.104,112 p

𝐷( = lim

.

q 0→X .0

y +zp

𝑟+ 𝑡 − 𝑟+ 0

M

(14)

Unscaled OPLS-2009IL ionic liquid simulations at 298 K were found to be mostly subdiffusive, that is, plotting the MSD as a function of time on a log-log scale gave a slope < 1.111,113 To ensure the ionic liquid systems were in the proper diffusive regime for extracting diffusivities, new 40 ns trajectories at 425 K were computed for both the unscaled and 0.8-scaled OPLS-2009IL systems. The increase in temperature considerably improved the linearity of the plot between MSD and time. The computed Ds at 425 K were compared to values derived from Vogel-Fulcher-Tamman (VFT) equation (eq. 15) and the best-fit parameters of ionic diffusivity, D0 (m2/s), B (K), and T0 (K), reported by Tokuda et al. (Table 7).114 𝐷( = 𝐷@ 𝑒𝑥𝑝 −𝐵/ 𝑇 − 𝑇@

(15)

The computed MAEs for the unscaled OPLS-2009IL derived self-diffusion coefficients were particularly large at 86.1% and 87.9% for the cations and anions, respectively. Scaling the charges for the OPLS-2009IL considerably improved prediction of the experimental self23 ACS Paragon Plus Environment


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 46

diffusion coefficients with MAEs of 24.2% and 26.6% for the cations and anions. The relative difference in magnitude for the increase in the diffusion coefficients from scaling are similar to comparable nonpolarizable force fields. For example, the computed Ds for the [BMIM][PF6] solvent using the Canongia Lopes and Pádua force field115 gave D+ and D– values of 3.5 x 10-11 and 2.9 x 10-11 m2/s, respectively, at 400 K; scaling the charges by 0.8 e increased the D+ and D– values to 24.4 x 10-11 and 20.5 x 10-11 m2/s.68

Table 7. Calculated Self-Diffusion Coefficients (10-11 m2/s) at 425 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

2009IL

0.8*2009IL

VFT eq.a

D+

D–

D+

D–

D+

D–

[BMIM][BF4]

7.3

6.6

43.1

42.9

40.0

47.6

[BMIM][NTf2]

6.2

5.7

25.4

23.5

46.9

42.6

[BMIM][PF6]

4.0

3.2

29.5

26.2

31.7

27.8

[BMIM][TfO]

5.2

4.0

27.4

22.7

43.7

42.2

MAE (%)

86.1

87.9

24.2

26.6

a

Vogel-Fulcher-Tamman (VFT) equation using Tokuda et al. experimental parameters.114

Surface Tension. Surface tension simulations were performed in a fashion similar to previous work,100 where the z-axis of a pre-equilibrated ionic liquid box was extended by a factor of 3 in order to create two liquid-vapor interfaces. Constant NVT trajectories of 10 ns were generated using the v-rescale algorithm and the pressure tensor, 𝑃€ , of the system was stored at every step. The surface tension, g, was then calculated from the directional components of the pressure tensor, i.e., Px, Py, and Pz, and Lz is the length of the box in the z direction (eq. 16).116


Page 25 of 46

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


1 1 𝛾 = 𝐿c 𝑃cc − 𝑃dd + 𝑃ƒƒ 2 2

(16)

Surface tensions were only computed for ion pairs in which experimental data was reported. Previous simulations by multiple researchers have emphasized the need for elevated temperatures to accurately reproduce ionic liquid surface tensions.117-119 Overestimated values at room temperature stemmed from poor equilibrium convergence at the interface. For example, Bhargava and Balasubramanian estimated a dramatically large standard deviation at 300 K for their computed [BMIM][PF6] surface tension, i.e., g = 73 ± 222 mN/m; however, elevating the simulation temperature to 400 K yielded a more reasonable surface tension of 38 mN/m.119 With this in mind, surface tension simulations were performed at 425 K for 10 ns. Experimental comparison was carried out by choosing publications that provided g at a range of temperatures, e.g., generally 298–355 K, and extrapolating the data to 425 K. The computed ionic liquid surface tensions yielded MAEs of 33.3% and 15.6% for the unscaled and 0.8-scaled OPLS2000IL force field, respectively (Table 8). Charge scaling improved the surface tensions, which correlates to reports of enhanced accuracy when using a polarizable force field as compared to fixed charges.120 Some ion combinations yielded particularly accurate results when employing the 0.8-scaled parameters, e.g., [EMIM][BF4], [BMIM][BF4], and [BMIM][PF6] gave individual errors as low as 1.8%, 2.2%, and 1.4%, respectively (Table 8).



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 46

Table 8. Calculated and Experimental Surface Tensions (mN/m) at 425 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. 2009IL

% error

0.8*2009IL

% error

Exptl.a

[BMIM][ACE]

39.2

24.5

37.4

18.6

31.5

[EMIM][BF4]

58.8

28.1

45.1

1.8

45.9

[BMIM][BF4]

45.1

25.0

36.9

2.2

36.1

[OMIM][BF4]

31.7

30.0

27.5

12.8

24.4

[BMIM][NTf2]

44.8

69.7

37.4

41.6

26.4

[BMIM][PF6]

42.3

17.9

36.4

1.4

35.9

[BMIM][SCN]

46.6

30.1

41.4

15.7

35.8

[BMIM][TfO]

41.8

40.9

38.9

30.9

29.7

Ionic Liquid

MAE (%)

33.3

15.6

a

Experimental value extrapolated from a range of temperatures; all individual values, uncertainties, methodology, and references given in the Supporting Information.

Radial distribution functions. Hydrogen bonding between [RMIM] cations and their corresponding anions were examined through the calculation of radial distribution functions (RDFs). Of particular importance are intermolecular interactions involving the most acidic proton on the [BMIM] ring, i.e. H2-C2 bisecting the nitrogen atoms (pKa = 21-23),121,122 as it has been shown to play a predominantly large role in dictating solvent properties.26-28 Many published ionic liquid parameters have been reported that underestimate interactions involving H2-C2 when compared to ab initio calculations,45 where the QM-derived charge set utilized has a large effect

upon structural properties.24 Schröder, in particular, noted the H2-C2

underestimation problem is exacerbated when scaled charges are employed.69 In order to examine the OPLS-2009IL force field’s ability to accurately reproduce local cation-anion interactions, RDFs were calculated for the [RMIM] ring hydrogen atoms and alkyl side chain hydrogen atoms with a corresponding H-bond accepting atom of the anion. The results were


Page 27 of 46

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


compared to RDFs derived from ab initio molecular dynamics (AIMD) and polarizable force field simulations. In recent work by Kirchner and coworkers, AIMD simulations were carried out on the ionic liquid pair [BMIM][TfO].123 The AIMD simulations were performed with a periodic box of 32 ion pairs using the BLYP-D3 method and MOLOPT-DZVP-SR-GTH basis set. Radial distribution functions were computed between the [BMIM] hydrogen atoms located on the ring and alkyl side chain, and the oxygen atoms of the [TfO] anion. The AIMD calculated RDFs are given in Figure 4(A), which was adapted from “Figure 3” of Kirchner’s publication.123 Visual examination of the 0.8-scaled OPLS-2009IL computed RDF plot, Figure 4(B), clearly shows an accurate reproduction of the AIMD results. The unscaled OPLS-2009IL RDF plot is given in the Supporting Information Figure S1. From AIMD, the H2 proton on the [BMIM] ring gave the greatest intensity with a value of 2.9 and an H2-O interaction distance of 230.6 pm.123 The 0.8scaled OPLS-200IL [BMIM][TfO] simulations found an identical H2-C2 intensity of 2.9, but a longer H2-O bond distance of 258.3 pm. The two other ring hydrogen atoms (H4 and H5) also gave similar RDF peaks with intensities of 2.2 and 2.3 for AIMD and 0.8*OPLS-2009IL, respectively, and H4/5-O bond distances of 228.4 and 255.0 pm (Figure 4). In work by Schröder, a comparison between polarizable and charge-scaled nonpolarizable force fields was detailed for the ionic liquid system [EMIM][TfO].69 Molecular dynamics were performed for 1000 ion pairs at 300 K utilizing the classical force field of Canongia Lopes and Pádua,115,124 but substituting the atomic partial charges with those reported by Hanke et al.32 The polarizable ionic liquid simulations used the same force field setup and induced dipoles were modeled using Drude oscillators as implemented in CHARMM.54 The conclusion of the study was that local interactions involving the H2-C2 group at distances of less



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 46

than 800 pm particularly suffered from charge scaling. For example, RDFs between the H2 imidazolium ring hydrogen of [EMIM] and the oxygen of [TfO] found the maximum height of the peaks located at 410 pm to be approximately 2.7, 2.4, and 2.2 for the polarizable, unscaled non-polarizable, and charge-scaled non-polarizable force fields, respectively.54 Computed RDF values from the unscaled and 0.8-scaled OPLS-2009IL gave maximum height values of 3.3 and 2.6 at peaks centered around 252.5 and 257.5 pm, respectively (Figure 5). Whereas, the 0.8scaled OPLS-2009IL gave a comparable peak intensity to Schroder’s polarizable model, the predicted H2-O interacting distance was much shorter. A crystal structure of [EMIM][TfO] reported the H2-O intermolecular distance to be 224 pm.125 In addition, a gas-phase wB97XD/6311++G(d,p) geometry optimization of an [EMIM][TfO] ion pair gave a similar H2-O distance of 223 pm.126 MD simulations of [BMIM][TfO] also gave an approximate distance of 220 pm between the H2 of the cation and the nearest neighboring O of the anion22 when employing the charge-scaled non-polarizable force field of Mondal and Balasubramanian.65 These results appear to be more consistent with the average H2-O distance computed with the OPLS-2009IL between all three oxygen atoms of [TfO] and the H2 [EMIM] ring proton (Figure 5).


Page 29 of 46

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


Figure 4. Radial distribution functions between the hydrogen atoms of [BMIM] and oxygen atom of [TfO] for (A) AIMD simulations reported by Kirchner and coworkers (figure adapted from ref. 123) and (B) the 0.8-scaled OPLS-2009IL force field. 29 ACS Paragon Plus Environment


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 30 of 46

Figure 5. Radial distribution functions between the ring H2 hydrogen atom of [EMIM] and oxygen atom of [TfO] for the unscaled and 0.8-scaled OPLS-2009IL force field.

A final ion-ion structure comparison was made between OPLS-2009IL and additional AIMD simulations performed by Weber and Kirchner between ionic liquids composed of [EMIM] and cyano-based anions [DCA] and [SCN].127 Similar to Kirchner’s AIMD study of [BMIM][TfO],123 the ionic liquids [EMIM][DCA] and [EMIM][SCN] were simulated with periodic boxes composed of 32 ion pairs at 298 K using the BLYP-D3/MOLOPT-DZVP-SRGTH level of theory. As expected, the most pronounced hydrogen bonding occurs between the acidic hydrogen atom H2-C2 and the terminal nitrogen atoms of the anions. Hydrogen bonding with the other ring hydrogens, H4 and H5, were qualitatively similar. The RDF plot computed using the OPLS-2009IL between the ring H2 hydrogen atom of [EMIM] and terminal nitrogen 30 ACS Paragon Plus Environment

Page 31 of 46

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


atoms of [DCA] found large intensities of 3.2 and 2.6 for the peaks centered around 258 and 265 pm for the unscaled and 0.8-scaled OPLS-2009IL ionic liquids simulations, respectively (Figure 6). The H2-N interaction between [EMIM] and [SCN] had closer interacting distances of 235 and 248 pm with peak intensities of 3.4 and 2.6 using the unscaled and 0.8-scaled force field, respectively. The shapes, peak intensities, and distances between H2 and N with the 0.8-scaled OPLS-2009IL potentials accurately reproduced the AIMD simulations performed by Weber and Kirchner (see “Figure 2” in ref. 127), e.g. peak intensities of approximately 2.5-2.6. Radial distribution functions between the hydrogen atoms located on methyl side-chain of [EMIM] and the terminal nitrogen atoms of [DCA] and [SCN] were also reported in the AIMD study; OPLS2009IL gave very similar RDF plots and are given in the Supporting Information Figure S2.



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 32 of 46

Figure 6. Radial distribution functions for the ionic liquids [EMIM][DCA] and [EMIM][SCN] between (A) the ring H2 hydrogen atom of [EMIM] and terminal nitrogen atoms of [DCA] or [SCN] and (B) the center of ring (CoR) for two [EMIM] cations, for the unscaled (solid line) and 0.8-scaled (dashed line) OPLS-2009IL force field.

Combined distribution functions. Combined distribution functions (CDFs) can relay additional information on hydrogen bonding by monitoring both the distance, d, and angle, a, of specific intermolecular interactions occurring between the cations and anions. The d and a interactions between the H2 acidic proton of [EMIM] and the terminal nitrogens of [DCA] and [SCN] are defined in Figure 7(A). A particularly strong peak intensity (colored in red in Figure 7(C)) was located around d = 200 – 300 pm and a = 120 – 150 degrees. Smaller secondary peaks at larger distances, ~600 – 650 pm, and low angles correspond to interactions between anions the other ring hydrogen atoms, H4 and H5, at the rear of the imidazolium cation and the alkyl sidechain hydrogen atoms. The shape and general intensities of the CDF plots are similar in nature to the AIMD simulations by Weber and Kirchner,127 although they reported the angle range for the stronger peak intensities to be located at a ≈ 135 – 180 degrees. Another important effect on the local ion-ion structure is the influence of cation-anion interactions upon p-p arrangements between the cations themselves. It has been noted by Matthews et al.128 and Brussel et al.129 when an anion, like [SCN], is primarily located in the cation ring plane (i.e., in-plane conformation) and interacts with the H2 ring proton it can help stabilize p-p stacking between the cations. Alternatively, bulkier anions, such as [DCA], that instead prefer to occupy the space on-top of the cation ring (i.e., above/below ring) consequently block the plane-parallel ring-ring conformation.127 The reproduction of cation-cation interactions 32 ACS Paragon Plus Environment

Page 33 of 46

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


by OPLS-2009IL was investigated by monitoring the RDF interactions between the center-ofring (CoR) of two [EMIM] cations (Figure 6). The calculations accurately reproduced the AIMD derived RDFs from simulations by Weber and Kirchner (see “Figure 2” of their publication).127 Combined distribution functions (CDFs) were also used to investigate the orientations and distributions of the cation-cation CoR-CoR interactions throughout the 0.8-scaled OPLS2009IL MD trajectories. Figure 7(B) illustrates the angle created by a normalized vector (RN) perpendicular to the plane of a given imidazolium ring and the vector that connects two [EMIM] rings by their geometric centers as a function of distance, d. Angles, a, close to 0 and 180 degrees denote that one ring is on top of the other, i.e., p-p stacking. The 0.8-scaled force field exhibited ideal p-p stacking for [EMIM][SCN] with strong peak intensities found primarily at distances of approximately 400 pm and a ≈ 0 and 180 degrees (Figure 7(D)). The results correlate well with previous ionic liquid simulations by Kirchner detailing the center-of-rings distance between two [BMIM] rings and two [EMIM] to be ca. 400 pm.130,131 In contrast, the [EMIM][DCA] simulations, which should have reduced p-p stacking due to the bulkier anion preferring the on-top cation-anion interaction, gave strong peak intensities at several locations on the CDF plot. For example, intense CoR-CoR peaks were found for [EMIM][DCA] at distances of 400 pm and 600 pm with angles close to 0 and 180 degrees and also at distances of approximately 800 pm when the angles were between ca. 50 and 120 degrees (Figure 7(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


Page 34 of 46

Page 35 of 46

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


Figure 7. Combined distribution functions (CDF) for the ionic liquids [EMIM][DCA] and [EMIM][SCN] from the 0.8-scaled OPLS-2009IL simulations. Illustrations are given of the plotted angle a versus the distance d for (A) the H2 hydrogen atom of [EMIM] and the terminal nitrogen atoms of [DCA] or [SCN] and (B) the center of ring (CoR) interaction between two [EMIM]. The CDF plot is given with a relative intensity color for (C) the H2-N interactions and (D) the CoR-CoR interactions.

Spatial distribution functions. Figure 8 shows the spatial distribution function (SDF) of preferential locations for the ions surrounding the cation in the [BMIM][Cl] and [BMIM][BF4] trajectories when using the 0.8-scaled OPLS-2009IL force field. The SDFs were built and analyzed with the TRAVIS program84 considering isosurfaces with values of 7.0 and 29.2 particles/nm3 for the cations and anions of [BMIM][Cl], respectively, and 5.3 and 17.7 particles/nm3 for [BMIM][BF4]. The unscaled version gave similar SDF plots and can be found in the Supporting Information Figures S3 and S4. The Cl- anions preferred to occupy the area between the H2 hydrogen atom and the hydrogen atoms bound to the alkyl side chains (anion region colored in yellow in Figure 8). This interaction may be considered an in-plane hydrogen bonding interaction, as the anion interacts directly with the H2-C2 bond.25 SAPT calculations of [EMIM][Cl] ion pairs reported by Izgorodina and MacFarlane found the on-top ion-ion interaction and in-plane hydrogen bond interactions to be energetically equivalent.25 Additional anion interactions were found at the H4 and H5 locations. For the [BMIM][BF4] simulations, the BF4- anions preferred to occupy the on-top position (colored orange in Figure 8), which reproduced local structure predictions in a recent study of spatial arrangements for ionic liquids using a nearest-neighbor approach.22



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 36 of 46

Figure 8. Top and side views of the spatial distribution function of the anions (in yellow for Cland orange for BF4-) and cations (in blue) around the cation for the [BMIM][Cl] and [BMIM][BF4] simulations using the 0.8-scaled OPLS-2009IL force field.

When comparing the favored locations for cation-cation interactions, the OPLS-2009IL force field gave a population located exclusively over both sides of the face of the ring. The preferred in-plane interaction for the Cl- anion relative to the imidazolium ring emphasizes an enhanced p-p stacking for smaller anions, which is also consistent with the p-p stacking found in the CDF analysis for [EMIM][SCN], (Figure 7(D)). Conversely, as the BF4- anion preferred the 36 ACS Paragon Plus Environment

Page 37 of 46

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


location over the [BMIM] ring, the p-p stacking potential was reduced in the spatial distribution function as expected (Figure 8).

Conclusions An OPLS-based force field has been reported for ionic liquids that feature the imidazolium-based [RMIM] cation, where R = E (ethyl), B (butyl), or O (octyl), in combination with 18 different anions. The OPLS-2009IL cation and anion parameters were validated versus an extensive amount of reported experimental literature, including solvent densities, heats of vaporization, viscosities, self-diffusivities, heat capacities, and surface tensions. A charge scaling of ± 0.8 e dramatically improved the reproduction of experimental DHvap, surface tension, and the self-diffusion coefficient values. In addition to bulk-solvent properties, the 0.8-scaled OPLS2009IL also accurately reproduced local ion-ion intermolecular interactions, e.g., cation-anion interactions and p-p stacking between the rings, when compared to both ab initio MD simulations and polarizable force field results. Spatial distribution functions highlighted preferential locations for the ions surrounding the cation to be similar to that of high-level ab initio gas phase results. Going forward, it is recommended that researchers employing using our original OPLS-2009IL parameters59 use instead the current 0.8-scaled [RMIM] version and the newly reported bonded and nonbonded anion parameters. All parameters can be downloaded at http://github.com/orlandoacevedo/IL

Acknowledgement. Gratitude is expressed to the National Science Foundation (CHE-1562205) and the Center for Computational Science at the University of Miami for support of this research.



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 38 of 46

Supporting Information Available: OPLS-AA bonded and nonbonded parameters for the [RMIM] cation and all anions discussed; additional figures for RDFs of [BMIM][TfO], [EMIM][DCA], and [EMIM][SCN]; additional figures for SDFs of [BMIM][Cl] and [BMIM][BF4] for the RMIM cation. This material is available free of charge via the Internet at http://pubs.acs.org. References (1) Freemantle, M. An introduction to ionic liquids; RSC Publishing: Cambridge, UK, 2009. (2) Rogers, R. D.; Seddon, K. R. Ionic liquids--solvents of the future. Science 2003, 31, 792793. (3) Hayes, R.; Warr, G. G.; Atkin, R. Structure and nanostructure in ionic liquids. Chem. Rev. 2015, 115, 6357–6426. (4) Moschovi, A. M.; Dracopoulos, V.; Nikolakis, V. Inter- and intramolecular interactions in imidazolium protic ionic liquids. J. Phys. Chem. B 2014, 118, 8673−8683. (5) Castner, E. W.; Wishart, J. F.; Shirota, H. Intermolecular dynamics, interactions, and solvation in ionic liquids. Acc. Chem. Res. 2007, 40, 1217-1227. (6) Iwata, K.; Okajima, H.; Saha, S.; Hamaguchi, H. Local structure formation in alkylimidazolium-based ionic liquids as revealed by linear and nonlinear raman spectroscopy. Acc. Chem. Res. 2007, 40, 1174-1181. (7) Ho, T. D.; Zhang, C.; Leandro W. Hantao; Anderson, J. L. Ionic liquids in analytical chemistry: Fundamentals, advances, and perspectives. Anal. Chem. 2014, 86, 262−285. (8) Olivier-Bourbigou, H.; Magna, L.; Morvan, D. Ionic liquids and catalysis: Recent progress from knowledge to applications. Appl. Catal., A 2010, 373, 1-56. (9) Weingaertner, H. Understanding ionic liquids at the molecular level: Facts, problems, and controversies. Angew. Chem. Int. Ed. 2008, 47, 654-670. (10) Anderson, J. L.; Ding, J.; Welton, T.; Armstrong, D. W. Characterizing ionic liquids on the basis of multiple solvation interactions. J. Am. Chem. Soc. 2002, 124, 14247–14254. (11) Carmichael, A. J.; Seddon, K. R. Polarity study of some 1-alkyl-3-methylimidazolium ambient-temperature ionic liquids with the solvatochromic dye, nile red. J. Phys. Org. Chem. 2000, 13, 591–595. (12) Niedermeyer, H.; Hallett, J. P.; Villar-Garcia, I. J.; Hunta, P. A.; Welton, T. Mixtures of ionic liquids. Chem. Soc. Rev. 2012, 41, 7780-7802. (13) Chatel, G.; Pereira, J. F. B.; Debbeti, V.; Wanga, H.; Rogers, R. D. Mixing ionic liquids – “simple mixtures” or “double salts”? Green Chem. 2014, 16, 2051-2083. (14) Seddon, K. R. Ionic liquids for clean technology. J. Chem. Tech. Biotech. 1997, 68, 351356. (15) Welton, T. Room-temperature ionic liquids. Solvents for synthesis and catalysis. Chem. Rev. 1999, 99, 2071-2083. (16) Hallett, J. P.; Welton, T. Room-temperature ionic liquids: Solvents for synthesis and catalysis. 2. Chem. Rev. 2011, 111, 3508–3576. 38 ACS Paragon Plus Environment

Page 39 of 46

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


(17) Plechkova, N. V.; Seddon, K. R. Applications of ionic liquids in the chemical industry. Chem. Soc. Rev. 2008, 37, 123-150. (18) Chen, S.; Zhang, S.; Liu, X.; Wang, J.; Wang, J.; Dong, K.; Suna, J.; Xu, B. Ionic liquid clusters: Structure, formation mechanism, and effect on the behavior of ionic liquids. Phys. Chem. Chem. Phys. 2014, 16, 5893-5906. (19) Katsyuba, S. A.; Vener, M. V.; Zvereva, E. E.; Fei, Z.; Scopelliti, R.; Laurenczy, G.; Yan, N.; Paunescu, E.; Dyson, P. J. How strong is hydrogen bonding in ionic liquids? Combined x-ray crystallographic, infrared/raman spectroscopic, and density functional theory study. J. Phys. Chem. B 2013, 117, 9094–9105. (20) Fumino, K.; Wittler, K.; Ludwig, R. The anion dependence of the interaction strength between ions in imidazolium-based ionic liquids probed by far-infrared spectroscopy. J. Phys. Chem. B 2012, 116, 9507–9511. (21) Roth, C.; Peppel, T.; Fumino, K.; Köckerling, M.; Ludwig, R. The importance of hydrogen bonds for the structure of ionic liquids: Single-crystal x-ray diffraction and transmission and attenuated total reflection spectroscopy in the terahertz region. Angew. Chem. Int. Ed. 2010, 49, 10221–10224. (22) Marekha, B. A.; Koverga, V. A.; Chesneau, E.; Kalugin, O. N.; Takamuku, T.; Jedlovszky, P. l.; Idrissi, A. Local structure in terms of nearest-neighbor approach in 1‑butyl-3methylimidazolium-based ionic liquids: MD simulations. J. Phys. Chem. B 2016, 120, 50295041. (23) Skarmoutsos, I.; Welton, T.; Hunt, P. A. The importance of timescale for hydrogen bonding in imidazolium chloride ionic liquids. Phys. Chem. Chem. Phys. 2014, 16, 3675-3685. (24) Kohagen, M.; Brehm, M.; Thar, J.; Zhao, W.; Müller-Plathe, F.; Kirchner, B. Performance of quantum chemically derived charges and persistence of ion cages in ionic liquids. A molecular dynamics simulations study of 1-n-butyl-3-methylimidazolium bromide. J. Phys. Chem. B 2011, 115, 693–702. (25) Izgorodina, E. I.; MacFarlane, D. R. Nature of hydrogen bonding in charged hydrogenbonded complexes and imidazolium-based ionic liquids. J. Phys. Chem. B 2011, 115, 14659– 14667. (26) Peppel, T.; Roth, C.; Fumino, K.; Paschek, D.; Kockerling, M.; Ludwig, R. The influence of hydrogen-bond defects on the properties of ionic liquids. Angew. Chem. Int. Ed. 2011, 50, 6661–6665. (27) Kempter, V.; Kirchner, B. The role of hydrogen atoms in interactions involving imidazolium-based ionic liquids. J. Mol. Struct. 2010, 972, 22-34. (28) Thar, J.; Brehm, M.; Seitsonen, A. P.; Kirchner, B. Unexpected hydrogen bond dynamics in imidazolium-based ionic liquids. J. Phys. Chem. B 2009, 113, 15129–15132. (29) Acevedo, O.; Jorgensen, W. L.; Evanseck, J. D. Elucidation of rate variations for a DielsAlder reaction in ionic liquids from QM/MM simulations. J. Chem. Theory. Comput. 2007, 3, 132-138. (30) Acevedo, O. Determination of local effects for chloroaluminate ionic liquids on DielsAlder reactions. J. Mol. Graphics Modell. 2009, 28, 95-101. (31) Chiappe, C.; Malvaldi, M.; Pomelli, C. S. Ab initio study of the Diels−Alder reaction of cyclopentadiene with acrolein in a ionic liquid by ks-dft/3d-rism-kh theory. J. Chem. Theory Comput. 2010, 6, 179–183. (32) Hanke, C. G.; Price, S. L.; Lynden-Bell, R. M. Intermolecular potentials for simulations of liquid imidazolium salts. Mol. Phys. 2001, 99, 801-809. 39 ACS Paragon Plus Environment


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 40 of 46

(33) Shah, J. K.; Brennecke, J. F.; Maginn, E. J. Thermodynamic properties of the ionic liquid 1n-butyl-3-methylimidazolium hexafluorophosphate from monte carlo simulations. Green Chem. 2002, 4, 112-118. (34) Morrow, T. I.; Maginn, E. J. Molecular dynamics study of the ionic liquid 1-n-butyl-3methylimidazolium hexafluorophosphate. J. Phys. Chem. B 2002, 106, 12807–12813. (35) Margulis, C. J.; Stern, H. A.; Berne, B. J. Computer simulation of a “green chemistry” room-temperature ionic solvent. J. Phys. Chem. B 2002, 106, 12017–12021. (36) Andrade, J. d.; Böes, E. S.; Stassen, H. A force field for liquid state simulations on room temperature molten salts: 1-ethyl-3-methylimidazolium tetrachloroaluminate. J. Phys. Chem. B 2002, 106, 3546–3548. (37) Andrade, J. d.; Böes, E. S.; Stassen, H. Computational study of room temperature molten salts composed by 1-alkyl-3-methylimidazolium cations force-field proposal and validation. J. Phys. Chem. B 2002, 106, 13344–13351. (38) Salanne, M. Simulations of room temperature ionic liquids: From polarizable to coarsegrained force fields. Phys. Chem. Chem. Phys. 2015, 17, 14270-14279. (39) Kirchner, B.; Hollóczki, O.; Lopes, J. N. C.; Pádua, A. A. H. Multiresolution calculation of ionic liquids. WIREs Comput. Mol. Sci. 2015, 5, 202-214. (40) Zahn, S.; Brehm, M.; Brüssel, M.; Hollóczki, O.; Kohagen, M.; Lehmann, S.; Malberg, F.; Pensado, A. S.; Schöppke, M.; Weber, H.; Kirchner, B. Understanding ionic liquids from theoretical methods. J. Mol. Liq. 2014, 192, 71-76. (41) Dommert, F.; Wendler, K.; Berger, R.; Site, L. D.; Holm, C. Force fields for studying the structure and dynamics of ionic liquids: A critical review of recent developments. ChemPhysChem 2012, 13, 1625-1637. (42) Wendler, K.; Dommert, F.; Zhao, Y. Y.; Berger, R.; Holmb, C.; Site, L. D. Ionic liquids studied across different scales: A computational perspective. Faraday Discuss. 2012, 154, 111132. (43) Zahn, S.; Kirchner, B. Uncovering molecular secrets of ionic liquids. Chem. Modell. 2012, 9, 1-24. (44) Yan, T.; Burnham, C. J.; Del-Pópolo, M. G.; Voth, G. A. Molecular dynamics simulation of ionic liquids: The effect of electronic polarizability. J. Phys. Chem. B 2004, 108, 11877-11881. (45) Zahn, S.; Cybik, R. Comparison of four ionic liquid force fields to an ab initio molecular dynamics simulation. Am. J. Nano. Res. Appl. 2014, 2, 19-26. (46) Dommert, F.; Schmidt, J.; Qiao, B.; Zhao, Y.; Krekeler, C.; Site, L. D.; Berger, R.; Holm, C. A comparative study of two classical force fields on statics and dynamics of [emim][bf4] investigated via molecular dynamics simulations. J. Chem. Phys. 2008, 129, 224501. (47) Maginn, E. J. Molecular simulation of ionic liquids: Current status and future opportunities. J. Phys.: Condens. Matter 2009, 21, 373101. (48) Krekeler, C.; Schmidt, J.; Zhao, Y. Y.; Qiao, B.; Berger, R.; Holm, C.; Site, L. D. Study of 1,3-dimethylimidazolium chloride with electronic structure methods and force field approaches. J. Chem. Phys. 2008, 129, 174503. (49) Zahn, S.; Thar, J.; Kirchner, B. Structure and dynamics of the protic ionic liquid monomethylammonium nitrate ([ch3nh3][no3])([ch3nh3][no3]) from ab initio molecular dynamics simulations. J. Chem. Phys. 2010, 132, 124506. (50) Bhargava, B. L.; Balasubramanian, S. Intermolecular structure and dynamics in an ionic liquid: A car–parrinello molecular dynamics simulation study of 1,3-dimethylimidazolium chloride. Chem. Phys. Lett. 2006, 417, 486-491. 40 ACS Paragon Plus Environment

Page 41 of 46

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


(51) Bühl, M.; Chaumont, A.; Schurhammer, R.; Wipff, G. Ab initio molecular dynamics of liquid 1,3-dimethylimidazolium chloride. J. Phys. Chem. B 2005, 109, 18591–18599. (52) Pópolo, M. G. D.; Lynden-Bell, R. M.; Kohanoff, J. Ab initio molecular dynamics simulation of a room temperature ionic liquid. J. Phys. Chem. B 2005, 109, 5895–5902. (53) Borodin, O. Polarizable force field development and molecular dynamics simulations of ionic liquids. J. Phys. Chem. B 2009, 113, 11463-11478. (54) Schroder, C.; Steinhauser, O. Simulating polarizable molecular ionic liquids with drude oscillators. J. Chem. Phys. 2010, 133, 154511. (55) Bedrov, D.; Borodin, O.; Li, Z.; Smith, G. D. Influence of polarization on structural, thermodynamic, and dynamic properties of ionic liquids obtained from molecular dynamics simulations. J. Phys. Chem. B 2010, 114, 4984–4997. (56) McDaniel, J. G.; Choi, E.; Son, C. Y.; Schmidt, J. R.; Yethiraj, A. Ab initio force fields for imidazolium-based ionic liquids. J. Phys. Chem. B 2016, 120, 7024-7036. (57) Son, C. Y.; McDaniel, J. G.; Schmidt, J. R.; Cui, Q.; Yethiraj, A. First-principles united atom force field for the ionic liquid bmim+bf4–: An alternative to charge scaling. J. Phys. Chem. B 2016, 120, 3560–3568. (58) Choi, E.; McDaniel, J. G.; Schmidt, J. R.; Yethiraj, A. First-principles, physically motivated force field for the ionic liquid [bmim][bf4]. J. Phys. Chem. Lett. 2014, 5, 2670−2674. (59) Sambasivarao, S. V.; Acevedo, O. Development of OPLS-AA force field parameters for 68 unique ionic liquids. J. Chem. Theory Comput. 2009, 5, 1038-1050. (60) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19-25. (61) Jiang, W.; Wang, Y. T.; Yan, T. Y.; Voth, G. A. A multiscale coarse-graining study of the liquid/vacuum interface of room-temperature ionic liquids with alkyl substituents of different lengths. J Phys Chem C 2008, 112, 1132-1139. (62) Yan, T.; Wang, Y.; Knox, C. On the dynamics of ionic liquids: Comparisons between electronically polarizable and nonpolarizable models ii. J. Phys. Chem. B 2010, 114, 6886-6904. (63) KoBmann, S.; Thar, J.; Kirchner, B.; Hunt, P. A.; Welton, T. Cooperativity in ionic liquids. J. Chem. Phys. 2006, 124, 174506. (64) Kirchner, B.; Malberg, F.; Firaha, D. S.; Hollóczki, O. Ion pairing in ionic liquids. J. Phys.: Condens. Matter 2015, 27, 463002. (65) 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. (66) Marin-Rimoldia, E.; Shahb, J. K.; Maginn, E. J. Monte carlo simulations of water solubility in ionic liquids: A force field assessment. Fluid Phase Equilibr 2016, 407, 117-125. (67) Chaban, V. Polarizability versus mobility: Atomistic force field for ionic liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055–16062. (68) Bhargava, B. L.; Balasubramanian, S. Refined potential model for atomistic simulations of ionic liquid [bmim][pf6]. J. Chem. Phys. 2007, 127, 114510. (69) Schroder, C. Comparing reduced partial charge models with polarizable simulations of ionic liquids. Phys. Chem. Chem. Phys. 2012, 14, 3089-3102. (70) Acevedo, O. Simulating chemical reactions in ionic liquids using QM/MM methodology. J. Phys. Chem. A 2014, 118, 11653-11666.



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 42 of 46

(71) Acevedo, O.; Jorgensen, W. L. Quantum and molecular mechanical (QM/MM) monte carlo techniques for modeling condensed-phase reactions. WIREs Comput. Mol. Sci. 2014, 4, 422-435. (72) Acevedo, O.; Jorgensen, W. L. Advances in quantum and molecular mechanical (QM/MM) simulations for organic and enzymatic reactions. Acc. Chem. Res. 2010, 43, 142-151. (73) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the opls allatom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. (74) Cornell, W. D.; 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 second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. (75) Sprenger, K. G.; Jaeger, V. W.; Pfaendtner, J. The general amber force field (gaff) can accurately predict thermodynamic and transport properties of many ionic liquids. J. Phys. Chem. B 2015, 119, 5882−5895. (76) Jaguar, version 7.5, Schrödinger, LLC, New York, NY, 2008. (77) Saebø, S.; Pulay, P. Local treatment of electron correlation. Ann. Rev. Phys. Chem. 1993, 44, 213-236. (78) Saebø, S.; Tong, W.; Pulay, P. Efficient elimination of basis set superposition errors by the local correlation method: Accurate ab initio studies of the water dimer. J. Chem. Phys. 1993, 98, 2170-2175. (79) Dunning Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007-1023. (80) Jorgensen, W. L.; Tirado-Rives, J. Molecular modeling of organic and biomolecular systems using boss and mcpro. J. Comput. Chem. 2005, 26, 1689-1700. (81) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 21572164. (82) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (83) Fernández-Navarro, J. J.; García-Álvarez-Coque, M. C.; Ruiz-Ángel, M. J. The role of the dual nature of ionic liquids in the reversed-phase liquid chromatographic separation of basic drugs. J. Chromatogr. A 2011, 1218, 398-407. (84) Brehm, M.; Kirchner, B. Travis - a free analyzer and visualizer for monte carlo and molecular dynamics trajectories. J. Chem. Inf. Model. 2011, 51, 2007-2023. (85) Dommert, F.; Wendler, K.; Berger, R.; Delle Site, L.; Holm, C. Force fields for studying the structure and dynamics of ionic liquids: A critical review of recent developments. ChemPhysChem 2012, 13, 1625-1637. (86) Lopes, J. N. C.; Pádua, A. A. H. Molecular force field for ionic liquids composed of triflate or bistriflylimide anions. J. Phys. Chem. B 2004, 108, 16893–16898. (87) Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A.; Paulechka, Y. U.; Tschersich, A.; Verevkin, S. P.; Heintz, A. Experimental vapor pressures of 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imides and a correlation scheme for estimation of vaporization enthalpies of ionic liquids. J. Phys. Chem. A 2006, 110, 7303–7306. (88) Swiderski, K.; McLean, A.; Gordon, C. M.; Vaughan, D. H. Estimates of internal energies of vaporisation of some room temperature ionic liquids. Chem. Commun. 2004, 2178-2179.


Page 43 of 46

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


(89) Kabo, G. J.; Paulechka, Y. U.; Zaitsau, D. H.; Firaha, A. S. Prediction of the enthalpies of vaporization for room-temperature ionic liquids: Correlations and a substitution-based additive scheme. Thermochim. Acta 2015, 609, 7-19. (90) Xu, A.; Wang, J.; Zhang, Y.; Chen, Q. Effect of alkyl chain length in anions on thermodynamic and surface properties of 1-butyl-3-methylimidazolium carboxylate ionic liquids. Ind. Eng. Chem. Res. 2012, 51, 3458-3465. (91) Chen, Y.; Zhang, Y.; Ke, F.; Zhou, J.; Wang, H.; Liang, D. Solubility of neutral and charged polymers in ionic liquids studied by laser light scattering. Polymer 2011, 52, 481-488. (92) Esperanca, J. M. S. S.; Lopes, J. N. C.; Tariq, M.; Santos, L. M. N. B. F.; Magee, J. W.; Rebelo, L. P. N. Volatility of aprotic ionic liquids — a review. J. Chem. Eng. Data 2010, 55, 312. (93) Verevkin, S. P. Predicting enthalpy of vaporization of ionic liquids: A simple rule for a complex property. Angew. Chem. Int. Ed. 2008, 47, 5071-5074. (94) Zaitsau, D. H.; Yermalayeu, A. V.; Emel’yanenko, V. N.; Butler, S.; Schubert, T.; Verevkin, S. P. Thermodynamics of imidazolium-based ionic liquids containing pf6 anions. J. Phys. Chem. B 2016, 120, 7949–7957. (95) Paulechka, Y. U. Heat capacity of room-temperature ionic liquids: A critical review. J. Phys. Chem. Ref. Data 2010, 39, 033108. (96) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Oxford University Press: UK, 1987. (97) Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water J. Chem. Phys. 1983, 79, 23752389. (98) Lin, S.-T.; Blanco, M.; Goddard III, W. A. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. J. Chem. Phys. 2003, 119, 11792. (99) Pascal, T. A.; Linc, S.-T.; Goddard III, W. A. Thermodynamics of liquids: Standard molar entropies and heat capacities of common solvents from 2PT molecular dynamics. Phys. Chem. Chem. Phys. 2011, 13, 169-181. (100) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force field benchmark of organic liquids: Density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory Comput. 2012, 8, 61-74. (101) Waheed, Q.; Edholm, O. Quantum corrections to classical molecular dynamics simulations of water and ice. J. Chem. Theory Comput. 2011, 7, 2903–2909. (102) Chen, M.; Pendrill, R.; Widmalm, G. r.; Brady, J. W.; Wohlert, J. Molecular dynamics simulations of the ionic liquid 1‑n‑butyl-3-methylimidazolium chloride and its binary mixtures with ethanol. J. Chem. Theory Comput. 2014, 10, 4465−4479. (103) Rey-Castro, C.; Vega, L. F. Transport properties of the ionic liquid 1-ethyl-3methylimidazolium chloride from equilibrium molecular dynamics simulation. The effect of temperature. J. Phys. Chem. B 2006, 110, 14426–14435. (104) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najaf, B. Molecular dynamics simulation of imidazolium-based ionic liquids. II. Transport coefficients. J. Chem. Phys. 2009, 130, 014703. (105) Atilhan, M.; Jacquemin, J.; Rooney, D.; Khraisheh, M.; Aparicio, S. Viscous behavior of imidazolium-based ionic liquids. Ind. Eng. Chem. Res. 2013, 52, 16774–16785.



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 44 of 46

(106) Zhang, Y.; Otani, A.; Maginn, E. J. Reliable viscosity calculation from equilibrium molecular dynamics simulations: A time decomposition method. J. Chem. Theory Comput. 2015, 11, 3537-3546. (107) Zhao, L.; Wang, X.; Wang, L.; Sun, H. Prediction of shear viscosities using periodic perturbation method and OPLS force field. Fluid Phase Equilib. 2007, 260, 212-217. (108) Ma, J.; Zhang, Z.; Xiang, Y.; Cao, F.; Sun, H. On the prediction of transport properties of ionic liquid using 1-n-butylmethylpyridinium tetrafluoroborate as an example. Mol. Simulat. 2017, 1-11 (dx.doi.org/10.1080/08927022.08922017.01321760). (109) Hess, B. Determining the shear viscosity of model liquids from molecular dynamics simulations. J. Chem. Phys. 2002, 116, 209-217. (110) Zhou, Y.; Miller, G. H. Green−kubo formulas for mutual diffusion coefficients in multicomponent systems. J. Phys. Chem. 1996, 100, 5516–5524. (111) Liu, X.; Schnell, S. K.; Simon, J.-M.; Krüger, P.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Diffusion coefficients from molecular dynamics simulations in binary and ternary mixtures. Int. J. Therophys. 2013, 34, 1169-1196. (112) Liu, H.; Maginn, E.; Visser, A. E.; Bridges, N. J.; Fox, E. B. Thermal and transport properties of six ionic liquids: An experimental and molecular dynamics study. Ind. Eng. Chem. Res. 2012, 51, 7242–7254. (113) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najafi, B. Molecular dynamics simulation of imidazolium-based ionic liquids. I. Dynamics and diffusion coefficient. J. Chem. Phys. 2008, 129, 224508. (114) Tokuda, H.; Hayamizu, K.; Ishii, K.; Susan, M. A. B. H.; Watanabe, M. Physicochemical properties and structures of room temperature ionic liquids. 1. Variation of anionic species. J. Phys. Chem. B 2004, 108, 16593–16600. (115) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. H. Modeling ionic liquids using a systematic all-atom force field. J. Phys. Chem. B 2004, 108, 2038-2047. (116) Alejandrea, J.; Tildesley, D. J.; Chapela, G. A. Molecular dynamics simulation of the orthobaric densities and surface tension of water. J. Chem. Phys. 1995, 102, 4574-4583. (117) Pensado, A. S.; Malfreyt, P.; Pádua, A. A. H. Molecular dynamics simulations of the liquid surface of the ionic liquid 1-hexyl-3-methylimidazolium bis(trifluoromethanesulfonyl)amide: Structure and surface tension. J. Phys. Chem. B 2009, 113, 14708–14718. (118) Lynden-Bell, R. M.; Pópolo, M. G. D. Simulation of the surface structure of butylmethylimidazolium ionic liquids. Phys. Chem. Chem. Phys. 2006, 8, 949-954. (119) Bhargava, B. L.; Balasubramanian, S. Layering at an ionic liquid−vapor interface: A molecular dynamics simulation study of [bmim][pf6]. J. Am. Chem. Soc. 2006, 128, 10073– 10078. (120) Yan, T.; Li, S.; Jiang, W.; Gao, X.; Xiang, B.; Voth, G. A. Structure of the liquid−vacuum interface of room-temperature ionic liquids: A molecular dynamics study. J. Phys. Chem. B 2006, 110, 1800–1806. (121) Amyes, T. L.; Diver, S. T.; Richard, J. P.; Rivas, F. M.; Toth, K. Formation and stability of n-heterocyclic carbenes in water: The carbon acid pka of imidazolium cations in aqueous solution. J. Am. Chem. Soc. 2004, 126, 4366–4374. (122) Dupont, J.; Spencer, J. On the noninnocent nature of 1,3-dialkylimidazolium ionic liquids. Angew. Chem. Int. Ed. 2004, 43, 5296–5297.


Page 45 of 46

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


(123) Thomas, M.; Sancho Sanz, I.; Holloczki, O.; Kirchner, B. In Nic symposium 2016; Binder, K., Muller, M., Kremer, M., Schnurpfeil, A., Eds. 2016; Vol. 48, p 117-124. (124) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. H. Modeling ionic liquids using a systematic all-atom force field. J. Phys. Chem. B 2004, 108, 11250–11250. (125) Choudhury, A. R.; Winterton, N.; Steiner, A.; Cooper, A. I.; Johnson, K. A. In situ crystallization of ionic liquids with melting points below 225 uc. CrystEngComm 2006, 8, 742– 745. (126) Singh, D. K.; Rathke, B.; Kiefer, J.; Materny, A. Molecular structure and interactions in the ionic liquid 1‑ethyl-3-methylimidazolium trifluoromethanesulfonate. J. Phys. Chem. A 2016, 120, 6274−6286. (127) Weber, H.; Kirchner, B. Complex structural and dynamical interplay of cyano-based ionic liquids. J. Phys. Chem. B 2016, 120, 2471−2483. (128) Matthews, R. P.; Welton, T.; Hunt, P. A. Hydrogen bonding and π–π interactions in imidazolium-chloride ionic liquid clusters. Phys. Chem. Chem. Phys. 2015, 17, 14437-14453. (129) Brüssel, M.; Brehm, M.; Pensado, A. S.; Malberg, F.; Ramzan, M.; Stark, A.; Kirchner, B. On the ideality of binary mixtures of ionic liquids. Phys. Chem. Chem. Phys. 2012, 14, 1320413215. (130) Weber, H.; Hollóczki, O.; Pensado, A. S.; Kirchner, B. Side chain fluorination and anion effect on the structure of 1-butyl-3-methylimidazolium ionic liquids. J. Chem. Phys. 2013, 139, 084502. (131) Pensado, A. S.; Brehm, M.; Thar, J.; Seitsonen, A. P.; Kirchner, B. Effect of dispersion on the structure and dynamics of theionic liquid 1-ethyl-3-methylimidazolium thiocyanate. ChemPhysChem 2012, 13, 1845-1853.



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 of Contents (TOC) Graphic


Page 46 of 46

Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations

Recommend Documents