A Force Field for Molecular Simulation of Tetrabutylphosphonium

Jun 7, 2007 - ... approaches to understanding reaction outcomes of organic processes in ionic liquids. Sinead T. Keaveney , Jason B. Harper , Anna K. ...
0 downloads 0 Views 124KB Size
Download PDF
7078

J. Phys. Chem. B 2007, 111, 7078-7084

A Force Field for Molecular Simulation of Tetrabutylphosphonium Amino Acid Ionic Liquids Guohui Zhou,†,‡ Xiaomin Liu,†,‡ Suojiang Zhang,*,† Guangren Yu,†,‡ and Hongyan He† State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, 100080, Beijing, China, and Graduate UniVersity of Chinese Academy of Sciences, 100049, Beijing, China ReceiVed: December 5, 2006; In Final Form: April 17, 2007

An all-atom force field was set up for a new class of ionic liquids (ILs), tetrabutylphosphonium amino acid, on the basis of the AMBER force field with determining parameters related to the phosphorus atom and modifying several parameters. Ab initio quantum chemical calculations were employed to obtain molecular geometries, infrared frequencies, and torsion energy profiles. Atom partial charges were obtained by using the one-conformation, two-step restraint electrostatic potential approach. Molecular dynamics simulation was carried out in the isothermal-isobaric ensemble for 14 tetrabutylphosphonium amino acid ILs at two temperatures to validate the force field against the experimental densities and heat capacities at constant pressure. Computed thermodynamic properties are in good agreement with available experimental values. Moreover, radial distribution functions were investigated to depict the microscopic structures of these ILs.

1. Introduction Ionic liquids (ILs) are a group of organic salts with melting points below 100 °C.1 ILs show negligible vapor pressures, high thermal stabilities and ionic conductivities, and a wide temperature range in the liquid state.1-6 They have been attracting a great deal of attention from both academic researchers and industrial engineers due to their designable nature for a variety of potential applications. To tailor the properties of ILs and design novel ILs for specific tasks, it is critically important to understand the relationship between microstructures and properties of ILs. Molecular dynamics (MD) simulation is a reliable approach to analyze the relationship between microstructures and macroscopic properties.7-12 ILs with a series of cationic derivatives, such as imidazoliumbased7,8,13,14 and pyridinium-based ILs,10 have been synthesized and simulated due to their convenient chemical modification, but there is much less information about the effect of anion structures.15 A new type of task specific ILs, tetrabutylphosphonium amino acid ([P(C4)4][AA]), was synthesized in our lab recently to capture CO2.16,17 This kind of ILs is composed of a series of anions, which is a good example for understanding the effects of anion structures on the microscopic and macroscopic properties by MD simulation, so we studied them. Such studies are also of vital importance for designing ILs. In this work, a refined force filed based on AMBER18 was set up for 14 [P(C4)4][AA] ILs composed of the tetrabutylphosphonium ([P(C4)4]+) cation and the following amino acid anions, such as glycine ([Gly]-), alanine ([Ala]-), β-alanine ([β-Ala]-), serine ([Ser]-), lysine ([Lys]-), leucine ([Leu]-), isoleucine ([Ile]-), phenylalanine ([Phe]-), proline ([Pro]-), methionine ([Met]-), aspartic acid ([Asp]-), glutamic acid ([Glu]-), glutamine ([Gln]-), and taurine ([Tau]-). Molecular dynamics simulation was carried out for these ILs at two different temperatures. * Address correspondence to this author. E-mail: [email protected]. † Institute of Process Engineering, Chinese Academy of Sciences. ‡ Graduate University of Chinese Academy of Sciences.

Densities and heat capacities at constant pressure were computed to validate the force field. The microscopic structures of the ILs were also analyzed. 2. Force Field Development 2.1. Form of the Potential. The AMBER force field of Cornell et al.18 has been used with the following functional form for the potential energy E:

Etotal )

Kφ (1 + 2 qiqj + (1) rij

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

cos(nφ - γ)) +

N

∑ ∑ i)1 j)i+1

{ [( ) ( ) ] } σij

4ij

12

rij

-

σij rij

6

where the potential parameters have their usual meaning. Lennard-Jones parameters for interactions between different atom types have been obtained by using the Lorentz-Berthelot mixing rule:19

ij ) xiijj σij ) (σii + σjj)/2

(2)

2.2. Atomic Partial Charges. The atomic charges, qi (eq 1), were fitted to reproduce the molecular electrostatic potential, calculated with the B3LYP wave function and the 6-311+G** basis set using the Gaussian 03 program. The fitting was carried out by using the one-conformation, two-stage standard restraint electrostatic potential fitting approach.20 The computed atomic partial charges are presented in Figure S1 (see the Supporting Information). 2.3. Atom Types and Van der Waals Parameters. Atom types were assigned by using the default AMBER parameters.18 Our choice of atom types for the cation and anions is presented in Figure S1 (Supporting Information). The intramolecular and

10.1021/jp068365e CCC: $37.00 © 2007 American Chemical Society Published on Web 06/07/2007

Tetrabutylphosphonium Amino Acid Ionic Liquids van der Waals parameters for these cation and anions have been perfectly defined after this choice of atom types. All of these parameters have been given in Table S1 (see the Supporting Information). 2.4. Bond and Angle Parameters. The bonding and bending parameters for the cation and anions were determined by fitting the IR vibration frequencies. Because the number of peaks in IR frequencies is far less than the quantity of freedom degree (3N - 6), the bond and angle force constants were also adjusted by performing a vibration analysis on quantum mechanics (QM) frequencies calculated from the Gaussian 03 program. However, the HF method overestimates the frequencies by approximately 10%.21 Therefore, a scaling factor of 0.9 is suitable for HF calculations.22 Molecular mechanics (MM) calculations were performed utilizing the TINKER package.23 The MINIMIZE and VIBRATE modules were employed for the normal-mode analysis. Four P-C stretching frequencies for [P(C4)4]+ calculated by QM are 638, 772, 780, and 786 cm-1, and the corresponding frequencies fitted by MM are 640, 755, 759, and 760 cm-1. For anions, take [β-Ala]- as an example, the C-O stretching frequencies of experimental and calculated data by QM are 1639 and 1578 cm-1,16 respectively, and the fitted value is 1616 cm-1. The experimental IR frequencies for N-H stretching are 3277 and 3354 cm-1,16 the corresponding QM frequencies are 3335 and 3406 cm-1, and the fitted ones are 3332 and 3394 cm-1. 2.5. Dihedral Parameters. There are no definitions and parameters for some dihedral angles in AMBER.18 The most important ones are the parameters for the atom P in the center of tetrabutylphosphonium, including P-CT-CT-HC, CT-PCT-HP, CT-CT-P-CT, and CT-CT-CT-P, and the torsions related to the -COO-, -NH2, and -SO3- end groups in amino acid anions. These corresponding parameters were optimized in this work. The ab initio torsion energy profiles were obtained by using the Gaussian 03 program at the MP2/6-311G**//HF/6-31+G* level. In the geometry optimization, the appropriate torsions for each of the torsional angles varied in a step of 10° grid while minimizing the energy with respect to all the other degrees of freedom. The optimized conformers from QM at various dihedral angles were used directly in MM energy calculations by the ANALYZE module in TINKER.23 Finally, the coefficients of the dihedrals were optimized by fitting the QM torsion energy profiles. The optimized parameters are listed in Table S1 (Supporting Information). The energy profiles from QM and MM calculations for CT-CT-P-CT and CT-CT-CT-P are shown in Figure 1.

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7079

Figure 1. Torsion profiles of the dihedrals in [P(C4)4]+, calculated by QM (lines) and MM (dots) in this work.

The equilibration of the simulations was extended to 400 ps and each production phase lasted for another 1.0 ns at 321.85 and 298.15 K under 1.0 atm. Simulations were carried out on Deepcomp 6800 running Linux. A 2964 ps simulation of 192 [P(C4)4][Ser] molecules (12672 atoms) took about 614 h (25.6 days) to complete on 4 nodes. Each node contains four Intel Itanium II processors operating at 1.3 GHz. 4. Results and Discussion

3. Simulation Details MD simulations for all 14 [P(C4)4][AA] ILs were performed with the M.DynaMix24 program version 4.3. Each system contained 192 tetrabutylphosphonium cations and 192 amino acid anions. The simulation started from an FCC lattice at a low density of 0.2 g/cm3. After the system was relaxed in the NVE ensemble for a few MD steps (2000 steps) to remove the possible overlapping in the initial configuration, the NoseHoover25 NpT ensemble was adopted with coupling constants of 700 and 100 fs, and the system density gradually increased to a reasonable value. The Tuckerman-Berne double time step algorithm26 was employed with long and short time steps of 2 and 0.2 fs, respectively. The intramolecule forces were cut off at 14 Å, while the long-range forces including LJ and Coulombic interactions were cut off at 18 Å, and the Ewald summation27 was implemented for the latter.

4.1. Liquid Densities. Liquid density is one of the most important sources that can be used to validate a proposed force field.8,28 In this work, experimental densities16,17 from the literature are compared with our calculated results. As is shown in Table 1, the predicted densities from our proposed force field are in good agreement with the experimental data. 4.2 Heat Capacities at Constant Pressure. The heat capacities at constant pressure were calculated with29

Cp(T,p) )

( ) [

]

H(T + ∆T) - H(T) ∂H ≈ ∂T p ∆T

p

(3)

where H is the enthalpy.

H ) ENB + EINT + K + pV

(4)

7080 J. Phys. Chem. B, Vol. 111, No. 25, 2007

Zhou et al.

Figure 2. Center of mass radial distribution functions of [P(C4)4][AA]: solid lines, +/-; dashed lines, +/+; dash-dotted lines, -/-.

Tetrabutylphosphonium Amino Acid Ionic Liquids

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7081

Figure 2. Continued.

The total potential energy E is split into nonbonded and intramolecular terms (E ) ENB + EINT), where ENB includes all Lennard-Jones and Coulombic interactions and EINT includes bond length, bond angle, dihedral angle, and improper angle terms. K is the kinetic energy, p is the pressure, and V is the molar volume. The heat capacities at constant pressure were computed from the simulations carried out at 1.0 atm and at varying temperatures (298.15 and 321.85 K). The computed heat capacities at 310 K are listed in Table 1. They are in reasonable agreement with the experimental values.17 4.3. Microstructures. To understand the structure of [P(C4)4][AA], the center-of-mass RDFs for the cation-cation, anionanion, and cation-anion were computed. The RDFs of the same IL at 298.15 and 321.85 K are qualitatively similar. The RDFs for each IL at 298.15 K are shown in Figure 3. The positions of the first minimum in cation-anion RDFs are quite different. For the anions like [Gly]- and [Ala]-, it locates at about 5 Å,

while for the anions just like [Glu]-, [Phe]-, and [Gln]-, it locates at about 6 Å. The differences are most likely caused by the anion volumes. The locations of the maximum and minimum in the cation-anion RDFs of [P(C4)4][AA] are listed in Table 2. By integrating g(r) out to the location of the first minimum, the coordination number for the first solvation shell, N, can be calculated via the following equation

N ) 4π

∫0R

min1

Fndg(r)r2 dr

(5)

where Fnd is the number density and Rmin1 refers to the first minimum in g(r). As is seen from Table 2, there are about three to four counterions surrounding each ion. A more detailed liquid structure can be described by the sitesite pair RDFs. Figure 3 shows the RDFs of the oxygen atom (O2) in -COO- or -SO3- of five typical anions and the hydrogen atoms (HP) in the tetrabutylphosphonium. The posi-

7082 J. Phys. Chem. B, Vol. 111, No. 25, 2007

Zhou et al.

TABLE 1: Simulated and Experimental Densities (G) and Heat Capacities at Constant Pressure (Cp) of [P(C4)4][AA]a ILs [P(C4)4][Gly] [P(C4)4][Ala] [P(C4)4][β-Ala] [P(C4)4][Leu] [P(C4)4][Ile] [P(C4)4][Ser] [P(C4)4][Lys] [P(C4)4][Asp] [P(C4)4][Glu] [P(C4)4][Gln] [P(C4)4][Pro] [P(C4)4][Phe] [P(C4)4][Met] [P(C4)4][Tau] a

Fsim (g cm-3) 321.85 K 298.15 K 0.9444 0.9295 0.9375 0.9114 0.9133 0.9784 0.9505 1.0134 1.0073 0.9915 0.9484 0.9674 0.9684 1.0214

0.9545 0.9425 0.9545 0.9243 0.9273 0.9934 0.9624 1.0284 1.0183 1.0036 0.9634 0.9803 0.9844 1.0374

Fexptl (g cm-3) 298.15 K

∆% 298.15 K

Cp(sim) (J g-1 K-1) 310 K

Cp(exptl) (J g-1 K-1) 310 K

0.9630 0.9500 0.9590 0.9269 0.9296 0.9910 0.9730 1.0173 1.0121 1.0519 0.9828 0.9524 0.9868 1.0199

-0.9 -0.8 -0.5 -0.3 -0.3 0.2 -1.1 1.1 0.6 -4.6 -2.0 2.9 -0.3 1.7

2.52506 2.58707 2.66508 2.56406 2.51507 2.46806 2.63405 2.45209 2.40510 2.41312 2.49310 2.41506 2.47605 2.46605

2.1422 2.4582 2.4577 2.3591 2.3186 2.3599 1.8434 / 2.2223 2.1748 2.2336 2.6112 2.4424 2.2403

All the experimental values are taken from refs 16 and 17.

TABLE 3: First Maximum and Minimum Positions of HP-O2 Radial Distribution Functions and the First Shell Coordination Numbers

Figure 3. Radial distribution functions between the HP atom of tetrabutylphosphonium and the O2 atom in five typical amino acid anions.

TABLE 2: First Maximum and Minimum Positions of Center of Mass Radial Distribution Functions for the Cation-Anion of [P(C4)4][AA] and the First Shell Coordination Numbers ILs

max position (Å)

min. position(Å)

peak height

coord no.

[P(C4)4][Gly] [P(C4)4][Ala] [P(C4)4][β-Ala] [P(C4)4][Leu] [P(C4)4][Ile] [P(C4)4][Ser] [P(C4)4][Lys] [P(C4)4][Asp] [P(C4)4][Glu] [P(C4)4][Gln] [P(C4)4][Pro] [P(C4)4][Phe] [P(C4)4][Met] [P(C4)4][Tau]

4.85 5.05 5.15 5.65 5.55 5.45 5.85 5.65 6.15 5.85 5.45 6.25 6.05 5.25

7.35 7.55 7.45 8.05 7.85 7.55 8.05 8.15 8.25 8.05 7.45 8.15 7.95 7.65

4.58 4.04 3.84 3.79 4.01 3.27 2.28 3.05 2.23 2.66 3.66 2.35 2.37 4.41

3.26 3.32 3.30 3.51 3.37 3.59 3.55 3.82 4.14 3.91 3.18 3.63 3.69 3.43

tions of the first peak and the first minimum for all the studied [P(C4)4][AA] are very similar. The first peaks and minimum occur at about 2.3 and 3.45 Å, respectively. Furthermore, [P(C4)4][Gly] and [P(C4)4][Ala] were studied by ab initio quantum chemical calculations,30 and similar structures are observed. The locations of the maximum and minimum in the HP-O2 RDFs are listed in Table 3. The heights of the first peak

ILs

max position (Å)

min. position (Å)

peak height

coord no.

[P(C4)4][Gly] [P(C4)4][Ala] [P(C4)4][β-Ala] [P(C4)4][Leu] [P(C4)4][Ile] [P(C4)4][Ser] [P(C4)4][Lys] [P(C4)4][Asp] [P(C4)4][Glu] [P(C4)4][Gln] [P(C4)4][Pro] [P(C4)4][Phe] [P(C4)4][Met] [P(C4)4][Tau]

2.35 2.35 2.35 2.25 2.35 2.35 2.35 2.35 2.35 2.35 2.35 2.35 2.35 2.35

3.45 3.45 3.45 3.45 3.35 3.45 3.45 3.45 3.45 3.45 3.45 3.45 3.45 3.45

4.22 4.26 4.30 4.68 4.71 3.43 4.53 3.61 3.69 3.67 4.52 4.73 4.59 3.19

0.50 0.48 0.50 0.46 0.42 0.46 0.45 0.42 0.41 0.40 0.48 0.45 0.47 0.40

are different: the peak heights are about 3.6 for [Ser]-, [Asp]-, [Gln]-, and [Glu]-, and the peak heights of others are greater than 4.2. To investigate the effects of the functional group of anions on the transport properties, the RDFs related to four typical groups are computed. The RDFs of O2-HX (HX is the hydrogen atom in the hydroxyl group of [Ser]-, the azyl group of [Lys]that is not bonded to the R-carbon atom, the carboxyl group of [Glu]-, and the acylamino group of [Gln]-) are shown in Figure 4a. The first peak height of carboxyl is 18.93, which is obvious higher than those of other groups. The first peaks occur at 1.65, 1.75, 1.85, and 1.95 Å for carboxyl, hydroxyl, acylamino, and azyl groups, respectively. It is revealed that the interaction between the O2 and the HX atom decreases in the order carboxyl > hydroxyl > acylamino > azyl. Figure 4b shows the RDFs of the O-HX (HX is the hydrogen atom in the carboxyl group of [Glu]- and the acylamino group of [Gln]-). Comparing with Figure 4a, it is found that the H atom in the acylamino group of [Gln]- prefers to distribute around O rather than O2, which is contrary to that of the carboxyl group of [Glu]-. The numbers of hydrogen bonds for O2-HX and O-HX were calculated. The criteria used to define a hydrogen bond were hydrogen-to-acceptor distance e2.4 Å and acceptor-hydrogendonor angle >135°.31 The hydrogen bond numbers for O2-HX and O-HX of the four anions are listed in Table 5. Figure 5 shows the experimental conductivities and viscosities as a function of the hydrogen bond numbers. It is obvious that the

Tetrabutylphosphonium Amino Acid Ionic Liquids

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7083

Figure 4. Radial distribution functions related to four typical groups in amino acid anions.

TABLE 4: First Maximum and Minimum Positions of O2-HX and O-HX Radial Distribution Functions and the First Shell Coordination Numbers anions

max position (Å)

[Ser][Lys][Glu][Gln]-

1.75 1.95 1.65 1.85

[Glu][Gln]-

1.85 1.95

min. position (Å)

peak height

coord no.

HO/H-O2 2.55 2.65 2.45 2.75

8.05 3.05 18.93 3.75

0.19 0.13 0.34 0.11

HO/H-O 2.85 2.85

1.58 4.95

0.05 0.18

TABLE 5: The Experimental Viscosities and Conductivities and Hydrogen Bonds Numbers of O2-HX and O-HXa anions

no. of hydrogen bonds

viscosities (mPa s)

conductivities (10-4 S cm-1)

[Ser][Lys][Glu][Gln]-

68 74 125 135

734.20 744.71 9499.68 9561.26

0.87 0.83 0.063 0.048

a

All the experimental values are taken from refs 16 and 17.

5. Conclusions An all-atom force field was set up for [P(C4)4][AA] ILs with determining parameters related to the atom P and modifications on several parameters. MD simulations of 14 [P(C4)4][AA] ILs were carried at two different temperatures. The densities and heat capacities at constant pressure were predicted. The agreement between the experimental densities and computed values is excellent, and the agreement for heat capacities at constant pressure is reasonable. The RDFs are employed to describe the local structures of these ILs. The RDFs of HP-O2 reveal that HP in the tetrabutylphosphonium can form hydrogen bonds with O2 in different anions. Hydrogen in azyl, hydroxyl, acylamino, and carboxyl groups of anions can form hydrogen bonds with O2 or O in different anions. The existence of hydrogen bonds explains to some extent the low electrical conductivities and high viscosities of the studied [P(C4)4][AA] ILs. Acknowledgment. This work was supported by the National Science Fund of China for Distinguished Young Scholars (20625618), the National Natural Science Foundation of China (20436050, 20603040), and the Supercomputing Center, CNIC, CAS. Supporting Information Available: Listing of proposed force field parameters and a figue showing atom types and partial charges of the tetrabutylphosphonium cation and 14 amino acid anions. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes

Figure 5. The experimental conductivities (square) and viscosities (point) as a function of hydrogen bond numbers.

conductivities decrease with hydrogen bond number, contrary to the viscosities. The relatively low electrical conductivities and high viscosities of these kinds of ILs may be due to the hydrogen bonds.

(1) Zhang, S. J.; Lu, X. M. Ionic Liquids: Form Fundamentals to Applications; Science Press: Beijing, 2006. (2) Blanchard, L. A.; Hancu, D.; Beckman, E. J.; Brennecke, J. F. Nature 1999, 399, 28. (3) Sheldon, R. Chem. Commun. 2001, 66, 2399. (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. (5) Rajagopal, R.; Jarikote, A.; Dilip, V.; Srinivasan, K. V. Chem. Commun. 2002, 67, 616. (6) Gao, H. X.; Han, B. X.; Li, J. C.; Jiang, T.; Liu, Z. M.; Wu, W. Z.; Chang, Y. H.; Zhang, J. M. Synth. Commun. 2004, 34, 3083. (7) Morrow, T. I.; Maginn, E. J. J. Phys. Chem. B 2002, 106, 12807. (8) Liu, Z. P.; Huang, S. P.; Wang, W. C. J. Phys. Chem. B 2004, 108, 12978. (9) Lopes, J. N. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2004, 108, 16893. (10) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. J. Phys. Chem. B 2006, 110, 2821. (11) Lopes, J. N. A. C.; Pa´dua, A. A. H. J. Phys. Chem. B 2006, 110, 3330.

7084 J. Phys. Chem. B, Vol. 111, No. 25, 2007 (12) Liu, X. M.; Zhang, S. J.; Zhou, G. H.; Wu, G. W.; Yuan, X. L.; Yao, X. Q. J. Phys. Chem. B 2006, 110, 12062. (13) Po´polo, M. G. D.; Voth, G. A. J. Phys. Chem. B 2004, 108, 1744. (14) Po´polo, M. G. D.; Voth, G. A. J. Phys. Chem. B 2006, 110, 14426. (15) Fukumoto, K.; Yoshizawa, M.; Ohno, H. J. Am. Chem. Soc. 2005, 127, 2398. (16) Zhang, J. M.; Zhang, S. J.; Dong, K.; Zhang, Y. Q.; Shen, Y.; Lu, X. Chem. Eur. J 2006, 12, 4021. (17) Zhang, J. M.; Zhang, S. J. Combinatorial chemistry reserch on amino acids and ionic liquids. Postdoctoral Thesis, Institute of Process Engineering, Chinese Academy of Sciences, 2005. (18) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (19) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, UK, 1987. (20) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (21) Lii, J. H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8566. (22) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-

Zhou et al. McCarthy, 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. (23) Pappu, R. V.; Hart, R. K.; Ponder, J. W. J. Phys. Chem. B 1998, 102, 9725. (24) Lyubartsev, A. P.; Laaksonen, A. Comput. Phys. Commun. 2000, 128, 565. (25) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (26) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (27) Deleeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, Ser. A 1983, 388, 177. (28) Zhang, S. J.; Sun, N.; He, X. Z.; Lu, X. M.; Zhang, X. P. J. Phys. Chem. Ref. Data 2006, 35, 1475. (29) Lagache, M.; Ungerer, P.; Boutina, A.; Fuchsa, A. H. Phys. Chem. Chem. Phys. 2001, 3, 4333. (30) Yu, G. R.; Zhang, S. J.; Yao, X. Q.; Zhang, J.; Dong, K.; Dai, W. B.; Mori, R. Ind. Eng. Chem. Res. 2006, 45, 2875. (31) Mark, P.; Nilsson, L. J. Phys. Chem. B 2002, 106, 9440.