Adsorption of Organic Electron Acceptors on Graphene-like Molecules

Article pubs.acs.org/JPCC

Adsorption of Organic Electron Acceptors on Graphene-like Molecules: Quantum Chemical and Molecular Mechanical Study Susanta Haldar,†,‡,∥ Michal Kolár,̌ †,‡,∥ Róbert Sedlák,†,‡ and Pavel Hobza*,†,§ †

Institute of Organic Chemistry and Biochemistry and Gilead Science Research Center, Academy of Sciences of the Czech Republic, Flemingovo nam. 2, 166 10 Prague 6, Czech Republic ‡ Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles University in Prague, Albertov 6, 128 43 Prague 2, Czech Republic § Regional Center of Advanced Technologies and Materials, Department of Physical Chemistry, Palacký University, Olomouc, 771 46 Olomouc, Czech Republic ABSTRACT: Various cluster models of graphene and a periodic graphene with two organic electron acceptors (tetracyanoethylene and tetracyanoquinodimethane) were investigated by means of several quantum chemical and molecular mechanical approaches. The benchmark interaction energies of the coronene complexes were calculated at the MP2.5/CBS/6-31G*(0.25) level of theory. The SCS-MI-MP2, BLYP-D3 and, surprisingly, also AMBER showed modest agreement in the absolute as well as relative interaction energies. Consequently, larger complexes were investigated at these lower levels of theory including also DFTB-D. Charge transfer was calculated on the basis of Mulliken and NBO analysis. A high correlation between the interaction energies and charge transfer was observed. Further, vibrational analysis of the complexes revealed the association free energies for the gas phase and aqueous environment at the DFTB-D and AMBER levels. Extensive potential of mean force molecular dynamics simulations were carried out for all of the graphene−organic acceptor complexes. The convergence with the graphene model size was observed for the interaction energies as well as for the association free energies, which justifies using a cluster graphene model when the periodic one is not accessible. The role of translational entropy loss upon binding and the solvent contribution were discussed thoroughly.

1. INTRODUCTION The graphene has attracted intense research interest since it was first experimentally discovered in 20041,2 as a fine semiconductor-based nanoelectronics.3 A large number of studies investigated the adsorption properties of the graphene surface, where the range of adsorbate particles spanned from metal ions to biomolecules such as DNA.4−6 It was soon realized that adding an organic electron acceptor molecule to graphene, which is an electron donor, can notably modify its electronic properties.7,8 P-type doping of the graphene might be achieved by adding such an electron acceptor onto the graphene surface; accordingly, the adsorption of an electron donor might result in n-type doping of the graphene. It was concluded that the extent of charge transfer depends on the quality of the adsorbate and the surface coverage.9−11 This tuning of the electronic transport properties is thus of fundamental importance, and deep insight into its nature is needed. Charge-transfer complexes are characterized by the transfer of electron density from an electron donor to an electron acceptor species, but the definition of this energy term is not unambiguous. Surprisingly, the charge-transfer energy term is missing among the various perturbation energy contributions as © 2012 American Chemical Society

defined in the symmetry adapted perturbation treatment (SAPT).12 Here, the total interaction energy between two subsystems consists of the systematically repulsive exchangerepulsion energy, the systematically attractive induction and the dispersion terms and their repulsive exchange parts, and the electrostatic term, which might be attractive as well as repulsive. The charge-transfer term is covered in the induction energy, and there is a technique to calculate this energy term.13 The charge-transfer energy term is included in the supersystem interaction energy, determined as the difference between the energy of a complex (supersystem) and the energies of interacting subsystems. When the supersystem interaction energy is determined accurately, the charge-transfer term is accurately described as well. The charge-transfer energy term is no doubt an important stabilization contribution, but the dispersion and electrostatic energy contributions are important as well. The theoretical methods applied should thus properly cover all these energy terms. Let us mention here that neither Received: July 18, 2012 Revised: October 4, 2012 Published: November 13, 2012 25328

dx.doi.org/10.1021/jp3071162 | J. Phys. Chem. C 2012, 116, 25328−25336

The Journal of Physical Chemistry C

Article

Figure 1. Top views on the structures investigated.

energies at the highest WFT level available, and on the basis of these benchmark calculations we used lower DFT levels of theory for larger graphene models. Besides the interaction energies, we also calculated the thermodynamic characteristics. These characteristics are obtained by using umbrella sampling for the potential of mean-force (PMF) evaluation and inherently contain information on the dynamic behavior of the complex.20−22 For such extended complexes, the PMF technique can be combined only with molecular mechanics (MM) because the higher levels of theory are computationally inaccessible. Hence, it is first necessary to estimate the accuracy of the MM force field. Then, in the second step, the PMF/MM calculations in the gas phase as well as in an aqueous solution are performed. The PMF calculations even with MM are tedious due to the sampling requirements; therefore, in the third step, we evaluate the suitability of a simpler rigid rotor−harmonic oscillator−ideal gas (RR-HO-IG) approximation.23 In the PMF case, apart from the cluster models the periodic graphene was also investigated. This study aims to contribute to the following issues: (i) the calculation of interaction energies at the highest level available with a reliable description of the relevant interactions between the graphene model and the organic electron acceptor; (ii) the elucidation of the role of charge transfer in calculations employing cluster graphene models of increasing size and consequently the description of the convergence behavior or the energies and CT features; and (iii) the free energy characteristics and the lateral mobility of organic molecules with the discussion on the accuracy of the empirical force field and the applicability of the RR-HO-IG approach.

Hartree−Fock nor (standard) density functional theory (DFT) techniques cover dispersion energy. Several papers on the interaction of small organic molecules with graphene have already appeared.10,14 Often, the results are based on periodic graphene models treated by such DFT functionals as local density approximation (LDA), which however lack many important contributions mentioned in the previous paragraph. In our previous studies, we performed accurate interaction energy calculations for various chargetransfer complexes created from small electron donors and acceptors (e.g., NH3...BH3 or NH3...SO2)15 as well as from carbon aromatic systems and metallic atoms (e.g., benzene, graphene...X; X = Ag, Au, Pd).16 The carbon aromatic systems in the latter study were modeled by benzene and coronene as well as by periodic graphene. The periodic calculations performed by the VASP program package17 employing the DFT were, however, not accurate enough, and sufficiently accurate results were obtained when the respective DFT calculations were calibrated to the wave function theory (WFT) cluster calculations. It may be questioned to what extent the finite size of coronene affects the interaction energies or other complex properties that are important for the understanding of the binding. The adsorption of tetracyanoethylene was investigated previously by Lu et al.10 using a periodic model of graphene. Anota et al. theoretically studied the adsorption of amine group on coronene.18 Herein, the graphene surface was thus modeled by carbon aromatic molecules of increasing size, starting from coronene (C2), circumcoronene (C3), dicircumcoronene (C4), and ending with the C5 system, which is one layer of benzene rings larger than C4 (see Figure 1). As electron acceptors, tetracyanoethylene (TCNE) and tetracyanoquinodimethane (TCNQ) were adopted. Extensive evidence has been accumulated in the past years that the coupled-cluster method, covering the single- and double-electron excitations iteratively and the triple-electron excitations perturbatively (CCSD(T)) calculated at the complete basis set limit (CBS), provides benchmark interaction energies differing from the accurate energies by less than about 0.5 kcal/mol.19 Nevertheless, this level of theory is not applicable even for the smallest graphene finite model used herecoronene. Therefore, we calculated the interaction

2. COMPUTATIONAL METHODS 2.1. Systems Investigated. In the cluster calculations, the graphene sheet was modeled by polycyclic aromatic molecules with an increasing number of the aromatic rings. The structures of the C2, C3, C4, and C5 systems were prepared manually. For each CX (X = 2−5) and TCNY (Y = E, Q), the complexes, abbreviated as CXY, were prepared by stacking the electron acceptor to the center of the carbon aromatic molecule in two different structural orientations. In the C2 case, only one orientation was prepared. All of the complexes are depicted in Figure 1. An infinite graphene layer used for the MM 25329



Article

The use of diffuse atomic basis sets for such an extended molecular cluster usually means that the basis set is overcomplete. This might be connected with some numerical problems. To overcome this problem, we have decided to use a combined basis set for the description of the donor−acceptor complexes.36 Nonaugmented basis sets were used for all of the hydrogen atoms. Augmented (aug-cc-pVDZ or aug-cc-pVTZ) basis sets were used only on every second carbon atom and the corresponding nonaugmented basis sets for the rest of the carbon atoms. This means that there was no bond between the two carbon atoms that were described by the augmented basis set. The heavy atoms of the electron-acceptor molecule were described by the corresponding augmented basis set. The MP2.5/CBS/6-31G*(0.25) interaction energies are among the most accurate ones, and the use of MP2.5 as a benchmark in this paper is fully justified.33 We note that using any higher level of theory providing even more accurate interaction energies (e.g., CCSD(T)/CBS or SCS(MI)-CCSD/ CBS methods)19 is unfeasible for the extended systems presented. All the MP as well as the scaled MP interaction energy calculations were performed on BLYP-D3 optimized structures. The BLYP-D3, AMBER, and DFTB-D interaction energy calculations were conducted on the optimized structures of the corresponding level of theory. The DFT calculations were performed in the Turbomole 6.3 program package,37 whereas all of the MP calculations were performed in the Gaussian09 package.38 The DFTB-D and MM calculations were carried out in the DFTB+ and AMBER packages, respectively.39 2.3. Charge-Transfer Characteristics. The net charges of the CX, TCNE, and TCNQ were calculated in two different manners. In previous studies, there was not much attention paid to the method of atomic charge calculations, with Mulliken charges being used commonly.10,14 Mulliken charges, however, suffer from large basis-set dependence and should be used with caution. We further performed a natural bond orbital (NBO) analysis40 of all of the complexes since the NBO charges are less basis set dependent and belong to more reliable ones. The CT energy is within the NBO theory,40 approximated by the second-order E2 energy:

calculations was prepared employing periodic boundary conditions in two dimensions. 2.2. Interaction Energy. The interaction energy of a CT complex R...T (ΔE(RT)) is defined by eq 1: ΔE(RT) = E(RT) − E(R) − E(T)

(1)

where E(RT) stands for the total electronic energy of the complex and E(R) and E(T) are the electronic energies of the electron donor R and electron acceptor T, respectively. Here, only the interaction energies based on the MP2 method were corrected for the basis-set superposition error using the standard counterpoise procedure.24 The subsystem deformation energies were not included.25 The interaction energies for all of the complexes were determined at the following levels of theory: (i) dispersioncorrected density functional theory (DFT-D3), (ii) dispersioncorrected density functional tight binding (DFTB-D), and (iii) empirical force field. The DFT-D3 calculations were performed at the BLYP/TZVPP level with empirical correction accounting for the dispersion attraction.26 The DFTB-D calculations were performed with the corrections proposed by Elstner et al.27,28 Finally, the empirical-potential calculations were based on the general Amber force field (GAFF).29,30 For the GAFF calculations, the atomic partial charges were calculated according to the restricted electrostatic potential (RESP) approach.31 The RESP fit was performed onto a grid of electrostatic potential points calculated at the HF/6-31G* level as recommended by the force-field designers.29,31 For C2Y and C3Y (Y = E, Q), the benchmark calculations were performed at the MP2.5/CBS/6-31G*(0.25), SCS-MP2, and SCS-MI-MP2 levels of theory. The former interaction energy was calculated according to eq 2: ΔE(MP2.5/CBS) = ΔE(MP2/CBS) + 0.5[ΔE(MP3/6‐31G*(0.25)) − ΔE(MP2/6‐31G*(0.25))]

(2)

The MP2/CBS interaction energy was calculated by the extrapolation of aug-cc-pVDZ and aug-cc-pVTZ basis sets (see below) while the MP3 correction term was determined with the 6-31G*(0.25) basis set, where the number in the parentheses represents the exponent of the d-function of the heavy (carbon and nitrogen) atoms.32 The MP2.5/CBS/6-31G*(0.25) interaction energies were found to be highly accurate, and their root-mean-square deviation toward S66 interaction energies is very small, 0.21 kcal/mol.33 The use of the spin-component scaled MP2 method (SCSMP2)34 is based on the separate scaling of the same and opposite components of the correlation energy. Originally, the respective coefficients were calibrated on the basis of reaction rates, and the overall performance of the method for noncovalent interactions was better than that of the plain MP2 method (see below). When however the calibration of these two coefficients was based on molecular interactions (MI), the performance of the SCS(MI)-MP2 method was considerably improved. 35 The root-mean-square errors (RMSE) of SCS-MP2 and SCS(MI)-MP2 toward the S66 data set amount to 0.87 and 0.38 kcal/mol, respectively.33 The SCS(MI)-MP2 calculations are less CPU-time intensive than the MP2.5 calculations and can thus be used for larger complexes.

E2 = −

nFσσ * eσ − eσ *

(3)

where n is the occupation number of the σ-orbital, F is the Fock matrix element between the σ and σ* orbitals, and eσ and eσ* are the respective orbital energies. According to eq 3, the energy of CT is directly connected with the matrix element F, which covers the overlap between the considered orbitals. Within the frontier molecular orbital theory (FMOT), 41 the CT interaction could be, to some extent, understood as the interaction between the HOMO of the electron donor and the LUMO of the electron acceptor. (To obtain a better understanding of these phenomena, a few more lower-lying occupied MOs can be considered. This is based on the idea that the MOs of graphene-like molecules are energetically close when the size of the molecule is increased.) Even though these two theories are independent, the charge transfer could be qualitatively described by a combination of both. 2.4. Free Energy Calculations. The association free energies were determined using two distinct methods: (i) the RR-HO-IG approximation in the case of aqueous solution combined with a continuum solvent model and (ii) the PMF 25330



Article

cubic box except for the periodic graphene...organic acceptor complex, which was done in a rectangular box. The size of the box varied with the size of the graphene-flake molecule. The box was prepared in a way that the distance of the solute from the box wall was at least 10 Å. The entire box containing the graphene flake, CT acceptor, and water molecules was minimized to avoid any possible close contacts and gradually heated to 300 K during a 20 ps long simulation, keeping the box volume constant. The initial velocities were chosen randomly from the Maxwell−Boltzmann distribution corresponding to the temperature of 10 K. Subsequently, the density was equilibrated by a 1 ns long simulation employing the biasing umbrella potential between the graphene flake and CT acceptor of a strength of 2.39 kcal/(mol Å2) with the equilibrium Z-distance 5.5 Å. A pressure of 1 bar and a temperature of 300 K were maintained by a Berensen barostat and thermostat, respectively.49 All of the 3 ns long umbrella windows started from the equilibrated coordinates and velocities, and the force constant of 2.39 kcal/(mol Å2) was used. Only the third nanosecond of each umbrella window was considered for the PMF evaluation while the first, second, and third nanoseconds were used for the block method error analysis. For the umbrella windows, Nosé−Hoover thermostat50,51 and Parrinello−Rahman barostat52 were employed to generate the isothermal−isobaric ensemble. The umbrella sampling simulations with periodic graphene complexes were performed in the same way. The x and y dimensions of the periodic box (i.e., the plane of the graphene) were 40.8 and 45.7 Å, respectively. The range of Deq values was extended, reaching the highest value of 17.0 Å, and the simulations were 4 ns long for each window. This provided suitably converged values with respect to all of the aspectsthe graphene size, the simulation length, and the extent of the dissociation. All of the explicit water simulations used the particle-mesh Ewald algorithm53 with a 10.0 Å direct-space cutoff for treating electrostatics and the plain 10.0 Å cutoff for the Lennard-Jones (LJ) interactions. The periodic graphene simulations in the gas phase were done with an identical setup, but the cutoffs were increased to 19.0 Å. For the gas-phase simulations of C2−C5 complexes, no cutoffs were used for LJ or for electrostatic interaction and no periodic boundary conditions were applied. The PMF simulations were performed in the Gromacs 4 program package.54

molecular dynamics simulation. The association free energies were calculated for the gas phase as well as for the aqueous environment. 2.4.1. RR-HO-IG Approximation. The complex and the respective monomers were optimized at a certain level of theory and the Hessian matrix was calculated providing a set of vibrational frequencies. The association free energies were calculated at 298.15 K and a pressure of 101 623 Pa on the basis of the equation ΔG = ΔE + ΔZPVE + ΔH(0 → T ) − T ΔS

(4)

where ΔE represents the change in the electronic energy, ΔZPVE is the change in the zero-point vibrational energy, ΔH is the change of enthalpy when passing from 0 to 298.15 K, and the last term in eq 4 stands for the change in the entropy at temperature T. The binding free energy in an aqueous environment ΔGw differs from the respective gas-phase values only by replacing the gas-phase interaction energy by interaction energy in the implicit solvation model (thus covering the solvation free energy change). The details of the construction of the free energy change can be found in e.g. ref 42. Various implicit solvation models were used depending on the level of the electronic energy calculations. For the BLYP-D3 calculations, the COSMO implicit solvation model was employed.43 For the DFTB-D calculations, the SMD model by Truhlar44 was used, whereas for the MM calculations, the generalized Born (GB) solvation model with parameters by Hawkins et al.45 was used. For the translational contributions to the association free energy, only the loss of one degree of freedom upon binging was assumed. This assumption is further discussed in section 3.3.2. 2.4.2. Potential of Mean Force. The PMF calculations rely on molecular dynamics simulations. Sufficient sampling of the configuration space is, however, computationally accessible at the empirical level of theory only. Hence, the complexes were described by GAFF, and the aqueous environment was modeled either by the GB implicit model or the TIP3P explicit water model.46 The bonded and nonbonded parameters were assigned by the Antechamber and Parmchk programs from the AMBER program suite.39 The binding free energies were calculated by the umbrella sampling technique,47 which introduces a biasing potential V of a quadratic form into the simulation (eq 5). V (d) = (1/2)k(d − Deq )2

(5)

3. RESULTS AND DISCUSSION 3.1. Interaction Energies. The interaction energies of the C2 complexes calculated at all levels of theory examined are summarized in Table 1. The MP2.5 values are claimed to be the

where d is the center-of-mass distance of the graphene flake and CT acceptor, k is a force constant, and Deq is the umbrella potential equilibrium distance. Here, only the distance in the Zdimension was biased, which yields the free energy changes of the association assuming the free movement of the CT acceptor molecule on the graphene-flake surface. Three atoms of the graphene flake were the subject of the positional restraints with a force constant of 2.39 kcal/(mol/Å2) to avoid undesirable movement of the flake during the simulation. We note that the CT acceptor was allowed to move freely in an xy-plane. All of the umbrella sampling simulations were evaluated by the weighted histogram analysis method (WHAM).48 Each of the 45 sampling windows varied Deq from 2.2 to11.0 Å in the solvent PMF while for gas-phase calculations 65 sampling windows with a range of Deq from 2.2 to15.0 Å were performed (i.e., 0.2 Å step window separation of both the solvent and gas phase). All the PMF calculations were done in a

Table 1. Interaction Energies of the C2 Complexes in (kcal/ mol)a methods

C2E

C2Q

E−Q diff

MP2.5/CBS/6-31G*(0.25) SCS-MP2/CBS SCS(MI)-MP2/CBS BLYP-D3 DFTB-D AMBER

−16.7 −17.1 −19.8 −14.5 −11.6 −13.4

−24.3 −25.0 −29.0 −20.6 −17.4 −22.1

7.5 7.8 9.2 6.1 5.8 8.7

a

There are more decimal places used for the E-Q difference calculation. 25331



Article

energy between the C5Ea and C5Qa complexes is 10.3, 9.8, and 10.7 kcal/mol for the BLYP-D3, DFTB-D, and AMBER methods, respectively. In terms of the benchmark data and the fact that the difference is well converged with respect to the graphene size, the value of 10.5 kcal/mol might represent a good estimate of the true difference in the stabilization of the TCNE and TCNQ complexes with graphene. The interaction energies calculated at the DFTB-D and AMBER levels converge slowly and are not too dependent on the orientation of the CT acceptor. However, the BLYP-D3 curves are not monotonic, and the energy difference in the two orientations is more pronounced (cf. Figure 2). 3.2. Charge Transfer. Table 3 shows the magnitude of charge transfer (in e) for all of the complexes. The NBO values

most accurate ones in the set and are used as the reference. The SCS-MP2 gives interaction energies which are within 1 kcal/ mol for both C2 complexes. Surprisingly, SCS(MI)-MP2 performs worse especially for the C2Q complex, where the complex stabilization is overestimated by about 4 kcal/mol. This might be attributed to the size of the complexes and their highly dispersion character, which were poorly supplied by the training S22 set55 used for the parametrization.34 Compared to the benchmark MP2.5 interaction energies, the BLYP-D3, DFTB-D, and molecular mechanics tend to underestimate the complex stabilization, although the AMBER molecular mechanics provide modest agreement of the absolute values and reasonable agreement in the relative difference between the two complexes. Table 2 shows the interaction energies of the C3, C4, and C5 complexes calculated at the BLYP-D3, DFTB-D, and MM

Table 3. Amount of Charge Transferred from C2, C3, C4, and C5 to the Acceptor Molecule

Table 2. Interaction Energies of the C3, C4, and C5 Complexes with TCNE and TCNQ (in kcal/mol)

a

complexa

BLYP-D3

DFTB-D

MM

C3Ea C3Eb C4Ea C4Eb C5Ea C5Eb C3Qa C3Qb C4Qa C4Qb C5Qa C5Qb

−20.1 −19.3 −19.7 −19.6 −22.9 −22.0 −26.9 −28.6 −31.2 −29.2 −32.2 −32.8

−13.8 −13.5 −14.2 −14.1 −14.5 −14.3 −21.9 −21.9 −23.8 −23.2 −24.3 −24.1

−15.7 −15.7 −16.5 −16.5 −17.0 −17.1 −24.6 −24.7 −26.6 −26.8 −27.7 −27.8

complex

NBO [e]

Mulliken [e]

complex

NBO [e]

Mulliken [e]

C2E C3Ea C3Eb C4Ea C4Eb C5Ea C5Eb

0.089 0.286 0.252 0.270 0.270 0.384 0.362

0.101 0.289 0.259 0.272 0.273 0.379 0.359

C2Q C3Qa C3Qb C4Qa C4Qb C5Qa C5Qb

0.046 0.212 0.292 0.337 0.258 0.379 0.423

0.071 0.217 0.299 0.337 0.272 0.375 0.426

Cf. Figure 1.

levels of theory. For each complex, two CT-acceptor orientations on the graphene surface (denoted as a and b) were investigated. The “a” orientation was identified as the more stable one in both the TCNE and TCNQ cases. The difference in the orientations is the least pronounced at the MM level. The stability of the TCNE complex is in agreement with the periodic DFT calculations performed earlier.10 The dependence of the interaction energy on the grapheneflake size is depicted in Figure 2. Following our expectations, the TCNQ complexes are more stable than the TCNE complexes. The stabilization of the complex is enormous, reaching about 20 and 30 kcal/mol in the TCNE and TCNQ cases, respectively. The average difference of the interaction

Figure 3. Dependence of the charge transfer on the size of the polycyclic aromatic molecule. The full lines stand for the “a” orientation while the dotted lines for the “b” orientation.

are plotted in Figure 3. It is shown that the Mulliken charge transfer agrees well with the (more accurate) NBO values and is generally comparable with the previous results.10,14 The charge transfer is comparable within the two CT acceptors studied with the fact that in the TCNQ case it is more orientation-dependent. It is known that the DFT functional chosen tends to slightly overestimate charge transfer (BLYP does not contain any exact exchange contribution);56 the actual absolute value is expected to be lower. The correlation coefficient R2 between the interaction energies and charge transfer calculated at the BLYP-D3 level is 1.00 for TCNE complexes and 0.97 for TCNQ complexes. The nonmonotonic behavior with increasing size is observed for charge-transfer as well as interaction energies. There may be various reasons for this. The frontier orbitals of the electron donors and electron acceptors considered in this paper are shown in Figure 4. The HOMO of C3 differs from that of C2 and C4, and this difference is quite large especially at the central benzene ring. Most likely, the symmetry of the highest occupied orbitals is quite sensitive to the molecular geometry. The larger the graphene model, the closer the orbitals. We note that the geometries were energy-minimized, resulting in a nonplanar shape of the polycyclic molecules, which may affect the orbital order. In conclusion, the different shape of the

Figure 2. Dependence of the interaction energy on the size of the polycyclic aromatic molecule. The full lines stand for the “a” orientation while the dotted lines for the “b” orientation. 25332



Article

which corresponds to the differences in the interaction energies. The gas-phase association free energies converge smoothly with the size of the complex. The average energies for C5E associated with both orientations were −3.1 and −8.2 kcal/mol for DFTB-D and AMBER, respectively, whereas for C5Q, they were −11.1 and −13.7 kcal/mol. The average difference between the interaction energy and association free energy for the DFTB-D and AMBER methods amounts to 3.2 and 2.2 kcal/mol, respectively. This difference is rather small, which indicates that the binding of the CT complexes studied is governed by the interaction energy, with the other terms (such as entropy) not being dominant or compensating for each other. When two molecules bind, they lose three degrees of freedom, and this loss is captured by the entropy decrease. However, when a molecule is adsorbed onto a surface (e.g., graphene), the number of degrees of freedom which are lost is reduced. In our calculations, only one degree of freedom was assumed to be lost upon binding. In other words, we supposed that the CT acceptor is able to move freely on the surface. Clearly, the larger the graphene model is, the better justified this assumption might be. The mobility is discussed below in section 3.3.2. The free energy of association in water has been calculated employing an implicit solvation model; in all instances, it has led to lower stabilization as compared to the gas phase. The MM binding free energies in a water environment are comparable with DFTB-D ones. These association free energies again converge smoothly with the size of the complexes, and the values for C5E averaged over both orientations were −10.1 and −12.1 kcal/mol for DFTB-D and AMBER, respectively, while for C5Q they were −18.2 and −17.3 kcal/mol. The association free energies for the C5 complexes differ from the gas-phase values by about 2 kcal/mol in the case of DFTB-D (SMD solvent). In the case of AMBER (GB solvent), the differences are larger and the solvent destabilizes the binding by about 5 kcal/mol for both C5 complexes. 3.3.2. Potential of Mean Force. Representative PMF profiles for C3 complexes are depicted in Figure 5. The PMF values were generated with the AMBER empirical potential. The

Figure 4. Frontier orbitals of the molecules investigated. Two energetically degenerated HOMO orbitals are shown for the C2, C3, and C4 molecules.

HOMO orbital is responsible for the unexpected shape of the curves as well as for the charge-transfer (nonmonotonous) values. 3.3. Free Energy Calculations. 3.3.1. RR-HO-IG Approximation. The binding free energies for C2−C5 with TCNE and TCNQ in the gas-phase and aqueous environment evaluated on the basis of molecular characteristics determined with the DFTB-D and AMBER methods are presented in Table 4. The gas-phase results will be discussed first. It is to be noted that the vibration analyses performed with both methods provided positive vibration frequencies for all of the complexes investigated. The only exception was the C2E complex, for which one imaginary frequency (−5.0 cm−1) was detected and discarded from further calculations. The association free energies determined at the DFTB-D level are less negative than those calculated at the MM level, Table 4. Association Free Energies (in kcal/mol) RR-HO-IG/DFTB-D

a

RR-HO-IG/MM

complex

gas

sol

gas

C2E C3Ea C3Eb C4Ea C4Eb C5Ea C5Eb periodicE C2Q C3Qa C3Qb C4Qa C4Qb C5Qa C5Qb periodicQ

−8.3 −9.6 −11.5 −10.5 −11.0 −10.7 −11.9

−6.2 −8.3 −10.2 −9.3 −9.7 −9.5 −10.7

−9.5 −12.6 −12.6 −13.4 −15.2 −16.8 −16.0

−12.6 −16.5 −17.3 −18.3 −19.2 −18.9 −20.2

−9.4 −14.2 −14.6 −16.7 −18.0 −17.4 −19.0

−15.6 −19.6 −19.6 −22.5 −23.3 −22.1 −22.2

PMF/MM sol

a

−3.6 −8.6 −8.6 −9.5 −10.3 −12.6 −11.7 −10.3 −14.1 −14.1 −17.4 −18.4 −17.3 −17.3

gas −10.5 −11.5 −11.8 −12.2 −12.5 −12.7 −12.7 −12.3 −16.8 −18.9 −19.0 −20.9 −20.8 −21.0 −21.6 −21.1

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

sol 0.14 0.28 0.11 0.40 0.16 0.28 0.26 0.98 0.50 0.29 0.23 0.17 0.29 0.27 0.14 0.77

−2.7 −4.6 −4.0 −5.2 −5.2 −4.9 −5.0 −7.1 −7.1 −10.2 −9.8 −11.4 −11.7 −12.4 −11.4 −13.7

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

0.16 0.25 0.50 0.39 0.17 0.10 0.05 0.15 0.14 0.06 0.37 0.20 0.69 0.14 0.37 0.69

One imaginary frequency of 5.0 cm−1 discarded. 25333



Article

and PMF values here lies in the solvent rather than in the different description of the translation. As shown in Figure 6 on the right, the mobility of TCNE on the C5 surface is quite high, and during the simulation time almost all of the surface is explored by the molecule. Thus, considering the loss of one degree of freedom upon binding in the RR-HO-IG approach, it is fully justified. One has to bear in mind that the MM approach allows the inclusion of CT effects only implicitly via Lennard-Jones parameters. An interesting feature of the graphene system is that it has very low absolute atomic charges on the surface-ring atoms when compared with the very small model of graphene systems like benzene. As we increase the size of the polycyclic molecule of the graphene flake, the charges of the carbon atoms are close to zero, from which we can conclude that for a graphene system with an infinite number of carbon atoms the average charge would go to zero. Regarding the benchmark C2E interaction energies calculations, the relative difference between TCNE and TCNQ PMF association free energies is expected to be more accurate than the absolute values, although still sufficient for insight into the thermodynamics. The PMF values are well converged with respect to the electron-donor size (Figure 7). The difference between the C5

Figure 5. Representative potential of mean-force curves.

average contact minima of the complexes were around 3.1 Å for the gas phase complexes and 3.3 Å for the solvent complexes. The TCNQ equilibrium distance was lower by about 0.2 Å when compared to TCNE. There were no solvent-separated minima found in the aqueous-environment PMF. Table 4 presents the determined gas-phase and water-association free energies. The errors were estimated by block averaging. Only one-dimensional biasing potential was introduced to allow the molecule to move naturally on the graphene surface. Contrary to RR-HO-IG, no assumptions about the mobility are needed, and consequently the PMF results should be considered as more reliable than the RR-HO-IG ones. The gas-phase results will be discussed first. All of the gas-phase PMF values are less negative than MM interaction energies by about 5 kcal/mol on average. The values are well converged and the association free energies with periodic graphene amount to −12.3 and −21.1 kcal/mol for TCNE and TCNQ, respectively. The lower absolute values of association free energies (in comparison with RR-HO-IG values) are caused by the loss of translational freedom. Figure 6

Figure 7. Dependence of the PMF association free energy on the size of the polycyclic aromatic molecule.

and periodic models in the gas phase is lower than 1 kcal/mol. In the case of water simulations, the difference is slightly higher, caused most likely by the higher mobility of the CT acceptor.

4. CONCLUSIONS By means of computational-chemistry methods, we have investigated two sets of charge-transfer complexes: TCNE and TCNQ adsorbed on polycyclic aromatic molecules of increasing size resembling graphene. The benchmark interaction energies in the gas phase were evaluated at the MP2.5/ CBS/6-31G*(0.25) level for coronene C2 complexes. SCSMP2 agreed better with the benchmark energies than SCS(MI)-MP2, which was attributed to the size of the complexes and/or the unreliable set used for the SCS(MI)MP2 parametrization. BLYP-D3, and surprisingly also molecular mechanical AMBER were shown to provide good agreement in the difference of the TCNE and TCNQ interaction energies. All of the non-MP2 methods provided slightly lower stabilization as compared to the benchmark MP2.5 values, whereas the DFTB-D performed worst among the methods in both the relative and absolute values of interaction energies. For the larger complexes, the interaction energies based on BLYP-D3 become more negative (i.e., more stabilizing) with the increasing polycyclic molecule size with the exception of

Figure 6. TCNE electron-acceptor mobility on the C5 surface. The positions of TCNE visited during 2 ns (gas phase) and 6 ns (solvent) long molecular dynamics simulations are shown.

shows the mobility of TCNE on the C5 molecule. Clearly, the molecule is not fixed on the surface, but it is not free either. The figure suggests that the actual number of the degrees of freedom lost upon the adsorption is not one, as was assumed in the RR-HO-IG calculations, but is higher. This gives an explanation to the difference between e.g. the C5E RR-HO-IG and PMF values, where the PMF is less negative according to the correct treatment of the TCNE-surface mobility. It should be noted that this emphasizes how valuable the data extracted from the molecular dynamics simulations might be, especially in the cases which are difficult to be tackled by the approximative RR-HO-IG method. The presence of a solvent lowers the stabilization substantially, and the periodic graphene PMF association free energies for TCNE and TCNQ complexes are −7.1 and −13.7 kcal/mol, respectively. The difference between the RR-HO-IG 25334



Article

(9) Chen, W.; Chen, S.; Qi, D. C.; Gao, X. Y.; Wee, A. T. S. J. Am. Chem. Soc. 2007, 129, 10418−10422. (10) Lu, Y. H.; Chen, W.; Feng, Y. P.; He, P. M. J. Phys. Chem. Lett. 2009, 113, 2−5. (11) Coletti, C.; Riedel, C.; Lee, D. S.; Krauss, B.; Patthey, L.; von Klitzing, K.; Smet, J. H.; Starke, U. Phys. Rev. B 2010, 81, 235401. (12) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887−1930. (13) Stone, J. A. Chem. Phys. Lett. 1993, 6, 101−109. (14) Zhang, Y. H.; Zhou, K. G.; Xie, K. F.; Zeng, J.; Zhang, L. H.; Peng, Y. Nanotechnology 2010, 21, 065201. (15) Munusamy, E.; Sedlák, R.; Hobza, P. ChemPhysChem 2011, 12, 3253−3261. (16) Granatier, J.; Lazar, P.; Otyepka, M.; Hobza, P. J. Chem. Theory Comput. 2011, 7, 3743−3755. (17) Sun, G.; Kurti, J.; Rajczy, P.; Kertesz, M.; Hafner, J.; Kresse, G. J. Mol. Struct. 2003, 624, 37−45. (18) Anota, E. C.; Juarez, A. R.; Castro, M.; Cocoletzi, H. H. J. Mol. Model. 2012, 18, 2147−2152. (19) Riley, K. E.; Pitoňaḱ , M.; Jurečka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023−5063. (20) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, 3742−3746. (21) Czaplewski, C.; Rodziewicz-Motowidło, S.; Liwo, A.; Ripoll, D .R.; Wawak, R. J.; Scheraga, H. A. Protein Sci. 2000, 9, 1235−1245. (22) Czaplewski, C.; Kalinowski, S.; Liwo, A.; Scheraga, H. A. Mol. Phys. 2005, 103, 3153−3167. (23) Podolsky, B. Phys. Rev. 1928, 32, 812−816. (24) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553−566. (25) Šponer, J.; Zgarbová, M.; Jurečka, P.; Riley, K. E.; Šponer, J. E.; Hobza, P. J. Chem. Theory Comput. 2009, 5, 1166−1179. (26) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (27) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114, 5149−5155. (28) Elstner, M.; Poreza, G. D.; Jungnickel, G.; Elstner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260− 7268. (29) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (30) Wang, J.; Wang, W.; Kolmann, P. A.; Case, D. A. J. Mol. Graphics Modell. 2006, 25, 247−260. (31) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (32) Kroon-Batenburg, L. M. J.; van Duijneveldt, F. B. J. Mol. Struct. 1985, 121, 185. (33) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. J. Chem. Theory Comput. 2011, 7 (8), 2427−2438. (34) Grimme, S. J. Chem. Phys. 2003, 118, 9095−9102. (35) Steele, R. P.; DiStasio, R. A.; Gordon, M. H. J. Chem. Theory Comput. 2009, 5, 1560−1572. (36) Janowski, T.; Ford, A. R.; Pulay, P. Mol. Phys. 2010, 108, 249− 257. (37) Ahlrichs, A.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Chem. Phys. Lett. 1989, 162, 165−189. (38) 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.1; Gaussian, Inc.: Wallingford, CT, 2009. (39) Case, D. A.; Darden, T. A.; Cheatham, III, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W. et al. AMBER 10; University of California: San Francisco, 2008. (40) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. Rev. 1988, 88, 899−926. (41) Fukui, K.; Yonezawa, T.; Shingu, H. J. Chem. Phys. 1952, 20, 722−725. (42) Cramer, C. J.; Truhlar, D. G. Rev. Comput. Chem. 1995, 6, 1−72. (43) Klamt, A.; Schüürmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 799−805.

C4E complex. This was attributed the unusual orbital shape caused by nonplanar optimized geometry. The simpler DFTBD and AMBER methods provided monotonically increasing stabilization with graphene model size. It was shown that the BLYP-D3 interaction energies highly correlate with the extent of charge transfer between the monomers calculated at the same level of theory. Mulliken charges were found to be very similar to NBO charges, reaching a value of about 0.4 e in the C5 complexes. The association free energies in the gas and aqueous phases were evaluated at the DFTB-D and AMBER levels of theory. The role of translational entropy was highlighted, just like the role of the solvent. Both lowered the stabilization of the complex as compared to the gas phase. An important finding was the fact that both CT acceptors are allowed to move on the graphene surface. In the gas phase, the movement is more restricted than in the solvent case, where the diffusive effects play a role. The RR-HO-IG approximation was shown to yield reasonable values only when a proper description of the translational entropy was used. A convergence behavior was observed for the interaction energies as well as for the association free energies. In the latter one, the C5 model of graphene provided sufficiently converged results with respect to those of the periodic model. Thus, the cluster simulations might provide authentic results in those cases where the periodic model is not accessible.

■

AUTHOR INFORMATION

Corresponding Author

*Tel +420 220 410311; e-mail [email protected]. Author Contributions ∥

These authors have contributed equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic [RVO: 61388963], and the Czech Science Foundation [P208/ 12/G016]. This work was also supported by the Operational Program Research and Development for Innovations−European Science Fund (CZ.1.05/2.1.00/03.0058). The support of Praemium Academiae of the Academy of Sciences of the Czech Republic awarded to P.H. in 2007 is also acknowledged.

■

REFERENCES

(1) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Science 2004, 306, 666−669. (2) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Katsnelson, M. I.; Grigorieva, I. V.; Dubonos, S. V.; Firsov, A. A. Nature 2005, 438, 197−200. (3) Tian, X.; Xu, j.; Wang, X . J. Phys. Chem. B 2010, 114, 11377− 11381. (4) Lu, C. H.; Yang, H. H.; Zhu, C. L.; Chen, X.; Chen, G. N. Angew. Chem., Int. Ed. 2009, 48, 4785−4787. (5) Shan, C. S.; Yang, H. F.; Song, J. F.; Han, D. X.; Ivaska, A.; Niu, L. Anal. Chem. 2009, 81, 2378−2382. (6) Baby, T. T.; Aravind, S. S. J.; Arockiadoss, T.; Rakhi, R. B.; Ramaprabhu, S. Sens. Actuators, B 2010, 145, 71−77. (7) Zhang, Z.; Huang, H.; Yang, X.; Zang, L. J. Phys. Chem. Lett. 2011, 2, 2897−2905. (8) Ishikawa, R.; Bando, M.; Morimoto, Y.; Sandhu, A. Nanoscale Res. Lett. 2011, 6, 111. 25335



Article

(44) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378−6396. (45) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 19824−19839. (46) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (47) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187−199. (48) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011−1021. (49) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (50) Nosé, S. Mol. Phys. 1984, 52, 255−268. (51) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (52) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182−7190. (53) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (54) Hess, B.; Kutzner, C.; Van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (55) Jurečka, P.; Šponer, J.; Č erný, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985−1993. (56) Steinmann, S. N.; Piemontesi, C.; Delachat, A.; Corminboeuf, C. J. Chem. Theory Comput. 2012, 9, 1626−1640.

25336


Adsorption of Organic Electron Acceptors on Graphene-like Molecules

Recommend Documents