Article pubs.acs.org/jced
Cite This: J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Extension of TEAM Force-Field Database to Ionic Liquids Zheng Gong and Huai Sun* School of Chemistry and Chemical Engineering, Materials Genome Initiative Center, and Key Laboratory of Scientific and Engineering Computing of Ministry of Education, Shanghai Jiao Tong University, Shanghai 200240, China
Downloaded via UNIV OF LOUISIANA AT LAFAYETTE on April 20, 2019 at 08:28:57 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: The TEAM force-field database is extended to cover roomtemperature ionic liquids (IL). The training set for parametrization includes 11 cations in five classes, imidazolium, pyridinium, pyrrolidinium, piperidinium, and tetraalkylammonium, and 12 anions, chloride, perchlorate, nitrate, tetrafluoroborate, hexafluorophosphate, thiocyanate, dicyanamide, tricyanomethanide, tetracyanoborate, bis(trifluoromethylsulfonyl)imide, triflate, and trifluoroacetate. The force field is compatible with other force fields of the TEAM database. The parameters are primarily derived from quantum-mechanical data, with a few nonbonded parameters optimized using experimental data of density and viscosity. The extrapolation of shear viscosity calculated using the nonequilibrium periodic perturbation (PP) method is examined by analyzing the Newtonian region in the calculated shear viscosities. With the identified extrapolation scheme, the PP method is about five times more efficient than the multiconfiguration-ensemble Green−Kubo (GK) method to reach similar convergence, and it is used for iterated calculations of viscosities as a part of the parametrization. The force field is validated using 46 ILs; the unsigned deviations from the experimental data are 1.4%, 12.8%, and 3.7% for density, viscosity, and isobaric heat capacity, respectively. The force field is used to analyze the interactions in terms of activation energies of viscosity and liquid structures in terms of distribution functions and free energy maps for common ILs. The ring−ring stacking is found to be the most stable local configuration for ILs of aromatic cation and chloride anion, which has a profound impact on the high viscosities of the ILs. field, all-atom force field (AAFF) or coarse-grained force field (CGFF).20,21 At different levels of approximation, the temporal and spatial scales of the simulation are obviously different. More importantly, the physical phenomena to be described are different. Using ab initio MD, the electron motions are taken into consideration; therefore, the polarization and chargetransfer phenomena, which have a profound effect on the local and long-range structural relaxation of ILs, are described on the first principles basis. The downside is that the simulation is too expensive to be carried out with sufficient sampling in the configuration space. Using a polarizable force field, the polarization and charge transfer can be represented by parametrized models, significantly advancing the temporal and spatial scales of the simulation.22−24 However, the computational cost is still a concern, and the parametrization must be carried out carefully, which is common for all force fields. Using AAFF, the polarization and charge-transfer effect are described by adjusting the force-field parameters.25−34 The fluidity of ILs is underestimated with the atomic partial charges derived from isolated ions.31,35 A common practice is to scale
I. INTRODUCTION Because of their unique properties, room-temperature ionic liquids (ILs) have been studied for a wide range of applications in which they are used as solvent, catalyst, adsorbent, electrolyte, and so on.1−5 ILs can be made from a number of cation classes, such as imidazolium, pyridinium, pyrrolidinium, piperidinium, ammonium, phosphonium, triazonium, guanidinium, and sulfonium. Different cation classes exhibit different properties. For example, pyridinium has the advantage of low cost and high thermal stability;6 pyrrolidinium, piperidinium, and tetraalkylammonium exhibit high electrochemical stability;7−9 phosphonium and ammonium share similar structure but different viscosity;10,11 and triazolium has a low melting point.12 The selection of anions is more diversified, providing variations in thermal stability,13 hydrolytic stability,14 toxicity,15 and adsorption capacity.16 The abundant combinations of cations and anions, coupled to chemical modifications, gives rise to the great potential of manipulating the physical properties of ILs. Molecular dynamics (MD) simulation is a valuable tool for investigating the structures and properties of ILs.17−19 With an increased level of approximation, MD simulations have been carried out using forces calculated by quantum mechanics (QM) density functional theory that may be generally referred to as ab initio MD, or based on the parametrized potential energy function (force field) in terms of the polarizable force © XXXX American Chemical Society
Special Issue: FOMMS Received: January 15, 2019 Accepted: April 8, 2019
A
DOI: 10.1021/acs.jced.9b00050 J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Journal of Chemical & Engineering Data
Article
the partial charges by a factor ranging from 0.6 to 0.9.26,34−37 The idea is that the charge transfer and dielectric screening occurring in the condensed phase weakens the electrostatic and the dispersion interactions.35,38−41 It has been shown that AAFF with the scaled atomic partial charges can be applied to predict energetic, structural, and kinetic properties under a wide range of conditions and ionic types.34,36,42 Using CGFF, the long-range structural correlations such as phase separation and self-assembly can be predicted, which is critical for applications where ILs are used alone or together with polymers as electrolytes in energy-storage devices,43−47 yet coarse graining removes the structural details of the molecule, and the transport properties cannot be predicted directly. Therefore, some kind of scaling scheme must be applied.48−50 In the multiscale simulation schemes summarized above, the AAFF plays a critical role in bridging the microscopic and mesoscopic levels. AAFF is simple enough that the simulations can be carried out in a system containing up to one million atoms, and the simulation time can be extended up to microseconds with the current computer power. It is sufficiently sophisticated to predict most thermodynamic and transport properties using MD simulations without any further calibration, except for the force-field parametrization. AAFF is also a stepping stone for further development. An ab initio CGFF can be developed using the bottom-up approach from data generated by AAFF simulations,51 analogous to an AAFF developed from data of QM calculations. In this work, we extend the TEAM force field41,52,53 to include common ILs. Despite several developments in this category,25,27,29 the necessity of this work is mainly for practical considerations: (a) The coverage must be sufficiently large and (b) the force field must be compatible with other force fields describing relevant but different materials. Both considerations come from the engineering perspective of real applications. We chose five classes of cations, imidazolium (Im), pyridinium (Py), pyrrolidinium (Pyrr), piperidinium (Pi), and tetraalkylammonium (N), and 12 anions, chloride (Cl), perchlorate (ClO4), nitrate (NO3), tetrafluoroborate (BF4), hexafluorophosphate (PF6), thiocyanate (NCS), dicyanamide (NC2N), tricyanomethanide (NC3C), tetracyanoborate (NC4B), bis(trifluoromethylsulfonyl)imide (TF2N), triflate (TFSO3), and trifluoroacetate (TFCO2), as the training set for parametrization (Figure 1). These cations and anions are the most common components of ILs; therefore, the method used in this work can be applied to development of other ILs, and the parameters obtained in this work can be transferred to ions containing the same chemical moieties. Because the ILs may be used together with other inorganic or organic (polymeric) materials, particularly in the applications of separation and energy storage, it is important to have compatible force fields of different molecules. The TEAM force-field database is designed for this purpose; it has multiple tables, each table is an independent force field covering a specific type of molecules, and all force fields share common functional forms and are parametrized using the same methodology. For a specific application, parameters taken from different force fields in the same functional forms are combined to a single force field by using parameter combination rules. In force-field applications, the interactions between different molecules are modeled as physical interactions only; therefore, the common combination rules of the Coulomb and Lennard-Jones (LJ) functions are used.
Figure 1. Structures, denotations, and atom types of ions considered in this work. The subscripts n, m, o, and p denote the numbers of carbon atom in alkyl tails. All hydrogen atoms connected to carbon share the same atom type of h_1 and are omitted in the figure for clarity.
The parametrization method used in this work is the same as that in other developments for TEAM force fields. Essentially, almost all parameters are derived by fitting the QM data (bottom-up); only a few LJ parameters are further optimized using high-quality experimental data (top-down). We chose to use the equilibrium density and viscosity for the top-down optimization because those are the most available experimental data of ILs. To calculate viscosities accurately for highly viscous fluids is nontrivial; we focused on how to calculate the viscosities accurately and efficiently. By comparing the extrapolation of shear viscosities calculated using the periodic perturbation (PP) method with the multiconfigurationensemble Green−Kubo (GK) method, we found that the PP method yields the same results as the GK method with about 1/5 of the computational resources. The resulting force field is tested on 46 ILs; the experimental densities and viscosities are well reproduced. The accurate prediction of isobaric heat capacity further validates this force field. As a preliminary application, we analyze the temperature dependency of viscosity and liquid structures of several representative ILs. In the following sections, we will first present the methods of parametrization and simulation and then discuss the results in terms of molecular property, liquid density, viscosity, and liquid structure, followed by conclusions and outlooks.
II. METHODS The function forms of TEAM-AMBER52 are used, which is written as B
DOI: 10.1021/acs.jced.9b00050 J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Journal of Chemical & Engineering Data
Article
Table 1. Equivalent Atom Types in This Force Fielda AAT
NB
ATC
BINC
BOND
A/C
A/S
T/C
T/S
O/C
O/S
b_4c_2nb c_2nc c_2nn c_2ns c_3c_35an c_35an2 c_3a c_3an c_3o2 c_45h2 c_45np c_4f3 c_4h2 c_4h3 c_4np c_4nph3 cl0cl4 f_1 f_1b f_1p h_1 h_1n n_1 n_2n_3+a n_35+da n_3o n_4+ n_45+ o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2nb c_2nc c_2nn c_2ns c_3c_35an c_35an2 c_3a c_3an c_3o2 c_45h2 c_45np c_4f3 c_4h2 c_4h3 c_4np c_4nph3 cl0cl4 f_1 f_1b f_1p h_1 h_1n n_1 n_2n_3+a n_35+da n_3o n_4+ n_45+ o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2nb c_2nc c_2nn c_2ns c_3c_35an c_35an2 c_3a c_3an c_3o c_45 c_45np c_4f3 c_4h2 c_4h3 c_4np c_4np cl0cl4 f_1 f_1 f_1 h_1 h_1n n_1 n_2n_3+a n_35+da n_3o n_4+ n_45+ o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2nb c_2nc c_2nn c_2ns c_3c_35an c_35an2 c_3a c_3an c_3o c_45 c_45np c_4 c_4 c_4 c_4np c_4np cl0cl4 f_1 f_1 f_1 h_1 h_1n n_1 n_2n_3+a n_35+da n_3o n_4+ n_45+ o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2n c_2n c_2n c_2n c_3c_35an c_35an2 c_3a c_3a c_3o c_45 c_45 c_4 c_4 c_4 c_4 c_4 cl0cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2n_3+a n_35+da n_3o n_4+ n_45+ o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2 c_2 c_2 c_2 c_3c_35an c_35an2 c_3a c_3a c_3o c_45 c_45 c_4 c_4 c_4 c_4 c_4 cl0cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2n_3a n_35+da n_3o n_4+ n_45 o_1 o_1-d o_1-q o_1-t p_6s_1s_4o
b_4c_2 c_2 c_2 c_2 c_3c_35an c_35an2 c_3a c_3a c_3o c_45 c_45 c_4 c_4 c_4 c_4 c_4 cl0cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2n_3a n_35+da n_3o n_4+ n_45 o_1 o_1o_1-q o_1p_6s_1s_4o
b_4c_2 c_2 c_2 c_2 c_3c_35an c_35an2 c_3a c_3a c_3o c_45 c_45 c_4 c_4 c_4 c_4 c_4 cl0cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2n_3a n_35+da n_3o n_4+ n_45 o_1 o_1o_1-q o_1p_6s_1s_4o
b_4c_2 c_2 c_2 c_2 c_3 c_35an c_35an2 c_3a c_3a c_3 c_4 c_4 c_4 c_4 c_4 c_4 c_4 cl0 cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2 n_3a n_35+da n_3 n_4+ n_45 o_1 o_1 o_1-q o_1 p_6s_1 s_4
b_4c_2 c_2 c_2 c_2 c_3 c_35an c_35an2 c_3a c_3a c_3o c_4 c_4 c_4 c_4 c_4 c_4 c_4 cl0cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2n_3a n_35+da n_3o n_4+ n_45 o_1 o_1o_1-q o_1p_6s_1s_4o
b_4c_2 c_2 c_2 c_2 c_3 c_35an c_35an2 c_3a c_3a c_3 c_4 c_4 c_4 c_4 c_4 c_4 c_4 cl0 cl4 f_1 f_1 f_1 h_1 h_1 n_1 n_2 n_3a n_35+da n_3 n_4+ n_45 o_1 o_1 o_1-q o_1 p_6s_1 s_4
a
AAT, apparent atom type; NB, van der Waals; ATC, atom charge; BINC, bond charge; BOND, bond; A/C, angle center; A/S, angle side; T/C, torsion center; T/S, torsion side; O/C, improper torsion center; O/S, improper torsion side.
E=
∑
kb(b − b0)2 +
bonds
+
∑
represented by two parameters: atomic charge and bond increment
kθ(θ − θ0)2
angles
∑
∑
kϕ , n(1 + cos(nϕ − ϕ0, n))
qi = qi ,0 +
dihedrals n = 1,2,3
∑
+
+
∑ i,j
Cqiqj εr
j
kχ (1 − cos(2χ ))
out‐of‐planes
+
∑ δij
i
6 i σ y yz − 2jjj zzz zzzz kr{ { kk r {
∑ ϵjjjjjijjj σ yzzz i,j
(2)
where j represents all atoms bonded to atom i. The bond increment δij represents the charge separation between two valence-bonded atoms i and j.54 The atom types are assigned using the hierarchical atomtype definition scheme.52 As shown in Figure 1, 39 atom types are defined, in which 19 are newly introduced to the TEAM database. An equivalence atom-type table52 is used to minimize the number of adjustable parameters, as shown in Table 1. For example, four sp-hybridized carbon atoms in nitrile-substituted ions, c_2nb, c_2nc, c_2nn, and c_2ns, have different vdW and charge parameters but share the same bond parameters identified by an atom type c_2n, and the angle, torsion, and
12
(1)
On the right side of the equation, the first four are bonded terms, and the last two are nonbonded terms. The bonded terms are functions of the bond length, bond angle, dihedral angle, and out-of-plane angle. The nonbonded terms are a function of the atom−atom separation, with the Coulomb function for the electrostatic and the LJ function for the van der Waals (vdW) interactions. The atomic partial charge is C
DOI: 10.1021/acs.jced.9b00050 J. Chem. Eng. Data XXXX, XXX, XXX−XXX
Journal of Chemical & Engineering Data
Article
analyzed by calculating the spatial distribution function (SDF) and radial distribution function (RDF), and free-energy maps were constructed from the radial-angle combined distribution function (CDF). Simulations were performed at 343 K unless otherwise specified. The uncertainties were estimated as standard errors calculated with a block average of five independent blocks. The GROMACS 201667 package was used for MD simulations. Open Babel68 and PACKMOL69 were used to build the simulation boxes. Gaussian 0970 was used for the QM calculations. SDF and CDF were calculated by utilizing the TRAVIS71 program.
improper torsion parameters are identified with the even more general atom type c_2. Although this force field shares atom types with other force fields in the TEAM database, the parameters in different force fields may or may not be the same because each force field is developed independently. This strategy provides the flexibility to have the best coverage of different types of molecules without using too many atom types. When multiple force fields are mixed, the parameters are taken following a user-defined precedence order of the force fields and are identified by a set of instant atom types. The instant atom types are made by combining the original atom types and the label of the force field from which the parameters are taken. For example, if we have force field A developed for polymers and force field B developed for ILs, and both force fields have the atom type c_4 (a generic sp3 hybrid carbon), then we would have two instant atom types, c_4A and c_4B, linked to the parameters in force fields A and B, respectively. The precedence order ensures that the correct parameters are called. However, if the coverage shares the same parameters (by transferring), then the precedence order does not matter. The parametrization strategy is similar to that used in TEAM force fields.52 First, the charge parameters were obtained by fitting to the CHelpG charges calculated on isolated ions optimized using the B3LYP/6-31G(d) method.55 Most parameters of the bonded terms were transferred from the TEAM-AMBER force field52 and the recently updated parameter set of alkyl groups,41 and the new parameters were optimized by fitting the structures of the ions. The dihedral parameters were further refined by fitting the conformational energies of the ions calculated at the M06/def2-SVP//madef2-TZVP level of theory.56−59 The LJ parameters were optimized by fitting the shear viscosity and equilibrium density. A total of 54 newly introduced LJ parameters were optimized using 90 pieces of experimental density and viscosity data in an iterative manner. Generally speaking, the radius parameters were optimized first to reproduce densities with energy parameters fixed; then, energy parameters were optimized targeting at viscosities with radius parameters fixed. After the LJ parameters were updated in each cycle, the torsion parameters were readjusted to count the couplings between torsion and nonbonded terms. The charge parameters were scaled by a factor of 0.8 in liquid simulations. MD simulations were carried out for parametrization and validation. The cutoff of vdW and electrostatic interactions was 1.2 nm, supplemented by the continuum correction for vdW energy60 and the smooth particle mesh Ewald (SPME)61 using a spline order of 4 and a Fourier spacing of 0.12 nm for electrostatic energy. The time step of MD integration was 1 fs. The Nosé−Hoover thermostat62,63 with a time constant of 0.5 ps was used to control the temperature. Berendsen64 and Parrinello−Rahman65 barostats with a time constant of 5 ps were used to control the pressure at the equilibration and production stages, respectively. The LINCS66 algorithm was used to constrain the length of the bond involving the hydrogen atom. The densities were calculated by averaging over trajectories of NPT simulations, with the simulation box containing ∼3000 atoms in the dimension of about 3.3 × 3.3 × 3.3 nm3. The calculation of viscosity was performed on doubled sizes. A rectangular box with dimensions of about 3.3 × 3.3 × 6.6 nm3 containing ∼6000 atoms was used for the PP calculation, and a cubic box containing the same number of atoms was used for the GK calculation. The structures were
III. RESULTS AND DISCUSSION Molecular Properties. The prediction of molecular structures and energies is a prerequisite for the prediction of condensed-phase properties because the molecular properties impact the molecular packing and interactions. Among the molecular properties, the conformational energies and structures are the most challenging to predict by using a force field. As examples, we list the energy profiles of four representative internal torsions of 1-ethylimidazolium, 1ethylimidazolium, 1-ethyl-1-methylpyrrodinium, and triflate in Figure 2. The force-field curves agree well with that
Figure 2. Comparison of energy profiles of four representative torsions calculated by quantum mechanics and force field. Hydrogen atoms are omitted for clarity.
calculated with M06/ma-def2-TZVP method.57−59 The rotation of the ethyl group about the five-member (1ethylimidazolium) or six-member (1-ethylimidazolium) aromatic ring exhibits a low energy barrier (5 mPa·s. We applied the linear extrapolation to estimate the zero-shear viscosities in this work. As shown in Figure 5, within the error bars, the GK and PP methods yield similar results in most cases. For low-viscosity fluid (viscosity
