The effect of various quantum mechanically derived partial atomic charges on the bulk properties of chloride-based ionic liquids

Accepted Manuscript The Effect of Various Quantum Mechanically Derived Partial Atomic Charges on the Bulk Properties of Chloride-Based Ionic Liquids Amin Reza Zolghadr, Mohammad Hadi Ghatee, Fatemeh Moosavi PII: DOI: Reference:

S0301-0104(16)30293-2 http://dx.doi.org/10.1016/j.chemphys.2016.05.022 CHEMPH 9566

To appear in:

Chemical Physics

Received Date: Accepted Date:

4 April 2016 25 May 2016

Please cite this article as: A.R. Zolghadr, M.H. Ghatee, F. Moosavi, The Effect of Various Quantum Mechanically Derived Partial Atomic Charges on the Bulk Properties of Chloride-Based Ionic Liquids, Chemical Physics (2016), doi: http://dx.doi.org/10.1016/j.chemphys.2016.05.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

The Effect of Various Quantum Mechanically Derived Partial Atomic Charges on the Bulk Properties of Chloride-Based Ionic Liquids Amin Reza Zolghadr,* ,1 Mohammad Hadi Ghatee,1 and Fatemeh Moosavi2 1

2

Department of Chemistry, Shiraz University, Shiraz, 71946-84795, Iran

Department of Chemistry, Ferdowsi University of Mashhad, Mashhad 91779, Iran

*Corresponding author Fax: +98 71 3646 0788. Tel: +98 71 3613 7157. E-mail address: [email protected]. 1

ABSTRACT Partial atomic charges using various quantum mechanical calculations for [Cnmim]Cl (

) ionic

liquids (ILs) are obtained and used for development of molecular dynamics simulation (MD) force fields. The isolated ion pairs are optimized using HF, B3LYP, and MP2 methods for electronic structure with 6311++G(d,p) basis set. Partial atomic charges are assigned to the atomic center with CHELPG and NBO methods. The effect of these partial charge sets on the static and dynamic properties of ILs is evaluated by performing a series of MD simulations and comparing the essential thermodynamic properties with the available experimental data and available molecular dynamics simulation results. In contrast to the general trends reported for ionic liquids with BF4, PF6, and iodide anions (in which restrained electrostatic potential (RESP) charges are preferred), partial charges derived by B3LYP-NBO method are relatively good in prediction of the structural, dynamical, and thermodynamic energetic properties of the chloride based ILs. Keywords: Partial charges, Molecular dynamics simulations, chloride ionic liquids, force fields

2

INTRODUCTION Room-temperature ionic liquids (ILs) define a new class of condensed matters because of their unique physicochemical properties. A molecular-based understanding of ILs poses a great challenge because the partial charge, the geometrical structure, and the electronic structure of the ions give rise to a complex interplay of molecular interactions.[1] Partial atomic charges, unlike the electron density, are not quantum mechanically observable, that is, they cannot be calculated from first principles. Hence, all methods for derivation of charge are ultimately arbitrary. They are, however, comparable based on the electrostatic performance of partial atomic charges produced for a given molecular system. Nevertheless, there is no universally accepted procedure, as the concept of partial atomic charge remains still arbitrary.[2] However, partial atomic charges are fundamental parameters for developing force fields in simulations and their accuracies play a substantial role in simulations of condensed phase properties. Under such circumstances, there is a significant debate about the best deal to determine effective atomic charges to simulate physical properties.[3] The most widely used methods for estimating atomic charges are derived from a least-square fit to the electrostatic potential (ESP).[4] One of these methods, CHELP (CHarges from ELectrostatic Potentials), was initially developed by Chirlian and Francl, and then modified by Breneman and Wiberg as the grid-oriented CHELPG method which is less dependent upon molecular orientation than the original CHELP method in which partial atomic charges are fitted to reproduce the molecular ESP at a number of points around the molecule. [5,6] ESP is obviously one of the most useful properties in order to acquire proper partial atomic charges to model shortrange as well as long-range molecule−molecule interactions. [7] On the other hand, natural bond orbital (NBO) analysis provides an accurate Lewis structure picture and charge density of a molecule using the highest percentage of orbital electron density representation. [8]

3

Generating the charges using ab initio methods may give the most accurate representation of the charge distribution, provided that a suitable atomic basis set is used. The values of the atomic charges obtained by this method can be successively improved as the representation of the basis set improves. The reason is that the total molecular charge distribution depends on the nuclear positions and the agreement of the optimized geometries with the experimental ones will vary with the basis set applied. Although the ab initio derived charges fluctuate significantly with small basis sets, after one reaches a basis set of 6-3lG* quality, the electrostatic potential is close to convergent with respect to improvements in the basis set. [3] Development of force fields that provide satisfying description of both static and dynamic properties of the system is up to now a major open problem. Some molecular dynamics simulations have used force fields with reduced ionic charges to achieve this aim. [9] Morrow and Maginn, in their study on 1-butyl-3-methylimidazolium hexafluorophosphate ([C4mim]PF6) were used a total charge of ±0.904 e on the basis of the quantum calculation of an ion pair in a vacuum. [10] Bhargava and Balasubramanian have explained that charge reduction leads to a faster dynamics for ionic liquids. [11] Holm and co-workers have presented detailed study of force field parameterization for 1,3-dimethylimidazolium chloride ionic liquid. [12,13] They performed MP2 electronic structure and DFT calculations on isolated ion pairs as well as Car-Parrinello MD simulation of 30 ion pairs. Their results indicated that absolute value of total ionic charge on cation and the anion are considerably less than unity. [14] Very recently, Mondal and Balasubramanian were carried out classical MD simulations of various ionic liquids with refined site charges obtained from condensed phase DFT calculations. [15,16] They indicated that calculated density, heat of vaporization, surface tension, diffusion coefficient and collective transport properties are in nearly quantitative agreement with experimental data. [17] Furthermore, we have carried out MD simulations with different atomic charges from various

4

models in order to clarify the structure and dynamics of the pure ILs, mixtures of ILs with water, and alkane/ILs/water interfaces. [18-20] The CHELPG and NBO methods represent reasonable atom-centered charges for use in molecular simulations. However, there are still some issues that remain to be explored in choosing most appropriate partial atomic charges for simulation of periodic systems. A major aim of this work is to propose the most efficient method that gives rise to the most reliable charge values for the simulation of bulk properties of [C1mim]Cl and [C4mim]Cl ionic liquids. The paper is organized as follows: this introduction is followed by details of the methodology employed to calculate partial atomic charges, force field parameters and details of simulations. The third section is devoted to the description of results obtained and comparison with literature, and followed by conclusions section. COMPUTATIONAL METHODS Ab initio calculations The structure of each IL was optimized at different level of theory: Hartree–Fock (HF) [21], hybrid density functional method (B3LYP) [22], and Möller–Plesset perturbation theory to the second-order (MP2) [23] using Gaussian 03 program [24], with 6-311++G(d,p) basis set. The basis set is selected large enough for the results to be independent of basis set. The computed structures of neutral ILs were checked for vibrational frequencies to assure the state of minimum energies is attained on the potential energy surface [25]. The structures of [C1mim]Cl and [C4mim]Cl ILs optimized with HF, B3LYP, and MP2 approaches are depicted in Figs. S1 and S2. (see Fig. 1 for atomic labels) Then, new sets of partial charges are assigned (to each IL) by two different procedures: (i) fitting to the electrostatic potential surface with the CHELPG procedure and (ii) carrying out population analysis using the NBO method, implemented in Gaussian 03. Force field parameters 5

For the cationic part of ILs, the force field is in the form of explicit fully flexible all-atom force field developed by Canongia Lopes et al. (CLaP) [26,27]. The electrostatic charges for dialkylimidazolium salts have been taken from our ab initio calculations (see Tables 1 and 2). The slow dynamics of ILs and lack of attaining true equilibrium essentially lead to a large deviation on liquid state properties especially on simulated viscosity and surface tension. Polarizable force field increases the diffusion constant by a factor of three compared to the nonpolarizable one [28]. On the other hand, due to high computational demands for the simulation of such complex liquids with ionic nature that dispersion interactions and hydrogen bonding play important roles, non-polarizable force fields are preferred [29]. The total interaction potential has the form V tot = ∑

bonds

3 V kb k ( r − req ) 2 + ∑ θ (θ −θ eq ) 2 + ∑ ∑ i [1 + (-1) i -1 cos( iΦ )] + 2 i 2 dihedrals i =1 2

  σ  ij ∑ ∑ 4ε ij   rij i =1 j 〉1  

N −1 N

12

   

 σ ij −  rij 

   

(1)

 qq  + i j   4πε
rij   

6

where all terms have their usual meaning. The values of the parameters of eq. 1 were taken from references [27] and [30]. The cross term parameters ε ij and σ ij are given by combining and mixing rules, ε ij = (ε ii ε jj )1 / 2 and σ ij = 1 / 2(σ ii + σ jj ) , respectively. Simulation details For each liquid system, [C1mim]Cl and [C4mim]Cl, we performed simulation studies at ambient pressure ( 1 .01325

× 10 5 Pa

) using the DL_POLY program version 2.17 [31]. The equations of

motion were solved using Verlet-Leapfrog integration algorithm, under the periodic boundary conditions. The Columbic long-range interactions were calculated using Ewald’s method with a precision of 1 × 10 −5 . The potential cut-off distance value of 16 Å was used for each simulation. For bulk simulations, a single ion pair with a geometrical structure optimized via ab initio

6

calculation was replicated to obtain an ensemble containing 512 ion pairs of each liquid system under study. The simulation had to begin with short time step between 1×10-4 and 5×10-4 ps for 500 ps, and then followed by time step of 1×10-3 ps. After the initial equilibration with short time step, the temperature of the system was increased at intervals of 100 K up to 800 K. At each temperature, the simulation was performed for 2ns under constant NPT conditions using the Nosé–Hoover thermostat/barostat algorithm [32,33] and the modification of Melchionna et al. [34], as implemented in the DL_POLY program. Then, the temperature of the system was decreased to 450 K at steps of 50 K. The relaxation times for the thermostat and barostat are 0.1 and 2.0 ps, respectively. After all these preliminary simulation adjustments, the simulations were continued at this temperature for 8 ns in the canonical (NVT) ensemble. The total energy of system indicates the state of equilibrium is reached satisfactorily. The obtained trajectory was considered equilibration and discarded. The NVT simulation was followed by 8 more ns collecting statistical data. The same simulation procedure is applied for each IL system while the quantum mechanically derived atomic charges were replaced in the force field. RESULTS AND DISCUSSION Atomic charges The six sets of atomic charges derived in the current work are shown in Tables 1 and 2 for [C1mim]Cl and [C4mim]Cl, respectively. Atomic labels follow those in Fig. 1. The sum of partial charges of cation obtained in this work are considerably less than unity in agreement with literature [12]. The absolute value of anion and cation charges are exactly the same in all cases. The amounts of total charges for different groups of atoms, as previously explored in literature [35,36], are shown in Tables S1 and S2 for [C1mim]Cl and [C4mim]Cl, respectively. While comparing the results calculated with different methods, it is interesting to note that the largest absolute value of total charge ( 0.939 e) was obtained by HF-NBO and the smallest value 7

( 0.719 e) was obtained by B3LYP-CHELPG for [C1mim]Cl (see Table 1 and Fig. S3). In the other words, the values of total charges obtained by B3LYP method are more distant from the unity (formal charge

) and the largest deviation was observed for CHELPG charges. The

same conclusion is true for the case of [C4mim]Cl (see Table 2 and Fig. S4). The partial charges of an ion pair determined for a single isolated [C1mim]Cl ion pair are visualized and compared in Fig. S3. The trend of partial charge variation is almost the same in all cases. Nitrogen atoms of the imidazolium ring are negatively charged and the CR carbon atom of the imidazolium ring has positive charge. Other carbon atoms are mostly negatively charged. All hydrogen atoms have positive charge. For [C4mim]Cl, the nitrogen atoms of the imidazolium ring possess negative charges and C2 atom has positive charge by the results of NBO calculations, contrary to CHELPG method. In addition, CR and H R atoms of the imidazolium ring have positive charges in line with [C1mim]Cl results. For both ion pairs, no striking differences in HR charge by different methods are noticed but the charge on chloride anion differs remarkably. Therefore, the larger negative charge on chloride could extend the hydrogen-bond interaction between anion and cation. As shown in Tables S1 and S2, each cation is divided into three different groups: imidazolium ring, methyl group(s), and butyl chain. According to the charge values of each group, although the positive charge of the cation (by NBO) is mostly localized on the alkyl chain, different trends are observed for CHELPG charges. Much smaller positive charges are observed for the alkyl chains specially in the case of [C4mim]Cl (see Fig. 2).

Temperature dependence of density A measuring utility for the atomic charge definition can be established by questioning how well it reproduces the macroscopic properties. One common quantity analyzed is the liquid density. 8

Here, the force fields involving atomic charges derived from quantum vacuum phase calculations and any scaling factor to the atomic charges were applied. The temperature dependence of liquid densities for both ionic liquids are simulated using the six sets of atomic charges in the temperature range 273-450 K. For [C1mim]Cl, MP2-NBO and HFNBO charges predict the highest densities among others, in good agreement with the experimentally measured densities of Fannin et al. [37], and with the MD simulation by Bhargava and Balasubramaninan [38]. The results based on Gaussian distributed multipole analysis (GDMA) and Gaussian electrostatic model-distributed multipole (GEM-DM) force fields are also shown in Fig. 3A for the sake of comparison [39]. The maximum deviation is observed for B3LYP-CHELPG results (see Fig. 3A). In the case of [C4mim]Cl, all six sets of partial charges generate slightly lower densities than the experimental data of Ghatee et al., which were measured after extracting and removing water content and dissolved atmospheric gases completely from liquid samples [40]. The results for [C4mim]Cl are also compared with the experimental data of He et al. [41], Machida et al. [42] and Govinda et. al. [43], as shown in Fig. 3B. Liquid structure: Radial distribution function To understand and verify the effect of electrostatic forces on the intermolecular liquid structure of target ILs, the radial distribution functions (g(r) plots) between different atomic sites are analyzed at the simulation temperature 450 K. Figs. 4A and 4B show g(r) plots between chloride anion and acidic hydrogen (HR) on imidazolium ring cation for [C1mim]Cl and [C4mim]Cl ILs, respectively. The g(r) plots of Cl −  H R for all six force fields appear to show two peaks: a sharp first peak at around 2.7 Å and a well-defined second peak at around 6.4 Å as Fig. 4 illustrates. All force fields show radial distribution functions that are qualitatively very similar except some differences in peak intensities imposed by the difference in the partial atomic charges. Sharper peaks are obtained using NBO based force fields. For [C1mim]Cl, the first peak intensity decreases in the

9

order of HF-NBO

MP2-NBO

HF-CHELPG

MP2-CHELPG

B3LYP-NBO

B3LYP-

CHELPG, consistent with the trend of chloride anion charge. Higher density of electrons on the chloride anions results in them interacting strongly with their neighbors. Another interesting observation is that the force fields based on CHELPG charges resulted in poorly structured N1 N1 correlations for [C1mim]Cl. Fig. 5 sheds light on this cation-cation structural aspect. The

force fields based on NBO charges show N1 N1 correlating with rather distinct peak at around 5.0 Å. The correlation between the terminal carbon atoms of the butyl tail of [C4mim]Cl is also studied by the radial distribution functions as shown in Fig. 6. It can be seen that the raising edge of the first peak is rather sharp from which the distance of closest contact can be determined to be 3.0 Å for all force fields. Quite interestingly, the position of the first peak is sharper with HF-CHELPG charges although the charge value of butyl group is lower compared to the NBO results. We also calculated the coordination number n by integrating the Cl −  H R radial distribution function up to the first solvation shell radius. The coordination number of cation-anion for each force field can be verified in Table 3. The coordination numbers are reduced by about 1.4 in the [C4mim]Cl system compared to the [C1mim]Cl system. This reduction could be due to the steric effect of butyl chain in [C4mim]Cl. The anion-cation structural relation of both ILs are further unraveled using spatial distribution functions (SDFs), which allow concluding the effect of various partial charges on the morphology of the anion-cation in this family of liquids. Each trajectory for SDF calculation is evaluated using TRAVIS program package [44]. Spatial distribution of chloride anion around cation is shown in Figs. 7 and 8 for [C1mim]Cl and [C4mim]Cl, respectively. SDFs are demonstrated by applying the same isosurface values for both systems, enabling to study the effect of partial charges on spatial distribution of anion around cation reliably.

Consequently, different quantum mechanically 10

derived partial charges produce relatively similar results that are in accordance with the results of ab initio molecular dynamics simulations obtained by Kirchner and co-workers [45]. In all cases, three preferred binding sites around the three ring hydrogens are observed as can be seen in Fig. 7 for [C1mim]Cl. Anions are also present next to the methyl groups. The results with CHELPG charges indicate that contact between the chloride and methyl hydrogens is of lower probability than with NBO charges. This could be due to the greater values of positive charges on the hydrogen atoms of methyl group. The three-dimensional organization of the chloride ions around the [C1mim]+ cation is shown in Fig. 7. The favor sites of interactions for chloride ions are centered on the three hydrogen atoms of the imidazolium ring as expected. It is also seen that anions favor the methyl side rather than the butyl side for [C4mim]Cl IL (see Fig. 8). Transport Properties: Diffusion and Viscosity While the different quantum mechanically derived partial atomic charges do not affect strongly the structural properties (density and radial distribution function) of chloride based ionic liquids, instead, major effects are found on the energetics and dynamics of these systems. Non-polarizable models have been remarkably successful in modeling many complex molecular systems. However, the simulation of polarization effects which was found to affect mainly the dynamics of ionic systems meet serious consideration. Despite the drastic simplifications, one issue of concern, however, is that the electrostatic interactions of ions are described in standard non-polarizable force fields, by their original integer charges (

), i.e. as if these ions were in vacuum, completely disregarding the effect of electronic dielectric

screening inherent to the condensed phase medium. Obviously, the electric potential of charged residues in such calculations should reflect the electronic screening of the medium. Therefore, we should account electronic polarization in non-polarizable force fields by choosing a suitable reduced atomic charge which could reproduce experimental results.

11

To evaluate the influence of different force fields on dynamical properties, we examined the corresponding simulated mean square displacements (MSDs). The ensemble-averaged MSD is calculated from the knowledge of trajectories of particles: N

MSD = ∑ ric (t ) - ric (0) i =1

2

2

= ∆ r(t ) ,

(2)

where ric (t ) is the location of the center of mass of ion i at time t. The MSDs of the cations in [C1mim]Cl and [C4mim]Cl ILs obtained from 8 ns simulations at 450 K are shown in Fig. 9. (for anions see Fig. S5 and S6). For both ILs, the dynamics of the ions significantly increased when CHELPG charges are used. Interestingly, MSDs tend to change in order of B3LYP-CHELPG > MP2-CHELPG > HF-CHELPG > B3LYP-NBO > MP2-NBO > HFNBO which is contrary to the spatial distribution of anion around cation as shown in Figs. 7 and 8. The main reason behind these observations lies in the nature of the method used for the calculation of atomic charges. For example, in the case of [C1mim]Cl, the charge values obtained by B3LYP-CHELPG are

0.719 e which show the largest deviation from unity of all schemes assessed. This charge reduction

leads to the higher dynamics of the systems modeled by B3LYP-CHELPG method. In the same manner, HF-NBO gives the largest absolute charge values and the lower dynamics. In the other word, the higher the spatial probability of anion around the cation the lower the MSD of both cations and anions. So there are strongest anion-cation associations in the case of HF-NBO charges. The diffusion coefficient, Di, was obtained from the linear regime of MSD curves at long simulation time using the Einstein relation,

Di =

2 1 d lim ric (t ) − ric (0) 6 t → ∞ dt

(3)

12

The calculated values of diffusion coefficient for [C1mim]Cl and [C4mim]Cl are in the order of 10-10 m2.s-1 as are shown in Tables 4 and 5, respectively. It should be noted that using CHELPG charges reduces the long-range interactions and thus increases the diffusion coefficient and ion mobility. Then, the viscosities were estimated using Stokes–Einstein equation, η=

k BT cπ Di rs,i

(4)

where, k B is the Boltzmann constant, c is a constant, and rs,i is the effective hydrodynamic (Stokes) radius of ion i. The value of c depends on the nature of the system and the interactions experienced by the solute. For large molecular size solutes in small molecular size solvent environments, c can be as high as 6. As the solute and solvent become more similar in size, especially in higher viscous liquids, the value of c reduces to ~ 4. [29] Tokuda et al. [46] have shown that ion diffusivity in the ILs obey Eq 4. While we applied c=4 in all calculations, values of rs

for

, [C1mim]+ and [C4mim]+ ions were adopted as 0.18, 0.287, and 0.335 nm,

respectively (ref. [47] for cations and ref. [48] for the anion). While ILs are highly hygroscopic by nature and moisture have substantial influence on viscosity, it is difficult to ascertain that simulated and experimental viscosities could be the same. The results of calculated viscosities for [C4mim]Cl IL with different potentials are shown in Fig. 10 and compared with the experimentally measured viscosities reported by Prausnitz and co-workers [49]. Accordingly, viscosities obtained by HF-NBO and MP2-NBO are higher substantially as compared with experimental values. The viscosities simulated by B3LYP-NBO charges are closer to experimental results (see Fig. 10). However, by using HF-CHELPG, MP2- CHELPG, and B3LYP-CHELPG charges, the viscosities are lower than the experimental data.

13

The transference numbers measure the relative ability of specific ions (cations or anions) to carry charge. Cationic and anionic transference numbers are defined based on diffusion coefficient of anion and cation:

t+ =

D+ D− and t− = ( D+ + D− ) ( D+ + D− )

(5)

where t+ + t− = 1 and 0 ≤ ti ≤ 1. The calculated values of cationic and anionic transference numbers are also shown in Tables 4 and 5 for [C1mim]Cl and [C4mim]Cl, respectively. In general, for ILs and especially for imidazolium-based ILs, the cationic transference number is greater than the anionic transference number. Therefore, the results reveal that transference number were properly represented by the B3LYP-NBO scheme. The molar electrical conductivity, Λ , was also calculated from diffusion coefficients for each ionic liquid (see Tables 4 and 5) using Nernst-Einstein equation, Λ = ( z cation Dcation + z anion Danion )F 2 / RT

(6)

where F is the Faraday constant, R is the gas constant, and z is the charges. In agreement with experiment, [50] the electrical conductivity of the ILs decrease with increasing chain length due to the higher van der Waals (vdW) interactions between the longer alkyl chains. The trends for electrical conductivities obtained from simulations with different force field charges are B3LYPCHELPG > MP2-CHELPG > HF-CHELPG > B3LYP-NBO > MP2-NBO > HF-NBO. From another viewpoint, the smaller electrical charge density obtained by B3LYP-CHELPG method leads to a greater mobility of ions. Therefore, the simulations performed by B3LYP-CHELPG charges possess the highest electrical conductivity. Thermodynamic properties

14

The enthalpy of vaporization is calculated by the following relationship to estimate the strength of non-bonding interactions for simulated systems [51]: (7)

∆HVap = H gas − H liquid = U gas + RT − U liquid

where all terms are self-defined. The PV term for the liquid phase was neglected in the relation of enthalpy and internal energy (H = U + PV). It has been assumed that PV = RT (ideal gas state) and the gas consists of non-interacting pairs of ions (cation + anion). Thus, the gas phase in a computer experiment was simulated in terms of a single ion pair in vacuum. The first molecular dynamics simulation of the enthalpy of vaporization of ILs have been performed by Maginn and coworkers [10,52]. They predicted the ∆H Vap in good agreement with experimental data [53,54]. ∆H Vap is available for [C4mim]Cl, [55] but not for [C1mim]Cl. The reported experimental values

of ∆HVap at 450 K are 138.5 and 137.2 kJ/mol measured by quartz crystal microbalance (QCM) and thermogravimetric analysis (TGA) methods, respectively [56]. The results of vaporization enthalpies obtained by different force field charges are shown in Tables 4 and 5 for [C1mim]Cl and [C4mim]Cl, respectively. The results obtained with different charges decrease in the order of HF-NBO > MP2-NBO > B3LYP-NBO > HF-CHELPG > MP2CHELPG > B3LYP-CHELPG which confirm with a reduction in cation-anion interaction. Energies of vaporization with HF charges are higher and is the highest with HF-NBO, which matches with the highest total charge obtained. The obtained value of the vaporization enthalpy for [C4mim]Cl IL with HF/CHELPG charges is 135 kJ/mol which agrees well with the experimental result. The vaporization enthalpies obtained from simulations with NBO charges slightly overestimate the experimental data. Dipole Distribution

15

Dipole moments were calculated for [C1mim] + and [C4mim]+ cations with the atomic charges obtained by different methods. The magnitude of the dipole moment on cation, µ , given by equation (8): 



(8)

µ = ∑ qi × ri − r0 i

  where qi and ri are the charge and the position of atom i, respectively, and r0 is some point of reference, typically the center of the imidazolium ring for [Cnmim]+. Fig. 11 shows the probability distribution of [C1mim]+ and [C4mim]+ dipole moments, respectively. In the case of [C1mim]+, the trend and the band width of broadening of the dipole distributions with NBO charges are completely different with that of CHELPG charges. These broadenings are due to greater fluctuations in the geometry of [C1mim] + as a result of method used for the partial atomic charges. Results for the dipole moments of [C1mim]+ with different charges are given in Fig. 11A for comparison.

The distribution width and shape from HF-

CHELPG charges are very similar to those of MP2-CHELPG charges. The predicted average dipole moments are in the range of 2.49 D to 2.96 D for [C1mim]+ cation. Therefore, the average field induced in this system due to different atomic charges changes the molecular dipole moment by about 0.47 D. To evaluate the results of these classical force fields, the comparison of the dipole moments by ab initio computation might be helpful. Lynden-Bell and co-workers obtained a histogram centered on 2.8D for dipole moment distribution in the [C1mim]+ cation at 450K and the liquid phase condition [57]. Krekeler et al. indicated that the average value for dipole moment of this cation is about 2.8D based on their CPMD calculations. [58] Wendler et al. observed similar distributions for the [C1mim]Cl cluster consisting of 30 ion pairs. In their case the dipole moment fell in the range between 0.0 and 5.0 D, with the maximum being around 3.0 D, which is in a good

16

agreement with our results from the NBO charges. [59] Fig. 11B reports the dipole distributions for [C4mim]+. Remarkably, the ion dipole moment distributions with NBO potential is shifted toward higher values (polarization) compared to CHELPG charges. CONCLUSIONS This paper examined the effect of various quantum mechanically derived partial atomic charges on the bulk properties of [C1mim]Cl and [C4mim]Cl ionic liquids that provided deeper insight into force field generation for these class of ILs. The partial charges for these ILs obtained by NBO and CHELPG analysis of ab initio calculations were performed by HF, B3LYP, and MP2 methods all with 6-311++G(d,p) basis set. The chloride ILs have presented an unusual trend in terms of their performance with respect to the charge schemes studied. It was observed that the RESP schemes are preferred in the case of

,

and iodide ionic liquids when producing atomic charges which can effectively capture electronic polarization effects and show an appreciable degree of charge transfer. In contrast, the chloride systems present an unusual strong inter-ionic bonds between the imidazolium ring and the chloride anion. This type of interaction elicits a most notable response from the NBO charge scheme. This findings are in line with the results of Rigby and Izgorodina which indicated that charge transfer in the [C1mim]Cl clusters converges to a lower value, showing an unusually strong inter-ionic bond with the chloride anion. [60]The general pattern claims that the charge on the chloride ion calculated from NBO method is larger in magnitude than CHELPG method. The highest absolute total charge values ( 0.939 e and

0.924 e for [C1mim]Cl and [C4mim]Cl,

respectively ) were obtained by HF-NBO calculations and the smallest values ( 0.719 e and 0.719 e for

[C1mim]Cl and [C4mim]Cl, respectively) were obtained by B3LYP-CHELPG

calculations.

17

The systematic change of the correlation function with partial charges in the order of HF-NBO > MP2NBO > HF-CHELPG> MP2-CHELPG > B3LYP-NBO > B3LYP-CHELPG is appealing as much as the electronic charges represent the electronic structure and hence the macroscopic density. The intensity and location of the first peak of g(r) is a function of charges generated. These structural changes play an important role on the dynamics of the system. The B3LYP-NBO force field charge is able to reproduce the experimental values of viscosity for [C4mim]Cl rather accurately. Simulated diffusion coefficients, MSDs and electrical conductivities are function of potentials in the order of B3LYP-CHELPG > MP2-CHELPG > HF-CHELPG > B3LYP-NBO > MP2-NBO > HF-NBO. A remarkable increase in the probability of presence of anion near the above of HR atom and below of HW and HZ atoms is concurrent with a significant decrease in the MSD values. This phenomenon leads to a reduction in the mobility of ions and therefore in their contribution to conductivity [61]. The predictions of electrical dipole moment of MP2-NBO results are in well agreement with ab initio MD dipole moment. The larger dipole moment distributions observed for simulations with NBO charges which indicate a stronger polarization induced in the system compared to CHELPG charges. Our results show that the static structures represented as the density and RDF are not affected strongly by the reorganization of atomic partial charges in IL, however, dynamical properties are affected due to the polarization effects. In addition, our benchmark calculations show that smaller differences between the results from B3LYP-NBO charges and available experimental data are observed when dynamical properties are considered. Overall, B3LYP-NBO calculations yield partial atomic charges being retained on the anions and cations, which could reproduce structural, dynamical and energetic properties of these ILs properly.

18

ACKNOWLEDGMENTS The authors are indebted to the research council of the Shiraz University for the financial supports.

REFERENCES [1] H. Weingärtner, Angew. Chem. Int. Ed. 47 (2008) 654. [2] C. J. Cramer, Essentials of computational chemistry: theories and models, John Wiley & Sons, 2013. [3] C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem. 97 (1993) 10269. [4] C. Chipot, B. Maigret, J. L. Rivail and H. A. Scheraga, J. Phys. Chem. 96 (1992) 10276. [5] L. E. Chirlian and M. M. Francl, J. Comput. Chem. 8 (1987) 894. [6] C. M. Breneman and K. B. Wiberg, J. Comput. Chem. 11 (1990) 361. [7] E. Sigfridsson and U. Ryde, J. Comput. Chem. 19 (1998) 377.

19

[8] A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys. 83 (1985) 735. [9] C. M. Tenney, M. Massel, J. M. Mayes, M. Sen, J. F. Brennecke and E. J. Maginn, J. Chem. Eng. Data 59 (2014) 391. [10] T. I. Morrow and E. J. Maginn, J. Phys. Chem. B 106 (2002) 12807. [11] B. Bhargava and S. Balasubramanian, J. Chem. Phys. 123 (2005) 144505. [12] F. Dommert and C. Holm, Phys. Chem. Chem. Phys. 15 (2013) 2037. [13] F. Dommert, K. Wendler, R. Berger, L. Delle Site and C. Holm, ChemPhysChem 13 (2012) 1625. [14] J. Schmidt, C. Krekeler, F. Dommert, Y. Zhao, R. Berger, L. D. Site and C. Holm, J. Phys. Chem. B 114 (2010) 6150. [15] A. Mondal and S. Balasubramanian, J. Phys. Chem. B 118 (2014) 3409. [16] A. Mondal and S. Balasubramanian, J. Phys. Chem. B 119 (2015) 11041. [17] A. Mondal and S. Balasubramanian, J. Chem. Eng. Data 59 (2014) 3061. [18] M. H. Ghatee, A. R. Zolghadr, F. Moosavi and Y. Ansari, J. Chem. Phys. 136 (2012) 124706. [19] M. H. Ghatee and A. R. Zolghadr, J. Phys. Chem. C 117 (2013) 2066. [20] A. R. Zolghadr, M. H. Ghatee and A. Zolghadr, J. Phys. Chem. C 118 (2014) 19889. [21] A. Szabo, N.S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Macmillan Publishing, New York, 1982. [22] C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37 (1988) 785. [23] C. Møller and M. S. Plesset, Phys. Rev. 46 (1934) 618. 20

[24] M. J. Frisch, G. W. Trucks and H. B. Schlegel, et al., Gaussian 03, Revision C.02, Gaussian Inc., Wallingford, CT, 2004. [25] M. A. Ratner, G. C. Schatz, Introduction to Quantum Mechanics in Chemistry, Prentice-Hall, Inc., 2001. [26] J. N. Canongia Lopes, J. Deschamps and A. A. Pádua, J. Phys. Chem. B 108 (2004) 2038. [27] J. N. Canongia Lopes and A. A. Pádua, J. Phys. Chem. B 110 (2006) 19586. [28] T. Yan, C. J. Burnham, M. G. Del Pópolo and G. A. Voth, J. Phys. Chem. B 108 (2004) 11877. [29] M. Kowsari, S. Alavi, M. Ashrafizaadeh and B. Najafi, J. Chem. Phys. 130 (2009) 014703. [30] J. N. Canongia Lopes, J. Deschamps and A. A. Pádua, J. Phys. Chem. B 108 (2004) 11250. [31] W. Smith, T. R. Forester, I. T. Todorov, The DL-POLY molecular simulation package, V. 2.17, Daresbury Laboratory, UK, 2006, www.cse.scitech.ac.uk/ccg/ software/DL_POLY/. [32] S. Nosé, J. Chem. Phys. 81 (1984) 511. [33] W. G. Hoover, Phys. Rev. A 1695 (1985) 31. [34] S. Melchionna, G. Ciccotti and B. Lee Holian, Mol. Phys. 78 (1993) 533. [35] F. Maier, T. Cremer, C. Kolbeck, K. Lovelock, N. Paape, P. Schulz, P. Wasserscheid and H.-P. Steinrück, Phys. Chem. Chem. Phys. 12 (2010) 1905. [36] P. A. Hunt, B. Kirchner and T. Welton, Chem. Eur. J.12 (2006) 6762. [37] A. A. Fannin Jr, D. A. Floreani, L. A. King, J. S. Landers, B. J. Piersma, D. J. Stech, R. L. Vaughn, J. S. Wilkes and L. Williams John, J. Phys. Chem. 88 (1984) 2614. [38] B. Bhargava and S. Balasubramanian, J. Chem. Phys. 123 (2005) 144505.

21

[39] O. N. Starovoytov, H. Torabifard and G. A. Cisneros, J. Phys. Chem. B 118 (2014) 7156. [40] M. H. Ghatee, F. Moosavi, A. R. Zolghadr and R. Jahromi, Ind. Eng. Chem. Res. 49 (2010) 12696. [41] R.-H. He, B.-W. Long, Y.-Z. Lu, H. Meng and C.-X. Li, J. Chem. Eng. Data 57 (2012) 2936. [42] H. Machida, R. Taguchi, Y. Sato and J. Smith, Richard L, J. Chem. Eng. Data 56 (2011) 923. [43] V. Govinda, P. Attri, P. Venkatesu and P. Venkateswarlu, Fluid Phase Equilibria 304 (2011) 35. [44] M. Brehm and B. Kirchner, J. Chem. Inf. Model. 51 (2011) 2007. [45] S. Zahn, M. Brehm, M. Bruessel, O. Holloczki, M. Kohagen, S. Lehmann, F. Malberg, A. S. Pensado, M. Schoeppke, H. Weber, and B. Kirchner, Journal of Molecular Liquids 192 (2014) 71. [46] a) H. Tokuda, K. Hayamizu, K. Ishii, M. A. B. H. Susan and M. Watanabe, J. Phys. Chem. B 108 (2004) 16593.; b) H. Tokuda, S. Tsuzuki, M. A. B. H. Susan, K. Hayamizu and M. Watanabe, J. Phys. Chem. B 110 (2006) 19593. [47] H. Tokuda, K. Hayamizu, K. Ishii, M. A. B. H. Susan and M. Watanabe, J. Phys. Chem. B 109, (2005) 6103. [48] J. Baxendale and P. Bevan, J. Chem. Soc. A (1969) 2240. [49] S. Fendt, S. Padmanabhan, H. W. Blanch and J. M. Prausnitz, J. Chem. Eng. Data 56 (2010) 31. [50] J. Leys, M. Wübbenhorst, C. P. Menon, R. Rajesh, J. Thoen, C. Glorieux, P. Nockemann, B. Thijs, K. Binnemans and S. Longuemart, J. Chem. Phys. 128 (2008) 064509. [51] W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225. [52] E. J. Maginn, J. Phys.: Condens. Matter 21 (2009) 373101.

22

[53] C. Cadena and E. J. Maginn, J. Phys. Chem. B 110 (2006) 18026. [54] K. Fumino, A. Wulf, S. P. Verevkin, A. Heintz and R. Ludwig, ChemPhysChem 11 (2010) 1623. [55] S. P. Verevkin, D. H. Zaitsau, V. N. Emel'yanenko, C. Schick, S. Jayaraman and E. J. Maginn, Chem. Commun. 48 (2012) 6915. [56] G. J. Kabo, Y. U. Paulechka, D. H. Zaitsau and A. S. Firaha, Thermochimica Acta 609 (2015) 7. [57] C. E. R. Prado, M. G. D. Pópolo, T. Youngs, J. Kohanoff and R. Lynden-Bell, Mol. Phys. 104 (2006) 2477. [58] C. Krekeler, F. Dommert, J. Schmidt, Y. Zhao, C. Holm, R. Berger and L. Delle Site, Phys. Chem. Chem. Phys. 12 (2010) 1817. [59] K. Wendler, S. Zahn, F. Dommert, R. Berger, C. Holm, B. Kirchner and L. Delle Site, J Chem Theory Comput. 7 (2011) 3040. [60] J. Rigby, E. I. Izgorodina, Phys.Chem. Chem. Phys. 15 (2013) 1632. [61] W. Zhao, F. Leroy, B. Heggen, S. Zahn, B. Kirchner, S. Balasubramanian and F. MüllerPlathe, J. Am. Chem. Soc. 131 (2009) 15825.

23

FIGURES (A)

(B)

Fig. 1. Structure and labeling of (A) [C1mim]Cl and (B) [C 4mim]Cl IL.

(A)

(B)

24

Fig. 2. Partial charges of different groups for single ion pair of (A) [C1mim]Cl and (B) [C4mim]Cl IL obtained by different methods.

(B)

(A)

Fig. 3. Temperature dependence of density for (A) [C1mim]Cl and (B) [C4mim]Cl.

(A)

Fig. 4. Comparison of H R  Cl radial distribution function in (A) [C1mim]Cl and (B) [C4mim]Cl systems using charges obtained by different methods.

25

(A)

(B)

Fig. 5. Comparison of N 1  N 1 radial distribution function in (A) [C1mim]Cl and (B) [C4mim]Cl systems using charges obtained by different methods.

Fig. 6. Comparison of methods.

C 4 C4 radial distribution function in [C4mim]Cl system using charges obtained by different

26

(a)

(d)

(c)

(b)

(e)

(f)

Fig. 7. Spatial distribution of anion around the cation in [C 1mim]Cl. (a) B3LYP-CHELPG, (b) HF-CHELPG, (c)

MP2-CHELPG, (d) B3LYP-NBO, (e) HF-NBO and (f) MP2-NBO. Cl− green; C, orange; H, white; N, blue.

(d)

(c)

(b)

(a)

(f)

(e)

Fig. 8. Spatial distribution of anion around the cation in [C 4mim]Cl. (a) B3LYP-CHELPG, (b) HF-CHELPG, (c)

MP2-CHELPG, (d) B3LYP-NBO, (e) HF-NBO and (f) MP2-NBO. Cl− green; C, orange; H, white; N, blue.

(A)

(B)

27

Fig. 9. Mean-square displacement of geometrical center of imidazolium cation for (A) [C1mim]Cl and (B) [C4mim]Cl.

Fig. 10. Viscosities of [C4mim]Cl IL as a function of temperature. Lines are only a guide to the eye.

(B)

(A)

+

+

Fig. 11. Simulated dipole moment distribution of (A) [C1mim] and (B) [C4mim] ensemble using charges obtained by

various methods.

28

TABLES

Table 1: Partial Atomic Charges for [C1mim]Cl Molecule Derived from HF, B3LYP and MP2 with 6311++G(d,p) Basis Set Using CHELPG and NBO Methods. For Atom Labels, See Fig. 1A. atom N1 CR N2 CW CZ CA C1 HR HW HZ HA H1 Cl

HF NBO CHELPG q/e q/e -0.451 -0.009 0.496 0.045 -0.451 -0.009 0.012 -0.123 0.012 -0.125 -0.253 0.054 -0.253 0.050 0.240 0.241 0.213 0.179 0.213 0.180 0.194 0.060 0.194 0.061 -0.939 -0.846

B3LYP NBO CHELPG q/e q/e -0.390 -0.071 0.337 0.100 -0.390 -0.056 -0.025 -0.129 -0.024 -0.08 -0.353 0.111 -0.356 0.071 0.250 0.18 0.224 0.165 0.224 0.154 0.218 0.040 0.219 0.052 -0.806 -0.719

NBO q/e -0.472 0.504 -0.472 0.015 0.015 -0.263 -0.263 0.231 0.215 0.215 0.196 0.196 -0.898

MP2 CHELPG q/e -0.089 0.071 -0.089 -0.100 -0.099 0.110 0.109 0.247 0.176 0.175 0.047 0.048 -0.797

Table 2: Partial Atomic Charges for [C4mim]Cl Molecule Derived from HF, B3LYP and MP2 with 6311++G(d,p) Basis Set Using CHELPG and NBO Methods. For Atom Labels, See Fig. 1B.

atom N1 CR N2 CW CZ CA C1 C2

NBO q/e -0.440 0.409 -0.447 -0.008 0.011 -0.255 -0.095 -0.337

HF CHELPG q/e 0.087 0.120 0.066 -0.174 -0.233 -0.05 -0.037 0.193

B3LYP NBO CHELPG q/e q/e -0.365 0.108 0.285 0.036 -0.375 0.11 -0.039 -0.17 -0.023 -0.189 -0.351 -0.082 -0.18 -0.016 -0.388 0.184

29

NBO q/e -0.454 0.394 -0.458 -0.015 -0.015 -0.255 -0.101 -0.343

MP2 CHELPG q/e 0.085 0.113 0.053 -0.204 -0.215 -0.049 -0.08 0.212

C3 C4 HR HW HZ HA H1 H2 H3 H4 Cl

-0.33 -0.492 0.299 0.214 0.215 0.182 0.235 0.155 0.184 0.162 -0.924

0.175 -0.202 0.174 0.192 0.223 0.079 0.053 -0.047 -0.032 0.042 -0.840

-0.384 -0.569 0.286 0.226 0.228 0.218 0.231 0.200 0.196 0.197 -0.850

0.182 -0.207 0.144 0.175 0.193 0.079 0.054 -0.072 -0.01 0.033 -0.765

-0.335 -0.494 0.304 0.216 0.218 0.192 0.204 0.178 0.172 0.172 -0.892

0.183 -0.207 0.191 0.202 0.218 0.077 0.063 -0.050 -0.032 0.040 -0.815

Table 3. The Coordination Number, n, for the chloride anions in the first coordination shell of each HR. The position of first minimum of the corresponding g(r), r, are also listed. Ionic Liquid B3LYP-CHELPG MP2-CHELPG HF-CHELPG B3LYP-NBO MP2-NBO HF-NBO

[C1mim]Cl r n 4.52 2.06 4.42 2.08 4.42 2.03 4.27 1.71 4.22 1.83 4.22 1.86

[C4mim]Cl r n 4.77 1.41 4.67 1.43 4.67 1.38 4.27 1.19 4.32 1.26 4.27 1.28

Table 4. Simulated Cation and Anion Diffusion Coefficients Di (in 10−10 m2.s-1), Molar Electrical Conductivity, (in 10-4 S.m2.mol-1), transference numbers, η (in mPa s), and Heat of Vaporization (in kJ/mol) of [C1mim]Cl IL at 450 K for Simulations with Various Charge Sets. Ionic Liquid

D+

D-

B3LYP-CHELPG MP2-CHELPG HF-CHELPG B3LYP-NBO MP2-NBO HF-NBO

14.6 6.79 4.44 4.67 2.68 1.16

13.5 7.32 4.44 4.55 2.00 1.00

Λ 50.0 31.5 18.7 18.5 10.4 4.87

t+

t-

η

∆HVap

0.519 0.481 0.500 0.507 0.572 0.547

0.481 0.519 0.500 0.493 0.428 0.453

0.84 1.69 2.66 2.56 5.83 12.93

108 122 129 138 147 162

Table 5. Simulated Cation and Anion Diffusion Coefficients Di (in 10−10 m2.s-1), Molar Electrical Conductivity, Λ (in 10-4 S.m2.mol-1), transference numbers, η (in mPa s), and Heat of Vaporization (in kJ/mol) of [C4mim]Cl IL at 450 K for Simulations with Various Charge Sets. B3LYP-CHELPG MP2-CHELPG HF-CHELPG B3LYP-NBO MP2-NBO HF-NBO

D+

D-

Λ

t+

t-

η

∆ HVap

11.6 7.31 4.92 4.26 1.32 0.57

11.2 5.84 5.05 3.29 1.35 0.54

43.3 26.6 20.8 16.0 5.81 2.59

0.508 0.556 0.494 0.565 0.496 0.512

0.492 0.444 0.506 0.435 0.504 0.488

1.22 2.08 2.81 3.62 8.89 21.26

119 129 135 140 151 167

30

GRAPHICAL ABSTRACT

31