Molecular Simulation of Guanidinium-Based Ionic Liquids - American

5658

J. Phys. Chem. B 2007, 111, 5658-5668

Molecular Simulation of Guanidinium-Based Ionic Liquids Xiaomin Liu,†,‡ Guohui Zhou,†,‡ Suojiang Zhang,*,† Guangwen Wu,§ and Guangren Yu†,‡ State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, 100080, Beijing, China, Graduate UniVersity of the Chinese Academy of Sciences, 100049, Beijing, China, and Center for Molecular Simulation, Swinburne UniVersity of Technology, Vic 3122, Australia ReceiVed: December 22, 2006; In Final Form: March 14, 2007

A new systematic all-atom force field was developed for cyclic guanidinium-based ionic liquids (ILs) based on the AMBER force field. Optimized molecular geometries and equilibrium bond lengths and angles were obtained by ab initio calculations, and charges were allocated to each atom center by fitting the ab initio electrostatic potential. Molecular dynamics simulations were performed for eleven kinds of ILs that are comprised of NO3- anions and cyclic guanidinium-based cations. Validation was carried out by comparing our simulated densities with experimental and calculated data from the literature. Transport properties such as self-diffusion coefficients, viscosities, and conductivities were calculated by molecular dynamic simulation, and their dependence on the length of the alkyl chains of cyclic guanidinium-based cations are discussed. Radial distribution functions and spatial distribution functions were investigated to depict the microscopic structures of the ILs, and the relationship between their properties and microstructures is also discussed.

E(rN) )

1. Introduction Ionic liquids (ILs) have emerged as promising alternative media for the replacement of conventional organic solvents.1-8 Guanidinium ILs have found many useful applications in dyesensitized solar cells,9 biological and molecular recognition,10,11 and as energetic materials.12,13 However, studies on such kinds of new energetic ILs are quite limited compared to the extensive studies on imidazolium ILs. There are extremely strong demands for exploring the relationship between microstructures and properties of guanidinium ILs in order to design task-specific ILs with great versatility of chemical and physical properties.14-17 To extend our previous studies on acyclic guanidinium-based ILs18 to a wider region covering cyclic ones, in this work we focused on molecular dynamic simulations on cyclic guanidinium-based ILs. In this work, a new force field for cyclic guanidinium-based ILs was exploited and molecular dynamics simulation was performed for eleven kinds of ILs that are comprised of NO3anions and cyclic guanidinium-based cations. Validation was carried out by comparing simulation densities with literature. Transport properties such as viscosity, conductivity, and selfdiffusion coefficients were studied for these cyclic guanidiniumbased ILs in order to fill up the void in both experimental and simulation studies. Radial distribution functions (RDFs) and spatial distribution functions (SDFs) were also studied to depict the microstructures and investigate the relationship between microstructure and macroproperties of these ILs. 2. Force Field Development In this work, the standard molecular mechanics force field was used with the functional form as follows:14,19 * Address correspondence to this author. E-mail: [email protected]. † Institute of Process Engineering, Chinese Academy of Sciences. ‡ Graduate University of the Chinese Academy of Sciences. § Swinburne University of Technology.

Kr(r - r0)2 + ∑ Kθ(θ - θ0)2 + ∑ bonds angles

∑

torsions

Kφ 2

N

(1 + cos(nφ - γ)) +

N

∑ ∑ i)1 j)i+1

{ [( ) ( ) ] } σij

4ij

rij

12

-

σij rij

6

+

qiqj rij

(1)

where the first three terms in eq 1 represent the bonded interactions: bonds, angles, and torsions. The nonbonded interactions include van der Waals (VDW, in the Lennard-Jones (LJ) 6-12 form) and Coulombic interactions of atom-centered point charges. Electrostatic and VDW interactions were calculated for the atoms in different molecules or in the same molecule separated by more than three bonds. The nonbonded interactions separated by exactly three bonds (1-4 interactions) are reduced by a scale factor, which was optimized as 1/2 for VDW and 1/1.2 for electrostatic interactions.20 The LJ parameters for unlike atoms were obtained from the Lorentz-Berthelot (LB) combining rule.14 The force field parameters were modified as that in the ref 19. The isolated ions were optimized by using Gaussian 03 package at the HF/6-31+G(d) level. The minimized structures were used to set equilibrium bond lengths (r0) and angles (θ0). Once the optimized geometries were obtained, the vibration frequencies were checked to ensure no negative frequencies existed and verify the ion energy of a true minimum. The structures and assignment of the atom types are shown in Figure 1. The force constants for bonds and angles were adjusted by fitting the vibration frequency data based on AMBER99. As the Hartree-Fock (HF) method overestimates the frequencies by approximately 10%, a scaling factor of 0.9 is applied for HF calculations. Torsion parameters were also adjusted by fitting torsion energy profiles. In this work, the ab initio torsion energy profiles were obtained at the MP2/6-31+G(d)//HF/6-31+G(d) level. The scans of energy profiles were started from the

10.1021/jp068849a CCC: $37.00 © 2007 American Chemical Society Published on Web 05/01/2007

Guanidinium-Based Ionic Liquids

J. Phys. Chem. B, Vol. 111, No. 20, 2007 5659

Figure 1. Structures and atom types for cations and anion.

minimum and the dihedral angle considered varied in a step of 10°. Two of the energy barrier fittings of CT-N2-CR-N* and N2-CR-N*-CK for F2 are presented in Figure 2. The one-conformation, two-step restraint electrostatic potential method21-25 was used to derive the atom charges by fitting the electrostatic potential generated from QM calculations at the

B3LYP/6-31+G(d) level. The force field parameters used in this work are listed in Table 1. 3. Simulations Details Molecular dynamics (MD) simulations for the eleven ILs were performed with the standard periodical boundary conditions with

5660 J. Phys. Chem. B, Vol. 111, No. 20, 2007

Liu et al. are higher than normal atmosphere temperature, only 4 kinds of ILs densities under the standard condition were measured in the literature.3 Our simulation densities, experimental data, and the calculated results at 298 K from ref 3 are presented in Table 2 for comparison. The simulated results are systematically smaller than the experimental and calculated ones. The fact that the densities are typically 5% below the reference data was possibly caused by the LJ potential not being realistic enough to describe interactions within these systems accurately. As a result, the agreement of simulated densities with the reference could be improved by fine-tuning the parameters for LennardJones parameters. Another reason may be that the charges derived from gas-phase calculations do not necessarily represent charges in the liquid phase, particularly when the liquid is composed of associative ions prone to polarizability and charge transfer. The simulation densities of these ILs at equilibrium temperatures for sampling are also listed in Table 2. It can be seen that the densities increase with reducing the length of alkyl chains for cations. Taking the A/A2/A3 series, for example, the density changes from 1.05 to 1.07, and finally 1.14 when the alkyl chain varies from butyl to propyl, and finally methyl. The same behavior is noticed when B/B2, C/C2, and F/F2 are taken into account. For imidazolium-based IL, similar behavior is also observed.28 The reason may be that the shorter the alkyl chains, the closer the packing of ionic pairs. 4.2. Transport Properties. A way to characterize the particle diffusion dynamics in a dense fluid is through the mean-square displacement (MSD). The coefficient of self-diffusion for a liquid can be calculated by using the Einstein relation17

d 1 Dself ) lim 〈|ri(t) - ri(0)|〉 6 tf∞ dt Figure 2. Energy barrier fitting of CT-N2-CR-N* and N2-CRN*-CK for F2. Quantum mechanic (lines) and molecular mechanic (dots).

use of the M.DynaMix package.26 Each system contained 192 ion pairs. A 2/0.2 fs multiple time-step algorithm was adopted. The intramolecule forces were cut off at 15 Å, while the longrange forces including LJ and Coulombic interactions were cut off at half of the simulation boxes, and Ewald summation was implemented for the latter. All simulations started from a facecentered cubic (fcc) lattice with a low density of 0.2 g/cm3. After the system was relaxed in the NVE ensemble for a few MD steps to remove the possible overlapping in the initial configuration, the Nose-Hoover27 NpT ensemble was adopted with coupling constants of 700 and 100 fs, and the system density gradually increased to a reasonable value. The system was equilibrated for at least 200 ps, making sure that the system configuration was stable. Each production phase lasted for 1 ns at 450, 420, or 320 K under 1.0 atm. For A/A2/A3, C/C2, and E the temperature was 450 K, for B/B2 and D 320 K, and for F under 420 K. The equilibrium temperatures were based on their melting points or glass phase transition points.3 For convenient comparison, the same temperature was used for the anions with different alkyl chain lengths. Trajectories were dumped in an intervals of 10 fs for further analysis. All of the simulations were performed at a 16-node Linux cluster. Each node contains four Itanium II processors operating at 1.3 GHz. 4. Results 4.1. Liquid Densities. As most of the melting points or glass phase transition points of these cyclic guanidinium-based ILs

(2)

the quantity in braces is the ensemble-averaged MSD of the molecules and ri is the vector coordinate of the center of mass of ion i. The start points and time regions for linear fitting of MSDs differ from each other in the literatures.17,29,30 It is reported that the error for estimating the self-diffusion coefficient (SDCs) by using MSD calculations is (1) the shorter the run length used for fitting, the greater will be the overestimation of self-diffusion coefficient and (2) fitting should avoid the initial and final regions. 29 In this work, trajectories were dumped for 1 ns and the range 200-600 ps is used for estimating the slope by linear regression. The log(MSD) - log(t) for FNO3 was plotted in Figure 3 c, a short-term rattling or caging in the first few picoseconds, and then a subdiffusive regime is observed. The proper selfdiffusion is achieved after 200 ps.15 The MSDs for the cation and anion of FNO3/F2NO3 at 420 K and 1.0 atm are also shown in Figure 3a,b. Overlaid in the figure is the best linear fit in the range 200-600 ps. From the slope of this line and according to eq 2, the SDCs of the cation and anion are 1.21 × 10-11 and 0.95 × 10-11 m2 s-1 for FNO3, and 0.59 × 10-11 and 0.72 × 10-11 m2 s-1 for F2NO3, respectively. The SDCs for all the ILs at sampling temperatures are shown in Table 3. It can be seen that nearly all the MSDs increase with the length of the alkyl chain for increasing cations. The unexpected SDCs for the A series may be due to the uncertainties in the calculation of MSDs. Comparing with the SDC for water at 2.30 × 10-9 m2 s-1 at 298 K,31 from the fact that the SDCs for these ILs are about 2 or 3



TABLE 1: Force Field Parameters for Guanidinium-Based ILs σi (Å)

atom CT CT* CA CR CK NA NB N* N N2 OS O2 H HC H1 H2 H5

i (kJ/mol)

3.3997 3.3997 3.3997 3.3997 3.3997 3.2500 3.2500 3.2500 3.2500 3.2500 3.0000 2.9603 1.0691 2.6495 2.4714 2.2932 2.4215

description sp3

0.4577 0.4577 0.3598 0.3598 0.3598 0.7133 0.7133 0.7133 0.7113 0.7133 0.7133 0.8786 0.0657 0.0657 0.0657 0.0657 0.0628

carbon sp3 carbon in the cycle of guanidine sp2 aromatic carbon sp2 aromatic in 5-membered ring next to two nitrogen atoms sp2 carbon in 5-membered aromatic between N and N-R sp2 nitrogen in aromatic rings with hydrogen attached sp2 nitrogen in 5-membered ring with lone pair sp2 nitrogen in 5-membered ring with carbon substituent sp2 nitrogen in amides sp2 nitrogen of aromatic amines and guanidinium ions sp3 oxygen in ethers sp2 oxygen in anionic acids H attached to N H attached to aliphatic carbon with no electron-withdrawing substituents H attached to aliphatic carbon with one electron-withdrawing substituent H attached to aliphatic carbon with two electron-withdrawing substituent H attached to aliphatic carbon with two electronegative neighbors

qi (e) atoma

A

A2

CR N*(m) N*(a) N2 CT*(m) CT*(a) H1(m*) H1(a*) CT(m) CT(a) H1(m) H1(a) CT(me) H1(me) CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3) CT(a4) HC(a4)

0.0660 0.0052 0.0161 -0.0058 -0.0373 -0.0436 0.0952 0.0954 -0.1608 -0.2143 0.1001 0.1159 -0.0893 0.0810 -0.0213 0.0521 0.0021 0.0174 0.0298 0.0129 -0.0639 0.0267

0.0604 0.0078 0.0067 -0.0132 -0.0473 -0.0296 0.0961 0.0948 -0.1540 -0.1906 0.0978 0.1112 -0.0772 0.0778 -0.0323 0.0570 0.0416 0.0209 -0.0827 0.0374

qi (e)

qi (e) A3

atomb

B

0.0643 0.0149 0.0150 0.0194 -0.0397 -0.0397 0.0958 0.0958 -0.2141 -0.2141 0.1180 0.1180 -0.1174 0.0897 -0.1174 0.0896

CR N*(h) N*(a) N2 CT*(h) CT*(a) H1(h*) H1(a*) CT(h) CT(a) H1(h) H1(a) H CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3)

0.2000 -0.0843 -0.0635 -0.2247 -0.0506 -0.0802 0.1127 0.1150 -0.1842 -0.1503 0.1206 0.1095 0.2383 -0.0669 0.0897 0.0235 0.0379 -0.0766 0.0396

qi (e) B2

atoma

C

C2

0.2047 -0.1070 -0.0308 -0.2438 -0.0246 -0.1164 0.1119 0.1224 -0.1490 -0.1825 0.1093 0.1189 0.3092 -0.2401 0.1425

CA NA(m) NA(a) N2 CT*(m) CT*(ma) CT*(a) H1(m*) HC(ma) H1(a*) CT(m) CT(a) H1(m) H1(a) CT(me) H1(me) CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3) CT(a4) HC(a4)

0.0550 -0.0339 0.0440 0.0122 -0.0392 -0.0232 -0.0934 0.0875 0.0675 0.0947 -0.1125 -0.1804 0.0874 0.0989 -0.0681 0.0685 -0.0142 0.0477 -0.0108 0.0172 0.0310 0.0127 -0.0678 0.0276

0.0608 -0.0259 0.0276 -0.0080 -0.0332 -0.0229 -0.0811 -0.0854 0.0659 0.0921 -0.1062 -0.1634 0.0846 0.0962 -0.0642 0.0686 -0.0322 0.0563 -0.0229 0.0231 -0.0738 0.0352

atom

NO3-

qi (e)

qi (e)

qi (e)

atomb

D

atoma

E

atoma

F

F2

CA NA(h) NA(a) N2 CT*(h) CT*(ha) CT*(a) H1(h*) HC(ha) H1(a*) CT(h) CT(a) H1(h) H1(a) H CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3)

0.0903 0.0204 0.0042 -0.2792 -0.0440 -0.0190 -0.0483 0.0866 0.0625 0.0837 -0.1793 -0.0945 0.1087 0.0820 0.2514 -0.0069 0.0743 0.0178 0.0299 -0.0882 0.0405

CA NA(m) NA(a) N2 CT*(m) CT*(a) OS H2(m*) H2(a*) CT(m) CT(a) H1(m) H1(a) CT(me) H1(me) CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3)

0.1473 0.0289 -0.1232 0.0046 0.0482 0.0716 -0.2890 0.1213 0.1173 -0.2208 -0.1225 0.1090 0.1093 -0.0868 0.0854 -0.1011 0.0621 0.0627 0.0135 -0.0656 0.0354

CR N*(m) N*(a) N2 NB CK H5 CT(m) CT(a) H1(m) H1(a) CT(me) H1(me) CT(a1) H1(a1) CT(a2) HC(a2) CT(a3) HC(a3) CT(a4) HC(a4)

0.2099 0.2864 0.0072 -0.1366 -0.4349 0.1460 0.1908 -0.0692 -0.2130 0.0895 0.1315 -0.2056 0.1183 -0.0873 0.0836 0.0300 0.0033 0.0521 0.0077 -0.0670 0.0280

0.1001 0.3196 0.0358 -0.1235 -0.4228 0.1369 0.1946 -0.1330 -0.1519 0.1080 0.1177 -0.0737 0.0868 -0.1275 0.1027

N O2

1.0335 -0.6778


Liu et al.

TABLE 1. (continued) Bonded Parameters 1c bond

2

Kr (kJ/mol/Å )

r0 (Å)

source

bond

Kr (kJ/mol/Å2)

r0 (Å)

source

1548.08 1548.08 2012.50 1422.56 1410.01 1297.04

1.336 1.322 0.996 1.082 1.463 1.529

this work this work AMBER AMBER AMBER AMBER

CT-CT CT-HC CT-H1 CT-N* CT-N2 CR-N*

1297.04 1422.56 1422.56 1255.20 1255.20 1464.40

1.528 1.086 1.081 1.458 1.470 1.333

structures A/A2/A3A, B/B2B AMBER CR-N2A AMBER CR-N2B AMBER N2-H this work CT*-H1 this work CT*-N* this work CT*-CT*

CT-CT CT-HC CT-H1 CT-NA CT-N2C CT-N2D CA-NA

1297.04 1422.56 1422.56 1410.01 1410.01 1410.01 1380.72

1.528 1.086 1.081 1.464 1.468 1.340 1.331

structures C/C2C, DD AMBER CA-N2 AMBER N2-H AMBER CT*-HC AMBER CT*-H1 AMBER CT*-NA AMBER CT*-CT* this work

1338.88 2020.87 1422.56 1422.56 1410.01 1297.04

1.348 0.995 1.084 1.083 1.470 1.519

this work this work AMBER AMBER AMBER AMBER

CT-CT CT-HC CT-H1 CT-H2 CT-NA

1297.04 1422.56 1422.56 1443.48 1410.01

1.529 1.085 1.081 1.081 1.467

structure E AMBER CT-N2 AMBER CT-OS AMBER CA-NA this work CA-N2 AMBER

1410.01 1338.88 1380.72 1338.88

1.470 1.382 1.336 1.337

AMBER AMBER this work this work

CT-CT CT-HC CT-H1 CT-N* CT-N2 CR-N*

1297.04 1422.56 1422.56 1410.01 1410.00 1380.72

1.529 1.090 1.080 1.464 1.464 1.333

structures F/F2 AMBER CR-N2 AMBER CK-H5 AMBER CK-NB AMBER CK-N* AMBER NB-N* this work

1338.88 1598.29 2213.34 1840.96 1853.51

1.340 1.070 1.270 1.370 1.358

this work this work AMBER AMBER this work

N-O2

2133.84

1.259

structure NO3this work Bonded Parameters 2c

angle

2

Kθ (kJ/mol/ rad )

θ0 (deg)

source

angle

Kθ (kJ/mol/ rad2)

θ0 (deg)

source

146.44 125.52 217.57 209.20 209.20 209.20 217.57 292.88 292.88 238.49 209.20

108.75 109.05 110.15 111.45 124.35 110.20 109.70 120.30 110.05 112.55 102.25

AMBER this work this work this work this work AMBER this work AMBER AMBER this work AMBER

CT-CT-CT CT-CT-HC CT-CT-H1 CT-CT-N2 CT-N*-CR CT-N2-CT CT-N2-H CT-N2-CR A CT-N2-CR B CR-N2-H HC-CT-HC

167.36 209.20 209.20 334.72 292.88 167.36 209.20 209.20 209.20 209.20 146.44

111.60 110.20 109.75 110.85 126.55 116.35 114.90 121.80 129.45 115.15 107.60

structures A/A2/A3A, B/B2B AMBER H1-CT-H1 AMBER H1-CT*-H1 AMBER N*-CT-H1 AMBER N*-CR-N* AMBER N*-CR-N2 this work N*-CT*-H1 AMBER N2-CT-H1 AMBER CT*-N*-CT AMBER CT*-N*-CR AMBER CT*-CT*-H1 AMBER CT*-CT*-N*

CT-CT-CT CT-CT-HC CT-CT-H1 CT-CT-N2 CT-NA-CA CT-N2-CT CT-N2-H CT-N2-CA C CT-N2-CA D CA-N2-H HC-CT-HC HC-CT*-HC H1-CT-H1

167.36 175.73 175.73 334.72 292.88 209.20 146.44 209.20 209.20 146.44 146.44 146.44 171.54

111.45 110.30 109.40 111.50 122.45 117.00 114.50 121.50 126.80 115.60 107.65 107.70 108.40

structures C/C2C, DD AMBER H1-CT*-H1 this work NA-CT-H1 this work NA-CA-NA AMBER NA-CA-N2 AMBER NA-CT*-H1 AMBER N2-CT-H1 this work CT*-NA-CT AMBER CT*-NA-CA AMBER CT*-CT*-HC this work CT*-CT*-H1 AMBER CT*-CT*-NA AMBER CT*-CT*-CT* this work

146.44 217.57 292.88 292.88 209.20 209.20 292.88 292.88 209.20 209.20 209.20 167.36

107.75 110.00 119.80 120.05 108.25 109.45 115.95 120.95 110.35 111.20 110.15 107.80

AMBER this work AMBER AMBER AMBER AMBER AMBER AMBER AMBER AMBER AMBER AMBER

CT-CT-CT CT-CT-HC CT-CT-H1 CT-CT-N2 CT-NA-CT CT-NA-CA CT-N2-CT CT-N2-CA CT-OS-CT HC-CT-HC

167.36 230.12 230.12 334.72 292.88 292.88 209.20 209.20 251.04 154.81

111.10 110.30 109.70 112.50 115.40 120.70 116.60 121.60 110.90 107.60

structure E AMBER H1-CT-H1 this work H2-CT-H2 this work NA-CT-H1 AMBER NA-CT-H2 AMBER NA-CT-OS AMBER NA-CA-NA AMBER NA-CA-N2 AMBER N2-CT-H1 AMBER OS-CT-H2 this work

154.81 167.36 230.12 238.49 209.20 292.88 292.88 230.12 238.49

108.80 109.70 109.90 109.10 109.70 117.70 121.10 109.60 109.60

this work this work this work this work AMBER AMBER AMBER this work this work



TABLE 1. (continued) Bonded Parameters 2 angle

Kθ (kJ/mol/rad2)

θ0 (deg)

CT-CT-CT CT-CT-HC CT-CT-H1 CT-CT-N2 CT-N*-CR CT-N*-CK CT-N*-NB CT-N2-CT CT-N2-CR CR-N*-CK CR-N*-NB

167.36 209.20 209.20 334.72 292.88 292.88 292.88 209.20 209.20 292.88 292.88

112.10 109.90 110.10 113.70 128.90 125.60 118.70 116.30 121.70 105.90 111.20

structures F/F2 AMBER HC-CT-HC AMBER H1-CT-H1 AMBER NB-CK-H5 AMBER NB-CK-N* AMBER N*-CT-H1 AMBER N*-CR-N* AMBER N*-CR-N2 AMBER N*-CK-H5 AMBER N*-NB-CK AMBER N2-CT-H1 AMBER

O2-N-O2

564.84

120.00

structure NO3this work

source

angle

Kθ (kJ/mol/ rad2)

θ0 (deg)

source

154.81 156.90 167.36 292.88 230.12 292.88 292.88 171.54 292.88 230.12

107.30 109.30 125.00 111.70 109.10 106.10 127.00 123.30 105.10 110.20

this work this work this work AMBER this work AMBER AMBER this work AMBER this work

Bonded Parameters 3 γ (deg)

Kφ (kJ/mol)

n

CT-CT-CT-CT CT-CT-CT-HC CT-CT-CT-H1 CT-CT-CT-N2 CT-CT-N2-CT CT-CT-N2-CR CT-CT-N2-H HC-CT-CT-HC HC-CT-CT-H1 H1-CT-N*-CR H1-CT-N*-CT* H1-CT-N2-CT H1-CT-N2-CR H1-CT-N2-H

0 0 0 0 0 180 180 0 0 0 0 0 0 180

0.753 0.669 0.653 0.653 3.849 2.092 2.092 0.628 0.653 0.837 0.628 1.548 0.209 -0.418

3 3 3 3 3 3 3 3 3 3 3 3 3 3

structures A/A2/A3, B/B2 AMBER H1-CT*-N*-CT AMBER H1-CT*-N*-CR AMBER H1-CT*-CT*-H1 AMBER N*-CR-N2-H this work N*-CR-N*-CT this work N*-CR-N*-CT* this work N*-CR-N2-CT AMBER N*-CT*-CT*-H1 AMBER N*-CT*-CT*-N* this work N2-CT-CT-HC this work N2-CR-N*-CT this work N2-CR-N*-CT* this work CT*-CT*-N*-CT this work CT*-CT*-N*-CR

CT-CT-CT-CT CT-CT-CT-HC CT-CT-CT-H1 CT-CT-CT-N2 CT-CT-N2-CT CT-CT-N2-CA CT-CT-N2-H HC-CT-CT-HC HC-CT-CT-H1 HC-CT*-CT*-H1 H1-CT-NA-CA H1-CT-NA-CT* H1-CT-N2-CT H1-CT-N2-CA H1-CT-N2-H

0 0 0 0 0 0 0 0 0 180 0 0 0 0 0

1.255 0.460 0.653 0.653 0.000 0.000 0.418 0.418 0.653 0.653 0.418 1.130 0.837 0.209 0.418

3 3 3 3 3 3 3 3 3 3 2 3 3 3 3

CT-CT-CT-HC CT-CT-CT-H1 CT-CT-CT-N2 CT-CT-N2-CT CT-CT-N2-CA HC-CT-CT-HC HC-CT-CT-H1 H1-CT-NA-CT H1-CT-NA-CA H1-CT-N2-CT H1-CT-N2-CA

0 0 0 0 0 0 0 0 0 180 180

0.464 0.653 1.004 2.092 0.418 0.628 0.653 0.753 0.753 0.418 0.418

3 3 3 3 3 3 3 3 3 3 3

torsion

γ (deg)

Kφ (kJ/mol)

n

source

0 0 180 180 180 0 180 0 180 0 180 0 0 0

2.301 2.092 1.071 5.858 11.715 7.113 4.184 6.234 6.276 0.653 11.715 7.531 8.368 6.276

2 2 2 2 2 2 2 3 3 3 2 2 2 2

this work this work this work this work this work AMBER this work this work this work AMBER this work this work this work this work

structures C/C2, D this work H1-CT*-NA-CT this work H1-CT*-NA-CA AMBER NA-CA-NA-CT AMBER NA-CA-NA-CT* AMBER NA-CA-N2-CT AMBER NA-CA-N2-H this work NA-CT*-CT*-HC this work N2-CT-CT-HC AMBER N2-CA-NA-CT AMBER N2-CA-NA-CT* this work CT*-CT*-NA-CT this work CT*-CT*-NA-CA this work CT*-CT*-CT*-H1 this work CT*-CT*-CT*-NA this work

0 0 180 180 180 180 0 0 180 180 0 0 0 0

-1.674 -2.092 11.715 10.042 14.644 6.276 0.695 0.653 10.042 10.042 0.418 0.418 0.611 2.745

3 3 2 2 2 2 3 3 2 2 3 3 3 3

this work this work this work AMBER this work this work this work AMBER AMBER AMBER this work this work this work this work

structure E this work H2-CT-NA-CT AMBER H2-CT-NA-CA this work H2-CT-OS-CT this work NA-CT-OS-CT this work NA-CA-NA-CT AMBER NA-CA-N2-CT AMBER N2-CT-CT-HC this work N2-CA-NA-CT this work OS-CT-NA-CT this work OS-CT-NA-CA this work

0 0 0 180 180 180 0 180 180 0

0.418 0.418 0.347 1.602 10.042 11.715 0.653 11.715 6.276 6.276

2 2 3 3 2 2 3 2 2 2

this work this work this work AMBER AMBER this work AMBER this work this work this work

180 180 180 180 180

15.481 2.929 15.481 2.929 29.706

2 2 2 2 2

this work this work this work this work this work

source

torsion

structures F/F2 CT-CT-CT-CT CT-CT-CT-HC CT-CT-CT-H1 CT-CT-CT-N2 CT-CT-N2-CT

0 0 0 0 0

0.753 0.669 0.653 0.653 2.092

3 3 3 3 3

proper torsions AMBER H5-CK-NB-N* AMBER H5-CK-N*-CT AMBER H5-CK-N*-CR AMBER NB-CK-N*-CT this work NB-CK-N*-CR


Liu et al.

TABLE 1. (continued) torsion

γ (deg)

Kφ (kJ/mol)

n

source

torsion

γ (deg)

Kφ (kJ/mol)

n

source

structures F/F2 CT-CT-N2-CR CT-N*-NB-CK CK-NB-N*-CR HC-CT-CT-HC HC-CT-CT-H1 H1-CT-N*-CR H1-CT-N*-CK H1-CT-N*-NB H1-CT-N2-CT H1-CT-N2-CR

0 0 180 0 0 0 0 0 0 0

2.092 10.460 29.706 0.628 0.653 0.000 0.837 0.502 0.000 1.255

3 2 2 3 3 3 3 3 3 3

proper torsions this work N*-CR-N*-CT this work N*-CR-N*-CK this work N*-CR-N*-NB AMBER N*-CR-N2-CT AMBER N*-NB-CK-N* AMBER N2-CT-CT-HC this work N2-CR-N*-CT this work N2-CR-N*-CK AMBER N2-CR-N*-NB this work

180 180 180 180 180 0 180 180 180

11.715 29.706 29.706 11.715 29.706 0.653 11.715 22.280 22.280

2 2 2 2 2 3 2 2 2

this work this work this work this work this work AMBER this work this work this work

CT-CR-CK-N* CT-CR-NB-N*

180 180

4.184 4.184

2 2

improper torsions this work H5-N*-NB-CK this work

180

10.460

2

this work

a The letter in parentheses indicates the atom is located close to methyl (m) or alkyl (a) side, (ma) means the distance for the CT*/HC is the same for methyl or alkyl, (m*) or (a*) stand for hydrogen connected with CT*, respectively, (me) indicates the atom for methyl, and (a1), (a2), (a3), and (a4) stand for the atom in alkyl. b The letter in parentheses indicate the atom is located close to atom hydrogen (h) or alkyl (a) side, (ha) means the distance for the CT*/HC is the same for methyl or alkyl, (h*) or (a*) stand for hydrogen connected with CT*, respectively, and (a1), (a2), (a3), and (a4) stand for the atom in alkyl. c Superscript A, B, C, and D stand for the corresponding parameters for A, B, C, or D series.

Figure 3. Average mean square displacement for (a) FNO3 and (b) F2NO3 and (c) log-log plot of MSD for FNO3

orders of magnitude smaller than that of water, it is obvious that the ILs are much “slow” than water. The lower diffusion constant value of the ILs could be explained by the strong interactions between the cation and anions, especially for the hydrogen-bonding interactions,29 which is discussed in more detail in part 4.3. Applying the Stokes-Einstein relation to these systems, the viscosities ηILfor the ILs can be estimated by the equation17

ηIL ) ηH2O

DH2O DIL

(3)

where DIL ) 1/2(Dcation + Danion) and ηH2O ) 0.9 mPa s.17,32 All the results listed in Table 3 are consistent with the fact that the viscosities for most of the ILs are much higher than those of water. Another crude estimation for these ILs is the molar conductivities ΛIL, which can be obtained from the Nernst-Einstein equation16,33

ΛIL )

NAe2 (D + Danion) kBT cation

(4)



Figure 4. Center-of-mass radial distribution functions (RDFs) of the cation-anion.

where NA is Avogadro’s number and e is the electron charge. The estimated molar conductivities ΛIL and conductivities κIL are all listed in Table 3. Together with the MSDs and viscosities, it is observed that the activity is reduced as the length of the alkyl chain is decreased for these cyclic guanidinium-based ILs. 4.3. Microstructures. The center-of-mass radial distribution functions (RDFs) were used to study the microstructure of the ILs. RDFs for cation-cation, cation-anion, and anion-anion of the eleven ILs are shown in Figure 4. Molten salts are characterized by strong oscillations in the radial distribution functions and are considered as strongly coupled ionic systems.34 As can be seen in Figure 4, all the oscillations extend beyond 18 Å, which is approximately the half-length of the simulation box. It is shown that electrostatic interactions are the major component of the system energy. All the peak heights of the RDFs are below 2.5, compared with that for imidazolium-based ILs of about 3.0.35 It may be that the interactions of ionic pairs for the cyclic guanidinium-based ILs are a litter weaker than that for imidazolium-based ones. It can be seen that some of the curves are not very smooth, for example, for F/F2-NO3the peaks even split into two, and the reason may be that anions most likely distribute around cations in different regions, and the distances for the regions to the center of cations are different.

The coordination numbers in the first solvent shell were obtained by integrating the RDFs from zero to the first minimum, which can be used for analyzing the organization of the bulk liquid. It is simply the average number of specific sites or atoms within a sphere of radius r about some other central site or atom. The coordination number, N, can be calculated by14

N ) 4π

∫0R

min1

Fg(r)r2 dr

(5)

where F is the number density. The first solvent shell locations and the coordination numbers for the eleven kinds of ILs are listed in Table 4, which shows that each cation is surrounded by five to seven anions. Although RDFs of different ILs with different alkyl chain lengths are somewhat similar, it is interesting to notice that the solvent shells locations a further apart from the cation centers with e decrease in the alkyl chain length for these ILs. However, the case is contrary for alkylimidazolium family.14 Another finding is that the coordination number increases with a decrease in the length of the alkyl chain, which may be explained by the space steric effect.


Liu et al.

Figure 5. Hydrogen-oxygen radial distribution functions between cation and anion. The site and structure of the hydrogen atoms (a, b, c, d, e, f) is shown in Figure 1.

A more detailed liquid structure can be described by the sitesite pair RDFs. The RDFs for cations with different alkyl chain lengths are very similar, only the RDFs for the oxygen atom in NO3- with different hydrogen atoms of (a) A, (b) B, (c) C, (d) D, (e) E, and (f) F are shown in Figure 5. The different hydrogen atoms involved in the RDFs are marked in Figure 1. In general, the oxygen atom in NO3- prefers to distribute around H1/H2/H5 at the “a” site for A, C, E, and F; however, when there is an H replacing the methyl in the “c” site (B and D), an obviously higher peak occurs, which implies that the activity of the H increases remarkably. The first maximum locations of the RDFs for H1/H2/H5 or H are 2.60, 1.96, 2.60, 2.04, 2.44, and 2.36 Å for A to F, respectively. No mater which systems, the first peak occurs at 2.6-2.7 Å for all the other site-site pair RDFs. It is also found that the interaction between O and HC at the end of the alkyl chain is rather weak. As is said by Marco Fioroni and Danilo Roccatano, hydrogenbonding interactions between the molecules will influence the transport properties,29 such as the diffusion constant and viscosity. The lower value of the diffusion constants and the higher viscosities for B, D, and F series could be explained by

the strong hydrogen-bonding interactions between H-O, H-O, and H5-O, respectively. To depicture more visual structures of these ILs, the spatial distribution functions (SDFs) that stand for the three-dimensional probability distributions of anion and cation were studied. In this work the SDFs for the anion around the cation were visualized by the software package gOpenMol36 and are shown in Figure 6. In each case, the coordinate lengths in the x, y, and z directions are chosen to be 20 Å, and the red and yellow contour surfaces are drawn at 10 and 6 times the average density, respectively. As is shown in Figure 6, the probability distributions of NO3- around cations with different alkyl chain lengths are very similar. For A and C series, the main regions for anions around the cation are above/below the cyclic guanidinium plane, but the anions are most likely distributed around the hydrogen atoms or alkyl chain for B and D. The anions in the first solvent shell are all distributed around the alkyl for E, which is caused by the repulsion of the oxygen atom of the cation to anions. Four regions are found for F and F2. It is interesting to notice that the SDFs between A2 and B


J. Phys. Chem. B, Vol. 111, No. 20, 2007 5667 TABLE 3: Dynamical Properties for Guanidinium-Based ILs

IL

Dcation (10-11 m2/s)

Danion (10-11 m2/s)

DIL (10-11 m2/s)

ηIL (10-3 Pa s)

ΛIL (10-4 S m2/mol)

κIL (S/m)

ANO3 A2NO3 A3NO3 BNO3 B2NO3 CNO3 C2NO3 DNO3 ENO3 FNO3 F2NO3

6.18 8.55 4.77 0.40 0.21 6.25 3.92 0.18 4.42 1.21 0.59

7.24 7.30 4.72 0.41 0.21 6.45 3.84 0.17 3.03 0.95 0.72

6.71 7.93 4.75 0.40 0.21 6.35 3.88 0.17 3.72 1.08 0.66

30.85 26.10 43.58 517.50 985.71 32.60 53.35 1217.65 55.65 191.67 316.03

5.03 5.94 3.56 0.30 0.16 4.76 2.91 0.13 2.79 0.81 0.49

2.14 2.75 1.99 0.16 0.10 1.92 1.27 0.06 1.29 0.37 0.30

TABLE 4: First Solvent Shell of Center of Mass RDFs IL

max position (Å)

min position (Å)

coord no.

ANO3 A2NO3 A3NO3 BNO3 B2NO3 CNO3 C2NO3 DNO3 ENO3 FNO3 F2NO3

5.18 5.18 5.63 4.73 5.00 5.54 5.90 5.00 5.90 4.01/5.81 3.74/5.72

7.97 7.97 7.70 7.52 7.43 8.15 8.15 7.52 7.70 7.70 7.43

5.9 6.3 6.9 6.4 7.4 5.9 6.3 6.1 6.0 5.7 6.9

5. Conclusions

Figure 6. Spatial distribution functions (SDFs) for anion around cation.

TABLE 2: Liquid Densities for Guanidinium-Based ILs IL ANO3 A2NO3 A3NO3 BNO3 B2NO3 CNO3 C2NO3 DNO3 ENO3 FNO3 F2NO3 c

Fexp (g/cm3)a,b

Fcal (g/cm3)a,b

FMD (g/cm3)b

FMD (g/cm3)

1.23

1.20

1.13

1.24 1.34 1.31

1.26 1.32 1.19

1.19 1.27 1.13

1.05c 1.07c 1.14c 1.18d 1.25d 1.05c 1.07c 1.18d 1.15c 1.14e 1.25e

a Calculated and experimental data in the ref 3. b Density at 298 K. Density at 450 K. d Density at 320 K. e Density at 420 K.

show significant differences, although the only difference between the structures of the two occurs at the “c” site. Anions that distribute above and below the ring are effectively removed and replaced with high-probability regions associated, presumably, with the amide hydrogen, which offers a far better point with which to form hydrogen bonds. A similar case was also observed for A3 and B2.

In this work, a force field based on AMBER was proposed for eleven cyclic guanidinium-based ionic liquids. The refinement of the force field was in the frame of AMBER, using a systematic method. Molecular dynamics simulations for these guanidinium-based ILs were performed. The simulation liquid densities are compared with the experimental results and the calculated data from reference. The calculation accuracy could be improved by fine-tuning the parameters for Lennard-Jones parameters. Transport properties including the coefficients of self-diffusion, viscosities, and conductivities were also studied. They are consistent with the fact that the viscosities for most of the ILs are much higher than those water. It is found that the activity reduces with decreasing the length of the alkyl chain for these cyclic guanidinium-based ILs. The microstructures are discussed by using radial distribution functions (RDFs) and spatial distribution functions (SDFs). It can be seen from integrating the RDFs that each cation is surrounded by five to seven anions. The anions are most likely distributed around cations in different regions because of their structural distinction. Acknowledgment. This work was supported by the National Science Fund of China for Distinguished Young Scholars (20625618), National Natural Science Foundation of China (20436050 and 20603040), and the Supercomputing Center, CNIC, CAS. References and Notes (1) Rogers, R. D.; Seddon, K. R. Ionic Liquids as Green Solvents: Progress and Prospects; American Chemical Society: Washington, DC, 2003. (2) Wasserscheid, P.; Welton, T. Ionic Liquids in Synthesis; WileyVCH: Weinheim, Germany, 2003. (3) Gao, Y.; Arritt, S. W.; Twamley, B.; Shreeve, J. M. Inorg. Chem. 2005, 44, 1704. (4) Cole, A. C.; Jensen, J. L.; Ntai, I.; Tran, K. L. T.; Weaver, K. J.; Forbes, D. C.; Davis, J. H., Jr. J. Am. Chem. Soc. 2002, 124, 5962.

5668 J. Phys. Chem. B, Vol. 111, No. 20, 2007 (5) Rajagopal, R.; Jarikote, A.; Dilip, V.; Srinivasan, K. V. Chem. Commun. 2002, 67, 616. (6) Sheldon, R. Chem. Commun. 2001, 66, 2399. (7) Blanchard, L. A.; Hancu, D.; Beckman, E. J.; Brennecke, J. F. Nature 1999, 399, 28. (8) Gao, H.; Han, B.; Li, J.; Jiang, T.; Liu, Z.; Wu, W.; Chang, Y.; Zhang. J. Synth. Commun. 2004, 34, 3083. (9) Wang, P.; Zakeeruddin, S. M.; Gra¨tzel, M.; Kantlehner, W.; Mezger, J.; Stoyanov, E. V.; Scherr, O. Appl. Phys. A 2004, 79, 73. (10) Tobey, S. L.; Anslyn, E. V. J. Am. Chem. Soc. 2003, 125, 14807. (11) Schmuck, C. Chem. Eur. J. 2000, 6, 709. (12) Chavez, D. E.; Hiskey, M. A.; Gilardi, R. D. Org. Lett. 2004, 6, 2889. (13) Neutz, J.; Grosshardt, O.; Schaeufele, S.; Schuppler, H.; Schweikert, W. Propellants Explos. Pyrotech. 2003, 28, 181. (14) Liu, Z.; Huang, S.; Wang, W. J. Phys. Chem. B 2004, 108, 12978. (15) Del Po´polo, M. G.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744. (16) Margulis, C. J.; Stern, H. A.; Berne, B. J. J. Phys. Chem. B 2002, 106, 12017. (17) Morrow, T. I.; Maginn, E. J. J. Phys. Chem. B 2002, 106, 12807. (18) Bhargava, B. L.; Balasubramanian, S. J. Chem. Phys. 2005, 123, 144505. (19) Liu, X.; Zhang, S.; Zhou, G.; Wu, G.; Yuan, X.; Yao, X. J. Phys. Chem. B 2006, 110, 12062. (20) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586.

Liu et al. (21) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (22) Fox, T.; Kollman, P. A. J. Phys. Chem. B 1998, 102, 8070. (23) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993, 115, 9620. (24) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1357. (25) Wang, J. M.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (26) Lyubartsev, A. P.; Laaksonen, A. Comput. Phys. Commun. 2000, 128, 565. (27) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (28) Zhang, S. J.; Sun, N.; He, X. Z.; Lu, X. M.; Zhang, X. P. J. Phys. Chem. Ref. Data 2006, 35, 1475. (29) Fioroni, M.; Burger, K.; Mark, A. E.; Roccatano, D. J. Phys. Chem. B 2000, 104, 12347. (30) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (31) Hertz, H. G.; Franks, F. Water, A comprehensiVe treatise; Plenum Press: New York, 1973; Vol. 3. (32) Yaws, C. L.; Miller, J. W.; Shah, P. N.; Schorr, G. R.; Patel, P. M. Chem. Eng. 1976, 83, 153. (33) Noda, A.; Hayamizu, K.; Watanabe, M. J. Phys. Chem. B 2001, 105, 4603. (34) Keblinshi, P.; Eggebrecht, J.; Woff, D.; Phillpot, S. R. J. Chem. Phys. 2000, 113, 282. (35) de Andrade, J.; Bo¨es, E. S.; Stassen, H. J. Chem. Phys. 2002, 106, 13344. (36) Laaksonen, L. J. Mol. Graphics 1992, 10, 33.

Molecular Simulation of Guanidinium-Based Ionic Liquids - American

Recommend Documents