Accurate and Efficient Calculation of the Desorption Energy of Small

Article pubs.acs.org/JPCC

Accurate and Efficient Calculation of the Desorption Energy of Small Molecules from Graphene Simone Conti and Marco Cecchini* Laboratoire d’Ingénierie des Fonctions Moléculaires ISIS, UMR 7006 CNRS, Université de Strasbourg, F-67083 Strasbourg Cedex, France S Supporting Information *

ABSTRACT: Understanding the energetics of molecular adsorption at nanostructured carbon materials opens up to outstanding applications in the postgraphene era. The accurate determination of the adsorbate binding affinity remains however a challenge both experimentally and computationally. Here, we report on the determination of the desorption energy of 25 chemically diverse compounds from graphite using computational methods at different levels of theory: empirical force fields (FF), semiempirical quantum mechanics (SQM), and density functional theory (DFT). By comparing the computational predictions with literature temperature-programmed desorption (TPD) experiments we found that the dispersion-corrected semiempirical method PM6-DH+ yields desorption energies in quantitative agreement with the experiments with an average error of 1.25 kcal mol−1. The discovery of a fast and accurate approach to molecular adsorption at surfaces prompted us to search the chemical space of the adsorbate to optimize the binding affinity for graphene. We found that polarizable groups containing sulfur and halogen atoms significantly enhance the interaction strength with the substrate. In particular, we predict that per-chlorination of aliphatic hydrocarbons doubles the desorption energy from graphene in vacuum. The efficiency of the PM6-DH+ calculations, which allows for screening libraries of compounds, guided the design of potentially improved graphene surfactants, which are commercially available.

I. INTRODUCTION Providing a chemical understanding of the energetics underlying the interaction of small molecules with nanostructured carbon materials is of fundamental importance and finds technological applications in nanoelectronics,1 catalysis,2 gas sensing,3,4 energy storage5 and conversion.6 For instance, the extraordinary surface area of graphene has been recently exploited for the design of smart materials for hydrogen storage7 or the extraction of aromatic pollutants from chemical waste.8,9 In addition, because the strength of the adsorbate/ substrate interaction determines the probability of binding at equilibrium, the energy of adsorption microscopically governs supramolecular processes such as molecular self-assembly at surfaces10 and the liquid-phase exfoliation of multilayered materials.11 Over the past decade, graphene has emerged as an exciting material with potential of impact in many areas of science and technology.12 The rise of graphene as a disruptive technology, however, strongly relies on the development of inexpensive strategies to produce high-quality material.13 Among the available methods, the liquid-phase exfoliation of graphite stands out as a versatile and potentially up-scalable approach to obtain high-quality graphene at low cost.11 Recently, it has been demonstrated that the yield of exfoliation in organic media is significantly enhanced by the addition of dispersion-stabilizing agents such as linear alkanes14 and long-chain fatty acids,15 © 2014 American Chemical Society

which form tightly packed self-assembled monolayers on the dispersed graphene sheets. The design of improved surfactants, whose binding affinity for the basal plane of graphene is optimized in a given solvent, appears as a promising strategy to shorten and simplify the production of graphene for technological applications. Despite its unquestionable importance, the accurate determination of the interaction energy of an isolated molecule with graphene remains a challenging problem both experimentally and computationally. Experimentally, the desorption energy from graphite can be accessed by temperatureprogrammed desorption (TPD) under ultrahigh vacuum conditions.16 Assuming that desorption kinetics are of first order in the adsorbate coverage, accurate determination of the adsorbate/substrate interaction strength can be obtained from a series of spectra recorded at different heating velocities17 or initial coverages.18 Multiple measurements are required to uncouple the height of the activation barrier from the preexponential factor of the rate constant, whose value is used for an indirect determination of the desorption energy.18 Although quantitative, this approach is indirect, time-consuming, and commonly limited to small-molecule adsorbates. Moreover, Received: October 17, 2014 Revised: November 28, 2014 Published: December 12, 2014 1867

DOI: 10.1021/jp5104774 J. Phys. Chem. C 2015, 119, 1867−1879

Article

The Journal of Physical Chemistry C

(MMFF94),54−58 which is able to handle natively most organic structural types. Applications of these methods to adsorption problems are not so common in the literature, perhaps due to a widespread skepticism toward the level of accuracy that can be achieved by the use of highly parametrized representations. In this work, the performance of various modeling approaches to molecular adsorption on graphene are benchmarked by comparing the desorption energies for a wide set of small molecules predicted using different levels of theory with TPD literature data. We have found that dispersion-corrected semiempirical methods provide the best compromise between accuracy and efficiency of the calculation with an average error of 1.25 kcal mol−1 on the desorption energy. The model based on the PM6-DH+ method was used to obtain chemical insight on how to optimize the binding affinity for graphene and used it to design synthetically accessible compounds with enhanced interaction energies with the substrate. Implications of these results on the design of improved dispersion stabilizing agents for graphite exfoliation are discussed.

special care must be taken to account for the coverage dependence of the desorption energy, which arises from both the presence of nonequivalent adsorption sites on the substrate (i.e., defects) and adsorbate−adsorbate interactions in the monolayer.19 Most of the problems encountered experimentally could be solved, in principle, by modeling the desorption process and accessing the adsorbate/substrate binding energy via high-level of theory calculations. Although the computational approach would allow for a direct determination of the activation barrier the size of the substrate to produce correct geometries and obtain accurate interaction energies is prohibitive for most ab initio methods and only a few coupled-cluster (CC) analysis of molecular adsorption have been reported so far.20,21 Due to the superior computational efficiency, approaches based on density functional theory (DFT) appear to be better suited and a number of successful applications have been reported.22−25 Nonetheless, the computational results are strongly dependent on the selection of the exchange-correlation functional21 and the basis set in use, and must include corrections for dispersive interactions.26,27 Most importantly, as the computational cost associated with a single desorption-energy determination scales with the number of basis functions, DFT analysis are bound to a restricted set of molecular compounds. Semiempirical quantum mechanics (SQM) offers a valuable alternative to molecular adsorption as it scales more favorably with the size of the system than ab initio or DFT methods. These simplified schemes treat the valence electrons explicitly and rely on a parametrization of the core electrons. Provided that a proper parametrization is available, a minimum basis set can be used, which results in faster converging calculations with no compromise on the accuracy. A wide variety of SQM methods are presently available. They include the original implementations AM128 and PM3,29,30 which were reparameterized in RM131 and PM6,32 respectively. More recently, empirical corrections for both dispersion interactions and Hbonding were included, leading to the development of the AM1-D and PM3-D,33 the PM6-DH,34 AM1-FS1,35 PM6DH2,36 PM6-DH+,37 and finally the D3H4 correction.38 The last-released method is PM7,39 which was parametrized using a large set of data and including corrections for both dispersion and hydrogen bonding natively. Recently, these dispersioncorrected methods were shown to provide an accurate description of the adsorbate/substrate interaction for a series of gas molecules and a few aromatic compounds on graphene.40 A last class of methods that can be used to probe molecular adsorption at surfaces includes empirical or force-fields methods, in which electrons are not treated explicitly. A clear advantage of these methods is the efficiency of the calculations, which makes them suitable to explore systems involving several thousands of atoms. Of course, the use of highly parametrized energy models is expected to introduce inaccuracies in the description of interactions and strictly requires that parameters for all atom types are known. In recent years, a great effort has been directed to develop generalized force fields. In this endeavor, the OPLS force field, which was originally designed for protein simulations,41 was gradually extended to the OPLSAA,42−48 which includes more chemotypes. Similarly the general Amber force field (GAFF)49,50 and the CHARMM general force field (CGenFF)51 were developed as an extension of the AMBER52 and the CHARMM53 force fields, respectively. Also, there have been efforts to implement general force field from scratch, such as the Merck molecular force field

II. MATERIAL Two sets of experimental desorption energies from highly oriented pyrolitic graphite (HOPG) were used in this work. These values were obtained from TPD measurements under ultra high-vacuum conditions and used to benchmark the performance of several computational approaches to molecular adsorption. The data include desorption energies for seven linear alkanes from C1 (methane) to C10 (decane), a few aromatic compounds and a series of inert gases and solvents; see Table 1. The first subset, which contains all linear alkanes, was obtained from Tait et al.19 In this work, desorption energies were accessed by multiple TPD measurements collected at different initial coverages using the inversion-optimization method.18 To account for substrate defects, the values reported in Table 1 were measured at 0.5 monolayer coverages. These data show that the desorption energy of n-alkanes scales linearly with the chain length (N) and produce a physically meaningful intercept of 1.7 kcal mol−1 when extrapolated to zero length (N = 0). Because conformational entropy is expected to play no role on the desorption energy of short-length alkanes,17 this data set represents a stringent test for modeling. In order to increase the chemical diversity of the benchmark, the TPD measurements of Ulbricht et al.59 were also considered. This second subset provides desorption energies from graphite for 23 additional compounds including polyaromatic hydrocarbons, inert gases, solvents, and hydrogen bonding molecules. Among them, the data for methane and fullerene were not considered. Coronene and ovalene were also excluded from the data set, because of the apparent disagreement between independent determinations of their desorption energies in the literature.59−61 Finally, to avoid complications with the treatment of radicals, nitrogen dioxide was also discarded. This provides 18 additional desorption energies, whose values are also given in Table 1. Unlike the data of Tait et al.,19 the desorption energies of Ulbricht et al.59 were determined using a Redhead maximum-peak analysis16 and a pre-exponential factor obtained from equilibrium vaporpressure constants and are expected to be of somewhat lower quality. 1868


Article


hydrocarbon mimicking a graphite substrate (Figure 1) was determined by subtracting the energies of the adsorbate and the graphite model in isolation from that of the complex. Because molecular adsorption is modeled here as a barrier-less exothermic reaction, the desorption energy is obtained from the calculated binding energy with opposite sign as

Table 1. Experimental Desorption Energies for a Series of Small Organic Compounds from Graphitea molecule

Edes [ kcal mol−1 ]

Linear Alkanes methane 3.37 ethane 5.88 propane 7.67 butane 9.75 hexane 15.05 octane 17.34 decane 21.83 Aromatic Compounds benzene 11.52 ± 1.92 naphtalene 18.48 ± 2.16 toluene 16.32 ± 1.68 ethylbenzene 18.96 ± 2.4 o-dichlorobenzene 16.56 ± 1.44 Gas and Solvents methanol 11.52 ± 0.72 1,1-dichloroethane 12.24 ± 0.72 N,N-dimethylformamide (DMF) 12.72 ± 0.96 ethanol 12.00 ± 0.72 water 11.04 ± 0.72 ammonia (NH3) 6.00 ± 0.48 trichloromethane (CHCl3) 12.96 ± 0.72 nitrogen (N2) 3.12 ± 0.24 carbon dioxide (CO2) 5.76 ± 0.48 carbon monoxide (CO) 3.12 ± 0.24 oxygen (O2) 2.88 ± 0.24 sulfur hexafluoride (SF6) 7.44 ± 0.48 xenon (Xe) 5.76 ± 0.48

ref Tait Tait Tait Tait Tait Tait Tait

et et et et et et et

al.19 al.19 al.19 al.19 al.19 al.19 al.19

Edes = −[Ea / s − (Ea + Es)]

Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht

et et et et et

al.59 al.59 al.59 al.59 al.59

Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht Ulbricht

et et et et et et et et et et et et et

al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59 al.59

(1)

with Ea/s the potential energy of the complex in vacuum at the minimum, and Ea and Es the potential energies of the adsorbate and the substrate separated at an infinite distance. In practice, the desorption energy of a test molecule on graphite was calculated as follows. First, the graphite model was energy minimized in vacuum fixing the z coordinate of its atoms to zero. Then, the adsorbate was placed at the center of the substrate and the geometry of the complex optimized with the position of the graphite atoms fixed in space to obtain Ea/s. The adsorbate and the substrate were then separated at an infinite distance and their single-point energies determined to obtained Ea and Es. Introducing these three energy values in eq 1 provides an estimate of the desorption energy of the test compound. To keep the procedure as general as possible, no particular attention was paid to the initial position of the adsorbate on the substrate. The error introduced by a suboptimal placement of the adsorbate was found to be as large as 0.38 kcal mol−1 for benzene, which corresponds to the difference between the AB staking and the sandwich configurations of the complex. Such an error is acceptable particularly when compared to the uncertainty on the experimental determinations of the desorption energy. Special care was taken only to verify that all adsorbates were lying flat on the surface, and for all linear alkanes that a zigzag configuration as the one depicted in Figure 1 was adopted. B. Computational Methods. 1. Density Functional Theory (DFT). In this work, the ωB97X-D dispersion corrected DFT functional62 in combination with the 6-31G(d) basis set was used. This functional was recently shown to be the best to model dispersion-dominated problems.63 The large computational effort involved in these calculations did not allow us to evaluate the performances of other functionals. Desorption energies were calculated using a small hexagonal graphene model including 216 carbon atoms and dangling hydrogens to

a

All values were obtained from temperature programmed desorption (TPD) experiments in ultra-high vacuum. For each determination the source is given.

III. THEORETICAL METHODS A. Calculation of the Desorption Energy. Desorption energies from graphene were computed for all compounds in Table 1 using three levels of theory: density functional theory (DFT), semiempirical quantum mechanics (SQM) and empirical force fields (FF). For this purpose, the binding energy of a small molecule adsorbed to a large poly aromatic

Figure 1. Graphite models for the calculation of the desorption energy. (a) Model used for both DFT and SQM (216 carbon atoms). (b) Model used with the CGenFF force field (610 carbon atoms). (c) Model used with the GAFF and OPLS-AA force fields (layer dimension: 103 × 110 Å) which include periodic boundary conditions. In MMFF94 calculations, a graphite model similar to that represented in part c, but with double size (layer dimension: 206 × 216 Å), was used. In these rendering the zigzag configuration of the linear alkanes adsorbed on the substrate is shown for decane. 1869


Article


Fletcher−Goldfarb−Shanno approach)76,77 until convergence to an energy gradient of 10−8 kcal mol−1 Å−1. OPLS-AA. Atom typing was performed via the MKTOP software78 that automatically produces topologies suitable for GROMACS. The determination of the desorption energy was done as described for the GAFF force field. A graphene layer of 103 × 110 Å with periodic boundaries was used with carbon atoms parametrized using the type opls_147 (naphtalene fusion carbon); see Figure 1c. C. Error Analysis. The performances of the various methods were evaluated from the correlation between the experimental desorption energies and the computational predictions. Specifically, the values of the slope, the correlation coefficient R2, and the intercept provide measures of the accuracy, the precision, and the systematic error of the calculation, respectively. The percent deviation of the slope from the idealized-correlation behavior, i.e., slope equals to unity, was used to quantify the accuracy of the methods. A quantitative measure of the predictivity, which is useful for ranking, was obtained from the root mean square errors (RMSE) on the predicted desorption energy

saturate the valence of the carbons on the edge; see Figure 1a. The energy of interaction between the adsorbate and the substrate was calculated using the counterpoise correction to reduce the basis set superimposition error (BSSE).64 All calculations were performed with Gaussian 09.65 Since the chosen basis set is defined from H to Zn, xenon was not considered. 2. Semiempirical Quantum Mechanics (SQM). The following SQM methods were tested: AM1, AM1-D3H4/ AM1, RM1, RM1-D3H4/RM1, PM6, PM6-DH2/PM6, PM6DH+, PM6-D3H4/PM6, PM6-D3H4/PM7, and PM7, which are all available in MOPAC201266 except for the D3H4 correction. Here, the notation methodA/methodB indicates that the desorption energy was evaluated using methodA on the geometry optimized with methodB. This is useful for instance with the PM6-DH2 method, which presents discontinuities in the energy gradient,37 and the D3H4 correction where geometry optimization was performed with AM1, RM1, PM6, or PM7 and the Rezac H438,67 and the Grimme D327,68 corrections were applied a posteriori. Similar to DFT, desorption energies were calculated using a hexagonal graphene model with 216 carbon atoms and dangling hydrogens; see Figure 1a. Analysis of the substrate size dependence of the desorption energy of benzene (Figure S1, Supporting Information) fully justified the use of this molecular model for the substrate. 3. Empirical Force Fields (FF). The performances of four classical force fields were evaluated using general procedures based on automatic atom typing. A short description of these calculations is given below. MMFF94. Atom typing was done by processing a sybyl mol2 file in CHARMM (version 36b2).69 Desorption energies were calculated using a graphene model of 206 × 216 Å in size with periodic boundaries; see Figure 1c. Carbon atoms were parametrized using the atom type 2 (generic sp2 carbon) with zero charge. For the calculation of the desorption energy, the adsorbate in complex with the substrate was energy minimized by 1000 steps of the steepest descent (SD) followed by an adopted basis Newton−Raphson (ABNR) optimization until convergence to an energy gradient of 10−8 kcal mol−1 Å−1. CGenFF. Atom typing was performed via the ParamChem Web server (version 0.9.7),70,71 which generates parameters for CGenFF version 2b8. The graphene substrate was modeled using a squared 610 carbons polyaromatic molecule with dangling hydrogens; see Figure 1b. The carbon atoms were parametrized using the atom type CG2R61 (aromatic carbon) with zero charge if internal and −0.115 e.u. on the edge. Hydrogens were parametrized with the atom type HGR61 (aromatic hydrogen) with a partial charge of +0.115 e.u.. The desorption energy was calculated with CHARMM as described for MMFF94. GAFF. Atom typing for GAFF is natively implemented in the ANTECHAMBER50 software, which is distributed with AmberTools. The calculation of the desorption energy was performed with GROMACS 4.5.572−74 and suitable topologies were created using the ACPYPE75 script. The graphene substrate was modeled as a bidimensional 103 × 110 Å periodic layer, see Figure 1c, with carbon atoms parametrized with the atom type CA (sp2 aromatic carbon). For the desorption energy calculations, the adsorbate on the substrate was energy minimized using the l-bfgs method (a quasiNewtonian algorithm based on the low-memory Broyden−

RMSE =

1 N

N

∑ (Eicalc − Eiexp)2 (2)

i

Eexp i

Ecalc i

where N is the number of data points, and are the experimental desorption energy and the computational prediction for the ith compound, respectively. Similarly, eq 2 was used to evaluate the predictivity of the models obtained by best-fitting the computational predictions to the set of experimental data. In this case, the predicted desorption energies Ecalc i , are corrected by the slope m and the intercept q of the linear fit, as Emodel = mEcalc + q. i i D. Polarizability. The isotropic polarizability of all molecules in Table 1 was evaluated and correlated with the calculated desorption energy. The calculation of the polarizability was performed with the PM6-DH+ semiempirical method on the conformation of the adsorbate optimized in the adsorption complex but in the absence of substrate using the finite-field method implemented in the MOPAC2012 software (POLAR module).79 The value of the isotropic polarizability was then obtained as an average of the diagonal elements of the diagonalized polarizability tensor.

IV. RESULTS AND DISCUSSION A. Modeling Adsorption on Graphene. Desorption energies calculated for all compounds in Table 1 with different levels of theory were compared with the experimental values from TPD. The semiempirical (SQM) results are presented first, the force field (FF) and DFT calculations follow. At the semiempirical level of theory, the interaction energies with the substrate were evaluated using the original parametrizations AM1 and RM1, the more recent PM6 and the latest-generation PM7 along with a series of a posteriori corrections to account for dispersion interactions and hydrogen bonding; see Methods. The results indicate that using the AM1 and RM1 parametrizations, the calculated desorption energies are nearly zero and the equilibrium distance from graphene is between 4.5 and 5.0 Å for all compounds; see Supporting Information. Clearly, these methods are unable to describe the interaction with graphene and were not considered further. In sharp contrast, the PM6 model was able to capture the adsorbate/substrate interaction and predicted structures for the 1870


Article

The Journal of Physical Chemistry C Table 2. Calculated Desorption Energies Using Various Semiempirical Methodsa PM6-DH2/PM6

a

PM6-DH+

methane ethane propane butane hexane octane decane

3.12 4.44 6.20 7.97 11.58 15.37 18.86

3.18 4.73 6.58 8.44 12.55 16.02 19.72

benzene naphtalene toluene ethylbenzene o-dichlorobenzene

10.99 17.04 12.82 14.55 15.53

11.88 18.70 13.86 15.53 16.74

methanol 1,1-dichloroethane DMF ethanol water NH3 CHCl3 N2 CO2 CO O2 SF6 Xe

4.91 8.97 9.55 6.47 3.63 4.21 6.07 2.83 3.61 2.56 1.76 0.28 5.08

4.74 9.89 10.29 6.54 3.64 4.23 9.07 2.97 3.83 2.92 2.26 3.75 5.20

PM6-D3H4/PM6 Alkanes 3.23 4.54 6.40 8.08 11.76 15.77 19.31 Aromatic Compounds 9.38 14.23 7.84 12.89 13.03 Gas and Solvent Molecules 5.08 8.02 9.13 6.57 3.89 4.30 5.37 2.71 3.30 2.48 1.91 0.14 5.16

PM7

PM6-D3H4/PM7

experimental

4.58 7.54 10.15 13.20 18.99 24.97 31.38

3.31 4.71 6.57 8.40 12.16 15.89 19.34

3.37 5.88 7.67 9.75 15.05 17.34 21.83

17.28 27.41 20.37 23.10 33.67

9.85 15.24 11.89 13.53 11.83

11.52 18.48 16.32 18.96 16.56

± ± ± ± ±

1.92 2.16 1.68 2.4 1.44

6.33 24.52 14.42 9.30 4.19 5.11 23.47 3.80 5.11 3.79 3.28 8.41 23.23

4.29 7.59 9.31 6.14 3.84 4.25 6.31 2.70 3.49 2.52 2.19 2.86 5.29

11.52 12.24 12.72 12.00 11.04 6.00 12.96 3.12 5.76 3.12 2.88 7.44 5.76

± ± ± ± ± ± ± ± ± ± ± ± ±

0.72 0.72 0.96 0.72 0.72 0.48 0.72 0.24 0.48 0.24 0.24 0.48 0.48

All values are in kcal mol−1. Reference values from TPD experiments are given in the last column.

data show that empirical methods yield desorption energies in remarkable agreement with the TPD data, although they tend to underestimate the strength of interaction with the substrate. Consistent with the SQM results, the calculated desorption energies for the hydrogen-bonding compounds, i.e. methanol, ethanol, and water, are significantly lower than the experimental ones. Computational results are missing for CHCl3, H2, CO2, CO, O2, H2O, SF6, and Xe. This is due to failures in the automatic atom typing, which was performed with the aid of preprocessing software (see Methods) or lack of support for a particular chemotype in the force field, which both undermine the robustness of the approach on screening of large libraries of compounds. In this respect, GAFF appears to be the best supported and the most general force field. Desorption energies obtained by DFT with the ωB97X-D functional and the 6-31G(d) basis set are also reported in Table 3. The results are puzzling. In fact, although very accurate desorption energies were obtained for the short-length linear alkanes as well as the entire set of aromatic compounds, the binding affinity for the two longest alkanes were overestimated and those for the small gases and solvents significantly underestimated. These results evidence a clear instability of the DFT predictions, perhaps arising from the use of a too small basis set, which makes the predictions unreliable. Alternatively, the ωB97X-D is biased toward a specific class of molecules. Based on the limited set of data, we cannot conclude on any of these hypotheses. B. Benchmarking the Computational Methods. The performance of the various methods were compared by analyzing the correlation of the computed predictions with

complex with an average distance of 3.4 Å from the graphene surface. However, the strength of the interaction was considerably underestimated and the predicted desorption energies were off by 1 order of magnitude relative to the TPD results; see Table S1. Remarkably, evaluating the energy of binding using the single-point DH2 or the D3H4 corrections on PM6 optimized geometries provided desorption energies in quantitative agreement with the experimental results; see Table 2. The same is true for geometries optimized with the PM6-DH + method; see Supporting Information. In sharp contrast, the most recent parametrization PM7 yielded reasonable geometries but significantly overestimated the interaction strength, as also reported elsewhere.80 Consistently, this problem disappears for desorption energies predicted by the PM6D3H4 method on PM7-optimized geometries. In all cases, methanol, ethanol and water were predicted to have desorption energies in striking disagreement with the TPD results; see Table 2 and Figure S4. In particular, the calculated desorption energies are systematically lower than the experimental ones. Considering the chemical nature of these compounds, which expose H-bond donor and acceptor groups, the apparent inconsistency could arise from the formation of 2D supramolecular structures connected via H-bonding interactions. If so, the experimental desorption energies could actually correspond to the desorption of dimers or trimers, or to the energy required to destroy the supramolecular structure before desorption. Given the uncertainty in the desorption mechanism, these molecules were not further considered. The results obtained with the force field methods (MMFF94, CGenFF, GAFF, and OPLS-AA) are reported in Table 3. The 1871


Article

The Journal of Physical Chemistry C Table 3. Calculated Desorption Energies Using Force Field and DFT Methodsa MMFF

CGenFF

methane ethane propane butane hexane octane decane

2.32 3.43 4.76 6.15 8.96 11.76 14.56

3.05 4.88 6.81 8.78 12.72 16.14 20.56

benzene naphtalene toluene ethylbenzene o-dichlorobenzene

9.15 14.60 10.75 11.12 13.28

10.40 16.43 12.54 13.27 14.67

methanol 1,1-dichloroethane DMF ethanol water NH3 CHCl3 N2 CO2 CO O2 SF6 Xe

3.15 7.06 7.33 4.18 1.85 2.25 8.10 2.28

4.51 8.51 10.25 6.23

GAFF LINEAR ALKANES 3.30 5.10 7.19 9.32 13.61 17.85 22.10 Aromatic Compounds 11.26 17.92 13.64 13.51 15.67 Gas and Solvent Molecules 4.85 8.71 11.19 6.72

1.91

4.37

2.45 8.48 3.94 6.19 4.31 3.05

OPLS-AA

ωB97X-D

3.15 4.82 6.73 8.67 12.60 16.51 20.41

3.63 5.68 8.24 10.77 15.79 20.73 25.64

3.37 5.88 7.67 9.75 15.05 17.34 21.83

10.15 16.07 12.34 12.49 14.41

11.35 18.47 14.16 16.46 14.82

11.52 18.48 16.32 18.96 16.56

± ± ± ± ±

1.92 2.16 1.68 2.4 1.44

4.40 7.63 10.70 6.05

4.53 6.60 10.67 6.32 2.31 3.10 7.75 2.38 2.82 2.15 1.82 4.20

11.52 12.24 12.72 12.00 11.04 6.00 12.96 3.12 5.76 3.12 2.88 7.44 5.76

± ± ± ± ± ± ± ± ± ± ± ± ±

0.72 0.72 0.96 0.72 0.72 0.48 0.72 0.24 0.48 0.24 0.24 0.48 0.48

3.64 8.82 3.73

experimental

All values are in kcal mol−1. Empty values in MMFF, CGenFF, GAFF and OPLS-AA correspond to molecules non supported by the force fields (see text). Reference values from TPD experiments are given in the last column. a

Figure 2. Correlation plots of experimental versus calculated desorption energies for five semiempirical methods (PM6-DH2/PM6, PM6-DH+, PM6-D3H4/PM6, PM6-D3H4/PM7, and PM7), four force fields (MMFF94, CGenFF, GAFF, and OPLS-AA), and one DFT functional (ωB97XD/ 6-31G(d)). Legend: alkanes in blue, aromatic compounds in red, and small molecules in orange. Error bars correspond to the uncertainty in the experimental determination of the desorption energy.

The data show that the PM6-DH2/PM6, PM6-DH+, PM6D3H4/PM6, and PM6-D3H4/PM7 are in good agreement with the experimental data. Among them, the most accurate is PM6-DH+ with an error on the slope of only 2.10% and on R2

the experimental desorption energies. The linear correlations plots are shown in Figure 2. The values of the slope, the intercept and the correlation coefficient R2 are given in Table 4. 1872


Article

The Journal of Physical Chemistry C Table 4. Linear Regressions of Computed Versus Experimental Desorption Energiesa % error

RMSE

method

level of theory

slope

intercept

R2

slope

R2

method

model

PM6-DH+ PM6-DH2/PM6 PM6-D3H4/PM7 PM6-D3H4/PM6 PM7 GAFF CGenFF OPLS-AA MMFF94 ωB97XD/6-31G(d)

SQM SQM SQM SQM SQM FF FF FF FF DFT

1.02 1.02 1.13 1.07 0.53 1.03 1.02 1.08 1.30 0.80

1.28 2.19 1.42 2.45 2.29 0.83 1.90 1.33 1.80 3.01

0.95 0.90 0.91 0.83 0.77 0.90 0.93 0.91 0.93 0.88

2.10 1.70 13.48 7.08 47.01 2.55 2.28 8.30 30.01 19.77

4.68 10.40 8.69 16.64 23.20 9.71 6.88 9.23 7.43 12.01

1.94 2.99 3.12 3.83 7.38 2.14 2.57 2.77 4.65 2.62

1.25 1.87 1.71 2.36 2.79 1.84 1.43 1.67 1.50 2.02

a

The values reported here correspond to the slope, the intercept and the determination coefficient (R2) of the linear regression, the percentage deviation of the slope and the determination coefficient from the ideal behavior, and the root mean square error (RMSE) calculated before (method) and after (model) best-fitting to the experimental results.

Figure 3. Root mean square error (RMSE) of the various computational approaches relative to the TPD experiments. The results for the semiempirical, empirical and DFT methods in kcal mol−1 are shown by red, blue, and green bars, respectively. RMSE values are given before (a) and after (b) best-fitting to the TPD results. In both cases, the methods are ranked by RMSE in ascending order.

of 4.68%. Very accurate results were also obtained with PM6DH2, but the spread was twice as large as evidenced by the error on R2 of 10.40%. The two methods based on PM6-D3H4 provided reasonable results. PM6-optimized complexes yielded a better slope but a poorer spread, whereas the PM7-optimized geometry increased the correlation at the expense of the accuracy. In this respect, it would be interesting to analyze the performance of the PM6-D3H4 method based on native geometry optimizations. Finally, the PM7 method yielded the worst results with a slope of 0.53 (47% error) and the largest spread (23%). All empirical methods showed surprisingly good correlation with the TPD data, with an error on the R2 between 6% and 10%. In particular, CGenFF and GAFF provided accurate results with an error on the slope of 2.28% and 2.55%, respectively. Remarkably, the quality of these results is comparable with those predicted by the best semiempirical method, while being computationally more efficient (see below). The other two force fields, MMFF94 and OPLS-AA, produced slopes much larger than one, which indicates a systematic underestimation of the desorption energy. Finally, the DFT results show that the ωB97X-D functional in combination with the 6-31G(d) basis set, poorly captures the adsorbate/substrate interactions for most molecules in the data set and show a linear correlation with an error in slope as large as 20%.

To quantify the predictability of the various methods and produce a ranking, the root mean square error (RMSE) on the computational predictions were calculated and compared; see Methods. RMSE values are given in Table 4 and shown in Figure 3a. The results show that PM6-DH+ is the most accurate semiempirical method with an average error of 1.9 kcal mol−1. Also the empirical GAFF and CGenFF present errors on the order of 2 kcal mol−1, with performances comparable to the best SQM method. Worse results are obtained for the other methods. Best-fitting the computational predictions to the experimental results is expected to reduce the systematic error on the calculated desorption energies and produce more accurate models. The RMSE of the best-fit models calculated as described in the Methods are given in Table 4 and compared pictorially in Figure 3b. Now that the computational predictions are corrected for both the slope and the intercept, most models show an average error smaller than 2 kcal mol−1. Strikingly, the model based on the PM6-DH+ method shows an average error of 1.25 kcal mol−1. Also, models based on GAFF and CGenFF force fields show an average error of about 1.5 kcal mol−1. C. Efficiency and CPU Time. Beside the accuracy of the methods, it is important to evaluate the efficiency of the calculations, which confers, or not, the ability to explore the chemical space by the computational approach. To this aim, the 1873


Article


Although the comparison in Figure 4 refers to the total CPU time, which for DFT can be reduced to days of calculations if properly parallelized, these results indicate that high-level of theory calculations are inherently expensive and can be used to study only a small set of molecules. By contrast both empirical and semiempirical methods provide results at a moderate computational effort, which grants an efficient translation of the computational predictions into chemical information. D. Exploring the Chemical Space. The model based on the PM6-DH+ semiempirical method, which was shown to yield accurate desorption energies, was used to search the chemical space to maximize the binding affinity for graphene. For this purpose, a library of chemically related compounds was designed by introducing heteroatoms as nitrogen, oxygen, and sulfur on a benzene scaffold or by 1,3,5 substitution with halogen, hydroxyl, amino, nitro, methyl, and carboxyl groups; see Figure 5a. The aim of this analysis is to explore the effect of both heteroatomic substitution inside the aromatic ring and the exposure of polarizable, electron-donating, and electron-withdrawing groups on the binding affinity for graphene. As TPD measurements were not available for these compounds, the comparison of their desorption energies was done in predictive mode. We stress that extensive analyses like the one presented here are possible only through the use of a fast and accurate computational strategy. Desorption energies predicted by the PM6-DH+ model are given in Table 5. By comparing the desorption energies predicted for benzene (C6H6) and hexane (C6H12), the calculations indicate that the latter binds stronger to graphene by 0.7 kcal mol−1. This

average computational time to process one molecule were measured and compared; see Figure 4. The values reported

Figure 4. Average process time (CPU time) required to compute the desorption energy of one small-molecule adsorbate using different levels of theory. In DFT calculations which are run in parallel, the total CPU time corresponds to the sum of the computing time spent by each processors.

here refer to the total processing or CPU time summed over all processes, or threads, generated by a single run. The results show that force-field methods accomplish the task in minutes, whereas SQM methods do so in less than 1 h. DFT calculations are 3−4 orders of magnitude more expensive and need on average four months of CPU time per simulation.

Figure 5. Molecular structures of the library of compounds used to explore the chemistry of adsorption on graphene. (a) The first data set includes the molecules used to study the effect of the substitution of the benzene aromatic ring with heteroatoms and to compare the difference between aliphatic and aromatic compounds. The second data set contains 1,3,5-trisubstituted benzene species, used to evaluate the effect of electronwithdrawing, electron donor and polarizable groups on the binding affinity for graphene. (b) The structures of commercially available partially chlorinated undecane compounds are shown and compared with the fully hydrogenated and the theoretically fully chlorinated undecanes. 1874


Article


symmetric nature of these compounds. In all cases, the calculations predict that halogen substitution increases the affinity for graphene. In particular, the introduction of three strongly polarizable groups at the periphery of the aromatic ring results in desorption energies that are enhanced by 40% with chlorine, 80% with bromine and are almost doubled with iodine substituents. As given in Table 5, substitution with fluorine results in a 0.40 kcal mol−1 increase per group, whereas substitution by chlorine, bromine and iodine results in affinity enhancement of 2.27, 4.04, and 5.00 kcal mol−1 per group, respectively. Consistent with previous observations, the increase in polarizability by moving down the group significantly enhances the affinity for graphene. Finally, the desorption energy predicted for five benzene derivatives 1,3,5-substituted by electron-donating and electronwithdrawing groups were compared to explore inductive effects on the binding affinity for graphene. The results in Table 5 suggest that substitution with electron-withdrawing groups (carboxyl, nitro) enhances the strength of adsorption more than the substitution with electron-donating groups (methyl, hydroxyl and amino); there is an affinity gain >4.0 kcal mol−1 per group with carboxyl and nitro substituents, whereas it is

Accurate and Efficient Calculation of the Desorption Energy of Small

Recommend Documents