Theoretical study of the adsorption of aromatic units on carbon allotropes including explicit (empirical) DFT dispersion corrections and implicitly dispersion-corrected functionals: the pyridine case

Abstract

The suitability of implicitly dispersion-corrected functionals, namely the M06-2X, for the determination of interaction energies and electron polarization densities in adsorption studies of aromatic molecules on carbon allotropes surfaces is analysed by comparing the results with those obtained using explicit dispersion through Grimme’s empirical corrections. Several models of increasing size for the graphene sheet together with one-dimensional curved carbon structures, (5,5), (6,6) and (7,7) armchair single-walled nanotubes, and two-dimensional curved carbon structures, C60 fullerene, have been considered as substrates in this work, whereas pyridine has been chosen as an example for the adsorbed aromatic molecule. Comparison with recent experimental estimations of the adsorption energy and calculations using periodic boundary conditions on a supercell of 72 carbon atoms indicates that a finite model containing ninety six carbon atoms (C96) approaches quite well the adsorption on a graphene sheet. Analysis of the interaction energy components reveals that the M06-2X functional accounts for most of the dispersion energy implicitly, followed far by wB97X and B3LYP, whereas B97 and BLYP do not differ too much from HF. It has been found that M06-2X corrects only the energy component associated to dispersion and leaves the rest, electrostatic, Pauli and induction “unaltered” with respect to the other DFT functionals investigated. Moreover, only the M06-2X functional reflects the effect of dispersion on the electron polarization density, whereas for the remaining functionals the polarization density does not differ too much from the HF density. This makes the former functional more suitable a priori for the calculation of electron density related properties in these adsorption complexes.

---

Introduction

The potential application of graphene materials^1,2 has motivated in the last few years a huge amount of experimental and theoretical studies.^3–5 This was favoured by the extraordinary features of this material consisting of a monolayer of carbon atoms packed into a two-dimensional honeycomb lattice. Its mechanical strength and optical and electronic properties have turned graphene into a key material for the development of multiple scientific fields. In particular, its unique electronic structure confers high electron mobility and an inherent ability to interact with other molecules, allowing it to be employed as a transistor or a chemical sensor.^6,7 Since graphene and other carbon allotropes can serve as substrates for different analytical techniques, the study of molecular adsorption on graphene has gained significant attention lately.⁸ For instance, recent studies have demonstrated the potentialities of graphene as a substrate for surface-enhanced Raman scattering (SERS),^9,10 a very sensitive and selective technique whose valuable applications to molecular detection and structural analysis are unquestionable nowadays.^11–15 The same studies have pointed out the lack of information about the enhancement mechanism involved in the process. Therefore, understanding the electronic factors responsible for the Raman enhancement requires a previous exhaustive study of the interaction of the adsorbed molecule with the carbon substrate.

In accordance with that mentioned above, several theoretical studies can be found in the literature tackling the adsorption of molecules on these surfaces, either pristine or modified.^16–23 In general, these studies conclude that for a precise quantum mechanical description of the interaction of molecules with graphene surfaces, nanotubes or fullerenes is necessary to account properly for the dispersion forces. Therefore, the choice of a suitable computational method should consider the correct description of long-range electron correlation. Several examples can be found in the literature. Thus, Voloshina et al.¹⁹ employed the PBE functional with the empirical dispersion correction proposed by Grimme²⁴ in the study of the adsorption of pyridine on graphene. Lazar et al.²³ investigated the interaction energies of different molecules with graphene using ab initio molecular dynamics (AIMD) with the optB88-vdW DFT functional and force field (FF) simulations, which include a contribution for non-local correlations. In the same work, they also carried out a comparative study of different methods and DFT functionals (MP2, CCSD(T), optB88-vdW and M06-2X) and an energy decomposition analysis using DFT-SAPT²⁵ for a single sheet of coronene as a model of the graphene surface. In spite of using ab initio methods including dispersion, the reduced size of a single sheet of coronene may not incorporate well all the contributions to long-range attractive forces in graphene. Therefore, it seems that an exhaustive and systematic search for an adequate molecular model of the graphene surface to be used in adsorption studies and subsequent spectroscopic analysis at reasonable computational cost is still to be done.

Even more important is the fact that dispersion corrections in DFT via explicit dispersion energy terms can be useless when the interest is focused on the electronic structure or electron density distribution instead of the adsorption energy. For instance, the chemical enhancement in SERS depends mainly on the electron density deformation and the resonance enhancement can be related to electron transitions within the surface (electromagnetic) or between surface and molecule (charge transfer). These electronic effects, however, cannot be well described using explicit dispersion corrections and a comparison of the effect of implicit and explicit dispersion corrections on the electronic structure of adsorption complexes must be carried out in these cases.

In the present work the adsorption of pyridine on graphene was investigated systemically using simplified carbon 2D-structures of increasing size. We have used pyridine as a model for aromatic units. The adsorption energy was obtained using different DFT functionals including dispersion corrections, either implicitly or explicitly. The nature of the interaction was examined by analysis of the interaction energy components obtained with a partitioning scheme based on the electron density.²⁶ The effect of the surface curvature on the binding strength and interaction energy components was also investigated by including in our study carbon allotropes of curved morphologies, namely, armchair single-walled carbon nanotubes and C₆₀ fullerene.

Methodology and computational details

The decomposition of the interaction energy was performed using a methodology based on the partition of the complex electron density into monomer unperturbed densities and deformation terms.²⁶ This method is intrinsically related to the perturbative treatment of intermolecular interactions, although the input electron density matrices can be calculated at any computational level. The interaction energy is then decomposed into the following components:

E_int = E_elec + E_rep + E_x + E_def + E_xc-def = E_elec + E_Pauli + E_pol (1)

where the first three terms are first order terms and correspond to the electrostatic energy, repulsion and intermolecular exchange. The sum of the repulsion, E_rep, and the first order intermolecular exchange, E_x, corresponds to the Pauli repulsion. These terms have exactly the same physical meaning as the first order terms obtained with symmetry adapted perturbation theory, SAPT, and their values coincide quantitatively.²⁶ The polarization energy, which corresponds to the last two terms, is decomposed into the deformation energy associated to changes in the electron density, E_def, and the deformation energy associated to changes in the exchange–correlation density, E_xc–def. Since the key quantity in the interaction of aromatic rings with carbon surfaces is the dispersion energy, for this work we have calculated the part of polarization energy corresponding to induction according to an expression derived from the second-order Rayleigh–Schröedinger perturbation theory (see the Appendix). Thus, the remaining part of the polarization energy contains the dispersion component. Although the induction energy calculated in this work is not the same as that obtained by SAPT due to the different treatment of the exchange corrections and intramolecular energy terms, it will be shown in the results and discussion section that this term becomes a revealing element for understanding the differences between implicitly and explicitly dispersion corrected DFT functionals. Thus, in this paper we will analyze the interaction energy according to the following decomposition scheme,

E_int = E_elec + E_Pauli + E_ind + E_res-pol (2)

where the total polarization energy is divided into induction, E_ind, and rest of polarization energy, E_res-pol, a term that mainly contains dispersion energy.

The structures of the complexes analyzed in this work are shown in Fig. 1, in all cases we have performed an exhaustive analysis to localize all the equilibrium structures and identify the most stable conformation. Only the most stable conformation has been considered in this work. Fig. 1 only shows the complex formed by pyridine with the largest graphene sheet model C150 (C150 means that the surface is formed with one hundred and fifty carbon atoms). Although the complexes formed with the remaining models (C24, C54 and C96) are not shown in the figure, all of them were optimized at the M06-2X/6-31G(d,p) level and its interaction energy, including the energy decomposition terms, were corrected for the basis set superposition error using the counterpoise method. For the sake of comparison with the finite graphene models, we have also calculated the structure and the interaction energy of the complex formed by pyridine adsorbed on a supercell of 72 carbon atoms using peridodic boundary conditions (C72-PBC). The structure of this complex is shown in Fig. 2. Calculation of total interaction energies and electron densities of the complex and the monomer was performed using the Gaussian 09 program package.²⁷ The energy decomposition analysis was performed for the adsorption complexes formed with the finite graphene models using an own code. This energy decomposition scheme allows for the calculation of the deformation electron densities associated to the Pauli and polarization components separately. These deformation densities were visualized using the Gaussview 5.0 interface program.²⁸ The empirical dispersion corrections of Grimme, denoted as GD3-corr in Tables 1 and 2, were computed using the DFT-D3 program.²⁹ A total of five functionals were compared; B97-D, BLYP-D, B3LYP-D, wB97XD and M06-2XD. For comparison, we have also performed Hartree–Fock calculations including the corresponding dispersion correction (HF-D). Finally, a Bader's topological analysis of the electron density³⁰ in the adsorption complexes was performed with the help of the AIM2000 software,³¹ characterizing the bond critical points (BCPs), ring critical points (RCPs) and cage critical points (CCPs) established between the atoms from pyridine and the carbon atoms from the surface.

Fig. 1 Optimized structures of the main adsorption complexes investigated.

Fig. 2 Optimized structure of the adsorption complex of pyridine on a supercell of 72 carbon atoms calculated using periodic boundary conditions.

Table 1 Interaction energy components (in kcal mol⁻¹) for the adsorption of pyridine on different graphene surface models obtained at different computational levels

	HF-D	BLYP-D	B97-D	B3LYP-D	wB97X-D^a	M062X-D
a This functional uses an own version of DFT-D2 correction instead of DFT-D3.
C24
Electrostatic	−7.75	−7.88	−7.01	−7.35	−6.47	−6.72
Pauli	21.31	22.52	21.43	21.85	20.74	20.78
Exchange	−31.66	−31.51	−29.35	−30.61	−28.95	−29.48
Repulsion	52.97	54.03	50.78	52.45	49.68	50.26
Polarization	−2.11	−3.61	−2.96	−5.70	−10.00	−19.15
Induction	−0.65	−1.13	−1.19	−1.05	−1.09	−1.07
Res-pol	−1.46	−2.48	−1.77	−4.65	−8.91	−18.08
GD3-corr	−16.92	−18.09	−17.49	−14.55	−11.45	−1.33
Total	−5.46	−7.06	−6.04	−5.76	−7.18	−6.43
C54
Electrostatic	−9.54	−9.27	−8.41	−8.81	−8.02	−8.31
Pauli	22.88	24.11	22.94	23.39	22.19	22.25
Exchange	−32.93	−32.89	−30.59	−31.92	−30.14	−30.77
Repulsion	55.82	57.00	53.53	55.31	52.33	53.02
Polarization	−2.28	−3.62	−2.99	−5.91	−10.52	−20.46
Induction	−0.75	−1.30	−1.37	−1.20	−1.25	−1.21
Res-pol	−1.53	−2.32	−1.62	−4.71	−9.27	−19.25
GD3-corr	−19.95	−21.18	−20.33	−17.34	−13.92	−2.14
Total	−8.90	−9.96	−8.79	−8.66	−10.27	−8.66
C96
Electrostatic	−9.65	−9.47	−8.59	−8.99	−8.19	−8.50
Pauli	22.98	24.23	23.05	23.51	22.28	22.34
Exchange	−32.83	−32.80	−30.48	−31.81	−30.03	−30.64
Repulsion	55.81	57.03	53.53	55.32	52.31	52.97
Polarization	−2.39	−3.59	−3.00	−5.91	−10.54	−20.42
Induction	−0.81	−1.35	−1.42	−1.26	−1.31	−1.28
Res-pol	−1.58	−2.24	−1.58	−4.65	−9.23	−19.14
GD3-corr	−20.29	−21.50	−20.61	−17.72	−14.33	−2.55
Total	−9.35	−10.32	−9.15	−9.10	−10.78	−9.13
C150
Electrostatic	−8.81	−8.75	−7.92	−8.27	−7.53	−7.82
Pauli	20.68	21.74	20.70	21.10	20.00	20.04
Exchange	−30.00	−29.91	−27.79	−29.00	−27.33	−27.91
Repulsion	50.68	51.66	48.49	50.10	47.34	47.95
Polarization	−2.08	−2.75	−2.37	−5.05	−9.54	−18.94
Induction	−0.72	−1.19	−1.25	−1.11	−1.16	−1.16
Res-pol	−1.36	−1.56	−1.12	−3.94	−8.38	−17.78
GD3-corr	−19.92	−21.06	−19.99	−17.42	−14.22	−2.70
Total	−10.13	−10.83	−9.57	−9.65	−11.29	−9.41

---

Table 2 Interaction energy components (in kcal mol⁻¹) for the adsorption of pyridine on C60 fullerene, and (5,5), (6,6) and (7,7) armchair single-walled carbon nanotubes obtained with the M06-2XD functional

	C60	(5,5)	(6,6)	(7,7)
Electrostatic	−4.88	−5.88	−5.98	−6.11
Pauli	13.09	16.16	16.40	16.88
Exchange	−18.28	−22.43	−22.67	−22.85
Repulsion	31.37	38.59	39.07	39.74
Polarization	−11.03	−14.37	−14.65	−15.24
Induction	−1.02	−0.85	−0.87	−0.87
Res-pol	−10.01	−13.52	−13.78	−14.37
GD3-corr	−1.30	−1.94	−1.98	−2.02
Total	−4.11	−6.02	−6.22	−6.48

---

Results and discussion

Graphene sheet modelling

As mentioned in the previous section, the graphene sheet models employed in this work comprise quasiplanar structures of carbon atoms. The smallest sheet considered corresponds to the coronene molecule (C24), which contains twenty-four carbon atoms disposed in a D_6h symmetrical structure. The following models were constructed by extending the perimeter of the coronene with additional benzene rings and keeping the D_6h symmetry. Thus, addition of benzene rings to the perimeter of coronene leads to the circumcoronene (C54) structure, subsequent additions leads to structures of ninety-six (C96) and one hundred and fifty (C150) carbon atoms. Only the optimized structure of the complex formed by pyridine and C150 is shown in Fig. 1.

As can be observed for the C150 model in the figure, graphene sheets are distorted with respect to planarity due to the interaction with pyridine. In order to account for the magnitude of the geometry distortion experienced by the four graphene models, we have obtained, by performing a multiple linear regression, the equation of the plane which fits better with the surface formed by the carbon atoms and calculated the residuals, ΔZ. The fitted equation in all cases has the following form,

z = ax + by + c (3)

For getting the equation above we first remove the pyridine molecule and then re-orientate the Cartesian axis so that the z axis lies perpendicular to the carbon surface and the origin of coordinates coincides with the mass center. To obtain a quantitative measure of the deviation of the carbon surface from the planarity, we represented the averaged ΔZ values calculated for the carbon atoms approximately equidistant from the z axis versus the distance to this axis (distance denoted by R in the figure). Thus, a representation whose slope is zero indicates similar residuals for all the atoms and would be found for a perfect planar surface, whereas a representation with a slope different from zero is found for a non-planar surface, that is, the larger the slope the larger the curvature. These representations for C54, C96 and C150 are shown in Fig. 3 and compared to that obtained for C72-PBC. As one can see the residuals and the distance display a very good linear correlation, as would be expected for quasi ‘bowl-shaped’ surfaces. Only in C150 the most internal carbons display an appreciable deviation from linearity. However, an apparently unexpected result is found. Thus, C150 appears to be the most distorted structure but C96 is a less distorted structure than C54 (see the slopes in Fig. 3). As will be shown later, this result has a direct effect on the interaction energy components. Comparing the finite models with the representation obtained for C72-PBC one can see that the slope, and then the curvature of the surface, found for C96 is the most similar to that obtained introducing periodic boundary conditions. However, this result must be taken with care as the inclusion of periodicity introduces the effect of adjacent molecules adsorbed on the surface. It is expected that between two adjacent molecules the graphene surface exhibits a change in its curvature giving rise to a ‘waveform’ of the surface. This is clearly appreciated in our model of 72 carbons, where distant atoms from the adsorbed molecule in the supercell show a negligible curvature (almost a slope zero in the representation of Fig. 3). This ‘waveform’ of the surface and its curvature will depend of course on the supercell size employed. To model experimental structures, the most appropriate size for the supercell will depend on the conditions of the experiment that lead to a specific concentration of pyridine on the surface. In this work, we have employed a supercell large enough to include the interaction with the most relevant atoms and sufficiently small to make possible the calculation within a reasonable time. The fact that truncated graphene sheets display a uniform curvature along the surface (a constant slope in the representation of Fig. 3) is mainly due to the absence of these near adsorbed molecules and also the introduction of terminal hydrogen atoms.

Fig. 3 Plots of the curvature of the fifty-four central carbon atoms in the C54 (upper left), C96 (upper right) and C150 (bottom left) graphene models. The scale represents the distance in Å of the surface with respect to the plane formed by the most external atoms in each model. Plot at the bottom right corner represents the residuals, ΔZ, obtained from multivariate linear regression performed with the Cartesian coordinates of the carbon atoms in each graphene model to get the approximated equation of the surface plane (see text for details), R represents the distance of the atoms to the Z axis (perpendicular to the plane). In this plot the representation obtained for the supercell of 72 carbon atoms and PBC (C72-PBC in the plot) is included for the sake of comparison. The positions of the carbon and nitrogen atoms from the pyridine molecule are also indicated by white and blue open circles, respectively.

In Fig. 3 we have also included plots that represent the depth of the central surface formed by the inner circumcoronene unit (54 carbon atoms). The darkness in the plots is proportional to the depth at each point and is referred to the most external carbon atoms. These plots give us an idea of the relative curvature of the different structures as well as the most curved regions. As one can see the most curved structure corresponds to C150, whereas the least curved one is C96. Moreover, these plots clearly reflect the “footprint” of the adsorbed pyridine. Thus, the darkest region involves the carbon atoms that interact more strongly with the molecule in all cases.

Interaction energies

The interaction energies of the graphene model surfaces and pyridine are collected in Table 1. The different energy components, electrostatic, Pauli, induction and remaining polarization, are also shown in the table. As one can see, the interaction energy (its absolute value) increases with the size of the carbon surface. The experimental adsorption energy of pyridine on graphite, −9.3 kcal mol⁻¹,³² is close to the interaction energy values obtained with C96 and C150 models. On the other hand, an estimated value of −8.5 kcal mol⁻¹ for the adsorption energy of pyridine on graphene was obtained by correlation of MM2 theoretical adsorption energies obtained for 118 molecules on 3-layer and 1-layer graphene with experimental adsorption energies of the same molecules on graphite.³³ M06-2XD adsorption energies, which also include the geometry deformation energy of the molecule and surface, for C96 (−7.0 kcal mol⁻¹) and C150 (−7.2 kcal mol⁻¹) models are reasonably close to this estimated value. It must be pointed out that no thermal corrections have been added to the adsorption energies given above. The use of the harmonic model on large complexes with many small intermolecular frequencies may lead to significant errors when calculating the binding energy.

On the other hand, the interaction energies obtained for C72-PBC with the M06-2X functional are −6.79 kcal mol⁻¹ and −9.09 kcal mol⁻¹ when the DFT-D3 correction is included. These values are very close to those obtained for the finite models of C96 (−6.58/−9.13 kcal mol⁻¹) and C150 (−6.71/−9.41 kcal mol⁻¹) at the same level of theory. The above results suggest that both C96 and C150 models are suitable to reproduce reasonably well the total interaction energy of pyridine on graphene. However, the interaction energy components reflect some significant differences in the complexes formed with C96 and C150. Thus, the largest curvature experienced by C150 with respect to C96, which was analysed in the previous section, leads to larger distances between the atoms of pyridine and the central atoms of the surface. As an example, the distance between the pyridine nitrogen and the geometrical mean of the central benzene ring is 3.17 Å in C54, 3.13 Å in C96 and 3.24 Å in C150 (this distance is a good reference since the nitrogen is placed over the center of this benzene ring). As a consequence, the attractive energy terms, electrostatic and polarization, decrease in C150, nevertheless this decrease is overcompensated by the Pauli repulsion energy, which is also smaller. Only the dispersion correction term (its absolute value) calculated with M06-2XD increases in C150 and this term accounts for most of the energy difference between C96 and C150 at this computational level (0.28 kcal mol⁻¹). In summary, excluding the dispersion correction, C54, C96 and C150 show very similar interaction energies but the relative strength of the different energy components differs due to the different distances to the surface. This also causes small changes in the electron density deformation when using different graphene models.

One of the most important points in this work is to compare the performance of the different functionals. Looking at the interaction energy components, M06-2XD is found to be the only functional that predicts an energetically stable complex without including the dispersion correction (GD3-corr). This is because the largest part of the dispersion energy is already accounted for by the functional itself. The second smallest value for the dispersion correction is displayed by the wB97XD functional, which is followed by B3LYP-D by more than 3 kcal mol⁻¹. Analysing term by term one discovers that the first order energy terms, electrostatic and Pauli, do not show significant differences among the functionals considered. Furthermore, as expected, the largest electrostatic energies are found for the HF-D method, followed by BLYP-D. These methods also display the largest repulsion energies, but in this case the values are higher in BLYP-D than in HF-D.

On the other hand, remarkable differences among the functionals are found in the polarization energies. The magnitude of these energies indicates the degree of dispersion included implicitly by the functionals. As expected, the smallest polarization is always found for HF-D as none of the dispersion is included within the HF method. Within the DFT functionals the magnitude of the polarization energy term follows the order: B97-D < BLYP-D < B3LYP-D < wB97XD < M06-2XD. The difference between M06-2XD and the remaining functionals is enormous; thus, the M06-2XD polarization energy for all the models studied is around twice that obtained with wB97XD, which in turn is also around twice that obtained by B3LYP-D. To understand the origin of these differences we have inspected the polarization energy by separating the part of induction from the rest of polarization, the latter including the part of dispersion. Although all the DFT functionals display quite similar induction energies, the largest values are always obtained for those displaying also the smallest total polarization energy (B97-D and BLYP-D). On the other hand, HF-D provides the smallest induction energy. All of this indicates that the large difference between M06-2XD and the other functionals analysed is concentrated on the term denoted by E_res-pol in eqn (2). In other words, this functional accounts for most of the dispersion energy implicitly, improving only the term directly related to this energy and leaving the remaining terms unaltered in comparison with the other functionals. This result deserves to be analysed in more detail using a larger set of adsorption complexes in order to reach a general conclusion. Also, a study of small size dispersion complexes for which high-level CCSDT calculations are available and perturbation methods for the calculation of dispersion energy may be applied would provide more solid arguments for a discussion about this question.

Even though the interaction energy components reflect a good trend of the M06-2XD functional at an energetic level, the effect that dispersion exerts on the electron density is quite important in the calculation of properties derived from this function. As mentioned in the Introduction section, the functional must be able to account for the dispersion effect on the electron density in order to be suitable for the determination of changes experienced in charge distributions, bond orders, polarizabilities or DFT reactivity indices upon adsorption. In Fig. 4 the electron polarization density of pyridine on C150 is represented using HF-D, B3LYP-D, wB97XD and M06-2XD. Only small differences are appreciated between B3LYP-D and wB97XD densities and even between these and the HF-D one. Notice that the electron density is not modified by the inclusion of the dispersion correction as it is introduced after the self-consistent field procedure. In contrast, the accumulation of electron density between pyridine molecule and the surface (represented using magenta colour) significantly increases using M06-2XD. This is the effect of the dispersion interaction, which tends to concentrate the electron density in the intermolecular region in a dispersion complex.²⁶ In order to shed light on this issue, we have calculated the electron polarization density for the smallest graphene model (C24) at the post-HF level using MP2 theory on the geometry obtained with the M06-2X functional. In addition, we have also optimized the complex formed between pyridine and benzene (C6) at the M06-2X level and calculated the polarization density using the different functionals and the MP2 and CCSD methods. The results obtained are plotted in Fig. 5. As one can see there are significant differences in the amount of electron density accumulated in the intermolecular region between HF-B3LYP-B97 and CCSD-MP2-M062X. Although the distribution of this electron density differs between MP2, CCSD and M06-2X, the three methods agree in magnitude in contrast to the remaining methods. These results indicate that electron density properties may be ill-defined in highly dispersive complexes or molecules adsorbed on carbon surfaces when dispersion is only accounted for via explicit empirical corrections, whereas functionals including dispersion implicitly (like M06-2X) may provide a proper solution to this problem.

Fig. 4 Electron polarization densities obtained using the C150 graphene model as a substrate with HF-D (upper left), B3LYP-D (upper right), wB97XD (bottom left) and M06-2XD (bottom right). The density is shown for an isosurface value of 2 × 10⁻⁴ au.

Fig. 5 Electron polarization density obtained using the C24 graphene model (bottom plots) and the benzene molecule (top plots) with different computational levels. The density is shown for an isosurface value of 2 × 10⁻⁴ au.

Adsorption on carbon nanotubes and fullerenes: effect of the surface curvature

We have finally explored the effect of the surface curvature on the adsorption strength. This is a quite important point when considering carbon structures as substrates for surface enhanced Raman scattering experiments. The chemical enhancement in the so-called graphene enhanced Raman scattering (GERS) technique may be significantly affected by the curvature of the substrate, both by changes in the molecule–surface interaction strength and resonance effects induced by the electromagnetic excitation.

In this work we have then explored the interaction of pyridine with the surface of cylindrical carbon structures (armchair single-walled carbon nanotubes) and spherical ones (fullerenes). In Table 2, the total interaction energy and its components are collected for calculations performed with the M06-2XD functional. As previously found in ref. 21 for carbon nanotubes the larger the curvature the smaller the interaction energy. The one-dimensional curvature in nanotubes makes the interaction energy decrease between 3 and 3.4 kcal mol⁻¹ with respect to C150. As expected the interaction energy increases when going from (5,5) to (7,7) since the curvature decreases. On the other hand, the two-dimensional curvature in C60 fullerene reduces even more the molecule–surface contact and the interaction energy decreases around 2 kcal mol⁻¹ with respect to the armchair (5,5) substrate. The interaction energy components less affected by the curvature of the surface correspond to the electrostatic and induction terms. In the case of the induction energy it depends on both the electrostatic field created by the unperturbed monomers and the magnitude of the multipoles induced by this electrostatic field (see Appendix A for description of this energy term). These induced multipoles depend strongly on the polarizability of the adsorbed molecule and the surface fragment interacting with it. The magnitude of the electrostatic field is reflected on the values of the electrostatic energy and depends on the distance and electronic distribution of the unperturbed molecules. This electrostatic energy is significantly larger in the graphene sheet models and consequently the magnitude of the induction energy also increases with respect to the curved surfaces. Among the curved surfaces, the induction energy is larger in C60 than in carbon nanotubes, even that it displays a smaller electrostatic energy. In this case the explanation must be found in the magnitude of the multipoles induced on the surface by the electrostatic field created by the molecule, which must be larger in the spherical (two-dimensional curved) geometry of C60.

The most significant difference among graphene, nanotubes and fullerene models is found in the polarization energy term that contains the dispersion energy (res-pol in Tables 1 and 2). These differences are also reflected in the electron polarization density (Fig. 6), whose magnitude decreases from C150 (Fig. 4) to C60 fullerene with nanotubes in between. The electron density topology (Fig. 7) also reflects a decrease in the number of intermolecular bond, ring and cage critical points when going from C150 to C60. Thus, only two BCPs are established between the ‘ortho’ carbon atoms of pyridine and the carbons from C60 fullerene, whereas up to three BCPs are found with the nanotubes, also involving the interaction of the ‘para’ carbon atom of pyridine. Although the topology obtained for C150 also shows three BCPs, the number of RCPs (six in this case but only three in nanotubes) and the number of CCPs (four in this case but only one in nanotubes) reveals that the interaction also involves the ‘meta’ carbon atoms and the nitrogen atom from pyridine. It must be remarked that all the topologies determined fulfil the Poincaré–Hopf relationship.³⁰ The number of ‘contacts’ between atoms from the pyridine and atoms from the surface characterized by the topology of the electron density contributes to explain the relative strength of the interaction along the different carbon substrates.

Fig. 6 Electron polarization densities obtained with M06-2XD using C60 fullerene (left) and armchair (7,7) single-walled carbon nanotube (right) as substrates. The density is shown for an isosurface value of 2 × 10⁻⁴ au.

Fig. 7 BCPs (red points), RCPs (yellow points) and CCPs (green points) linking pyridine with the C150 graphene model (top), C60 fullerene (bottom left) and armchair (7,7) single-walled carbon nanotube (bottom right).

Summary

In view of the interest emerged in the last few years on the use of carbon substrates as chemical sensors we have performed in this work a theoretical study of the adsorption of an aromatic molecule (pyridine) on different carbon allotrope models, graphene, nanotubes and fullerenes. As the main source of stability for the adsorption complexes formed is the dispersion interaction, we have investigated the performance of different DFT functionals that incorporates this component implicitly or via empirical corrections. An interaction energy decomposition scheme has allowed analysing the differences in the interaction energy components among the different functionals. Thus, we have found that the implicitly dispersion-corrected functional M06-2X provides similar values of the electrostatic, Pauli and induction energies to those obtained with other functionals, but the M06-2X polarization energy term that includes dispersion is substantially different to the others. The wB97X functional, and in less degree the B3LYP one, also correct part of the dispersion energy implicitly, but it is not enough to predict a stable complex without empirical corrections. In contrast, B97 and BLYP functionals only correct marginally the HF result in terms of dispersion. The inclusion of an empirical dispersion correction improves the calculated interaction energy but does not change the electron polarization density. In contrast, M06-2X displays significant changes in the polarization density, predicting a larger accumulation of electron density in the intermolecular region than other functionals. As a consequence, problems in the calculation of electron density derived properties such as charge distributions, bond orders or electron polarizabilities may emerge in highly dispersive adsorption complexes when dispersion is only corrected explicitly. Finally, the study of the pyridine adsorption on nanotubes and fullerenes indicates that one-dimensional curvatures (nanotubes) and in larger degree two-dimensional curvatures (fullerenes) mainly decrease the polarization energy, destabilizing the adsorption complex with respect to quasi-planar graphene model surfaces. Polarization density plots and topological analysis of the electron density reflect that the surface curvature reduces the number and strength of the ‘contacts’ established between the pyridine atoms and the surface atoms.

Appendix

A clear-cut interpretation of the interaction energy in intermolecular force theory is given by the Rayleigh–Schrödinger perturbation theory, RSPT,^34–37 and symmetry adapted perturbation theories, SAPT,³⁸ the latter including exchange terms. In this section we will compare the interaction energy partitioning obtained using deformation densities²⁶ to that obtained using 2nd-order RSPT.³⁷ In order to simplify the task we will neglect the exchange energy terms by considering the systems at large distance. Since in the deformation density partitioning there is a clear correspondence between energy terms and exchange–correlation corrections the conclusions obtained in this section for non-overlapped systems can be extrapolated to the overlapped case.

In most of the cases 1st- and 2nd-order energies in the RSPT are enough to account for most of the interaction energy. They are associated to electrostatic and polarization energies, respectively, and the polarization term can be in turn split up into pure induction and pure dispersion energies.³⁷ These are represented by eqn (A1)–(A3) for the interaction of two monomers A and B, where the super index 0 denotes the ground state.

The perturbation energy is written in terms of the ground state one-electron densities, ρ⁰⁰_A and ρ⁰⁰_B, and the induced transition one-electron densities, ρⁿ⁰_A and ρ^m0_B. The summations in eqn (A2) and (A3) run over all the states obtained by intramonomer excitations of the electrons from occupied to unoccupied molecular orbitals in molecules A and B. The operators in the first two terms of eqn (A1) correspond to the nuclear potential of nuclei A and B, respectively, whereas those of eqn (A2) also incorporate the unperturbed electron densities of monomers A and B, respectively.

It is obvious that the electrostatic terms in eqn (A1) and eqn (12) of ref. 26 are identical since both represent the intermolecular coulombic interaction between the nuclei and electrons of unperturbed monomers. So, the polarization terms (induction plus dispersion in the 2nd-order RSPT) must be equivalent to the polarization energy of eqn (15) in ref. 26, excluding the intermolecular exchange–correlation part.

The key point here is to determine how induction and dispersion energies are accounted for in the deformation density partitioning scheme. As deformation density partitioning contains intermolecular and intramolecular terms and the RSPT gives rises to only intermolecular terms, these two methodologies could seem to be a priori non-comparable. However, the only requirement for the deformation density scheme to accurately describe the polarization energy is the knowledge of the deformation density at a good precision independently of the quantum chemical method employed for its calculation. Thus, one can calculate the deformation density up to 2nd-order (in energy) using the RSPT and then introduce the result within the expressions obtained in the previous section in order to merge both methodologies.

In order to get the deformation density induced by the intermolecular interaction using 2nd-order RSPT first we need to establish the order of correction required for the wave function. Since we intend to go up to 2nd-order in energy, the corresponding wave function requires only 1st order correction. Using 1st order correction to the wave function the perturbed density, ρ′, at a point r, can be calculated through eqn (A4),

where Ψ⁽¹⁾ and Ψ⁽⁰⁾ represent the 1st-order correction and the unperturbed wave functions, and [small rho, Greek, circumflex](r⃑) is the charge density operator, which is given by:where n_A and n_B are the number of electrons of monomers A and B, respectively.

The normalization factor 1/N² is given by eqn (A6).

We have employed the intermediate normalization of the wave function, 〈Ψ⁽⁰⁾|Ψ⁽⁰⁾〉 = 1 and 〈Ψ^(k)|Ψ⁽⁰⁾〉 = 1 for ∀k ≠ 0. The first term in eqn (A4) is the unperturbed density so that the remaining terms constitute the deformation density.

The unperturbed wave function is constructed with the product of the wave functions for the isolated molecules. For the case of two non-overlapped systems A and B the unperturbed wave function is just given by the Hartree product of the wave functions of A and B in their ground states.

Ψ⁽⁰⁾ = ψ⁰_Aψ⁰_B (A7)

On the other hand, the 1st-order correction to the wave function is constructed by summation of all the accessible excited states obtained by intramonomer excitations of electrons from occupied to unoccupied unperturbed MOs, where each excited state is multiplied by a coefficient that mainly depends on the energy gap between the ground and the excited state.

The operator V̂ represents the perturbation, which in this case corresponds to the intermolecular terms included in the Hamiltonian given by:

It must be noticed that the third term in eqn (A4) comes from the self-product of Ψ⁽¹⁾. However, looking at the form of Ψ⁽¹⁾ one can realize that this term contains only 3th-order corrections to the energy. Since only 1st- and 2nd-order energies need to be considered, this term can be removed from the expression of ρ′. By subtracting the unperturbed density, the deformation density up to 2nd-order (in energy) is obtained.

Δρ(r⃑) = [〈Ψ⁽⁰⁾|[small rho, Greek, circumflex](r⃑)|Ψ⁽¹⁾〉 + 〈Ψ⁽¹⁾|[small rho, Greek, circumflex](r⃑)|Ψ⁽⁰⁾〉] (A10)

It must be noticed that eqn (A10) is actually the 1st-order correction to the electron density derived by Ángyán.³⁹ Expanding eqn (A10) using eqn (A7) and (A8) one obtains

Now, perturbation theory and deformation density schemes can be merged introducing eqn (A11) into eqn (15) of ref. 26. However, it is better to split first the polarization energy, E_pol, of ref. 26, excluding the exchange–correlation part, into that corresponding to only deformation density, E_def, and that merging deformation density with unperturbed densities, E_elec-def, obtaining the following expressions,

By comparison with eqn (A2) the induction energy is located within the term E_elec-def.

It must be noticed that the last two terms of eqn (A12) and those of eqn (A13) represent intramonomer energy changes in the deformation density partitioning. Since eqn (A11) can be unambiguously split into deformation densities of monomers A and B by separating their corresponding induced transition one-electron densities, E_ind can be expressed as follows.

eqn (A15) relates the energy term E_elec-def of the deformation density partitioning with the induction energy obtained from 2nd-order RSPT. Extrapolating this equation to any computational level one can calculate the contribution of the induction energy provided that deformation densities are known.

In addition, eqn (A15) can also be rewritten as,

allowing connection with the classical view of the induction energy. Thus, comparing the first two terms of the electrostatic energy of eqn (12) in ref. 26 (electron–nuclei electrostatic interaction) with eqn (A16), one can see how the expression of the induction energy is half that of the electrostatic energy. It is assumed in classical theory of intermolecular forces that half of the energy involved in the electrostatic interaction between induced multipoles and the electric field responsible of the induction is employed in the induction process. Thus, eqn (A16) is in agreement with this classical view of induction energy.

The final requirement for the application of eqn (A16) is to define the monomers’ domains for the decomposition of the deformation density into A and B contributions. This can be done by introducing a real space atomic partitioning of the deformation density or the Mulliken-like atomic partitioning of the Hilbert space of the basis functions.⁴⁰ Real space atomic partitionings are computationally much more expensive and for that reason the Mulliken-like scheme was chosen in this work.

Acknowledgements

The authors thank ‘Ministerio de Economía y Competitividad’ for funding their research through project CTQ2010-20229 and Centro de Supercomputacion de Galicia (CESGA) for providing access to its computational facilities.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
2. D. R. Dreyer, R. S. Ruoff and C. W. Bielawski, Angew. Chem., Int. Ed., 2010, 49, 9336–9344.
3. A. K. Geim, Science, 2009, 324, 1530–1534.
4. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
5. K. S. Novoselov, V. Fal, L. Colombo, P. Gellert, M. Schwab and K. Kim, Nature, 2012, 490, 192–200.
6. S. Stankovich, D. A. Dikin, G. H. Dommett, K. M. Kohlhaas, E. J. Zimney, E. A. Stach, R. D. Piner, S. T. Nguyen and R. S. Ruoff, Nature, 2006, 442, 282–286.
7. F. Schedin, A. K. Geim, S. Morozov, E. Hill, P. Blake, M. Katsnelson and K. Novoselov, Nat. Mater., 2007, 6, 652–655.
8. V. Georgakilas, M. Otyepka, A. B. Bourlinos, V. Chandra, N. Kim, K. C. Kemp, P. Hobza, R. Zboril and K. S. Kim, Chem. Rev., 2012, 112, 6156–6214.
9. X. Ling, L. Xie, Y. Fang, H. Xu, H. Zhang, J. Kong, M. S. Dresselhaus, J. Zhang and Z. Liu, Nano Lett., 2009, 10, 553–561.
10. W. Xu, N. Mao and J. Zhang, Small, 2013, 9, 1206–1224.
11. M. Fleischmann, P. Hendra and A. McQuillan, Chem. Phys. Lett., 1974, 26, 163–166.
12. S. Nie and S. R. Emory, Science, 1997, 275, 1102–1106.
13. K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan, R. R. Dasari and M. S. Feld, Phys. Rev. Lett., 1997, 78, 1667.
14. H. Xu, E. J. Bjerneld, M. Käll and L. Börjesson, Phys. Rev. Lett., 1999, 83, 4357.
15. A. M. Michaels, M. Nirmal and L. Brus, J. Am. Chem. Soc., 1999, 121, 9932–9939.
16. A. Aguiar, S. Fagan, L. da Silva, J. M. Filho and A. Souza Filho, J. Phys. Chem. C, 2010, 114, 10790–10795.
17. J. D. Wuest and A. Rochefort, Chem. Commun., 2010, 46, 2923–2925.
18. O. V. Ershova, T. C. Lillestolen and E. Bichoutskaia, Phys. Chem. Chem. Phys., 2010, 12, 6483–6491.
19. E. Voloshina, D. Mollenhauer, L. Chiappisi and B. Paulus, Chem. Phys. Lett., 2011, 510, 220–223.
20. C. Thierfelder, M. Witte, S. Blankenburg, E. Rauls and W. G. Schmidt, Surf. Sci., 2011, 605, 746–749.
21. D. Umadevi and G. N. Sastry, J. Phys. Chem. Lett., 2011, 2, 1572–1576.
22. Y.-n. Guo, X. Lu, J. Weng and Y. Leng, J. Phys. Chem. C, 2013, 117, 5708–5717.
23. P. Lazar, F. E. Karlický, P. Jurečka, M. S. Kocman, E. Otyepková, K. R. Šafářová and M. Otyepka, J. Am. Chem. Soc., 2013, 135, 6372–6377.
24. S. Grimme, WIREs Comput. Mol. Sci., 2011, 1, 211–228.
25. K. Szalewicz, WIREs Comput. Mol. Sci., 2012, 2, 254–272.
26. M. Mandado and J. M. Hermida-Ramón, J. Chem. Theory Comput., 2011, 7, 633–641.
27. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009.
28. R. Dennington, T. Keith and J. Millam, GaussView, Version 5, Semichem Inc., Shawnee Mission KS, 2009.
29. S. Grimme, J. Antony, S. Ehrlich and H. J. Krieg, Chem. Phys., 2010, 132, 154104.
30. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, 1990.
31. F. Biegler-König, J. Schönbohm and D. Bayles, J. Comput. Chem., 2001, 22, 545–559.
32. S. N. Yashkin, D. A. Svetlov and A. K. Buryak, Russ. Chem. Bull., 2003, 52, 344–353.
33. J. H. Son and T. R. Rybolt, Graphene, 2013, 2, 18–34.
34. E. Schrödinger, Ann. Phys., 1926, 80, 437–490.
35. J. O. Hirschfelder, Chem. Phys. Lett., 1967, 1, 325.
36. J. O. Hirschfelder, Chem. Phys. Lett., 1967, 1, 363.
37. A. J. Stone, The theory of intermolecular forces, Clarendon Press, Oxford University Press, Oxford, U. K, 1997.
38. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887.
39. J. G. Ángyán, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 78, 022510.
40. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833.

---