Improved Force Field Model for the Deep Eutectic Solvent Ethaline

Article pubs.acs.org/JPCB

Improved Force Field Model for the Deep Eutectic Solvent Ethaline: Reliable Physicochemical Properties Elisabete S. C. Ferreira,*,†,‡ Iuliia V. Voroshylova,†,‡ Carlos M. Pereira,*,‡ and M. Natália D. S. Cordeiro*,† †

Departamento de Química e Bioquímica, LAQV@REQUIMTE, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal ‡ Departamento de Química e Bioquímica, CIQ(UP), Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal S Supporting Information *

ABSTRACT: In this work, we combined various parameters found in the literature for the choline cation, chloride anion, and ethylene glycol to set up force field models (FFMs) for a eutectic mixture, namely, ethaline (1:2 choline chloride/ ethylene glycol (ChCl:2EG)). The validation of these models was carried out on the basis of physical and chemical properties, such as the density, expansion coefficient, enthalpy of vaporization, self-diffusion coefficients, isothermal compressibility, surface tension, and shear viscosity. After the initial evaluation of the FFMs, a refinement was found necessary and accomplished by taking into account polarization effects in a mean-field manner. This was achieved by rescaling the electrostatic charges of the ions based on partial charges derived from ab initio molecular dynamics (MD) simulations of the bulk system. Classical all-atom MD simulations performed over a large range of temperatures (298.15−373.15 K) using the refined FFMs clearly showed improved results, allowing a better prediction of experimental properties. Specific structural properties (radial distribution functions and hydrogen bonding) were then analyzed in order to support the adequacy of the proposed refinement. The final selected FFM leads to excellent agreement between simulated and experimental data on dynamic and structural properties. Moreover, compared to the previously reported force field model (Perkins, S. L.; Painter, P.; Colina, C. M. Experimental and Computational Studies of Choline Chloride-Based Deep Eutectic Solvents. J. Chem. Eng. Data 2014, 59, 3652−3662), a 10% improvement in simulated transport properties, i.e., self-diffusion coefficients, was achieved. The isothermal compressibility, surface tension, and shear viscosity for ethaline are accessed in MD simulations for the first time.

1. INTRODUCTION Since the 1990s, room-temperature ionic liquids (RTILs) have attracted a great deal of interest because of their potential to provide new ways to solve many problems in chemistry and industry. Being known as green solvents, their potential applications vary from electrochemical, analytical, synthesis, and engineering processes.1−6 RTILs’ properties include low volatility, high conductivity, large potential windows, excellent thermal and chemical stability, and a high solubility of organic and inorganic materials.4,7 However, the drawback for largescale applications is their high cost. Deep eutectic solvents (DESs) are a special group of ionic liquids (ILs), which differ from RTILs mainly because they are easy to synthesize and cheaper.6,8−10 DES is a liquid typically composed of at least one salt and a neutral specie that form a eutectic mixture, often through hydrogen bonding interactions. The melting point of DES is lower than either of its individual components.11−13 One of the most interesting features of using DESs is their tunability because the components can be picked up according to the final application.12,14 A large number of combinations are © XXXX American Chemical Society

possible for this type of mixture. Compared to conventional RTILs, there is a lack of experimental and theoretical information on fundamental thermodynamic and transport properties of DESs,15 although some choline chloride-based eutectic mixtures have been explored experimentally.2−4,16 Among all types of DESs, the most promising ones are those obtained by mixing a quaternary ammonium salt, typically choline chloride (ChCl) and a hydrogen bond donor (HBD) molecule.4 For example, ethaline is a DES obtained from mixing choline chloride (ChCl) with the HBD ethylene glycol in a molar ratio of 1:2 (ChCl:2EG). Ethaline melts at 213.15 K, a temperature much lower than that of most RTILs.17 It is nonflammable, nontoxic, completely miscible with water,17 and readily biodegradable.10,18 All of these properties make ethaline a good substitute for volatile organic solvents and RTILs. Received: July 19, 2016 Revised: August 29, 2016

A

DOI: 10.1021/acs.jpcb.6b07233 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 1. Description of the Force Field Model Combinations Studieda chloride (Cl−)

choline (Ch+)

a

ethylene glycol (EG)

FFM

charges

parameters

charges

parameters

charges

parameters

FFM1 FFM2 FFM3 FFM4 FFM5 FFM6 FFM7 FFM8

30 31

OPLS-AA OPLS-AA

OPLS-AA OPLS-AA

30

OPLS-AA OPLS-AA OPLS-AA

OPLS-AA OPLS-AA

30

OPLS-AA OPLS-AA OPLS-AA

32

32

31

31

35

35

36

36

31

31

35

35

32

32

30

30

35

35

OPLS-AA

OPLS-AA

30

30

35

35

36

36

30

30

35

35

32

32

Charges and remaining parameters were taken from the included references.

procedure using the restrained electrostatic charge potential (RESP) fit procedure. The authors opt, however, to use reduced ±0.9e charges for the ionic species. In so doing, the final FFM proposed by Perkins et al.29 showed very good agreement for the ethaline’s density and other thermodynamic properties (volume expansivity and heat capacity), but transport properties were ∼20−30% underestimated, compared to experimental values. The main goal of the present work is to set up an all-atom force field model that best reproduces the structural, dynamic, and thermodynamic properties of ethaline. To develop such a FFM, we combined various models available in the literature30−36 for individual components of ethaline and analyzed the performance of the resulting force fields. Thus, MD simulations were carried out to obtain a set of physicochemical properties of ethaline, such as the density, molar volume, isothermal compressibility, shear viscosity, surface tension, self-diffusion coefficients, and enthalpy of vaporization. As a result, a parametrization of the force fields was deemed necessary to cope properly with the transport properties of ethaline. Therefore, ab initio MD calculations coupled with an RESP fitting procedure were performed to guide the rescaling of charges of the studied force fields. Additionally, radial distribution functions and hydrogen bonding were analyzed to probe the structure and interactions of ethaline and finally settle the best refined force field for this deep eutectic solvent. The article is organized as follows. In Section 2, the details of the theoretical methods employed, including the force field model parameters, are provided, as is the setup for classical and ab initio MD simulations, and the details for the calculation of the different physicochemical properties are described. The results are presented and discussed in Section 3, and the conclusions are given in Section 4.

Reports on physical and chemical properties of ethaline start to appear in the literature. Density is the most studied property, being obtained both experimentally2,15,17,19−21 and through predictive models.15,22 Regarding other properties, the experimentally measured isothermal compressibility,19 surface tension,11,21 viscosity20,23−25 and self-diffusion23 can be mentioned. D’Agostino et al.23 showed that the interactions between the choline cation and the hydrogen bond donor become weaker with increasing temperature. According to the authors, the choline cation diffuses more slowely than ethylene glycol because of its larger size and higher molecular weight compared to those of the HBD molecule.26 Three main advantages of substituting commercially available solvents with ethaline in processes of electropolishing stainless steel have been shown by Abbott et al.27 These are (i) high current efficiencies; (ii) negligible gas evolution at the anode/ solution interface during polishing; and (iii) benign properties and corrosiveness of ethaline compared to those of the currently used aqueous acid solutions. Popescu et al.28 showed, in a series of electrodepositions of Cu from CuCl2 dissolved in choline chloride-based DES with different HBDs (e.g., urea, malonic acid, oxalic acid, and ethylene glycol), that better deposits are attained when using the ChCl:2EG solvent. In previous work,4 our group described the study of the electrodeposition mechanism of tin onto glassy carbon electrodes, and excellent results were also obtained for ChCl:2EG. However, a more extended study is required in order to understand the processes governing the deposit growth in different DESs. Computer modeling of such electrochemical processes may contribute to the elucidation of the phenomena controlling the nucleation rate constant and thereby reveal possible methods of improvement. Indeed, despite extensive efforts, it is extremely challenging to experimentally probe the interactions involved between the different components of these complex mixtures.26 In particular, MD simulations can both predict the physicochemical properties and help in understanding the interactions and interfacial processes in these liquids. The first step in any classical MD simulation is the validation of available force field models, the refinement of these if necessary, or the eventual development of a new one that is appropriate for reproducing the system’s properties under study. To the best of our knowledge, only Perkins et al.29 have recently developed and published an FFM for ethaline. As a starting point for this FFM, the authors chose the general Amber force field (GAFF) and adopted all the parameters for bonded terms of intramolecular potentials from it. Each ion and molecule was optimized in an isolated state. The final charges in EG molecules were obtained from a multiconformation

2. METHODS 2.1. Force Field Models. Various potential models previously proposed for the simulation of choline, chloride, and ethylene glycol were combined with each other to test possible force field models for ethaline. All are all-atom potentials, where each atom is represented by a Lennard-Jones (LJ) center and a fixed-point partial charge, with the latter leading up to integer charges (±1) on the Ch+ and Cl− ions. For the choline cation, two sets of partial charges were found in the literature. The first set is due to Sambasivarao and Acevedo,30 who determined the charges by averaging and weighting those initially obtained from an electrostatic potential fit based on ab initio energy-minimized conformations. The second set of partial charges used was proposed by Perkins et B


Article

The Journal of Physical Chemistry B al.31 The charges were obtained using the RESP charge derivation method. Both of these sets of charges were combined with bonded and nonbonded parameters from the classical OPLS-AA force field33,34 as well as with those from original papers.30,31 Regarding the chloride anion, two potentials were tested: the OPLS-AA and the one reported by Canongia Lopes et al.,35 where the intermolecular potentials were determined by fitting to the Born-Mayer potential developed for crystalline and molten salts. Three FFMs describing the ethylene glycol molecule were tested, namely, the OPLS-AA FF33,34 and the ones proposed by Szefczyk and Cordeiro32 and by Górny et al.36 In the FFM due to Szefczyk and Cordeiro,32 the charges initially generated using the CHELPG algorithm were then hand-tuned to reproduce the physical properties of liquid glycol under ambient conditions. The rest of the parameters were based on the OPLS-AA force field. Its LJ parameters were mostly maintained, except for the hydroxyl group oxygen van der Waals radius, whereas the bonded interactions were modified for several dihedrals. In the FFM by Górny et al.,36 the EG was described with a flexible model based on the CHARMM 27 force field, including intramolecular harmonic stretching, harmonic bending, torsional, van der Waals, and Coulombic terms. Table 1 summarizes all of the combinations of FFMs for Ch+, − Cl , and EG, which will be denoted as FFMn, where n is the number of each FFM. FFM1 differs from FFM2 only in the partial charges of the Ch+ cation, and all the rest of the parameters for all other particles in the mixture were taken from the classical OPLS-AA force field.33,34 In the case of FFM3, both partial charges and bonded and nonbonded parameters for Ch+ proposed by Sambasivarao and Acevedo30 were used. Cl− was described with the potential from OPLSAA, whereas for the EG molecule the force field model published by Szefczyk and Cordeiro32 was used. For FFM4 and FFM5, the partial charges and parameters (bonded and nonbonded) utilized for Ch+ were taken from Perkins et al.,31 and for Cl−, the Canongia Lopes et al.35 potential was used. These two FFMs differ solely with respect to the parameters applied for EG, in which the LJ and intramolecular potentials employed in FFM4 were established by Górny et al.,36 whereas the parameters used in FFM5 were proposed by Szefczyk and Cordeiro.32 In the last three force field models, FFM6−FFM8, the partial charges and parameters used for the cation were those reported by Sambasivarao and Acevedo,30 and for the anion, the parameters described by Canongia Lopes et al.35 were employed. Regarding the EG molecule, the potential from the OPLS-AA force field was used in FFM6, the force field described by Górny et al.36 was used in FFM7, and the parameters proposed by Szefczyk and Cordeiro32 were used in FFM8. All of he parameters for the first five FFMs (FFM1− FFM5) are shown in the Supporting Information (SI) in Tables S1−S17, as well as the schematic diagrams for choline cations and ethylene glycol molecules with the respective nomenclature used for each force field (Figures S1−S4). A schematic diagram of the choline cation, chloride anion, and ethylene glycol molecule, showing the entire atom types used in this work, is given in Figure 1. 2.2. Computational Details. MD simulations were carried out for the system composed of choline chloride and ethylene glycol (molar ratio of 1:2) using the GROMACS 5.0.4 software package.37

Figure 1. Geometries of the choline cation, chloride anion, and ethylene glycol molecule, along with the atom types used throughout this work. Nitrogen, carbon, oxygen, hydrogen, and chlorine atoms are represented by blue, cyan, red, white, and pink spheres, respectively.

The system containing 512 choline chloride ion pairs and 1024 ethylene glycol molecules was placed in a cubic box (dimensions of approximately 6 × 6 × 6 nm3). Initial configurations of the system were built with the PACKMOL package,38 ensuring a random distribution of molecules. To explore the temperature dependence of transport properties, a series of MD simulations were carried out at 313.15, 323.15, 333.15, 348.15, 358.15, and 373.15 K. The temperatures were selected by taking into account the available experimental data. Each temperature was sampled for 60 ns in an isothermal− isobaric (NpT) ensemble. The first 10 ns was regarded as equilibration, whereas the following 50 ns was used to derive properties. The temperature was kept constant by applying the Nosé−Hoover thermostat39 with a coupling constant of 0.5 ps, and the pressure was held fixed using the Parrinello−Rahman barostat40 with a coupling constant of 2.0 ps. In all runs, the Verlet leapfrog algorithm41 was used to integrate the equations of motion, with a time step of 2 fs. Constraints were enforced on all bond lengths using the LINCS algorithm.42 Periodic boundary conditions were applied in all three directions with long-range van der Waals tail corrections for both energy and pressure applied for distances longer than 1.1 nm. The same cutoff was used for the long-range electrostatic interactions in combination with the particle mesh Ewald (PME) technique.43 The neighbor list was updated using a Verlet cutoff scheme. The simulation boxes for surface tension and viscosity calculations were elongated in the z direction, and simulations were carried out in the NVT ensemble. 2.3. Calculation of Physicochemical Properties. Thermodynamic, transport, and structural properties for the ethaline system were computed from MD simulations to assess the quality of the applied force fields. All simulated properties were compared, whenever possible, with experimental data. 2.3.1. Density. ρ at different temperatures was obtained by averaging the results from 50 ns NpT simulations using GROMACS built-in tools. Pressure in these simulations was maintained at 1 bar. 2.3.2. Thermal Expansion Coefficient. αp was estimated from the variation of the density with the temperature (from 298.15 to 373.15 K), using the numerical derivative at constant pressure ⎛ ρ2 ⎞ ⎜ ln ρ1 ⎟ 1 ⎛ ∂ρ ⎞ 1 ⎛⎜ ∂V ⎞⎟ ⎟ ≈ −⎜ αp = − ⎜ ⎟ ≈ ρ ⎝ ∂T ⎠ V ⎝ ∂T ⎠ p ⎜ T2 − T1 ⎟ ⎝ ⎠p

()

3

(1)

−1

where V is the molar volume (cm mol ). 2.3.3. Isothermal Compressibility. βT was calculated from the following equation: C


Article

The Journal of Physical Chemistry B

Figure 2. Experimental (symbols) and simulated (lines) density for ethaline 1:2 (ChCl:2EG) at 1 bar. Experimental: black circles.2,10,15,17,20−24 Simulated lines: (a) black, FFM1; red, FFM2; green, FFM3; blue, FFM4; cyan, FFM5; magenta, FFM6; dark yellow, FFM7; and gray, FFM8; (b) black, 0.8FFM1; red, 0.8FFM2; green, 0.8FFM3; blue, 0.8FFM4; cyan, 0.8FFM5; magenta, 0.8FFM6; dark yellow, 0.8FFM7; and gray, 0.8FFM8.

⎛ ρ2 ⎞ ⎜ ln ρ1 ⎟ ⎟ βT = ⎜ ⎜ P2 − P1 ⎟ ⎝ ⎠T

in which Lz is the length of the box in the z direction, Pij is the ij component of the pressure tensor, and the factor of 1/2 stems from the presence of two interfaces in the system.32 2.3.7. Shear Viscosity. η was calculated in the framework of the nonequilibrium periodic perturbation method.46 A 100 ns NVT simulation was performed using the same box as for surface tension, with the same volume of vacuum, and a periodic acceleration ax(z) applied along the x axis

()

(2)

For this purpose, additional NpT runs of 10 ns each were performed at all considered temperatures and under a pressure of 100 bar. 2.3.4. Enthalpy of Vaporization. ΔHvap was calculated as follows44 ΔH vap = Epot(gas) − Epot(liquid) + RT

⎛ 2π ⎞ ax(z) = A cos⎜ z⎟ ⎝ Lz ⎠

(3)

with A being an arbitrary amplitude parameter. This parameter was fine-tuned manually. After a number of tests, the A value was set to 0.001 nm ps−2. The shear viscosity, η, was then calculated from the generated velocity profile (ν) computed during the simulation, as follows:

where Epot (gas) and Epot (liquid) are the total energies in the gas phase and in the liquid phase, respectively, R is the universal gas constant, and T = 298.15 K. 2.3.5. Self-Diffusion Coefficients. D was estimated from a 50 ns NpT simulation, following the Einstein relation31 D=

1 d lim ⟨|ri(t ) − ri(0)|2 ⟩ 6 t →∞ dt

η= (4)

d log(Δr 2) d log(t )

(5)

where Δr2 is the average MSD of the species. Self-diffusion coefficients were calculated in the region where β = 1, i.e., in the diffusive regime. 2.3.6. Surface Tension. γ was estimated from a simulation performed with the initial box, extended two times in the z direction, followed by vacuum of the same length. Such a system was equilibrated and simulated in the NVT ensemble for 10 ns, and the surface tension was calculated by the so-called virial route γ=

⎛ Pxx + Pyy 1 ⎜ Lz ⟨Pzz⟩− 2 ⎜ 2 ⎝

⎞ ⎟ ⎟ ⎠

Aρ ⎛ Lz ⎞ ⎜ ⎟ v ⎝ 2π ⎠

2

(8)

2.4. Ab Initio Calculations. All ab initio MD calculations, used to refine the nonpolarizable FFMs, were carried out with the CP2K program package,47 which provides a general framework for density functional theory (DFT) using a mixed Gaussian and plane waves (GPW) approach and classical pair and many-body potentials.48−51 The gradient-corrected exchange-correlation functional of Becke−Lee−Yang−Parr49,50 (BLYP) was employed. The wave functions were expanded in a double-ζ valence polarized (DZVP) Gaussian basis set52 for all atomic kinds. Analytical dual-space Goedecker−Teter−Hutter pseudopotentials53 were used to treat the core electrons. After a number of convergence tests, the electron density cutoff was set to 310 Ry and the plane wave cutoff of a reference grid covered by a Gaussian was set to 70 Ry. The pre-equilibrated configurations of 10 ethaline units (10 choline cations + 10 chloride anions + 20 ethylene glycol molecules) obtained in classical MD simulations using the original force field models were placed in a cubic cell in order to investigate the bulk behavior of the considered systems. The cell dimensions were calculated from the experimental density of ethaline at 313.15 K.15,19 Two short additional sequential equilibrations of 1 ps at elevated temperature (373.15 K), for the sake of faster dynamics, and at 313.15 K were conducted for each system. Then, five consecutive production runs of 2 ps each were performed at 313.15 K. During equilibration and production runs, the motion of nuclei follows Newton’s equation of motion and is propagated on the basis of the velocity Verlet algorithm

where ri is the vector coordinate of particle i at time t. D was obtained from the slope of the linear portion of the meansquare displacement (MSD) of atomic positions at long times. To ensure that the simulations were in the diffusive regime, the autocorrelation functions for two vectors in the studied systems were calculated. As suggested in the work of Perkins et al.,31 the N−O vector for the choline cation and the Og−Og vector for the ethylene glycol molecule (Figure 1) were selected as references. The β parameter was calculated to determine the location of the diffusive region, as given by45 β (t ) =

(7)

(6) D


Article


overcome this issue, the annealing method was applied, but no improvement was achieved. From this point on, the discussion related to the original FFMs will consider only the ones that reach diffusive regime. These are FFM1−FFM5 (Table 1). The simulated self-diffusion coefficients at 313.15 K are listed in Table 3, along with their relative deviations. As can be seen, the

with time steps of 0.5 and 1 fs. All of these ab initio MD simulations were performed in the NVT ensemble. A Nosé− Hoover thermostat54 with a chain length of 3 and a time constant of 50 fs fixed the temperature. Periodic boundary conditions were applied in all directions. Using the resulting wave function, the restrained electrostatic charge potential fit was generated for each system. For the fitting procedure of the RESP charges, the total charge of the system is constrained to zero, as is the sum of charges of ethylene glycol molecules. Then the charges on the cation and anion were optimized, starting from the following restriction: charges on each atom of the choline cation were restrained to their partial charges in the original FFM, and the charge on the chloride anion was restrained to −1. The obtained partial atomic charges were summed up on the cation and anion separately and averaged over the whole system. These total charges on the ions were further used as scaling factors so that only the absolute values of the original charges were changed, whereas all of their ratios remain constant.

Table 3. Self-Diffusion Coefficients Simulated with FFM1− FFM5 at 313.15 Ka

FFM1 FFM2 FFM3 FFM4 FFM5 FFM6 FFM7 FFM8

5.80 6.21 6.37 5.59 5.48 6.06 5.81 5.75

± ± ± ± ± ± ± ±

0.8FFM1 0.8FFM2 0.8FFM3 0.8FFM4 0.8FFM5 0.8FFM6 0.8FFM7 0.8FFM8

0.02 0.02 0.03 0.01 0.02 0.01 0.03 0.01

αp/10−4 K−1 7.22 7.28 7.67 6.29 6.90 6.51 6.26 7.07

± ± ± ± ± ± ± ±

FFM1 FFM2 FFM3 FFM4 FFM5

0.4 (±0.2) 0.5 (±0.1) 0.32 (±0.07) 0.19 (±0.08) 0.21 (±0.09)

91.60 90.13 93.28 96.01 95.59

DEG/1011 m2 s−1 1.05 1.01 0.82 0.45 0.39

(±0.3) (±0.2) (±0.1) (±0.07) (±0.06)

error/% 87.39 87.88 90.16 94.60 95.32

self-diffusion coefficients resulting from the FFM1−FFM5 models are systematically underestimated. The self-diffusion coefficient of the choline cation (DCh+) is underestimated at least 10 times (the case of FFM2) and sometimes even 20 times (in the case of FFM4) when compared to experimental values reported by D’Agostino et al.,23 who measured DCh+ = 4.76 × 10−11 m2 s−1 at 313.15 K. As for the self-diffusion coefficient of the ethylene glycol molecule (HBD donor), its self-diffusion coefficient in the best case is ∼8 times smaller (for FFM1) than the experimental value at the studied temperature (DHBD = 8.33 × 10−11 m2 s−1). The experimental values reported in the literature present a relative error not exceeding 2%. At this stage of analysis, it became clear that in order to predict dynamic properties accurately these FFMs require refinement to take into account the effect of electronic polarization. One known simple way of fighting the overestimation of electrostatic potential in the simulated system is to resort to scaled charge models in which the ions possess noninteger charges. As it has been shown recently,29,31,56−60 this scaling can considerably improve the dynamics of RTILs and DESs. To account for electronic polarization in the liquid ethaline, DFT calculations in the bulk phase were performed. A RESP fitting was employed and partial charges on ethaline atoms were generated so that they properly reproduce the potential obtained in DFT calculations. For all investigated FFMs, the average electrostatic charges were found to be very similar. Thus, we present here only charges for the system described by the FFM1 and FFM3 models because they represent the smallest and the largest values. The average electrostatic charge was found to be ±0.777 (±0.003) on both the cation and anion when the charges from FFM3 were used as a starting point for the RESP fit, whereas it was ±0.793 (±0.003) on the ions when the charges from FFM1 were used. On the basis of these results, we scaled the initial charges of each FFM, used in following study, by the rounded scaling factor of 0.8. Charges on EG were left unchanged. Thus, another electrostatic potential was assumed for ethaline at this point with the aim of decreasing the pairwise electrostatic interaction energy of the simulated system. As has been shown,56−61 this approach allows us to capture on average specific cation−anion interactions, including polarization. However, there are some known limitations of these FFMs, such as the nonpredictability of

temperature range 298.15−373.15 K FFM

error/%

Time interval used was 40−50 ns. Numbers in parentheses refer to standard deviations.

Table 2. Thermal Expansion Coefficients Obtained by Means of the Original and Refined FFMs αp/10−4 K−1

DCh+/1011 m2 s−1

a

3. RESULTS AND DISCUSSION 3.1. Original FFM Performance. The density values obtained for each tested FFM at the studied temperatures are compared to those obtained experimentally2,10,15,17,20−24 in Figure 2. As it can be seen from Figure 2a, the simulated densities obtained using the original FFMs are in good agreement with the experimental ones. The largest standard deviation of calculated densities is about 0.0003 g cm−3. The best reproducibility of the experimental measurements is observed for FFM7, at all temperatures, with an average error not exceeding 4%. These and the following reference to errors stand for (|Xexp − Xsim|)/Xexp × 100%, where X is the property obtained by experimental measurements, Xexp, or by simulation, Xsim. The thermal expansion coefficients (αp), calculated from the densities determined in simulations using the FFM1−FFM8 models, are shown in Table 2. The results obtained for the

FFM

FFM

0.01 0.02 0.02 0.01 0.01 0.02 0.03 0.01

original FFMs appear to be in a good agreement with the available experimental data because the experimental values of αp for ionic liquids are usually in the range of (5−7) × 10−4 K−1.2,15,17,19,20,22,55 Although densities and thermal expansion coefficients are in good agreement with experimental data, that is clearly not sufficient criteria for the validation of any FFM. Thus, selfdiffusion coefficients were analyzed. Unfortunately, for some FFMs the diffusive regime, which is necessary to calculate meaningful self-diffusion coefficients, was not attained, even after increasing the production time to 100 ns. In an attempt to E


Article


Table 4. Self-Diffusion Coefficients for Ethaline Simulated with Refined Force Field Models at Temperatures of 298.15−373.15 Ka DCh+/1011 m2 s−1

DEG/1011 m2 s−1

FFM

T/K

exp

sim

error/%

sim

sim

error/%

0.8FFM1

298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 313.15 313.15 313.15

2.62 4.76 7.50 8.95 ± 0.1 14.8 ± 0.8 19.5 ± 0.8 30.2 ± 0.5 2.62 4.76 7.50 9.6 ± 0.4 15.1 ± 0.4 19 ± 1 28 ± 1 2.62 4.76 7.50 8.5 ± 0.9 16 ± 1 21.1 ± 0.4 29.2 ± 0.7 2.62 4.76 7.50 5.5 ± 0.2 8.1 ± 0.5 12.0 ± 0.3 17.1 ± 0.1 2.62 4.76 7.50 4.8 ± 0.2 7.5 ± 0.3 10.4 ± 0.9 15 ± 1 4.76 4.76 4.76

1.76 ± 0.1 3.83 ± 0.02 6.4 ± 0.2 13.5 ± 0.7 22 ± 1 29.9 ± 0.3 39.7 ± 0.5 1.89 ± 0.06 4.2 ± 0.2 6.0 ± 0.2 13.63 ± 0.06 22.9 ± 0.9 28.3 ± 0.1 41 ± 1 2.08 ± 0.02 4.60 ± 0.06 6.5 ± 0.3 14.9 ± 0.7 21.9 ± 0.8 29.6 ± 2 42 ± 1 1.0 ± 0.1 2.57 ± 0.02 3.9 ± 0.2 8.5 ± 0.3 14 ± 1 18 ± 1 27.5 ± 0.7 0.84 ± 0.05 1.91 ± 0.08 3.1 ± 0.1 7.2 ± 0.5 12.5 ± 0.3 16.2 ± 0.8 25.1 ± 0.7 1.77 (±0.06) 1.46 (±0.1) 1.26 (±0.08)

32.8 19.5 14.7 17.0 ± 0.6 27.5 ± 0.4 34.9 ± 0.3 50.0 ± 2 27.9 11.8 20.0 16.3 ± 0.3 27.4 ± 0.8 35 ± 1 51.5 ± 0.9 20.6 3.4 13.3 17.46 ± 0.08 28.6 ± 0.3 36.5 ± 0.8 52 ± 1 62 46 48 9.8 ± 0.2 15.7 ± 0.2 20.2 ± 0.2 29.2 ± 0.3 68 60 59 9.1 ± 0.2 13.5 ± 0.5 18.29 ± 0.01 28 ± 2 62.82 69.33 73.53

3.0 ± 0.3 6.8 ± 0.5 10.2 ± 0.4

4.77 8.33 13.2

3.8 ± 0.3 8.3 ± 0.2 12.0 ± 0.4

20.3 0.4 9.1

3.06 ± 0.01 6.8 ± 0.4 9.1 ± 0.9

4.77 8.33 13.2

3.85 ± 0.09 8.1 ± 0.2 12.2 ± 0.6

19.3 2.8 7.6

3.2 ± 0.1 6.7 ± 0.5 9.5 ± 0.2

4.77 8.33 13.2

4.15 ± 0.04 8.6 ± 0.1 12.7 ± 0.2

13.0 3.2 3.8

1.90 ± 0.05 4.4 ± 0.4 6.4 ± 0.3

4.77 8.33 13.2

2.14 ± 0.01 4.38 ± 0.07 7.12 ± 0.08

55.1 47.4 46.1

1.51 ± 0.08 3.3 ± 0.1 5.28 ± 0.06

4.77 8.33 13.2

1.84 ± 0.02 3.7 ± 0.1 5.84 ± 0.01

61.4 55.6 55.8

8.33 8.33 8.33

3.07 (±0.1) 3.15 (±0.1) 2.69 (±0.08)

0.8FFM2

0.8FFM3

0.8FFM4

0.8FFM5

0.8FFM6 0.8FFM7 0.8FFM8 a

DCl−/1011 m2 s−1 exp

63.15 62.18 67.71

Time interval used was 10−50 ns. Numbers in parentheses refer to standard deviations.

phase separation in IL mixtures.62 This seems to be due to an underestimation of the IL cohesive energy, which in turn may lead to poor performance whenever these are used to predict the phase behavior of IL-based mixtures.63 Nevertheless, a good performance of FFMs refined in this way for the simulation of bulk properties of pure ionic liquids is obtained.29,31,56−62 This, together with the higher computational affordability of FFM without an explicit interaction term responsible for electronic polarization effects, leads to nonpolarizable FFMs with a refined electrostatic potential being widely used. 3.2. Refined FFMs Performance. The density values obtained for all of the refined FFMs (e.g., with partial atomic charges on both the cation and anion scaled by 0.8) over a wide range of temperature are shown in Figure 2b. As expected, the resulting densities are slightly lower compared to those obtained with the original FFMs because the reduced charges lead to a decrease in attractive interactions. However, the

maximum standard deviation of those does not exceed 0.0001 g cm−3. As can be seen in Figure 2b, results closest to the experimental density are obtained with 0.8FFM5. The largest deviation in the simulated data from experimental values is observed for 0.8FFM3 (average error of 11%). In addition, the thermal expansion coefficients (αp) derived from MD simulations with the refined FFMs (Table 2) are in good to excellent agreement with experimental data. Self-diffusion, although sometimes tricky to calculate, allows a general perception of dynamics in the model system and therefore is useful for FFM validation. The self-diffusion coefficients obtained by using refined FFMs at 298.15−373.15 K (for 0.8FFM6−0.8FFM8 at 313.15 K only) are gathered in Table 4, along with their standard deviations and respective errors. Clearly, the self-diffusion coefficients are much higher and closer to experimental data than the ones obtained with FFMs before refinement. The diffusive regime attainment was F


Article

The Journal of Physical Chemistry B reached for all of the refined FFMs; therefore, D values were calculated in all cases. The diffusion regime attainment was controlled by autocorrelation functions, as described previously. Examples of autocorrelation curves, used to determine whether the diffuse regime was achieved, are shown in Figures S5−S14 of the SI for several refined FFMs. By analyzing the results summarized in Table 4, one can see that the first three FFMs reproduce the experimental diffusion coefficients better. On the other hand, the performance of the 0.8FFM6−0.8FFM8 models is quite poor compared to that of the other ones. Therefore, only models 0.8FFM1−0.8FFM5 are considered for further discussion. Generally, a slight underestimation of self-diffusion coefficients by all of these refined FFMs can be observed. 0.8FFM1−0.8FFM3 present a somewhat better reproduction of experimental D±, with a maximum average difference of 16% for 0.8FFM1. The minimal difference between the experimental and calculated diffusion is registered for 0.8FFM3, with an average difference of only 10% (average of 12% for choline, Ch+, and 7% for ethylene glycol, EG). To the best of our knowledge, there is no reported selfdiffusion coefficient for Cl− in deep eutectic solvents. Here, it should be noticed that the values for self-diffusion coefficients obtained in this work are closer (average errors of 10%) to experimental ones than those calculated with the previously developed FFM (average errors of 20%) for ethaline.29 As expected, the choline cation diffuses more slowly than the hydrogen bond donor molecule, reflecting the trends in molecular size and weight. This observation is in agreement with the NMR findings reported by D’Agostino et al.23 The ratio between the self-diffusion coefficients of HBD/ cation and anion/cation decreases with temperature, displaying average decreases of 16 and 13%, respectively. This reduction indicates a decrease in the interactions among these species (Figures S15 and S16 in the SI). Interestingly, the ratio between the self-diffusion coefficients of anion/HBD remain constant with the increase in temperature (Figure S17 in the SI), suggesting that the anion/HBD interactions are not affected by temperature. As shown by D’Agostino et al.,23 the temperature dependence of self-diffusion coefficients can be described by an Arrhenius exponential function. To verify if our results can be described with the same type of function, Arrhenius dependences for DCh+ and DEG obtained from 0.8FFM3 were plotted and are depicted in Figure 3. Because the dependence is linear in the range of temperature studied (298.15−373.15 K), the activation energies for the diffusion of the choline cation and ethylene glycol can be derived. The calculated activation energies (Ea) for each species are presented in Table 5, along with the experimental NMR values measured for ethaline 1:2 (ChCl:2EG).23 Evidently the Ea values calculated in this work are very similar to those found experimentally, with a maximum difference of 11%. Also, it can be seen that the activation energy for the diffusion of the choline cation is higher than that of HBD, which is in accordance with D± values, as discussed above. The simulated enthalpies of vaporization, ΔHvap, of ethaline are listed in Table 6 for three possible compositions in the gas phase. As a first hypothesis, the gas phase constituted of one choline chloride ion pair and two ethylene glycol molecules (ChCl:2EG, one unit of DES ethaline) was assumed. Then, second and third guesses were devised in which the gas phase was composed by one molecule of ethylene glycol or one ion pair of choline chloride, respectively. As can be seen from Table

Figure 3. Arrhenius plot of the self-diffusion coefficients for choline (Ch+) (black) and ethylene glycol (EG) (red) in ethaline, derived from MD simulations using 0.8FFM3.

Table 5. Diffusion Activation Energies for Choline (Ch+) and Ethylene Glycol (EG) FFM 0.8FFM1 0.8FFM2 0.8FFM3

Ea/kJ mol−1a

Ea/kJ mol−1b

error/%

± ± ± ± ± ±

37.1 34.6 37.1 34.6 37.1 34.6

5 10 11 8 11 10

+

35 31 33 32 33 31

Ch EG Ch+ EG Ch+ EG

1 1 1 1 1 1

a Obtained in this work using the refined FFMs. bObtained from the experimental data reported in ref 23.

Table 6. Enthalpies of Vaporization for Ethaline Simulated with the Refined Force Field Models at 298.15 K ΔHvap/(kJ mol−1) FFM 0.8FFM1 0.8FFM2 0.8FFM3 0.8FFM4 0.8FFM5

ChCl:2EG 175.5 173.6 167.5 179.5 170.9

± ± ± ± ±

0.6 0.3 0.3 0.5 0.5

EG 205.7 234.6 133.7 314.29 324

± ± ± ± ±

ChCl 0.6 0.6 0.2 0.07 1

245.8 232.5 289.6 284.9 307.9

± ± ± ± ±

0.2 0.3 0.1 0.3 0.2

6, the enthalpy of vaporization for the ChCl ion pair is generally higher than the ΔHvap values for the other species, independent of the FFM used. This is an expected result because ChCl is a nonvolatile solid salt at room temperature. Regarding the vaporization enthalpies of ChCl:2EG, one would expect these values to be higher than the enthalpy of vaporization of pure EG. This is not observed for all of the tested FFMs, except for 0.8FFM3, which output values of ΔHvap for EG smaller than for ChCl:2EG. However, ΔHvap(EG) and ΔHvap(ChCl:2EG) are not too different, which leads us to expect that ethaline will be composed of both species in the gas phase. The results obtained with the remaining force fields seem to be less consistent. Unfortunately, no experimental values for the enthalpy of vaporization of ethaline have been reported in the literature. Nevertheless, for RTILs the vaporization process is well studied and ΔHvap values are known. Generally, the ΔHvap values for RTILs are in the interval of 100 and 200 kJ mol−1.64,65 Although there are no experimental data to compare with, the computed values for the enthalpy of vaporization of ethaline G


Article


0.8FFM2), the overall differences are not large, considering that the isothermal compressibility is a second-order property. As for the three FFMs providing better self-diffusion coefficients (0.8FFM1−0.8FFM3, Table 4), model 0.8FFM3 reproduces an isothermal compressibility with the smallest difference of 17%. As can be seen, surface tension values reproduce the experimental data very accurately, with a lower mean difference of 5% obtained for model 0.8FFM3. Viscosity is a property that is difficult to estimate, both experimentally and computationally, because many ionic liquids have relatively high viscosity compared to organic solvents. Thus, an accurate prediction of the viscosity of deep eutectic mixtures is computationally challenging.32,66,67 Our simulated results obtained from refined force field models are presented in Figure 4 as a function of temperature. For the refined FFMs performing better for all of the other properties (0.8FFM1− 0.8FFM3), the viscosity values are on average a factor of 2.3 (for both 0.8FFM2 and 0.8FFM3) to 2.5 (for 0.8FFM1) lower than the corresponding experimental values,23,25 which is in agreement with the simulated data for other ionic liquids when using nonpolarizable force field models.67 To the best of our knowledge, there are no reported results of viscosity for any DES from MD simulations. To quantitatively compare the trends in the temperature dependence of the viscosity of ethaline determined from simulation and experiment,25 both sets of values were normalized by the viscosity at the lower temperature (298.15 K). As can be seen in Figure 4b, the normalized trends in the temperature dependence of the viscosity resulting from both setsexperimental and simulationsare in very good agreement. Considering the overall properties obtained with refined force field models 0.8FFM1−0.8FFM5, it can be concluded that 0.8FFM3 better reproduces the experimental data (i.e., presents lower mean difference for most of the properties: 10% for self-diffusion coefficients, 17% for isothermal compressibility, and 5% for surface tension). This conclusion is still valid when considering the shear viscosity results. Additionally, 0.8FFM3 is the only FF from all of those tested here that provides values of ΔHvap for EG and ChCl:2EG within the expected order. Unfortunately, there are no experimental studies on the enthalpy of vaporization of ethaline, which can support or disprove our hypothesis of the presence of both EG and ChCl:2EG in ethaline in the gas phase. This type of experimental study would be very welcome. 3.3. Structural Analysis. To probe the underlying structure of ethaline, radial distribution functions (RDFs) and corresponding coordination numbers (CNs) were computed and hydrogen bonds (H-bonds) were analyzed. Taking into account the results obtained for physical and chemical properties, discussed above, for our further study, only 0.8FFM1−0.8FFM3 will be discussed. Aiming to prove that the suggested refinement does not cause significant changes in the liquid structure, we compared the center of mass (COM) cation−cation RDFs obtained before and after the refinement. As shown in Figure S18 (SI), the electrostatic charge scaling only slightly affects the positions (∼0.02 nm) and heights (∼0.5) of the RDF peaks, attributable to a marginally looser packing arrangement due to the lower density. Thus, it is safe to conclude that the structure of the DES is preserved after applying the scaling factor of ±0.8e. A comparison of three best-performing FFMs was initiated by analyzing the dihedral angle distributions for the species present in the system. Note that only the most interesting

obtained with the refined FFMs and considering the composition of the gas phase to be the ChCl ion pair plus two molecules of EG, appear to be reasonable because they are a similar order of magnitude (Table 6). Table 7 contains the isothermal compressibility (β) and surface tension (γ) for ethaline derived from MD simulations at Table 7. Isothermal Compressibility and Surface Tension for Ethaline Simulated Using the Refined FFMs at Temperatures of 298.15−373.15 K, with the Respective Standard Deviations FFM

T/K

0.8FFM1

298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15 298.15 313.15 323.15 333.15 348.15 358.15 373.15

0.8FFM2

0.8FFM3

0.8FFM4

0.8FFM5

βT/105 bar−1

error/%

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

28.6a 12.5a 22.2a

3.7 3.4 4.0 3.9 4.39 4.35 4.63 3.3 3.7 4.1 4.1 4.2 4.47 4.55 3.3 3.44 4.0 4.17 4.6 4.76 5.01 2.47 2.76 2.9 3.26 3.31 3.2 3.63 2.57 2.87 3.02 3.3 3.5 3.42 3.7

0.1 0.1 0.2 0.2 0.07 0.08 0.07 0.3 0.1 0.1 0.1 0.1 0.06 0.06 0.2 0.08 0.1 0.09 0.1 0.06 0.06 0.03 0.08 0.1 0.09 0.09 0.1 0.07 0.02 0.02 0.02 0.1 0.1 0.06 0.1

14.1a 25.9a 26.0a

17.2a 14.9a 18.9a

14.5a 5.3a 3.8a

11.0a 1.5a 0.1a

γ/mN m−1

error/%

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

18.3b 1.7b 5.0b

58 46.7 49 47 45 43 40.9 54 50 48 47 44 44 42 48 47 45 44 42.8 41.3 40 55 53 56 56 51 54 53 55 57 46 52 48 46

5 0.4 3 1 1 2 0.8 6 1 2 2 2 2 1 3 3 6 2 0.7 0.5 1 8 3 4 2 4 2 2 4 4 3 2 2 1

9.4b 5.5b 1.9b

3.5b 4.9b 5.7b

12.2b 12.4b 19.4b

12.1b 19.3b 1.3b

a Calculated using as reference values βT(298.15) = (2.89 ± 0.3) × 10−5, βT(313.15) = (2.91 ± 0.3) × 10−5, and βT(323.15) = (3.02 ± 0.3) × 10−5 bar−1, taken from ref 19. bCalculated using as reference values γ(298.15) = 48.91 ± 0.1, γ(313.15) = 47.50 ± 0.1, and γ(323.15) = 46.67 ± 0.1 mN m−1, taken from ref 21.

298.15−373.15 K, using the refined 0.8FFM1−0.8FFM5. The obtained isothermal compressibility values are slightly overestimated compared to the experimental values. Because β is derived from the density, one would expect to obtain it with lower relative error using FFMs that better predict density. Indeed, as it can be seen from Table 7 and Figure 2b, 0.8FFM5 allows better reproducibility of both the density and isothermal compressibility. Although discrepancies in β values for other FFMs are higher (with the highest average difference of 22% for H


Article


Figure 4. (a) Experimental and simulated viscosity for ethaline 1:2 (ChCl:2EG) at 1 bar. (b) Experimental and simulated normalized viscosity (normalized by the viscosity at temperature 298.15 K). Experimental: black circles.23,25 Simulated: red squares, 0.8FFM1; green up triangles, 0.8FFM2; blue down triangles, 0.8FFM3; cyan right triangles, 0.8FFM4; magenta left triangles, 0.8FFM5.

Figure 5a, the distribution, produced by 0.8FFM1 (red line), indicates the coexistence of two main EG conformations in the ethaline medium, cis (∼0°) and trans (±180°), whereas the distribution, resulting from 0.8FFM3 (black line), suggests the prevalence of gauche conformations (±50°). As has been demonstrated in ref 68, the ethylene glycol fluid, when composed mainly by either trans or gauche conformers, yields simulated self-diffusion coefficients closer to experimental values. Our results are consistent with this observation because the gauche conformation is predominant for 0.8FFM3 and this model permits us to obtain the closest to experimental selfdiffusion coefficients (Table 4). Interestingly, the ethylene glycol Og−Cg−Cg−Og dihedral angle distribution resulting from 0.8FFM3 and the one reported in the original study of pure EG in the condensed phase32 are quite different. In fact, Szefczyk et al.32 obtained a prevalence for angles of 0 and ±150°, distributed almost equivalently. This difference in the angle distribution for pure EG and for ethaline is a direct consequence of the influence of the medium because in our study the EG molecules are surrounded by strongly charged choline and chloride ions. As for the choline cation, the obtained N−Cs−Cw−O angle distributions are quite similar for 0.8FFM1 and 0.8FFM3 (Figure 5b), presenting peaks at ±180 and ±90°. However, for the 0.8FFM1 model, the peaks corresponding to ±180° angles are much sharper and the ±90° ones are very weakly pronounced, which is clear evidence of a reduced rotational mobility for the Cs−Cw−O−Ho tail in this model. Indeed, the N−O distance distribution, shown in the inset of Figure 5b, for the 0.8FFM1 model consists only of two lengths, whereas for 0.8FFM3 the length distribution is much broader. The completely extended N−Cs−Cw−O angle (predicted in the case of 0.8FFM1) would be an expected behavior in gas-phase calculations without solvation effects. Thus, the flexible N−Cs− Cw−O angle (the case of 0.8FFM3) aligns better with simulations in the ethaline medium. The radial distribution functions obtained for the COM of the different species in the system at 298.15 K, by using 0.8FFM3, are shown in Figure 6. The RDFs obtained for anion−anion and anion−HBD, with the first three FFMs at two temperatures, 298.15 and 373.15 K, are also presented in Figure 6. As can be seen, the increase in temperature results in a decrease in molecular interactions between the various species present in the DES. This is reflected, for example, in the increase in anion−anion distances (positions of maxima, Figure 6b) and the decreasing intensity of the RDF curves (Figure 6b,c), related to looser contacts between the species.

conformation distributions are discussed. Figure 5 contrasts the distributions for the Og−Cg−Cg−Og angle of the HBD

Figure 5. Probability distribution for the (a) Og−Cg−Cg−Og angle of ethylene glycol and the (b) N−Cs−Cw−O angle of the choline cation in ethaline obtained with 0.8FFM1 (red) and 0.8FFM3 (black).

molecule (Figure 5a) and the N−Cs−Cw−O angle of the choline cation (Figure 5b), obtained by means of 0.8FFM1 and 0.8FFM3 models. Results for 0.8FFM2 are not shown because they are similar to those found with 0.8FFM1. Notoriously, the EG molecule possesses two stable configurations of the angle between hydroxymethyl groups: trans, with an Og−Cg−Cg−Og dihedral angle of 180°, and gauche, where the angle is equal to ±60°.68,69 As can be seen in I


Article


Figure 6. COM RDFs for ethaline derived from MD simulations using 0.8FFM3 (a): cation−cation (black), anion−cation (red), cation−HBD (green), and HBD−HBD (blue). Temperature dependence of COM RDFs for (b) anion−anion and (c) anion−HBD, derived from MD simulations using different FFMs: black, 0.8FFM1 at 298.15 K; red, 0.8FFM2 at 298.15 K; green, 0.8FFM3 at 298.15 K; blue, 0.8FFM1 at 373.15 K; cyan, 0.8FFM2 at 373.15 K; and magenta, 0.8FFM3 at 373.15 K.

Force field models 0.8FFM1 and 0.8FFM2 have structural properties very close to each other, at both temperatures for all of the COM RDFs analyzed (Figure 6b,c). As it can be noticed in Figure 6b, the anion−anion curves for these two FFMs, unlike the same curve for 0.8FFM3, present a shoulder on the first peak. The form of Cl−Cl RDFs calculated using 0.8FFM3 seems to be more logical because these curves correspond to an interaction of one-atom anions, and as such, one would expect a more equidistant distribution. In summary, it can be concluded that 0.8FFM3 presents more plausible structural behavior and reproduces experimental observations better compared not only to all FFMs tested here but also to the FFM for ethaline reported previously.29 This, allied with 0.8FFM3’s compatibility with any other popular force field models, leads us to confirm that it is the best force field model for DES ethaline from the FFMs evaluated in this work. Therefore, the following discussion will be focused on exploring the structural properties resulting from 0.8FFM3. The RDFs depicted in Figure 6 for the 0.8FFM3 model show important anion−cation, anion−HBD, and HBD−HBD contacts because the distances between the COMs of these species are less than 0.5 nm. The chloride anion displays a bimodal distribution with respect to the choline cation and HBD molecule. For the anion−cation RDF curve, two peaks appear at distances of 0.47 and 0.54 nm, and for the anion− HBD RDF curve, the two peaks are located at 0.38 and 0.45 nm. The first maximum peaks of cation−HBD are located at 0.55 nm (i.e., further than the anion−cation separation (0.47 nm) but closer than that of the cation−cation (0.67 nm)). These results suggest that ethylene glycol molecules enter only into the second layer around the choline cations, and the cation−cation interactions are insignificant. As for COM EG− EG RDFs, one can see that the EG molecule presents two contact peaks with respect to another HBD molecule, with the first peak located at 0.45 nm and the second located at 0.50 nm.

These distances are positioned between the anion−HBD (0.38 nm) and the cation−HBD (0.55 nm) distances, which leads us to conclude that in the first layer around the EG are the chloride anions, followed by the EG molecules in the second layer. This positioning suggests that interactions between cation−HBD are weaker than those between anion−HBD. Taking into account the close contact between the anion and the HBD molecule (Figure 6c), it is important to discuss these interactions in detail. Figure 7 depicts possible relative

Figure 7. Snapshots with possible conformations and corresponding distances (in nm) between anion−HBD pairs. Cyan spheres represent anion (Cl−); pink spheres represent center of mass of the ethylene glycol (HBD) molecule.

orientations of the ethylene glycol molecule with respect to the anion. Various possible conformations of ethylene glycol at two main distances from the chloride anion, matching the ones obtained in the RDF (Figure 6c), are shown. One chloride anion can be observed at distances of, for example, 0.37 and 0.41 nm from two EG molecules in the gauche conformation as well as at distances of 0.36 and 0.45 nm from two EG molecules in the trans conformation. These distance values are close to those observed in the anion−HBD COM RDF curve (0.38 and 0.45 nm). Thus, no correlation between Cl−HBD distances and HBD conformations can be established. Moreover, it is observed that there are many more distances larger and smaller than 0.4 nm, which is consistent with the COM RDFS presented. J


Article

The Journal of Physical Chemistry B Additional information about the structure of this deep eutectic solvent can be provided by atom−atom radial distribution functions. As mentioned previously, DESs form a eutectic mixture mainly through H-bonding interactions, so the analysis of such close contacts is of crucial importance. As can be seen in Figure 8, ethaline presents two close contacts between the chloride anion and the hydroxyl hydrogen

Figure 9. RDFs for oxygen atoms of ethaline, derived from MD simulations using 0.8FFM3. Black line corresponds to Og−Ho and red to Og−Hg atom pairs (intermolecular interactions). The inset shows the total contributions for Og−Hg (green) and relative inter- (red) and intramolecular (blue) contributions for the same atom pair.

When looking at the Og−Hg interactions, it becomes evident that they can occur in the same molecule. These intramolecular interactions are accounted for in the hydrogen bonding analysis because they are not covalently bonded as the distance is much smaller (0.094 nm). Although the increase in temperature does not affect significantly the average number of neighbors for the Cl−N pair, this is not the case for all atom pairs. Figure S20 in the SI shows the average coordination number (CN) for the oxygen atoms and for the anion around hydroxyl groups. For all RDFs, there is in general a slightly reduction of the peak heights with the increase in temperature, reflecting the weakening of the interatomic interactions (Figure S20). Concerning CNs, those for pairs Og−Ho and Cl−Ho remain constant, regardless of the temperature change. This does not occur for the Og−Hg and Cl−Hg pairs, for which the respective CN decreases and increases (>15%) with the temperature increase. It is interesting that the decrease in CN with temperature occurs only when a greater decrease in the RDF peak heights is observed, for instance, for the Og−Hg interaction. The above indicates the existence of a somewhat rigid solvation shell for the Cl−Hg interaction. Further hydrogen bond analysis was made using the criterion proposed by Luzar and Chandler.70 According to this criterion, a hydrogen bond exists if the distance between the participating atoms is less than 3.5 Å and simultaneously the X−H····Y angle is less than 30°, where X is the hydrogen donor atom, H is the hydrogen atom, and Y is the hydrogen acceptor atom. Table 8 contains the mean relative contribution of hydrogen bonds that are present in ethaline along the trajectory analyzed. It can be seen from Table 8 that, contrary to Perkins et al.’s29 findings, the higher contribution to hydrogen bond interactions occurs between ethylene glycol molecules. This is not surprising because these molecules are present in a larger quantity. Nevertheless, it is important to notice that the interaction between the anion and HBD also presents large relative contributions for hydrogen bonds (24.8% at 298.15K and 29.3% at 373.15 K), even larger than those for anion− cation. In principle, this would be contradictory to what is obtained by RDFs because Cl−Ho presents a higher first peak than that obtained for Cl−Hg. However, we have to be aware that there is another parameter, in addition to the distance, that

Figure 8. RDFs for ethaline derived from MD simulations using 0.8FFM3. Black line corresponds to anion−cation, red to Cl−Ho, green to anion−HBD, and blue to Cl−Hg. The inset shows a comparison between the RDFs obtained for anion−cation (black) and Cl−N (magenta).

atoms. One interaction, Cl−Ho (Figure 1), is referred to the interaction with choline (Ch+), and the other, Cl−Hg (Figure 1), is established with either hydrogen atom of the two hydroxyl groups of ethylene glycol. In addition, it can be seen that the differences between the COM and atom−atom RDFs are mainly due to the position of the center of mass, which may coincide with the atom positions. For instance, in the case of the anion−cation and anion−HBD RDFs, the closest contact occurs only at a distance of 0.38 nm (Cl−HBD), whereas when considering Cl−Ho and Cl−Hg the interactions occur at a much smaller distance of 0.23 nm in both cases. At the same time, when one compares the anion−cation and the Cl−N no significant differences between these two RDFs can be noticed because the center of mass of choline is very close to the nitrogen atom. This is also confirmed by the average number of neighbors in the first coordination shell that can be estimated by integrating the corresponding RDFs up to their first minimum. As can be seen in Figure S19 of the SI, the RDF Cl−N yields an average coordination number of about 4, independent of the temperature. Let us discuss the closest contacts of the oxygen atoms of ethylene glycol (Og) with the hydroxyl hydrogen atoms of both ethylene glycol (Hg) and choline (Ho) (Figure 9). The Og− Hg and Og−Ho RDFs exhibit strong first peaks at 0.17 and 0.18 nm, respectively, indicating strong interactions between these atoms. The inset in Figure 9 presents the relative contributions from inter- and intramolecular interactions between the oxygen atoms and hydroxyl hydrogen atoms of ethylene glycol. It can be observed that at 0.17 nm both contributions are present, although it appears that there is a slightly larger fraction of intermolecular interactions than intramolecular ones. The same tendency is observed at larger distances: the second smaller peak located at 0.33 nm presents stronger intermolecular interactions than intramolecular ones. K


Article


Table 8. Average Number and Mean Relative Contribution of Hydrogen Bonds Derived from MD Simulations Using 0.8FFM3 at 298.15 and 373.15 K no. of hydrogen bonds

% hydrogen bonds

interaction

type of hydrogen bond

298.15 K

373.15 K

Δno. of H‑bonds

298.15 K

373.15 K

cation−cation cation−HBD

O−Ho···O O−Ho···Og Og−Hg···O O−Ho···Cl Og−Hg···Cl Og−Hg···Og

2.3 161.0 5.9 157.1 303.4 594.3 1224.0

3.2 134.0 7.5 130.8 303.5 458.2 1037.2

0.9 −27.0 1.6 −26.3 0.1 −136.1 −186.8

0.2 13.1 0.5 12.8 24.8 48.6 100

0.3 12.9 0.7 12.6 29.3 44.2 100

anion−cation anion−HBD HBD−HBD total

4. CONCLUSIONS In this work, we report a detailed analysis of existing force field models for choline, chloride, and ethylene glycol, which are components of the deep eutectic mixture ethaline. A series of parameters for ethaline were tested, and electrostatic charges for the respective force fields were explored with ab initio MD calculations of the bulk system. As a result, refined classical nonpolarizable FFMs, which indirectly account for charge transfer in the system, were proposed. After the initial validation of the refined FFMs by a large number of transport and static properties, such as density, expansion coefficients, self-diffusion coefficients, enthalpies of vaporization, isothermal coefficients, surface tensions, and viscosity over a wide range of temperature (298.15−373.15 K), particular attention was paid to the structural properties of the solvent. The structure of the DES was characterized via radial distribution functions, the running coordination number, and hydrogen bond analysis. Such an analysis allowed us to establish the most important anion−cation, anion−HBD, and HBD−HBD interactions in bulk ethaline. The average coordination number for chloride in its well-defined first solvation shell with choline was found to be equal to 4 at all studied temperatures. The hydrogen bond structure in ethaline system is formed mainly by EG−EG Hbonds (∼50% of all H-bonds) and Cl−EG H-bonds (∼25%). The temperature dependence of the liquid’s structure reveals that most of the interactions become weaker with the increase in temperature, with the exception of the close contact of anion−HBD, which has a somewhat rigid structure leading to an independent temperature steady state. In summary, the refined force field model proposed here for ethaline (0.8FFM3) provides excellent agreement between simulated and experimental data on physicochemical and structural properties, and it is fully compatible with the OPLS all-atom force field. Comparison with other available force fields shows that in fact progress was achieved and that the presented FFM can predict experimental data better, whereas existing FFMs of the same kind underestimate the transport properties. Simulated values of shear viscosity for the eutectic mixture ethaline, as well as surface tensions and isothermal compressibility values, are presented for the first time.

should be considered when studying hydrogen bonds, that is, the cutoff angle. Thus, hydrogen bonds cannot be evaluated only from RDF curve analyses. As the temperature increases, a clear decrease in the following hydrogen bond interactions can be observed (Table 8): cation−HBD, anion−cation, and HBD−HBD, suggesting that these interactions are less pronounced at higher temperatures. This observation is consistent with our abovementioned findings regarding the temperature dependence of the RDFs and CNs (Figure S20), where the larger difference in RDF peaks, upon increasing temperature, led to a larger CN decrease (Og−Hg). These observations are also in agreement with the self-diffusion coefficient behavior. Indeed, the ratio of the self-diffusion coefficients of ethylene glycol/choline cation and chloride anion/choline cation decreases with increasing temperature (Table 4), indicating a decrease in molecular interactions between these two species. The hydrogen bond interactions for HBD−anion remain constant with increasing temperature, which can explain the invariance of the ratio of the self-diffusion coefficients with temperature. A general overview of the variations observed upon increasing temperature is shown in Table 9 for the most important interactions. Table 9. Variations of the Different Type of Interactions Involving Hydrogen upon Increasing the Temperature from 298.15 to 373.15 K interactions method

Og−Ho

Og−Hg

Cl−Ho

Cl−Hg

RDF peaks CN Δno H‑bonds

decrease no change decrease

decrease decrease decrease

decrease no change decrease

no change increase no change

Interestingly, the Og−Ho and Cl−Ho hydrogen bonds behave in a similar manner: even the numbers of hydrogen bonds destroyed with the increase in temperature are virtually the same (i.e., 27 and 26, respectively (Table 8)). In the case of the Og−Hg interaction, where a decrease in CN is observed, the decrease in the number of hydrogen bonds is the most significant. Predictably, the most remarkable results are obtained for the Cl−Hg interaction: although RDF peaks do not represent any significant variation, an increase in the coordination numbers is evident, as is a constant number of hydrogen bonds. This evidence reinforces the idea that complex interactions involving the HBD are present in deep eutectic mixtures such as ethaline.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b07233. The parameters of the force fields tested in this study and presenting better results (FFM1−FFM5) and autocorrelation functions in the 298.15−373.15 K range, used to verify the diffusive regime attainment (PDF) L


Article


■

Chloride/Malonic Acid) Systems. J. Chem. Thermodyn. 2014, 68, 216−220. (14) Nkuku, C. A.; LeSuer, R. J. Electrochemistry in Deep Eutectic Solvents. J. Phys. Chem. B 2007, 111, 13271−13277. (15) Shahbaz, K.; Baroutian, S.; Mjalli, F. S.; Hashim, M. A.; AlNashef, I. M. Densities of Ammonium and Phosphonium Based Deep Eutectic Solvents: Prediction Using Artificial Intelligence and Group Contribution Techniques. Thermochim. Acta 2012, 527, 59−66. (16) Figueiredo, M.; Gomes, C.; Costa, R.; Martins, A.; Pereira, C. M.; Silva, F. Differential Capacity of a Deep Eutectic Solvent Based on Choline Chloride and Glycerol on Solid Electrodes. Electrochim. Acta 2009, 54, 2630−2634. (17) Yadav, A.; Kar, J. R.; Verma, M.; Naqvi, S.; Pandey, S. Densities of Aqueous Mixtures of (Choline Chloride + Ethylene Glycol) and (Choline Chloride + Malonic Acid) Deep Eutectic Solvents in Temperature Range 283.15−363.15 K. Thermochim. Acta 2015, 600, 95−101. (18) Haerens, K.; Matthijs, E.; Chmielarz, A.; Van der Bruggen, B. The Use of Ionic Liquids Based on Choline Chloride for Metal Deposition: A Green Alternative? J. Environ. Manage. 2009, 90, 3245− 3252. (19) Leron, R. B.; Li, M.-H. High-Pressure Volumetric Properties of Choline Chloride−Ethylene Glycol Based Deep Eutectic Solvent and Its Mixtures with Water. Thermochim. Acta 2012, 546, 54−60. (20) Ghosh, S.; Ryder, K.; Roy, S. Electrochemical and Transport Properties of Ethaline Containing Copper and Tin Chloride. Trans. Inst. Met. Finish. 2014, 92, 41−46. (21) Mjalli, F. S.; Vakili-Nezhaad, G.; Shahbaz, K.; AlNashef, I. M. Application of the Eötvos and Guggenheim Empirical Rules for Predicting the Density and Surface Tension of Ionic Liquids Analogues. Thermochim. Acta 2014, 575, 40−44. (22) Shahbaz, K.; Mjalli, F. S.; Hashim, M. A.; AlNashef, I. M. Prediction of Deep Eutectic Solvents Densities at Different Temperatures. Thermochim. Acta 2011, 515, 67−72. (23) D’Agostino, C.; Harris, R. C.; Abbott, A. P.; Gladden, L. F.; Mantle, M. D. Molecular Motion and Ion Diffusion in Choline Chloride Based Deep Eutectic Solvents Studied by 1h Pulsed Field Gradient Nmr Spectroscopy. Phys. Chem. Chem. Phys. 2011, 13, 21383−21391. (24) Mjalli, F. S.; Naser, J. Viscosity Model for Choline ChlorideBased Deep Eutectic Solvents. Asia-Pac. J. Chem. Eng. 2015, 10, 273− 281. (25) Mjalli, F. S.; Ahmed, O. U. Physical Properties and Intermolecular Interaction of Eutectic Solvents Binary Mixtures: Reline and Ethaline. Asia-Pac. J. Chem. Eng. 2016, 11, 549. (26) D’Agostino, C.; Gladden, L. F.; Mantle, M. D.; Abbott, A. P.; Ahmed, E. I.; Al-Murshedi, A. Y. M.; Harris, R. C. Molecular and Ionic Diffusion in Aqueous - Deep Eutectic Solvent Mixtures: Probing InterMolecular Interactions Using Pfg Nmr. Phys. Chem. Chem. Phys. 2015, 17, 15297−15304. (27) Abbott, A. P.; Capper, G.; McKenzie, K. J.; Ryder, K. S. Voltammetric and Impedance Studies of the Electropolishing of Type 316 Stainless Steel in a Choline Chloride Based Ionic Liquid. Electrochim. Acta 2006, 51, 4420−4425. (28) Popescu, A. M. J.; Constantin, V.; Olteanu, M.; Demidenko, O.; Yanushkevich, K. Obtaining and Structural Characterization of the Electrodeposited Metallic Copper from Ionic Liquids. Rev. Chim. 2011, 62, 626−632. (29) Perkins, S. L.; Painter, P.; Colina, C. M. Experimental and Computational Studies of Choline Chloride-Based Deep Eutectic Solvents. J. Chem. Eng. Data 2014, 59, 3652−3662. (30) 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. (31) Perkins, S. L.; Painter, P.; Colina, C. M. Molecular Dynamic Simulations and Vibrational Analysis of an Ionic Liquid Analogue. J. Phys. Chem. B 2013, 117, 10250−10260.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work had the financial support of Fundaçaõ para a Ciência e a Tecnologia (FCT/MEC) through national funds and was cofinanced by FEDER under partnership agreement PT2020 (LAQV@REQUIMTE: projects UID/QUI/50006/2013 and POCI/01/0145/FEDER/007265). E.S.C.F. and I.V.V. further acknowledge FCT for postdoctoral grants SFRH/BPD/90343/ 2012 and SFRH/BPD/97918/2013, respectively, cofinanced by the European Social Fund. I.V.V. and C.M.P. are grateful for funding from COST Action CM1206 − EXIL − Exchange on Ionic Liquids.

■

REFERENCES

(1) 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. (2) Leron, R. B.; Soriano, A. N.; Li, M.-H. Densities and Refractive Indices of the Deep Eutectic Solvents (Choline Chloride + Ethylene Glycol or Glycerol) and Their Aqueous Mixtures at the Temperature Ranging from 298.15 to 333.15 K. J. Taiwan Inst. Chem. Eng. 2012, 43, 551−557. (3) Ferreira, E. S. C.; Pereira, C. M.; Silva, A. F. Electrochemical Studies of Metallic Chromium Electrodeposition from a Cr(Iii) Bath. J. Electroanal. Chem. 2013, 707, 52−58. (4) Salomé, S.; Pereira, N. M.; Ferreira, E. S.; Pereira, C. M.; Silva, A. F. Tin Electrodeposition from Choline Chloride Based Solvent: Influence of the Hydrogen Bond Donors. J. Electroanal. Chem. 2013, 703, 80−87. (5) Abbott, A. P.; Harris, R. C.; Ryder, K. S.; D’Agostino, C.; Gladden, L. F.; Mantle, M. D. Glycerol Eutectics as Sustainable Solvent Systems. Green Chem. 2011, 13, 82−90. (6) Durand, E.; Lecomte, J.; Villeneuve, P. Deep Eutectic Solvents: Synthesis, Application, and Focus on Lipase-Catalyzed Reactions. Eur. J. Lipid Sci. Technol. 2013, 115, 379−385. (7) Koller, T. M.; Rausch, M. H.; Ramos, J.; Schulz, P. S.; Wasserscheid, P.; Economou, I. G.; Fröba, A. P. Thermophysical Properties of the Ionic Liquids [Emim][B(Cn)4] and [Hmim][B(Cn)4]. J. Phys. Chem. B 2013, 117, 8512−8523. (8) Wu, S.-H.; Caparanga, A. R.; Leron, R. B.; Li, M.-H. Vapor Pressure of Aqueous Choline Chloride-Based Deep Eutectic Solvents (Ethaline, Glyceline, Maline and Reline) at 30−70 °C. Thermochim. Acta 2012, 544, 1−5. (9) Yadav, A.; Pandey, S. Densities and Viscosities of (Choline Chloride + Urea) Deep Eutectic Solvent and Its Aqueous Mixtures in the Temperature Range 293.15 to 363.15 K. J. Chem. Eng. Data 2014, 59, 2221−2229. (10) Smith, E. L.; Abbott, A. P.; Ryder, K. S. Deep Eutectic Solvents (Dess) and Their Applications. Chem. Rev. 2014, 114, 11060−11082. (11) Shahbaz, K.; Mjalli, F. S.; Hashim, M. A.; AlNashef, I. M. Prediction of the Surface Tension of Deep Eutectic Solvents. Fluid Phase Equilib. 2012, 319, 48−54. (12) Leron, R. B.; Li, M.-H. Molar Heat Capacities of Choline Chloride-Based Deep Eutectic Solvents and Their Binary Mixtures with Water. Thermochim. Acta 2012, 530, 52−57. (13) Lin, C.-M.; Leron, R. B.; Caparanga, A. R.; Li, M.-H. Henry’s Constant of Carbon Dioxide-Aqueous Deep Eutectic Solvent (Choline Chloride/Ethylene Glycol, Choline Chloride/Glycerol, Choline M


Article

The Journal of Physical Chemistry B (32) Szefczyk, B.; D. S. Cordeiro, M. N. Physical Properties at the Base for the Development of an All-Atom Force Field for Ethylene Glycol. J. Phys. Chem. B 2011, 115, 3013−3019. (33) Jorgensen, W. L. Optimized Intermolecular Potential Functions for Liquid Alcohols. J. Phys. Chem. 1986, 90, 1276−1284. (34) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the Opls All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (35) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. J. Phys. Chem. B 2004, 108, 2038−2047. (36) Górny, K.; Dendzik, Z.; Pabiszczak, M.; Gburski, Z. Non-Debye Dipolar Relaxation of Ethylene Glycol Embedded in Zsm-5 Zeolite Host Matrix  Computer Simulation Study. J. Non-Cryst. Solids 2013, 364, 28−33. (37) Abraham, M.; van der Spoel, D.; Lindahl, E.; Hess, B. Gromacs User Manual version 5.0.4; 2014. (38) 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. (39) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (40) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (41) Hockney, R. W.; Goel, S. P.; Eastwood, J. W. Quiet HighResolution Computer Models of a Plasma. J. Comput. Phys. 1974, 14, 148−158. (42) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. Lincs: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (43) Lindahl, E.; Hess, B.; van der Spoel, D. Gromacs 3.0: A Package for Molecular Simulation and Trajectory Analysis. Molecular modeling annual 2001, 7, 306−317. (44) Bedrov, D.; Borodin, O. Thermodynamic, Dynamic, and Structural Properties of Ionic Liquids Comprised of 1-Butyl-3Methylimidazolium Cation and Nitrate, Azide, or Dicyanamide Anions. J. Phys. Chem. B 2010, 114, 12802−12810. (45) Del Pópolo, M. G.; Voth, G. A. On the Structure and Dynamics of Ionic Liquids. J. Phys. Chem. B 2004, 108, 1744−1752. (46) Ciccotti, G.; Jacucci, G.; McDonald, I. R. Thought-Experiments” by Molecular Dynamics. J. Stat. Phys. 1979, 21, 1−22. (47) The Cp2k Developers Group; Cp2k is Freely Available from http://www.Cp2k.org/, accessed November 2015. (48) Kinz-Thompson, C.; Conwell, E. Proton Transfer in Adenine− Thymine Radical Cation Embedded in B-Form DNA. J. Phys. Chem. Lett. 2010, 1, 1403−1407. (49) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (50) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (51) Mantz, Y. A.; Gervasio, F. L.; Laino, T.; Parrinello, M. Solvent Effects on Charge Spatial Extent in DNA and Implications for Transfer. Phys. Rev. Lett. 2007, 99, 058104. (52) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127, 114105. (53) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (54) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (55) Ciocirlan, O. I.; Iulian, O.; Croitoru, O. Effect of Temperature on the Physico-Chemical Properties of Three Inic Liquids Containing Choline Chloride. Rev. Chim (Bucharest) 2010, 61, 721−723.

(56) Voroshylova, I. V.; Chaban, V. V. Atomistic Force Field for Pyridinium-Based Ionic Liquids: Reliable Transport Properties. J. Phys. Chem. B 2014, 118, 10716−10724. (57) Chaban, V. V.; Voroshylova, I. V. Systematic Refinement of Canongia Lopes−Pádua Force Field for Pyrrolidinium-Based Ionic Liquids. J. Phys. Chem. B 2015, 119, 6242−6249. (58) 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. (59) Vergadou, N.; Androulaki, E.; Hill, J.-R.; Economou, I. G. Molecular Simulations of Imidazolium-Based Tricyanomethanide Ionic Liquids Using an Optimized Classical Force Field. Phys. Chem. Chem. Phys. 2016, 18, 6850−6860. (60) Ishizuka, R.; Matubayasi, N. Self-Consistent Determination of Atomic Charges of Ionic Liquid through a Combination of Molecular Dynamics Simulation and Density Functional Theory. J. Chem. Theory Comput. 2016, 12, 804−811. (61) Fileti, E. E.; Chaban, V. V. The Scaled-Charge Additive Force Field for Amino Acid Based Ionic Liquids. Chem. Phys. Lett. 2014, 616−617, 205−211. (62) Choi, E.; Yethiraj, A. Entropic Mechanism for the Lower Critical Solution Temperature of Poly(Ethylene Oxide) in a Room Temperature Ionic Liquid. ACS Macro Lett. 2015, 4, 799−803. (63) Son, C. Y.; McDaniel, J. G.; Schmidt, J. R.; Cui, Q.; Yethiraj, A. First-Principles United Atom Force Field for the Ionic Liquid Bmim +Bf4−: An Alternative to Charge Scaling. J. Phys. Chem. B 2016, 120, 3560−3568. (64) Schröder, B.; Coutinho, J. A. P. Predicting Enthalpies of Vaporization of Aprotic Ionic Liquids with Cosmo-Rs. Fluid Phase Equilib. 2014, 370, 24−33. (65) Verevkin, S. P.; Zaitsau, D. H.; Emel’yanenko, V. N.; Yermalayeu, A. V.; Schick, C.; Liu, H.; Maginn, E. J.; Bulut, S.; Krossing, I.; Kalb, R. Making Sense of Enthalpy of Vaporization Trends for Ionic Liquids: New Experimental and Simulation Data Show a Simple Linear Relationship and Help Reconcile Previous Data. J. Phys. Chem. B 2013, 117, 6473−6486. (66) Borodin, O.; Smith, G. D.; Kim, H. Viscosity of a Room Temperature Ionic Liquid: Predictions from Nonequilibrium and Equilibrium Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 4771−4774. (67) Zhang, Y.; Xue, L.; Khabaz, F.; Doerfler, R.; Quitevis, E. L.; Khare, R.; Maginn, E. J. Molecular Topology and Local Dynamics Govern the Viscosity of Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2015, 119, 14934−14944. (68) Saiz, L.; Padró, J. A.; Guàrdia, E. Structure of Liquid Ethylene Glycol: A Molecular Dynamics Simulation Study with Different Force Fields. J. Chem. Phys. 2001, 114, 3187−3199. (69) Murli, C.; Lu, N.; Dong, Z.; Song, Y. Hydrogen Bonds and Conformations in Ethylene Glycol under Pressure. J. Phys. Chem. B 2012, 116, 12574−12580. (70) Luzar, A.; Chandler, D. Effect of Environment on Hydrogen Bond Dynamics in Liquid Water. Phys. Rev. Lett. 1996, 76, 928−931.

N


Improved Force Field Model for the Deep Eutectic Solvent Ethaline

Recommend Documents