First-principles parameterization of polarizable coarse-grained force

Subscriber access provided by READING UNIV

Article

First-principles parameterization of polarizable coarse-grained force fields for ionic liquids Frank Uhlig, Johannes Zeman, Jens Smiatek, and Christian Holm J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00903 • Publication Date (Web): 22 Jan 2018 Downloaded from http://pubs.acs.org on January 24, 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 free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

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

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

Journal of Chemical Theory and Computation

First-principles parameterization of polarizable coarse-grained force fields for ionic liquids Frank Uhlig,∗ Johannes Zeman, Jens Smiatek, and Christian Holm Institute for Computational Physics, University of Stuttgart, D-70569 Stuttgart, Germany E-mail: [email protected] Phone: +49(0)711/685-63607. Fax: +49(0)711/685-63658

Abstract We present an ab initio parameterization scheme for explicitly dipole-polarizable force fields for simulation of molecular liquids. The scheme allows for, in principle, arbitrarily coarse-grained representations. All parameters in the force field are derived from first principles, based on simple physical arguments. Only one fit parameter enters the parameterization, a global scaling factor for the size of the particles, which is adjusted to reproduce the experimental mass density. As important examples and for the first time, polarizable coarse-grained force fields are derived for 1-alkyl-3methylimidazolium cations with varying alkyl-chain lengths (alkyl=ethyl,butyl,hexyl), and hexafluorophosphate and tetrafluoroborate anions. Our findings are in good agreement with experimental results and results of further atomistic simulations. Hence, the force fields can be faithfully used where polarizability is expected to play a significant role, such as simulations of energy storage devices.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

Molecular dynamics (MD) simulations have become an indispensable tool for modern computational research. Accurate interaction potentials are readily available in many program packages for MD. However, parameter derivation for these potentials is neither obvious nor unique. They can be determined by performing extensive fits to experimental data, or deriving them from first principles. Notable large-scale efforts for non-polarizable force fields are OPLS, 1 PARM, 2 GROMOS, 3 and CHARMM. 4 Parameterizing force fields from first principles is becoming more and more a focus of recent force field developments. 5–7 The hope in using first-principles-based force fields is that they could alleviate the drawback of having to perform extensive parameter sweeps and rely on effective representations that in the end might only correspond to a certain state point. This implies that the potentials chosen have to cover all relevant physical effects driving the molecular motion and thus reproduce thermophysical properties. Hence, these potentials should also include the explicit treatment of polarization effects, i.e., changes in the charge distribution, or multipoles, respectively. However, the explicit treatment of many-body effects comes at an elevated computational cost. The many-body polarization problem has to be solved self-consistently, and one needs to include additional interaction sites and potentials. The polarizability is commonly restricted to inducible dipoles, whereas also force fields with higher-order inducible multipoles exist. 8 Inducible dipoles are then either treated as mathematical dipoles 9,10 or approximated as Drude dipoles. 11 The self-consistency requirement can be circumvented by dual-thermostatting schemes 11 or explicit propagation of the dipoles using time-reversible predictor-corrector methods. 12 Corresponding, efficient implementations have become widespread in major MD packages. 10,11,13–15 If many-body effects are not explicitly treated in molecular simulations, one can still account for them via effective pair potentials. Polarizability can be neglected for systems with little polarizability and polarization. However, for the simulation of highly polar, and highly polarizable liquids, polarizability needs to be included, even if only in a mean-field 2


Page 2 of 55

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


manner. 16,17 In a like-minded fashion, effective pair-potentials can be parameterized for the interactions of coarse-grained representations of bundles of atoms. Systematic ways of parameterizing such force fields have been described in the literature 18–22 and commonly rely on fitting to either experimental or theoretical reference data. Attempts at effectively including missing reorientational polarizability have been successfully used to improve existing coarsegrained models. 23 We have recently applied a coarse-graining scheme to the molecular polarizability for an already exisiting force field, yet still relying on empirical parameterizations of van-der-Waals and electrostatic interactions. 24 However, to the best of our knowledge, no explicitly dipole-polarizable coarse-grained force fields directly derived from ab initio data exist to date. One system class in which polarizability has a significant influence, particularly on dynamical properties, is formed by ionic liquids. 25–27 Room-temperature ionic liquids (RTILs) are salts that are molten below 100 ◦ C. They have attracted a lot of attention in recent years, both from experimental 28–30 as well as theoretical sides (ref. 31 and references therein). RTILs draw their appeal from being (close to) non-volatile, having a large electrochemical window, and being barely combustible. Hence, they are prime candidates for potential applications in fuel cells, 28,29 catalysis, 30 and energy storage devices (e.g., as electrolytes in supercapacitors). 32–34 RTILs have been investigated quite extensively using methods of molecular simulation. 31,35 Both ab initio MD as well as force-field-based MD simulations have been used to determine dynamics and structure in ionic liquids. 27 Accurate molecular simulations of ionic liquids are, however, hampered by substantial computational costs due to slow dynamics in these systems. To obtain converged properties from molecular simulations, one needs to calculate hundreds of nanoseconds for system sizes of at least 500 ion pairs of substantial individual size (tens of atoms). 36,37 Computational treatment of RTILs requires careful treatment of polarization effects due to strong internal electric fields. Both polarizable and non-polarizable force fields exist for ILs, 25,27 and many combinations of anions and cations are possible to create a plethora of different ionic liquids. 31 Hence,

3



researchers have aimed at providing generic, majorly non-polarizable, force fields to cover a wide range of possible chemical structures. 38–51 The mentioned force fields use general atom types and interaction potentials, and thus aim at transferability, which often comes at a loss of accuracy. Besides these united- or all-atom force fields, current studies focus more on specific model representations. These representations can be either significantly more complicated models, with explicit description of the electronic structure, 52 or also simpler descriptions as charged hard spheres. 53 In ionic liquids, the Coulombic forces lead to strong electric fields that distort the electronic clouds of the ions and induce significant molecular dipole moments. 54–56 This makes them an interesting test bed for the development of accurate polarizable force fields. 25,27,57 One of the most efficient approaches to include polarizability is to scale down the magnitude of the charge on the ions, 16,17,25 corresponding to an implicit, mean-field treatment of electronic polarizability. Reduced-charge models have also been justified semi-empirically from effective ion charges obtained from ab initio molecular simulations, 58 and have been successfully employed in reparameterizations of existing non-polarizable force fields. 59–61 Researchers have gone beyond treating polarizability in a mean-field manner. To this end, accurate polarizabilities can be obtained from quantumchemical calculations for specific compounds, 62–64 from empirical fits to experimental 65,66 or theoretical data. 67–69 Other important non-bonded interactions are exchange repulsion and dispersion attraction, together commonly referred to as van-der-Waals (vdW) interactions. Recently, successful schemes to derive complete parameter sets for atomistic force fields, including those for vdW interactions, from first principles have been reported. 5,6 Notably, therein, also correlations between the molecular volume and distributed dispersion coefficients have been exploited. 6,70 Alternatively, symmetry-adapted perturbation theory (SAPT) has been suggested for force field parameterizations 7 and successfully used to parameterize force fields for ionic liquids. 71,72 However, SAPT parameterizations require many different input structures to obtain reasonable force field parameters, akin to fitting potentials to reproduce forces in ab initio MD simulations. 73 To achieve transferability in these

4


Page 4 of 55

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


tailored force fields large data sets are needed. 74,75 Due to the high computational cost of modeling ionic liquids, already over a decade ago, coarse-grained force fields were developed for their simulation. 76,77 In recent years, further coarse-grained models of ionic liquids were created based on atomistic representations. 78–81 None of these coarse-grained force fields explicitly include electronic polarizability. For instance, the coarse-grained force field in ref. 79 uses charge-scaling to implicitly account for effects of electronic polarizability and thus overcome issues with slow diffusion. However, as mentioned before, this charge-scaling limits the applicability of these force fields. Only recently, we introduced an explicitly polarizable force field for the ionic liqiud 1-butyl-3methylimidazolium hexafluorophosphate. 24 In this work, we show how to generally derive force field parameters from quantumchemical calculations of isolated molecules. Accurate quantum-chemical calculations provide input data, which are transformed into a force field description using simple and physically intuitive arguments. Using volume-scaling arguments, these non-transferable, yet accurate, pair-potentials can be generated for arbitrary sub-groups of atoms in a molecule. This way, we parameterize force fields for coarse-grained representations of the 1-alkyl-3-methylimidazolium (alkyl=ethyl,butyl,hexyl) cations in combination with hexafluorophosphate and tetrafluoroborate anions. We demonstrate good agreement of our molecular simulation results with experimental and theoretical data. The resulting force fields are, to the best of our knowledge, the first ab initio derived, polarizable, and coarse-grained force fields for ionic liquids. These force fields provide an accurate representation of system properties at large length and long time scales.

5



2

Page 6 of 55

Methods

2.1

General parameterization strategy

For parameterizing force fields, we need to determine a couple of key components: constituents of our system, the functional form of the force field, and parameters to go with this force field. For now, we envision our system’s components simply as a set of points that interact with each other via some force field. We then choose intuitive interaction potentials that are able to represent physical reality with reasonable accuracy. However, only with accurate values for the involved parameters will we be able to run meaningful molecular dynamics simulations. The interaction potentials are divided into bonded and non-bonded terms. The former describe chemical bonds, bond-angles, and torsions. We ignore bonded contributions for our models and treat them as rigid bodies. Bonded terms could in general still be included and could be parameterized using ab initio data and the non-bonded interactions described below. The non-bonded terms define the internal energy U in a system of N interacting sites as:

U (rij ) =

N X N X i

j>i

"

( 4εij

σij rij

12

−

σij rij

6 # (1)

) 2 (µi · µj ) − 3 (µi · rij ) (µj · rij ) qi q j qi (µj · rij ) rij + + + . 3 5 4πε0 rij 4πε0 rij 4πε0 rij The double-sum runs over all pairs of interaction sites i 6= j. The distance between two sites is given by rij = |ri − rj |. ε0 is the vacuum permittivity. Each interaction site is characterized by a charge qi , a van-der-Waals (vdW) diameter σi , the well depth of the vdW potential εi , the dipole-polarizability αi and a correspondingly induced dipole moment µi . The first term in eq. 1 is a 12-6-Lennard-Jones potential, which is split into Pauli repulsion and dispersion attraction. The remaining terms in eq. 1 describe electrostatics, namely charge-charge, charge-dipole, and dipole-dipole interactions. Dipoles µi on each site are 6


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


induced by the local electric field (assuming linear response) according to: " µi = α i E +

# X

(T (rij )µj + G(rij )qj ) ,

(2)

j6=i

where E is an external electric field, T (r) = (3rr − 1r2 )/(4πε0 r3 ) is the dipole-field tensor, and G(rij )qj the electric field due to the point charges qj , with G(r) = r/(4πε0 r3 ). Equation 2 needs to be solved self-consistently to obtain the induced dipole moments µi , which then enter eq. 1. The cost to induce such a dipole moment (µ2i /(2αi )) needs to be added to the interaction energy in eq. 1 to obtain the potential energy of the system. We approximate these dipoles by so-called “Drude oscillators”, which are pairs of charges of opposite sign and equal magnitude, where each pair is connected by a harmonic spring. 11 The inducible dipoles diverge if they are spatially too close to each other. 68 This artifact can be avoided by using an additional screening of the dipole-dipole interaction: 68,82

Fij = 1 −

s r ij ij exp (−sij rij ) , 2

(3)

with sij = tij /(αij )1/3 , where αij is the geometric average of the polarizabilities of the respective sites, and tij is a semi-empirical parameter. Instead of the tensor αi in eq. 2, we assume an isotropic polarizability αi . This isotropic site-polarizability is then also used throughout all of the force-field-based simulations. Depending on the symmetry of the molecule, the coupling between the individually isotropic site-polarizabilities will in general lead to an overall anisotropic molecular polarizability. Up to now, we abstractly referred to interaction sites. These sites should represent our particular molecules in some manner. Atoms are an obvious, and usually good choice as interaction sites. However, we do not want to restrict ourselves to an atomistic representation. Rather, our representation combines several atoms to form a coarse-grained (CG) interaction site. With a reasonable choice of parameters, this allows to maintain both accuracy and computational feasibility. We then approximate the molecular charge distribution by a set 7



Page 8 of 55

of point charges on these interaction sites. However, there is no unique mapping between the real charge density and these point charges. The charge-fitting scheme we employ aims at reproducing the charge distribution and its multipole expansion by a reciprocal-space fit of a model charge density to a quantum-chemically determined one. 83 First, we fit a set of point charges to an atomistic representation of our molecules. Second, charges of the actual interaction sites of the models are calculated from the respective centers of charge of the atoms comprising a particular site. For a coarse-grained representation of a molecule of at least two charged beads, this procedure guarantees that total charge and dipole are conserved between the atomistic and a coarse-grained representation. For a one-bead representation only the total charge can be conserved. The next ingredient for our force field are the polarizabilities of the respective interaction sites. Those polarizabilities mimic the electronic response due to the surrounding particles and their electrostatic interaction with them. Accurate values for these polarizabilities can be obtained from linear-response calculations. 63 In general, these calculations yield molecular polarizability tensors α. We now need to distribute per-site polarizabilities in order to reproduce the molecular polarizability α. Here, we can make use of known correlations between polarizabilities and other properties (for details see ref. 84 and references therein). One particular, and recurring observation is the correlation of an isotropic molecular polarizability α = Tr(α)/3 with the molecular volume vtot of the corresponding electronic cloud. 65,85 Hence, we chose to calculate the static polarizabilities of an interaction site i in our models from the corresponding ratio of its volume vi to the total molecular volume vtot and the isotropic molecular polarizability α:

αi =

vi α. vtot

(4)

The volume of individual interaction site i is obtained from Bader’s atoms-in-molecules (AIM) analysis. 86,87 We define arbitrary interaction sites by concatenating respective atoms,

8


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


and their volumes. Because the volumes in an isolated molecule are not strictly bound, one needs to define a density cutoff beyond which the volume is not summed up anymore. We have chosen a value for the electron density contour of 0.0026 e/a30 , with the elementary charge e and the Bohr radius a0 . The cutoff provides reasonable parameters for the van-derWaals potential and yields a good guess for molecular volumes in our coarse-grained representation. This value is comparable to the one used in parameterizations of semi-empirical dispersion corrections. 88 The overall procedure provides a simple and intuitive alternative to more complicated procedures that are also based on the AIM analysis. 62 Nevertheless, the static polarizability tensor is well reproduced with this simple method (vide infra). Polarizabilities obtained from the AIM analysis are then transformed into parameters for the Drude oscillators. The polarizability αi of a Drude oscillator is given by: 11

αi =

qi2 , k

(5)

where qi is the charge on the Drude particle i, and k is the harmonic force constant connecting the Drude pair. It is formally only valid to approximate a mathematical dipole by a Drude oscillator for small separations. The overall separation of the Drude pair is governed through the force constant, and deviations between simulations with either inducible point-dipoles or Drude dipoles have been observed for high polarizabilities. 89 The force constant k was set to 1000 kcal/mol/Å2 for all Drude oscillators, as successfully employed in biomolecular force fields. 90 Increasing the value by an order of magnitude did not lead to any significant changes in structural and dynamical quantities. Lastly, we describe how to obtain the parameters in the vdW potential. A 12-6-LennardJones potential (first term in eq. 1) approximates the van-der-Waals interactions. The attractive term proportional to r−6 describes the dispersion interaction. The corresponding interaction strength can be calculated using quantum-chemical methods. The C6 dispersion coefficients, also called Hamaker constants, are generally calculated from the following

9



integral: 91 3 C6 = π

Page 10 of 55

Z∞ α(iν)α(iν)dν,

(6)

0

where α(iν) is the isotropic dynamic polarizability at the imaginary frequency ν. Hamaker constants for molecules can be obtained from both experiment, e.g., from dipole-oscillator strength distributions 92,93 and theory, 94,95 but need to be distributed to individual interaction sites for use in force fields. From eq. 6 one can see that C6 coefficients correlate with the square of the polarizability, and thus also with the square of the volume. Hence, we adopt the same argument as before and calculate C6 coefficients of individual sites from the molecular C6 coefficient and the squared ratio of the volume vi of site i to the total molecular volume vtot : C6i

=

vi vtot

2 C6 .

(7)

Similar arguments have been brought forward in the calculation of atomic dispersion coefficients for pairwise corrections to density functional theory 70 which in turn have also been used in the parameterization of force fields. 6 In these parameterizations, dispersion coefficients of atoms were rescaled according to the square of the volume ratio of the bound atom (in a molecule) compared to its free volume. Importantly, by virtue of the volume scaling, both methods to obtain distributed C6 coefficients reproduce the decrease in dispersion interaction when transferring free atoms to chemically bound systems. Furthermore, the individual volumes obtained from the AIM analysis are connected to each other and encompass continuous regions in space for neighboring atoms. Hence, arbitrary groups of connected atoms can be coarse-grained into respective interaction sites. Parameters σ and ε of the Lennard-Jones-Potential cannot be obtained from quantumchemical calculations directly. However, the C6 dispersion coefficient relates σ and ε through C6 = 4εσ 6 . Hence, we need to determine either σ or ε. In our AIM analysis, we have already obtained the quantum-mechanical volume of an interaction site. Assuming a sphere with this volume, we obtain a reasonable guess for the vdW diameter σ, and thus also 10


Page 11 of 55 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 depth of the Lennard-Jones potential well ε = C6 /(4σ 6 ). The σ’s are then adjusted to obtain the correct bulk liquid density at one thermodynamic state point while keeping the C6 coefficients constant. The heterogeneous parameters σij and εij between different sites i and j are obtained from the homogeneous analogues through geometric combination √ √ rules, i.e., σij = σi σj and εij = εi εj . This choice guarantees that the total dispersion coefficient stays the same also after distributing it among interaction sites. The second part of the vdW interaction, the exchange repulsion, is described by the term proportional to r−12 . There are functional forms that describe the exchange repulsion better, 96 but we chose this representation because of sufficient accuracy and common availability. The involved parameters are the same as for the dispersion attraction. All of the aforementioned potentials are commonly implemented in molecular simulation packages and hence widely available. Simulation software often differs in the available algorithms for the explicit treatment of polarizability. Details of our simulation setup follow in the next section.

2.2

Computational methods and algorithms

We performed both static, quantum-chemical calculations, and molecular dynamics simulations. Different methods and program packages were used due to differences in available implementations. We start by describing the calculations needed for the force field parameterization, followed by the methodology of the molecular dynamics simulations. Static density functional theory calculations on isolated ions were performed with Orca 97 and the CP2K 98 program suites. We used density functional theory (DFT) as implemented in the Quickstep module 99 of CP2K. Quickstep employs a dual basis-set representation. The Kohn-Sham one-electron wavefunctions were expanded into atom-centered dual- or triple-ζ Gaussian basis sets with polarization functions (“DZVP-MOLOPT-SR”, “TZV2P-MOLOPT”, and “aug-TZV2P”). 100 The electronic density was also expanded into a plane-wave basis set using a kinetic energy cutoff of 400 Ry. Goedecker-Teter-Hutter-type 11



pseudopotentials 101 replaced the core-electrons. Kohn-Sham equations were then solved in the generalized-gradient approximation to DFT, with the revised Perdew-Burke-Ernzerhof (revPBE) density functional. 102 Periodic images were decoupled using density-derived atomic point charges (DDAPC). 83 Calculations of polarizabilities of isolated molecules used variational density functional perturbation theory. 103,104 Polarizabilities were then obtained with various basis sets. With Orca, density functional theory calculations employed as well the revPBE density functional approximation, and coupled-perturbed self-consistent field equations. Furthermore, we performed second-order perturbation theory calculations with Orca. In these Møller-Plesset perturbation theory (second order, MP2) calculations, the dipole moment was differentiated numerically to obtain the polarizability. 63 Again, polarizabilities were then calculated using dual- and triple-ζ basis sets. 105,106 Polarizabilities for use in the force fields were calculated on isolated molecules as the dielectric response of the solvent is treated explicitly in our molecular dynamics simulations. The change of dipole moment in the bulk liquid, in response to an external electric field, was probed using an applied electric field. The applied electric field was spatially homogeneous, with a field strength of 0.01 a.u. (atomic units), which is low enough to stay in the linear-response regime. 63 We performed three one-sided finite-difference calculations with applied fields in x-, y-, and z-direction to obtain the respective responses. These responses can be calculated equivalently in the bulk DFT as well as in the polarizable force field calculation. By inverting eq. 2, this procedure yields polarizabilities of our system. The so-obtained polarizability can be split into molecular contributions using maximally localized Wannier functions 107–109 (MLWF) in the DFT calculations, and according to the molecular topology in the force field calculations. Dispersion coefficients for all species were calculated with real-time time-dependent density functional theory calculations using the Octopus code. 110 The real-time response of the density to an external perturbation was recorded and integrated numerically to obtain the frequency-dependent, dynamic polarizabilities. 111 The initial perturbation to the elec-

12


Page 12 of 55

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


tronic structure was 0.00037 a.u. 111 in magnitude, which provided stable electronic dynamics. For calculation of dispersion coefficients we used self-interaction-corrected 112 density functional theory calculations in the local density approximation, which provide results in good agreement with experimental values. 111,113 The density was represented on a grid composed of individual spherical grids centered on the atoms. The grid spacing was set to 0.3 a.u. Core-electrons were replaced using norm-conserving pseudopotentials of the HartwigsenGoedecker-Hutter type. 114 All the obtained parameters and simulations were used to parameterize force fields for ionic liquids containing the cations 1-ethyl-3-methylimidazolium, 1-butyl-3-methylimidazolium, or 1-hexyl-3-methylimidazolium, and anions hexafluorophosphateor tetrafluoroborate. The final force field parameters are all given in the supporting information (SI). We performed both ab initio molecular dynamics (AIMD) simulations, and force-fieldbased molecular simulations (FFMD). FFMD simulations were performed using either atomistic or coarse-grained models. The parameters for the atomistic force fields were obtained from available force fields for ionic liquids. 38,115 The coarse-grained simulations used force fields as described in the preceding section, as well as previously reported non-polarizable force fields. 78,79 All FFMD simulations were performed with LAMMPS 116 (http://lammps. sandia.gov/index.html). Simulations were prepared with the FFTOOL 117 program and adapted manually. Newton’s equations of motion were integrated with a timestep of 1 fs, or 5 fs 78 for simulations using an atomistic, or those using a coarse-grained representation, respectively. Pairwise interactions were cut off at 1.2 nm. Long-range electrostatics were calculated using 3D-periodic boundary conditions, and the particle-particle-particle-mesh method. 118 The relative accuracy on the forces was set to 0.0001 because tighter accuracies did not have a quantitative effect on the dynamics. The long-range contribution to the dispersion energy and virial was computed assuming an isotropic distribution beyond the pairwise cutoff. Temperature and pressure were kept constant using Nosé-Hoover thermostats and Parrinello-Rahman barostats. 119–121 The coupling constants were 100 timesteps

13



for the thermostat, and 1000 timesteps for the barostat. All simulations were run at ambient pressure (1 bar), and different temperatures ranging from 300 to 450 K. In all simulations, chemical bonds including hydrogen atoms and all the internal degrees of freedom of the coarse-grained models were constrained using the RATTLE algorithm. 122,123 The systems consisted of 512 ion pairs in all the force-field-based molecular simulations, if not stated otherwise. In explicitly polarizable force fields, the motion of the Drude particles on all sites was integrated using a time-reversible always-stable predictor-corrector (ASPC) algorithm. 12 We implemented this algorithm using the infrastructure of the recent Drude-polarizability implementation 13 into LAMMPS. The ASPC routines are currently not included in the official LAMMPS distribution, but can be downloaded from GitHub. 124 The force constants for all Drude oscillators were set to 1000 kcal/mol/Å2 . The damping factor used in the ASPC integration was set 0.6, which provided stable and accurate time integration. The coulombic interaction between Drude-oscillators was screened at short range using a modified 82 version of the damping proposed in ref. 68, to avoid a polarization catastrophy. The screening factor was set to 2.0 as successfully used in ionic liquid simulations before. 72 Initial structures were prepared using Packmol 125 by randomly placing molecules into a box of dimensions corresponding to the respective experimental mass density. As explained in the previous section, the σ parameters of the Lennard-Jones potential in our models were adjusted to reproduce the experimental mass density of the solutions. Initial test runs used the guess for the σ’s as obtained from the quantum-mechanical volume of the atoms in the molecule. The parameters were then rescaled using the third root of the ratio between current mass density and experimental mass density. One to two trial runs sufficed to accurately reproduce the experimental mass density. The systems were then equilibrated for about 20 ns, and final production runs were at least 100 ns long. For the AIMD simulations, 10 ion pairs of 1-butyl-3-methylimidazolium hexafluorophosphate were randomly placed in a (1.507 nm)3 cubic box, again using Packmol, 125 which

14


Page 14 of 55

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


yields a mass density of 1.38 g/cm3 , close to the experimental density at 300 K. 126,127 The structure was then equilibrated in the isothermal, canonical ensemble with the CL&P force field 38,115 using the methodology described above. The resulting structure was then used as input to the AIMD simulations. The DFT setup is the same as described above, using the MOLOPT-DZVP-SR basis set. Additionally, a pairwise-additive term 128,129 was included to account for the lack of dispersion interaction within the revPBE density functional. Molecular dynamics were performed for 30 ps, from which the last 20 ps were used for analysis. Molecular dipole moments along the trajectory were obtained from maximally-localized Wannier functions.

3

Results and discussion

In this section we will describe the parameterization of our force fields using as explicit example 1-butyl-3-methylimidazolium hexafluorophosphate ([C4 C1 Im][PF6 ]). We first validate the choice of methods for our parameter calculations, and proceed with validating our force field by explicitly comparing our data to ab initio molecular dynamics simulations, available experimental data, and other molecular simulations. The parameterization process is then applied to ionic liquid species with ethyl and hexyl side chains, and the BF4 – anion.

Force field parameterization at the example of [C4 C1 Im][PF6 ] One of the first issues in generating a force field is the determination of an appropriate molecular structure. The chemical structure of our model system, [C4 C1 Im][PF6 ] is shown in Fig. 1. Also displayed is the coarse-grained representation of the respective ions, in which butyl-, methyl-, imidazolium-, and hexafluorophosphate-groups are collected into respective coarse-grained beads (interaction sites). This coarse-grained representation has been used successfully before in non-polarizable force fields. 78,79 If not stated otherwise, this coarsegraining was employed in all force-field-based simulations. Throughout this work, we use

15



coarse-grained bead representations where bead positions are obtained from centers of charge of the constituting group of atoms (see also Methods section). The coarse-grained, polarizable simulations are then generally referred to as “CGPol”. A more general mapping scheme from atomistic to coarse-grained representation for ionic liquids could be inferred from common functional groups and used to parameterize other species as well.

Figure 1: Chemical structure of the C4 C1 Im+ cation and PF6 – anion (left), and their respective coarse-grained representations (right).

In this work, we aim at reproducing properties of bulk liquid solutions. As all the coarsegrained models are treated as being rigid, this means we need to use a molecular structure in some way representative of the bulk liquid. Hence, a 3D-periodic liquid of [C4 C1 Im][PF6 ] was simulated using ab initio molecular dynamics, which does not make any assumptions on the molecular structure. The extracted internal, bonded degrees of freedom were then compared with those in a widely used atomistic force field (CL&P). The explicit values are given in the supporting information. The force field values closely agree with ab initio simulations. Hence, we opted to use the structures in the CL&P force field, because it allows us to quickly, yet faithfully, investigate and parameterize force fields for a large set of chemically similar species, without having to perform computationally expensive AIMD simulations. 16


Page 16 of 55

Page 17 of 55 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 structures are included in the supporting information. Except for the comparison of molecular structures, we do not extract any other quantities for the parameterization from the bulk liquid AIMD simulations. Extracting, e.g., charges and polarizabilities from ab initio liquid phase calculations would mean to implicitly account for induction and manybody effects present in the liquid. 58 All of our final force fields include explicit treatment of many-body polarization effects. Thus, we rather extract molecular properties from gasphase calculations and investigate how inclusion of many-body effects in the bulk FFMD simulations influences properties, and compare to corresponding ab initio results. Arguably, the most important ingredient needed to faithfully reproduce many-body effects in our simulations are the dipole-polarizabilities. We obtain them from quantum chemical calculations on the individual cations and anions comprising our ionic liquids. Resulting diagonal components of the dipole-polarizability tensors of C4 C1 Im+ and PF6 – (in the geometry given in the SI) are shown in Tab. 1. These polarizabilities in vacuum strongly depend Table 1: Static polarizabilities obtained with different quantum-chemical methods and basis sets for C4 C1 Im+ . All MP2 calculations, and DFT calculations with the (d-)aug-cc-pVD/TZ basis set were performed with Orca, all other DFT calculations with CP2K. The geometry of the ions is given in the supporting information. All values in Å3 . method

basis set C4 C1 Im+ RI-MP2 aug-cc-pVDZ revPBE aug-TZV2P Drude heavy atoms Drude coarse-grained PF6 – RI-MP2 aug-cc-pVDZ revPBE aug-TZV2P Drude heavy atoms Drude coarse-grained

αxx 19.24 19.62 19.21 16.8

αyy

αzz

15.43 10.98 15.33 10.67 13.78 9.31 14.59 13.9 4.67 4.65 4.53 4.65

on the description of the diffuse nature of the electronic cloud. Hence, RI-MP2 calculations employing augmented, diffuse basis sets (aug-cc-pVDZ) serve as reference calculations. Further augmentation, or using a triple-ζ basis set, did not lead to significant changes. Table 1 also includes polarizabilities obtained with density functional theory as implemented 17



in CP2K. As CP2K is our workhorse for bulk DFT calculations, we want to verify that polarizabilities are well represented with available methods and basis sets. Accurate values are obtained, compared to the MP2 reference, using triple-ζ basis sets (aug-TZV2P) and the revPBE density functional. The relative error of the DFT-calculated diagonal components of the polarizability tensor is below 3 %. We thus use this method and basis set throughout all force field parameterizations. This combination does not just provide good accuracy in the gas phase calculations, it also gives access to calculations of polarizabilities of large molecules. As explained in the methods section, calculated molecular polarizabilities need to be converted to the force field description and hence distributed among the interaction sites. We calculated distributed polarizabilities for an atomistic representation and our coarse-grained representation. The polarizable, atomistic representation was only used to investigate the polarizability tensor and not further in molecular simulations. For atomistic calculations of the polarizability tensor we attach Drude oscillators to all but hydrogen atoms. The volume of the hydrogen atoms was added to the volume of the connecting heavy atom. In the coarse-grained representation, the volumes of the interaction sites are obtained from the sum of the volumes of constituent atoms (cf. Fig. 1). The polarizability tensor of the molecule is reproduced well using Drude-oscillators and the atomistic representation (Tab. 1). The coarse-grained representation qualitatively reproduces the polarizability tensor, but the diagonal components are much closer to each other in magnitude than in the quantumchemical calculations. The coarse-grained representation is likely missing too many chemical details to quantitatively describe the anisotropy of the polarizability. In Fig. 2, we show the change in diagonal components of the polarizability tensor with respect to the short-range screening parameter tij introduced in eq. 3 (solid lines). The screening parameter dampens the short-range coulombic interaction to avoid polarization catastrophes. It quite strongly influences the split of the components of the polarizability tensor. Reference values from revPBE//aug-TZV2P calculations are shown as dashed lines. The change in the diagonal

18


Page 18 of 55

Page 19 of 55

αxx αyy αzz

30

25

α ˚ 3] αii [A

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

15

10

5 1.0

1.5

2.0

2.5

1.0

1.5

tintramol

2.0

2.5

tintramol

Figure 2: Diagonal components αii (i = {x, y, z}) of the molecular polarizability tensor of C4 C1 Im+ from DFT (dashed lines) and from Drude dipoles (lines with points). Atomistic representation with Drude dipoles on heavy atoms (left), and coarse-grained (right) in dependence of the Thole screening parameter tintralmol , αxx -component in blue, αyy in green, and αzz in red. Corresponding isotropic values α are shown in black. components of the polarizability tensor is stronger in the atomistic representation than in the coarse-grained one. As the screening parameter varies, the atomistic representation both over- and underestimates the polarizability. However, the screening parameter barely influences the polarizability in the CG representation, but the overall agreement with the reference polarizability is worse. The CG representation separates the Drude-oscillators farther than the atomistic representation and thus the screening has a smaller impact. In the end, we selected a screening parameter of 2.0 for our force fields. This value is close to the originally suggested value of 2.1, 68 and has already been used successfully in ionic 19



liquid force fields, 72 and yields reasonable accuracy with respect to our quantum-chemical calculations. To further elucidate the validity of our parameterization, we examine the isotropic polarizabilities and dispersion coefficients created by volume-scaling of the molecular properties (for details see Methods section). We construct a set of reference molecules by splitting the C4 C1 Im+ molecule into three parts: butyl-, methyl-, and imidazolium groups. The broken bonds were saturated with hydrogen atoms and used for separate calculations of polarizabilities and dispersion coefficients. Properties of these molecules are compared to their counterparts in the coarse-grained representation in Tab. 2. The polarizabilities of these coarse-grained fragments agree well with those of the reference molecules. The discrepancies between dispersion coefficients are larger than for the polarizabilities. However, a small change in polarizability will manifest itself more strongly in the dispersion coefficient as the latter is proportional to the integral of the square of the dynamic polarizability (cf. eq. 6). We derive parameters for the Lennard-Jones potential from values given in Tab. 2. The σ Table 2: Polarizabilities α and dispersion coefficients C6 on individual sites, derived from scaling with volume ratio, and from equivalent, hydrogen-saturated, molecular fragments. +

C4 C1 Im Im butane methane

αtot [Å] 15.8 6.0 7.5 2.3

vpart /vtot · αtot [ Å] 15.2 5.7 7.4 2.1

C6 [a.u.] 5883.8 807.6 1335.5 137.9

(vpart /vtot )2 · C6 [a.u.] 5428.2 760.9 1288.9 103.9

parameter is calculated as the diameter of a sphere with volume v, and the ε parameter from the relation C6 = 4εσ 6 . The final values are then obtained by scaling the σ values to reproduce the experimental mass density (at T=350 K), while keeping the dispersion coefficients fixed. The volume-scaling procedure, for both dispersion coefficients and polarizabilities, is in particular useful for coarse-graining approaches. It allows to dissect molecules into arbitrary, connected sub-groups of atoms, and still obtain physically reasonable parameters for the resulting (coarse-grained) beads. The missing ingredient for our force field are the charges on the interaction sites. As men20


Page 20 of 55

Page 21 of 55 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


tioned previously, we employ so-called “density-derived atomic point charges” (DDAPC). 83 These charges are straightforwardly calculated from both ab initio and force field molecular dynamics simulations and have been successfully employed in force field parameterizations for ionic liquids. 58 As a sensitivity analysis and a posteriori verification of the coarse-grained representation, we calculated the DDAPCs with respect to rotation around the torsional angle connecting the methyl-imidazolium and the butyl parts of C4 C1 Im+ . The average charges are 0.63 e on the imidazolium, 0.15 e on the methyl, and 0.22 e on the butyl group. Their respective standard deviations are of the order of 0.01 e or less. Hence, the structurally important, torsional degree of freedom 130,131 only marginally influences the DDAP charge distribution, and can thus, in the CG representation, be captured by a rotation of the molecule around the methyl-imidazolium axis. The butyl chain itself is quite flexible. In atomistic molecular dynamics simulations, the distribution of distances between the centers of mass of the butyl chain and the imidazolium group are bimodal, with the main peak close to our model parameters (see SI). The final parameters for the [C4 C1 Im][PF6 ] model are given in Tab. 3. The table shows non-bonded and bonded parameters. As explained earlier, there are no bonded degrees of freedom. Hence, the bonded parameters are used only to define the structure of the molecules. Table 3: Parameters for [C4 C1 Im][PF6 ] model

site Im C1 C4 PF6

Im–C1 bond Im–C4 bond C4 –Im–C1 angle

nonbonded parameters mass [amu] charge [e] σ [Å] 67.070 0.6344 4.585 15.035 0.1508 3.290 57.115 0.2148 5.006 144.960 -1.0000 5.152 bonded parameters value 3.489 Å 4.352 Å 130.36◦ 21


[kJ/mol] α [Å3 ] 1.180 5.693 1.180 2.103 1.180 7.409 0.775 4.653


Model validation We now compare properties extracted from simulations of our newly parameterized [C4 C1 Im][PF6 ] model to both experimental as well as other theoretical data. Fig. 3 shows the mass densities of [C4 C1 Im][PF6 ] solutions at varying temperature.

1400

ref. 117 ref. 118 CGPol RoyMa CL&P

1350 ρ [kg/m3]

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

Page 22 of 55

1300

1250

1200 300

320

340

360

380 T [K]

400

420

440

Figure 3: Mass density ρ of [C4 C1 Im][PF6 ] as a function of temperature as calculated in simulations with explicit (“CGPol”) and implicit polarizability (coarse-grained: RoyMa, 79 all-atom: CL&P 38,115 ) and two experimental sources. 126,127

Also included are values for the mass density obtained with two other, non-polarizable force fields. The CL&P 38,40 force field uses an atomistic representation, and the RoyMa 79 force field uses a coarse-grained representation of [C4 C1 Im][PF6 ]. The latter force field uses the same general abstraction of atoms into coarse-grained interaction sites as our model.

22


Page 23 of 55 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


All force fields reproduce the trend in the experimental mass density well. However, both the RoyMa and the CGPol force fields were fitted to reproduce experimental mass densities. Hence, the good comparison to experimental data is not surprising. Nevertheless, the CGPol force field was only fitted to reproduce the density at one particular state point (T = 350 K), and it is not necessarily given that the mass density at different state points is reproduced accurately as well. It is noteworthy that although the CL&P force field underestimates the experimental data, deviations of the different force fields are within each others’ errors. The dependence of the density on the temperature differs between the various models. The density of the CGPol models decreases more steeply with increasing temperature than with the other models. However, none of the theoretical curves is parallel to the experimental data. We exclude the data points at 300 K from further analyses, because it is close to the melting point of [C4 C1 Im][PF6 ] (283 K 132 ) and the CGPol model showed sudden jumps in the potential energy. The density of this, possibly supercooled, solution was not affected. We examine the structure within the solution using radial distribution functions (RDF) between centers of mass of various groups. In Fig. 4 the RDFs between ions, i.e., anion-anion (AA), cation-cation (CC), and anion-cation (AC) at a temperature of 350 K are shown. Data are shown for different force field simulations and our ab initio simulations.

23



CGPol RoyMa CLaP AIMD

4

g(r)

3

AC

2 1 0 4

g(r)

3 AA

2 1 0 4 3

g(r)

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

CC

2 1 0

5

10

15 20 ˚ r [A]

25

Figure 4: Radial distribution functions g (r) between anions and cations (AC, top), anions and anions (AA, middle), and cations and cations (CC, bottom). Results obtained with AIMD in light blue (only for AC), CGPol in dark blue, CL&P in red, and RoyMa in green. Individual groups are represented by their respective centers of mass.

The range of distances is limited for the AIMD data, as the unit cell contained only ten 24


Page 24 of 55

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


ion pairs. Hence, we only included data for the anion-cation RDF, where the first ion shell is clearly visible. From the AIMD data we can see that the ion shells, and thus sizes, of the ions are well represented in all force field simulations. The agreement between the CGPol model and other force field simulations is best for the anion-anion distribution, while the cation size appears to be underestimated. This underestimation is visible in both anion-cation, and cation-cation RDFs. However, the structure of the cation-cation RDF for small values of r differs quite significantly between all data. The anion in the CGPol simulation is also smaller than in all other data, but to a lesser extent than the cation. Hence, the RDFs with CGPol generally appear compressed compared to the other data. Furthermore, both CL&P and RoyMa force fields predict a significant shoulder to the first peak in the anion-anion RDF due to the asymmetric coordination of the cation. This shoulder is less pronounced in the CGPol force field. The smaller shoulder further leads to an apparent compression of the RDF. The long-range-structure decays similarly for all RDFs and fluctuates with similar magnitude at long distances. Overall, all simulations agree reasonable well with each other. This can also be seen from RDFs between the individual sites (imidazolium, methyl, and butyl), which are included in the SI. To go beyond simple structural properties, we calculate both polarization and polarizability in our bulk liquids. Bulk-phase dipole moments in the AIMD simulations are calculated from maximally-localized Wannier functions. These localized functions allow us to dissect the total dipole moment into individual contributions of separate molecules. The distributions of C4 C1 Im+ ’s and PF6 – ’s dipole moments, from both AIMD and coarse-grained molecular dynamics simulations, are shown in Fig. 5.

25



3.0

CGPol C4C1Im+ CGPol PF− 6

2.5 normalized frequency

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

Page 26 of 55

AIMD C4C1Im+ AIMD PF− 6

2.0 1.5 1.0 0.5 0.0 0

1

2

3

4

5

6

7

8

9

µ [D] Figure 5: Distributions of dipole moments µ in CGPol (solid lines) and AIMD (dashed lines) simulations of C4 C1 Im+ (blue) and PF6 – (red) in bulk liquid at 350 K. Upon solvation in the liquid phase, the dipole moment of C4 C1 Im+ increases by about a third of its gas-phase value of 4.5 D to about 6 D. This strong polarization is predicted by both AIMD and CGPol simulations. The large induction signifies the strong polarization from the surrounding (counter)ions. The dipole moment of the anion increases as well, but dipole moments of the anions differ more between AIMD and CGPol descriptions. The difference is in the magnitude of the polarization, as well as in the width of the dipole moment distributions. This discrepancy can be traced back to the coarse description of the PF6 – anion. However, it is not a priori clear whether the description of electronic polarizability or the missing internal degrees of freedom are responsible for the too small polarization of the anion. Hence, we also investigate the response function α0 of the system with respect to an applied external field. The induced dipole moments just described signify the response of

26


Page 27 of 55

the molecules to the internal electric field in the bulk solutions. This response differs from its gas-phase equivalent. To quantify and validate our description of the electronic part of this response, we recorded the change in dipole moments in the liquid phase to an external, homogeneous, electric field E. Three one-sided differentiations of the dipole moments, one for each spatial direction, were carried out on snapshots of the CGPol and AIMD trajectories. From the response of the dipole moments, one can calculate the dipole polarizability of the (species in the) solution. However, the external field E not just induces a dipole locally in the sample, but this induction also leads to polarization of the other molecules within the sample, resulting in a changed local electric field Eresp . In principle, to get accurate bulk dipole polarizabilities, one needs to account for both these effects. The corresponding relation for the induced dipole moment is given in eq. 2. The inversion of eq. 2 to obtain α

3.0

CGPol C4C1Im+ CGPol PF− 6

2.5 normalized frequency

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


AIMD C4C1Im+ AIMD PF− 6

2.0 1.5 1.0 0.5 0.0 0

5

10

15

20

25

˚ 3] α 0 [A Figure 6: Distributions of the response-function α0 to an external electric field E in CGPol (solid lines) and AIMD (dashed lines) simulations of C4 C1 Im+ (blue) and PF6 – (red) in bulk liquid at 350 K.

27



can be performed by the finite difference calculations mentioned above. Here, the focus lies on the comparison of the total response of the system to an external perturbation. Hence, in all that follows, the re-radiated field of the induced dipoles Eresp is ignored, and we define a response δµ = α0 E. The distributions of the isotropic susceptibilities α0 for both C4 C1 Im+ and PF6 – are displayed in Fig. 6. The response within the force field description compares favorably with the ab initio data. Compared to the gas phase values, both the AIMD and CGPol simulations show an increase of the polarizability values for both anion, and the cation distributions. For the anion, only about 90 % of the response in the AIMD data is captured. In case of the cation, about 96 % of the ab initio response is recovered. Response of anion and cation within the liquid are tightly coupled to each other, and also to the employed charge distributions and representation. Hence, we see that the missing polarization of the anion is partly due to the lack in the (coupled) response to an electric field. However, overall more than 90 % of the response is recovered. The coarse description of the anion as a single sphere causes a loss in degrees of freedom that could lead to a significant polarization. The anion in its atomistic representation has a high quadrupolar moment due to a strong polarization of the bonds between the phosphorous atom and the six highly electronegative fluorine atoms. Thus, the anion can easily gain a significant dipole moment upon symmetry-breaking of the ideally octahedral coordination of fluorine atoms around the phosphorous atom. We recalculate dipole moments along our AIMD trajectory using the respective structures, but with arbitrary sets of point charges. This allows us to estimate the effect on the dipole moment by the distortion of the PF6 – geometry. Assuming a fixed point-charge distribution on the PF6 – anion, with average charges obtained from DDAPC fits on the bulk phase charge distribution, a dipole moment of about 0.4 D alone is induced by changes to the geometry of the anion. The missing reorientational, vibrational, and electronic polarizability in the anion lead to a substantially smaller response in the CGPol force field compared to the ab initio data. An accurate description of polarizability and polarization is important as it allows to

28


Page 28 of 55

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


correctly capture the intermolecular interactions in liquids. One can infer the interaction strength within a liquid from its enthalpy of vaporization (experimentally and theoretically). These enthalpies do not just measure the interaction, but also include the work required to change the volume at constant pressure. The latter is simplified due to the large difference of molar volumes in gaseous and liquid phases, and the assumption of an ideal gas. Furthermore, we neglect changes in molecular structure, as well as associated differences in ground-state energy. The enthalpy of vaporization Hvap is then given as: 1

∆Hvap (T ) = hEvap i −

hEliq i + RT. NIP

(8)

Evap/liq are average potential energies in the liquid and gas phase, NIP is the number of ion pairs in the liquid, R is the universal gas constant and T the temperature. The gas phase of ionic liquids consists predominantly of ion pairs. 49,133–135 This means that one needs to explicitly calculate the interaction energy in the gas phase and the difference in polarization energy between gaseous phase and bulk liquid. In explicitly polarizable simulations, this difference is inherently accounted for. This is in contrast to non-polarizable models, where this difference in polarization energy can only be accounted for effectively. 136 The polarization energy is given by: N

Epol =

1 X (µi − µ0i )2 , 2 i=1 αi

(9)

where µi − µ0i is the induced dipole moment of the ith Drude oscillator, αi is the corresponding dipole polarizability. In our simulations, the polarization energy differs between gaseous and liquid phase by about 21 kJ/mol. The associated gain in electrostatic attraction due to the increased dipole-dipole interaction is -40 kJ/mol, i.e., the average dipole-dipole interaction is significantly larger in magnitude in the gaseous phase. In a non-polarizable force field, the polarization is the same in both gaseous and liquid phase. This would imply a too high enthalpy of vaporization for an implicitly polarizable model. For our CGPol model, we calculate a heat of vaporization of 114 kJ/mol. The estimated, experimental value at 298 K 29



Page 30 of 55

is 155 kJ/mol, 137 and 129 kJ/mol at 353 K from other, explicitly polarizable, molecular dynamics simulations. 72 With our force field, the binding energy of an ion pair in the gas phase is about -337 kJ/mol, which is close to the value obtained from atomistic, polarizable force field simulations of -344 kJ/mol. 71 Similarly, for the liquid phase the values are -457 kJ/mol in our simulations, and -470 kJ/mol in ref. 71. Coulombic interactions make up about 99 % of the interaction energy in the liquid phase. Of this attraction, about 8 % are due to induced dipoles. In case of the gas-phase ion pair, the polarization is stronger as there is no screening due surrounding ions. Hence, the dipole-dipole contribution to the electrostatic interaction rises to about 20 %. Altogether, this demonstrates that interactions are dominated by strong Coulomb forces, and can be well reproduced by our force field. The description of diffusivity in ionic liquids has reportedly shown a strong influence on the inclusion of polarizability. 57,138 We thus calculate the self-diffusion coefficients from simulations with our CGPol force field and the aforementioned, corresponding non-polarizable force field. Self-diffusion coefficients D were obtained from the mean squared displacement of the center-of-mass trajectories of individual particles in the solution according to:

D=

N 1 1 X

1 lim |Ri (t) − Ri (0)|2 . 6 t→∞ t N i=1

(10)

Here, the average over all N particles of a particular species is already included. Figure 7 depicts diffusion coefficients for [C4 C1 Im][PF6 ] from experimental measurements as well as from our molecular simulations. The lines in Fig. 7 represent fits to experimental data, and the points correspond to molecular simulations. The simulational data comes from the CGPol model and an equivalent model with the same parameters but the explicit polarizabilities removed. In the latter model, the polarizability is neither explicitly nor implicitly accounted for. However, simply removing the polarizability from the CGPol model leads to a nonpolarizable model which still reproduces the mass density of the original model to within 1% without refitting the van-der-Waals parameters. Hence, this model appears to be well suited

30


Page 31 of 55

to directly study the influence of polarizability. All models provide diffusion coefficients that are close to the experimental data. An example plot showing the mean square displacement of the ions is given in the SI.

10−5

D [cm2/s]

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

exp. CGpol nonpol 10−7

2.2

2.3

2.4

2.5 2.6 1000/T [K]

2.7

2.8

2.9

3.0

Figure 7: Comparison between experimental (full lines) diffusion coefficient and theoretical (dots) ones. Each data set consists of diffusion coefficients of the anions (red) and the corresponding cations (blue). Circles correspond to data from the explicitly polarizable simulations, squares to the non-polarizable model generated by removing the polarizability from the CGPol model. Experimental results are from a fit to the Vogel-Fulcher-Tamman equation 139–141 in ref. 142. Errors in the experimental data are indicated by the shaded areas. Error bars for the computational data are not shown. Fitting errors for the computational data are smaller than the symbols. A conservative estimate for the relative error in the diffusion coefficients can be obtained from a couple of independent trajectories and is about 5 % (of the order of the size of the symbols).

Diffusion coefficients from the CGPol model are larger than experimental values by less

31



than a factor of two. The ordering between anion and cation diffusion is not reproduced with the CGPol force field. According to the experimental data, cations should diffuse faster than the anions. However all values are very close to each other and within the significant uncertainties of the experimental data. The only exception is at 450 K where the anion diffusion is faster even beyond the experimental error bars. The order of anion and cation diffusion coefficients is switched when using the non-polarizable force field, in qualitative agreement with the experiment. Diffusion coefficients in the non-polarizable force field simulations with full charges are lower than the ones in the explicitly polarizable simulations. These non-polarizable simulations are incidentally also closer to the experimental values. Removing the explicit polarizability has a larger effect on the anion diffusion than on the cation diffusion. It has been argued before, that inclusion of charge-scaling is actually necessary to reproduce the correct bulk diffusivity for these kinds of coarse-grained models. 79 This scaling is de facto not necessary for our parameter sets, which points to the fact that the used van-der-Waals parameters might play a larger role than anticipated. However, using not-explicitly polarizable force fields limits the applicability of the respective force fields to bulk solutions, or rather, periodic solutions with a homogeneous, fast component of the dielectric response.

Varying alkyl-chain length and PF6 – and BF4 – anions We want to validate our force field parameterization procedure further by generating force fields for various, chemically similar, species. Hence, we parameterized force fields for 1alkyl-3-methylimidazolium cations, with varying alkyl-chain length (ethyl, buytl, and hexyl) and two spherical anions of different size (PF6 – and BF4 – ). In the following, when we talk about changes in alkyl chain length, this implies changes in the larger alkyl-chain bead, meaning varying size, dispersion coefficient, and polarizability. All force field parameters are given in the SI. We start our comparison with RDFs as a structural property. Nanodomain formation has 32


Page 32 of 55

Page 33 of 55

been observed with increasing alkyl-chain length. 143 This domain formation is reflected in a stronger association of the long, apolar side chains. Figure 8 depicts the radial distribution functions of the long alkyl side-chain in our coarse-grained simulations. With increasing alkyl-chain length, there is increased association of the long alkyl chains.

2.5

C2C1Im C4C1Im C6C1Im

2.0

1.5 g(r)

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


1.0

0.5

0.0

5

10

15 ˚ r [A]

20

25

Figure 8: Site-site radial distribution functions g (r) of C2 (blue), C4 (green), and C6 (red) site in [C2 C1 Im][PF6 ], [C4 C1 Im][PF6 ], and [C6 C1 Im][PF6 ], respectively. The difference in association is clearly visible when going from ethyl to butyl side-chain. One can observe a shift to larger alkyl-chain distances due to the increased size of the butyl-side-chain, and an increased attraction due to the stronger dispersion interaction. The change in peak-height not as pronounced when changing from butyl to hexyl side chain. This already demonstrates a limit of our choice of coarse-grained representation, which might be

33



a bit too coarse to represent very large alkyl chains. We furthermore determined enthalpies of vaporization for all our model systems (see SI). These enthalpies are all in good agreement with available experimental data. Subtle trends with respect to alkyl-chain length and anions are, however, not well reproduced. These appear to require more detailed models for a full, quantitative description. 72,144,145 Much clearer trends with respect to alkyl-chain length are observed for dynamical properties. Diffusion coefficients for the various ionic liquids are listed in Tab. 4, together with available experimental values, and simulational data. Table 4: Diffusion coefficients for ionic liquids with varied alkyl-chain length in the cation, and different anions. Data is from our CGPol simulations, polarizable all-atom atomistic (AApol) simulations, 72 and available experimental data at 350 K. All data in 10−5 cm2 /s.

compound [C2 C1 Im][BF4 ] [C2 C1 Im][PF6 ] [C4 C1 Im][BF4 ] [C4 C1 Im][PF6 ] [C6 C1 Im][BF4 ] [C6 C1 Im][PF6 ]

CGPol AApol exp cation anion cation anion cation anion 0.18 0.17 0.155 0.102 0.2 146 0.18 146 0.24 0.19 0.060 0.035 142 0.058 0.073 0.113 0.093 0.105 0.11 142 142 0.079 0.086 0.059 0.041 0.07 0.056 142 0.040 0.057 0.051 0.054 0.051 0.064 0.029 0.024 -

All diffusion coefficients agree well with the experimental and theoretical reference values. The ordering between cation and anion is reproduced for [C2 C1 Im][BF4 ] and [C4 C1 Im][BF4 ]. However, differences in diffusion coefficients are rather small, and at least for [C2 C1 Im][BF4 ] within the error of the analysis. The force field simulations reproduce the slow-down with respect to increasing alkyl-chain length. Not reproduced is the order of diffusion coefficients with respect to the anion. For [C4 C1 Im][BF4 ] and [C4 C1 Im][PF6 ] the latter has larger diffusion coefficients, while the experiment predicts an opposite trend. As mentioned before, error bars in both the experimental, and computational data can be as large as to obfuscate the apparent ordering within the ion species. The wrong ordering between different anions can be attributed to their very coarse descriptions as single, polarizable spheres. Diffusion coefficients are related to the electrical conductivity through the Nernst-Einstein 34


Page 34 of 55

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


equation. However, the so-obtained conductivities do not contain any correlation effects between charge carriers. We calculate the full conductivities κ using the mean squared displacement between all ions, including all cross-correlations, as follows: n X n X 1 h(qi [Ri (t) − Ri (0)]) (qj [Rj (t) − Rj (0)])i . t→∞ 6tV kB T i j

κ = lim

(11)

Ri are the center-of-mass positions of the ions, qi the respective charges, V the volume of the unit cell, kB Boltzmann’s constant, t time, and T the temperature. Table 5 lists the conductivities of different ionic liquids with different alkyl chain lengths and compares them to experimental and other simulation data, from a recent study using explicitly polarizable, atomistic force fields. 72 Table 5: Conductivities for ionic liquids with different alkyl-chain lengths and different anions. Data from our CGPol simulations, polarizable all-atom atomistic (AApol) simulations, 72 and available experimental data. Errors in the computational conductivity calculations are large, up to 50 %. All data in S/m. species [C2 C1 Im][BF4 ] [C4 C1 Im][BF4 ] [C6 C1 Im][BF4 ] [C2 C1 Im][PF6 ] [C4 C1 Im][PF6 ] [C6 C1 Im][PF6 ]

CGPol AApol 72 4.1 5.65 1.7 2.30 0.4 0.81 6.9 1.40 2.8 0.89 0.9 0.49

expt 3.65, 4.35, 147 5.7 148 2.20, 142 2.17 148 1.00 148 147

1.24 142 0.64 149

The conductivities calculated with our CGPol force field agree well with the reference values. The CGPol force field captures the decrease in electrical conductivity with increasing alkyl-chain length. However, as for the diffusion coefficients, ordering of conductivities is not reproduced with respect to the different anions. Errors in the conductivity calculations can be large. However, the qualitative trends can be clearly observed in the full mean-square displacements, given in the SI. The CGPol force fields for the alkyl-imidazolium ionic liquids are currently constructed using specific, non-transferable parameter sets. However, charge distributions as well as 35



parameters for the van-der-Waals potential are fairly similar to each other. Hence, a simple averaging might yield reasonable parameters for here parameterized and other ionic liquids. We aim at exploring this and other approaches to create transferable parameters in future work.

4

Conclusions

In this article we have suggested a new method to consistently derive accurate polarizable, coarse-grained force fields from first principles. This procedure was demonstrated on a set of commonly used and widely investigated alkyl-imidazolium ionic liquids. Our approach has two main advantags: first, little computational effort is needed for the parameterization, and the methods employed are general enough to be applicable to a broader range of systems. Second, only one fitting parameter for the ion size was involved in the parameterization, in order to reproduce the correct mass density of the respective solutions. No other experimental data entered the parameterization. With our newly parameterized force fields, we achieved good correspondence to other simulational work and available experimental data. In particular, diffusion coefficients are well reproduced, as well as trends with respect to varying alkyl-chain length in the imidazolium cation. Furthermore, electronic properties like induced dipole moment and response to an applied electric field correspond well to ab initio results. One limiting factor in our current study is the focus on rigid, coarse-grained models. Hence, some details cannot be resolved within this representation, e.g., for different, spherical anions of similar size, or large alkyl chains. The force field parameterization procedure is, however, quite general and can also be applied to all-atom force field parameterizations that would inherently have those details built in. Our choice of dispersion coefficients is beneficial for the direct investigation of the influence of explicit inclusion of polarizability in simulations of ionic liquids. Without any additional effort, non-polarizable models appear readily available by simply removing the

36


Page 36 of 55

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


dipole-polarizability from our current force fields. Upon removal of the polarizability from the force field for [C4 C1 Im][PF6 ], we observed only little change in the mass density and in the self-diffusion coefficients. We plan to further investigate the influence of explicit polarizability and the construction of implicitly polarizable models in future work. The presented combination of explicitly polarizable force fields for coarse-grained systems offers direct access to accurate simulations of large systems on long time scales, such as investigations of ionic liquids in electrochemical applications. Because of the explicit treatment of polarizability, the newly parameterized force fields are straightforwardly usable in a broader range of system geometries, including porous interfaces, e.g., in supercapacitor investigations. The quick and straightforward parameterization of these force fields makes them good candidates in research of ionic liquids, where a high number of possible combinations of ions exists that form potentially relevant ionic liquids.

Acknowledgement FU thanks Jiří Kolafa for discussion of the ASPC implementation and providing data for its validation in LAMMPS. FU also thanks Jeanne Miriam Kohagen for careful reading of the manuscript and helpful suggestions. The authors gratefully acknowledge resources provided from bwHPC, DFG: INST 40/467-1 FUGG, and funding from DFG SFB 716 and the cluster of excellence “Simulation Technology” (EXC 310).

Supporting Information Available Supporting information is available free of charge on the ACS Publication website. Details on nomenclature and molecular structures. All parameter sets for the newly parameterized force fields. Sensitivity analysis of coarse-grained representation. Comparison of liquid structure between different force fields. Data on diffusion, enthalpy

37



of vaporization, and electrical conductivity. An example input file for the simulations.

References (1) 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. (2) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. (3) Scott, W. R. P.; Hünenberger, P. H.; Tironi, I. G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Krüger, P.; van Gunsteren, W. F. The GROMOS Biomolecular Simulation Program Package. J. Phys. Chem. A 1999, 103, 3596–3607. (4) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiórkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616, PMID: 24889800. (5) Grimme, S. A General Quantum Mechanically Derived Force Field (QMDFF) for Molecules and Condensed Phase Simulations. J. Chem. Theory Comput. 2014, 10, 4497–4514, PMID: 26588146. (6) Cole, D. J.; Vilseck, J. Z.; Tirado-Rives, J.; Payne, M. C.; Jorgensen, W. L. Biomolecu-

38


Page 38 of 55

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


lar Force Field Parameterization via Atoms-in-Molecule Electron Density Partitioning. J. Chem. Theory Comput. 2016, 12, 2312–2323. (7) McDaniel, J. G.; Schmidt, J. Physically-Motivated Force Fields from SymmetryAdapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053–2066, PMID: 23343200. (8) Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549–2564, PMID: 20136072. (9) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: New York, NY, USA, 1989. (10) Toukmaji, A.; Sagui, C.; Board, J.; Darden, T. Efficient particle-mesh Ewald based approach to fixed and induced dipolar interactions. J. Chem. Phys. 2000, 113, 10913– 10927. (11) Lamoureux, G.; Roux, B. Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm. J. Chem. Phys. 2003, 119, 3025–3039. (12) Kolafa, J. Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules. J. Comput. Chem. 2004, 25, 335–342. (13) Dequidt, A.; Devémy, J.; Pádua, A. A. H. Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator. J. Chem. Inf. Model. 2016, 56, 260–268. (14) Lemkul, J. A.; Roux, B.; van der Spoel, D.; MacKerell, A. D. Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model. J. Comput. Chem. 2015, 36, 1473–1479. 39



(15) Todorov, I. T.; Smith, W.; Trachenko, K.; Dove, M. T. DL_POLY_3: new dimensions in molecular dynamics simulations via massive parallelism. J. Mater. Chem. 2006, 16, 1911–1918. (16) Leontyev, I. V.; Stuchebrukhov, A. A. Electronic continuum model for molecular dynamics simulations. J. Chem. Phys. 2009, 130, 085102. (17) Leontyev, I.; Stuchebrukhov, A. Accounting for electronic polarization in nonpolarizable force fields. Phys. Chem. Chem. Phys. 2011, 13, 2613–2626. (18) Lyubartsev, A. P.; Laaksonen, A. Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach. Phys. Rev. E 1995, 52, 3730–3737. (19) Meyer, H.; Biermann, O.; Faller, R.; Reith, D.; Müller-Plathe, F. Coarse graining of nonbonded inter-particle potentials using automatic simplex optimization to fit structural properties. J. Chem. Phys. 2000, 113, 6264–6275. (20) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. J. Comput. Chem. 2003, 24, 1624–1636. (21) Izvekov, S.; Voth, G. A. A Multiscale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109, 2469–2473, PMID: 16851243. (22) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812–7824, PMID: 17569554. (23) Michalowsky, J.; Schäfer, L. V.; Holm, C.; Smiatek, J. A refined polarizable water model for the coarse-grained MARTINI force field with long-range electrostatic interactions. J. Chem. Phys. 2017, 146, 054501.

40


Page 40 of 55

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


(24) Zeman, J.; Uhlig, F.; Smiatek, J.; Holm, C. A coarse-grained polarizable force field for the ionic liquid 1-butyl-3-methylimidazolium hexafluorophosphate. J. Phys.: Condens. Matter 2017, 29, 504004. (25) 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. (26) Paek, E.; Pak, A. J.; Hwang, G. S. On the influence of polarization effects in predicting the interfacial structure and capacitance of graphene-like electrodes in ionic liquids. J. Chem. Phys. 2015, 142 . (27) Salanne, M. Simulations of room temperature ionic liquids: from polarizable to coarsegrained force fields. Phys. Chem. Chem. Phys. 2015, 17, 14270–14279. (28) Yasuda, T.; Watanabe, M. Protic ionic liquids: Fuel cell applications. MRS Bulleti 2013, 38, 560–566. (29) Watanabe, M.; Thomas, M. L.; Zhang, S.; Ueno, K.; Yasuda, T.; Dokko, K. Application of Ionic Liquids to Energy Storage and Conversion Materials and Devices. Chem. Rev. 2017, 117, 7190–7239, PMID: 28084733. (30) Dai, C.; Zhang, J.; Huang, C.; Lei, Z. Ionic Liquids in Selective Oxidation: Catalysts and Solvents. Chem. Rev. 2017, 117, 6929–6983, PMID: 28459547. (31) Hayes, R.; Warr, G. G.; Atkin, R. Structure and Nanostructure in Ionic Liquids. Chem. Rev. 2015, 115, 6357–6426. (32) Sato, T.; Masuda, G.; Takagi, K. Electrochemical properties of novel ionic liquids for electric double layer capacitor applications. Electrochim. Acta 2004, 49, 3603 – 3611.

41



(33) Armand, M.; Endres, F.; MacFarlane, D. R.; Ohno, H.; Scrosati, B. Ionic-liquid materials for the electrochemical challenges of the future. Nat. Mater. 2009, 8, 621 – 629. (34) Breitsprecher, K.; Košovan, P.; Holm, C. Coarse-grained simulations of an ionic liquidbased capacitor: I. Density, ion size, and valency effects. J. Phys.: Condens. Matter 2014, 26, 284108. (35) Dong, K.; Liu, X.; Dong, H.; Zhang, X.; Zhang, S. Multiscale Studies on Ionic Liquids. Chem. Rev. 2017, 117, 6636–6695, PMID: 28488441. (36) Gabl, S.; Schröder, C.; Steinhauser, O. Computational studies of ionic liquids: Size does matter and time too. J. Chem. Phys. 2012, 137, 094501. (37) Lesch, V.;

Heuer, A.;

Holm, C.;

Smiatek, J. Solvent effects of 1-ethyl-3-

methylimidazolium acetate: solvation and dynamic behavior of polar and apolar solutes. Phys. Chem. Chem. Phys. 2015, 17, 8480–8490. (38) 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. (39) Canongia Lopes, J. N.; Pádua, A. A. H. Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. J. Phys. Chem. B 2004, 108, 16893–16898. (40) Canongia Lopes, J. N.; Pádua, A. A. H. Molecular Force Field for Ionic Liquids III:âĂĽ Imidazolium, Pyridinium, and Phosphonium Cations; Chloride, Bromide, and Dicyanamide Anions. J. Phys. Chem. B 2006, 110, 19586–19592, PMID: 17004824. (41) Canongia Lopes, J. N.; Pádua, A. A. H.; Shimizu, K. Molecular Force Field for Ionic Liquids IV:âĂĽ Trialkylimidazolium and Alkoxycarbonyl-Imidazolium Cations; Alkylsulfonate and Alkylsulfate Anions. J. Phys. Chem. B 2008, 112, 5039–5046, PMID: 18380506. 42


Page 42 of 55

Page 43 of 55 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


(42) 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, PMID: 26609613. (43) Chaban, V. Polarizability versus mobility: atomistic force field for ionic liquids. Phys. Chem. Chem. Phys. 2011, 13, 16055–16062. (44) Chaban, V. V.; Prezhdo, O. V. A new force field model of 1-butyl-3-methylimidazolium tetrafluoroborate ionic liquid and acetonitrile mixtures. Phys. Chem. Chem. Phys. 2011, 13, 19345–19354. (45) 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. (46) Shah, J. K.; Brennecke, J. F.; Maginn, E. J. Thermodynamic properties of the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate from Monte Carlo simulations. Green Chem. 2002, 4, 112–118. (47) Cadena, C.; Maginn, E. J. Molecular Simulation Study of Some Thermophysical and Transport Properties of Triazolium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110, 18026–18039, PMID: 16956294. (48) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. Molecular Modeling and Experimental Studies of the Thermodynamic and Transport Properties of Pyridinium-Based Ionic Liquids. J. Phys. Chem. B 2006, 110, 2821–2832, PMID: 16471891. (49) Maginn, E. J. Atomistic Simulation of the Thermodynamic and Transport Properties of Ionic Liquids. Acc. Chem. Res. 2007, 40, 1200–1207, PMID: 17953449. (50) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463–11478, PMID: 19637900. 43



(51) Starovoytov, O. N.; Torabifard, H.; Cisneros, G. A. Development of AMOEBA Force Field for 1,3-Dimethylimidazolium Based Ionic Liquids. J. Phys. Chem. B 2014, 118, 7156–7166, PMID: 24901255. (52) 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. (53) Lee, A. A.; Kondrat, S.; Vella, D.; Goriely, A. Dynamics of Ion Transport in Ionic Liquids. Phys. Rev. Lett. 2015, 115, 106101. (54) Krekeler, C.; Dommert, F.; Schmidt, J.; Zhao, Y.; Holm, C.; Berger, R.; Delle Site, L. Electrostatic properties of liquid 1,3-dimethylimidazolium chloride: role of local polarization and effect of the bulk. Phys. Chem. Chem. Phys. 2010, 12, 1817–1821. (55) Wendler, K.; Zahn, S.; Dommert, F.; Berger, R.; Holm, C.; Kirchner, B.; Site, L. D. Locality and Fluctuations: Trends in Imidazolium-Based Ionic Liquids and Beyond. J. Chem. Theory Comput. 2011, 7, 3040–3044. (56) Wendler, K.; Dommert, F.; Zhao, Y.; Berger, R.; Holm, C.; Delle Site, L. Ionic liquids studied across different scales: A computational perspective. Faraday Discuss. 2012, 154, 111–132. (57) 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. (58) Schmidt, J.; Krekeler, C.; Dommert, F.; Zhao, Y.; Berger, R.; Delle Site, L.; Holm, C. Ionic Charge Reduction and Atomic Partial Charges from First-Principles Calculations of 1,3-Dimethylimidazolium Chloride. J. Phys. Chem. B 2010, 114, 6150–6155.

44


Page 44 of 55

Page 45 of 55 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


(59) Dommert, F.; Wendler, K.; Berger, R.; Site, L. D.; Holm, C. Force Fields for Studying the Structure and Dynamics of Ionic Liquids: A Critical Review of Recent Developments. ChemPhysChem 2012, 13, 1625–1637. (60) Dommert, F.; Holm, C. Refining classical force fields for ionic liquids: theory and application to [MMIM][Cl]. Phys. Chem. Chem. Phys. 2013, 15, 2037–2049. (61) Dommert, F.; Wendler, K.; Qiao, B.; Site, L. D.; Holm, C. Generic force fields for ionic liquids. J. Mol. Liq. 2014, 192, 32 – 37, Fundamental Aspects of Ionic Liquid Science. (62) Ángyán, J. G.; Jansen, G.; Loss, M.; Hättig, C.; Heß, B. A. Distributed polarizabilities using the topological theory of atoms in molecules. Chem. Phys. Lett. 1994, 219, 267 – 273. (63) Elking, D. M.; Perera, L.; Duke, R.; Darden, T.; Pedersen, L. G. A finite field method for calculating molecular polarizability tensors for arbitrary multipole rank. J. Comput. Chem. 2011, 32, 3283–3295. (64) Millot, C.; Chaumont, A.; Engler, E.; Wipff, G. Distributed Polarizability Models for Imidazolium-Based Ionic Liquids. J. Phys. Chem. A 2014, 118, 8842–8851, PMID: 25133873. (65) Bica, K.; Deetlefs, M.; Schröder, C.; Seddon, K. R. Polarisabilities of alkylimidazolium ionic liquids. Phys. Chem. Chem. Phys. 2013, 15, 2703–2711. (66) Bernardes, C. E. S.; Shimizu, K.; Lopes, J. N. C.; Marquetand, P.; Heid, E.; Steinhauser, O.; Schroder, C. Additive polarizabilities in ionic liquids. Phys. Chem. Chem. Phys. 2016, 18, 1665–1670. (67) Miller, K. J.; Savchik, J. A new empirical method to calculate average molecular polarizabilities. J. Am. Chem. Soc. 1979, 101, 7206–7213. 45



(68) Thole, B. Molecular polarizabilities calculated with a modified dipole interaction. Chem. Phys. 1981, 59, 341 – 350. (69) van Duijnen, P. T.; Swart, M. Molecular and Atomic Polarizabilities: Thole’s Model Revisited. J. Phys. Chem. A 1998, 102, 2399–2407. (70) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (71) Son, C. Y.; McDaniel, J. G.; Schmidt, J. R.; Cui, Q.; Yethiraj, A. First-Principles United Atom Force Field for the Ionic Liquid BMIM+ BF− 4 : An Alternative to Charge Scaling. J. Phys. Chem. B 2016, 120, 3560–3568. (72) 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. (73) Salanne, M.; Rotenberg, B.; Jahn, S.; Vuilleumier, R.; Simon, C.; Madden, P. A. Including many-body effects in models for ionic liquids. Theor. Chem. Acc. 2012, 131, 1143. (74) Lindorff-Larsen, K.; Maragakis, P.; Piana, S.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. Systematic Validation of Protein Force Fields against Experimental Data. PLoS One 2012, 7, 1–6. (75) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern

46


Page 46 of 55

Page 47 of 55 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


Free-Energy Calculation Protocol and Force Field. J. Am. Chem. Soc. 2015, 137, 2695–2703, PMID: 25625324. (76) Wang, Y.; Izvekov, S.; Yan, T.; Voth, G. A. Multiscale Coarse-Graining of Ionic Liquids. J. Phys. Chem. B 2006, 110, 3564–3575. (77) Bhargava, B. L.; Devane, R.; Klein, M. L.; Balasubramanian, S. Nanoscale organization in room temperature ionic liquids: a coarse grained molecular dynamics simulation study. Soft Matter 2007, 3, 1395–1400. (78) Roy, D.; Patel, N.; Conte, S.; Maroncelli, M. Dynamics in an Idealized Ionic Liquid Model. J. Phys. Chem. B 2010, 114, 8410–8424. (79) Roy, D.; Maroncelli, M. An Improved Four-Site Ionic Liquid Model. J. Phys. Chem. B 2010, 114, 12629–12631. (80) Merlet, C.; Salanne, M.; Rotenberg, B. New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations. J. Phys. Chem. C 2012, 116, 7687–7693. (81) Wang, Y.; Feng, S.; Voth, G. A. Transferable Coarse-Grained Models for Ionic Liquids. J. Chem. Theory Comput. 2009, 5, 1091–1098, PMID: 26609619. (82) Noskov, S. Y.; Lamoureux, G.; Roux, B. Molecular Dynamics Study of Hydration in Ethanol-Water Mixtures Using a Polarizable Force Field. J. Phys. Chem. B 2005, 109, 6705–6713. (83) Blöchl, P. E. Electrostatic decoupling of periodic images of plane-wave-expanded densities and derived atomic point charges. J. Chem. Phys. 1995, 103, 7422–7428. (84) Blair, S. A.; Thakkar, A. J. Relating polarizability to volume, ionization energy, electronegativity, hardness, moments of momentum, and other molecular properties. J. Chem. Phys. 2014, 141 . 47



(85) Brinck, T.; Murray, J. S.; Politzer, P. Polarizability and volume. J. Chem. Phys. 1993, 98, 4305–4306. (86) Laidig, K. E.; Bader, R. F. W. Properties of atoms in molecules: Atomic polarizabilities. J. Chem. Phys. 1990, 93, 7213–7224. (87) Yu, M.; Trinkle, D. R. Accurate and efficient algorithm for Bader charge integration. J. Chem. Phys. 2011, 134, 064111. (88) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. (89) Schmollngruber, M.; Lesch, V.; Schröder, C.; Heuer, A.; Steinhauser, O. Comparing induced point-dipoles and Drude oscillators. Phys. Chem. Chem. Phys. 2015, 17, 14297–14306. (90) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 116, 4983–5013. (91) Hamaker, H. The London-van der Waals attraction between spherical particles. Physica 1937, 4, 1058 – 1072. (92) Margenau, H. Note on the Calculation of van der Waals Forces. Phys. Rev. 1931, 37, 1425–1430. (93) Zeiss, G.; Meath, W. J. Dispersion energy constants C 6(A, B), dipole oscillator strength sums and refractivities for Li, N, O, H2, N2, O2, NH3, H2O, NO and N2O. Mol. Phys. 1977, 33, 1155–1176. (94) London, F. The general theory of molecular forces. Trans. Faraday Soc. 1937, 33, 8b–26.

48


Page 48 of 55

Page 49 of 55 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) Casimir, H. B. G.; Polder, D. The Influence of Retardation on the London-van der Waals Forces. Phys. Rev. 1948, 73, 360–372. (96) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Science 1938, 168, 264–283. (97) Neese, F. The ORCA program system. Wiley Interdiscip Rev Comput Mol Sci 2012, 2, 73–78. (98) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k: atomistic simulations of condensed matter systems. Wiley Interdiscip Rev Comput Mol Sci 2014, 4, 15–25. (99) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103 – 128. (100) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (101) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. (102) Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple”. Phys. Rev. Lett. 1998, 80, 890–890. (103) Putrino, A.; Sebastiani, D.; Parrinello, M. Generalized variational density functional perturbation theory. J. Chem. Phys. 2000, 113, 7102–7109. (104) Luber, S.; Iannuzzi, M.; Hutter, J. Raman spectra from ab initio molecular dynamics and its application to liquid S-methyloxirane. J. Chem. Phys. 2014, 141, 094503. (105) Jr., T. H. D. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. 49



(106) Woon, D. E.; Jr., T. H. D. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358–1371. (107) Wannier, G. H. The Structure of Electronic Excitation Levels in Insulating Crystals. Phys. Rev. 1937, 52, 191–197. (108) Resta, R.; Vanderbilt, D. Physics of Ferroelectrics: A Modern Perspective; Springer Berlin Heidelberg: Berlin, Heidelberg, 2007; pp 31–68. (109) Marzari, N.; Mostofi, A. A.; Yates, J. R.; Souza, I.; Vanderbilt, D. Maximally localized Wannier functions: Theory and applications. Rev. Mod. Phys. 2012, 84, 1419–1475. (110) Marques, M. A.; Castro, A.; Bertsch, G. F.; Rubio, A. octopus: a first-principles tool for excited electron-ion dynamics. Comput. Phys. Commun. 2003, 151, 60 – 78. (111) Marques, M. A. L.; Castro, A.; Malloci, G.; Mulas, G.; Botti, S. Efficient calculation of van der Waals dispersion coefficients with time-dependent density functional theory in real time: Application to polycyclic aromatic hydrocarbons. J. Chem. Phys. 2007, 127 . (112) Legrand, C.; Suraud, E.; Reinhard, P.-G. Comparison of self-interaction-corrections for metal clusters. J. Phys. B: At., Mol. Opt. Phys. 2002, 35, 1115. (113) Chu, X.; Dalgarno, A. Linear response time-dependent density functional theory for van der Waals coefficients. J. Chem. Phys. 2004, 121, 4083–4088. (114) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. Phys. Rev. B 1998, 58, 3641–3662. (115) Canongia Lopes, J. N.; Pádua, A. A. H. CL&P: A generic and systematic force field for ionic liquids modeling. Theor. Chem. Acc. 2012, 131, 1129. (116) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys., 1995, 117, 1 – 19. 50


Page 50 of 55

Page 51 of 55 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


(117) Padua, A. fftool v1.0.0. 2015; https://doi.org/10.5281/zenodo.18618. (118) Hockney, R. W.; Eastwood, J. W. Particle-Particle-Particle-Mesh (P3M) Algorithms; Taylor & Francis, 1988. (119) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied Physics 1981, 52, 7182–7190. (120) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177–4189. (121) Shinoda, W.; Shiga, M.; Mikami, M. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Phys. Rev. B 2004, 69, 134103. (122) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys., 1977, 23, 327 – 341. (123) Andersen, H. C. Rattle: A "velocity" version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys., 1983, 52, 24 – 34. (124) Uhlig, F. Implementation of ASPC algorithm for LAMMPS. https://github.com/ fuulish/ASPC.git, 2016–2017. (125) 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, 2157–2164. (126) Qiao, Y.; Yan, F.; Xia, S.; Yin, S.; Ma, P. Densities and Viscosities of [Bmim][PF6] and Binary Systems [Bmim][PF6] + Ethanol, [Bmim][PF6] + Benzene at Several Temperatures and Pressures: Determined by the Falling-Ball Method. J. Chem. Eng. Data 2011, 56, 2379–2385.

51



(127) Harris, K. R.; Woolf, L. A.; Kanakubo, M. Temperature and Pressure Dependence of the Viscosity of the Ionic Liquid 1-Butyl-3-methylimidazolium Hexafluorophosphate. J. Chem. Eng. Data 2005, 50, 1777–1782. (128) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (129) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (130) Hunt, P. A.;

Gould, I. R. Structural Characterization of the 1-Butyl-3-

methylimidazolium Chloride Ion Pair Using ab Initio Methods. J. Phys. Chem. A 2006, 110, 2269–2282, PMID: 16466265. (131) Koßmann, S.; Thar, J.; Kirchner, B.; Hunt, P. A.; Welton, T. Cooperativity in ionic liquids. J. Chem. Phys. 2006, 124, 174506. (132) Huddleston, J. G.; Visser, A. E.; Reichert, W. M.; Willauer, H. D.; Broker, G. A.; Rogers, R. D. Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation. Green Chem. 2001, 3, 156–164. (133) 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. (134) Chaban, V. V.; Prezhdo, O. V. Ionic Vapor: What Does It Consist Of? J. Phys. Chem. Lett. 2012, 3, 1657–1662. (135) Neto, B. A. D.; Meurer, E. C.; Galaverna, R.; Bythell, B. J.; Dupont, J.; Cooks, R. G.;

52


Page 52 of 55

Page 53 of 55 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


Eberlin, M. N. Vapors from Ionic Liquids: Reconciling Simulations with Mass Spectrometric Data. J. Phys. Chem. Lett. 2012, 3, 3435–3441. (136) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269–6271. (137) Zaitsau, D. H.; Kabo, G. J.; Strechan, A. A.; Paulechka, Y. U.; Tschersich, A.; Verevkin, S. P.; Heintz, A. Experimental Vapor Pressures of 1-Alkyl-3methylimidazolium Bis(trifluoromethylsulfonyl)imides and a Correlation Scheme for Estimation of Vaporization Enthalpies of Ionic Liquids. J. Phys. Chem. A 2006, 110, 7303–7306, PMID: 16737284. (138) 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. (139) Vogel, H. The law of the relation between the viscosity of liquids and the temperature. Physikalische Zeitschrift 1921, 22, 645. (140) Fulcher, G. S. Analysis of recent measurements of the viscosity of glasses. J. Am. Ceram. Soc. 1925, 8, 339–355. (141) Tammann, G.; Hesse, W. Die Abhängigkeit der Viskosität von der Temperatur bei unterkühlten Flüssigkeiten. Z. Anorg. Allg. Chem. 1926, 156, 245–257. (142) 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. (143) Canongia Lopes, J. N. A.; Pádua, A. A. H. Nanostructural Organization in Ionic Liquids. J. Phys. Chem. B 2006, 110, 3330–3335.

53



(144) Köddermann, T.; Paschek, D.; Ludwig, R. Ionic Liquids: Dissecting the Enthalpies of Vaporization. ChemPhysChem 2008, 9, 549–555. (145) Zaitsau, D. H.; Fumino, K.; Emel’yanenko, V. N.; Yermalayeu, A. V.; Ludwig, R.; Verevkin, S. P. Structure-Property Relationships in Ionic Liquids: A Study of the Anion Dependence in Vaporization Enthalpies of Imidazolium-Based Ionic Liquids. ChemPhysChem 2012, 13, 1868–1876. (146) Noda, A.; Hayamizu, K.; Watanabe, M. Pulsed-Gradient SpinâĹŠEcho 1H and 19F NMR Ionic Diffusion Coefficient, Viscosity, and Ionic Conductivity of NonChloroaluminate Room-Temperature Ionic Liquids. J. Phys. Chem. B 2001, 105, 4603–4610. (147) Zhang, S.; Sun, N.; He, X.; Lu, X.; Zhang, X. Physical Properties of Ionic Liquids: Database and Evaluation. J. Phys. Chem. Ref. Data 2006, 35, 1475–1517. (148) Stoppa, A.; Zech, O.; Kunz, W.; Buchner, R. The Conductivity of Imidazolium-Based Ionic Liquids from (-35 to 195) ◦ C. A. Variation of Cation’s Alkyl Chain. J. Chem. Eng. Data 2010, 55, 1768–1773. (149) Ning, H.; Hou, M.; Mei, Q.; Liu, Y.; Yang, D.; Han, B. The physicochemical properties of some imidazolium-based ionic liquids and their binary mixtures. Science China Chemistry 2012, 55, 1509–1518.

54


Page 54 of 55

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


Graphical TOC Entry

55


First-principles parameterization of polarizable coarse-grained force

Recommend Documents