Development of a ReaxFF Reactive Force Field for ... - ACS Publications

Article pubs.acs.org/JPCB

Development of a ReaxFF Reactive Force Field for Tetrabutylphosphonium Glycinate/CO2 Mixtures Bo Zhang,†,‡ Adri C. T. van Duin,†,§ and J Karl Johnson*,†,‡ †

United States Department of Energy, National Energy Technology Laboratory, Pittsburgh, Pennsylvania 15236, United States Department of Chemical & Petroleum Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15261, United States § Department of Mechanical and Nuclear Engineering, Pennsylvania State University, University Park, Pennsylvania 16802, United States ‡

S Supporting Information *

ABSTRACT: Carbon dioxide interacts with the ionic liquid tetrabutylphosphonium glycinate, [P(C4)4][Gly], through both physical and chemical absorption. We present a parametrization of the ReaxFF force field for this system that accounts for both chemical and physical interactions. The parametrization was developed from an extensive training set including periodic density functional theory (DFT) calculations of reaction pathways between CO2 and the anion [Gly]− in the condensed phase, condensed-phase molecular dynamics (MD) configurations, gas-phase CO2−anion and CO2−cation interactions, and gas-phase cluster calculations for intra-ion interactions. The optimized ReaxFF parameters capture the essential features of both physical and chemical interactions between CO2 and [P(C4)4][Gly] as compared with experiments, van der Waals-corrected DFT calculations, or, in the case of physical interactions, classical force field calculations. The probability distributions of the distance between C (from CO2) and N (from the anion) and the CO2 bend angles calculated from MD simulations with the optimized ReaxFF force field are in good general agreement with those from DFT-based MD simulations. We predict that the density of CO2/[P(C4)4][Gly] mixtures increases with increasing CO2 concentration up to at least 50 mol % CO2. We attribute the significant increase in density to the small effective volume occupied by chemically bound CO2 in the mixture. The predicted increase in density may be tested experimentally.

1. INTRODUCTION The increase of CO2 concentration in the atmosphere due to anthropogenic activity has become a major concern due to its link with global climate change.1,2 While renewable energy sources and nuclear power offer the potential to reduce carbon emissions over the long-term, fossil fuels will very likely remain the major source of energy for several decades. Therefore, CO2 capture and sequestration is one of the most promising strategies to mitigate CO2 emissions in the near future.3 Chemisorption (scrubbing) of CO2 with aqueous monoethanolamine (MEA) is a commercially available technique that is able to capture 90% of the CO2 from flue gas.3,4 However, this approach suffers from thermal and chemical degradation of MEA5,6 and high energy costs needed to regenerate the solvent.7 Despite research and development on many different materials, there are currently no CO2 capture technologies available that meet all of the DOE targets for cost and efficiency.3 Therefore, significant effort is being devoted to the development of new sorbents, either in solid form, such as oxides,8−11 hydroxides,9−12 carbonates,10,13 and metal organic frameworks,14−16 or in liquid form, such as amine-based solvents7,17,18 or ionic liquids (ILs).19−21 As demonstrated by many studies from experiments and theoretical calculations, the © 2014 American Chemical Society

interaction strength between CO2 and ILs can be tailored by using different combinations of cations and anions. The interaction can be tuned to range from physical interactions (∼10−20 kJ/mol),22 to complexation interactions (∼30−40 kJ/mol),23 to chemical interactions (∼60−80 kJ/mol).24 This ability to tune the interaction strength between CO2 and the sorbent makes ILs promising candidates for CO2 sorbents suitable for specific capture conditions. However, there are a tremendous number of different cations and anions that can be chosen to make an IL, making it impractical to synthesize even a small fraction of all possible ILs and test their CO2 capture capability. The process of identifying promising ILs for CO2 capture can be accelerated by screening ILs through calculations of thermodynamic and transport properties from molecular simulations. Most classical force fields employ simple empirical potentials to account for bonded interactions, such as bond lengths, angles, and torsions, and nonbonded interactions, such as van der Waals (vdW) and Coulomb forces. This methodology of using simple functional forms makes the Received: June 2, 2014 Revised: September 22, 2014 Published: October 6, 2014 12008

dx.doi.org/10.1021/jp5054277 | J. Phys. Chem. B 2014, 118, 12008−12016

The Journal of Physical Chemistry B

Article

devoted to force field optimization. Finally, the optimized force field is validated for the CO2/[P(C4)4][Gly] system by comparing various properties, such as density, heat capacity, and reaction barriers, with density functional theory (DFT) calculations.

molecular simulations with classical force fields computationally very efficient. There are a number of different classical force fields for different ILs that exhibit physical interactions with CO2.19,25 However, these potentials cannot capture the complicated interactions between ILs and CO2 resulting from charge-transfer complexation or chemical reactions. Although simulations of reactive ILs with CO2 have been performed by ad hoc mixing of different fractions of unreacted ILs and reacted ILs,26−28 this approach lacks information on how reactions between ILs and CO2 take place and how the fraction of reacted CO2 changes with the process conditions and therefore lacks real predictive power. Thus, there is a need for a force field that will accurately describe both physical and chemical interactions between CO2 and the IL. The ReaxFF formalism proposed by van Duin et al. is one of the force fields that can provide such functionality.29,30 ReaxFF is bond-orderdependent and thus naturally allows for bond breaking/forming during simulation. ReaxFF has been successfully applied to many different chemical systems, such as crystals and surfaces,31 glycine tautomerization,32 and hydrocarbon oxidations,30 to name a few. In this work, tetrabutylphosphonium glycinate, [P(C4)4][Gly], is chosen as a model IL because it has both −CO−2 and − NH2 functional groups that exhibit complex formation and reactive interactions with CO2, respectively. Thus, physical, complexation, and chemical interactions will be expected to occur simultaneously in this IL/CO2 system. The structure of [P(C4)4][Gly] is shown in Figure 1. Zhang et al. have prepared

2. COMPUTATIONAL METHODS 2.1. ReaxFF Reactive Force Field Formalism. ReaxFF is a reactive force field that can describe bond-breaking and bondformation events.29,30 Unlike traditional force fields, each element is described by only one single atom type, even within different chemical environments, in the ReaxFF formalism. The reaction site or connectivity information is not needed beforehand. Instead, this information is derived from bond orders calculated from interatomic distances at every molecular dynamics (MD) step. Similar to traditional force fields, the vdW and Coulomb interactions are also accounted for by pair interactions from all atoms, irrespective of connectivity. A shielding term is also included to avoid excessively close range nonbonded interactions. Polarization effects are considered through a geometry-dependent charge equilibration scheme.35,36 A detailed description of the ReaxFF potential function can be found elsewhere.29,30 The force field parameters were determined by starting from a previous glycine force field with new training sets relevant to [P(C4)4][Gly]/CO2 mixtures. The single-parameter search optimization method37 was used to minimize the following sum of squares n

Error =

⎛ xi ,DFT − xi ,ReaxFF ⎞ ⎟ σi ⎝ ⎠

∑⎜ i=1

(1)

where xi,DFT is the energy value calculated from DFT for some configuration i, xi,ReaxFF is the value computed from ReaxFF, and σi is the weighting factor for that configuration. 2.2. Training Set Generation. An extensive set of DFT calculations was performed in order to obtain accurate transition states and reaction energies for processes involving CO2 reacting with [P(C4)4][Gly] in the condensed phase. The Vienna ab initio simulation package (VASP),38−40 a periodic plane wave DFT code, was used for most calculations. Core− electron interactions were described by projector augmented wave (PAW) potentials.41,42 The generalized gradient approximation exchange−correlation functional of Perdew, Burke, and Ernzerhof (PBE) was employed.43,44 It is well-known that DFT does not typically describe the long-range electron correlation that represents the dispersion interactions due to the local or semilocal character of most functionals.45 Various methods have been proposed to overcome this shortcoming of DFT.46−49 DFT-D2 is one of the simplest methods; it introduces an additional dispersion energy term to the Kohn−Sham energy to empirically account for the dispersion interactions,48 and we use this approach here due to its computational efficiency. A detailed description of the DFT-D2 method can be found elsewhere.48 A plane wave cutoff energy of 520 eV was used for all calculations. The Brillouin zone was sampled with only the Gamma point due to the large system size. Geometry relaxation of atomic positions with the conjugate gradient algorithm was used to obtain the local minima of the systems. The stopping criteria used were energy differences less than 0.1 meV and forces less than 0.01 eV/Å. The minimum-energy pathway and reaction barriers for reactions between [P(C4)4][Gly] and CO2 were obtained from

Figure 1. Structure of [P(C4)4][Gly].

several new ILs based on tetrabutylphosphonium amino acid, [P(C4)4][AA], where AA = glycine, L-alanine, L-β-alanine, Lserine, or L-lysine.21 Due to high viscosity, these ILs were impregnated into silica gel supports in order to decrease masstransfer limitations for CO2 capture. The supported ILs show faster CO2 absorption than the bulk ILs. They proposed that a CO2 molecule attacks the nitrogen site of the amine group and forms an −NHCO−2 group, as measured by IR and 13C NMR. They also proposed that the leaving proton will attach to the newly formed −NHCO−2 group or the original −CO−2 group in the anion to form a new CO2H group, which will form a hydrogen bond O−H···N with the electron pair of NH2 in another IL.21 In contrast, Bates et al. proposed that a primary amine-functionalized IL will form a carbmate salt, as the result of a proton leaving the amine group to which the CO2 is bound and attaching to the amine on a different ion.20 Zhou et al. have developed an AMBER force field for [P(C4)4][Gly] and validated it by comparing with experimental density and heat capacities of the pure IL.33 Kowsari et al. used the same force field to compute the transport properties of [P(C4)4][Gly].34 To the best of our knowledge, simulations of chemical reactions in [P(C4)4][Gly]/CO2 mixtures have not been reported. We present the development and validation of ReaxFF for a mixture of [P(C4)4][Gly] and CO2. The reaction mechanism and its energetics are discussed in section 2. Section 3 is 12009



Article

Figure 2. Structures of [P(C4)4][Gly] with CO2 at various stages of reaction. (a) Initial reactant local minimum structure of [P(C4)4][Gly] and CO2. (b) Final product structure of [P(C4)4][Gly] and CO2. (c) Transition state 1 and (d) transition state 2 of the reaction between [P(C4)4][Gly] and CO2 as computed from DFT NEB methods. The imaginary vibrational modes are shown as black arrows in (c) and (d).

geometry parameters were constrained according to the specified reaction coordinates. For example, for the P−C bond dissociation energy, the key geometry parameter was the distance between the P atom and C atom. Thus, this distance was constrained while other geometry variables were relaxed to their local minima in the ReaxFF calculations. 2.3. MD Simulation Details. MD simulations were performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) MD package.55 The isothermal− isobaric (NpT) ensemble was employed for equilibration and to measure the density. A time step of 0.5 fs was used, and the temperature was maintained using a Nosé−Hoover thermostat with a damping parameter of 50 fs. Pressure was maintained at 1 atm for all NpT simulations using by Nosé−Hoover barostat with a damping parameter of 500 fs. The temperature was first ramped up to 350 K for 1 ns to quickly equilibrate the system and relax high-energy configurations. Then, the temperature was lowered to 300 K to carry out production runs of 1 ns to determine the equilibrium density. We also used constant volume NVT ensemble simulations, again with the Nosé− Hoover thermostat with a damping parameter of 50 fs. We also employed constant energy NVE simulations for some calculations to avoid any influence of the thermostat. The time step was lowered to 0.25 fs for NVE simulations to ensure good energy conservation.

the climbing-image nudged elastic band (NEB) method as implemented within VASP.50−53 The stopping criterion for NEB calculations was forces less than 0.01 eV/Å. The transition states were all verified by frequency analysis, where only one significant imaginary frequency was found. We developed a training set specifically for describing the intra-ion interactions of the [P(C4)4]+ cation by carrying out gas-phase cluster calculations with Gaussian 09.54 The calculations were carried out at the PBEPBE/aug-cc-pVTZ or MP2/aug-cc-pVTZ levels of theory to ensure that the quality of the calculations was as good as or better than our VASP calculations. The vdW energy contribution was also added post hoc for all PBE functional calculations with DFT-D2 via Grimme’s code.49 P−C bond dissociation energies, C−P−C angle distortion energies, and C−C−P−C rotational barriers were calculated by fixing the bond/angle of interest at a specified value and optimizing all of the other coordinates. We compared the P−C bond dissociation energy of [P(C4)4][Gly] with VASP and [P(C2)4]+ with Gaussian 09 and found no significant difference, giving us confidence in using data from these two different computational approaches to construct our training set. The C−O bond dissociation and CO2 bend angles for the CO2 molecule and the N−C bond dissociation of CO2 interacting with the gas-phase [Gly]− anion were also calculated within Gaussian 09 and included in the training sets. In the parametrization process, the structures in the training set were relaxed with the ReaxFF force field, while the key 12010



Article

3. RESULTS AND DISCUSSION 3.1. Reaction Pathway Analysis of [P(C4)4][Gly] with CO2. The reaction pathway configurations for CO2 binding to the [Gly]− anion are shown in Figure 2. We start the calculations for a reactant state (see Figure 2a) that already shows a strong interaction between CO2 and the N atom on [Gly]− but for which we could find no barrier. The distance between the carbon atom in CO2 and the N atom in [Gly]− is 1.62 Å, while the CO2 bend angle (O−C−O) is about 137.4°. This implies that CO2 interacts with [Gly]− very strongly. However, this structure still has an energy 10 kcal/mol higher than the product state shown in Figure 2b, as computed from DFT. In the product state, the C−N distance is reduced to 1.38 Å, while the CO2 bend angle is about 123°. In addition, the proton previously attached to nitrogen in the anion is now shared between the oxygen from CO2 and the carboxylic group from [Gly]− and with preference for CO2 over the carboxylate, as indicated by the shorter bond length between the O in CO2 and the H (about 1.11 Å) compared with the carboxylate O−H distance of 1.35 Å. These all suggest stronger interaction between CO2 and [Gly]− in the product state than that in the reactant state; the ring structure appears to further stabilize the structure. On the basis of the reactant and product structures, a linear interpolating scheme was used to construct the pathway used in NEB calculations. The final minimum-energy pathway from NEB calculations is shown in Figure 3. We found two

frequency analysis of transition state 1, and the associated vibrational mode vector is also indicated in Figure 2c, which is the vibrational mode of the proton pointing to the nitrogen atom or the oxygen atom. The barrier height is about 0.7 kcal/ mol. After the proton is transferred to the oxygen atom from the carboxylic group in the anion, the next step in the reaction involves the proton approaching the oxygen atom in CO2. The second barrier can be observed during this process in the minimum-energy pathway, with a barrier height of about 2 kcal/mol, as shown in Figure 3. The corresponding structure of transition state 2 is shown in Figure 2d, which shows that the proton is oscillating between the oxygen atom from the carboxylic group and the oxygen atom from CO2, based on vibrational mode analysis. The transition state is also verified by frequency analysis, and only one significant imaginary frequency is found. That mode corresponds to the librational mode and is shown in Figure 2d. A ReaxFF training set was constructed from all of the VASP and Gaussian 09 calculations, along with some of configurations from the previous training set for glycine.32 3.2. ReaxFF Force Field Parametrization. The ReaxFF force field was parametrized against the training set. The optimized parameters are given in the Supporting Information. Figure 3 shows the relative energies calculated from DFT and the ReaxFF force field for the structures taken from the reaction pathway. Generally, the trend from the ReaxFF energies agrees with the DFT energies, although some discrepancies do exist. For example, the first barrier is overestimated by ReaxFF by about 1 kcal/mol, while the second barrier is slightly overestimated by about 0.2 kcal/mol. We also parametrized the ReaxFF force field to better describe the [P(C4)4]+ cation. The relative P−C bond dissociation energy for [P(C4)4][Gly] computed from VASP within periodic boundary conditions is almost identical to the P−C bond dissociation energy computed for the [P(C2)4]+ cation within Gaussian 09, as can be seen from Figure 4. The

Figure 3. Comparison of the DFT and ReaxFF relative energies of the configurations from the reaction pathway of condensed-phase [P(C4)4][Gly] and CO2 obtained from NEB calculation. Blue circles: DFT energy from the VASP NEB calculation. Red squares: ReaxFF energy. The reactant (at 0 on the abscissa) geometry is given in Figure 2a. The first transition state, shown in Figure 2c, appears at the first maximum in the curves (at about 0.7 Å). The second transition state, Figure 2d, occurs at about 3.2 Å.

energy barriers along the minimum-energy pathway. We denote the structure corresponding to the first barrier (from left to right) as transition state 1 (shown in Figure 2c) and denote the structure corresponding to the second barrier (from left to right) as transition state 2 (shown in Figure 2d). We can see that as the CO2 approaches the amine group in [Gly]−, the proton bonded to the nitrogen atom gradually leaves it and approaches one of the oxygens in the carboxylic group. This results in the formation of transition state 1, where the proton is about midway between the nitrogen atom and the oxygen atom in the anion. There is only one significant imaginary frequency (three other modes are very close to zero, and these are assumed to be associated with translational motion) from

Figure 4. Comparison of the DFT and ReaxFF relative energies of P− C bond dissociation in [P(C4)4]+. Blue circles: DFT energies from VASP for [P(C4)4][Gly]. Red squares: ReaxFF energies for [P(C4)4][Gly]. Green triangles: DFT energies from Gaussian 09 for the [P(C2)4]+ cation.

P−C bond stretching energies from ReaxFF agree fairly well with DFT energies from 1.4 to 1.8 Å. The energies computed from ReaxFF are too large from about 1.9 to 2.5 Å and are too small beyond about 2.6 Å. Because there is no bond breaking for the cation when [P(C4)4][Gly] reacts with CO2, the discrepancy between ReaxFF and DFT energies beyond 2.6 Å will not be a major concern. 12011



Article

The C−P−C angle distortion and C−P−C−C torsion angle distortion relative energies were calculated using the isolated cation [P(C2)4]+ in Gaussian 09. Figure 5 shows the

Figure 6. Comparison of the DFT and ReaxFF relative energies of (a) the C−O bond dissociation and (b) the bending angle for the CO2 molecule. Blue circles: DFT energies. Red squares: ReaxFF energies. The DFT calculations were performed with one CO2 molecule in Gaussian 09. Figure 5. Comparison of the DFT and ReaxFF relative energies for (a) the C−P−C angle distortion in [P(C2)4]+ and (b) the C−P−C−C torsion angle distortion in [P(C2)4]+. Blue circles: DFT energies computed from Gaussian 09. Red squares: ReaxFF energies.

from DFT and ReaxFF are qualitatively similar. The energy as a function of the CO2 bend angle is shown in Figure 6b. ReaxFF reproduces the general trend but overestimates the energy of the CO2 angle distortion from 135 to 160° by at most 10 kcal/ mol while reproducing the energies at higher fidelity for all of the other angles. This error in the bending energy is not as significant as it might appear because this bend energy shown in Figure 6b is for the gas-phase or unreacted condensed-phase CO2 molecule, where it is extremely unlikely that one would observe bend angles < 160°, and not for the species that is bound to the [Gly]− anion, as seen in Figure 2, where large bend angles are observed. The ReaxFF description of the CO2− [Gly]− compound is fundamentally different from that of the gas-phase CO2 because the former is identified by ReaxFF as a very different chemical species due to the changes in bond order between free and bound CO2. The bend angle for the CO2−[Gly]− complex is actually reproduced accurately, as seen by comparing VASP and ReaxFF MD simulations of CO2/ [P(C4)4][Gly] (vide infra). 3.3. Validation and Predictions. The training set included data to give reasonable bulk densities for bulk [P(C4)4][Gly]. We included calculations of the energy as a function of the density, where we artificially manipulated the density by increasing the distance between the center of mass of the ions, holding the intra-ion bond angles and bond lengths fixed. Single-point energies were computed from both VASP DFTD2 and the published classical force field33 for [P(C4)4][Gly]. We used the fitted ReaxFF parameters to predict the density of CO2/IL mixtures as a function of CO2 concentration. Densities

comparison of the relative energies calculated from Gaussian 09 and ReaxFF. Overall, there is generally good agreement for the C−P−C angle distortion and C−P−C−C rotational barrier between DFT and ReaxFF calculations. The C−P−C equilibrium angle is adequately predicted by ReaxFF, although it overestimates the energies for angles significantly larger than the equilibrium value. The C−P−C−C torsion barriers are overestimated from ReaxFF by only about 0.2 kcal/mol, but the equilibrium angles are somewhat off, and the shape of the curve is qualitatively different. As we noted above, the bend angle for a CO2 molecule changes depending on its chemical interactions with the anion. For instance, when there is no reaction between CO2 and the IL, the CO2 bend angle should be around 180°. As the product state is formed, this angle will gradually close to about 120°. Therefore, it is necessary to include the C−O bond dissociation and CO2 bending energies in the training set and optimize the force field to reproduce those energetics from DFT calculations. The optimized ReaxFF force field can reproduce the relative energies of C−O bond stretching calculated from Gaussian 09 illustrated in Figure 6a. The curve of C−O bond dissociation energies from Gaussian 09 calculation is a combination of singlet states and triplet states of CO2. The singlet CO2 has a lower energy up to a C−O bond distance of 2.0 Å, while the energy of the triplet CO2 is lower at larger O− C distances. The C−O bond dissociation curves computed 12012



Article

at ambient temperature and pressures computed from ReaxFF and the classical force field at various concentrations of CO2 are plotted in Figure 7, and the values are given in Table 1. The

Figure 8. Probability density distribution of the C(CO2)−N([Gly]−) distances from MD simulations with the classical force field (red), the ReaxFF force field (blue), and DFT-D2 (green). The simulation cells for classical and ReaxFF calculations contained 10 CO2 with 50 ILs. There were four CO2 molecules with four ILs in the DFT-D2 simulation cell.

Figure 7. Density of [P(C 4 ) 4 ][Gly] as a function of CO 2 concentration as computed from ReaxFF and classical potential33 NpT simulations at ambient conditions. Blue circles: ReaxFF. Red squares: classical potential. Green triangle: experimental density of the pure IL.21

are predictions. The initial structures were generated from the package Packmol.56 The classical simulations give a relatively sharp peak centered at 3.5 Å and a shoulder at about 5.5 Å in the C−N distance, indicating that the interaction between CO2 and the anion is weak and mainly physical in nature. This is expected because only physical interactions are considered within the classical force field framework. Similar peaks around the same C−N distances can also be seen with relatively lower intensity in the ReaxFF and DFT-D2 simulations shown in Figure 8. This indicates the presence of these same physical interactions in these simulations. The striking difference between the classical simulations and the ReaxFF and DFTD2 calculations is the existence of a sharp peak at about 1.7 Å for the latter two cases. This peak is evidence of chemical interaction between CO2 and [Gly]− that is missing in the classical potentials. We have observed similar trends for simulations with different CO2 concentrations using DFT-D2 and ReaxFF, as shown in Figure S1 of the Supporting Information. The probability distributions of the CO2 bend angle computed from classical, ReaxFF, and DFT-D2 methods are plotted in Figure 9. The CO2 bend angles for the classical force field calculations lie within the range of 170−180°, as expected for the absence of chemical interactions. We likewise see a similar peak in the bend angle in the ReaxFF and DFT-D2 simulations, but in addition, there are larger and broader peaks at about 137° for ReaxFF and DFT-D2 calculations. These smaller bend angles indicate CO2 molecules that are chemically bound to the anions, complementing the peaks at small C−N distances shown in Figure 8. The ReaxFF results are predictions because the bend angle of the chemically bound CO2 was not included in the training set. The agreement between the bend angle distributions computed from DFT-D2 and ReaxFF indicates that the description of the energetics of the chemically induced bond bending computed from ReaxFF is more accurate than one would expect from the gas-phase CO2 bend angle energetics shown in Figure 6b. Similar trends are again seen for various concentrations of CO2, as seen in Figure S2 of the Supporting Information. The prefactor for the proton-transfer reaction can be estimated from transition-state theory from the vibrational modes of the transition-state and reactant-state complexes. This

Table 1. Comparison of Density Calculated from Classical and ReaxFF Potenitials for 50 ILs with Increasing Numbers of CO2 within NpT Simulations at 300 K and a Pressure of 1 atma number of CO2

mole % CO2

classical

ReaxFF

0 10 20 25 30 40 50

0 16.7 28.6 33.3 37.5 44.4 50

0.962(0.002) 0.963(0.005) 0.969(0.005) 0.972(0.004) 0.974(0.006) 0.976(0.003) 0.981(0.003)

0.976(0.001) 1.000(0.001) 1.013(0.001) 1.017(0.002) 1.024(0.002) 1.036(0.002) 1.049(0.001)

a

The experimental density of the pure IL at 298.15 K is 0.963 g/cm3.21 The standard deviations of the calculated densities are given in parentheses. Density units are g/cm3.

densities computed from ReaxFF and the classical force field for the pure IL33 are 0.976 and 0.962 g/cm3, respectively. The experimental density is reported to be 0.963 g/cm3.21 The ReaxFF density is larger than the experimental value by about 1%. We note that the density increases with the concentration of CO2 up to at least 50 mol %, as predicted from both the classical and ReaxFF calculations. However, the density increase predicted by ReaxFF is considerably larger than that predicted by the classical potential. We attribute the difference to reactions between CO2 in the [Gly] anion, which naturally would reduce the volume (and thus increase the density) of the system relative to the unreacted (classical) case. This constitutes a prediction that needs to be tested experimentally because, to the best of our knowledge, no experimental data for density as a function of CO2 concentration in [P(C4)4][Gly] have been published. Figure 8 shows the probability distribution of distances between C in CO2 and the nearest N in a [Gly] anion as computed from simulations with the classical and ReaxFF potentials. Also shown is the same quantity computed from VASP DFT-D2 simulations. We note that the VASP data were not included in the training set; therefore, the ReaxFF results 12013



Article

Figure 9. Probability distribution of the CO2 bend angle from MD simulations at 300 K using the classical force field (red), the ReaxFF force field (blue), and the DFT-D2 method (green). The simulation cells for classical and ReaxFF calculations contained 25 CO2 with 50 ILs. There were four CO2molecules with four ILs in the DFT-D2 simulation cell.

Figure 10. RDFs, g(r), for the P atom in the cation and the N atom in the anion for the pure IL (blue line) and a 50 mol % mixture of CO2 and IL (red line).

adduct provides a more symmetric environment for the N atom by sharing of the negative charge between the carboxylic group and the bound CO2, making the charge more diffuse. This means that there are many more orientations of the CO2/ [Gly]− complex that give favorable electrostatic interactions with the P atom in [P(C4)4]+ than are available to the unreacted glycinate anion. These favorable orientations have P−N distances that are close to 6 Å. Results from NpT simulations are qualitatively similar (see Supporting Information Figure S4).

calculation gives a prefactor of 2.84 × 1011 s−1 at 300 K based on DFT-D2 calculations (see the Supporting Information). The rate coefficient calculated using this prefactor and the ReaxFF barrier height (Figure 3) is 1.47 × 1010 s−1. In comparison, the rate coefficient using the DFT-D2 barrier (Figure 3) is 8.34 × 1010 s−1. Thus, we expect to see the proton-transfer reaction occurring roughly on the nanosecond time scale. This estimate has been verified in our ReaxFF MD simulations for a system with 50 CO2 and 50 IL pairs. We observed the proton originating from an anion transferring to a chemically bound CO2, although only in a few cases during a 3 ns total simulation. In contrast, our DFT-D2MD simulations were only carried out for a few picoseconds, and we did not observe the protontransfer event. This shows an advantage of the ReaxFF in accessing longer time scales and observing events that cannot be probed with DFT methods. One would expect that the interactions between the anion and cation of the IL pairs would change significantly upon reaction of the anion with CO2. We have quantified this change by computing the radial distribution function (RDF) for P−N pairs in the pure IL and also for a system where every [Gly] anion is bound to a CO2 molecule. We have carried out these calculations in both the NVT and NpT ensembles (details provided in the Supporting Information). Results for the RDF calculations in the NVT ensemble are presented in Figure 10. The simulations use the same volume for both the pure IL and the equimolar CO2/IL mixture, so that any differences in the RDF are not due to an increase in the average distance between the IL pairs. The RDF for the pure IL shows a small peak at a distance of 3.6 Å that is completely absent in the mixture system. This is to be expected because when CO2 is bound to the anion at the N site, steric hindrance prevents P−N interaction at such small distances. The first large peak for pure [P(C4)4][Gly] occurs at about 5.6 Å and is due to the favorable interaction between the carboxylic group in [Gly]−, which carries most of the charge. This peak corresponds to the carboxylic group in the glycinate oriented toward the P atom in [P(C4)4]+, which causes the N atom to be further away from the phosphorus. Unexpectedly, a new peak appears in the P−N RDF at a distance of about 6.3 Å for the CO2/[P(C4)4][Gly] mixture that is significantly larger than the first peak in the pure IL system. The explanation of this is that the CO2−glycinate

4. CONCLUSION We have obtained a minimum-energy pathway for the reaction of CO2 with tetrabutylphosphonium glycinate at the PBE/ DFT-D2 level of theory with the NEB method, as implemented in VASP. There are two energy barriers of 0.7 and 2 kcal/mol, respectively. The first step involves proton transfer from the N atom in [Gly]− to one of the O atoms in the carboxylic group in [Gly]−. The next step involves movement of the proton to a position where it is shared by the O atom from CO2 bound to the N atom in [Gly]− and the O atom from the carboxylic group in the anion. A ReaxFF force field was developed based on an extensive training set, including the minimum-energy pathway. The optimized force field has been used in NpT MD simulations to determine the equilibrium density at 300 K. The resulting density of 0.976 g/cm3 is within about 1% of the experimental value. We have used ReaxFF to predict the density of CO2/[P(C4)4][Gly] mixtures, and we found a significantly larger density increase with increasing CO2 concentration than predicted from classical, nonreactive simulations. Note that these calculations could not have been carried out with DFT-D3 because the system size is too large and the time scales required to equilibrate the system (several ns) are not accessible to DFT MD simulations. We predict that the enhanced density is due to the compact CO2−[Gly] adduct observed in the ReaxFF calculations. The significant increase in density with CO2 loading can be tested experimentally.

■

ASSOCIATED CONTENT

S Supporting Information *

Optimized ReaxFF parameters, rate coefficient calculation details, probability density distributions for C−N distances and CO2 bend angles, and radial distribution function calculation details. This material is available free of charge via the Internet at http://pubs.acs.org. 12014



■

Article

Theory and Lattice Phonon Dynamics Study. J. Chem. Phys. 2010, 133, 074508−11. (10) Duan, Y.; Zhang, B.; Sorescu, D. C.; Johnson, J. K. CO2 Capture Properties of M−C−O−H (M=Li, Na, K) Systems: A Combined Density Functional Theory and Lattice Phonon Dynamics Study. J. Solid State Chem. 2011, 184, 304−311. (11) Zhang, B.; Duan, Y.; Johnson, K. Density Functional Theory Study of CO2 Capture with Transition Metal Oxides and Hydroxides. J. Chem. Phys. 2012, 136, 064516/1−064516/13. (12) Siriwardane, R. V.; Stevens, R. W. Novel Regenerable Magnesium Hydroxide Sorbents for CO2 Capture at Warm Gas Temperatures. Ind. Eng. Chem. Res. 2008, 48, 2135−2141. (13) Lee, S. C.; Choi, B. Y.; Lee, T. J.; Ryu, C. K.; Ahn, Y. S.; Kim, J. C. CO2 Absorption and Regeneration of Alkali Metal-Based Solid Sorbents. Catal. Today 2006, 111, 385−390. (14) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. High-Throughput Synthesis of Zeolitic Imidazolate Frameworks and Application to CO2 Capture. Science 2008, 319, 939−943. (15) Banerjee, R.; Furukawa, H.; Britt, D.; Knobler, C.; O’Keeffe, M.; Yaghi, O. M. Control of Pore Size and Functionality in Isoreticular Zeolitic Imidazolate Frameworks and Their Carbon Dioxide Selective Capture Properties. J. Am. Chem. Soc. 2009, 131, 3875−3877. (16) Liu, J.; Keskin, S.; Sholl, D. S.; Johnson, J. K. Molecular Simulations and Theoretical Predictions for Adsorption and Diffusion of CH4/H2 and CO2/CH4 Mixtures in ZIFs. J. Phys. Chem. C 2011, 115, 12560−12566. (17) Xie, H.-B.; Johnson, J. K.; Perry, R. J.; Genovese, S.; Wood, B. R. A Computational Study of the Heats of Reaction of Substituted Monoethanolamine with CO2. J. Phys. Chem. A 2010, 115, 342−350. (18) Xie, H.-B.; Zhou, Y.; Zhang, Y.; Johnson, J. K. Reaction Mechanism of Monoethanolamine with CO2 in Aqueous Solution from Molecular Modeling. J. Phys. Chem. A 2010, 114, 11844−11852. (19) Maginn, E. J. Atomistic Simulation of the Thermodynamic and Transport Properties of Ionic Liquids. Acc. Chem. Res. 2007, 40, 1200− 1207. (20) Bates, E. D.; Mayton, R. D.; Ntai, I.; Davis, J. H. CO2 Capture by a Task-Specific Ionic Liquid. J. Am. Chem. Soc. 2002, 124, 926−927. (21) Zhang, J.; Zhang, S.; Dong, K.; Zhang, Y.; Shen, Y.; Lv, X. Supported Absorption of CO2 by Tetrabutylphosphonium Amino Acid Ionic Liquids. Chem.Eur. J. 2006, 12, 4021−4026. (22) Cadena, C.; Anthony, J. L.; Shah, J. K.; Morrow, T. I.; Brennecke, J. F.; Maginn, E. J. Why Is CO2 So Soluble in ImidazoliumBased Ionic Liquids? J. Am. Chem. Soc. 2004, 126, 5300−5308. (23) Yokozeki, A.; Shiflett, M. B.; Junk, C. P.; Grieco, L. M.; Foo, T. Physical and Chemical Absorptions of Carbon Dioxide in RoomTemperature Ionic Liquids. J. Phys. Chem. B 2008, 112, 16654−16663. (24) Gurkan, B. E.; de la Fuente, J. C.; Mindrup, E. M.; Ficke, L. E.; Goodrich, B. F.; Price, E. A.; Schneider, W. F.; Brennecke, J. F. Equimolar CO2 Absorption by Anion-Functionalized Ionic Liquids. J. Am. Chem. Soc. 2010, 132, 2116−2117. (25) Hunt, P. A. The Simulation of Imidazolium-Based Ionic Liquids. Mol. Simul. 2006, 32, 1−10. (26) Gutowski, K. E.; Maginn, E. J. Amine-Functionalized TaskSpecific Ionic Liquids: A Mechanistic Explanation for the Dramatic Increase in Viscosity Upon Complexation with CO2 from Molecular Simulation. J. Am. Chem. Soc. 2008, 130, 14690−14704. (27) Wu, H.; Shah, J. K.; Tenney, C. M.; Rosch, T. W.; Maginn, E. J. Structure and Dynamics of Neat and CO2-Reacted Ionic Liquid Tetrabutylphosphonium 2-Cyanopyrrolide. Ind. Eng. Chem. Res. 2011, 50, 8983−8993. (28) Wu, H.; Maginn, E. J. Water Solubility and Dynamics of CO2 Capture Ionic Liquids Having Aprotic Heterocyclic Anions. Fluid Phase Equilib. 2014, 368, 72−79. (29) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. Reaxff: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 412-624-5644. Notes

This project was funded by the Department of Energy, National Energy Technology Laboratory, an agency of the United States Government, through a support contract with URS Energy & Construction, Inc. Neither the United States Government nor any agency thereof, nor any of their employees, nor URS Energy & Construction, Inc., nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe upon privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This technical effort was performed with the support of the National Energy Technology Laboratory (NETL) through its ongoing research program into CO2 capture (RES Contract DE-FE0004000). Calculations were performed at the University of Pittsburgh Center for Simulation and Modeling. A.C.T.v.D. acknowledges the support of the Fluid Interface Reactions, Structures and Transport (FIRST) Center, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences.

■

REFERENCES

(1) Stern, N. The Economics of Climate Change: The Stern Review; Cambridge University Press: New York, 2007. (2) IPCC Fourth Assessment Report (Ar4); IPCC: Geneva, Switzerland, 2007. (3) Figueroa, J. D.; Fout, T.; Plasynski, S.; McIlvried, H.; Srivastava, R. D. Advances in CO2 Capture TechnologyThe U.S. Department of Energy’s Carbon Sequestration Program. Int. J. Greenhouse Gas Control 2008, 2, 9−20. (4) Astarita, G.; Marrucci, G.; Gioia, F. The Influence of Carbonation Ratio and Total Amine Concentration on Carbon Dioxide Absorption in Aqueous Monoethanolamine Solutions. Chem. Eng. Sci. 1964, 19, 95−103. (5) Supap, T.; Idem, R.; Tontiwachwuthikul, P.; Saiwan, C. Kinetics of Sulfur Dioxide- and Oxygen-Induced Degradation of Aqueous Monoethanolamine Solution During CO2 Absorption from Power Plant Flue Gas Streams. Int. J. Greenhouse Gas Control 2009, 3, 133− 142. (6) Chi, S.; Rochelle, G. T. Oxidative Degradation of Monoethanolamine. Ind. Eng. Chem. Res. 2002, 41, 4178−4186. (7) Abu-Zahra, M. R. M.; Niederer, J. P. M.; Feron, P. H. M.; Versteeg, G. F. CO2 Capture from Power Plants: Part II. A Parametric Study of the Economical Performance Based on Mono-Ethanolamine. Int. J. Greenhouse Gas Control 2007, 1, 135−142. (8) Mosqueda, H. A.; Vazquez, C.; Bosch, P.; Pfeiffer, H. Chemical Sorption of Carbon Dioxide (CO2) on Lithium Oxide (Li2O). Chem. Mater. 2006, 18, 2307−2310. (9) Duan, Y.; Sorescu, D. C. CO2 Capture Properties of Alkaline Earth Metal Oxides and Hydroxides: A Combined Density Functional 12015



Article

(30) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. Reaxff Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (31) LaBrosse, M. R.; Johnson, J. K.; van Duin, A. C. T. Development of a Transferable Reactive Force Field for Cobalt. J. Phys. Chem. A 2010, 114, 5855−5861. (32) Rahaman, O.; van Duin, A. C. T.; Goddard, W. A.; Doren, D. J. Development of a Reaxff Reactive Force Field for Glycine and Application to Solvent Effect and Tautomerization. J. Phys. Chem. B 2010, 115, 249−261. (33) Zhou, G.; Liu, X.; Zhang, S.; Yu, G.; He, H. A Force Field for Molecular Simulation of Tetrabutylphosphonium Amino Acid Ionic Liquids. J. Phys. Chem. B 2007, 111, 7078−7084. (34) Kowsari, M. H.; Alavi, S.; Najafi, B.; Gholizadeh, K.; Dehghanpisheh, E.; Ranjbar, F. Molecular Dynamics Simulations of the Structure and Transport Properties of Tetra-Butylphosphonium Amino Acid Ionic Liquids. Phys. Chem. Chem. Phys. 2011, 13, 8826− 8837. (35) Mortier, W. J.; Ghosh, S. K.; Shankar, S. ElectronegativityEqualization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (36) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358−3363. (37) van Duin, A. C. T.; Baas, J. M. A.; van de Graaf, B. Delft Molecular Mechanics: A New Approach to Hydrocarbon Force Fields. Inclusion of a Geometry-Dependent Charge Calculation. J. Chem. Soc., Faraday Trans. 1994, 90, 2881−2895. (38) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558−561. (39) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal−Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251−14269. (40) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169−11186. (41) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953−17979. (42) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758− 1775. (43) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396−1396. (45) Kristyán, S.; Pulay, P. Can (Semi)Local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175−180. (46) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. Optimization of Effective Atom Centered Potentials for London Dispersion Forces in Density Functional Theory. Phys. Rev. Lett. 2004, 93, 153004. (47) Langreth, D. C.; Dion, M.; Rydberg, H.; Schröder, E.; Hyldgaard, P.; Lundqvist, B. I. Van Der Waals Density Functional Theory with Applications. Int. J. Quantum Chem. 2005, 101, 599−610. (48) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (49) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H−Pu. J. Chem. Phys. 2010, 132, 154104/1−154104/19. (50) Mills, G.; Jónsson, H. Quantum and Thermal Effects in H2 Dissociative Adsorption: Evaluation of Free Energy Barriers in Multidimensional Quantum Systems. Phys. Rev. Lett. 1994, 72, 1124−1127. (51) Mills, G.; Jónsson, H.; Schenter, G. K. Reversible Work Transition State Theory: Application to Dissociative Adsorption of Hydrogen. Surf. Sci. 1995, 324, 305−337.

(52) Jónsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific: Singapore, 1998; pp 385−404. (53) Henkelman, G.; Jonsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978−9985. (54) Frisch, M. J.; et al. Gaussian 09, revision B.01; Gaussian Inc.: Wallingford, CT, 2009. (55) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (56) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164.

12016


Development of a ReaxFF Reactive Force Field for ... - ACS Publications

Recommend Documents