Quantum Chemical Estimation of Acetone Physisorption on Graphene

Publication Date (Web): April 4, 2017. Copyright © 2017 American Chemical Society. *E-mail: [email protected]. Cite this:J. Phys. Chem. C 121...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Quantum Chemical Estimation of Acetone Physisorption on Graphene Using Combined Basis Set and Size Extrapolation Schemes Yoshifumi Nishimura,†,‡ Takao Tsuneda,§ Takeshi Sato,∥ Michio Katouda,⊥ and Stephan Irle*,†,# †

Department of Chemistry, Graduate School of Science, and #Institute of Transformative Bio-Molecules (WPI-ITbM), Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan § Fuel Cell Nanomaterials Center, University of Yamanashi, Kofu 400-0021, Japan ∥ Photon Science Center, University of Tokyo, Tokyo 113-8656, Japan ⊥ Computational Molecular Science Research Team, RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan S Supporting Information *

ABSTRACT: The physisorption of an acetone molecule on hexagonal graphene nanoflakes with increasing size has been investigated using a variety of quantum chemical methods capable of describing weak intermolecular interactions: coupledcluster theory (CCSD and CCSD(T)), second-order Møller−Plesset perturbation theory with and without spin-component scaling (SCS-MP2 or standard MP2), longrange corrected density functional theory combined with a van der Waals functional (LC-BOP+ALL), meta-generalized gradient approximation functionals (M06-2X and M05-2X), and the dispersion augmented self-consistent-charge density functional tight-binding (SCC-DFTB-D) method. Our benchmark results for model systems as large as dicircumcoronene C96H24 have confirmed the suitability of the SCS-MP2 method for this specific system and the satisfactory performance of the computationally much more economical semiempirical SCC-DFTB-D method. The latter delivers a qualitatively accurate description of physisorption for flakes containing more than 800 carbon atoms. We predict accurate interaction energies of acetone with an infinitely large, defect-free graphene monolayer by combining extrapolation approaches for both increasing ab initio basis sets and graphene flake size in a two-dimensional extrapolation scheme.

1. INTRODUCTION Nowadays, carbon-based nanostructures such as fullerenes, carbon nanotubes (CNTs), and graphenes are widely recognized as the key materials of technological applications in various research fields. One of the most active subjects is to elucidate the adsorption of small molecules with nanoscale carbonaceous surfaces, related to the development of gas storage,1,2 chemical sensors,3−5 or drug delivery systems.6−8 It is thus of crucial importance to obtain accurate estimates for such weak intermolecular interactions from quantum chemical calculations. However, due to the large size of molecular models required for this kind of system, highly accurate ab initio electronic structure methods including electron correlation [e.g., Møller−Plesset perturbation theory (MP)9 and coupled-cluster theory (CC)10,11] are computationally too demanding. They can only be practically employed for small model systems of graphite, such as the benzene dimer or benzene−methane complex,12−14 to accurately describe π−π or CH/π interactions in terms of binding energies and equilibrium intermolecular distances. Application of second-order MP (MP2)15−17 or quadratic configuration interaction16 methods to larger graphene dimers have been reported in calculations highlighting the performance of parallel program implementa© 2017 American Chemical Society

tions, but there are no scientific applications of such computationally highly demanding studies. First-principles density functional theory (DFT)18,19 is computationally far less demanding than the above-mentioned ab initio wave function based methods. However, traditional exchange-correlation functionals generally perform poorly in the description of weak interactions, specifically in the case of van der Waals attractions such as dispersion interactions.20,21 Several approaches have been proposed to overcome this severe drawback. Adding a dispersion correction term to the Kohn−Sham total energy (vdW-DFT) has enjoyed much popularity, since it was shown that its use with nonlocal DFT functionals has produced encouraging results for many systems.22−25 Another strategy was chosen in a series of DFT-Dn (up to n = 3) methods, developed by Grimme and coworkers,26−28 where the use of a London dispersion term with precomputed atom dependent coefficients and scaling factors for a variety of functionals is able to deliver good accuracy with negligible additional computational effort. Zhao and Truhlar Received: December 28, 2016 Revised: March 27, 2017 Published: April 4, 2017 8999

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

packages (including a private development version for LC-BOP +ALL calculations). Ahlrichs’ SVP,58 TZVPP,59 QZVPP,60 and corresponding auxiliary basis sets for the resolution-of-identity (RI) approach61,62 were applied to perform geometry optimizations and to evaluate intermolecular interaction energies. The equilibrium geometries were optimized by imposing the following symmetry constraints given in parentheses at the RI-MP2/SVP level of theory63,64 under default convergence criteria, as implemented in Turbomole: acetone (C2v), hexagonal graphene flakes (C6n2H6n, n = 1−4) (D6h), and acetone-graphene flake complexes (Cs or C2v). For benzene (C6H6, n = 1) and coronene (C24H12, n = 2) model systems, RI-MP2/QZVPP geometry optimizations were performed additionally. Higher-level, single-point energies for these optimized geometries were computed using the following methods: ab initio coupled-cluster singles and doubles (CCSD), including perturbative triples [CCSD(T)], RI-MP2, and its spin-component-scaled version (RI-SCS-MP2),65 DFT using two meta-GGA functionals (M05-2X29 and M06-2X30) and the so-called LC-BOP+ALL functional.24,66,67 “LC” denotes the long-range correction for exchange functionals.68 The LC scheme was applied to the Becke 1988 exchange functional69 plus one-parameter progressive (OP) correlation functional70 (LC-BOP). For the parameter μ controlling the ratio between short- and long-range exchange, we used two different values: μ = 0.330 which is accepted to be optimum for response properties including excitation energies,71 and μ = 0.470 which was found to be more appropriate for ground state property calculations.72 “+ALL” represents the AndersonLangreth-Lundqvist (ALL) dispersion functional.22 The applicability of LC-BOP+ALL to a wide range of weakly bound systems is documented elsewhere.73 The SCC-DFTB-D calculations were carried out with the Paderborn DFTB/“Dylax” code and the standard mio parameter set.35 The set of constants for the London dispersion term was taken from ref 36. For geometry optimizations, the conjugated gradient method was adopted with convergence thresholds of 1 × 10−5 a.u. and 1 × 10−12 a.u. for maximum force and convergence of the SCC-DFTB Mulliken charge, respectively. The hexagonal graphene flake size was extended up to n = 12 (C864H72). For comparison, the interaction energies of an acetone molecule with an infinite, one-layer graphite sheet were also computed using periodic boundary condition (PBC) including 160 carbon atoms in a unit cell and used as the infinite size limit of the graphene flake model. The intermolecular interaction energy (INT) of the physisorbed acetone on a graphene flake is defined as follows:

have designed multiple meta generalized-gradient-approximation (GGA) functionals such as M05-2X29 and M06-2X,30 in which noncovalent interaction is accounted for by the balanced incorporation of spin densities, kinetic energy densities, and other variables along with careful parametrization. But although current dispersion-augmented DFT methods are capable to cope with molecular systems containing several hundreds of atoms, it remains still difficult to converge basis set and size effects in the study of small molecule/nanomaterials complexes. Approximate semiempirical methods were developed to allow the quantum chemical treatment of systems containing hundreds or a few thousands of atoms. The most commonly used semiempirical techniques based on molecular orbital theory are MNDO,31 AM1,32 and PM3,33 and the selfconsistent-charge density-functional tight-binding (SCCDFTB) method.34,35 The latter is firmly rooted in the formalism of DFT. Several sensible, yet very efficient approximations enable this method to reach computational speeds that are 2−3 orders of magnitude faster than traditional DFT methods with qualitative and often quantitative accuracy. SCC-DFTB inherited the problems from traditional DFT regarding the prediction of van der Waals interactions, but it was shown that this limitation can be overcome empirically by two different schemes; one is the R−6 function attenuated for short distances with damping function where C6 coefficient is calculated using atomic hybridization dependent polarizabilities;36 the other consists of Lennard-Jones type functional form in which the short-range interaction is given by polynomial function with parametrized coefficients.37 Earlier publications have demonstrated the suitability of the dispersion-augmented SCC-DFTB (SCC-DFTB-D) method toward molecular physisorption on carbon nanomaterials for relatively simple molecules: hydrogen,37 water and its clusters,38,39 ammonia,40 methanol,41 and others.42,43 In this work, we focus on the theoretical investigation of the interaction between a single acetone molecule and a graphite (0001) surface, modeled using small to large hexagonal graphene flakes. Acetone is the simplest organic compound of ketones and typically used as volatile organic solvent. The adsorption of acetone on a carbonaceous surface has been extensively studied experimentally,44−48 however, fundamental knowledge of its mode of interaction is still lacking. Only few theoretical studies have tried to estimate the interaction energy of acetone on single-walled CNTs (SWCNTs), mainly using standard DFT.49,50 We have previously studied acetone− coronene interactions and compared MP2 energies and geometries with those from SCC-DFTB-D,51−53 yet, coronene is too small a model system to fully describe the interaction of acetone on an infinitely large graphene sheet. The present benchmark study aims to obtain deep insights into the structural and energetic properties between acetone and a series of hexagonal graphene flakes with systematically increasing size to assess the validity of representative quantum chemical methods for the description of their intermolecular interactions. In the present study, we report an accurate estimate of the interaction energy using a novel twodimensional extrapolation scheme that addresses both the convergence behavior of ab initio calculations with increasing basis set size and that of the graphene model system size.

INT = E(complex opt) − E(acetoneopt) − E(graphene flakeopt) + BSSE − DEF

(1)

In eq 1, E(complexopt), E(acetoneopt), and E(graphene flakeopt) denote the calculated electronic energy of acetone− graphene flake complex, acetone molecule, and graphene flake at the respective fully optimized geometry in vacuum. BSSE is an abbreviation of basis set superposition error, which is corrected using the Boys−Bernardi counterpoise (CP) method,74 except for the SCC-DFTB-D method, for which CP methods cannot be applied. DEF stands for the magnitude of structural deformation of each constituent by forming the complex in terms of electronic energy E given by

2. COMPUTATIONAL DETAILS All ab initio and DFT calculations were performed using the Turbomole,54,55 Gaussian 03,56 and Gaussian 0957 program 9000

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Figure 1. RI-MP2/SVP optimized structures of acetone-dicircumcoronene (C96H24) complexes. The calculated intermolecular distance r is indicated by a dashed line.

atoms. See Tables S1 and S2 in the Supporting Information for the comparison of optimized bond lengths and bond angles of monomers. The graphene flakes in the acetone-flake complexes were slightly warped due to the attraction between acetone and the nearest carbon atoms of the flake at all levels of theory. The PP complexes tend to induce the largest deformations compared to the UP or DP types, which is a result of the larger contact surface area. The degree of distortion becomes larger as the size of the graphene flakes (n) increases, which is reflected in their deformation energy (DEF) as well. The decomposition of DEF into acetone and graphene flake components shows that the latter dominates with a contribution greater than 90%. While a DEF of a few kJ/mol is computed for the SCC-DFTB-D and RI-MP2/QZVPP optimized structures, the RI-MP2/SVP optimized structures exhibit a DEF larger than 10 kJ/mol. This is a consequence of the BSSE since the RI-MP2/SVP gradients were not CP-corrected. Although the usual supermolecular approach discusses the total interaction energy ΔE = DEF + INT, the larger DEFs for larger graphene flakes in RIMP2/SVP geometries prevented a smooth convergence behavior of ΔE with respect to graphene size. For enabling the analysis and discussion of all benchmark results in the same manner, we therefore decided to use the quantity INT, as defined in eq 1. Figure 2 displays the evolution of the intermolecular distance r between acetone and the graphene flake defined in Figure 1 as a function of graphene size n. Note that the surface normal of flakes given in Figure 1 was defined as the normal vector of a plane specified in Cartesian coordinates of carbon atoms with the first three short distances from a reference atom of acetone molecule (oxygen atom for PP and DP complexes and hydrogen atom for UP complexes). We find that the distance becomes shorter as n increases, and the decrease in r is largest for the size change from n = 1 to 2. An exception is found for the case of the RI-MP2 optimized PP2 geometries, caused by changes of the angle between the acetone CCCO plane and the graphene flake surface in going from n = 1 to 2. This angle is 28.5 and 20.0 degrees for SVP and QZVPP in the case of n = 1,

DEF = {E(acetonecomplex ) − E(acetoneopt)} + {E(graphene flakecomplex ) − E(graphene flakeopt)} (2)

where the subscript “complex” corresponds to the geometry of monomer in acetone−graphene flake complex. We will later elaborate on the role of DEF. To estimate INT at the complete basis set (CBS) limit, different extrapolation schemes are used in regard to Ahlrichs’ basis sets. For the wave function based methods, the Hartree−Fock (HF), and electron correlation energies are extrapolated separately following the procedure by Feller75,76 and Helgaker and co-workers,77 respectively. For the DFT energies, we applied the formula proposed by Kahn and Bruice.78 For the estimation of CCSD(T) energies with a large basis set, the technique reported by Tsuzuki and co-workers79 was additionally employed.

3. RESULTS AND DISCUSSION 3.1. Optimized Structures of Physisorption Complexes. Following our previous work on acetone−coronene interactions,51 we take a total of six possible acetone−graphene flake van der Waals complex geometries into account. The acetone molecule is placed above the central hexagonal ring of graphene flakes where the CCCO plane is roughly parallel or the CO axis is perpendicular with respect to the graphite (0001) surface. The former is referred to as planar-parallel (PP) structure possessing Cs symmetry. The latter has C2v symmetry and depends on the direction of the CO axis; if the oxygen atom is pointing away from the central hexagonal ring, this type of complex is called up-perpendicular (UP); if it points toward the central ring, we name the complex downperpendicular (DP). Furthermore, each complex type is distinguished with labels “1” and “2”, which specifies the position of methyl groups relative to the central hexagonal ring. As an example, the RI-MP2/SVP optimized structures of acetone-dicircumcoronene complexes are shown in Figure 1. Additionally, we also optimized the smaller complexes using RIMP2/QZVPP. SCC-DFTB-D geometry optimizations were carried out for complexes with flakes containing more than 800 9001

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Figure 2). As the flakes become larger, the relative orientation of the acetone molecule with respect to the two possible orientations labeled “1” and “2” on the hexagonal flakes does not make a significant difference. The converged intermolecular distances are shorter than the sum of van der Waals radii,80 indicating the existence of π−π and/or CH/π interactions between the two partner molecules in the complexes. In summary, the computationally much more economical SCC-DFTB-D method is capable to describe the important structural information on acetone-graphene flake physisorption complexes without severe deterioration in comparison to RIMP2, which is one of the most common methods for the qualitative treatment of van der Waals interactions for small to medium-sized molecules. 3.2. Benchmark of Acetone−Benzene Interaction Energies. Here we assess the performance of various quantum chemical methods and their basis set dependence for the description of interaction energies (INT). Figure 3 shows CC, MP2, DFT, and SCC-DFTB-D single-point INT energies at RIMP2/QZVPP optimized acetone−benzene complexes using increasingly larger basis sets up to respectively extrapolated CBS limits. The corresponding results of RI-MP2/SVP and SCC-DFTB-D optimized geometries are presented in Figure S1 in the Supporting Information. It should be mentioned that the correlation energies of RI-MP2 and RI-SCS-MP2 at the CBS limit were extrapolated based on TZVPP and QZVPP energies, while those of CCSD and CCSD(T) were derived using only SVP and TZVPP basis sets, due to computational cost. All employed methods lead to the following order of complex stability: PP > UP > DP, which we had already reported earlier.51 The behavior of INT variation is virtually the same with respect to the relative orientation of acetone to the hexagonal ring, labeled “1” and “2”. It is clear that the value of INT predicted by the wave function-based methods (solid lines in Figure 3) strongly depends on the employed basis set. We see that the “error” of RI-SCS-MP2, defined as the difference relative to CCSD(T) in combination with the respective basis set, is smaller than that of CCSD, up to TZVPP. The RI-MP2 method tends to overestimate INT especially for PP and UP structures, which is a well-documented overstabilization of π−π

Figure 2. Intermolecular distance r (Å) between acetone and graphene flake surface with respect to flake size n. The results of SCC-DFTB-D PBC calculations are displayed as black dotted line together with its value.

and significantly reduced to 13.7 and 10.5 degrees for SVP and QZVPP in the case of n = 2. We note that the SCC-DFTB-D method gives much smaller values (2.9 and 0.4 degrees for n = 1 and 2, respectively). The complete list of this angle for all optimized PP structures is given in Table S3 in the Supporting Information. For n = 3 and 4, the deviation in the prediction of r between the SCC-DFTB-D and the RI-MP2 methods is less than 0.3 Å for PP, and less than 0.1 Å for UP and DP structures. The SCC-DFTB-D distance essentially levels off after n ≥ 5, and a similar value of r is obtained using periodic boundary condition (PBC) calculations (black dotted line in

Figure 3. Calculated INT (kJ/mol) of various methods at RI-MP2/QZVPP optimized acetone−benzene model systems. 9002

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Figure 4. Calculated BSSE (kJ/mol) of various methods at RI-MP2/QZVPP optimized acetone-benzene model systems.

Table 1. Calculated and Estimated BSSE-Corrected-INT (kJ/mol) at the CCSD(T) Level of RI-MP2/QZVPP Optimized Acetone−Benzene Model Systems basis set TZVPP

QZVPP

CBS with eq 5

a

extrapolation

PP1

PP2

UP1

UP2

DP1

DP2

calcd eq 3 large-basis, TZVPP; small-basis, SVP eq 4 large-basis, TZVPP; small-basis, SVP Δ[eq 3]a Δ[eq 4]a eq 3 large-basis, QZVPP; small-basis, TZVPP eq 4 large-basis, QZVPP; small-basis, TZVPP eq 5 SVP (calcd) TZVPP (calcd) SVP (calcd) TZVPP (calcd) TZVPP (calcd) QZVPP [eq 3] TZVPP (calcd) QZVPP [eq 4] TZVPP (calcd) QZVPP [eq 5]

−11.92 −13.25 −11.35 −1.33 0.57 −14.21 −13.88 −13.76 −19.25 −15.06 −15.41 −18.39

−12.21 −13.56 −11.62 −1.35 0.59 −14.48 −14.16 −14.11 −20.34 −15.24 −15.69 −19.41

−9.12 −10.02 −8.53 −0.90 0.59 −10.56 −10.36 −10.55 −14.20 −11.16 −11.57 −13.62

−9.33 −10.19 −8.74 −0.85 0.59 −10.68 −10.49 −10.72 −14.36 −11.25 −11.66 −13.80

3.24 2.46 3.34 −0.78 0.10 1.98 2.23 2.06 −2.14 1.26 0.96 −1.63

3.24 2.46 3.34 −0.78 0.10 1.98 2.23 2.06 −2.14 1.26 0.95 −1.63

Δ = INT(extrapolated) − INT(calculated).

and CH−π interactions by this method.14,79,81,82 The DFT methods (dashed lines in Figure 3) do not exhibit significant changes of INT with increased basis set size as expected, thus we only focus on the comparison of first principle results with ab initio CBS limits. Concerning the LC-BOP+ALL method, the difference of INT between μ = 0.330 and 0.470 is less than 3 kJ/mol. The M06-2X functional shows closer agreement with the CBS limit INT of RI-SCS-MP2 for PP and UP structures, yet the interaction of DP structures is weakened by approximately 4 kJ/mol. The M05-2X method underestimates INT by approximately 1−3 kJ/mol with respect to INT of M06-2X. Interestingly, the CP-uncorrected SCC-DFTB-D INT is comparable to the high-level quantum chemical methods. For example, the largest difference of the SCC-DFTB-D INT relative to the INT of RI-SCS-MP2/CBS is 2.99 kJ/mol in the UP1 structure. This finding encourages us to apply the SCCDFTB-D method to the study of the graphene model size dependence of INT using larger model systems. Figure 4 displays the BSSE according to CP calculations for the various quantum chemical methods for the RI-MP2/ QZVPP optimized acetone benzene complexes as a function of basis set size. See Figure S2 in the Supporting Information for the same results calculated using other optimized geometries.

As expected, wave function-based methods have a twice larger BSSE than DFT results for PP and UP structures; on the other hand, the BSSE of DP structures is similar between wave function and DFT methods. Since the contact surface area is reduced in this particular conformation, the perpendicular orientations experience smaller BSSEs than PP structures. As can be seen from BSSE values larger than 5 kJ/mol in Figure 4, the SVP basis set seems to be too small for wave function-based methods to quantitatively predict physisorption energies even after CP corrections. Therefore, the quality of CCSD and CCSD(T) results at the CBS limit in Figure 3 is somewhat doubtful, because our CBS extrapolation in this case was based only on SVP and TZVPP data. To check the CBS limit more carefully in this case, INT at the CCSD(T) level along with large basis sets are estimated using different approximation schemes. The results of RI-MP2/QZVPP model systems are shown in Table 1. See Tables S4 and S5 in the Supporting Information for the corresponding results with other optimized structures. In previous theoretical works,14,79 an MP2-based extrapolation of the influence of increased basis set on energies obtained using CCSD(T) with smaller basis sets was applied to obtain estimates for CCSD(T) energies with larger basis sets. In such an extrapolation scheme, 9003

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Table 2. Dispersion (D) and Nondispersion (non-D) Energy Contributions to BSSE-Corrected-INT (kJ/mol) of RI-MP2/ QZVPP Optimized Acetone−Benzene Model Systems with LC-BOP+ALL and SCC-DFTB-D Methods LC-BOP+ALL basis set SVP

μ

energy

PP1

PP2

UP1

UP2

DP1

DP2

0.470

INT non-D D INT non-D D INT non-D D INT non-D D INT non-D D INT non-D D

−14.69 0.34 −15.03 −16.56 −1.41 −15.15 −17.00 −0.36 −16.65 −19.26 −2.30 −16.96 −17.00 −0.50 −16.50 −19.38 −2.45 −16.94

−14.71 0.30 −15.01 −16.58 −1.41 −15.16 −17.34 −0.60 −16.75 −19.54 −2.55 −17.00 −17.36 −0.72 −16.65 −19.54 −2.68 −16.86

−12.59 0.50 −13.08 −14.07 −0.82 −13.26 −13.69 0.46 −14.15 −15.53 −1.02 −14.51 −13.51 0.47 −13.98 −15.27 −1.00 −14.27

−12.85 0.11 −12.97 −14.21 −1.11 −13.10 −14.03 0.08 −14.11 −15.60 −1.32 −14.29 −13.86 0.09 −13.94 −15.58 −1.31 −14.28

3.19 10.91 −7.72 1.76 9.47 −7.71 0.50 10.05 −9.55 −1.34 8.41 −9.75 0.07 10.00 −9.93 −2.07 8.33 −10.41

3.21 10.91 −7.70 1.77 9.47 −7.69 0.47 10.05 −9.57 −1.38 8.41 −9.79 0.08 10.00 −9.92 −2.09 8.33 −10.42

0.330

TZVPP

0.470

0.330

QZVPP

0.470

0.330

energy INT non-D D

PP1 −15.94 −3.36 −12.58

SCC-DFTB-D PP2 UP1 −15.51 −13.30 −3.34 −2.52 −12.17 −10.78

UP2 −13.19 −2.54 −10.65

DP1 0.22 5.09 −4.87

DP2 0.22 5.09 −4.87

Figure 5. Calculated INT (kJ/mol) with TZVPP basis set of RI-MP2/QZVPP optimized acetone−coronene model systems. Note that CCSD(T) ≈ ECCSD(T) + ERI‑SCS‑MP2 − ERI‑SCS‑MP2 . results are estimated using extrapolation:ECCSD(T) TZVPP SVP TZVPP SVP

TZVPP calculation reveals that the error is slightly reduced in eq 4 for all orientations and therefore preferable over eq 3 in this study. For the QZVPP energies, we additionally examined the extrapolation scheme as follows

we have applied both RI-MP2 as well as RI-SCS-MP2 based extrapolations, that is, RI‐MP2 RI‐MP2 CCSD(T) CCSD(T) E large ‐basis ≈ Esmall‐basis + E large‐basis − Esmall‐basis CCSD(T) E large ‐basis



CCSD(T) Esmall ‐basis

+

RI‐SCS‐MP2 E large ‐basis



RI‐SCS‐MP2 Esmall ‐basis

(3)

HF corr CCSD(T) E large ‐basis ≈ E large‐basis + E large‐basis

(4)

(5)

where large-basis is QZVPP, the Hartree−Fock (HF) energy EHF large‑basis is explicitly calculated, and the corresponding electron correlation energy Ecorr large‑basis is estimated following the inverse relationship of SVP and TZVPP electron correlation energies

where the relationship between large-basis and small-basis shows the consecutive cardinal number. The comparison of both basis set extrapolation schemes with explicit CCSD(T)/ 9004

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Figure 6. Calculated INT (kJ/mol) with TZVPP basis set of RI-MP2/SVP optimized acetone−graphene flake model systems. The fitted curves using eq 6 are drawn with solid lines for the wave function-based methods, dashed lines for DFT methods, and dotted line for SCC-DFTB-D method.

3.3. Benchmark of Acetone−Coronene Interaction Energies. Our previous study51 assessed the validity of the SCC-DFTB-D method by making a comparison for the acetone−coronene (n = 2) interaction energies, based on RIMP2 results alone. Here we reinvestigate INT, employing various quantum chemical methods with the TZVPP basis set. As shown in Figure 5 for RI-MP2/QZVPP optimized model systems, irrespective of the method, the coronene graphene flake experiences 10−20 kJ/mol stronger interactions with acetone than benzene. See Figure S3 in the Supporting Information for the corresponding results with other optimized geometries. In the following, we summarize the main observations from this data. Similar to acetone−benzene systems, RI-SCS-MP2 results are in excellent agreement with extrapolated CCSD(T) values, and they predict the order of stability as PP ≈ UP > DP. INT of RI-MP2 is almost twice as large as those of the other wave function-based methods for PP and DP structures. Interestingly, it is found that LC-BOP+ALL (μ = 0.470) predicts that the UP structures are more stable than PP structures, although the deviation of INT from CCSD(T) is only 2−3 kJ/mol for PP structures. This indicates that a larger value of μ emphasizes CH/π interactions. For structures with acetone aligned in perpendicular to the graphene plane, the difference between LC-BOP+ALL (μ = 0.330) and RI-MP2 is less than 3 kJ/mol. INT of DP structures is much underestimated using M06-2X and M05-2X methods. M05-2X produces about 10 kJ/mol smaller INT than M06-2X for PP structures. SCC-DFTB-D overestimates INT by 7−8 kJ/mol for PP and UP structures and by 4−5 kJ/mol for DP structures relative to CCSD(T) results. 3.4. Graphene Model System Size Dependence of Acetone−Graphene Interaction Energies. For model systems larger than acetone−circumcoronene (n = 3), even CCSD(T)/SVP calculations were prohibitively expensive. However, the above benchmarks with smaller flakes have demonstrated CCSD(T)-like performance of the RI-SCS-MP2 method for acetone−graphene flake interactions, and our

with the third power of cardinal number as proposed by Helgaker and co-workers.77 The deviation of predicted INT with QZVPP basis set is less than 0.5 kJ/mol among the three schemes, thus, eq 5 is acceptable for the prediction of electron correlation energies obtained with basis sets of the next cardinal number. The values of INT at the CBS limit are determined from the summation of CBS limit HF and electron correlation energies. The former is calculated with the exponential fitting of SVP, TZVPP, and QZVPP HF energies with respect to cardinal number, as argued by Feller. 75,76 The straightforward extrapolation model by Helgaker and co-workers77 is again used for the latter quantity. We note that the fitting of a sequence of basis set levels TZVPP and QZVPP leads to a doubly extrapolated value. The results definitely depend on the quality of basis set and extrapolation schemes. Considering the close asymptotic behavior of electron correlation energies for basis sets larger than triple-ζ levels77 and the observation of reduced errors for prediction of CCSD(T)/TZVPP energies, the reliable extrapolation to CBS is most likely achieved by combining TZVPP (calculated) and QZVPP (eq 4) for the acetone−graphene flake model systems. Table 2 reports two separable energy contributions to INT from the LC-BOP+ALL and SCC-DFTB-D methods, namely, the DFT and DFTB electronic energy (non-D) and dispersion energy (D) contribution. See Tables S6 and S7 in the Supporting Information for the same results computed using other geometries. It is apparent that the dispersion energy governs the main part of the interaction energy. The triple-ζ quality basis set seems to be sufficiently large to evaluate both contributions accurately with the LC-BOP+ALL method. The SCC-DFTB-D empirical dispersion term is about 5 kJ/mol smaller than the value calculated using the ALL functional. Both methods show roughly 10 kJ/mol larger nondispersion contribution of DP structures than PP and UP, likely a manifestation of the stronger Pauli repulsion between oxygen lone-pairs and the graphene flake π-system relative to other orientations. 9005

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C

Figure 7. Two dimensional contour plot of RI-SCS-MP2 INT for RI-MP2/SVP optimized acetone−graphene flake model systems. The contour lines are highlighted by a magenta color and pointed arrow every 10 kJ/mol for easier reading.

computational resources allowed to calculate RI-SCS-MP2/ TZVPP single point energies for system sizes up to n = 4. RI-MP2/TZVPP-quality values of INT at RI-MP2/SVP optimized geometries for hexagonal graphene flake sizes up to n = 4 are shown in Figure 6. As already mentioned, RI-SCS-MP2 yields a relative stability order of PP ≈ UP > DP for n = 2, however, increasing the graphene size results in a slightly different order, PP > UP > DP. The energy differences between the relative acetone orientations to the hexagonal graphene flake, denoted by “1” and “2” structures, are only on the order of 2−3 kJ/mol, which is even smaller than for the smaller graphene flakes. The increment of INT from n = 3 to 4 is about 2 kJ/mol, thus, the effect of hydrogen atoms at the rim of graphene flakes is pretty much negligible at this size. The overbinding of PP structures with RI-MP2 reaches about 20 kJ/ mol relative to RI-SCS-MP2. As before, LC-BOP+ALL (μ = 0.470) predicts the order of stability PP ≈ UP > DP for large graphene flakes; on the other hand, μ = 0.330 predicts PP > UP

> DP similar to MP2 calculations. M05-2X and M06-2X functionals indicate very weak binding of DP structures even if the flake size is increased beyond dicircumcoronene. The SCCDFTB-D results are close to the behavior of LC-BOP+ALL (μ = 0.330) and reproduce the associated order of stability qualitatively. For extrapolation of INT to n → ∞ (corresponding to defect-free monolayer graphene), we use the following exponential form: INT(n) = a0 +

1 exp(a1n + a 2)

(6)

where n is size of graphene flake (C6n2H6n) and ai (i = 0−2) are parameters determined by solving a least-squares fitting problem with INT at n = 1−4. Similar fitting techniques of interaction energies with graphene flake sizes have been applied for MP2 binding energies of water−graphene systems by Feller and Jordan.83 The shape of the fitting curves are plotted by 9006

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C solid, dashed, and dotted lines in Figure 6 for RI-MP2/SVP optimized model systems. See Figure S4 in the Supporting Information for the corresponding data of SCC-DFTB-D optimized structures. The fitted parameters ai are listed in Tables S8 and S9 in the Supporting Information. The correlation coefficient is generally larger than 0.99, indicating a reasonable asymptotic behavior. It is found that a small number of fitting data points is not a significant problem to estimate the interaction between acetone and pristine singlelayer graphite sheet using this extrapolation technique. We confirmed using eq 6 that at n = 1−4 shows less than 2 kJ/mol deviation in comparison to the interaction energies for model systems with PBC at the SCC-DFTB-D level. The details of this analysis are shown in the Supporting Information. 3.5. Basis Set and Graphene Flake Size Dual Extrapolation for Accurate Acetone−Graphene Interactions. The results presented above verified that (i) RI-SCSMP2 is capable of describing acetone−graphene flake interactions similar to CCSD(T) quality and (ii) the graphene size extrapolation allows to predict the interaction of acetone with an infinitely large graphene sheet. There is little doubt that the RI-SCS-MP2/TZVPP//RI-MP2/SVP energies and the corresponding flake size extrapolation, shown as the fitted blue solid curve in Figure 6, is already fairly accurate, but as the previous discussion of the acetone-coronene system showed, it would be preferable to have RI-SCS-MP2 INT data obtained with a larger size of basis set. For this purpose we have performed a two-dimensional (2D) fit in terms of the size of the basis set (measured by cardinal numbers, 3 for TZVPP, 4 for QZVPP, 5 for a hypothetical polarized quintuple zeta basis set, and so on) and graphene flake size n. Regarding the basis set, the extrapolation of correlation energy77 is performed by taking TZVPP and QZVPP results for n = 1 and 2, and SVP and TZVPP results for n = 3 and 4. For the model systems with n = 1 and 2, it is possible to determine the HF energy for each basis set using an exponential formula.75,76 For n = 3 and 4, we used the HF/TZVPP energy as reference value for the Hartree−Fock energy, since HF/QZVPP calculations were not feasible for graphene flakes of this size. Note that the error arising from this simplification is almost negligible since the difference of INT between HF/TZVPP and HF/QZVPP for n = 2 is less than 0.3 kJ/mol in all orientations. Figure 7 exhibits the 2D maps of RI-SCS-MP2 INT for RI-MP2/SVP optimized geometries where contour lines are plotted at intervals of 2 kJ/ mol. These plots suggest that TZVPP basis sets underestimate INT by at most 8 kJ/mol relative to sufficiently large basis sets. From this 2D-extrapolation, we can finally obtain an accurate estimate of INT(n → ∞) at the RI-SCS-MP2 in conjunction with basis set of cardinal number 12 as follows: −42 to −43 kJ/ mol for PP, −30 kJ/mol for UP, and −20 kJ/mol for DP structures. The comparison of predicted interaction energies of acetone with pristine single-layer graphite sheet using extrapolation techniques with selected available experimental measurements44−46 is tabulated in Table 3. On the basis of acetone− benzene benchmarks, we concluded that graphene size extrapolations with results computed from TZVPP basis set are reliable enough for DFT methods. Even though our theoretical data do not include zero-point vibrational energies or thermal corrections and neglect contributions due to structural deformation and possible interactions like physisorption of acetone cluster assemblies on a carbonaceous surface, the RI-SCS-MP2, LC-BOP+ALL, and SCC-DFTB-D

Table 3. Theoretical Interaction Energies of Acetone with Defect-Free Monolayer Graphene (Extrapolation of INT to n → ∞ According to Eq 6) via Extrapolation Schemes for RI-MP2/SVP Optimized Geometries (Upper Table) and Experimentally Available Data (Lower Table) Theory |INT| (kJ/mol) a

RI-SCS-MP2 RI-MP2a LC-BOP+ALLb (μ = 0.470) LC-BOP+ALLb (μ = 0.330) M06-2Xb M05-2Xb SCC-DFTB-Dc

ref

energy (kJ/mol)

44

31 ± 2

45

31.8

46

40

PP1

PP2

UP1

UP2

DP1

DP2

43.2 66.9 28.8

42.0 65.2 27.5

30.3 45.5 27.4

29.9 44.6 27.4

20.0 30.4 17.4

20.5 30.7 18.5

35.1

34.2

30.9

30.9

20.7

20.8

28.3 24.8 19.4 17.1 14.7 15.0 38.6 38.5 31.6 Experiments

18.1 14.1 32.0

6.5 6.4 17.5

7.0 7.8 17.5

comment desorption of acetone multilayer from graphite using temperature-programmed desorption technique heat of adsorption of acetone on graphitized carbon black surface using gas chromatography heat of adsorption of acetone on activated carbon on the basis of adsorption isotherms

2D-extrapolation at the basis set with cardinal number 12 and n → ∞. bGraphene flake size extrapolation using TZVPP basis set for n = 1−4. cGraphene flake size extrapolation for n = 1−4. a

INTs of PP and UP orientations are in a reasonable order of magnitude. RI-MP2 tends to overestimate by 15−20 kJ/mol, and meta-GGA functionals produce weak interactions by more than 10 kJ/mol in comparison to the RI-SCS-MP2 values.

4. CONCLUSIONS We have theoretically investigated the interactions of the acetone molecule with D6h hexagonal graphene flakes, including six different orientations of the physisorbed molecule. The interaction clearly originates from dispersion interactions since electrostatic contributions are 1 order of magnitude smaller at best. It is found that the RI-SCS-MP2 yields interaction energies comparable to CCSD(T) in benchmark systems for which we could perform the expensive coupled cluster calculations. The difference between RI-SCS-MP2 and CCSD(T) is only a few kJ/mol in conjunction with the triple-ζ quality of basis set. The performance of dispersion-augmented DFT methods is also assessed. The LC-BOP+ALL functional overestimates INT by several kJ/mol relative to CCSD(T) and RI-SCS-MP2 at the same basis set size, and the μ-value of LC-BOP of 0.470 yields a slightly different order of relative stability for acetone−coronene conformations. In particular, we find that a larger μ-value emphasizes CH/π interactions. The meta-GGA functionals underestimate the interaction energy of DP-type structures considerably and has a nonconsistent basis set dependence (smaller basis set seems to be better). The semiempirical SCC-DFTB-D method is able to capture the qualitative trends of geometrical and energetics properties and is in excellent agreement with RI-SCS-MP2, which is probably somewhat fortuitous. In this work, we took advantage of the low computational cost of SCC-DFTB-D and applied this method to size extrapolations using very large model systems, up to C864H72. We have demonstrated that acetone−graphene interaction energies are accurately fitted to an exponentially 9007

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C decaying function with respect to the increase in graphene flake size, allowing an extrapolation to an infinitely large graphene sheet. It should be noted that the calculations up to dicircumcoronene C96H24 complexes are sufficient to estimate the converged interaction energies. The precise interaction energy of acetone with a pristine single-layer graphite sheet has been determined via a two-dimensional basis set and graphene flake size extrapolation. At the convergence limit, we determine a clear acetone orientation dependence (PP > UP > DP). In particular, the RI-SCS-MP2 interaction energy of acetone on graphene sheet was determined to be 42.6 kJ/mol for PP, 30.1 kJ/mol for UP, and 20.3 kJ/mol for DP structures within chemical accuracy limits. The double-extrapolation scheme presented here demonstrates that SCC-DFTB-D can be used for size extrapolation in combination with high-level ab initio RI-SCS-MP2 calculations to accurately estimate binding energies of small molecules on nanosized systems.



(3) Kong, J.; Franklin, N. R.; Zhou, C. W.; Chapline, M. G.; Peng, S.; Cho, K.; Dai, H. Nanotube Molecular Wires as Chemical Sensors. Science 2000, 287, 622−625. (4) Snow, E. S.; Perkins, F. K.; Houser, E. J.; Badescu, S. C.; Reinecke, T. L. Chemical Detection with a Single-Walled Carbon Nanotube Capacitor. Science 2005, 307, 1942−1945. (5) Snow, E. S.; Perkins, F. K. Capacitance and Conductance of Single-Walled Carbon Nanotubes in the Presence of Chemical Vapors. Nano Lett. 2005, 5, 2414−2417. (6) Bianco, A.; Kostarelos, K.; Partidos, C. D.; Prato, M. Biomedical Applications of Functionalised Carbon Nanotubes. Chem. Commun. 2005, 571−577. (7) Kam, N. W. S.; O’Connell, M.; Wisdom, J. A.; Dai, H. Carbon Nanotubes as Multifunctional Biological Transporters and NearInfrared Agents for Selective Cancer Cell Destruction. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 11600−11605. (8) Bianco, A.; Kostarelos, K.; Prato, M. Applications of Carbon Nanotubes in Drug Delivery. Curr. Opin. Chem. Biol. 2005, 9, 674− 679. (9) Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (10) Č ížek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256−4266. (11) Watts, J. D.; Gauss, J.; Bartlett, R. J. Coupled-Cluster Methods with Noniterative Triple Excitations for Restricted Open-Shell Hartree-Fock and Other General Single Determinant Reference Functions. Energies and Analytical Gradients. J. Chem. Phys. 1993, 98, 8718−8733. (12) Sinnokrot, M. O.; Valeev, E. F.; Sherrill, C. D. Estimates of the Ab Initio Limit for π-π Interactions: The Benzene Dimer. J. Am. Chem. Soc. 2002, 124, 10887−10893. (13) Pitoňaḱ , M.; Neogrády, P.; Ř ezác,̌ J.; Jurečka, P.; Urban, M.; Hobza, P. Benzene Dimer: High-Level Wave Function and Density Functional Theory Calculations. J. Chem. Theory Comput. 2008, 4, 1829−1834. (14) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. An Assessment of Theoretical Methods for Nonbonded Interactions: Comparison to Complete Basis Set Limit Coupled-Cluster Potential Energy Curves for the Benzene Dimer, the Methane Dimer, Benzene-Methane, and Benzene-H2S. J. Phys. Chem. A 2009, 113, 10146−10159. (15) Katouda, M.; Nagase, S. Efficient Parallel Algorithm of SecondOrder Møller-Plesset Perturbation Theory with Resolution-of-Identidy Approximation (RI-MP2). Int. J. Quantum Chem. 2009, 109, 2121− 2130. (16) Janowski, T.; Pulay, P. A Benchmark Quantum Chemical Study of the Stacking Interaction Between Larger Polycondensed Aromatic Hydrocarbons. Theor. Chem. Acc. 2011, 130, 419−427. (17) Katouda, M.; Naruse, A.; Hirano, Y.; Nakajima, T. Massively Parallel Algorithm and Implementation of RI-MP2 Energy Calculation for Peta-Scale Many-Core Supercomputers. J. Comput. Chem. 2016, 37, 2623−2633. (18) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864−B871. (19) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (20) Kristyán, S.; Pulay, P. Can (semi)Local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175−180. (21) Pérez-Jordá, J. M.; Becke, A. D. A Density-Functional Study of van der Waals Forces: Rare Gas Diatomics. Chem. Phys. Lett. 1995, 233, 134−137. (22) Andersson, Y.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Interactions in Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 102−105.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b13002. Optimized geometrical parameters of acetone and graphene flakes, angle between acetone CCCO plane and graphene flake in PP structures, calculated interaction energy data for different optimized structures, and list of optimized parameter sets for graphene flake size extrapolation (PDF).



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Michio Katouda: 0000-0001-7980-5386 Stephan Irle: 0000-0003-4995-4991 Present Address

‡ Research Institute for Science and Engineering, Waseda University, Tokyo 169−8555, Japan.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Prof. Keiji Morokuma for helpful discussions. S.I. acknowledges support by the Program for Improvement of Research Environment for Young Researchers from the Ministry of Education, Culture, Sports, Science and Technology Agency/Core Research for Evolutional Science and Technology Grant in the area of “High Performance Computing for Multiscale and Multiphysics Phenomena”, and by a special grant from the Japan Society for the Promotion of Science (JSPS) in Priority area “Molecular Theory for Real Systems”.



REFERENCES

(1) Dillon, A. C.; Jones, K. M.; Bekkedahl, T. A.; Kiang, C. H.; Bethune, D. S.; Heben, M. J. Storage of Hydrogen in Single-Walled Carbon Nanotubes. Nature 1997, 386, 377−379. (2) Liu, C.; Fan, Y. Y.; Liu, M.; Cong, H. T.; Cheng, H. M.; Dresselhaus, M. S. Hydrogen Storage in Single-Walled Carbon Nanotubes at Room Temperature. Science 1999, 286, 1127−1129. 9008

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C (23) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (24) Sato, T.; Tsuneda, T.; Hirao, K. Van der Waals Interactions Studied by Density Functional Theory. Mol. Phys. 2005, 103, 1151− 1164. (25) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103, 063004. (26) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (27) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (28) 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. (29) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364−382. (30) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (31) Dewar, M. J. S.; Thiel, W. Ground States of Molecules. 38. The MNDO Method. Approximations and Parameters. J. Am. Chem. Soc. 1977, 99, 4899−4907. (32) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (33) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods I. Method. J. Comput. Chem. 1989, 10, 209−220. (34) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of Tight-Binding-Like Potentials on the Basis of Density-Functional Theory: Application to Carbon. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 12947−12957. (35) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260−7268. (36) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen Bonding and Stacking Interactions of Nucleic Acid Base Pairs: A Density-Functional-Theory Based Treatment. J. Chem. Phys. 2001, 114, 5149−5155. (37) Zhechkov, L.; Heine, T.; Patchkovskii, S.; Seifert, G.; Duarte, H. A. An Efficient a Posteriori Treatment for Dispersion Interaction in Density-Functional-Based Tight Binding. J. Chem. Theory Comput. 2005, 1, 841−847. (38) Xu, S.; Irle, S.; Musaev, D. G.; Lin, M. C. Water Clusters on Graphite: Methodology for Quantum Chemical A Priori Prediction of Reaction Rate Constants. J. Phys. Chem. A 2005, 109, 9563−9572. (39) Lin, C. S.; Zhang, R. Q.; Lee, S. T.; Elstner, M.; Frauenheim, T.; Wan, L. J. Simulation of Water Cluster Assembly on a Graphite Surface. J. Phys. Chem. B 2005, 109, 14183−14188. (40) Feng, X.; Irle, S.; Witek, H.; Morokuma, K.; Vidic, R.; Borguet, E. Sensitivity of Ammonia Interaction with Single-Walled Carbon Nanotube Bundles to the Presence of Defect Sites and Functionalities. J. Am. Chem. Soc. 2005, 127, 10533−10538. (41) Ganji, M. D.; Goodarzi, M.; Nashtahosseini, M.; Mommadinejad, A. Theoretical Studies on Interaction Between Methanol and Functionalized Single-Walled Carbon Nanotubes. Commun. Theor. Phys. 2011, 55, 365−370.

(42) Xu, S. C.; Irle, S.; Musaev, D. G.; Lin, M. C. Quantum Chemical Prediction of Reaction Pathways and Rate Constants for Dissociative Adsorption of COx and NOx on the Graphite (0001) Surface. J. Phys. Chem. B 2006, 110, 21135−21144. (43) Xu, S. C.; Irle, S.; Musaev, D. G.; Lin, M. C. Quantum Chemical Study of the Dissociative Adsorption of OH and H2O on Pristine and Defective Graphite (0001) Surfaces: Reaction Mechanisms and Kinetics. J. Phys. Chem. C 2007, 111, 1355−1365. (44) Kwon, S.; Russell, J.; Zhao, X.; Vidic, R. D.; Johnson, J. K.; Borguet, E. Combined Experimental and Theoretical Investigation of Polar Organic Adsorption/Desorption from Model Carbonaceous Surfaces: Acetone on Graphite. Langmuir 2002, 18, 2595−2600. (45) Elkington, P. A.; Curthoys, G. Heats of Adsorption on Carbon Black Surfaces. J. Phys. Chem. 1969, 73, 2321−2326. (46) Gales, L.; Mendes, A.; Costa, C. Hysteresis in the Cyclic Adsorption of Acetone, Ethanol and Ethyl Acetate on Activated Carbon. Carbon 2000, 38, 1083−1088. (47) Dinger, A.; Lutterloh, C.; Biener, J.; Küppers, J. Reactions of Hydrogen Atoms with Acetone Monolayers Adsorbed on Graphite Monolayer Covered Pt(111) Surfaces. Surf. Sci. 1999, 437, 116−124. (48) Shih, Y.-h.; Li, M.-s. Adsorption of Selected Volatile Organic Vapors on Multiwall Carbon Nanotubes. J. Hazard. Mater. 2008, 154, 21−28. (49) Chakrapani, N.; Zhang, Y. M.; Nayak, S. K.; Moore, J. A.; Carroll, D. L.; Choi, Y. Y.; Ajayan, P. M. Chemisorption of Acetone on Carbon Nanotubes. J. Phys. Chem. B 2003, 107, 9308−9311. (50) Guirado-López, R. A.; Sánchez, M.; Rincón, M. E. Interaction of Acetone Molecules with Carbon-Nanotube-Supported TiO2 Nanoparticles: Possible Applications as Room Temperature Molecular Sensitive Coatings. J. Phys. Chem. C 2007, 111, 57−65. (51) Kazachkin, D.; Nishimura, Y.; Irle, S.; Morokuma, K.; Vidic, R. D.; Borguet, E. Interaction of Acetone with Single Wall Carbon Nanotubes at Cryogenic Temperatures: A Combined Temperature Programmed Desorption and Theoretical Study. Langmuir 2008, 24, 7848−7856. (52) Kazachkin, D. V.; Nishimura, Y.; Irle, S.; Feng, X.; Vidic, R.; Borguet, E. Temperature and Pressure Dependence of Molecular Adsorption on Single Wall Carbon Nanotubes and the Existence of an “Adsorption/Desorption Pressure Gap”. Carbon 2010, 48, 1867− 1875. (53) Kazachkin, D. V.; Nishimura, Y.; Witek, H. A.; Irle, S.; Borguet, E. Dramatic Reduction of IR Vibrational Cross Sections of Molecules Encapsulated in Carbon Nanotubes. J. Am. Chem. Soc. 2011, 133, 8191−8198. (54) TURBOMOLE, V5.10 2008, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com. (55) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (56) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C., et al. Gaussian 03, Revision E.01; Gaussian, Inc.: Wallingford, CT, 2004. (57) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (58) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (59) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (60) Weigend, F.; Furche, F.; Ahlrichs, R. Gaussian Basis Sets of Quadruple Zeta Valence Quality for Atoms H-Kr. J. Chem. Phys. 2003, 119, 12753−12762. 9009

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010

Article

The Journal of Physical Chemistry C (61) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (62) Hättig, C. Optimization of Auxiliary Basis Sets for RI-MP2 and RI-CC2 Calculations: Core-Valence and Quintuple-ζ Basis Sets for H to Ar and QZVPP Basis Sets for Li to Kr. Phys. Chem. Chem. Phys. 2005, 7, 59−66. (63) Häser, M.; Ahlrichs, R. Improvements on the Direct SCF Method. J. Comput. Chem. 1989, 10, 104−111. (64) Weigend, F.; Häser, M. RI-MP2 First Derivatives and Global Consistency. Theor. Chem. Acc. 1997, 97, 331−340. (65) Grimme, S. Improved Second-Order Møller−Plesset Perturbation Theory by Separate Scaling of Parallel- and Antiparallel-Spin Pair Correlation Energies. J. Chem. Phys. 2003, 118, 9095−9102. (66) Kamiya, M.; Tsuneda, T.; Hirao, K. A Density Functional Study of van der Waals Interactions. J. Chem. Phys. 2002, 117, 6010−6015. (67) Sato, T.; Tsuneda, T.; Hirao, K. A Density-Functional Study on π-Aromatic Interaction: Benzene Dimer and Naphthalene Dimer. J. Chem. Phys. 2005, 123, 104307. (68) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A Long-Range Correction Scheme for Generalized-Gradient-Approximation Exchange Functionals. J. Chem. Phys. 2001, 115, 3540−3544. (69) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (70) Tsuneda, T.; Suzumura, T.; Hirao, K. A New One-Parameter Progressive Colle-Salvetti-Type Correlation Functional. J. Chem. Phys. 1999, 110, 10664−10678. (71) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-Range-Corrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425−8433. (72) Song, J.-W.; Hirosawa, T.; Tsuneda, T.; Hirao, K. Long-Range Corrected Density Functional Calculations of Chemical Reactions: Redetermination of Parameter. J. Chem. Phys. 2007, 126, 154105. (73) Sato, T.; Tsuneda, T.; Hirao, K. Long-Range Corrected Density Functional Study on Weakly Bound Systems: Balanced Descriptions of Various Types of Molecular Interactions. J. Chem. Phys. 2007, 126, 234114. (74) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (75) Feller, D. Application of Systematic Sequences of Wave Functions to the Water Dimer. J. Chem. Phys. 1992, 96, 6104−6114. (76) Feller, D. The Use of Systematic Sequences of Wave Functions for Estimating the Complete Basis Set, Full Configuration Interaction Limit in Water. J. Chem. Phys. 1993, 98, 7059−7071. (77) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639−9646. (78) Kahn, K.; Bruice, T. C. Systematic Convergence of Energies with Respect to Basis Set and Treatment of Electron Correlation: Focal-Point Conformational Analysis of Methanol. Theor. Chem. Acc. 2004, 111, 18−24. (79) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. Origin of Attraction and Directionality of the π/π Interaction: Model Chemistry Calculations of Benzene Dimer Interaction. J. Am. Chem. Soc. 2002, 124, 104−112. (80) Bondi, A. Van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (81) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. The Magnitude of the CH/π Interaction Between Benzene and Some Model Hydrocarbons. J. Am. Chem. Soc. 2000, 122, 3746−3753. (82) Tsuzuki, S. In Intermolecular Forces and Clusters; Wales, D. J., Ed.; Springer: Berlin, 2005; Vol. 115, pp 149−193. (83) Feller, D.; Jordan, K. D. Estimating the Strength of the Water/ Single-Layer Graphite Interaction. J. Phys. Chem. A 2000, 104, 9971− 9975.

9010

DOI: 10.1021/acs.jpcc.6b13002 J. Phys. Chem. C 2017, 121, 8999−9010