OPLS Force Field for Choline Chloride-Based ... - ACS Publications

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

OPLS Force Field for Choline Chloride-Based Deep Eutectic Solvents Brian Doherty, and Orlando Acevedo J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b06647 • Publication Date (Web): 20 Aug 2018 Downloaded from http://pubs.acs.org on August 21, 2018

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

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

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

The Journal of Physical Chemistry

Revised version 8/16/18

OPLS Force Field for Choline Chloride-Based Deep Eutectic Solvents Brian Doherty and Orlando Acevedo* Department of Chemistry, University of Miami, Coral Gables, Florida 33146 E-mail: [email protected] Submitted July 11, 2018 Abstract: Deep eutectic solvents (DES) are a class of solvents frequently composed of choline chloride and a neutral hydrogen bond donor (HBD) at ratios of 1:1, 1:2, or 1:3, respectively. As cost-effective and eco-friendly solvents, DESs have gained considerable popularity in multiple fields, including materials, separations, and nanotechnology. In the present work, a comprehensive set of transferable parameters have been fine-tuned to accurately reproduce bulkphase physical properties and local intermolecular interactions for 8 different choline chloridebased DESs. This nonpolarizable force field, OPLS-DES, gave near quantitative agreement at multiple temperatures for experimental densities, viscosities, heat capacities, and surface tensions yielding overall mean absolute errors (MAEs) of ca. 1.1%, 1.6%, 5.5%, and 1.5%, respectively. Local interactions and solvent structuring between the ions and HBDs, including urea, glycerol, phenol, ethylene glycol, levulinic acid, oxalic acid, and malonic acid, were accurately reproduced when compared to radial distribution functions and coordination numbers derived from experimental liquid-phase neutron diffraction data and from first-principles molecular dynamics simulations. The reproduction of transport properties presented a considerable challenge and behaved more like a supercooled liquid near room temperature; higher temperature simulations, e.g., 400–500 K, or an alternative polarizable force field is recommended when computing self-diffusion coefficients.

1 ACS Paragon Plus Environment

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

Page 2 of 38

Revised version 8/16/18 Introduction An emerging family of ionic fluids called deep eutectic solvents (DES) have been shown to be an environmentally safer and a relatively inexpensive alternative to ionic liquids.1-7 DESs possess many favorable physical properties, such as, low vapor pressures, non-flammability, wide liquid range, and excellent biodegradability, but compared to ionic liquids, provide a more facile synthesis, a simpler purification process, and a 100% atom utilization rate.8 Accordingly, multiple fields have benefitted from their use, including biodiesel production,9 gas separation,10,11 materials,12,13 extractions,14 nanotechnology,15 and electrochemistry.16 As a medium for chemical reactions, DESs also possess numerous advantages,17 but are only beginning to emerge as potential alternatives to organic solvents.18-21 DESs are made by mixing together a solid quaternary ammonium salt and a neutral hydrogen bond donor (HBD), which like ionic liquids, often have melting points below 100 °C. While many combinations are available, DESs are frequently composed of the cost-effective choline chloride ion pair and an organic hydrogen bond donor such as urea or glycerol at ratios of 1:1, 1:2, or 1:3 (Figure 1 and Table 1).22,23

Figure 1. Deep eutectic solvent hydrogen bond donors and acceptors.


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


Revised version 8/16/18 Table 1. Deep eutectic solvents composed of choline chloride (ChCl) and a hydrogen bond donor (HBD) at specific ratios. Abbreviation CCEtg CCGly CCLev CCMal CCOx CCPhe(1:2) CCPhe(1:3) CCU

HBD Ethylene Glycol Glycerol Levulinic Acid Malonic Acid Oxalic Acid Phenol Phenol Urea

ChCl:HBD 1:2 1:2 1:2 1:1 1:1 1:2 1:3 1:2

Name Ethaline Glyceline Maline Oxaline Reline

Computational efforts employing quantum mechanics,24-26 QM/MD,27 ab initio molecular dynamics (AIMD),28-30 and first principles molecular dynamics (FPMD),31 have challenged earlier suggestions that simple delocalization of negative charge via strong hydrogen-bond interactions from the HBD prevents the high-melting point solid organic compounds from crystallizing.13 However, accurate predictions of thermodynamic and transport properties for bulk phase DESs requires hundreds to thousands of explicit solvent molecules which places the aforementioned theoretical methods out of reach in terms of computational cost. Molecular dynamics (MD) is well-suited to describe such properties as evidenced by the success of previous ionic liquid simulations.32-34 Force field development for DESs began with work by Colina and coworkers through the modification of generalized amber force field (GAFF) parameters.35,36 Despite a lack of experimental data available for parameter refinement, good agreement was found with measured densities and heat capacities; however, reproducing transport properties, i.e., self-diffusion coefficients, proved to be more difficult. For example, CCGly was overestimated by 17–27% at 330 K and underestimated by 14–20% at 298 K.35,36 3 ACS Paragon Plus Environment


Page 4 of 38

Revised version 8/16/18 Recent work by Mainberger et al. further modified the GAFF force field for DESs, but still overestimated the self-diffusivities by 7–16%, whereas their modification to the Merck molecular force field (MMFF) found the self-diffusion coefficients to be underestimated.37 The inconsistencies reported in these nonpolarizable force fields and in other MD studies38,39 highlight the need for a systematic developmental effort dedicated to DESs. For example, Ferreira et al. reported40 a 10% improvement in simulated transport properties compared to Colina’s work36 when specifically focusing upon a single DES, i.e., CCEtg, and by utilizing a combination of OPLS-AA parameters originally intended for ionic liquid simulations.41,42 In this work, nonpolarizable OPLS-AA parameters have been created for a transferrable choline chloride ion pair and multiple organic hydrogen bond donors (Figure 1) which have been fine-tuned to accurately reproduce bulk-phase physical properties and local intermolecular interactions. Overall, the novel partial charges and Lennard-Jones (LJ) terms yielded close agreement with experimental DES data that included densities, heat capacities, surface tensions, and viscosities. However, reproduction of self-diffusion coefficients at multiple temperatures still presented a considerable challenge. Radial distribution functions (RDFs) and solvent coordination numbers were computed and compared favorably to those derived from experiment and first-principles molecular dynamics. The new force field reported here for deep eutectic solvents is called OPLS-DES.


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


Revised version 8/16/18 Computational Methods OPLS-AA Force Field. The OPLS-AA force field calculates the total energy as the sum of individual intra- and intermolecular energies. Specifically, the harmonic bond stretching and angle bending terms, a Fourier series for dihedral angles, and the intermolecular energies are computed using Coulomb and 12-6 Lennard-Jones terms (Eq. 1-4). The adjustable parameters are the force constants k, the ro and qo equilibrium bond and angle values, Fourier coefficients V, partial atomic charges, q, and Lennard-Jones radii and well-depths, s and e.

;

𝐸"#$%& = ∑+ 𝑘",+ (𝑟+ − 𝑟/,+ )1

(1)

𝐸2$345& = ∑+ 𝑘6,+ (𝜃+ − 𝜃/,+ )1

(2)

;

;

;

𝐸8#9&+#$ = ∑+ :1 𝑉;,+ (1 + 𝑐𝑜𝑠 𝜙+ ) + 1 𝑉1,+ (1 − cos 2𝜙+ ) + 1 𝑉H,+ (1 + cos 3𝜙+ ) + 1 𝑉J,+ (1 − cos 4𝜙+ )L

𝐸$#$"#$% = ∑+ ∑S[+ M

NONP 5 Q 9OP

V

;1

OP + 4𝜀+S TU 9 W OP

V

X

OP − U 9 W YZ OP

(3)

(4)

Geometric combining rules were used for the Lennard-Jones coefficients, i.e., sij = (siisjj)1/2 and eij = (eiiejj)1/2. Nonbonded interactions were calculated intermolecularly and for intramolecular atom pairs separated by three or more bonds. To apply the same parameters for both intra- and intermolecular interactions, the 1,4-intramolecular interactions were scaled by 0.5. Molecular Dynamics. Molecular dynamic (MD) simulations were carried out using the GROMACS 5.0.7 software package.43 Cubic boxes containing 500 choline chloride pairs and the proper ratio of HBDs from Table 1 for each DES were constructed using the program Packmol.44 Periodic boundary conditions were employed with the short-range electrostatic cutoff set to 16



Page 6 of 38

Revised version 8/16/18 Å. Long-range electrostatics were treated with Particle-Mesh Ewald summations. Minimization of each box was performed using a steepest-descent algorithm for at least 5000 steps. Temperature ranges were held constant using a velocity rescaling with a stochastic term (vrescale)45 and pressure was maintained at 1.0 bar using the Berendsen coupling method during an isothermal-isobaric ensemble (NPT) for 10 ns of equilibration. Equations of motion were integrated with the leap-frog algorithm using a time step of 1 fs. The LINCS method constrained all covalently bound hydrogen atoms. Production runs were carried out for an additional 50 ns. Combined distribution functions (CDF) and spatial distribution functions (SDF) were calculated utilizing the TRAVIS program.46

Results and Discussion OPLS-DES Parameters. Inter- and intramolecular parameters for the choline chloride ion pair were developed starting with potentials reported in our OPLS-2009IL ionic liquid force field.33,41 Molecular charge scaling of DESs by noninteger values, i.e., ±0.8e and ±0.9e, has been previously investigated by Colina as a means to improve predictions of bulk-phase and transport properties.35,36 A ±0.8e scaling was reported to yield the best experimental agreement for CCU, whereas ±0.9e provided better results for other DESs, e.g., CCGly. Additional DES parameterization efforts40 and multiple ionic liquid studies32,33 have reported ±0.8e to provide the best overall agreement in both bulk-phase properties and local intermolecular interactions over a broad range of solvents. Accordingly, parameters for the OPLS-DES were developed specifically towards a ±0.8e charge scaling for the cations, anions, and HBDs. DES simulations have been shown to be particularly sensitive to the methodological choice of computing partial atomic


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


Revised version 8/16/18 charges.47 As such, the choline chloride nonbonded parameters, i.e., charges and Lennard-Jones terms, were instead adjusted to match radial distribution functions derived by Hammond et al. from fitting liquid-phase CCU neutron diffraction data with an empirical potential structure refinement (ND/EPSR) model.48 Choline geometries and torsional parameters were adjusted to fit conformational energy minima from LMP2/cc-pVTZ(-f) calculations. Potentials for the organic hydrogen donors were taken directly from the OPLS-AA force field49 and the nonbonded terms were fine-tuned towards the choline chloride pair to reproduce experimental and/or ab initio MD derived RDFs and bulk physical properties. As such, the new parameters for the organic HBD molecules should only be used for DES calculations and not in pure liquid simulations. All the OPLS-DES parameters are provided in the Supporting Information Tables S1-S6. Density. The OPLS-DES computed densities and their corresponding literature values are given in Table 2 for a range of temperatures. A weighted experimental average was derived by determining the weight of each data point from the inverse of its uncertainty. Reported in the Supporting Information are all of the individual experimental values, uncertainties, and references. Overall mean absolute errors (MAE) of 1.1–1.3% for 8 unique DESs over 5 different temperatures emphasize a close agreement with the measured densities. However, it should be noted that the parameters gave higher errors of 3.9–4.2% for the CCU solvent. This deviation is the result of a compromise as urea required additional scaling of the LJ terms to improve the agreement with other bulk-phase properties, including the heat capacity and surface tension, and local intermolecular interactions (discussed later).



Page 8 of 38

Revised version 8/16/18 Table 2. Calculated and Experimental Liquid Densities (g/cm3) Across Multiple Temperatures for Choline Chloride-Based Deep Eutectic Solvents. Temperature (K)

298.15

303.15 313.15 323.15 333.15 CCEtg OPLS-DES 1.132 1.129 1.121 1.115 1.107 a Exptl. 1.117 1.114 1.108 1.103 1.097 % Error 1.3 1.3 1.2 1.0 0.9 CCGly OPLS-DES 1.193 1.192 1.186 1.179 1.173 a Exptl. 1.192 1.189 1.183 1.177 1.172 % Error 0.1 0.2 0.2 0.1 0.1 CCLev OPLS-DES 1.135 1.132 1.125 1.118 1.110 a Exptl. 1.138 1.134 1.128 1.121 1.115 % Error 0.2 0.2 0.3 0.3 0.4 CCMal OPLS-DES 1.232 1.231 1.230 1.228 1.223 a Exptl. 1.231 1.224 1.221 1.215 1.209 % Error 0.1 0.5 0.7 1.1 1.1 CCOx OPLS-DES 1.237 1.233 1.231 1.227 1.223 a Exptl. 1.259 1.255 1.248 1.241 1.234 % Error 1.8 1.8 1.3 1.1 0.9 CCPhe(1:2) OPLS-DES 1.088 1.084 1.076 1.068 1.061 a Exptl. 1.097 1.093 1.087 % Error 0.8 0.8 1.0 CCPhe(1:3) OPLS-DES 1.086 1.082 1.074 1.066 1.058 a Exptl. 1.092 1.089 1.083 % Error 0.5 0.6 0.8 CCU OPLS-DES 1.148 1.147 1.141 1.135 1.129 a Exptl. 1.198 1.195 1.189 1.183 1.178 % Error 4.0 3.9 4.0 4.1 4.2 MAE (%) 1.1 1.2 1.2 1.3 1.3 a Weighted experimental average; all the individual values, uncertainties, and references are given in the Supporting Information. 8 ACS Paragon Plus Environment

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


Revised version 8/16/18 Viscosity. To treat the highly viscous nature of the DESs, the nonequilibrium periodic perturbation method was used to calculate the shear viscosity.50,51 Briefly, an external force, a, is applied that will create a velocity field, u, for the solvent (eq. 5). Three-dimensional periodic cells are utilized during the simulation with an external force chosen in the x direction, such that ax is a function of z only, and ay and az are equal to zero. The result is that uy and uz become zero as well. ]𝒖

𝑝 ]8 + 𝑝(𝒖 ⋅ 𝛁)𝐮 = p𝐚 − ∇p + η∇1 𝒖

(5)

A lack of pressure gradient in the x-direction allows ux to be simplified to equation 6. The velocity field is calculated during a simulation and is defined by a velocity profile, V. Consequently, the viscosity, h, of the system can be obtained from equation 7, where k = 2p/lz, and lz is the height of the box. 𝑝

]fg (h) ]8

= 𝑝𝑎j (𝑧) + 𝜂

] Q fg (h) ]h Q

m o

𝜂 = n pQ

(6)

(7)

An acceleration amplitude, L, of the external force ax(z) is chosen that provides a smooth velocity profile with small local shear rates (eq. 8). The L must be small enough so that the perturbation does not disturb the system equilibrium, but large enough to ensure proper sampling. 𝑎j (𝑧) = Λcos (kz)

(8)

Each system had 6 to 12 acceleration amplitudes varying from 0.04–0.24 nm/ps2. Each amplitude point was simulated by extending the production run by 10 ns. Plotting the viscosities versus the 9 ACS Paragon Plus Environment


Page 10 of 38

Revised version 8/16/18 amplitudes then allows for the undisturbed system where L = 0 to be extrapolated, and the intercept is taken as the calculated viscosity. The OPLS-DES gave excellent agreement with experiment for 14 different simulations at 3 different temperatures with an overall MAE of 1.6%. All individual percent errors were in the low single digit range. The largest error of 5.4% belonged to CCPhe(1:2) at 303.15 K; however, the accuracy is considerably improved to 1.3% at the lower temperature of 298.15 K. Table 3. Calculated and Experimental Viscosities (cP) at Various Temperatures for Choline Chloride-Based Deep Eutectic Solvents. DES CCEtg CCGly CCLev CCPhe(1:2) CCPhe(1:3) CCU CCEtg CCGly CCLev CCPhe(1:2) CCPhe(1:3) CCU

OPLS-DES 298.15 K 39.5 258.8 220.8 89.1 44.4 753.1 303.15 K 35.0 246.8 164.0 64.7 36.5 520.5 348.15 K 94.9 205.2

Exptl.a

% error

39.7 259.0 226.8 90.3 44.6 749.9

0.6 0.1 2.6 1.3 0.4 0.4

35.0 238.9 164.5 68.4 35.2 514.8

0.0 3.3 0.3 5.4 3.7 1.1

CCMal 96.7 1.9 CCOx 208.3 1.5 MAE (%) 1.6 a Weighted experimental average; all the individual values, uncertainties, and references are given in the Supporting Information.


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


Revised version 8/16/18 Heat Capacity. Heat capacities at constant pressure were computed by using classical statistical mechanics (Cpclass) corrected to account for neglected quantum mechanical vibrational Nv

contributions (𝛿𝐶n ) via the two-phase thermodynamics density of states (DoS) method.52

𝐶ow42&& = x

Nv

𝛿𝐶n

〈𝜕𝐻1 〉 ~ 𝑘" 𝑇 1 •

(9)

†

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

(10)

The Cpclass equation examines enthalpy and temperature fluctuations throughout the simulation (eq. 9), where H is the enthalpy, T is the temperature, and kb is the Boltzmann’s constant. Quantum mechanical effects are then accounted for by considering the distribution of normal modes as a function of the normal-mode frequency n, which is calculated53 using the Fourier transform of the mass-weighted atomic velocity autocorrelation functions.52,54,55 The distribution is normalized to the total number of degrees of freedom present in the system. Each normal mode is estimated to be a quantum mechanical harmonic oscillator with frequency n. The quantum mechanical correction to the classical heat capacity can then be approximated by eq. 10, where W is a weighting function. The final expression for the quantum-corrected heat capacity is given by eq. 11, which has provided accurate heat capacities in prior ionic liquid studies.32,33,56 Nv

𝐶ow#99 = 𝐶ow42&& + 𝛿𝐶n

(11)

The heat capacities, Cpcorr, were computed using the last 40 ns of a 50 ns trajectory. As the simulations utilized covalent hydrogen bond constraints, and each harmonic bond potential contributes kb to the heat capacity, a term Nckb was added to Cpclass where Nc is the number of constraints. The g_dos utility program in GROMACS was used to calculate the normal mode



Page 12 of 38

Revised version 8/16/18 distributions from the autocorrelation functions.53 The simulations were extended by 100 ps for the DoS calculations and employed the NVT ensemble with fully-flexible bonds. Atomic velocities were stored every 4 ns for the velocity autocorrelation function with a basic time step of 0.5 fs. All Cpclass and dCVqm values computed for each solvent are provided in the Supporting Information Table S7. Table 4 provides the computed heat capacities at 353.15 K for three DESs, i.e., CCEtg, CCGly, and CCU. The simulations provided reasonably accurate results with an overall MAE of 5.5%. As a point of comparison, the accurate ionic liquid force field OPLSVSIL gave an overall MAE of 9.5% for the Cp’s of 11 different ionic liquids32 and Perkins et al. reported 7% overestimations for the DESs CCEtg and CCGly.36

Table 4. Calculated (Quantum Corrected) and Experimental Heats Capacities (J/mol.K) at 353.15 K for Choline Chloride-Based Deep Eutectic Solvents. DES CCEtg CCGly CCU MAE (%) a Ref. 57.

Exptl.a 205.6 ± 0.2 254.3 ± 0.4 190.8 ± 0.8

OPLS-DES 215.8 240.1 201.9

% error 4.9 5.6 5.8 5.5

Surface Tension. Surface tensions were computed using DES boxes in which the z-axis was elongated by a factor of 3 to generate two liquid-vapor interfaces. Simulations in the NVT ensemble were carried out for 10 ns using the v-rescale algorithm and the pressure tensor, 𝑃$ , of the system was stored at every step. The surface tension, g, was computed using the directional components of the pressure tensor, and Px, Py, and Pz, and the length of the box in the z direction, 12 ACS Paragon Plus Environment

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


Revised version 8/16/18 Lz (eq. 12). Previous studies on ionic liquids reported overestimated surface tension predictions with strikingly large standard deviations when performed near room temperature.58-60 The high viscosities and slow dynamics for the DES systems also require elevated temperatures for proper equilibrium convergence at the interface. As such, surface tension simulations were performed at 425 K for 10 ns. 1 1 𝛾 = 𝐿h Š𝑃hh − ‹𝑃jj + 𝑃ŒŒ •Ž 2 2 (𝑇w − 𝑇)

(12) ;;⁄•

Y 𝛾 (𝑇) = 𝛾95• T ‹𝑇w − 𝑇95• •

(13)

As experimental values for DESs at elevated temperatures are unavailable, in the present work the surface tensions at 425 K have been estimated using the Othmer equation (eq. 13).61 Equation 13 requires the critical temperature, Tc, of the liquid, i.e., where the surface tension becomes zero, and a reference surface tension value, gref, at another temperature, Tref. Shahbaz et al. examined 9 DESs in a temperature range of 298–328 K and found the Othmer equation to accurately reproduce experimental surface tensions giving an average error of 2.57%.61 In this work, gref and Tref were determined at 298 K by using a weighted average of multiple reported experimental surface tensions. Critical temperatures were derived using a modified LydersonJoback-Reid method,62 which has been reported to predict highly accurate critical properties for 200 ionic liquids.63 Table 5 gives the OPLS-DES computed surface tensions for 4 DESs and their experimentally estimated values at 425 K. The simulations gave an overall MAE of 1.5%. These highly accurate results significantly surpass previous ionic liquid MD simulations that gave MAEs of 15.2% and 7.3% for 7 different solvents using the OPLS-2009IL and OPLS-VSIL force fields, respectively.32,33 13 ACS Paragon Plus Environment


Page 14 of 38

Revised version 8/16/18 Table 5. Calculated and Experimental Surface Tensions (mN/m) at 425 K for Choline ChlorideBased Deep Eutectic Solvents. DES CCEtg CCGly CCMal CCU

Exptl.a 35.4 44.1 52.3 38.7

OPLS-DES 35.9 43.2 51.3 38.9

MAE (%)

% error 1.4 2.0 1.9 0.5 1.5

a

Experimental estimate at 425 K using the Othmer equation and reference data at 298 K.

Self-Diffusion Coefficients. Accurate predictions of self-diffusivity are of particular importance due to widespread interest in transport properties for the DESs.8,64,65 Self-diffusion coefficients, Ds, for choline and the organic hydrogen bond donors of the DESs were computed by applying the Einstein relation (eq. 14) and a mean-square displacement (MSD) of the center of mass position for each ion/molecule, ri.66 To sample in the proper diffusive regime for extracting diffusivities, the MSD for each simulation trajectory was used to calculate the betaparameter (b) as previously discussed by Del Pópolo and Voth to determine the location of the diffusive regime in ionic liquids.67 Equation 15 shows how to calculate b as a function of time, t, where Dr(t)2 is the average MSD of the system. The b parameter was also used by Colina and coworkers in the calculation of self-diffusivities for DESs.35,36 At low temperatures and short time scales b < 1 is often observed as the system dynamics are in the subdiffusive regime, e.g., motions are similar to supercooled liquids.68 The appropriate region for computing self-diffusion coefficients is at b = 1, where the calculated diffusion coefficient is in the diffusive regime.


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


Revised version 8/16/18 ;

%

𝐷& = X lim %8 〈∑› ⃗+ (𝑡) − 𝑟⃗+ (0)]1 〉 +œ;[𝑟

(14)

8→†

𝛽 (𝑡 ) =

𝑑 log;/ 〈Δ𝑟(𝑡)1 〉 𝑑 log;/ 𝑡

(15)

A plot of b versus time found that DES simulations near room temperature, e.g., 298–328 K, lie exclusively within the subdiffusive regime (b < 1) over the entire trajectory and possess large fluctuations of the parameter (Supporting Information Figure S1). Increasing the temperature to over 400 K yielded trajectories that were more consistently within the diffusive regime (b = 1) (Supporting Information Figures S2-S4). This dependence is a consequence of a reduced viscosity at higher temperatures and was also found to be the case in computing selfdiffusion coefficients for ionic liquids.33,69 Table 6 provides the OPLS-DES computed selfdiffusion coefficients for choline (D+) and the HBD (DHBD) at the temperatures that encompassed the proper diffusive regime, i.e., 400.15–500.15 K. The trajectory ranges utilized for analysis at each temperature are given in the Supporting Information Table S8.

Table 6. Calculated Self-Diffusivity Coefficients (10-11 m2 s-1) Between 400.15–500.15 K for Choline Chloride Deep Eutectic Solvents.

T (K) 400.15 420.15 440.15 460.15 480.15 500.15

CCEtg D DHBD +

15.78 25.53 37.08 51.06 69.01 87.83

30.99 47.37 67.51 93.99 123.73 160.61

CCGly D DHBD +

3.67 5.60 9.33 13.69 22.29 27.66

8.61 13.86 21.10 30.61 41.72 50.93


CCU +

D

DHBD

11.72 20.03 31.07 44.90 61.58 80.25

24.28 39.68 61.25 86.79 114.40 145.90


Page 16 of 38

Revised version 8/16/18 Experimental self-diffusivity values for DESs, i.e., CCU, CCGly, and CCEtg, were reported by D’Agostino et al. within a temperature range of 298–328 K.64 Their pulse field gradient nuclear magnetic resonance (PFG-NMR) study of diffusion in DESs found that an Arrhenius relationship between the diffusion coefficient and temperature was uniformly followed.64 Examining the ability of the OPLS-DES parameters to reproduce experimental selfdiffusion coefficients at lower temperatures proved to be technically challenging. Comparison of simulations at 298–328 K to experiment would be invalid as they sampled exclusively within the subdiffusive regime. Alternatively, a linear fit of our calculated Ds values at 400–500 K was extrapolated down to temperatures of 298–328 K by plotting ln(D) versus 1000/T.64 However, the extrapolation introduced considerable error despite an R2 > 0.995 in all cases; reproducing the calculated Ds values between 400–500 K using the linear fit equations gave overall MAE errors of ~4% with individual errors as high as 8.2% (Supporting Information Table S9). In addition, comparison of the extrapolated self-diffusivities at 298–328 K to experiment yielded wildly different errors that were dependent upon the chosen temperature. For example, the CCU DES had errors of 31.4%, 0.0%, and 23.2% for D+ (choline cation) at temperatures of 298.15 K, 308.15 K, and 323.15 K, respectively (Supporting Information Table S10). This error pattern is not unique to our parameters as Ferreira et al. also found a parabola-like temperature dependence in their nonpolarizable force field MD study of CCEtg.40 They reported self-diffusion coefficients errors in CCEtg of 20.6%, 3.4%, and 13.3% for D+ at 298.15 K, 313.15 K, and 323.15 K, respectively.40 Colina and coworkers also computed a temperature dependence during their parameterization efforts of 3 different DESs; for example, CCGly was overestimated by 17–27% at 330 K and underestimated by 14–20% at 298 K.35,36 These results indicate that adequate sampling for self-diffusivity calculations requires, at minimum, higher temperatures,


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


Revised version 8/16/18 i.e., >400 K. However, given the reported inconsistencies from multiple groups, nonpolarizable force field do not appear to be well-suited for reproducing DES transport properties. Higher-level models, e.g., polarizable force fields, are likely required for improved accuracy. CCU solvent structure. Previous studies featuring computationally-intensive methods, such as first-principles molecular dynamics (FPDM),31 ab initio molecular dynamics (AIMD),29 and quantum mechanical molecular dynamics (QM/MD),27 gave generally accurate radial distribution function peak positions, heights, and shapes for the interactions between the CCU solvent ions and HBDs when compared to results from liquid-phase CCU neutron diffraction data with an empirical potential structure refinement (ND/EPSR) model.48 Some discrepancies were noted; for example, QM/MD indicated a preferential structuring of urea moieties through the bulk liquid, where proximal hydrogens interacted more strongly with chloride anions than distal hydrogen, which is the reverse trend of that observed experimentally.27 In addition, AIMD simulations of CCU had some issues quantitatively reproducing RDF peak heights. For example, the interaction between the hydroxyl hydrogen atom of choline and chloride is reported to have an approximated peak height of 6.5 at 2.1 Å from ND/EPSR, whereas AIMD predicted a value of ~9 at 2.10 Å.29 This highlights the large technical challenge of quantitatively reproducing RDFs for CCU, even with decidedly accurate theoretical methods. The OPLS-DES derived liquid structure of CCU was investigated by comparing RDFs and coordination numbers between the center-of-mass for the different ion and molecule combinations. RDFs provide the probability of finding an atom at a certain distance from a reference atom compared to the statistical average. Table 7 compares the OPLS-DES to results from Hammond et al.’s ND/EPSR48 and Fetisov et al’s FPMD31 studies. Agreement between the maximum and minimum of the first peaks (rmax and rmin) and the average coordination numbers 17 ACS Paragon Plus Environment


Page 18 of 38

Revised version 8/16/18 (Ncoord) compared to ND/EPSR and FPMD was excellent with the values primarily falling within statistical uncertainties. Supporting Information Figure S5 provides the center-of-mass RDF peak positions, heights, and shapes. The center-of-mass mean coordination number also predicted that each chloride is coordinated by 2 urea molecules (given the 2:1 ratio of urea:ChCl), which agrees with experimental observations.48

Table 7. Average Coordination Number (Ncoord) and Positions (Å) of the First Maximum and Minimum in Center-of-Mass RDFs between Choline Cation (Ch), Chloride Anion (Cl), and Urea in the CCU Deep Eutectic Solvent. OPLS-DES (303 K) rmax rmin Ncoord

ND/ESPR (303 K)a rmax rmin Ncoord

FPMD (333 K)b rmax rmin Ncoord

center

shell

urea Ch

Cl Cl

4.3 4.1

5.4 6.4

1.90 3.49

4.0 4.2

5.5 6.7

2.08 ± 1.01 4.35 ± 1.30

4.1 4.2

5.3 6.5

1.9 ± 0.4 3.1 ± 0.6

Ch

urea

4.7

7.2

8.76

5.4

6.9

5.91 ± 2.84

5.1

7.1

8.6 ± 0.7

Ch

Ch

6.5

8.2

5.41

6.3

8.0

6.74 ± 2.16

-

-

-

urea

urea

4.8

6.6

6.00

4.3

6.1

6.77 ± 3.05

4.7

6.3

4.9 ± 0.5

a

Ref. 48. bRef. 31.

Site-site RDFs were also analyzed between the hydrogen bonding atoms of choline, chloride, and urea (Figure 2 and Table 8). The strongest interaction was between the hydrogen atom on the hydroxyl group of choline (HY) and chloride with a peak height of 7.2 centered at 2.20 Å (Figure 2a), which compares fairly well to the approximated peak height of 6.5 at 2.1 Å reported from ND/EPSR. This intense HY–Cl interaction gave a computed mean coordination number of 0.63 (exptl. Ncoord = 0.66 ± 0.50). The remaining hydrogen atoms on choline, i.e., HA, 18 ACS Paragon Plus Environment

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


Revised version 8/16/18 HS, and HW in Figure 2a, had significantly lower intensity peaks centered around 2.90 Å and the OPLS-DES accurately reproduced the ND/EPSR peak shapes and positions (see Hammond et al.’s fig. 3 in their publication48). No distinct structuring preference was found between chloride and any one particular choline hydrogen atom (HA, HS, or HW) as the mean coordination numbers and RDF structuring were all comparable (Table 8). The second most dominate interaction was between urea and chloride, as the anion primarily interacted with the HC and HT atoms bound to the nitrogen atoms. Consistent with experimental findings, the HC–Cl interaction was preferred over HT–Cl with mean coordination numbers of 2.35 and 1.64, respectively (Figure 2b and Table 8). RDF analysis for HC–Cl found intensities and distances of 4.03 at 2.42 Å and 1.57 at 3.74 Å for the first and second peaks, respectively. In contrast, the HT–Cl had weaker intensities at longer distances, i.e., 3.05 at 2.44 Å and 1.35 at 4.08 Å for the first and second peaks, respectively. The most prominent peak computed between urea and chloride was Cl–N and was centered at 3.38 Å with an intensity of 4.5, which is nearly indistinguishable from the ND/EPSR values of approximately 3.4 Å and 4.5, respectively. The comparatively weaker interactions between choline-urea, choline-choline, and ureaurea were also monitored (Figure 2c-e) and the structural ordering was very similar to the published ND/EPSR RDF plots.48 The coordination of urea around choline was found to be particularly weak, as low RDF intensities of g(r) < 2, were computed for all interactions (Figure 2c). The most significant interaction was that of the choline hydroxy HY atom and a nitrogen atom from urea with a computed coordination number of 2.96, which implies the participation of roughly one and half urea molecules. Another one and half urea molecules approximately coordinate with the remaining non-hydroxyl hydrogen atoms on choline. The RDF interactions between two choline molecules were also weak and possessed long interacting distances (Figure 19 ACS Paragon Plus Environment


Page 20 of 38

Revised version 8/16/18 3d). The ND/EPSR study found the most distinct interaction to be between the N atoms, with a peak centered at approximately 6 Å,48 which was well-reproduced by the OPLS-DES parameters. Finally, the RDF analysis between two urea molecules found an expected hydrogen bond occurring between the urea hydrogen atoms (HC and HT) and the oxygen of urea (Figure 2e). The hydrogen atoms in the “trans” orientation to the oxygen atom (HT) gave a higher intensity peak of 1.20 compared to 0.75 for the “cis” orientation (HC) and both were centered at 2.3 Å. The OPLS-DES computed Ncoord values were 0.86 and 0.49 urea molecules for the O–HT and O– HC, respectively; however, the ND/ESPR RDFs give a much larger Ncoord difference between HT and HC, i.e., 2.33 ± 0.95 and 0.62 ± 0.84. Curiously, visual inspection of the ND/ESPR RDF plots for O–HT and O–HC finds both curves to be nearly identical (see Hammond et al.’s fig. 3(e) in their publication48), which does not appear to correlate with a four-fold integration difference. Zahn also reported AIMD simulations of CCU at 375 K that gave Ncoord values, 0.78 and 0.93 for O–HT and O–HC, that were closer in magnitude to OPLS-DES predictions, but reversed the H-bond preference.29 Lastly, the O–N interaction gave a sharp peak centered around 3.26 Å with a Ncoord of 2.19 (exptl. = 2.25 ± 1.50), which is indicative of an ordered head-to-tail urea-urea preference.


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



Figure 2. Radial distribution functions (site-site) between the different species present in the CCU deep eutectic solvent.



Page 22 of 38

Revised version 8/16/18 Table 8. Calculated Coordination Numbers (Ncoord) for the CCU Deep Eutectic Solvent by Integrating Each Radial Distribution Function Peak Over a Radius Range Spanning Each Minimum (Rcoord given in Å). RDF Ch-Cl

Cl-urea

Ch-urea

Ch-Ch

urea-urea

a

Atom Aa HY HA HS HW Cl Cl Cl Cl HS HA HW HY HY NA NA NA NA NA O O O N N

Atom Ba Cl Cl Cl Cl C N HC HT O O O O N NA CS CW OY HY N HC HT HC HT

Rcoord 1.7-3.1 2.3-3.9 2.3-4.0 2.3-3.9 3.5-5.4 2.8-4.6 1.9-3.2 1.9-3.2 2.2-3.5 2.1-3.4 2.1-3.6 1.5-2.6 1.8-4.2 5.1-7.7 4.1-7.8 3.9-6.3 3.3-6.0 2.9-5.8 2.7-4.2 1.7-2.9 1.7-3.1 2.0-4.9 2.0-4.8

See Figure 2 for atom names of CCU. bRef. 48.


Ncoord 0.63 0.68 0.64 0.69 3.80 4.92 2.35 1.64 0.65 0.72 0.74 0.18 2.96 4.37 4.47 1.87 1.82 2.87 2.19 0.49 0.86 4.25 4.11

ND/ESPRb 0.66 ± 0.50 0.70 ± 0.66 0.73 ± 0.67 0.70 ± 0.66 3.76 ± 2.27 4.10 ± 2.63 1.73 ± 1.58 1.25 ± 1.15 0.41 ± 0.60 0.48 ± 0.66 0.45 ± 0.64 0.16 ± 0.38 2.08 ± 1.68 3.48 ± 1.72 5.34 ± 2.03 2.05 ± 1.34 2.31 ± 1.11 2.73 ± 1.31 2.25 ± 1.50 0.62 ± 0.84 2.33 ± 0.95 5.62 ± 2.65 5.77 ± 2.50

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


Revised version 8/16/18 In order to monitor the hydrogen bond network within the CCU DES, combined distribution functions (CDFs) were used to correlate radial and angular distributions to elucidate the average orientation for each molecule within the solvent (Figure 3). A quantitative description of hydrogen bonding in DESs was fashioned by Perkins et al. as a set of separate criteria that recommended a distance < 3.5 Å between the donor (X) and acceptor (Y) heavy atoms and an X–H…Y angle of >150°.35 However, an alternative hydrogen bond definition proposed by Wernet et al. for liquid water coupled the distances and angles of X and Y to define an ellipse.70 Accordingly, Fetisov et al. utilized a similar DES H-bond analysis based on elliptical dependencies that gave a more accurate description relative to a rectangular boundary.31 Their analysis of CCU found the >150° criteria to be too restrictive and the < 3.5 Å distance to vary largely based on the interacting atoms, i.e., N, O, and Cl. Instead, Figure 3 highlights the dominant hydrogen bonding region between (A) chloride and choline (Cl–OY), (B) choline and urea (OY–O), and (C) chloride and urea (Cl–N). The OPLS-DES parameters accurately reproduced the FPMD H-bond analysis reported in fig. 6 of Fetisov et al.’s publication.31 For example, the majority of chloride anions were bound to the HC and HT atoms of urea (Figure 3C), whereas choline seldomly formed substantial H-bond interactions with the urea hydrogen atoms. In addition, the O from urea was a very weak hydrogen bond acceptor with the hydroxy group of choline (Figure 3B) and instead favored forming hydrogen bonds with a N from a second urea molecule (Table 8).



Page 24 of 38


Figure 3. Combined distribution functions for the CCU deep eutectic solvent from the OPLSDES simulations. Illustrations are given of the plotted angle a versus the distance d for (A) the Cl- anion and the oxygen atom (OY) of choline, (B) the OY atom of choline and O atom of urea, and (C) the Cl- anion and a N atom from urea. The CDF plots are given with a relative intensity color.

CCEt and CCOx solvent structure. The DESs CCU, CCEtg, and CCOx were studied by Zahn et al. using ab initio molecular dynamics.28 A follow-up AMID study by Zahn gave detailed insight into CCU liquid structuring and hydrogen bonding dynamics.29 Given their high accuracy, the AIMD calculations provide an excellent benchmark to gauge the local structuring predictions from the OPLS-DES parameters for CCEtg and CCOx. The OPLS-DES did a very good job of reproducing the AIMD RDF peak shapes and positions for CCEtg and CCOx 24 ACS Paragon Plus Environment

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


Revised version 8/16/18 (compare Figures 4 and 5 here to figs. 2a and 2b of Zahn et al.’s publication,28 respectively). The agreement in terms of distances between the hydrogen bond donor atoms and chloride anions were particularly excellent (Table 9); a correlation between hydrogen bond strength and donoracceptor distance is well-known.71 Integration up to a distance of 4.50 Å of the first solvation shell peak between chloride and the choline oxygen atom (Cl–OY) for CCOx and CCEtg revealed an average coordination number (Ncoord) of 1.73 and 1.31 anions, respectively, surrounding the choline nitrogen, which compared favorably to 1.47 and 1.22 reported from the AIMD simulations.28

Table 9. Comparison of Radial Distribution Function Maximum Peak Heights, g(r), for the First Solvation Shell and Their Distances (rmax in Å) for the CCEtg and CCOx Deep Eutectic Solvents at 375 K. OPLS-DES rmax g(r) CCEtg CCOx

4.28 4.26

CCEtg CCOx

3.08 3.08

CCEtg CCOx

3.06 2.94

Cl-NA 4.6 4.8 Cl-OY 2.5 2.6 Cl-Ob 4.6 6.6

AIMDa rmax

g(r)

4.20 4.20

3.5 3.6

3.05 3.05

3.9 4.6

3.10 2.95

4.0 6.5

a

Ref. 28. bO = OG in CCEtg or OH in CCOx. See Figures 4 and 5 for atom names.



Page 26 of 38


Figure 4. Radial distribution functions (site-site) between the different species present in the CCEtg deep eutectic solvent.


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



Figure 5. Radial distribution functions (site-site) between the different species present in the CCOx deep eutectic solvent.

CCGly Structure. The liquid structure of CCGly has been previously studied by Korotkevich et al. using AIMD to monitor SO2 solvation into the DES.30 The strongest hydrogen bond predicted by AIMD was between chloride and the OH groups on glycol, i.e., HM and HO in Figure 6. The RDFs indicated that Cl- had the most dominant interaction with the middle hydroxyl group on glycerol (HM–Cl) with a peak height of ~9, whereas interaction with the 27 ACS Paragon Plus Environment


Page 28 of 38

Revised version 8/16/18 peripheral OH groups (HO–Cl) gave a peak height at ~7.5. While OPLS-DES predicted the correct trend of HM > HO when interacting with chloride, the peak heights were significantly lower at 4.6 and 3.9, respectively. Similar RDF peak heights and distances were reported by Mainberger et al. with a nonpolarizable force field utilizing GAFF charges scaled by ±0.75e for CCGly.37 However, the AIMD computed RDF between the chloride and the hydroxyl group of the choline cation (HY–Cl) gave a peak height of ~7.5, which was well-reproduced with a value of 7.4 from the OPLS-DES (Figure 6).

Figure 6. Radial distribution functions (site-site) between the different species present in the CCGly deep eutectic solvent.


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


Revised version 8/16/18 The interaction between the hydrogen atom on the hydroxyl group on choline (HY) and the terminal (OH) and center (OM) oxygen atoms on glycerol were also computed in Figure 6. The RDF analysis found both HY–OH and HY–OM possessed particularly small peak heights of approximately 0.5. In contrast, the AIMD simulations had larger peak heights of around 2.8-3.2 for the same interactions. Mainberger et al.’s nonpolarizable force field also gave very small peak heights similar to that of OPLS-DES when scaling the GAFF charges by ±0.75e.37 The unscaled GAFF charges increased the peak heights of HY–OH and HY–OM to approximately 1.5, but were still considerably below that of AIMD. Table 10 gives the average coordination number (Ncoord) between the CCGly ions and molecules and the OPLS-DES provided good agreement with the reported AIMD simulations.

Table 10. Calculated Coordination Numbers (Ncoord) for the CCGly Deep Eutectic Solvent by Integrating Each Radial Distribution Function Peak Over a Radius Range Spanning Each Minimum (Rcoord given in Å). RDF Ch-Cl Gly-Cl Ch-Gly Gly-Gly a

Atom Aa HY HO HM HY HY HO HM

Atom Ba Cl Cl Cl OH OM OH OM

Rcoord 1.8-3.0 1.9-3.0 1.9-3.0 1.7-2.5 1.7-2.5 1.7-3.0 1.7-2.6

See Figure 6 for atom names of CCGly. bRef. 30.


Ncoord 0.56 0.31 0.37 0.17 0.07 0.88 0.14

AIMDb 0.50 1.11 0.68 0.26 0.16 0.53 0.10


Page 30 of 38

Revised version 8/16/18 Spatial Distribution Functions. Spatial distribution functions (SDFs) can supply additional insight as to the preferential positioning of the ions surrounding the HBDs of the solvents. The TRAVIS program46 was used to construct and analyze the SDFs for all DESs listed in Table 1. The SDF isosurface values for the systems are provided in the Supporting Information Table S11. Figure 7 shows the most probable location of finding the chloride anion (given in orange) and choline cation (shown in green) relative to each DES HBD. For example, the CCU solvent was found to coordinate two urea molecules per chloride from the center-ofmass average coordination number (Table 7). Site-site RDF analysis determined that the anion strongly interacted with the HC and HT atoms bound to the urea nitrogen atoms (Table 8) and, accordingly, the SDF analysis in Figure 7A confirmed the chloride population density clustered around the hydrogen atoms. Choline, as expected, preferentially occupied the position above the carbonyl oxygen. In a similar fashion, the other HBD hydrogen atoms also preferentially interacted with the chloride anions, whereas the oxygen atoms of the HBDs hydrogen bonded to the hydroxyl group of choline. This ability to form specific hydrogen bonding intermolecular order, while maintaining favorable electrostatic interactions that impedes crystallization has been discussed as the geometric driving force for the deep eutectic solvent formation.48


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



Figure 7. Spatial distribution functions for (A) urea, (B) ethylene glycol, (C) levulinic acid, (D) malonic acid, (E) oxalic acid, (F) glycerol, and (G) phenol in choline chloride deep eutectic solvents. Orange denotes the position of the chloride anion and green depicts the choline cation. 31 ACS Paragon Plus Environment


Page 32 of 38

Revised version 8/16/18 Conclusions New force field potentials, called OPLS-DES, for choline chloride-based deep eutectic solvents have been developed and validated on 8 unique DESs. The physical properties, including densities, heat capacities, surface tensions, and viscosities, gave highly accurate reproduction of the experimentally measured values, i.e., typical MAEs of approximately 1-5%. Transport properties were more difficult to reproduce as they were highly sensitive to simulation temperature and gave errors comparable to other nonpolarizable force fields. More computationally intensive methodology, such as polarizable force fields, are likely required for accurate reproduction of self-diffusion coefficients. Nonbonded parameters, i.e., charges and Lennard-Jones terms, were adjusted to match radial distribution functions reported by Hammond et al. that fit liquid-phase DES neutron diffraction data with an empirical potential structure refinement (ND/EPSR) model.48 The solvent structure and local intermolecular interactions were monitored by using radial distribution functions, combined distribution functions, and spatial distribution functions. Accordingly, computed peak heights and intensities for each DES plot and solvent coordination numbers were comparable to those of experiment and QM-based methods, e.g., FPMD and AIMD simulations. All of the OPLS-DES parameters can be downloaded from http://github.com/orlandoacevedo/DES.

Acknowledgement. Gratitude is expressed to the National Science Foundation (CHE-1562205) and the Center for Computational Science at the University of Miami for support of this research.


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


Revised version 8/16/18 Supporting Information Available: OPLS-AA bonded and nonbonded parameters for choline chloride and all HBDs discussed; heat capacity adjustment values; b parameter plots for selfdiffusivity calculations; SDF isovalues; cubic lengths for equilibrated boxes; solvent experimental values, uncertainties, measurement methodologies, and references. The Supporting Information is available free of charge on the ACS Publications website at http://pubs.acs.org.

References (1) Clarke, C. J.; Tu, W.-C.; Levers, O.; Bröhl, A.; Hallett, J. P. Green and sustainable solvents in chemical processes. Chem. Rev. 2018, 118, 747-800. (2) Smith, E. L.; Abbott, A. P.; Ryder, K. S. Deep eutectic solvents (DESs) and their applications. Chem. Rev. 2014, 114, 11060–11082. (3) Khandelwal, S.; Tailor, Y. K.; Kumar, M. Deep eutectic solvents (DESs) as eco-friendly and sustainable solvent/catalyst systems in organic transformations. J. Mol. Liq. 2016, 215, 345-386. (4) Tang, B.; Row, K. H. Recent developments in deep eutectic solvents in chemical sciences. Monatsh. Chem. 2013, 144, 1427-1454. (5) Francisco, M.; van den Bruinhorst, A.; Kroon, M. C. Low-transition-temperature mixtures (lttms): A new generation of designer solvents. Angew. Chem. Int. Ed. 2013, 52, 3074-3085. (6) 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. (7) Abbott, A. P.; Boothby, D.; Capper, G.; Davies, D. L.; Rasheed, R. K. Deep eutectic solvents formed between choline chloride and carboxylic acids: Versatile alternatives to ionic liquids. J. Am. Chem. Soc. 2004, 126, 9142-9147. (8) Zhang, Q.; Vigier, K. D. O.; Royer, S.; Jerome, F. Deep eutectic solvents: Syntheses, properties and applications. Chem. Soc. Rev. 2012, 41, 7108-7146. (9) Zhao, H.; Zhang, C.; Crittle, T. D. Choline-based deep eutectic solvents for enzymatic preparation of biodiesel from soybean oil. J. Mol. Catal. B Enzym. 2013, 85-86, 243-247. (10) Leron, R. B.; Li, M.-H. Solubility of carbon dioxide in a choline chloride–ethylene glycol based deep eutectic solvent. Thermochim. Acta 2013, 551, 14-19. (11) Yang, D.; Hou, M.; Ning, H.; Zhang, J.; Ma, J.; Yang, G.; Han, B. Efficient SO2 absorption by renewable choline chloride-glycerol deep eutectic solvents. Green Chem. 2013, 15, 22612265. (12) Wagle, D. V.; Zhao, H.; Baker, G. A. Deep eutectic solvents: Sustainable media for nanoscale and functional materials. Acc. Chem. Res. 2014, 47, 2299–2308. (13) Carriazo, D.; Serrano, M. C.; Gutiérrez, M. C.; Ferrer, M. L.; del Monte, F. Deep-eutectic solvents playing multiple roles in the synthesis of polymers and related materials. Chem. Soc. Rev. 2012, 41, 4996-5014.



Page 34 of 38

Revised version 8/16/18 (14) Gu, T.; Zhang, M.; Tan, T.; Chen, J.; Li, Z.; Zhang, Q.; Qiu, H. Deep eutectic solvents as novel extraction media for phenolic compounds from model oil. Chem. Commun. 2014, 50, 11749-11752. (15) Abo-Hamad, A.; Hayyan, M.; AlSaadi, M. A.; Hashim, M. A. Potential applications of deep eutectic solvents in nanotechnology. Chem. Eng. J. 2015, 273, 551-567. (16) Sebastián, P.; Gómez, E.; Climent, V.; Feliu, J. M. Copper underpotential deposition at gold surfaces in contact with a deep eutectic solvent: New insights. Electrochem. Commun. 2017, 78, 51-55. (17) Alonso, D. A.; Baeza, A.; Chinchilla, R.; Guillena, G.; Pastor, I. M.; Ramón, D. J. Deep eutectic solvents: The organic reaction medium of the century. Eur. J. Org. Chem. 2016, 2016, 612-632. (18) Liu, P.; Hao, J.-W.; Mo, L.-P.; Zhang, Z.-H. Recent advances in the application of deep eutectic solvents as sustainable media as well as catalysts in organic reactions. RSC Adv. 2015, 5, 48675-48704. (19) Ruß, C.; König, B. Low melting mixtures in organic synthesis – an alternative to ionic liquids? Green Chem. 2012, 14, 2969-2982. (20) Ghumro, S. A.; Alharthy, R. D.; al-Rashida, M.; Ahmed, S.; Malik, M. I.; Hameed, A. Nalkylated 1,4-diazabicyclo[2.2.2]octane–polyethylene glycol melt as deep eutectic solvent for the synthesis of fisher indoles and 1H-tetrazoles. ACS Omega 2017, 2, 2891-2900. (21) Marset, X.; Guillena, G.; Ramón, D. J. Deep eutectic solvents as reaction media for the palladium-catalysed C−S bond formation: Scope and mechanistic studies. Chem. Eur. J. 2017, 23, 10522-10526. (22) López-Porfiri, P.; Brennecke, J. F.; Gonzalez-Miquel, M. Excess molar enthalpies of deep eutectic solvents (DESs) composed of quaternary ammonium salts and glycerol or ethylene glycol. J. Chem. Eng. Data 2016, 61, 4245-4251. (23) Abbott, A. P.; Capper, G.; Davies, D. L.; Rasheed, R. K.; Tambyrajah, V. Novel solvent properties of choline chloride/urea mixtures. Chem. Commun. 2003, 70-71. (24) García, G.; Atilhan, M.; Aparicio, S. An approach for the rationalization of melting temperature for deep eutectic solvents from DFT. Chem. Phys. Lett. 2015, 634, 151-155. (25) Ashworth, C. R.; Matthews, R. P.; Welton, T.; Hunt, P. A. Doubly ionic hydrogen bond interactions within the choline chloride-urea deep eutectic solvent. Phys. Chem. Chem. Phys. 2016, 18, 18145-18160. (26) Rimsza, J. M.; Corrales, L. R. Adsorption complexes of copper and copper oxide in the deep eutectic solvent 2:1 urea–choline chloride. Comput. Theor. Chem. 2012, 987, 57-61. (27) Stefanovic, R.; Ludwig, M.; Webber, G. B.; Atkin, R.; Page, A. J. Nanostructure, hydrogen bonding and rheology in choline chloride deep eutectic solvents as a function of the hydrogen bond donor. Phys. Chem. Chem. Phys. 2017, 19, 3297-3306. (28) Zahn, S.; Kirchner, B.; Mollenhauer, D. Charge spreading in deep eutectic solvents. ChemPhysChem 2016, 17, 3354-3358. (29) Zahn, S. Deep eutectic solvents: Similia similibus solvuntur? Phys. Chem. Chem. Phys. 2017, 19, 4041-4047. (30) Korotkevich, A.; Firaha, D. S.; Padua, A. A. H.; Kirchner, B. Ab initio molecular dynamics simulations of SO2 solvation in choline chloride/glycerol deep eutectic solvent. Fluid Phase Equilibria 2017, 448, 59-68.


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


Revised version 8/16/18 (31) Fetisov, E. O.; Harwood, D. B.; Kuo, I. F. W.; Warrag, S. E. E.; Kroon, M. C.; Peters, C. J.; Siepmann, J. I. First-principles molecular dynamics study of a deep eutectic solvent: Choline chloride/urea and its mixture with water. J. Phys. Chem. B 2018, 122, 1245-1254. (32) Doherty, B.; Zhong, X.; Acevedo, O. Virtual site OPLS force field for imidazolium-based ionic liquids. J. Phys. Chem. B 2018, 122, 2962-2974. (33) Doherty, B.; Zhong, X.; Gathiaka, S.; Li, B.; Acevedo, O. Revisiting OPLS force field parameters for ionic liquid simulations. J. Chem. Theory Comput. 2017, 13, 6131-6145. (34) Salanne, M. Simulations of room temperature ionic liquids: From polarizable to coarsegrained force fields. Phys. Chem. Chem. Phys. 2015, 17, 14270-14279. (35) 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. (36) 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. (37) Mainberger, S.; Kindlein, M.; Bezold, F.; Elts, E.; Minceva, M.; Briesen, H. Deep eutectic solvent formation: A structural view using molecular dynamics simulations with classical force fields. Mol. Phys. 2017, 115, 1309-1321. (38) Sun, H.; Li, Y.; Wu, X.; Li, G. Theoretical study on the structures and properties of mixtures of urea and choline chloride. J. Mol. Model. 2013, 19, 2433–2441. (39) Ullah, R.; Atilhan, M.; Anaya, B.; Khraisheh, M.; Garcıa, G.; ElKhattat, A.; Tariq, M.; Aparicio, S. A detailed study of cholinium chloride and levulinic acid deep eutectic solvent system for CO2 capture via experimental and molecular simulation approaches. Phys. Chem. Chem. Phys. 2015, 17, 20941-20960. (40) Ferreira, E. S. C.; Voroshylova, I. V.; Pereira, C. M.; Cordeiro, M. N. D. S. Improved force field model for the deep eutectic solvent ethaline: Reliable physicochemical properties. J. Phys. Chem. B 2016, 120, 10124−10137. (41) 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. (42) Canongia Lopes, J. N.; Deschamps, J.; Pádua, A. H. Modeling ionic liquids using a systematic all-atom force field. J. Phys. Chem. B 2004, 108, 2038-2047. (43) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19-25. (44) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 21572164. (45) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (46) Brehm, M.; Kirchner, B. Travis - a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories. J. Chem. Inf. Model. 2011, 51, 2007-2023. (47) García, G.; Atilhan, M.; Aparicio, S. The impact of charges in force field parameterization for molecular dynamics simulations of deep eutectic solvents. J. Mol. Liq. 2015, 221, 506-514. (48) Hammond, O. S.; Bowron, D. T.; Edler, K. J. Liquid structure of the choline chloride-urea deep eutectic solvent (reline) from neutron diffraction and atomistic modelling. Green Chem. 2016, 18, 2736-2744.



Page 36 of 38

Revised version 8/16/18 (49) 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. (50) Zhao, L.; Wang, X.; Wang, L.; Sun, H. Prediction of shear viscosities using periodic perturbation method and opls force field. Fluid Phase Equilib. 2007, 260, 212-217. (51) Hess, B. Determining the shear viscosity of model liquids from molecular dynamics simulations. J. Chem. Phys. 2002, 116, 209-217. (52) Lin, S.-T.; Blanco, M.; Goddard III, W. A. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. J. Chem. Phys. 2003, 119, 11792-11805. (53) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force field benchmark of organic liquids: Density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory Comput. 2012, 8, 61-74. (54) Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water. J. Chem. Phys. 1983, 79, 23752389. (55) Pascal, T. A.; Lin, S.-T.; Goddard III, W. A. Thermodynamics of liquids: Standard molar entropies and heat capacities of common solvents from 2pt molecular dynamics. Phys. Chem. Chem. Phys. 2011, 13, 169-181. (56) Chen, M.; Pendrill, R.; Widmalm, G.; Brady, J. W.; Wohlert, J. Molecular dynamics simulations of the ionic liquid 1-n-butyl-3-methylimidazolium chloride and its binary mixtures with ethanol. J. Chem. Theory Comput. 2014, 10, 4465-4479. (57) 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. (58) Pensado, A. S.; Malfreyt, P.; Pádua, A. A. H. Molecular dynamics simulations of the liquid surface of the ionic liquid 1-hexyl-3-methylimidazolium bis(trifluoromethanesulfonyl)amide: Structure and surface tension. J. Phys. Chem. B 2009, 113, 14708–14718. (59) Lynden-Bell, R. M.; Del Pópolo, M. G. Simulation of the surface structure of butylmethylimidazolium ionic liquids. Phys. Chem. Chem. Phys. 2006, 8, 949-954. (60) Bhargava, B. L.; Balasubramanian, S. Layering at an ionic liquid−vapor interface: A molecular dynamics simulation study of [BMIM][PF6]. J. Am. Chem. Soc. 2006, 128, 10073– 10078. (61) Shahbaz, K.; Mjalli, F. S.; Hashim, M. A.; AlNashef, I. M. Prediction of the surface tension of deep eutectic solvents. Fluid Phase Equilibria 2012, 319, 48-54. (62) Joback, K. K.; Reid, R. Estimation of pure component properties from group contribution. Chem. Eng. Commun. 1987, 57, 233-247. (63) Valderrama, J. O.; Sanga, W. W.; Lazzús, J. A. Critical properties, normal boiling temperature, and acentric factor of another 200 ionic liquids. Ind. Eng. Chem. Res. 2008, 47, 1318-1330. (64) 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. (65) Shen, Y.; He, X.; Hung, F. R. Structural and dynamical properties of a deep eutectic solvent confined inside a slit pore. J. Phys. Chem. C 2015, 119, 24489−24500.


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


Revised version 8/16/18 (66) Liu, X.; Schnell, S. K.; Simon, J.-M.; Krüger, P.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Diffusion coefficients from molecular dynamics simulations in binary and ternary mixtures. Int. J. Therophys. 2013, 34, 1169-1196. (67) Del Pópolo, M. G.; Voth, G. A. On the structure and dynamics of ionic liquids. J. Phys. Chem. B 2004, 108, 1744-1752. (68) Qian, J.; Hentschke, R.; Heuer, A. Dynamic heterogeneities of translational and rotational motion of a molecular glass former from computer simulations. J. Chem. Phys. 1999, 110, 45144522. (69) Kowsari, M. H.; Alavi, S.; Ashrafizaadeh, M.; Najafi, B. Molecular dynamics simulation of imidazolium-based ionic liquids. I. Dynamics and diffusion coefficient. J. Chem. Phys. 2008, 129, 224508. (70) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Näslund, L. Å.; Hirsch, T. K.; Ojamäe, L.; Glatzel, P.; et al. The structure of the first coordination shell in liquid water. Science 2004, 304, 995-999. (71) Wendler, K.; Thar, J.; Zahn, S.; Kirchner, B. Estimating the hydrogen bond energy. J. Phys. Chem. A 2010, 114, 9529–9536.



Page 38 of 38

Revised version 8/16/18 TOC graphic


OPLS Force Field for Choline Chloride-Based ... - ACS Publications

Recommend Documents