A Virtual Site OPLS Force Field for Imidazolium ... - ACS Publications

Subscriber access provided by - Access paid by the | UCSB Libraries

Article

A Virtual Site OPLS Force Field for Imidazolium-Based Ionic Liquids Brian Doherty, Xiang Zhong, and Orlando Acevedo J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b11996 • Publication Date (Web): 23 Feb 2018 Downloaded from http://pubs.acs.org on February 26, 2018

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

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

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

The Journal of Physical Chemistry

Revised version 2/22/18

A Virtual Site OPLS Force Field for Imidazolium-Based Ionic Liquids Brian Doherty, Xiang Zhong, and Orlando Acevedo* Department of Chemistry, University of Miami, Coral Gables, Florida 33146 E-mail: [email protected] Submitted December 5, 2017 Abstract: Molecular simulations of ionic liquids can provide deeper insight into the relationship between intermolecular interactions and macroscopic measurements for the solvents. However, many existing force fields have multiple shortcomings, including poor solvent dynamics, the underestimation of hydrogen bonding strength, and errors in solvent interactions/organization. A new force field, called OPLS-VSIL, has been developed for imidazolium-based ionic liquids featuring a novel topology incorporating a virtual site bisecting the nitrogen atoms that offloads negative charge to inside the plane of the ring. Guided by free energy of hydration calculations, an empirically-derived set of partial charges and nonbonded Lennard-Jones terms for both 1alkyl-3-methylimidazolium [RMIM] and 11 different anions provided accurate bulk phase ionic liquid properties and produced radial distribution functions nearly indistinguishable from ab initio molecular dynamics simulations. For example, overall mean absolute errors (MAE) of 3.1– 3.4% were computed for the density, heat of vaporization, and viscosity of approximately 20 different ion pair combinations. Additional physical properties, such as, self-diffusion coefficients, heat capacity, and surface tension also gave significant MAE improvements using OPLS-VSIL compared to existing fixed charge ionic liquid force fields. Local interactions, including cation-anion hydrogen bonding and p-p stacking between the imidazolium rings, were also accurately reproduced.

ACS Paragon Plus Environment

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

Introduction Imidazolium-based ionic liquids possess a molecular asymmetry built into the cation that impedes the strong charge ordering with anions found in the crystalline phase.1-4 Instead, the widely prevalent 1-alkyl-3-methylimidazolium [RMIM] cations, where R = E (ethyl), B (butyl), etc., form liquids with inorganic or organic anions at ambient temperatures.5,6 The attractive physical and chemical properties of these ionic liquids, including negligible vapor pressures at 298 K, wide electrochemical windows, and large thermal stabilities, has led to an ever increasing number of technological applications.6-12 Fundamentally, these desirable ionic liquid characteristics are derived from the solvent structure,1,13 including the local organization stemming from ion-ion interactions and specific hydrogen-bonding occurring between the cations and anions.14-16 Molecular simulations have provided a deeper insight into the relationship between intermolecular interactions and macroscopic measurements.17-24 However, shortcomings in existing ionic liquid force fields, including poor solvent dynamics25 and the underestimation of hydrogen bonding strength,26 have led to critical evaluations26,27 and a desire for more accurate potentials.16,28,29 Much of the error stems directly from the partial charges, as predictions of the ionic liquid thermodynamic and transport properties can be significantly improved through charge scaling specific to each ion pair.30 However, charge scaling can potentially degrade local intermolecular interactions at short ranges.31 An ideal charge set would provide the correct local interactions and solvent dynamics, while mimicking both the average charge screening caused by polarization and charge transfer effects.32,33

2 ACS Paragon Plus Environment

Page 2 of 42

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


According to traditional resonance theory, the imidazolium protonated cation can be represented by structure I, where only II and III are believed to contribute (Scheme 1). However, this traditional hybrid representation neglects a potential resonance contributor that includes positive charge on the 2-C atom as described by López and Méndez using Fukui function calculations and the local hard and soft acids/bases principle.34 In their work, they proposed an alternative resonance structure (Scheme 2) which may explain some of the more unusual features of imidazole chemistry, such as the ease with which the proton can be removed from the 2-C atom in neutral or basic conditions to produce imidazolium ylides.35-37 In our own OPLS-2009IL ionic liquid force field the 2-C atom was assigned a negative partial charge,38,39 which was also the case for other popular imidazolium-based force fields by the groups of Canongia Lopes and Pádua,40 Wang,41, and Balasubramanian.30

Scheme 1. Traditional resonance structures for imidazolium protonated cation.

Scheme 2. López and Méndez’s postulated resonance structure for imidazolium protonated cation. 3 ACS Paragon Plus Environment


In this work, a novel OPLS-AA force field parameterization effort has been undertaken where [RMIM] cation parameters were created that utilized a virtual site bisecting the ring nitrogen atoms. The creation of empirical charges (not QM derived) were guided by free energy of hydration calculations of an imidazolium cation in aqueous solution. Lennard-Jones terms were fine-tuned for both [RMIM] and the anions (Figure 1) to accurately reproduce bulk phase ionic liquid properties and improve local cation-anion intermolecular interactions. Overall, the new parameters yielded near quantitative reproduction of experimental data in many cases, including densities, heats of vaporization, viscosities, diffusion coefficients, surface tensions, and heat capacities. Reproducing radial distribution functions (RDFs) derived from ab initio molecular dynamics was also emphasized. This new imidazolium-based ionic liquid virtual site force field is called OPLS-VSIL and is reported herein.

Figure 1. Ionic liquid-forming ions. R = E (ethyl), B (butyl), and O (octyl). 4 ACS Paragon Plus Environment

Page 4 of 42

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


Computational Methods OPLS-AA Force Field. The OPLS-AA force field calculates the total energy as a combination of intramolecular energies, that is, the harmonic bond stretching and angle bending terms, a Fourier series for dihedral angles, and the intermolecular energies are computed using Coulomb and 12-6 Lennard-Jones terms (Eq. 1-4).42 The adjustable 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.

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)

The Lennard-Jones coefficients used the geometric combining rules, i.e., sij = (siisjj)1/2 and eij = (eiiejj)1/2. Nonbonded interactions were calculated 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. Molecular Dynamics. Molecular dynamic (MD) simulations were performed using the GROMACS 5.0.7 software package.43 Cubic boxes containing 500 ion pairs were created using the Packmol program.44 Periodic boundary conditions and Particle-Mesh Ewald summations 5 ACS Paragon Plus Environment


were applied. The systems were minimized using a steepest descent algorithm for 5000 steps. Equations of motion were integrated using the leap-frog algorithm with a time step of 1 fs. A temperature value of 298 K was kept constant using velocity rescaling with a stochastic term (vrescale)45 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. All covalent hydrogen-containing bonds were constrained using the LINCS algorithm and a cutoff range for the short-range electrostatics was set to 13 Å. Production runs were carried out for an additional 40 ns. Radial distribution functions (RDF), combined distribution functions (CDF), and spatial distribution functions (SDF) were computed using the TRAVIS program.46 Results and Discussion [RMIM] Virtual Site Topology. López and Méndez postulated an alternative resonance structure for imidazolium (Scheme 2) based on condensed Fukui function values obtained from a Mulliken charge analysis.34 Drawing inspiration from their work, a new empirically-derived set of partial charges for imidazolium was developed that incorporated a virtual site bisecting the nitrogen atoms which offloads negative charge to inside the plane of the ring. Condensed-phase simulation improvements with a charged dummy atom bisecting real atoms has precedent in the TIP4P versus TIP3P water models.47 The virtual site in the [RMIM] maintained a relatively constant in-plane position by using 3 reference atoms (X, N1, and C2) with a fixed N1-X bond distance of 1.0847 Å and fixed X-N1-C2 angle of 47.869 degrees (Figure 2). To ensure the virtual site did not appreciably drift throughout the simulations, the atomic coordinates were monitored for [BMIM][BF4] using an extended 100 ns MD simulation. For example, N1-X, C2-X, and N3-X distances were computed for all 500 [BMIM] cations in the simulation box and averaged over 6 different snapshots, i.e., at 5, 10, 20, 40, 80 and 100 ns. The resultant distances were 1.0847 ± 6 ACS Paragon Plus Environment

Page 6 of 42

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


0.0041, 0.99 ± 0.01, and 1.11 ± 0.04 for N1-X, C2-X, and N3-X, respectively. While the standard deviations in the distances were relatively minor, the slightly larger variations for the C2-X and N3-X distances can be attributed to atom movements within the imidazolium ring itself.

Figure 2. The OPLS-VSIL [RMIM] cation employing the GROMACS 3fad virtual site construction type, where d = 1.0847 Å is the fixed distance and q = 47.869 degrees is the fixed angle. The traditionally accurate predictions of conventional liquid properties48 by OPLS-AA stem largely from Jorgensen’s emphasis on the correct reproduction of the experimental densities, heats of vaporization (DHvap), and free energies of hydration (DGhyd) during development.49 While many ionic liquid parameterization efforts considered densities and DHvap, none to our knowledge used DGhyd as a criterion. Not surprisingly, a recent assessment of water solubility in ionic liquids found that all tested nonpolarizable force fields gave overestimated results.50 In this work, the OPLS-VSIL explicitly considered the reproduction of density, DHvap, and DGhyd for guidance in developing parameters for imidazolium and [RMIM] (Table 1). In addition to bulk-phase properties, development of the [RMIM] partial charges and LennardJones terms emphasized the reproduction of radial distribution functions (RDF) derived from ab initio molecular dynamics (AIMD) simulations reported by Kirchner and coworkers.51,52



Atom types, charges, Lennard-Jones terms, and the topology for OPLS-VSIL are given for both the imidazolium cation and the [RMIM] cations in Figure 3 and Table 1. The nonbonded parameters were developed specifically towards a ± 0.8 e charge scaling for both the cations and anions. A non-inter molecular charge has been reported to dramatically improve experimental ionic liquid-phase properties including enthalpy of vaporization, surface tension, and ion diffusion coefficients.30,39,50,53,54 Correct solute-solvent interactions, e.g., solubility, are also retained when utilizing charge scaled solvents. For example, Maginn computed a Henry’s Law constant of 51.4 ± 1.6 bar for CO2 in a 0.8-charge scaled [BMIM][NTf2] model55 that compared favorably to experimental measurements of 46.5 bar56 and 50.76-51.08 bar57 at 333K. However, it has been previously noted that discrepancies can exist at the local level between the scaled charges and polarizable models.31 As such, the OPLS-VSIL parameters took care in fine-tuning the Lennard-Jones parameters of the anions listed in Figure 1 towards the virtual site [RMIM] model to accurately reproduce the local intermolecular interactions (see the radial distributions discussion). The [RMIM] bonded parameters, i.e., bonds, angles, and torsions, were kept intact from our OPLS-2009IL ionic liquid force field (values are provided in the Supporting Information Tables S1 and S2).38,39 All anion parameters are given in the Supporting Information Tables S3–S5. Vibrational frequencies were not parameterized, as is characteristic of the OPLSAA force field.


Page 8 of 42

Page 9 of 42 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. Atom types for (A) 1-alkyl-3-methylimidazolium [RMIM] and (B) imidazolium cations. Table

1.

OPLS-VSIL

Charges

and

Lennard-Jones

Parameters

for

the

1-Alkyl-3-

methylimidazolium [RMIM] and Imidazolium Cations.a Atom

q (e)

s (Å)

e (kcal/mol)

CR type NA

0.040

3.55

0.07

0.176

3.25

0.17

CW

-0.056

3.55

0.07

CM

-0.200

3.50

0.066

CA

-0.064

3.50

0.066

CS

-0.096

3.50

0.066

CT

-0.192

3.50

0.066

X

-0.088

0

0

H

0.104

0

0

HR

0.128

1.85

0.03

HW

0.136

1.85

0.03

HM

0.104

2.15

0.03

HA

0.080

2.25

0.03

HS

0.048

2.50

0.03

HT

0.064

2.50

0.03

a

Atom types for [RMIM] and imidazolium given in Figure 3. Overall cation charge is +0.8 e.

Free Energy of Hydration. To compute the DGhyd for imidazolium, the BOSS program58 was utilized and a procedure similar to that reported by Jensen and Jorgensen for halide ions59 and charged amino acids60 was followed. Briefly, a TIP4P water molecule sphere of 15 Å radii was used with an effectively infinite (100 Å) nonbonded cutoff and a corresponding Born correction of 10.93 kcal/mol for size-consistent results. An individual imidazolium cation was 9 ACS Paragon Plus Environment


placed in the center of the sphere and free energies were computed by annihilating in one calculation the charge while keeping LJ parameters unaltered, and in a second perturbation annihilating LJ parameters of the uncharged molecules. Monte Carlo simulations used 10 windows of double-wide sampling, with 5 million configurations of equilibration and 10 million configurations of averaging per window, which has been reported to be sufficient for accurate results.60 The DGhyd value computed is the difference between the free energy of annihilating all nonbonded interactions in both water and gas. The partial charges for the imidazolium atom types in this case were scaled up to give an overall +1 cation charge. Calculating the free energy of hydration yielded a value of -56.7 kcal/mol, which was reasonably close to experimental estimates of -62 and -64 kcal/mol when considering the reported error bars of ± 5 kcal/mol (Table 2).61,62 Table 2. Computed Free Energy of Hydration, DGhyd (kcal/mol), for Imidazolium Cation using OPLS-VSIL. Calc. DGhyd

Exptl. DGhyd62

-56.7 ± 0.5

-62, -64 ± 5

Density. A comparison between the OPLS-VSIL calculated densities and corresponding literature values at 298 K is given in Table 3. A weighted experimental average was computed by determining the weight of each data point from the inverse of its uncertainty. All individual experiments, uncertainties, measurement methodologies, and references are given in the Supporting Information. The calculated deviation from reported densities was relatively small with an overall mean absolute error (MAE) of 3.1% for the OPLS-VSIL. The results compared favorably to our 0.8-charge scaled OPLS-2009IL force field38 that gave an MAE of 2.8%. 10 ACS Paragon Plus Environment

Page 10 of 42

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


Individual percent errors compared to the experimental weighted mean are also reported for each ion combination. Notably, the [RMIM][BF4] and [RMIM][PF6] ion combinations, which are extensively used in the literature both experimentally11,12,63,64 and computationally,19,22 gave near quantitative agreement, with errors of less than 1%, when using the OPLS-VSIL force field. Heats of Vaporization. Accurately computing vaporization enthalpies for ionic liquids can be very challenging for nonpolarizable force fields.22,25 While substantial improvement in DHvap predictions can be obtained by employing polarizable65-67 or first-principles68-70 force fields, the extensive computational resources required by both methods can present a significant barrier for many researchers. Alternatively, polarization effects can be emulated by simply scaling the atomic charges to non-integer values. 30,50,53 For example, when our OPLS-2009IL force field was recently used to compute the DHvap for 25 different ionic liquids, the overall MAE improved from a 33.8% error when using unscaled charges to 11.2% for 0.8-charge scaling.39 The reduced electrostatics mimic charge transfer effects and the average charge screening caused by polarization.32,33 In this work, the DHvap was calculated at 298 K for both the 0.8*OPLS-2009IL and OPLS-VSIL force fields using eq. 5. As ionic liquids vaporize exclusively as neutral contact ionpairs,71-73 the Etotal (gas) simulations were carried out on a single ion-pair for 10 ns using the same setup described for liquid simulations. The Etotal (liquid) term is the corresponding value in the condensed phase. The RT term is used in place of a PV-work term in the enthalpy. In the specific cases of [EMIM][PF6], [BMIM][PF6], and [OMIM][PF6] an elevated temperature method developed by Zaitsau et al. was employed.74 This is necessary, in part, because [EMIM][PF6] is solid up to temperatures of 333 K. Briefly, MD was performed for the [RMIM][PF6] solvents at the average temperature, Tav, of the elevated temperature range 11 ACS Paragon Plus Environment


Page 12 of 42

measured from quartz-crystal microbalance (QCM) experiments. Equation 6 was used to adjust #

( the DHvap(Tav) values down to 298.15 K, where ∆" 𝐶%,' is the molar heat capacity difference #

#

" ( between the gas and liquid phases, i.e., 𝐶%,' − 𝐶%,' . Experimental values of ∆" 𝐶%,' from the

Zaitsau et al. paper74 were utilized and are given along with Tav in the Supporting Information Table S6. ∆𝐻+,% = ∆𝐻#,. − ∆𝐻"/01/2 = 𝐸454," (𝑔𝑎𝑠) − 𝐸454," (𝑙𝑖𝑞𝑢𝑖𝑑) + 𝑅𝑇 #

( ∆𝐻+,% (298.15 𝐾 ) = ∆𝐻+,% (𝑇,+ ) + ∆" 𝐶%,' (298.15 − 𝑇,+ )

(5) (6)

The OPLS-VSIL yielded highly accurate vaporization enthalpies at 298 K with an overall MAE of 3.4%, a significant improvement over the 6.6% predicted from 0.8*OPLS-2009IL (Table 4). It should be noted that experimental measurement of liquid-vapor equilibrium thermodynamic data, i.e., DHvap, for ionic liquids can be extremely difficult as room temperature ionic liquids often possess low and even undetectable vapor pressures.75 Consequently, most DHvap are measured at elevated temperatures of approximately 400–500 K and adjusted down to #

( 298 K. However, experimental measurements of ∆" 𝐶%,' for correcting temperature to 298 K are

still uncommon in the literature and using a crude general value, such as -100 J/K.mol, has been shown to give poor results.75,76 As a final point of comparison, ionic liquid simulations were performed at the same elevated temperatures used in the DHvap experiments. The OPLS-VSIL yielded a slighter higher MAE of 4.5% at the higher temperatures, but the individual results were generally comparable to the room temperature simulations (Table 5).


Page 13 of 42 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 3. Calculated and Experimental Liquid Densities (g/cm3) at 298 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

0.8*2009IL

% error

VSIL

% error

Exptl.b

[EMIM][BF4]

1.212

5.2

[BMIM][BF4]

1.150

4.9

1.278 1.201

0.0 0.6

1.278 1.209

[OMIM][BF4]

1.074

3.0

1.109

0.2

1.107

[BMIM][Br]

1.277

1.5

1.345

3.8

1.296

[BMIM][Cl]

1.048

3.0

1.106

2.4

1.080

[EMIM][DCA]

1.069

3.1

1.149

4.1

1.103

[BMIM][DCA]

1.039

2.0

1.099

3.7

1.060

[EMIM][MS]

1.229

4.1

1.298

1.2

1.282

[BMIM][MS]

1.177

2.9

1.225

1.1

1.212

[BMIM][NO3]

1.155

0.1

1.206

4.3

1.156

[EMIM][NTf2]

1.584

4.4

1.579

4.0

1.517

[BMIM][NTf2]

1.504

4.7

1.509

5.1

1.436

[EMIM][PF6]a

1.363

-

1.474

-

-

[BMIM][PF6]

1.315

3.5

1.373

0.7

1.363

[OMIM][PF6]

1.207

3.5

1.247

0.3

1.251

[EMIM][SCN]

1.130

1.2

1.177

5.4

1.117

[BMIM][SCN]

1.068

0.4

1.115

4.0

1.072

[EMIM][TCM]

1.055

2.4

1.143

5.7

1.081

[BMIM][TCM]

1.028

1.8

1.098

4.9

1.046

[EMIM][TfO]

1.419

2.6

1.456

5.2

1.384

[BMIM][TfO]

1.333

2.6

1.357

4.5

1.299

MAE (%) a

2.8

3.1

b

Simulation at 333 K. Weighted experimental average; all individual values, uncertainties, methodology,

and references given in the Supporting Information.



Page 14 of 42

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

0.8*2009IL

% error

VSIL

% error

Exptl.b

[EMIM][BF4]

135.1

2.8

137.2

1.3

[BMIM][BF4]

140.5

3.0

144.5

5.9

139.0 136.4

[OMIM][BF4]

159.0

1.9

163.3

0.8

162.0

[BMIM][Br]

157.7

3.8

158.7

4.5

151.9

[BMIM][Cl]

157.7

4.4

158.7

5.1

151.0

[EMIM][DCA]

137.1

12.3

154.0

1.5

156.4

[BMIM][DCA]

140.2

11.0

160.7

2.1

157.5

[EMIM][MS]

158.1

5.7

163.6

9.4

149.6

[BMIM][MS]

164.8

-

170.6

-

-

[BMIM][NO3]

160.0

1.4

161.4

0.6

162.3

[EMIM][NTf2]

146.9

9.6

124.1

7.4

134.0

[BMIM][NTf2]

154.8

14.9

139.4

3.5

134.7

[EMIM][PF6]a

139.1

0.6

142.9

2.1

140.0

a

[BMIM][PF6]

144.3

1.5

148.2

1.2

146.5

[OMIM][PF6]a

159.0

2.6

163.5

5.6

154.9

[EMIM][SCN]

125.9

17.7

140.3

8.2

152.9

[BMIM][SCN]

147.7

0.3

147.6

0.3

148.1

[EMIM][TCM]

121.6

12.3

134.4

3.0

138.6

[BMIM][TCM]

130.5

8.8

140.3

2.0

143.1

[EMIM][TfO]

148.1

7.4

134.7

2.3

137.9

[BMIM][TfO]

154.4

9.2

142.3

0.7

141.4

MAE (%) a

6.6

3.4

b

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


Page 15 of 42 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 and Experimental Heat of Vaporization (kJ/mol) at Elevated Temperatures for 1-Alkyl-3-methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

Temp (K)a

VSIL

% error

Exptl.b

[EMIM][BF4]

431.6

126.0

3.1

[BMIM][BF4] [OMIM][BF4]

500.0 520.0

124.3 134.7

6.5 3.1

122.2 133.0 139.0

[BMIM][Br]

447.2

142.6

6.3

152.2

[BMIM][Cl]

431.3

144.3

4.5

138.0

[BMIM][DCA]

474.2

140.9

0.9

139.6

[EMIM][MS]

442.0

149.0

10.2

135.2

[EMIM][NTf2]

463.0

108.8

8.4

118.8

[BMIM][NTf2]

477.6

118.9

0.5

118.3

[EMIM][PF6]

435.2

132.8

2.2

129.9

[BMIM][PF6]

425.7

138.8

1.2

137.1

[OMIM][PF6]

433.7

152.0

6.0

143.4

[EMIM][SCN]

413.2

130.1

8.5

142.2

[EMIM][TCM]

423.2

120.7

4.2

126.0

[EMIM][TfO]

412.8

124.3

1.7

126.4

MAE (%)

4.5

a

Simulation and experimental temperature. bAll individual values, uncertainties, methodology, and references given in the Supporting Information.

Heat Capacity. The heat capacity of the system at constant pressure can be calculated using classical statistical mechanics (Cpclass) by examining enthalpy and temperature fluctuations throughout the simulation (eq. 7). In equation 7, H is the enthalpy, T is the temperature, and kb is the Boltzmann’s constant. However, neglecting the quantum mechanical vibrational contributions has been shown to substantially overestimate Cp values at room temperature.77

𝐶%K",..

〈𝜕𝐻O 〉 =L S 𝑘R 𝑇 O T


(7)


Page 16 of 42

Quantum mechanical effects can be accounted for by using a two-phase thermodynamics method called density of states (DoS). This DoS method considers the distribution of normal modes as a function of the normal-mode frequency n, which is computed48 using the Fourier transform of the mass-weighted atomic velocity autocorrelation functions.77-79 The distribution is normalized to the total number of degrees of freedom present in the system and each normal mode is estimated to be a quantum mechanical harmonic oscillator with frequency n. The quantum mechanical correction to the classical heat capacity can then be approximated by eq. 8, where W is a weighting function.48 0'

𝛿𝐶V

]

= 𝑘R ∫( (𝑊 (𝜈) − 1)𝐷𝑜𝑆(𝜈)𝑑𝜈

(8)

The final expression for the quantum-corrected heat capacity is given by eq. 9, which has provided accurate heat capacities for various systems,77,80 including ionic liquids.39,81 0'

𝐶%K5^^ = 𝐶%K",.. + 𝛿𝐶V

(9)

The quantum mechanically corrected heat capacities, Cpcorr, were computed using the final 40 ns of a 50 ns trajectory. Since the simulations utilized 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.48 The simulations were extended by 100 ps for the DoS calculations and employed the NVT ensemble with fully-flexible bonds. A basic time step of 0.5 fs was utilized with all atomic velocities saved every 4 fs. The individual Cpclass and dCVqm values computed for each solvent are provided in the Supporting Information Table S7. The final quantum corrected heat capacities, Cpcorr, are reported within a temperature range of 358.2-403.2 K (Table 6). The computed heat capacity MAE of 9.5% from 16 ACS Paragon Plus Environment

Page 17 of 42 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 OPLS-VSIL is comparable to both the 0.8-charge scaled OPLS-2009IL force field (MAE = 10.6%) and the popular TIP3P and SPC/E water models, i.e., 9.9% and 18.7% error, respectively.82

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

T (K)a

0.8*2009IL

% error

VSIL

% error

Exptl.b

[EMIM][BF4] [BMIM][BF4]

358.2 370.0

321.1 427.0

2.9 6.6

323.1 397.4

2.3 0.8

330.7 400.5

[OMIM][BF4]

370.0

643.8

18.4

595.8

9.6

543.8

[BMIM][Br]

403.2

497.5

35.9

359.3

1.8

366.0

[BMIM][Cl]

403.2

337.3

0.9

343.8

1.0

340.4

[EMIM][DCA]

372.2

441.2

17.1

308.0

18.2

376.7

[BMIM][DCA]

370.0

415.3

2.5

434.9

7.4

405.0

[EMIM][SCN]

372.2

411.3

1.6

322.7

20.3

404.7

[BMIM][SCN]

372.2

412.2

9.9

380.6

16.8

457.3

[EMIM][TCM]

372.2

458.3

13.7

373.3

7.4

403.1

[BMIM][TfO]

403.2

438.2

6.5

381.9

18.6

468.9

MAE (%)

10.6

9.5

a

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

Viscosity. The non-equilibrium periodic perturbation method has been shown to properly handle the highly viscous nature of ionic liquids.83 An outstanding methodological review of the perturbation approach for calculating shear viscosity is available.84 To summarize, an external force, a, is applied that will create a velocity field, u, of the liquid (eq. 10). Simulations are



Page 18 of 42

carried out in 3-D periodic cells with an external force chosen in the x direction such that ax is a function of z only, and ay and az are equal to zero. Consequently, uy and uz become zero as well. `𝒖

𝑝 `4 + 𝑝(𝒖 ⋅ 𝛁)𝐮 = p𝐚 − ∇p + η∇O 𝒖

(10)

With no pressure gradient in the x-direction, ux can be simplified to equation 11. The velocity field can be defined by a velocity profile, V, which is calculated during a simulation. As a result, the viscosity, h, of the system can be obtained from equation 12, where k = 2p/lz, and lz is the height of the box. 𝑝

`1i (j) `4

= 𝑝𝑎k (𝑧) + 𝜂

` n 1i (j) `j n

o %

𝜂 = V pn

(11)

(12)

To obtain a smooth velocity profile with small local shear rates, an acceleration amplitude, L, of the external force ax(z) is chosen that satisfies both conditions (eq. 13). The acceleration amplitude must be small enough so that the perturbation does not disturb the equilibrium of the system, but large enough to allow for proper sampling. 𝑎k (𝑧) = Λcos (kz)

(13)

Accordingly, 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 provides 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.


Page 19 of 42 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 OPLS-VSIL force field dramatically improved the predicted viscosities with a computed MAE of 3.2% compared to 10.4% for the 0.8-charge scaled OPLS-2009IL (Table 7). All individual percent errors were in the single digit range, with the exception of the solvent [EMIM][NTf2] that gave a 29.0% error compared to experiment. Removal of this single data point brings the overall MAE down to 1.8% for the OPLS-VSIL simulations. Table 7. Calculated and Experimental Viscosities (cP) at 298 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. 0.8*2009IL

% error

VSIL

% error

Exptl.a

[EMIM][BF4] [BMIM][BF4]

37.7 97.8

8.0 10.9

39.9 108.9

2.7 0.9

41.0 109.8

[OMIM][BF4]

337.5

0.6

348.2

3.7

335.6

[BMIM][Br]

231.1

0.5

237.7

3.4

230.0

[EMIM][DCA]

17.9

14.5

15.7

0.8

15.6

[BMIM][DCA]

24.6

19.3

30.5

0.2

30.5

[EMIM][MS]

74.9

9.9

83.3

0.3

83.1

[BMIM][MS]

194.7

1.4

195.1

1.1

197.4

[BMIM][NO3]

157.7

1.7

162.0

1.0

160.5

[EMIM][NTf2]

24.9

26.3

24.0

29.0

33.8

[BMIM][NTf2]

49.9

2.5

51.1

0.2

51.2

[BMIM][PF6]

264.4

6.2

279.4

0.9

282.0

[OMIM][PF6]

734.3

1.6

720.9

0.2

722.6

[EMIM][SCN]

18.9

23.0

24.5

0.0

24.5

[BMIM][SCN]

51.2

7.9

56.4

1.4

55.6

[EMIM][TCM]

16.3

14.9

13.9

2.2

14.2

[BMIM][TCM]

21.1

23.6

28.2

2.1

27.6

[EMIM][TfO]

40.1

8.7

40.2

8.5

43.9

[BMIM][TfO]

73.7

15.8

85.3

2.5

87.5

Ionic Liquid

MAE (%)

10.4

3.2

a

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



Page 20 of 42

Surface Tension. Surface tension was calculated by elongating the z-axis of a preequilibrated ionic liquid box by a factor of 3 in order to generate two liquid-vapor interfaces. Simulations of 10 ns were performed in the NVT ensemble using the v-rescale algorithm with the pressure tensor, 𝑃x , of the system stored at every step. The surface tension, g, was then calculated using equation 14, where Lz is the length of the box in the z direction and Px, Py, and Pz, are the directional components of the pressure tensor (eq. 14).85 1 1 𝛾 = 𝐿j {𝑃jj − |𝑃kk + 𝑃}} ~• 2 2

(14)

Simulations performed near room temperature can lead to poor equilibrium convergence at the interface and yield overestimated surface tension predictions.86-88 For example, in their surface tension calculation of [BMIM][PF6], Bhargava and Balasubramanian estimated a strikingly large standard deviation of g = 73 ± 222 mN/m at 300 K, which was considerably improved when the simulation was carried out at 400 K.88 Consequently, surface tension simulations were performed here at 425 K for 10 ns. The experimental data listed in Table 8 was determined by selecting experiments that provided g at a range of temperatures, approximately 298–355 K, and extrapolating the results to 425 K. In all cases, linear fits gave an r2 > 99%. The OPLS-VSIL force field gave an overall MAE value of 7.3%, which cut the error in half compared to the 0.8charge scaled OPLS-2009IL potentials (Table 8). The majority of the improvement stemmed primarily from significant corrections for two specific ionic liquids, i.e., [BMIM][NTf2] and [BMIM][TfO].


Page 21 of 42 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 8. Calculated and Experimental Surface Tensions (mN/m) at 425 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

0.8*2009IL

% error

VSIL

% error

Exptl.a

[EMIM][BF4]

45.1

1.8

45.2

1.5

[BMIM][BF4]

36.9

2.2

35.7

1.1

45.9 36.1

[OMIM][BF4]

27.5

12.8

27.9

14.2

24.4

[BMIM][NTf2]

37.4

41.6

29.1

10.2

26.4

[BMIM][PF6]

36.4

1.4

39.1

8.9

35.9

[BMIM][SCN]

41.4

15.7

40.7

13.7

35.8

[BMIM][TfO]

38.9

30.9

30.1

1.4

29.7

MAE (%)

15.2

7.3

a

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

Self-diffusion Coefficients. Self-diffusion coefficients, Ds, for the cations and anions of the ionic liquids were computed by applying the Einstein relation (eq. 15) and a mean-square displacement (MSD) of the center of mass position for each ion, ri.89 The Einstein relation has been found to yield more accurate results for ion self-diffusivity in the molten salts as compared to alternative methods, such as the Green-Kubo expression.90,91 To sample in the proper diffusive regime for extracting diffusivities, 40 ns trajectories were computed at an elevated temperature of 425 K to ensure a linear plot between MSD and time. The computed D+ and D– given in Table 9 were compared to values derived from Vogel-Fulcher-Tamman (VFT) equation (eq. 16) and the best-fit parameters of ionic diffusivity, D0 (m2/s), B (K), and T0 (K), reported by Tokuda and co-workers.92 €

2

𝐷. = • lim 24 〈∑• ⃗/ (𝑡) − 𝑟⃗/ (0)]O 〉 /Ž€[𝑟 4→]


(15)


𝐷. = 𝐷( 𝑒𝑥𝑝[−𝐵/(𝑇 − 𝑇( )]

Page 22 of 42

(16)

The calculated self-diffusivity coefficient MAEs for the 0.8-charge scaled OPLS-2009IL force field were 24.4% and 26.6% for the cations and anions, respectively, which are comparable to reported values from other popular charge-scaled ionic liquid force fields.54 Remarkably, the OPLS-VSIL improved the D+ prediction to a computed MAE of 9.6%. The improvement in D– from OPLS-VSIL was more modest with a MAE of 21.1% (Table 9). It should be noted that the velocity rescale thermostat was employed and under certain circumstances stochastic thermostats may influence solvent dynamics. However, given proper sampling these effects can be minimized.93 For example, the velocity rescale thermostat has accurately reproduced transport properties in multiple ionic liquid studies.94-96 Alternative thermostats, including the NoséHoover temperature coupling method, have been utilized in prior ionic liquid calculations.97 As Nosé-Hoover dynamics are nonergodic,98 a series of Nosé-Hoover chains are typically used to impart ergodicity on the thermostatted system.99 Problematically, the number of chains is limited to one in the GROMACS software when using the leap-frog algorithm. As a test of the effect of thermostat choice upon the self-diffusion coefficients, simulations employing the Nosé-Hoover method were performed at 425 K in a similar fashion to Table 9. An additional 10 ns of NVT equilibration was performed prior to the 40 ns trajectory to ensure sufficient sampling. The Nosé-Hoover computed values were underestimated compared to the velocity rescale method giving calculated self-diffusivity coefficient MAEs of 16.7% and 22.7% for the cations and anions, respectively (Supporting Information Table S8).


Page 23 of 42 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 9. Computed Self-Diffusivity Coefficients (10-11 m2/s) at 425 K for 1-Alkyl-3methylimidazolium [RMIM]-Based Ionic Liquids. Ionic Liquid

0.8*2009IL

VFT eq.a

VSIL

D+

D–

D+

D–

D+

D–

[BMIM][BF4]

43.1

42.9

43.0

30.6

40.0

47.6

[BMIM][NTf2]

25.4

23.5

47.6

45.7

46.9

42.6

[BMIM][PF6]

29.5

26.2

23.2

20.1

31.7

27.8

[BMIM][TfO]

27.4

22.7

44.8

36.5

43.7

42.2

MAE (%)

24.4

26.6

9.6

21.1

a

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

Bulk Solvent Summary. Overall, the virtual site OPLS-VSIL force field consistently reproduced experimental physical properties, i.e., density, DHvap, heat capacity, viscosity, surface tension, and self-diffusivity, with a high degree of accuracy. Notable improvements in the ionic liquid simulations compared to the traditional QM-charge derived OPLS-2009IL can been seen in Table 10. For example, the MAE for DHvap at 298 K went from 6.6% with the 0.8-charge scaled OPLS-2009IL to 3.4% for the OPLS-VSIL and viscosity improved from 10.4% to 3.2%, respectively. As the bulk-phase properties studied are primarily governed by intermolecular interactions, this suggests that the OPLS-VSIL force field is correctly describing the local interactions occurring between the cations and anions.



Page 24 of 42

Table 10. Summarized Mean Absolute Errors (%) for Bulk-Phase Physical Properties of 1Alkyl-3-methylimidazolium [RMIM]-Based Ionic Liquids.

Density Heat of vaporization Heat capacity Viscosity Surface tension Self-diffusion coefficient (D+/D-)

0.8*OPLS-2009IL 2.8 6.6 10.6 10.4 15.2 24.4/26.6

OPLS- VSIL 3.1 3.4 9.5 3.2 7.3 9.6/21.1

Hydrogen Bonding. Along with improved accuracy in the reproduction of condensed phase properties, parameter development for both the cations and anions also emphasized the reproduction of radial distribution functions (RDF) derived from ab initio molecular dynamics (AIMD) simulations. Direct comparisons were made to reported AIMD calculations performed by Kirchner and coworkers for the ion pair combinations of [BMIM][TfO],51 [EMIM][DCA],52 and [EMIM][SCN].52 Of particular importance are intermolecular interactions involving the most acidic proton on the [RMIM] ring, i.e., H2-C2 bisecting the nitrogen atoms, as it has been shown to play a predominantly large role in dictating solvent properties.16,100,101 Problematically, many ionic liquid force fields underestimate the hydrogen-bonding ability of the imidazolium H2-C2 ring proton.26,102 It should be noted that some parameterization efforts have paid special attention to the H2-C2 hydrogen. For example, Liu et al. adjusted the van der Waals diameter of the acidic hydrogen to match gas-phase molecular structures computed using ab initio methods.41 Even so, a poor cation charge set103 or improper charge scaling31 can still have a large negative effect on the prediction of structural properties. To examine the OPLS-VSIL force field’s ability to reproduce local cation-anion intermolecular interactions, RDFs were calculated for the [RMIM] 24 ACS Paragon Plus Environment

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


ring hydrogen atoms and alkyl side chain hydrogen atoms with corresponding H-bond accepting atoms from the [TfO], [DCA], and [SCN] anions. Radial distribution functions were first calculated for the ionic liquid pair [BMIM][TfO] and compared to Kirchner’s simulations.51 In their work, the AIMD RDFs were computed using a periodic box of 32 ion pairs and employed the BLYP-D3/MOLOPT-DZVP-SR-GTH theory level. The OPLS-VSIL gave RDF intensities and interacting distances between the oxygen atoms of the [TfO] anion and the [BMIM] protons that were remarkably similar to the AIMD-derived RDFs (compare Figure 4 to “Figure 3” from Kirchner’s publication51). A prominent cation-anion interaction was found between the H2-C2 proton and the O atom of the [TfO] anion (Figure 4), as expected.51,104 While the height of the peak is determined by the number of neighbors, a tall peak with a short interacting distance is generally indicative of strong hydrogen bonding. The often-criticized poor hydrogen bonding description common to fixed charge ionic liquid force fields is improved here using the OPLS-VSIL. For example, a computed H2–O distance of 258.3 pm found with the 0.8*OPLS-2009IL39 was improved to 228.3 pm for OPLS-VSIL, which compared favorably to the 230.6 pm value from AIMD (Table 11). Even polarizable force fields can yield overestimated hydrogen bonding distances, e.g., a reported [EMIM][TfO] CHARMM simulation using Drude oscillators gave a H2–O peak height of 2.7 centered at a distance of 410 pm.31 See Table 11 for additional [BMIM][TfO] hydrogen bonding interaction comparisons.



Figure 4. Computed radial distribution functions between the hydrogen atoms of [BMIM] and oxygen atom of [TfO] using the OPLS-VSIL force field.


Page 26 of 42

Page 27 of 42 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 5. Radial distribution functions for the ionic liquids [EMIM][DCA] and [EMIM][SCN] between the ring H2 and methyl side-chain H10 hydrogen atoms of [EMIM] and the terminal 27 ACS Paragon Plus Environment


Page 28 of 42

nitrogen or sulfur atoms of [DCA] or [SCN] from (a) AIMD simulations reported by Weber and Kirchner (figure adapted from ref. 52) and (b) the OPLS-VSIL force field.

Table 11. Interaction Distances (pm) and g(r) From Radial Distribution Functions Computed Using OPLS-VSIL and Ab Initio MD (AIMD) for the Ionic Liquid [BMIM][TfO].a Distance (pm) VSIL AIMDb H2 – O H4-5 – O H6 – O H10 – O a

228.3 225.0 255.0 258.3

230.6 228.4 256.0 257.1

g(r) VSIL

AIMDb

2.9 2.3 1.5 1.7

2.9 2.2 1.6 1.7

See Fig. 4 for atom definitions. bRef. 51.

Radial distribution functions were also computed for the [EMIM][DCA] and [EMIM][SCN] ionic liquids using the OPLS-VSIL force field and compared to AIMD simulations reported by Weber and Kirchner.52 The AIMD calculations were comprised of 32 ion pairs in a periodic box and applied the same density functional theory (DFT) methodology as their [BMIM][TfO] work.51 The AIMD calculated RDFs are shown in Figure 5(a), which was adapted from “Figure 4” of Weber and Kirchner’s publication.52 Comparison to the OPLS-VSIL simulations in Figure 5(b) finds the virtual site force field’s results to be nearly indistinguishable from AIMD. The H2-C2 imidazolium ring proton gave large intensities when interacting with the terminal nitrogen atoms of the [DCA] and [SCN] anions. For example, OPLS-VSIL computed intensity peaks of 2.7 and 2.6 centered around 241 and 238 pm for [EMIM][SCN] and [EMIM][DCA], respectively. Both values compared well to the peak intensities of approximately 28 ACS Paragon Plus Environment

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


2.6-2.7 from Weber and Kirchner’s AIMD simulations. The OPLS-VSIL computed interactions between the methyl side-chain hydrogen atoms (H10-C10) and the cyano-based anions’ terminal nitrogen atoms also reproduced the AIMD shapes, intensities, and distances. Hydrogen bonding between the [EMIM] protons and the sulfur atom of [SCN] have been reported to have weaker intensities and longer distances as compared to the nitrogen atom.52,105 Accordingly, these H2–S and H10–S interactions are again correctly modeled using the virtual site force field (Figure 5). Ring–Ring Stacking. The ability to properly model p–p cation-cation arrangements and simulate the influence of cation-anion interactions upon p–p stacking was another important criterion in developing the OPLS-VSIL parameters. Combined distribution functions (CDFs) were used to analyze the interaction between two [RMIM] cations by monitoring their distance, d, and the angle, a, between their center-of-ring (CoR) throughout the MD trajectories. Figure 6 shows the angle created by a normalized vector (RN) perpendicular to the plane of a given imidazolium ring and the vector that connects two [RMIM] rings by their geometric centers as a function of distance. Angles near 0 and 180 degrees are indicative of p–p stacking. Focusing on the [BMIM][Cl] simulation finds ideal p–p stacking between the imidazolium rings with peak intensities located primarily at distances close to 400 pm and with angles of approximately 0 and 180 degrees. A typical equilibrium dCoR–CoR for two [BMIM] rings has been reported to be around 400 pm.106 In addition, a computed dCoR–CoR of approximately 650-800 pm at an angle of 90 degrees indicates that the second cation has migrated into the plane of the ring. A decrease in the population of stacking imidazolium rings represents increased competition with an anion for space above and below the center of the ring, such as for [BMIM][BF4] (Figure 6).107,108 Highlevel QM calculations, e.g., symmetry-adapted perturbation theory (SAPT) and DFT, of gas phase [RMIM][Cl] ion pairs have found that on-top cation-anion interactions (i.e. above/below 29 ACS Paragon Plus Environment


the ring) and in-plane hydrogen bonding are energetically equivalent.15,109 The anions can stabilize p–p stacking if they primarily lie in-plane ([Cl] or [NO3]), as pointed out previously by Matthews et al.110 Substitution with weakly coordinating anions, e.g., [TfO] or [BF4], leads to shift from p–p stacking to an alternating cation-anion conformation.111 Accordingly, the strong intensities localized around 400 pm for the [BMIM][Cl] solvent computed here using CDF have moved to a dCoR–CoR of 600 pm with angles between 0 and 45º for [EMIM][SCN] and [BMIM][BF4]. Spatial distribution functions (SDF) can provide further insight into the preferential positions of the ions surrounding the [RMIM] cation throughout the course of the MD trajectory. The TRAVIS program46 was used to construct and analyze the SDFs derived from the OPLSVSIL simulations of [BMIM][Cl], [BMIM][BF4], [BMIM][TfO], [EMIM][DCA], and [EMIM][SCN]. The SDF isosurface values for all of the cations and anions are provided in the Supporting Information Table S9. As expected, the anions preferred to occupy the area over the H2 atom and the H6 and H10 hydrogen atoms on the alkyl side chains, as well as at the H4 and H5 locations (Figures 7 and 4). The smaller Cl- anion favored direct interaction with the H2-C2 bond, often referred to as an “in-plane hydrogen bonding interaction”,15,112 whereas the larger anions, e.g., BF4- and TfO-, occupied the location over the [RMIM] ring, or an “on-top,” ion-ion interaction.107,108


Page 30 of 42

Page 31 of 42 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 6. Combined distribution functions (CDF) for the ionic liquids [BMIM][Cl], [BMIM][BF4], and [EMIM][SCN] from the OPLS-VSIL simulations. Illustrations are given of the plotted angle a versus the distance d for the center of the ring (CoR) interaction between two [RMIM]. The CDF plot is given with a relative intensity color for the CoR-CoR interactions.



Figure 7. Spatial distribution function of the anions and cations around the cation in the ionic liquid simulations using OPLS-VSIL, where blue equals the [RMIM] cations and the other colors represent the corresponding anions.

Conclusions In summary, new OPLS force field parameters have been developed for imidazoliumbased ionic liquids that yielded accurate dynamics, thermochemical properties, and correct local cation-anion intermolecular interactions. The OPLS-VSIL provides an empirically-derived set of partial charges for [RMIM] with a new topology that incorporates a virtual site atom bisecting 32 ACS Paragon Plus Environment

Page 32 of 42

Page 33 of 42 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 nitrogen atoms to offload negative charge to inside the plane of the ring. Guided by free energy of hydration calculations, the atom charges and Lennard-Jones terms were fine-tuned for both [RMIM] and 11 different anions to accurately compute bulk phase ionic liquid properties and radial distribution functions derived from ab initio molecular dynamics simulations. Overall, the OPLS-VSIL reproduced experimental densities, DHvap, heat capacities, viscosities, surface tensions, and self-diffusivities with a high degree of accuracy and properly replicated cationanion hydrogen bonding and p-p cation-cation stacking. 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.

Supporting Information Available: OPLS-AA bonded and nonbonded parameters for the [RMIM] cation and all anions discussed; heat capacity and heat of vaporization adjustment values; SDF isovalues; cubic lengths for equilibrated boxes; solvent experimental values, uncertainties, measurement methodologies, and references. The Supporting Information is available free of charge on the ACS Publications website at http://pubs.acs.org.

References (1) Hayes, R.; Warr, G. G.; Atkin, R. Structure and nanostructure in ionic liquids. Chem. Rev. 2015, 115, 6357–6426. (2) Moschovi, A. M.; Dracopoulos, V.; Nikolakis, V. Inter- and intramolecular interactions in imidazolium protic ionic liquids. J. Phys. Chem. B 2014, 118, 8673−8683.



(3) Castner, E. W.; Wishart, J. F.; Shirota, H. Intermolecular dynamics, interactions, and solvation in ionic liquids. Acc. Chem. Res. 2007, 40, 1217-1227. (4) 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. (5) Freemantle, M. An introduction to ionic liquids; RSC Publishing: Cambridge, UK, 2009. (6) Rogers, R. D.; Seddon, K. R. Ionic liquids--solvents of the future. Science 2003, 31, 792793. (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) Seddon, K. R. Ionic liquids for clean technology. J. Chem. Tech. Biotech. 1997, 68, 351-356. (10) Plechkova, N. V.; Seddon, K. R. Applications of ionic liquids in the chemical industry. Chem. Soc. Rev. 2008, 37, 123-150. (11) Welton, T. Room-temperature ionic liquids. Solvents for synthesis and catalysis. Chem. Rev. 1999, 99, 2071-2083. (12) Hallett, J. P.; Welton, T. Room-temperature ionic liquids: Solvents for synthesis and catalysis. 2. Chem. Rev. 2011, 111, 3508–3576. (13) Chen, S.; Zhang, S.; Liu, X.; Wang, J.; Wang, J.; Dong, K.; Sun, 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. (14) 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. (15) 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. (16) Kempter, V.; Kirchner, B. The role of hydrogen atoms in interactions involving imidazolium-based ionic liquids. J. Mol. Struct. 2010, 972, 22-34. (17) Dong, K.; Liu, X.; Dong, H.; Zhang, X.; Zhang, S. Multiscale studies on ionic liquids. Chem. Rev. 2017, 117, 6636–6695. (18) Izgorodina, E. I.; Seeger, Z. L.; Scarborough, D. L. A.; Tan, S. Y. S. Quantum chemical methods for the prediction of energetic, physical, and spectroscopic properties of ionic liquids. Chem. Rev. 2017, 117, 6696–6754. (19) Salanne, M. Simulations of room temperature ionic liquids: From polarizable to coarsegrained force fields. Phys. Chem. Chem. Phys. 2015, 17, 14270-14279. (20) 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. (21) 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. (22) 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.


Page 34 of 42

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


(23) 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. (24) Zahn, S.; Kirchner, B. Uncovering molecular secrets of ionic liquids. Chem. Modell. 2012, 9, 1-24. (25) 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. (26) 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. (27) 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. (28) Maginn, E. J. Molecular simulation of ionic liquids: Current status and future opportunities. J. Phys.: Condens. Matter 2009, 21, 373101. (29) 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. (30) 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. (31) Schroder, C. Comparing reduced partial charge models with polarizable simulations of ionic liquids. Phys. Chem. Chem. Phys. 2012, 14, 3089-3102. (32) KoBmann, S.; Thar, J.; Kirchner, B.; Hunt, P. A.; Welton, T. Cooperativity in ionic liquids. J. Chem. Phys. 2006, 124, 174506. (33) Kirchner, B.; Malberg, F.; Firaha, D. S.; Hollóczki, O. Ion pairing in ionic liquids. J. Phys. Condens. Matter 2015, 27, 463002. (34) López, P.; Méndez, F. Fukui function as a descriptor of the imidazolium protonated cation resonance hybrid structure. Org. Lett. 2004, 6, 1781-1783. (35) Arduengo, A. J. Looking for stable carbenes: The difficulty in starting anew. Acc. Chem. Res. 1999, 32, 913–921. (36) Hollóczki, O.; Firaha, D. S.; Friedrich, J.; Brehm, M.; Cybik, R.; Wild, M.; Stark, A.; Kirchner, B. Carbene formation in ionic liquids: Spontaneous, induced, or prohibited? J. Phys. Chem. B 2013, 117, 5898-5907. (37) Thomas, M.; Brehm, M.; Hollóczki, O.; Kirchner, B. How can a carbene be active in an ionic liquid? Chem. Eur. J. 2014, 20, 1622–1629. (38) 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. (39) Doherty, B.; Zhong, X.; Gathiaka, S.; Li, B.; Acevedo, O. Revisiting OPLS force field parameters for ionic liquid simulations. J. Chem. Theory Comput. 2017, 13, 6131-6145. (40) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. Modeling ionic liquids using a systematic all-atom force field. J. Phys. Chem. B 2004, 108, 2038-2047. (41) Liu, Z.; Huang, S.; Wang, W. A refined force field for molecular simulation of imidazolium-based ionic liquids. J. Phys. Chem. B 2004, 108, 12978-12989. (42) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. 35 ACS Paragon Plus Environment


(43) 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. (44) 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. (45) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (46) 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. (47) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. (48) 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. (49) Jorgensen, W. L.; Tirado-Rives, J. Potential energy functions for atomic-level simulations of water and organic and biomolecular systems. Proc. Nat. Acad. Sci. USA 2005, 102, 66656670. (50) Marin-Rimoldi, E.; Shah, 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. (51) Thomas, M.; Sancho Sanz, I.; Holloczki, O.; Kirchner, B. Ab initio molecular dynamics simulations of ionic liquids. In Nic symposium 2016; Binder, K., Muller, M., Kremer, M., Schnurpfeil, A., Eds. 2016; Vol. 48, p 117-124. (52) Weber, H.; Kirchner, B. Complex structural and dynamical interplay of cyano-based ionic liquids. J. Phys. Chem. B 2016, 120, 2471−2483. (53) Chaban, V. Polarizability versus mobility: Atomistic force field for ionic liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055–16062. (54) Bhargava, B. L.; Balasubramanian, S. Refined potential model for atomistic simulations of ionic liquid [bmim][pf6]. J. Chem. Phys. 2007, 127, 114510. (55) Budhathoki, S.; Shah, J. K.; Maginn, E. J. Molecular simulation study of the solubility, diffusivity and permselectivity of pure and binary mixtures of CO2 and CH4 in the ionic liquid 1-n-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide. Ind. Eng. Chem. Res. 2015, 54, 8821-8828. (56) Carvalhoa, P. J.; Álvarez, V. H.; Marruchoa, I. M.; Aznar, M.; Coutinho, J. A. P. High pressure phase behavior of carbon dioxide in 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide and 1-butyl-3-methylimidazolium dicyanamide ionic liquids. J. Supercrit. Fluids 2009, 50, 105-111. (57) Jacquemin, J.; Husson, P.; Majer, V.; Gomes, M. F. C. Influence of the cation on the solubility of CO2 and H2 in ionic liquids based on the bis(trifluoromethylsulfonyl)imide anion. J. Solution Chem. 2007, 36, 967-979. (58) Jorgensen, W. L.; Tirado-Rives, J. Molecular modeling of organic and biomolecular systems using BOSS and MCPRO. J. Comput. Chem. 2005, 26, 1689-1700. (59) Jensen, K. P.; Jorgensen, W. L. Halide, ammonium, and alkali metal ion parameters for modeling aqueous solutions. J. Chem. Theory Comput. 2006, 2, 1499-1509.


Page 36 of 42

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


(60) Jensen, K. P. Improved interaction potentials for charged residues in proteins. J. Phys. Chem. B 2008, 112, 1820–1827. (61) Pearson, R. G. Ionization potentials and electron affinities in aqueous solution. J. Am. Chem. Soc. 1986, 108, 6109-6114. (62) Florián, J.; Warshel, A. Langevin dipoles model for ab initio calculations of chemical processes in solution: Parametrization and application to hydration free energies of neutral and ionic solutes and conformational analysis in aqueous solution. J. Phys. Chem. B 1997, 101, 55835595. (63) Wasserscheid, P.; Welton, T. Ionic liquids in synthesis; Wiley-VCH: Germany, 2008; Vol. 1. (64) Reichardt, C.; Welton, T. Solvents and solvent effects in organic chemistry; 4th ed.; WileyVCH: Weinheim, Germany, 2011. (65) Borodin, O. Polarizable force field development and molecular dynamics simulations of ionic liquids. J. Phys. Chem. B 2009, 113, 11463-11478. (66) Schroder, C.; Steinhauser, O. Simulating polarizable molecular ionic liquids with drude oscillators. J. Chem. Phys. 2010, 133, 154511. (67) 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. (68) 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. (69) 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. (70) 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. (71) Armstrong, J. P.; Hurst, C.; Jones, R. G.; Licence, P.; Lovelock, K. R. J.; Satterley, C. J.; Villar-Garcia, I. J. Vapourisation of ionic liquids. Phys. Chem. Chem. Phys. 2007, 9, 982-990. (72) Strasser, D.; Goulay, F.; Kelkar, M. S.; Maginn, E. J.; Leone, S. R. Photoelectron spectrum of isolated ion-pairs in ionic liquid vapor. J. Phys. Chem. A 2007, 111, 3191-3195. (73) Leal, J. P.; Esperança, J. M. S. S.; Minas da Piedade, M. E.; Canongia Lopes, J. N.; Rebelo, L. P. N.; Seddon, K. R. The nature of ionic liquids in the gas phase. J. Phys. Chem. A 2007, 111, 6176-6182. (74) 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. (75) Esperanca, J. M. S. S.; Canongia Lopes, J. N.; Tariq, M.; Santos, L. M. N. B. F.; Magee, J. W.; Rebelo, L. P. N. Volatility of aprotic ionic liquids — a review. J. Chem. Eng. Data 2010, 55, 3-12. (76) Cervinka, C.; Pádua, A. A. H.; Fulem, M. Thermodynamic properties of selected homologous series of ionic liquids calculated using molecular dynamics. J. Phys. Chem. B 2016, 120, 2362-2371. (77) 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.



(78) 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-11805. (79) Pascal, T. A.; Lin, 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. (80) Waheed, Q.; Edholm, O. Quantum corrections to classical molecular dynamics simulations of water and ice. J. Chem. Theory Comput. 2011, 7, 2903–2909. (81) 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. (82) 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. (83) 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. (84) Hess, B. Determining the shear viscosity of model liquids from molecular dynamics simulations. J. Chem. Phys. 2002, 116, 209-217. (85) Alejandre, 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. (86) 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. (87) Lynden-Bell, R. M.; Del Pópolo, M. G. Simulation of the surface structure of butylmethylimidazolium ionic liquids. Phys. Chem. Chem. Phys. 2006, 8, 949-954. (88) 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. (89) 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. (90) 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. (91) 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. (92) 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. (93) Berendsen, H. J. C. Transport properties computed by linear response through weak coupling to a bath. In Computer simulation in materials science; NATO ASI Series (Series E: Applied Sciences); Meyer, M., Pontikis, V., Eds.; Springer: Dordrecht, 1991; Vol. 205. (94) Kalugin, O. N.; Riabchunova, A. V.; Voroshylova, I. V.; Chaban, V. V.; Marekha, B. A.; Koverga, V. A.; Idrissi, A. Transport properties and ion aggregation in mixtures of room temperature ionic liquids with aprotic dipolar solvents. In Modern problems of molecular physics; Bulavin, L., Chalyi, A., Eds.; Springer: Cham, 2018. 38 ACS Paragon Plus Environment

Page 38 of 42

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


(95) Fileti, E. E.; Chaban, V. V. The force field for imidazolium-based ionic liquids: Novel anions with polar residues. Chem. Phys. Lett. 2015, 633, 132-138. (96) Chaban, V. V.; Voroshylova, I. V.; Kalugin, O. N. A new force field model for the simulation of transport properties of imidazolium-based ionic liquids. Phys. Chem. Chem. Phys. 2011, 13, 7910-7920. (97) Köddermann, T.; Paschek, D.; Ludwig, R. Molecular dynamic simulations of ionic liquids: A reliable description of structure, thermodynamics and dynamics. ChemPhysChem 2007, 8, 2464–2470. (98) Cooke, B.; Schmidler, S. C. Preserving the boltzmann ensemble in replica-exchange molecular dynamics. J. Chem. Phys. 2008, 129, 164112. (99) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635-2643. (100) 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. (101) 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. (102) 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. (103) 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. (104) Gehrke, S.; Domaros, M. v.; Clark, R.; Hollóczki, O.; Brehm, M.; Welton, T.; Luzar, A.; Kirchner, B. Structure and lifetimes in ionic liquids and their mixtures. Faraday Discuss. 2018, 206, 219-245. (105) Batista, M. L. S.; Kurnia, K. A.; Pinho, S. P.; Gomes, J. R. B.; Coutinho, J. A. P. Computational and experimental study of the behavior of cyano-based ionic liquids in aqueous solution. J. Phys. Chem. B 2015, 119, 1567–1578. (106) 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. (107) Marekha, B. A.; Kaluginb, O. N.; Idrissi, A. Non-covalent interactions in ionic liquid ion pairs and ion pair dimers: A quantum chemical calculation analysis. Phys. Chem. Chem. Phys. 2015, 17, 16846-16857. (108) 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. (109) Hunt, P. A.; Gould, I. R. Structural characterization of the 1-butyl-3-methylimidazolium chloride ion pair using ab initio methods. J. Phys. Chem. B 2006, 110, 2269–2282. (110) 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.



(111) Matthews, R. P.; Villar-Garcia, I. J.; Weber, C. C.; Griffith, J.; Cameron, F.; Hallett, J. P.; Hunt, P. A.; Welton, T. A structural investigation of ionic liquid mixtures. Phys. Chem. Chem. Phys. 2016, 18, 8608-8624. (112) Del Pópolo, M. G.; 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.


Page 40 of 42

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


TOC graphic



TOC graphic 82x29mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 42 of 42

A Virtual Site OPLS Force Field for Imidazolium ... - ACS Publications

Recommend Documents