First-Principles United Atom Force Field for the ... - ACS Publications

Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI

Article

First Principles United Atom Force Field for the Ionic Liquid [BMIM][BF] : An Alternative to Charge Scaling 4

Chang Yun Son, Jesse Gatten McDaniel, Jordan R. Schmidt, Qiang Cui, and Arun Yethiraj J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b12371 • Publication Date (Web): 22 Feb 2016 Downloaded from http://pubs.acs.org on February 27, 2016

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

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

Page 1 of 30

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

First Principles United Atom Force Field for the Ionic Liquid BMIM+BF− 4 : An Alternative to Charge Scaling Chang Yun Son, Jesse G. McDaniel, J. R. Schmidt, Qiang Cui, and Arun Yethiraj∗ Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin-Madison, 1101 University Avenue, Madison, WI 53706, USA E-mail: [email protected]

Abstract Molecular dynamics study of ionic liquids (IL) is a challenging task. While accurate fully polarizable atomistic models exist, they are computationally too demanding for routine use. Most non-polarizable atomistic models predict diffusion constants that are much lower than experiment. Scaled charge atomistic models are cost-effective and give good results for single component ILs, but are in qualitative error for the phase behavior of mixtures, due to inaccurate prediction of the IL cohesive energy. In this work, we present an alternative approach for developing computationally efficient models that importantly preserves both the correct dynamics and cohesive energy of the IL. Employing a “top-down” approach, a hierarchy of coarse-grained models for BMIM+ BF− 4 are developed by systematically varying the polarization/atomic resolution of the distinct functional groups. Parameterization is based on symmetry-adapted perturbation theory (SAPT) calculations involving the homo-molecular species; all

1

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

cross interactions are obtained from mixing rules and there are no adjustable parameters. We find that enhanced dynamics from an united-atom description counteracts the effect of reduced polarization, enabling computationally-efficient models that exhibit quantitative agreement with experiment for both static and dynamic properties. We give explicit suggestions for reduced-description models that are computationally more efficient, more accurate, and more fundamentally sound than existing non-polarizable atomistic models.

2


Page 2 of 30

Page 3 of 30

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

Introduction

Ionic liquids (IL) have many promising applications as solvents due to their low vapor pressure and formation of both polar and hydrophobic domains, leading to their ability to solvate both polar and non-polar molecules. In addition, their ionic nature leads to further applications as conductive electrolytes. There exist many different types of ILs, originating from different choices of the organic cations and anions, suggesting the possibility to specifically choose and tailor an IL for a particular solvent application. The success of such tailoring, however, requires a detailed understanding of properties of the ILs, as well as all potential IL/solute interactions; such an understanding can be advanced with the aid of molecular dynamics (MD) simulations. In this work, we develop a hierarchy of models of varying complexity and computational efficiency for use in MD simulations. There have been numerous MD simulation studies on ILs, as well as parallel work on model development for these systems, and we refer the reader to two comprehensive review articles on this subject. 1,2 It has become apparent that the standard trade-offs between accuracy, transferability, and computational expense inherent to all models, are especially relevant for the simulation of ILs. For example, it is generally recognized that explicit polarization of the constituent ions plays an important role in determining both the dynamics and heat of vaporization of the liquid. 3 However, due to the high cohesive energy of ILs, their dynamics at room temperature is very slow, requiring long simulations for accurate equilibration and statistical sampling. The computational expense of these simulations becomes even more problematic when attempting to study phase behavior of mixtures involving ILs. It is therefore of practical use to have a collection of IL models, encompassing different levels of accuracy and sophistication as well as computational expense. There have been much effort 3–5 dedicated to developing explicitly polarizable models that target high accuracy and correct physics in order to predict macroscopic liquid properties. Other groups have worked to develop more simplistic models that attempt to achieve good accuracy at a much lower computational cost, which would allow for the practical use of these models in a wider 3



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

range of applications that demand more extensive simulations. 6,7 Care must be taken in the latter approach, however, as such simplistic models may sacrifice important physics in the description of intra- and inter-molecular interactions, which may limit the transferability of the model beyond the study of the pure ILs. Some of the most successful examples of computationally efficient IL models involve nonpolarizable scaled or reduced charge models, in which the (monovalent) cations and anions possess absolute charges less than unity. It has been empirically demonstrated that such models exhibit enhanced dynamics and lower heats of vaporization, in both cases bringing the predicted properties in closer agreement with experiment. The apparent success of these models in predicting bulk IL properties has now led to their routine use in simulations. 2 However, in their recent work on the phase behavior of BMIM+ BF− 4 and poly(ethylene oxide) (PEO) mixtures, Choi and Yethiraj 8 found that scaled charge models for BMIM+ BF− 4 failed to predict phase separation of PEO in BMIM+ BF− 4 at any temperature, in qualitative disagreement with experiment. 9 We speculate that this failure is not system specific, but rather that primarily due to their significant and artificial underestimation of the IL cohesive energy (vide infra), scaled charge models may generally perform poorly for predicting phase behavior in IL mixtures. It therefore seems necessary to develop alternative approaches to charge scaling for deriving computationally efficient IL models, especially for use in the simulation of the phase behavior of IL-based mixtures; this is the motivation for the present work. In the present work, we employ a top-down approach to develop a hierarchy of reducedcomplexity IL models, focusing on the commonly studied IL BMIM+ BF− 4 . Starting with 4 an accurate, explicitly polarizable, all-atom model for BMIM+ BF− 4 , we systematically

coarse-grain the polarization and all-atom resolution, resulting in several different models employing different amounts of polarization and united-atom character. Following Choi et al., 4 the inter-molecular interactions in these models are parameterized entirely based on ab initio, symmetry-adapted perturbation theory calculations. The end result of our study is the

4


Page 4 of 30

Page 5 of 30

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


development of several different BMIM+ BF− 4 models that are potentially useful in different applications, due to their different compromises in accuracy and computational expense. We choose our specific coarse-graining schemes motivated by the potential application of such models to study phase-equilibria in ionic-liquid based mixtures, for which we desire to maintain accuracy in the prediction of the cohesive energy, dynamics, structure, and density (vide infra), but at reduced computational expense. The rest of the paper is organized as follows: We first present a discussion of the ideal attributes of a reduced-complexity IL model, from which we motivate our specific coarsegraining scheme. We then discuss the methodology of our approach, and present predicted bulk properties for six different resolution BMIM+ BF− 4 models that we develop. Finally, we highlight both the advantages and disadvantages of each model, and applications for which each model is potentially useful.

2

Coarse-graining philosophy

There is an inevitable tradeoff between computational cost and physical accuracy inherent to any coarse grained model development; the challenge lies in determining the essential physics to include in the model for a specific application, while omitting non-essential details to enhance computational efficiency. For ILs, the most fundamental physical interactions are the electrostatic attractions and repulsions between ions. In fact, it has been shown that in the bulk liquid, short-range exchange repulsion interactions approximately cancel the long-range van-der-Waals dispersion attraction, so that the total energy mirrors the electrostatic energy-both of which are very large and net attractive. 4 Polarization, while playing an important role in speeding up the dynamics, and contributing to differences in ion-pair energies in the gas and liquid phase (enthalpy of vaporization), only contributes a few percent to the cohesive energy of the liquid. 4 Therefore, altering the electrostatic cohesive energy correspondingly alters the total cohesive energy of the liquid by the roughly

5



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 same extent. The cohesive energy of ILs is one of the fundamental determinants of their properties. It gives rise to high heats of vaporization and low vapor pressure, as well as high viscosity and slow, glassy dynamics. Reducing this cohesive energy therefore will dramatically change these properties, speeding up dynamics and lowering heats of vaporization; this is exactly what is found in scaled/reduced charge models of ILs. 10 However, there is no physical basis for such reduced electrostatics for neighboring ion pairs; both polarization and charge transfer effects are strictly attractive interactions 11 , and can only increase the magnitude of the cohesive energy, not lower it. As mentioned previously, the finding that reduced charge models are unable to predict phase separation in BMIM+ BF− 4 -PEO mixtures can be attributed to the reduced cohesive energy of the IL phase. We note the following important caveat in our discussion of charge scaling–our discussion entirely pertains to pure ionic solutions such as ILs and not dilute aqueous electrolyte solutions. For the latter systems, charge scaling of (dilute) ions may effectively compensate for missing polarization in non-polarizable water models, 12–14 as motivated by electronic continuum models. 15,16 However, such cases should be kept distinct, due to the important difference that “solvent” interactions are dramatically altered in the former but not the latter case. Strong and local polarization in the neighboring ions are distinctive in ILs, 17 and electronic continuum models may not be appropriate as in the case of electrolytes near their solubility limit. 14 The central goal of this work is therefore to develop a BMIM+ BF− 4 model of greater computational efficiency that maintains the same dynamics and cohesive energy as our previous model. 4 To achieve this, we attempt to both reduce the extent of polarization in our model as well as employ united atom sites to absorb explicit hydrogen atoms. We anticipate that such effects may offset each other, with the former lowering and the later enhancing the dynamics of the model while maintaining a relatively constant cohesive energy of the liquid. We rely on previous model development, based on symmetry-adapted perturbation

6


Page 6 of 30

Page 7 of 30


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

evaluation of the importance of the level of detail in the description of each group. The functional group force fields are developed through a “top-down” approach, starting with an explicitly polarizable, all-atom model, and systematically reducing the resolution through reduction/omission of polarization, removal of explicit hydrogen atoms (UA-description), and coarse graining the BF− 4 anion to a single-site model (CG). We note that parameters of the coarser models (UA or CG) are obtained directly from ab initio calculations. We expect that the accuracy of the UA-description may depend on the functional group; the alkyl hydrogens may have minimal impact on the system properties, while the hydrogen atoms of the imidazolium ring have been suggested to form important hydrogen-bond like interactions with anions or solutes. 18 Atomistic simulations of PEO / BMIM+ BF− 4 mixtures also suggested the important role of specific interactions between the ring hydrogens and the oxygen atoms of the polymer. 8 Thus it is interesting to separately test the robustness of these various united atom cation models. The approximate treatment of the tetrahedral anion as a single spherical bead is also interesting to evaluate, as for imidazolium-based ILs, the influence of different anions is generally considered to be dependent on size, rather than structure. On the other hand, the atomistic detail of the anion may be important to describe the specific interaction with the imidazolium ring. By varying the resolution of each functional groups independently, the importance of the above-mentioned effects are systematically investigated. We investigate six models with increasing levels of coarse-graining, and these are shown in Figure 1 and listed in Table 1. The first, fully atomistic and polarizable model is taken from recently developed model (named as AA AP, All Polarizable) of Choi et al. 4 The second, AA RP (Ring Polarizable) model, has the same atomic resolution as the AA AP model but does not have drude particles on alkyl chain carbons. The polarizability of the BF− 4 anion is described with a single drude (SD) particle on the central boron instead of drude particles on each of the 5 atoms. The next model (ringH) removes hydrogen atoms from alkyl groups of cation but preserves the hydrogen atoms on imidazolium ring carbons. Two versions of

8


Page 8 of 30

Page 9 of 30

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


this model are tested, with explicit polarization on imidazolium ring and anion (ringH pol model) and without any drude particles with same atomic description (ringH nopol model). In the non-polarizable united atom models (UA AA and UA CG) the hydrogen atoms on the ring carbons are further removed; UA AA keeps fluorine atoms on the anion while UA CG model uses a single bead for the anion. Table 1: The resolution of each functional groups in tested models of BMIM+ BF− 4 . AA: all atom, UA: united atom(hydrogen is removed), CG:coarse grained(single bead anion), SD:single drude particle, pol :drude particles on each non-hydrogen atoms,none:no drude particles Model abbr. AA AP AA RP ringH pol ringH nopol UA AA UA CG

3.2

cation ring atoms drude AA pol AA pol AA pol AA none UA none UA none

cation atoms AA AA UA UA UA UA

alkyl drude pol none none none none none

anion atoms drude AA pol AA SD AA SD AA none AA none CG none

Force field parameterization

The force field parameters are optimized to have good agreement with the SAPT calculated homo-molecular interaction energy in four individual energy components - electrostatics, exchange-repulsion, induction, and dispersion interactions. For the details on the physical basis of functional form and parameterization scheme of the force field, we refer the reader to previous papers 4,19 and supporting information. Three sets of parameters are developed in this work - united atom parameters for butane, united atom parameters for the imidazolium ring and coarse-grained single bead parameters for BF− 4 . An exception to the functional group-based parameterization was fitting point charges for the BMIM+ cation, in which two new sets of atomic charges were developed for the united atom models(ringH and UA models) by explicitly fitting the ab initio electric field calculated at many points on different 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

van der Waals surfaces of the BMIM+ cation (the single bead anion in UA CG model has full charge of -1.0). For the UA models, the polarization is identical to the all atom model 4 since the all-atom model had drude particles only on non-hydrogen atoms. The new polarizability for the single drude(SD) BF− 4 model as well as dispersion parameters for the UA models were parameterized by fitting coupled-perturbed Kohn-Sham (CPKS) 20 linear response calculations of the functional groups. Short-range, charge penetration terms are fit to DFT-SAPT 21 calculations of homonuclear dimers, namely, imidazolium dimers, butane dimers and BF− 4 dimers. We note that no additional computational cost was required for the present model development, as all electronic structure calculations utilized for fitting were taken from our previous work. 4

3.3

Simulation details

Molecular dynamics simulations are carried out using GROMACS version 4.6.5. 22 Initially, 200 BMIM+ BF− 4 ion pairs are randomly placed in a large box and then were energy minimized. The system is equilibrated for 1 ns at three different temperatures: 300, 350 and 400 K, and at 1 bar with a Nosé-Hoover 23,24 thermostat and Berendsen 25 barostat. After the equilibration, 5 ns NPT simulation is performed with Parrinello-Rahman barostat to ensure the correct volume fluctuation in NPT ensemble. For the calculation of dynamic properties, a 20 ns final production run is carried out in NVT ensemble using Nosé-Hoover thermostat. The van der Waals interactions are shifted to zero at 1.4 nm. For long-range electrostatic interactions, the Particle Mesh Ewald(PME) method 26,27 with a spacing of 0.12 nm and a real space cutoff distance 1.4 nm is used. To incorporate explicit polarization effect, the positions of the shell particles are optimized each time step until the root-mean-squared force on the shells becomes less than 0.1kJ mol−1 nm−1 . The charges of the Drude oscillators are given √ by q = αk, with a spring constant k of 0.1 atomic units and atomic polarizability α. Thole screening function 28 is used for intra-molecular shell-shell interactions with a Thole parameter of 2.0. Bonded potentials are adopted from the OPLS-AA model of Sambasivarao and 10


Page 10 of 30

Page 11 of 30

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


Acevedo 29 with two exceptions for computational convenience: bond stretching potential parameters for boron and fluorine are taken from De Andrade and coworkers’ model 30 instead of using rigid bonds, and intramolecular 1-4 nonbonded interactions are ignored instead of scaling them down to 50%.

3.4

Calculation of thermodynamic properties

We define the cohesive energy, Ecohes , as the total intermolecular potential energy per ion pair. Although cohesive energy is normally assumed equivalent to the enthalpy of vaporization for neutral molecules, this is not true for an ionic liquids, since the enthalpy of vaporization is lowered by the strong attraction energy of gas-phase ion pairs 4 (see eq 1). We obtain this quantity by re-calculating the total potential energy of the system from the production NVT simulation trajectories with all bonded and intra-molecular non-bonded interactions turned off. Other thermodynamic and dynamic properties are calculated as follows. The enthalpy of vaporization is calculated from gas liquid ∆Hvap = Epotential − Epotential + RT

(1)

gas where Epotential is the total potential energy. Note that Epotential reflects the ion pair binding

energy in the gas phase. The self diffusion coefficient, D, is calculated using Einstein relation

2 1 ∂ |rcom,i (t)| D = lim t→∞ 6 ∂t

(2)

where |rcom,i (t)| is the position of the center-of-mass of particle i, and the average is over all initial times and particles of the same species (cation or anion). The dielectric constant, ǫ, is calculated from the fluctuation of rotational dipole moment

ǫS = 1 +

2 1 MD − hMD i2 , 3ǫ0 V kB T 11


(3)


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 30

where MD =

XX i

α

qi,α ri,α −

X

qi rcom,i

(4)

i

where qi,α and ri,α are the charge and position of particle i of species α. We further use the approximation that for a system of pure IL, the second term of Eq. 3 vanishes. For the simulation of scaled charge models, it has been suggested to reverse the scaling factor when calculating the dielectric constant and dipole moment of the system. 16 To compare this corrected calculation with Eq.3 we also calculated the corrected dielectric constant and dipole moment by ǫcorr = ǫel · ǫM D

µcorr =

√

ǫel · µM D

(5)

, respectively. The electronic dielectric constant ǫel is related to the charge scaling factor Qscale with ǫel = 1/Q2scale . For OPLS-AAQ084 model of which the atomic charges of OPLSAA model 29 is scaled with scaling factor of Qscale = 0.84, 31 ǫel = 1.42.

4

Results and Discussion

4.1

Cation-anion interaction energy

The energy of an ion pair is a good test of whether a reduced model can describe the properties in the gas phase. Of course, a reproduction of the ion pair energies is not equivalent to the performance of the model in the liquid state because non-pairwise-additive contributions may play important role. Nevertheless, a good description of pairwise interaction energies are essential for the transferability of the model both in gas and liquid phases. Figures 2 and 3 compare the total and exchange energy for a BMIM+ BF− 4 ion pair to the DFT-SAPT energy 4 for 1300 configurations. We emphasize that no parameters are adjusted in this comparison; combination rules are used on the homo-molecular interactions. For comparison results from the OPLS-AA 29 and OPLS-AAQ084 29,31 models are also shown. For the total energy, all of the models are asymptotically in semi-quantitative agreement 12


Page 13 of 30

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


with the ab initio calculations, except for the scaled charge (OPLS-AAQ084) model. Also, the inclusion of polarization is crucial for an accurate description of the strong attraction between ions at short separations. In Figure 2, first three models with some extent of polarization (AA AP, AA RP, ringH pol) perform equally well and closely match the DFTSAPT calculated energies, while the non-polarizable models (ringH nopol, UA AA, UA CG) and the OPLS-AA/OPLS-AAQ084 models underestimate the attraction. This clearly shows that the effect of polarization is strictly attractive and scaling the charge cannot capture the correct physics in highly charged media. For the non-polarizable models, we do not refit the force field to compensate for the residual energy difference due to neglect of polarization; this is intentional in our force field design in order to isolate the impact of polarization. Figure 2 also shows the effect of atomic description; the coarser the models, the more scattered the energies. The scattering in energies arises primarily from the short ranged exchange interaction (Figure 3), which is highly anisotropic and cannot be accurately reproduced with a spherical isotropic potential on an united atom site. Because the ion pair configurations are biased towards the imidazolium ring-anion (and not alkyl tail-anion) interactions, the UA resolution of the ring introduces more error than the UA resolution of the alkyl chain. The UA CG model suffers from this the most, suggesting that the single-site coarse graining of BF− 4 might not be ideal.

4.2

Static thermodynamic properties

All the models predict the liquid density quite accurately when compared to experiment (see Table 2), within the range of experimental values, except for the OPLSAA-Q084 model. The two ringH models predict almost identical densities when compared to the AA AP model, implying that retaining the ring hydrogen atoms is more important than retaining the alkyl hydrogen atoms, for the prediction of the density. The deviation of the UA AA and UA CG models’ densities from that of AA AP model highlights the enhanced description 13




Page 14 of 30

Page 15 of 30




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 30

Table 2: Comparison of thermodynamic properties from models of different resolution.

exp. AA AP AA RP ringH pol ringH nopol UA AA UA CG OPLS-AA OPLS-AAQ084

ρ (kg/m3 ) 300K 400K 1202 32 1130 32 1171 33 1180 1107 1182 1107 1187 1113 1185 1111 1202 1132 1172 1104 1177 1103 1138 1053

∆Hvap (kJ/mol) 300K 400K

Ecohes (kJ/mol) 300K 400K

128 34

-

-

-

128 128 135 159 151 158 184 144

118 120 127 150 143 150 171 132

-491 -489 -485 -481 -481 -467 -494 -360

-477 -476 -473 -469 -469 -457 -479 -345

The calculated liquid phase cohesive energy shows striking consistency throughout all tested models except the coarsest UA CG model, where the deviation from AA AP model is 5%. Since the enthalpy of mixing depends strongly on the cohesive energy, we expect that nonpolarizable models may be suitable for simulations of phase equilibria of large systems such as polymers in ILs. As expected, the OPLS-AA model predicts the density and cohesive energy of liquid phase relatively well, but predicts too high an enthalpy of vaporization, due to neglect of polarization in the gas phase. On the other hand, the scaled charge OPLSAAQ084 model predicts the enthalpy of vaporization more accurately than the full charge OPLS-AA model (∼15 kJ/mol higher than polarizable models) but predicts a liquid density that is too low, due to the significantly underestimated cohesive energy. It is interesting that full charge OPLS-AA model shows almost equivalent cohesive energy to AA AP model, and this explains why OPLS-AA model is able to simulate the correct phase behavior of PEO/IL mixture but the scaled charge model fails. 8 Observed failure of scaled charge models in describing the energetics of ILs is general; because of the reduced electrostatics, they will always underestimate the net IL attraction both in gas phase and in condensed phase. To demonstrate this, we computed the gas-phase interaction energy, condensed phase cohesive energy, density, and enthalpy of vaporization of three additional scaled charge models. 31,35,36 Figure S7 and Table S4 in the supporting 16


Page 17 of 30




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 30

model resolution for anions and cations are different. A similar trend is seen going from the ringH pol to the ringhH nopol to the UA AA model. Interestingly, the UA CG model has lower diffusion constants than the UA AA model, which again indicates the limitation of modeling the anion as a single sphere. The qualitative trend, i.e., removing polarization results in slower dynamics and removing hydrogen atoms results in faster dynamics, also holds for T=300K (see table 3). There are, however, quantitative differences, e.g., in the AA RP model D+ /D− < 1 at 400K D+ /D− > 1 at 300K. We speculate that this temperature dependent enhancement of dynamics in the united atom description can be explained with the change in diffusion mechanism of BMIM+ BF− 4 at different temperatures. At low temperature, each imidazolium ring is trapped in the surrounding counter ions by strong electrostatic interactions and may diffuses through a hopping-like mechanism; 37 thus the united atom description on alkyl groups does not enhance the dynamics significantly. On the other hand, there are more free ions at higher temperature 38 thus united atom alkyl group can enhance the dynamics of these free cation effectively. Table 3: Comparison of dynamic properties from different models of BMIM+ BF− 4 . Properties at 300K suffers from poor statistics and has higher uncertainty due to the slow dynamics of ILs. * Directly calculated from MD using Eq.3

exp. AA AP AA RP ringH pol ringH nopol UA AA UA CG OPLS-AA OPLS-AAQ084

D+ ( 10−11 m2 s−1 ) 300K 400K

D− (10−11 m2 s−1 ) 300K 400K

1.6 32

27.5 32

1.5 32

31.6 32

1.8 0.9 1.0 0.8 3.2 1.7 0.1 1.3

26.6 22.7 29.6 19.1 35.2 24.4 7.7 32.2

1.0 0.5 0.6 0.3 1.6 1.2 0.04 0.8

28.7 19.9 23.6 14.0 22.7 16.3 5.9 25.1

18


300K 14.6 39 11.7 40 11.0 ± 1 12.4 9.2 7.4 7.7 6.7 3.9 3.2*

ǫS 350K

400K

-

-

8.7 8.0 7.6 6.8 7.3 7.0 4.4 2.8*

8.2 7.3 6.8 6.3 6.3 6.0 3.5 2.6*

Page 19 of 30




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

charge model has an order of magnitude higher translational diffusion constant compared to the full charge OPLS-AA model, but its static dielectric constant is slightly lower. If one uses correction factor of Eq. 5 as suggested in the MDEC model, 16 calculated values of OPLS-AAQ084 model is very close to that of OPLS-AA model (ǫcorr = 3.7, µcorr = 3.3 at 400K). The inability of either OPLS-AA model (reduced or full charge) to correctly predict the dielectric response is due to the incorrect description of the dipole moment on the butyl group of the cation. This is illustrated in Figure 5b; for the AA AP model, removing the alkyl group charges reduces the average dipole moment significantly, while the trend is opposite for the OPLS models. We verified that this effect was not due to ensemble differences, as conclusions were similar when trajectories were switched and dielectric recalculated.

4.4

Structural analysis

Comparison of the spatial distribution function of anions around the cation as predicted by AA AP model and UA CG model sheds lights on the existence of hydrogen bonding in BMIM+ BF− 4 . As shown in Figure 6, there is strong association between the anion and the C1 carbon in all of the cation models. For the UA CG model the distribution of both anion and the alkyl group is more diffuse, while AA AP model shows strongly preferred orientations of both the anion and the alkyl tail. It is interesting that both with or without the hydrogen on the imidazolium ring, anions preferred to locate close to the carbon atoms of the ring. In the literature, there has been significant discussion about whether the C1-H hydrogen exhibits a hydrogen bond-like interaction with highly electronegative atoms in other groups. The BF− 4 association is not very directional (i.e. close to linear C1-H-anion angle) and only slightly stronger than the C2-H-anion interaction (as predicted by the models with ring hydrogens). Indeed, there is greater population of linearly associated anion to the C1 carbon in UA CG model, which is counter intuitive if the hydrogen bonding exists. The fact that this orientational organization of anions is observed from united atom description of the 20


Page 20 of 30

Page 21 of 30





Page 22 of 30

Page 23 of 30

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


shown by the prediction of the translational diffusion constant and static dielectric constant of the reduced charge OPLS-AAQ084 model. In this work, we have outlined a new approach to obtain cost effective and reliable united atom model based on first-principle SAPT calculations. A hierarchy of increasingly coarsegrained models for BMIM+ BF− 4 are explicitly parameterized to quantum calculations, i.e., there are no adjustable parameters. The predicted properties of each model are compared to experiments (where possible) and a fully polarizable all atom model. All the six developed models (to varying degrees) accurately predict the liquid density, cohesive energy, structure, and dynamics. Inherently, the non-polarizable models show deficiencies in the enthalpy of vaporization, which we attribute to the failure of these models to accurately describe the vapor state.

Remarkable equivalence of the cohesive energy in all

developed models indicate that the polarization effect in ILs is largely local and short-ranged; thus invalidating electronic continuum models in ILs. Interestingly, removal of polarizability and removal of explicit hydrogen atoms have contrasting effects on the ion diffusion constants, allowing for a useful cancellation of errors in the development of united atom models. The process of hierarchical coarse-graining provides physical insight of its own. For example, we find a strong correlation between the anion and the C1 hydrogen of the BMIM+ cation, as is inferred from experiments, and is usually explained in terms of hydrogen bonding. This correlation, however, also exists when the relevant hydrogen atoms are removed, thus emphasizing that the correlation is purely electrostatic in origin. We also found that coarse-graining the tetrahedral anion to a single spherical particle was rather drastic approximation, inducing large errors in the description of short range repulsion interaction, which is highly anisotropic in nature. This is a surprising result because it seems intuitively appealing to replace a symmetric BF− 4 anion with a single sphere. We have several recommendations regarding a choice of model for computer simulations of this IL. If the system is expected to be heterogeneous, as happens with interfaces, we suggest a polarizable model because polarizability is important for the description of a low

23



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

density phase. In this case we suggest the ringH pol model where the hydrogen atoms and polarizability in the cation ring are retained, but both are removed in the alkyl tail and the anion is atomistic and has a single drude particle. For the liquid phase we suggest the UA AA model, where the cation groups are treated using a united atom description and the anion is atomistic. This model provides a good description of the static and dynamic properties while also being computationally convenient, about 38 times faster than the fully atomistic polarizable model. We see no reason to recommend the OPLS-AA or scaled charge OPLS-AAQ084 models which are significantly less accurate than the UA AA model and also twice as expensive computationally.

6

Acknowledgement

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1111835. This work was partially supported by Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy, under Award DE-FG02-09ER16059. Computational resources were provided by the Center for High Throughput Computing at the University of Wisconsin.

7

Supporting information

Definition of atomic sites in different models. Ab initio electric field fits for ringH and UA cation models. Functional form of intermolecular force-fields and the complete force field parameters for each models. DFT-SAPT force-field fits for homo-dimer UA models (butane, BF− 4 , imidazolium dimer). Comparison of three additional scaled charge models to AA AP model.

24


Page 24 of 30

Page 25 of 30

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


References (1) Maginn, E. J. Molecular Simulation of Ionic Liquids: Current Status and Future Opportunities. J. Phys. Condens. Matter 2009, 21, 373101–373117. (2) 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. (3) Salanne, M. Simulations of Room Temperature Ionic Liquids: From Polarizable to Coarse-Grained Force Fields. Phys. Chem. Chem. Phys. 2015, 17, 14270–14279. (4) 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. (5) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463–11478. (6) Urahata, S. M.; Ribeiro, M. C. C. Structure of Ionic Liquids of 1-Alkyl-3methylimidazolium Cations: A Systematic Computer Simulation Study. J. Chem. Phys. 2004, 120, 1855–1863. (7) Roy, D.; Maroncelli, M. An Improved Four-Site Ionic Liquid Model. J. Phys. Chem. B 2010, 114, 12629–12631. (8) Choi, E.; Yethiraj, A. Entropic Mechanism for the Lower Critical Solution Temperature of Poly(ethylene oxide) in a Room Temperature Ionic Liquid. ACS Macro Lett. 2015, 4, 799–803. (9) Lee, H. N.; Newell, N.; Bai, Z. F.; Lodge, T. P. Unusual Lower Critical Solution Temperature Phase Behavior of Poly(ethylene oxide) in Ionic Liquids. Macromolecules 2012, 45, 3627–3633. 25



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

(10) Chaban, V. Polarizability versus Mobility: Atomistic Force Field for Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055–16062. (11) Szefczyk, B.; Sokalski, W. A. Physical Nature of Intermolecular Interactions in [BMIM][PF6 ] Ionic Liquid. J. Phys. Chem. B 2014, 118, 2147–2156. (12) Pegado, L.; Marsalek, O.; Jungwirth, P.; Wernersson, E. Solvation and Ion-Pairing Properties of the Aqueous Sulfate Anion: Explicit versus Effective Electronic Polarization. Phys. Chem. Chem. Phys. 2012, 14, 10248–10257. (13) Mason, P. E.; Wernersson, E.; Jungwirth, P. Accurate Description of Aqueous Carbonate Ions: An Effective Polarization Model Verified by Neutron Scattering. J. Phys. Chem. B 2012, 116, 8145–8153. (14) Kann, Z. R.; Skinner, J. L. A Scaled-Ionic-Charge Simulation Model that Reproduces Enhanced and Suppressed Water Diffusion in Aqueous Salt Solutions. J. Chem. Phys. 2014, 141, 104507–104513. (15) Leontyev, I. V.; Vener, M. V.; Rostov, I. V.; Basilevsky, M. V.; Newton, M. D. Continuum Level Treatment of Electronic Polarization in the Framework of Molecular Simulations of Solvation Effects. J. Chem. Phys. 2003, 119, 8024–8037. (16) Leontyev, I.; Stuchebrukhov, A. Accounting for Electronic Polarization in NonPolarizable Force Fields. Phys. Chem. Chem. Phys. 2011, 13, 2613–2626. (17) Schröder, C. Comparing Reduced Partial Charge Models with Polarizable Simulations of Ionic Liquids. Phys. Chem. Chem. Phys. 2012, 14, 3089–3102. (18) Liu, Z.; Wu, X.; Wang, W. A Novel United-Atom Force Field for Imidazolium-Based Ionic Liquids. Phys. Chem. Chem. Phys. 2006, 8, 1096–1104. (19) McDaniel, J. G.; Schmidt, J. R. Physically-Motivated Force Fields From SymmetryAdapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053–2066. 26


Page 26 of 30

Page 27 of 30

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


(20) Colwell, S. M.; Handy, N. C.; Lee, A. M. Determination of Frequency-Dependent Polarizabilities Using Current Density-Functional Theory. Phys. Rev. A 1996, 53, 1316–1322. (21) Hesselmann, A.; Jansen, G.; Sch¨ utz, M. Density-Functional Theory-Symmetry-Adapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. J. Chem. Phys. 2005, 122, 14103–14119. (22) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845–854. (23) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (24) Hoover, W. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697. (25) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (26) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. (27) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. (28) Thole, B. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341–350. (29) 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. 27



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

(30) Andrade, J. D.; Bes, E. S.; Stassen, H.; Bo, E. S. Computational Study of Room Temperature Molten Salts Composed by 1-Alkyl-3-methylimidazolium Cations—ForceField Proposal and Validation. J. Phys. Chem. B 2002, 106, 13344–13351. (31) 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. (32) Tokuda, H.; Hayamizu, K.; Ishii, K.; Bin, A.; Susan, 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. (33) Liu, W.; Zhao, T.; Zhang, Y.; Wang, H.; Yu, M. The Physical Properties of Aqueous Solutions of the Ionic Liquid [BMIM][BF4 ]. J. Solution Chem. 2006, 35, 1337–1346. (34) Verevkin, S. P. Predicting Enthalpy of Vaporization of Ionic Liquids: A Simple Rule for a Complex Property. Angew. Chem. Int. Ed. Engl. 2008, 47, 5071–5074. (35) Liu, H.;

Maginn, E. A Molecular Dynamics Investigation of the Struc-

tural and Dynamic Properties of the Ionic Liquid 1-n-Butyl-3-methylimidazolium Bis(trifluoromethanesulfonyl)imide. J. Chem. Phys. 2011, 135, 124507–124522. (36) 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. (37) Pal, T.; Biswas, R. Slow Solvation in Ionic Liquids: Connections to Non-Gaussian Moves and Multi-Point Correlations. J. Chem. Phys. 2014, 141, 104501–104512. (38) Gebbie, M. a.; Dobbs, H. a.; Valtiner, M.; Israelachvili, J. N. Long-Range Electrostatic Screening in Ionic Liquids. Proc. Natl. Acad. Sci. USA 2015, 112, 7432–7437. 28


Page 28 of 30

Page 29 of 30

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


(39) Hunger, J.; Stoppa, A.; Schrödle, S.; Hefter, G.; Buchner, R. Temperature Dependence of the Dielectric Properties and Dynamics of Ionic Liquids. ChemPhysChem 2009, 10, 723–733. (40) Wakai, C.; Oleinikova, A.; Ott, M.; Weingärtner, H. How Polar Are Ionic Liquids? Determination of the Static Dielectric Constant of an Imidazolium-Based Ionic Liquid by Microwave Dielectric Spectroscopy. J. Phys. Chem. B 2005, 109, 17028–17030.

29




Page 30 of 30

First-Principles United Atom Force Field for the ... - ACS Publications

Recommend Documents