2 BindingbyAlcoholamines AssessmentofDispersion-ImprovedExchange-CorrelationFunctionalsfortheSimulationofCO

(1)

Assessment of Dispersion-Improved Exchange-Correlation Functionals for the Simulation of CO ₂ Binding by

Alcoholamines

Hsueh-Chien Li,

^[a]

Jeng-Da Chai,*

^[b]

and Ming-Kang Tsai*

^[a]

In this study, 12 bound complexes were selected to construct a database for testing 15 dispersion-improved exchange-correlation (XC) functionals, including hybrid generalized gradient approximation (GGA), modified using the Grimme’s pairwise strategy, and double hybrid XC functionals, for specifically characterizing the CO₂ binding by alcoholamines. Bound complexes were selected based on the characteristics of their hydrogen bonds, dispersion, and electrostatic (particularly between the positive charge of CO₂ and the lone pair of N of alcoholamines) interactions. The extrapolated binding energy from the aug-cc-pVTZ (ATZ) to aug-cc-pVQZ (AQZ) basis set at the CCSD(T)/CBS(MP21DZ) level was used as the reference for the XC functional comparison. M06-2X produced the optimal agreement if the optimized geometries at MP2/ATZ level were chosen for all the test bound complexes. However, M06-L,

xB97X, and xB97, and were preferred if the corresponding density functional theory (DFT) optimized geometries were adapted for the benchmark. Simple bimolecular reaction between CO₂ and monoethanolamine simulated using polarizable continuum solvation model confirmed thatxB97, xB97X, and xB97XD qualitatively reproduced the energetics of MP2 level. The inconsistent performance of the tested XC functionals, observed when using MP2 or DFT optimized geometries, raised concerns regarding using the single-point ab initio correction combined with DFT optimized geometry, particularly for determining the nucleophilic attack by alcoholamines to CO₂.V^C 2014 Wiley Periodicals, Inc.

DOI: 10.1002/qua.24670

Introduction

Capturing and storing CO₂ is crucial for ameliorating global climate change. Using alcoholamines, such as monoethanolamine (MEA), diethanolamine (DEA), and methyldiethanol- amine (MDEA), is one of the commercially available technologies involving solvent-based processes to immobilize the CO₂ produced by fossil-fuel power stations.^[1] Alcohol- amines are typically mixed with water to form carbamate solu- tions in the capturing process. The fundamental physics of CO₂ binding by alcoholamines is to take advantage of the electrostatic interaction between the lone pair on the N of alcoholamines, which undergoes a nucleophilic attack, and the partial-positively charged C of CO₂. Many experimental and theoretical studies have focused on elucidating the reaction mechanism of the MEA/CO₂system to understand the dynamics and kinetics of the CO₂ capturing process, and ultimately reduce the operational cost of the solvent-based CO₂-capture technology.^[2–14]

A schematic representation of the CO₂ binding moiety is shown in Figure 1. From the classical intermolecular interaction perspective, the bound complex (or a physically absorbed complex) is formed by the quadrupole–dipole interaction. The alpha-C of alcoholamines can be further polarized as CO₂ approaches, in addition to the intrinsic polarization resulting from the electron-withdrawing effect of N due to the formation of zwitterion. The structures of the transition state and zwitterion (Fig. 1) can be intuitively considered as the chemical bond formation through the

orbital overlapping between the p orbital of MEA and the low-unoccupied molecular orbital (LUMO) of CO₂, based on the quantum perspective. Accurately describing the CO₂ capture process by alcoholamines is therefore a challenging theoretical task. The physical-absorbed complex can be simply described by weak-electrostatic and dispersion interactions. As the alkylation of the N of alcoholamine increases, the significance of the dispersion interaction also increases.

Moreover, the chemical-bound zwitterionic species is formed by a bonding orbital consisting of the bent-p^* orbital of CO₂ and ther-character lone-pair electron of N. The lack of appropriate description for the dispersion interactions, as well as the low-lying unoccupied orbital along the CO₂ absorption reaction coordinate, could produce controversial results when interpreting the experimental observations.

[a] H.-C. Li, M.-K. Tsai

Department of Chemistry, National Taiwan Normal University, Taipei 11677, Taiwan

E-mail: [email protected] [b] J.-D. Chai

Department of Physics, Center for Theoretical Sciences, and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan

E-mail: [email protected]

Contract grant sponsor: National Science Council; contract grant numbers:

101-2113-M-003-003-MY2, 101-2112-M-002-017-MY3.

Contract grant sponsor: National Taiwan University; contract grant numbers: 99R70304, 101R891401, 101R891403 (J.D.C.).

Contract grant sponsor: National Center for Theoretical Sciences of Taiwan.

V^C2014 Wiley Periodicals, Inc.

(2)

Many theoretical works have focused on outlining the reaction mechanism of the capturing process. Three types of mechanisms were summarized by Xie et al., including the (1) zwitterion mechanism,^[2,4,13] (2) single-step mechanism,^[9,10]

and (3) carbamic acid reaction mechanism.^[12] A sophisticated effort to differentiate the various proposed mechanisms by simulating the explicit solvation effects, dynamics of the intermolecular proton transfer process, and the various substitution effects to N of alcoholamines has been attempted by several studies in the past,^[13,15–20] but is beyond the scope in this study. In this study, we aimed to search for a simple and accessible method at the density functional theory (DFT) level to describe CO2–alcoholamine interactions based on the design of a three-category database (denoted as C1). The par- tition of the C1 was meant to examine the selected XC functionals separately in the category of NH3 and MEA, DEA, and CO2-binding. The presence of NH3, MEA, and DEA also accounted for the basicity resulted from the alkylation effect on N of alcoholmines, and triethanolamine was omitted due to its slow CO2uptake.^[21]

The ideal description of CO2–alcoholamines interactions using the ab initio methods and triple-zeta quality basis set could be accessible for the bimolecular reaction like CO2 and MEA. More substitution on N of amines or even accounting the microsolvation effect would be computationally impractical at such levels. Therefore, finding suitable XC functionals to provide qualitative description of the important intermediate and minimum structures would be useful for the future design the CO2 scrubbers computationally. The assessments of XC function for the noncovalent interactions have been docu- mented in the past,^[22–25]and the S22 and S66 databases were commonly used for the benchmark evaluation.^[26,27] These databases consist of a wide range of noncovalent and biologi- cally important complexes, and show most of the essential physics of the intermolecular interactions in nature. However, in this study, we aimed to search an optimal method for the specific application of CO2 binding by alcoholamines. As a consequence, we systematically developed a testing database including 12 bound complexes ranging from pure hydrogen- bond (HB) complexes, mixed HB, and dispersion complexes, to dispersion-dominated complexes. Fourteen functionals reported as being able to improve dispersion interaction were chosen for the benchmark tests; these functionals included two generalized gradient approximation (GGA) modified by the Grimme’s semiempirical function, 11 hybrid GGA functionals, and two perturbed GGA functionals. The details of the database design, as well as the ab initio reference data, are

presented in Section Methodology. The results and discussion are presented in Section Discussion. The conclusion is summarized in Section Conclusions.

Methodology

C1 database

Nine HB interacting complexes were included in the C1 database (Fig. 2). These complexes included pure HB complexes such as (H₂O)₂ and (H₂O)NH₃, mixed HB, and dispersion-bound complexes of MEA and DEA with H₂O and NH₃, and dispersion- dominant bound complexes such as (MEA)₂, (MEA)DEA, and (DEA)₂. In addition, the binding interactions of CO₂ with NH₃, MEA, and DEA were also included. The C1 database was con- structed chemical-intuitively, accounting for the electronic effect resulting from the various levels of alkylation on the N of alcoholamines, the consequent intermolecular HB interactions, and the CO₂ binding. Moreover, the structures of the testing complexes were chosen based on their potential as stable intermedi- ates along the reaction coordinate of the capturing process. In contrast to the conventional S22 database, the dispersion interaction resulting from the p electrons was omitted in the C1 database. Moreover, the r electron was well represented and well suited for studying CO₂ binding resulting from the nucleophilic attack of alcoholamines. The largest complex in C1 was 36 atoms, whereas the largest complex in S22 was 28 atoms.

Ab initio electronic structure calculation

MP2/aug-cc-pVTZ (ATZ) optimized geometries were used to conduct all the ab initio reference calculations, where ATZ included 4s3p2d1f basis functions for the second row ele- ments. Both couple cluster singles and doubles (CCSD)(T)/

CBS(MP21DZ) and MP2 at the complete basis set limit (CBS) were adopted to calculate the binding energies of complexes a–l listed in C1 database. The CBS was approximated using the extrapolation from ATZ to aug-cc-pVQZ (AQZ). The formula of the CCSD(T)/CBS(MP21DZ) scheme adopted in this study was listed as Eq. 1^[28–30]:

DE_CBS^CCSD_ðMP21DZÞ^ðTÞ 5DECBS^MP2ðTZ!QZÞ1ðDE^CCSD^ðTÞ2DE^MP2Þsmall basis set (1)

where each DE represents the binding energy with the basis set superposition error (BSSE) correction. The superscript ofDE denotes the theory level of the BSSE calculations. Despite the counterpoise approach for correcting the binding energy for basis set incompleteness generally results in underestimated values, in particular for hydrogen bonded systems, ^[25,31] the tested cases in this study actually covered species from pure HB to dispersion-dominate complexes. The CBS limit at the CCSD(T)/CBS(MP21DZ) level, denoted as DE_CBS^CCSDðTÞ_{ðMP 21DZ Þ} was approximated byDE_CBS^{MP 2}_{ðTZ !QZ Þ} in addition to a correction term of evaluating the difference between CCSD(T) and MP2 energy with using aug-cc-pVDZ basis set. As being introduced by Hobza et al., the DE_CBS^CCSDðTÞ_{ðMP 21DZ Þ} scheme was based on the assumption that the difference between the binding energies at CCSD(T) and MP2 levels converges faster than the binding Figure 1. Schematic representation of CO2 capturing process by MEA.

[Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

(3)

energies themselves in respect to the quality of the basis set.

[28,29]DE^{MP 2}_CBS_{ðTZ !QZ Þ}adopted the formula in Eq. 2, where X and Y represented ATZ and AQZ in the extrapolation, respectively.

DE^MP2_XY 5DE^HF_Y 1DEX^corr3n³_X2DEY^corr3n³_Y

n³_X2n³_Y (2)

The first term,DE^HF_Y , of Eq. 2 was calculated at the HF/AQZ level, where DE^corr_{Xor Y} terms were the difference of the correlation energy in the binding energy calculation at the MP2 level using basis set X and Y, respectively. The correlation contribu- tion was further weighted by n, and n 5 2–4 if the basis set was aug-cc-pV(T or Q)Z, respectively. Equation 2 was also adopted to calculateDE^CCSD_TQ ^ðTÞ andDE^CCSD_DT ^ðTÞ using different ab initio method and basis set extrapolation, and DE_{Xor Y}^corr was adjusted consequently.

DFT functionals

Fifteen exchange-correlation (XC) functionals were chosen for testing the statistics error. Five functionals (CAM-B3LYP,^[32] LC- xPBE,^[33–35] xB97,^[36] xB97X,^[36] and xB97XD^[37]) were long- range corrected by adding the Gauss error function in the exchange term. Three dispersion-improved functionals (BLYP- GD2,^[38,39] B97-GD2,^[40] B3LYP-GD2,^[38,39,41] and B3LYP-GD3) were corrected using the Grimme’s D2 and D3 dispersion, respectively.^[42,43] Four Minnesota series hybrid-GGA functionals were included (M05-2X,^[44] M06-2X,^[45] M06-HF,^[46,47] and M06-L^[48]). For two double hybrid GGA functionals (B2PLYP and B2PLYPD), the PT2 scheme was included in the correlation calculation.^[49] All DFT calculations were performed using the Gaussian 09 package^[50] and NWChem.^[51] The geometries

used for the BSSE calculations were optimized using the corresponding XC functional and ATZ basis set.

Discussion

Converging ab initio references using MP2/ATZ optimized geometries

Four types of calculations including CCSD(T),DCCSD(T), and MP2, coupled with different basis set extrapolations (i.e., ATZ ! AQZ and aug-cc-pVDZ (ADZ)! ATZ) were adopted to determine the selection of ab initio reference data, and the results are summarized in Table 1. DE_CBS^CCSD_{ðMP 21DZ Þ}^ðTÞ systematically but marginally underestimated the binding energy of the small complexes (a–d) in comparison with DE^{MP 2}_CBS_{ðTZ !QZ Þ}. However, D E^CCSD_CBS_{ðMP 21DZ Þ}^ðTÞ provided the optimal agreement with DE_CBS^{MP 2}_{ðTZ !QZ Þ} for the CO2 alcoholamine bound complexes (j–k). For the non-CO2 complexes (a–

i), DE_CBS^{MP 2}_{ðTZ !QZ Þ} predicted slightly stronger binding energy thanDE_CBS^CCSD_{ðMP 21DZ Þ}^ðTÞ andDE_CBS^CCSD_{ðDZ !TZ Þ}^ðTÞ . Given the interest of studying CO2binding by alcoholamines, we, therefore, selectedDE_CBS^CCSDðTÞ_{ðMP 21DZ Þ}as the reference for comparison with the DFT results.

Sensitivity of the selection of MP2/ATZ optimized geometries To determine the accuracy of the MP2/ATZ optimized geometries, we also calculatedDE^{MP 2}_CBS_{ðTZ !QZ Þ} using the MP2/ADZ and CCSD/ADZ optimized geometries (Table 2), and compared several selected interatomic distances. Geometrically, the MP2/ATZ optimized geometries were close to the CCSD/ADZ optimized Figure 2. Graphical representation of C1 database, that is, a) is (H2O)2, b) is NH3(H2O), c) is MEA(H2O), d)

is MEA(NH3), e) is DEA(H2O), f ) is DEA(NH3), g) is (MEA)2, h) is MEA(DEA), i) is (DEA)2, j) is NH3(CO2), k) is MEA(CO2), and l) is DEA(CO2). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

Table 1. Binding energy (kcal/mol) of complexes a–l^[a]at the various basis set extrapolations.

Complexes DECBSðTZ !QZ Þ^CCSDðTÞ DE^CCSDCBSðDZ !TZ Þ^ðTÞ DE^CCSDCBSðMP 21DZ Þ^ðTÞ DECBSðTZ !QZ Þ^{MP 2}

A 25.06 24.95 24.97 25.02

B 26.60 26.52 26.48 26.66

C 27.93 27.85 27.80 28.07

D 24.80 24.77 24.73 24.85

E – 28.07 28.01 28.25

F – 24.91 24.87 25.07

G – 24.50 24.47 24.56

H – – 25.62 25.84

I – – 212.22 212.92

J 23.15 23.08 23.10 22.96

K 24.32 24.24 24.26 24.18

L – 24.25 24.26 24.39

Dashes denote the data are unavailable due to the computational expense. [a] MP2/ATZ optimized geometries were used.

(4)

geometries, without any relative interatomic distance being greater than 0.07 A˚ , whereas a negligible difference was observed when using the ADZ and ATZ basis sets at the MP2 level. However, the binding energy extrapolated using the MP2/ATZ optimized geometries was slightly closer to the energy extrapolated using the CCSD/ADZ optimized geometries. The energies extrapolated using the MP2/ADZ geometries were also in good agreement with the ATZ basis set.

Therefore, using MP2/ATZ optimized geometries provided a reasonable accuracy for the CCSD level, geometrically and energetically. Although an additional examination of the CCSD/ATZ optimized geometries would be preferred, computationally it would be too expensive or impractical.

XC functional comparison using MP2/ATZ optimized geometry

Table 3 summarizes the overall performance of the 15 functionals tested, using the complexes a–l listed in Figure 2. The root mean square deviations (RMS) of 10 out of 15 functionals were observed to be less than 1.0 kcal/mol. The optimal approach was using M06-2X. The three functionals corrected using the Grimme’s dispersion approach (B3LYP-GD3, B97-GD2,

and BLYP-GD2) also produced a good agreement, having RMS values of 0.56, 0.58, and 0.81, respectively. Nevertheless, the two double-hybrid functionals (B2PLYPD and B2PLYP) did not rank highly among these 15 cases.

Complexes a–l can be further categorized into three subgroups, as summarized in Table 4. The first group, labeled

“small molecule and MEA series,” included a–d, g, j, and k; the second subgroup, labeled “DEA series,” included complexes e, f, h, i, and l, where the dispersion interaction played a crucial role; and the last subgroup, called the CO₂-bound series, con- tained complexes j–l. Both M06-2X and M06-L were effective in all three subcategories. Grimme’s dispersion correction showed the optimal agreement in the DEA series, particularly for B97-GD2 and B3LYP-GD3. The double hybrid functionals performed poorly in all subgroups, and their performance worsened in the dispersion-dominated DEA series. In the CO₂- bound series, 10 functionals (M05-2X, B3LYP-GD3, BLYP-GD2, B3LYP-GD2,xB97XD, xB97X, M06HF, M06-2X, xB97, and M06- L) were found to predict CO₂ binding energy with RMS error less than 0.6 kcal/mol. Seven out of these ten functionals overestimated the binding energy of CO2 with alcoholamines, as shown by the negative value of the corresponding mean signed error (MSE). Accounting DE_CBS^CCSDðTÞ_{ðMP 21DZ Þ} gave slightly under-bound character of NH3(CO2) and MEA(CO2) with respect toDE_CBS^CCSDðTÞ_{ðTZ !QZ Þ}(Table 1), the over-bound character of these seven functionals may provide an accurate description of the energy required for the CO2capture by alcoholamines.

XC functional comparison using DFT optimized geometry Tables 5 and 6 list the overall and categorized statistical analysis of the relative binding energy between DFT and DE^CCSD_CBS_{ðMP 21DZ Þ}^ðTÞ of complexes a–l, respectively. The binding energy at the DFT level was calculated using the geometry optimized for each corresponding functional whereas the ab initio reference used the MP2/ATZ optimized geometry. The performance of the tested functionals was substantially different from the previous analysis using MP2/ATZ geometries (Table 3). The reordering raised concerns about combing DFT optimized geometries and a single-point ab initio energetic correction. Only four functionals (M06-L, xB97X, xB97, and M05-2X) were predicted to have RMS < 0.6 kcal/mol. Surpris- ingly, M06-2X, the optimal functional listed in Table 3, dropped to almost the last position, with an RMS of 3.02 kcal/mol due Table 2. Comparison of the optimized geometries at MP2/ATZ, MP2/ADZ,

and CCSD/ADZ.

Complexes MP2/ATZ MP2/ADZ CCSD/ADZ

r_ONof (H2O)(NH3) 2.923 2.938 (0.014) 2.980 (0.057) r_NCof (NH3)(CO2) 2.937 2.933 (20.004) 2.950 (0.013) r_NN of MEA(NH3) 3.208 3.201 (20.007) 3.269 (0.062) r_NOof MEA(H2O) 2.868 2.882 (0.013) 2.930 (0.062) r_NCof MEA(CO2) 2.857 2.861 (0.003) 2.896 (0.038) r_NN of (MEA)2 3.204 3.208 (0.005) 3.270 (0.064) Complexes DECBSðTZ !QZ Þ^{MP 2} //

MP2/ATZ^[a]

DECBSðTZ !QZ Þ^{MP 2} //

MP2/ADZ^[a]

DE^{MP 2}CBSðTZ !QZ Þ//

CCSD/ADZ^[a]

(H2O)(NH3) 26.66 26.72 (20.05) 26.61 (0.05) (NH3)(CO2) 22.96 22.96 (0.00) 22.97 (20.01)

MEA(NH3) 24.85 24.93 (20.08) 24.88 (20.03)

MEA(H2O) 28.07 28.15 (20.08) 27.97 (0.10)

MEA(CO2) 24.18 24.22 (20.04) 24.04 (0.14)

(MEA)2 24.56 24.63 (20.08) 24.61 (20.05)

The extrapolated binding energy was in kcal/mol and the intermolecular distance was in A˚. The parentheses denoted the relative binding energy and interatomic distance in respect to the values listed in the first column. [a] Each column used MP2/ATZ, MP2/ADZ, and CCSD/ADZ optimized geometries, respectively.

Table 3. Overall performance (kcal/mol) of the tested 15 functionals.^[a]

Methods RMS MSE MAE Methods RMS MSE MAE

M06-2X 0.30 20.08 0.25 xB97XD 0.64 20.33 0.50

M06-L 0.45 0.41 0.43 BLYP-GD2 0.81 20.46 0.58

B3LYP-GD3 0.56 20.45 0.45 B3LYP-GD2 1.07 20.84 0.84

xB97 0.57 20.47 0.51 B2PLYPD 1.46 1.13 1.13

B97-GD2 0.58 0.06 0.45 CAM-B3LYP 2.06 1.37 1.47

xB97X 0.59 20.18 0.50 LC-xPBE 2.54 2.02 2.02

M06HF 0.60 0.20 0.48 B2PLYP 3.72 2.90 2.90

M05-2X 0.61 0.22 0.39

RMS is the root mean square deviation. MSE stands for the mean signed error. MAE is mean absolute error. [a] MP2/ATZ optimized geometries were used for complexes a–l for all functionals.

(5)

to its different optimized geometry in the DEA series. The other two functionals treated using Grimme’s correction (B3LYP-GD3 and B97-GD2) also showed an acceptable agreement regarding the CCSD(T)/CBS(MP21DZ) level. The double hybrid functional still showed a considerable statistical deviation, even when using the corresponding optimized geometries.

Similarly, dispersion-dominated series regardless using MP2 or DFT geometries were generally troublesome for the tested functionals as shown in Table 6. Only three functionals (xB97, M06-L, and xB97X) had RMS < 0.6 kcal/mol, and the double hybrid performed poorly in each of the subcategories. Ten functionals (BLYP-GD2, M05-2X, B3LYP-GD3, xB97X, B3LYP-GD2, xB97, M06-2X, xB97XD, M06-L, and B97-GD2) described the CO2 interaction with alcoholamines effectively as

listed in CO2 bound series of Table 6. B97-GD2, as being performed strongly in small molecule and MEA subcategories, produced RMS of 0.59 kcal/mol with the under-bound character.

Table 7 summarized the statistical results of the various interatomic measurements, to show the consistency of the optimized geometries in each theoretical method. The DEA dimer was selected because the dominance of the dispersion interaction. The reference used the MP2/ATZ optimized geometry for the DEA dimer, and three interatomic distances (N13- N20, O35-N13, and O15-O33) were selected for the comparison (Fig. 3). Although the three interatomic measurements could be simply considered as hydrogen-bonding interactions, the mutual interaction between the backbone alcohol groups was believed to substantially affect these interatomic distances.

Thus, the description of the dispersion interaction of each Table 4. Statistical analysis of the relative binding energy in kcal/mol of the subcategories using MP2/ATZ optimized geometries.^[a]

Small Molecule and MEA series^[b]

M06-2X 0.27 20.13 0.24 xB97X 0.54 20.48 0.48

M05-2X 0.29 20.01 0.25 B3LYP-GD2 0.65 20.55 0.55

xB97XD 0.32 20.12 0.29 xB97 0.68 20.66 0.66

B3LYP-GD3 0.39 20.33 0.33 B2PLYPD 0.72 0.64 0.64

BLYP-GD2 0.43 20.19 0.39 CAM-B3LYP 0.88 0.55 0.71

M06HF 0.45 0.13 0.36 LC-xPBE 1.38 1.23 1.23

M06-L 0.48 0.46 0.46 B2PLYP 1.86 1.68 1.68

B97-GD2 0.51 0.32 0.37

DEA series^[c]

M06-2X 0.33 20.02 0.26 xB97XD 0.92 20.64 0.79

xB97 0.36 20.20 0.30 BLYP-GD2 1.14 20.83 0.85

M06-L 0.41 0.33 0.38 B3LYP-GD2 1.48 21.23 1.23

xB97X 0.66 0.25 0.52 B2PLYPD 2.10 1.81 1.81

B97-GD2 0.68 20.31 0.56 CAM-B3LYP 3.02 2.52 2.52

B3LYP-GD3 0.74 20.63 0.63 LC-xPBE 3.58 3.12 3.12

M06HF 0.76 0.30 0.65 B2PLYP 5.33 4.62 4.62

M05-2X 0.87 0.55 0.59

CO2bound series^[d]

M05-2X 0.18 20.17 0.17 xB97 0.53 20.50 0.50

B3LYP-GD3 0.24 20.21 0.21 M06-L 0.59 0.57 0.57

BLYP-GD2 0.27 0.24 0.24 B97-GD2 0.76 0.75 0.75

B3LYP-GD2 0.29 20.25 0.25 B2PLYPD 1.07 0.95 0.95

xB97XD 0.34 0.33 0.33 CAM-B3LYP 1.54 1.34 1.34

xB97X 0.34 20.12 0.32 LC-xPBE 2.10 1.95 1.95

M06HF 0.36 20.35 0.35 B2PLYP 2.63 2.38 2.38

M06-2X 0.49 20.48 0.48

[a] MP2/ATZ optimized geometries were used for complexes a–l for all functionals. [b] Including complexes a–d, g, j, and k. [c] Including complexes e, f, h, i, and l. [d] Including complexes j, k and l.

Table 5. Statistical analysis of the relative binding energy in kcal/mol of complexes a–l using the optimized geometry at the corresponding theory level.

M06-L 0.46 0.36 0.44 BLYP-GD2 1.27 20.67 0.91

xB97X 0.51 20.20 0.45 B3LYP-GD2 1.28 20.99 0.99

xB97 0.54 20.47 0.47 CAM-B3LYP 1.68 1.15 1.26

M05-2X 0.59 0.20 0.39 B2PLYPD 1.74 0.76 1.34

M06HF 0.65 0.06 0.54 LC-xPBE 2.21 1.81 1.81

B3LYP-GD3 0.71 20.57 0.57 M06-2X 3.02 21.30 1.36

xB97XD 0.71 20.34 0.53 B2PLYP 3.03 2.43 2.43

B97-GD2 0.76 20.16 0.54

DE^CCSDCBSðMP 21DZ Þ^ðTÞ //MP2/ATZ was used as the reference to calculate the errors of DFT.

(6)

functional could be examined in terms of geometrical variance from the MP2 reference. The results in Table 7 were sorted with respect to the size of the RMS of the DFT geometrical data, and all the DFT optimizations were started from the MP2 optimized geometry. B2PLYPD exhibited the optimal agreement with the MP2 geometry, but its binding energy between two DEA molecules was considerably underestimated. The

three functionals treated using Grimme’s correction were within a 0.1 A˚ deviation from the reference, whereas the corresponding DEA dimer binding energy was also in reasonable agreement. M06-2X (Supporting Information Fig. S1) was observed to locate at a slightly different minimum in comparison with MP2 and other XC functionals, where additional HB was formed between O15-O33 and N13-N20 was elongated.

Table 6. Statistical analysis of the relative binding energy in kcal/mol of the subcategories using the optimized geometries of the corresponding theory level.^[a]

Small Molecule and MEA series^[b]

M06-2X 0.23 20.07 0.19 xB97 0.66 20.60 0.60

M05-2X 0.30 20.04 0.26 BLYP-GD2 0.71 20.18 0.60

xB97XD 0.39 20.06 0.33 B3LYP-GD2 0.73 20.64 0.64

B3LYP-GD3 0.44 20.38 0.38 CAM-B3LYP 0.79 0.49 0.67

B97-GD2 0.44 0.20 0.36 LC-xPBE 1.27 1.16 1.16

M06HF 0.46 0.06 0.37 B2PLYPD 1.47 0.03 1.02

M06-L 0.48 0.46 0.46 B2PLYP 1.61 1.49 1.49

xB97X 0.53 20.47 0.47

DEA series^[c]

xB97 0.33 20.27 0.27 BLYP-GD2 1.78 21.35 1.35

M06-L 0.45 0.23 0.43 B3LYP-GD2 1.79 21.48 1.48

xB97X 0.49 0.18 0.42 B2PLYPD 2.06 1.78 1.78

M05-2X 0.84 0.53 0.59 CAM-B3LYP 2.44 2.08 2.08

M06HF 0.85 0.06 0.76 LC-xPBE 3.08 2.73 2.73

B3LYP-GD3 0.97 20.82 0.82 B2PLYP 4.28 3.76 3.76

xB97XD 1.00 20.73 0.81 M06-2X 4.67 23.01 3.01

B97-GD2 1.05 20.65 0.79

CO2bound series^[d]

BLYP-GD2 0.22 0.09 0.21 M06-L 0.58 0.56 0.56

M05-2X 0.26 20.25 0.25 B97-GD2 0.59 0.57 0.57

B3LYP-GD3 0.32 20.28 0.28 M06HF 0.64 20.54 0.54

xB97X 0.35 20.15 0.33 B2PLYPD 1.03 0.92 0.92

B3LYP-GD2 0.37 20.32 0.32 CAM-B3LYP 1.39 1.22 1.22

xB97 0.45 20.35 0.35 LC-xPBE 1.91 1.79 1.79

M06-2X 0.45 20.35 0.38 B2PLYP 2.21 2.03 2.03

xB97XD 0.46 0.39 0.39

[a]DE^CCSDCBSðMP 21DZ Þ^ðTÞ //MP2/ATZ was used as the reference to calculate the errors of DFT. [b] Including complexes a-d, g, j, and k. [c] Including complexes e, f, h, i, and l. [d] Including complexes j, k, and l.

Table 7. Geometry analysis (A˚ ) and binding energy (kcal/mol) using the corresponding optimized geometries.

N13-N20^[a] O35-N13^[a] O15-O33^[a] BE RMS MSE MAE

MP2^[b] 2.996 2.846 4.694 212.22

B2PLYPD 2.987 2.841 4.698 28.49 6.21E-03 23.23E-03 5.83E-03

B97-GD2 2.961 2.830 4.720 214.36 2.69E-02 28.35E-03 2.58E-02

M06HF 3.032 2.804 4.663 211.04 3.67E-02 21.23E-02 3.64E-02

xB97 3.078 2.910 4.726 212.28 6.27E-02 5.91E-02 5.91E-02

BLYP-GD2 2.940 2.820 4.569 215.75 8.02E-02 6.89E-02 6.89E-02

B3LYP-GD2 2.951 2.821 4.564 215.62 8.06E-02 26.67E-02 6.67E-02

xB97XD 3.032 2.865 4.879 214.13 1.10E-01 8.00E-02 8.00E-02

M05-2X 3.110 2.917 4.860 210.54 1.24E-01 21.17E-01 1.17E-01

M06-L 3.091 2.922 4.419 212.71 1.73E-01 20.35E-02 1.48E-01

xB97X 3.095 2.905 4.991 211.30 1.84E-01 1.52E-01 1.52E-01

B3LYP-GD3 3.015 2.849 5.050 213.99 2.05E-01 1.27E-01 1.27E-01

B2PLYP 3.104 2.898 5.448 24.43 4.41E-01 3.05E-01 3.05E-01

LC-xPBE 3.153 2.895 5.451 26.80 4.47E-01 3.21E-01 3.21E-01

CAM-B3LYP 3.125 2.895 5.614 27.86 5.37E-01 3.66E-01 3.66E-01

M06-2X 3.203 2.897 2.967 221.08 1.00E100 24.90E-01 6.61E-01

[a] The interatomic distances are labeled in Figure 3. [b] MP2/ATZ optimized geometry and binding energy atDECBSðMP 21DZ Þ^CCSD^ðTÞ level.

(7)

Bimoelcular and trimolecular reaction mechanism test Nine XC functionals (B3LYP-GD2, B3LYP-GD3, M05-2X, M06-2X, M06-HF, xB97, xB97X, xB97XD, and BLYP-GD2) were selected to calculate the bimolecular CO2AMEA and trimolecular CO₂A(MEA)₂ binding mechanism as described in Figure 1 and Supporting Information Figure S1 using ATZ basis set, compared with the results of MP2/ATZ. The solvent effect was approximated by the polarizable continuum solvation model (PCM) and the dielectric constant used 24.852 as ethanol. The calculated barriers and reaction energies were shown in Tables 8 and 9.xB97, xB97X, and xB97XD were found as the optimal choices and gave qualitatively agreement with the MP2 results in describing the transition state with less than 0.7 kcal/mol deviation, and acceptable accuracy at <1.6 kcal/mol in estimating the reaction energy using the PCM approximation for the bimolecular and trimolecular reactions. Combining implicit solvent description and BLYP-GD2 was not able to locate the transition state and zwitterion along the bimolecular reaction pathway but underestimated the barrier by more than 2.5 kcal/mol for the trimolecular reaction pathway. The inconsistent results of BLYP-GD2 may be inferred to the underestimated kinetics of CO2 absorption in dilute MEA solution but overestimated kinetics in liquid MEA in the condensed phase simulations.

Conclusions

We built a specific but systematic database for describing the potential interactions involved in the CO2 capture by alcoholamines, including hydrogen bond, dispersion, and electrostatic (positive charge of CO2vs. the lone pair of N of alcoholamines) interactions. Using the proposed database, we tested 15 XC dispersion-improved XC functionals to determine the optimal choice for the future theoretical characterization of the CO2

capturing process. We used the binding energy at the CCSD(T)/CBS(MP21DZ) level extrapolated from ATZ to AQZ level as a reference because it had the smallest deviation, compared with the extrapolated results of CCSD(T), for small bound complexes such as (H2O)2, NH3(H2O), and NH3(CO2).

Using MP2/ATZ-optimized geometry, we determined that M06- 2X was the optimal choice, whereas the other eight XC functionals also had RMS < 0.6 kcal/mol in the overall comparison.

Using the corresponding DFT optimized geometries com- pletely changed the ranking order, indicating that combining DFT optimized geometries with an ab initio single-point correction could be questionable.xB97XD was the optimal option, followed by B97-GD2 and xB97; all three functionals had RMS < 0.6 kcal/mol. xB97X and BLYP-GD2 were the next two optimal options, having RMS of 0.76 and 0.87 kcal/mol, respectively. Three out of 15 functionals (xB97XD, B97-GD2, and xB97) described the dispersion-dominated DEA series effectively and had RMS < 0.6 kcal/mol. Additionally, using DFT optimized geometries, we identified six functionals (M05-2X, BLYP- GD2, B3LYP-D, xB97X, M06-2X, and xB97) as being suitable (RMS < 0.6 kcal/mol) for describing CO2-alcoholamine interactions. Furthermore, by filtering out the XC functional selection through comparison of the crucial geometrical descriptors listed in Table 7 (the criteria was chosen as <10²³A˚ ), we identified B97-GD2, BLYP-GD2, and xB97 as the recommended choices for characterizing the process of CO2 capture by alcoholamine.xB97 is available in most computational chemistry packages and is suitable for simulating intermolecular interactions. B97-GD2 and BLYP-GD2 are available for the condensed- Figure 3. Local minimum structure of DEA dimer at MP2 level. [Color figure

can be viewed in the online issue, which is available at wileyonlinelibrary.

com.]

Table 8. The calculated relative barrier and reaction energy for the bimolecular CO2AMEA reaction (in kcal/mol) using PCM model.

Theory DETS DERec DDE_TS^{MP 2} DDE_Rec^{MP 2}

MP2 4.86 3.72 0.0 0.0

M06-HF 0.68 27.75 24.18 211.47

M05-2X 2.52 20.35 22.34 24.07

B3LYP-GD2 3.11 3.09 21.75 20.63

B3LYP-GD3 3.43 21.43 21.43 20.39

M06-2X 3.61 1.45 21.26 22.27

xB97XD 4.46 2.80 20.40 20.92

xB97X 4.87 2.58 0.01 21.14

xB97 5.15 2.44 0.28 21.28

BLYP-GD2^[a] – – – –

DETS andDERec are the barrier and reaction energy predicted at each corresponding theory level. DDE_TS^{MP 2} and DDE^{MP 2}_Rec are the errors in respect to MP2 results. [a] TS and zwitterion minimum structures of BLYP-GD2 using PCM model were not able to locate.

Table 9. The calculated barrier and reaction energy for the trimolecular CO2A(MEA)2reaction (in kcal/mol) using PCM model as shown in Sup- porting Information Figure S2.^[a]

Theory DETS 1 DE122 DE123 DDE^{MP 2}_{TS 1} DDE^{MP 2}₁₂₂ DDE^{MP 2}₁₂₃

MP2 2.70 23.30 23.52 0.00 0.00 0.00

M06-HF 0.10 216.61 219.00 22.60 213.31 215.48

M05-2X 1.45 26.63 28.17 21.25 23.32 24.65

B3LYP-GD2 0.90 23.64 23.80 21.80 20.34 20.28

B3LYP-GD3 1.26 22.80 22.80 21.44 0.50 0.72

M06-2X 2.05 25.16 25.63 20.65 21.86 22.12

xB97XD 2.49 24.34 24.84 20.22 20.71 21.32

xB97X 3.19 23.81 24.73 0.49 20.51 21.21

xB97 3.37 24.01 25.03 0.67 21.04 21.52

BLYP-GD2 0.14 22.45 21.48 22.56 0.85 2.04

DETS 1,DE122, andDE123are the barrier and reaction energy in respect to the bound complex (min1).DDE^{MP 2}_x , x 5 TS1, 1–2, and 1–3 are the errors in respect to MP2 results. [a] The results of TS2 listed in Support- ing Information Figure S2 were omitted due to the difficulty of converging PCM optimizations.

(8)

phase simulations. B97-GD2 performed strongly in the interal- coholamine interaction, and BLYP-GD2 produced an effective description of the CO2–alcoholamine interaction. Simple bimolecular reaction between CO₂ and MEA simulating using PCM model also confirmedxB97, xB97X, and xB97XD qualitatively reproduced the energetics at MP2 level.

Considering the different ranking in terms of root-mean- square error for the two comparisons (using MP2 or DFT optimized geometries), the inconsistent results of the tested XC functionals raised concerns about using the conventional pro- cedure of single-point ab initio electronic structure correction in combination with DFT optimized geometries, particularly for describing the CO₂binding by alcoholamines.

Acknowledgment

Computational resources were provided by the National Center for High Performance Computing.

Keywords: CO₂ capture

#

^density ^functional ^theory

#

alcoholamines

How to cite this article: H.-C. Li, J.-D. Chai, M.-K. Tsai Int. J.

Quantum Chem. 2014, 114, 805–812. DOI: 10.1002/qua.24670 Additional Supporting Information may be found in the online version of this article.

[1] J. T. Murphy, A. P. Jones, DOE/NETL Annual Report, 2009. Available at:

http://www.netl.doe.gov. Accessed on September 17, 2013.

[2] M. Caplow, J. Am. Chem. Soc. 1968, 90, 6795.

[3] H. Hikita, S. Asai, Y. Katsu, S. Ikuno, AIChe J. 1979, 25, 793.

[4] P. V. Danckwerts, Chem. Eng. Sci. 1979, 34, 443.

[5] E. Alper, Ind. Eng. Chem. Res. 1990, 29, 1725.

[6] S. H. Ali, Int. J. Chem. Kinet. 2005, 37, 391.

[7] D. E. Penny, T. J. Ritter, J. Chem. Soc., Faraday Trans. 1983, 79, 2103.

[8] H. Hikita, S. Asai, H. Ishikawa, M. Honda, Chem. Eng. J. 1977, 13, 7.

[9] J. -G. Shim, J. -H. Kim, Y. H. Jhon, J. Kim, K. -H. Cho, Ind. Eng. Chem.

Res. 2009, 48, 2172.

[10] E. F. da Silva, H. F. Svendsen, Ind. Eng. Chem. Res. 2004, 43, 3413.

[11] P. M. M. Blauwhoff, G. F. Versteeg, W. P. M. Vanswaaij, Chem. Eng. Sci.

1983, 38, 1411.

[12] B. Arstad, R. Blom, O. Swang, J. Phys. Chem. A 2007, 111, 1222.

[13] H. -B. Xie, Y. Zhou, Y. Zhang, J. K. Johnson, J. Phys. Chem. A 2010, 114, 11844.

[14] W. Conway, X. Wang, D. Fernandes, R. Burns, G. Lawrance, G. Puxty, M.

Maeder, J. Phys. Chem. A 2011, 115, 14340.

[15] C. A. Guido, F. Pietrucci, G. A. Gallet, W. Andreoni, J. Chem. Theory Com- put. 2013, 9, 28.

[16] H. Yamada, Y. Matsuzaki, T. Higashii, S. Kazama, J. Phys. Chem. A 2011, 115, 3079.

[17] S. Gangarapu, A. T. M. Marcelis, H. Zuilhof, Chemphyschem 2012, 13, 3973.

[18] P. Jackson, A. Beste, M. Attalla, Struct. Chem. 2011, 22, 537.

[19] K. R. Jorgensen, T. R. Cundari, A. K. Wilson, J. Phys. Chem. A 2012, 116, 10403.

[20] H. B. Xie, J. K. Johnson, R. J. Perry, S. Genovese, B. R. Wood, J. Phys.

Chem. A 2011, 115, 342.

[21] S. -W. Park, H. -B. Cho, I. -J. Sohn, H. Kumazawa, Sep. Sci. Technol. 2002, 37, 639.

[22] Y. Zhao, D. G. Truhlar, J. Chem. Theory Comput. 2006, 2, 1009.

[23] K. E. Riley, M. Pitonak, P. Jurecka, P. Hobza, Chem. Rev. 2010, 110, 5023.

[24] O. A. Vydrov, T. Van Voorhis, J. Chem. Theory Comput. 2012, 8, 1929.

[25] G. A. DiLabio, E. R. Johnson, A. Otero-de-la-Roza, Phys. Chem. Chem.

Phys. 2013, 15, 12821.

[26] P. Jurecka, J. Sponer, J. Cerny, P. Hobza, Phys. Chem. Chem. Phys. 2006, 8, 1985.

[27] J. Rezac, K. E. Riley, P. Hobza, J. Chem. Theory Comput. 2011, 7, 3466.

[28] P. Jurecˇka, P. Hobza, Chem. Phys. Lett. 2002, 365, 89.

[29] P. Hobza, J. !Sponer, J. Am. Chem. Soc. 2002, 124, 11802.

[30] I. Dabkowska, P. Jurecka, P. Hobza, J. Chem. Phys. 2005, 122, 204322.

[31] I. D. Mackie, G. A. DiLabio, J. Chem. Phys. 2011, 135, 134318.

[32] T. Yanai, D. P. Tew, N. C. Handy, Chem. Phys. Lett. 2004, 393, 51.

[33] O. A. Vydrov, G. E. Scuseria, J. Chem. Phys. 2006, 125, 234109.

[34] O. A. Vydrov, J. Heyd, A. V. Krukau, G. E. Scuseria, J. Chem. Phys. 2006, 125, 074106.

[35] O. A. Vydrov, G. E. Scuseria, J. P. Perdew, J. Chem. Phys. 2007, 126, 154109.

[36] J. D. Chai, M. Head-Gordon, J. Chem. Phys. 2008, 128, 084106.

[37] J. D. Chai, M. Head-Gordon, Phys. Chem. Chem. Phys. 2008, 10, 6615.

[38] A. D. Becke, Phys. Rev. A 1988, 38, 3098.

[39] C. T. Lee, W. T. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785.

[40] A. D. Becke, J. Chem. Phys. 1997, 107, 8554.

[41] A. D. Becke, J. Chem. Phys. 1993, 98, 5648.

[42] S. Grimme, J. Comput. Chem. 2006, 27, 1787.

[43] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104.

[44] Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2005, 109, 5656.

[45] Y. Zhao, D. G. Truhlar, Theor. Chem. Acc. 2008, 120, 215.

[48] Y. Zhao, D. G. Truhlar, J. Chem. Phys. 2006, 125, 194101.

[49] T. Schwabe, S. Grimme, Phys. Chem. Chem. Phys. 2007, 9, 3397.

[50] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R.

Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H.

Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G.

Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J.

Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T.

Vreven, J. J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J.

Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J.

Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J.

Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O.

Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J.

Dannenberg, S. Dapprich, A. D. Daniels, €O. Farkas, J. B. Foresman, J. V.

Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09; Gaussian, Inc.: Wallingford CT, 2009.

[51] M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J.

Van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus, W. de Jong, Comput. Phys. Commun. 2010, 181, 1477.

Received: 17 September 2013 Revised: 18 February 2014 Accepted: 26 February 2014 Published online 19 March 2014