Structure and Properties of Double-Sandwich Complexes at Graphene

Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

C: Physical Processes in Nanomaterials and Nanostructures

Structure and Properties of Double-Sandwich Complexes at Graphene Surface: A Theoretical Study Vijay Madhav Miriyala, Rabindranath Lo, Susanta Haldar, Amrit Sarmah, and Pavel Hobza J. Phys. Chem. C, Just Accepted Manuscript • Publication Date (Web): 15 May 2019 Downloaded from http://pubs.acs.org on May 15, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Structure and Properties of Double-Sandwich Complexes at Graphene Surface: A Theoretical Study Vijay Madhav Miriyala,1, 3 Rabindranath Lo,1, 3 Susanta Haldar,2 Amrit Sarmah1, 3 and Pavel Hobza*, 1, 3 AUTHOR ADDRESS 1

Department of Computational Chemistry, Institute of Organic Chemistry and Biochemistry, Czech Academy of

Sciences, Flemingovo Náměstí 542/2, 16610 Prague, Czech Republic. 2

School of Chemistry, University of Bristol, Cantock’s Close, Clifton, Bristol, BS8 1TS, United Kingdom.

3

Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacký University, Šlechtitelů 27, 78371 Olomouc, Czech Republic. *E-mail: [email protected]

ABSTRACT

Graphene and its derivatives are useful building blocks for the bottom-up assembly of advanced functional materials. Non-covalently functionalized graphene networks offer a wide range of applications. We investigated the formation of sandwich-like three layered nanostructures with graphene. Novel architectures have been generated by stacking selected suitable organic molecules, based on their characterized donor and acceptor strengths, vertically on the graphene surface. This paper describes the adsorption of electron-acceptor and electron-donor molecules on the graphene layer through non-covalent interactions. Cluster as well as crystal models of the

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 41

graphene surface have been selected to design sandwich-like two layered and three-layered structures, and their stabilities have been verified using DFT calculations. Further, stability of the complexes has been confirmed on the basis of factors such as interaction energy and charge transfer. The stability of the macrostructures has been tested by metadynamics simulations. We have found that the most stable complex C4...HAT-CN...TAB prefers the double sandwich state over the dissociated state.

Introduction Ever since it was made,1 graphene (G), due to its special honeycomb-like structure, has been known to possesses some very special properties such as the intersection of the valence band and the conduction band at the Dirac point, where the energy-momentum dispersion is linear. This makes graphene a very good charge carrier.2 The interactions between the graphene and various chemical species such as organic biomolecules and electro-active molecules, where graphene acts as an electron-donor and the latter systems as electron-acceptor, have had many potential applications in the field of graphene-based devices.3–11 Therefore, both experimental and theoretical investigations have been performed intensively to assess the effect of the interaction of graphene with different molecules on the electronic properties of graphene. To this end, a large number of adsorption studies have been carried out so far on the graphene surface, where the adsorbate particles ranged from small organic molecules to biomacromolecules such as proteins, amino acids, ss-RNA, ds-DNA, inorganic molecules and also monovalent and divalent ions.12–16 Many of those adsorption studies mainly focused on the adsorption pattern, stability and the origin of the stability.17,18 In most of the cases, it was confirmed that the molecules bind via non-covalent interactions, mainly via stacking [For more details and information please read Ref. 16 and the references within the article]. The adsorption of a conjugated electro-active donor or acceptor on the graphene surface has been a promising


2

Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


way for the design of new functionalized nanostructures.19–24 The work done by Chen et al. on the charge carrier mobility in a noncovalently functionalized graphene has revealed new prospects in the design of conductive graphene electrodes with a tailored work function.25 Among the various molecules studied for adsorption on the graphene surface, conjugated organic donors and acceptors are very important. These are used to functionalize graphene because of their compatibility and binding ability via π-π stacking, electrostatic force, and hydrogen bonds. In the past, several studies dealt with the adsorption of such small molecules as CO, NO2, NH3 etc. on the pristine graphene and nanotubes.16 With recent advances in science and increased computational power, it has been reported that several organic electron-acceptors, donors, and even DNA fragments could play a significant role in changing the electronic structure of the system, including tailoring the electronic properties of graphene via adsorption.26-28 Graphene is important in the field of nanoengineering because of the presence of highly mobile πelectrons, which is widely used in semiconductor-based nano-electronics. Graphene is thus an electron-donor. It has been confirmed that an electron-acceptor molecule can be adsorbed at the graphene surface and besides dispersion energy, the complex is stabilized by charge-transfer energy. Upon the adsorption of an electron-acceptor molecule, an increase in its electronic properties has been detected.21 The present study mainly focuses on the construction of graphene-based sandwich and doublesandwich (DS) complexes. In the former case, an electron-acceptor (EA) molecule is sandwiched at the graphene surface, while in the latter case another molecule, an electron-donor (ED), is attached at the sandwich complex facing an electron-acceptor system. The sandwich complex (G...EA) is formed by depositing an electron-acceptor at the graphene surface. Additional deposition of an electron-donor leads to the formation of either sandwich complexes between graphene and an electron-donor (G...ED) or DS complexes (G...EA...ED), where an electrondonor


3


Page 4 of 41

faces the electron-acceptor part of the G...EA sandwich. The formation of the G...ED complex is not favored because both subsystems, G and ED, are electron-donors and the only stabilization originates from the omnipresent dispersion energy. The G...EA...ED DS complex is preferentially formed if the binding free energy of the ED to the G...EA complex is higher than that of the ED to G.

Strategy The adsorption of several organic molecules on the graphene surface will be investigated by using either the cluster or crystal models. Graphene in the former model is represented by circumcircumcoronene (referred to as C4) while in the latter model by an infinite graphene surface. The selection of such a finite-size-cluster model as a graphene sheet depends on at least two complementary factors. First, the model should be large enough to avoid edge–edge interactions. Second, the model should not be too large to restrict the use of large/flexible basis sets. The C4 system that has been widely used fulfills both requirements.29 In both models, the electron-acceptor molecule is placed horizontally with respect to the center of the graphene. Finally, the electrondonor molecules are placed on top of these complexes. An infinite graphene layer is used to perform all the periodic calculations employing the periodic boundary conditions in two dimensions. The stability of different sandwich complexes will be estimated on the basis of interaction energies calculated by both cluster and crystal models at the DFT level. DFT methods augmented by dispersion energy (DFT-D) are known to provide reliable estimates of interaction energies for different types of noncovalent complexes including the presently investigated complexes with dominant dispersion and charge-transfer energies. The present paper will preferentially use the cluster model, whose comparison with the physically more justified crystal model, albeit with


4

Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


some limitations, makes it possible to estimate its accuracy. The preferential binding of different EDs at the G...EA complexes cannot be based on interaction energy – binding free energy should be utilized instead. This study has used metadynamics simulations (MetaD) along with the standard molecular mechanics (MM) based molecular dynamics (MD) simulations. Predictions made by MetaD thus critically depend on the quality of MM interaction energies, which can be evaluated by comparing MM and DFT-D interaction energies. Finally, the following computation strategy will be used: i) the properties of isolated subsystems (graphene, EA, ED) will be investigated; ii) the structure and properties of G...EA complexes will be studied at cluster and crystal levels; iii) the structure and properties of G...EA...ED complexes will be evaluated at cluster and crystal levels as well as using MetaD simulations; the performance of MM used in MetaD will first be assessed by its comparison with DFT-D results.

Systems investigated Graphene. In the cluster model, graphene is represented by circumcircumcoronene (C4), having 96 carbon and 24 hydrogen atoms (Figure 1). In the crystal model of graphene, the Brillouin zone is mapped with in-plane k-point densities corresponding to an 8x8 mesh in a primitive cell. Electron-acceptors. The following electron-acceptors, shown in Figure 1(b–e), have been investigated: HAT-CN (1,4,5,8,9,11-hexaazatriphenylenehexacarbonitrile),30 which is a planar molecule with a triangular π-deficient area formed due to the presence of three pyrazine rings attached to the central benzene ring and six -CN substituents that lie around the periphery of the molecule. We have also chosen hexacyanobenzene (HCB), 1,3,5-triazine-2,4,6-tricarbonitrile (TACN) and 2,4,6-trinitro-1,3,5-triazine (TANO), which are also classified as strong electronacceptors. The main criterion that makes these monomers suitable candidates for electron


5


Page 6 of 41

acceptors is due to the presence of deactivating substituents on the aromatic rings that induces withdrawal of electrons making the aromatic rings electron deficient.

electron-donors. To complete the formation of DS complexes, the following well-known small electron-donors, shown in Figure 1(f–i), with different donating abilities have been considered: 1,3,5-trimethoxybenzene (TMOB), hexamethylbenzene (HMB), 1,3,5-triaminobenzene (TAB) and tetrathiafulvalene (TTF). The first three monomers contain activating substituents on the aromatic ring. Electron donation occurs here due to resonance effect, as the atoms connected to the aromatic ring N and O possess non-bonding valence shell electron pairs. Therefore, electrons flow into the aromatic ring by p-π conjugation. Methyl substituents in HMB also increase the nucleophilicity of the aromatic ring in the same fashion. TTF is a unique imidazole-derived neutral organic super-electron-donor that falls under the category of doubly-bridged donors (DBDs).31 It is a neutral stable compound with four sulfur atoms attached to the central double bond. The sulfur atoms are able to donate electrons to the π system, hence making TTF an electron-rich donor.

Computations Quantum Chemical Calculations The cluster model. The structures of all isolated subsystems as well as G...EA and G...EA...ED complexes have been optimized at the DFT-D3/B97D/TZVPP level of theory with the TURBOMOLE 6.6 suite of program.32 In recent studies, B97-D functional has been successfully utilized for the adsorption of aromatic molecules on graphene.33 Grimme’s advanced dispersioncorrected approach34 (DFT-D3) with the TZVPP basis set has been incorporated into the model.


6

Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The interaction energies of the complexes have been calculated using the following equation: E(interaction energy) = E(complex)-(E(molecule 1)+E(molecule 2))………(1), where, E(complex) stands for the total energy of the complex and E(molecule 1) and E(molecule 2) represent the energies of all isolated subsystems. The basis set superposition error is not considered here. NBO analysis has been performed using the NBO 3.1 program35 as implemented in Gaussian 09,36 employing optimized complex geometries at the B97D/def2-SVP level. Electrostatic potentials have been generated on 0.001a.u. molecular surfaces of selected systems at the B97D/ def2-SVP level. Each subsystem is characterized by its ionization potential (IP) and electron affinity (EA) calculated as the difference between the energy of a neutral subsystem and the respective cation or anion radical, where the structure of each subsystem is fully optimized. IP = the energy of the cationic state − the energy of the neutral state ………..(2) EA = the energy of the neutral state − the energy of the anionic state…………(3)

The DFT-D3 calculations have been automated using the Cuby framework.37 Constrained DFT calculations performed with NWChem38 have been used for the proper estimation of the charge transfer in (C4...EA)...ED DS. The constraint in this case is to force a certain region or a fragment of the complex to carry a predefined charge by introducing an auxiliary potential. In constrained DFT the auxiliary potential implements a charge constraint where we can localize a charge on part of the complex i.e. one of the monomer fragments in our case. By doing so we tell the program implicitly that the partial charges on the particular fragment we constraint


7


Page 8 of 41

add up to zero. Therefore, the program looks for the lowest energy state that satisfies the constraint we have imposed.39-42 Hence, two calculations have been performed for each complex. The first one is a standard DFT calculation, while the other is a constrained DFT where the charge on the (C4...EA) is constrained to zero, hence preventing any charge transfer between (C4...EA) and ED. Both calculations have been performed at the B97D/6-31G* level of theory. Dispersion has not been included in this case. Interaction energies for DFT and cDFT as well as the difference between them have been computed. Based on the difference in interaction energies between DFT and cDFT (a non-zero value implies the occurrence of charge transfer), one can predict the charge transfer between the monomers of the complex. Lowdin population analysis has also been performed to quantify the charge transfer. Bader charges have been calculated using a fast algorithm developed by Henkelman group at the University of Texas at Austin.43 The program reads charge densities in the Gaussian CUBE format. The program outputs the total charge associated with each atom. The crystal model based on periodic boundary calculations. All DFT-based ab initio calculations have been performed using the open-source Quantum Espresso 5.4.0 code.44 Troullier and Martins45 norm-conserving pseudopotentials have been used along with exchangecorrelation (XC) effects incorporated at the Becke, Lee, Yang, Parr (BLYP)46 level. The van der Waals correction to the calculations has been incorporated through semi-empirical Grimme’s DFT-D2 methodology47 as implemented in Quantum Espresso. This approach has been successfully implemented to describe the graphene-based structures.48 In the previous reports, this method performed reasonably well to describe the electronic structures of graphene nanoribbon/h-BN49-51 and graphene nanoribbon/AlN52,53 interfaces with the ambient accuracy. The standard kinetic energy cutoff for wavefunctions and the corresponding cutoff for charge density and potential are set at 40 Ry and 160Ry, respectively. The reciprocal space integrations have been performed at the “Tao” point. The model systems are considered inside a simple cubic super-cell with the 8×8


8

Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


graphene layer with the periodic images separated by at least 10 Å of vacuum space in both x and y directions in order to maintain negligible interactions between the repetitive units. The forces have been minimized up to < 10-3 Ry/a.u. for all the atoms to incorporate the fully relaxed atomic geometries. Molecular Dynamics Simulation. We have carried out a standard molecular dynamics (MD) simulation using the molecular mechanics force field to check the stability of the different complexes formed at room temperature. For the MD simulation, we prepared a periodic graphene with the dimensions of 38.0 x 38.0 Å in the x and y directions and later built sandwich complexes in the direction of the z axis, like in the cluster model. In the next step, we created a rectangular box by maintaining the same dimensions in the x and y directions while extending the z direction to approximately 40.0 Å to avoid interactions between the graphene and its periodic image. General AMBER force-field (GAFF)54,55 parameters were used for the MD simulation of all the complexes. The charges were calculated using the restrained electrostatic potential (RESP)56 fitting procedure. The RESP fit was performed onto a grid of electrostatic potential points calculated at the HF/6-31G* level. The MD simulation was performed with and without explicit solvent water molecules. For the solvent-phase calculations, complexes were solvated with 810 explicit TIP3P57 water molecules. Herein, the stability of all the complexes was investigated using a well-tempered version of the metadynamics (WTMetaD) simulation,58 which was run using the same simulation protocol as adopted in the standard MD simulation. Metadynamics (MetaD) simulation is a well-established enhanced sampling technique to overcome free-energy barriers; it provides a more robust exploration of free-energy surfaces (FESs). This technique accelerates the sampling through the introduction of a history-dependent biasing potential at selected degrees of freedom, called collective variables (CVs). The last frame found in the gas-phase MD equilibration was taken as


9


Page 10 of 41

the starting structure for the WTMetaD simulation. The plumed 2.2.3 plugin59 was used to carry out the simulation with the Gromacs-5.1.2 code.60 We have also applied a recently developed reweighting procedure61 to reconstruct a different FES. This algorithm allowed us to rebuild the Boltzmann distribution relative to CVs different from those biased in the actual WTMetaD simulation. We did not run a new WTMetaD simulation to obtain the reweighted FES.

Results and Discussion

I. The cluster model

Isolated subsystems. All the subsystems, i.e. electron-acceptors and electron-donors, are characterized by molecular electrostatic potential (ESP), its maximum value Vs,max (in kcal/mol), ionization potential and electron affinity (cf. Table 1). The characteristics of electron-acceptors are discussed first. The molecular electrostatic potentials of selected electron-acceptor molecules are shown in Figure 2, where the magnitude of the π-hole, Vs.max, of the acceptors is given below their ESPs. In HAT-CN, the π-hole is delocalized over the three six-membered rings attached to the central six-membered ring of the molecule, whereas in HCB, TACN and TANO, the π-hole is concentrated at the center of the only six-membered ring. Investigating the entries below the ESPs in Figure 2, we found the highest Vs,max for TANO followed by HAT-CN. However, upon further investigation, the electron affinities of the EAs (cf. Table 1) were found to be in the reverse order, with the highest electron-affinity value for HAT-CN. From this, it can be concluded that these two EAs, HAT-CN and TANO, will be the most efficient candidates. This statement is also supported by the fact that, due to the full or partial positive charge on the nitrogen atoms directly attached to the central aromatic ring, they all have a moderate to strong electron-withdrawing inductive effect,


10

Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-I. Thus, making the aromatic ring very electron poor in turn making them strong electron acceptors. Electron-donor ability is estimated on the basis of the magnitude of ionization potential. The lower the ionization potential value, the better the electron-donating ability of the system. In this study, we have used four donor systems. Investigating entries in Table 1, we found that TTF and TAB are the best electron-donors. While TTF contains four sulfur atoms that possess strong donor ability, TAB has three -NH2 substituents that activate the aromatic ring with their resonance donating effect. groups A plot of Ionization Potentials and electron Affinities of the all the subsystems have been presented in Figure 3.

Sandwich complexes. C4 (electron-donor) was combined with four electron-acceptors. The optimized structures of the resulting complexes are shown in Figure 4. Calculated characteristics (stabilization energy ΔE, dispersion energy, the distance between the centroids of theE, dispersion energy, the distance between the centroids of the subsystems, the charge transfer between the donor and acceptor from the NBO and HOMOLUMO gap) are collected in Table 2. The intermolecular distances, as measured between the centroids of the subsystems, were found to be highly uniform, lying between 3.3 and 3.5 Å. Evidently, the intermolecular distance is mainly determined by the van der Waals (vdW) repulsion, which is proportional to the sum of the vdW radii of the respective atoms. The highest stabilization energy of 46.5 kcal/mol was found for the C4…HAT-CN complex, which is in full agreement with the highest electron affinity found for HAT-CN. Investigation of the energy components has revealed that the complex is stabilized not only by charge-transfer energy (see NBO CT) but also by dispersion energy. The stabilization energy of C4...HCB is about ~16 kcal/ mol weaker, but this difference comes exclusively from dispersion energy, which is weaker by ~24 kcal/mol. Chargetransfer energy and the HOMO-LUMO gap for this complex are comparable to those of C4...HAT-


11


Page 12 of 41

CN. The stabilization energies of the remaining complexes are lower, with the main contribution to the stabilization coming from dispersion. With the highest ΔE, dispersion energy, the distance between the centroids of theE, dispersion energy and charge transfer, C4...HAT-CN and C4...HCB form stable sandwich complexes and are potential candidates for DS complexes. In DS complexes, EDs are attached to C4...EA sandwich complexes, specifically to their EA parts. Before the modelling of DS complexes, first we have analyzed the characteristics of the C4...ED sandwich (two-body) complexes. The results are summarized in Table 2, while their structures are shown in Figure 5. Evidently, the stabilization energies of all these C4…ED sandwich complexes are lower than the C4…EA complexes presented in Table 2, and the strongest complexes are formed between the subsystems with the highest donor and acceptor abilities. Double-sandwich complexes. In the previous section, sandwich complexes between C4 and four EAs were investigated. Complexes with HAT-CN and HCB were shown to be the strongest. This is due to the presence of very strong electron withdrawing -CN groups attached directly or in conjugation with the central aromatic ring. Double-sandwich (DS) complexes were thus constructed from (C4...HAT-CN) and (C4...HCB) combined with four electron-donors. The efficient formation of C4...EA...ED DS complexes is conditioned by the higher stabilization energy of the (C4...EA)...ED complex in comparison with the C4...ED one. Table 2 summarizes structural and energy characteristic of these complexes. Table 2 shows that the EA...ED intermolecular distances are very similar to those found for the C4...EA complexes, all being between 3.1 and 3.6 Å as measured between the centroids of the subsystems. The differences in the stabilization energies of the (C4...EA)...ED and C4...ED complexes are presented in the 7th column of Table 2 as ΔE, dispersion energy, the distance between the centroids of the E’. This value is decisive for the successful formation of the DS complex. The more negative the difference, the higher the probability of the formation of the DS


12

Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


complex. This difference was the most negative for C4...HAT-CN...TTF followed by C4...HATCN...TAB and C4...HCB...TAB. In TAB, density of the π-cloud over the benzene ring increases due to the presence of lone-pair containing -NH2 functional groups. In other words, the aromatic ring is activated by the increased electron density over the benzene ring through a resonance donating effect also known as the +R effect due to presence of substituent groups such as -NH2. However, when the alkyl substituents such as -CH3 are attached to the benzene ring, the electron density on the aromatic ring increases due to weaker inductive donating effect. TTF on the other hand possesses strong electron donating ability due to high-lying HOMO that arises from core sulfur atoms. The alignment of TTF parallel to an acceptor molecule also leads to charge transfer, which we show in the next section. Therefore, donors such as HMB are weaker compared to TAB and TTF. Hence, these results suggest that the EDs, TTF and TAB form stable DS complexes with C4…HAT-CN or C4…HCB units. Consequently, these EDs prefer to interact with EA in C4…EA unit compared to C4 surface. Optimized structures of various double sandwich complexes are presented in Figures 6 & 7. A plot of interaction energies in kcal/mol of stable sandwich and double sandwich complexes is presented in Figure 8. Charge-transfer analysis based on cDFT All the complexes investigated are supposed to be partially stabilized by charge-transfer interaction. To elucidate the role of charge-transfer in the formation of stable DS complexes, we have performed constrained DFT for all the (C4...HAT-CN)... DS complexes. We constrained the total charge over the (C4...HAT-CN)... system to zero. As the overall DS complex is neutral, the charge over the donor species also became zero. Next, we performed a standard DFT, where the flow of charge between the monomers is not restricted. The difference between the cDFT and DFT interaction energies, ΔE’’ (Table 3), indicates the charge-transfer contribution to the overall


13


Page 14 of 41

stabilization energy. ΔE’’ energies for (C4...HAT-CN)...TTF and (C4...HAT-CN)...TAB are 6.4 kcal/mol and 5.3 kcal/mol, respectively, showing that the charge transfer contribution for these complexes is relatively the largest. Hence, the cDFT results emphasize the role of charge-transfer in the formation of stable DS complexes. Now, we have analyzed the cDFT results for two less preferable DS complexes formed by two EDs, HMB and TMOB, respectively. Two DS complexes, (C4...HAT-CN)...HMB and (C4...HAT-CN)...TMOB, exhibit a very low ΔE’’ of 0.6 kcal/mol and 0.8 kcal/mol, respectively, which implies that the charge-transfer contribution in these DS complexes is considerably smaller. To support this argument further, we added the charge transfer calculated from the Lowdin charges obtained (from the outputs) for DFT calculations. The calculated values fully confirmed previous results. The largest relative charge transfer was again found for (C4...HAT-CN)...TAB and (C4...HAT-CN)...TTF complexes. To further testify the charge transfer, we have calculated the Bader Charges using the Henkelman algorithm over the Gaussian density cube. The results, presented in Table 3 follow similar trend and are comparable to that of the Lowdin charges discussed previously. As mentioned in the earlier section, in case of TAB there is an extra electron density over the aromatic ring due to presence of three -NH2 group that activate the electron density through a resonance donating effect leading to an effective charge transfer between the donor and the acceptor. In case of TTF the high lying HOMO, lower intermolecular distance that leads to a strong π-π overlap and the molecular fastening effect of the alkylthio group leads to an easier donation and charge transfer. The charge transfer is not prominent in case of HMB and TMOB as the -CH3 groups of HMB increases the electron density through a mere inductive donating effect, which does not induce any charge transfer between the donor and acceptor. The -OCH3 on the other hand, due to the presence of electronegative oxygens, pulls the electron cloud towards itself by exerting a negative inductive effect (-I effect). There is an overlap of the oxygen p-orbital with the extended π-electron cloud of the benzene ring, hence delocalizing the electrons into the aromatic ring. This in turn leads to not only poor


14

Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


donating ability but also hinders charge transfer. Whatever stabilization occurs in the later systems is only due to the omnipresent dispersion. Hence, it is confirmed that the complexes (C4...HATCN)...TTF and (C4...HAT-CN)...TAB are stabilized by charge transfer interactions along with the omnipresent dispersion. II. The crystal model Three DS complexes based on graphene (G) surface, specifically (G...HAT-CN)...TAB, (G...HCB)...TAB and (G...HAT-CN)...TTF, have been considered. The EDs in these complexes are strategically placed over the graphene(G)...EA, which forms the bottom layer and interacts directly with the ED. The thermodynamic stability of the model system is an effective measure to predict the structural stability of the proposed hypothetical DS systems. Based on the computed ∆E value of 43.8 kcal/mol, the DS formation is highly favorable for the (G...HAT-CN)...TAB complex. Due to the strong acceptor properties of HAT-CN and the strong donating ability of graphene, the DS complex has been stabilized by the charge-transfer interaction. The (G...HAT-CN)...TTF also found to be a good candidate for DS formation. The computed ∆E of 32.1 kcal/mol is also relatively high. It is worth mentioning here that the electron-donating ability of TTF is relatively lower than that of TAB. Although dispersion is the predominant interaction between graphene and HAT-CN, charge transfer is important for energy minimization in the interaction between G...EA and ED. The higher the amount of charge transfer, the greater the stability. A critical analysis of the stabilization process also indicates the possible existence of a dynamic charge-transfer equilibrium between the interacting systems. Investigation of the (G...HCB)...TAB DS complex has revealed that despite the fact that HCB is a good EA and TAB is a good ED, the resulting ΔE, dispersion energy, the distance between the centroids of the E is the least favorable. There are two possible


15


Page 16 of 41

explanations for this. Firstly, due to its smaller size, HCB does not have an effective dispersion interaction with the graphene layer (cf. Table 2). The second explanation is associated with the dynamic charge-transfer aspect. The graphene sheet has a significantly high amount of mobile πelectrons on the surface. While HCB has easily accumulated the electrons from the graphene surface, the lower level of electron delocalization brings the system to a saturation point, which limits its inherent electron-accepting ability. The placement of TAB above the G...HCB structure leads to a negligible charge transfer between G...EA and ED, and the overall DS formation becomes thermodynamically unfavorable. It can also be argued that a synchronized electronic charge transfer between the individual layers is a prerequisite for the DS formation. To shed some light on the chemical properties of the newly proposed DS complexes, it is essential to develop a comprehensive understanding of the electronic behavior of the system. The computed density of states (DOS) plots for the systems are presented in Figure 9, from which the different electronic structure of each system is evident. The characteristic strong charge-transfer interactions for (G...HAT-CN)...TAB and (G...HAT-CN)...TTF are attributed to the intense peak at the Fermi level. This also shows the induced metallic nature of the two systems. Similarly, the insulating nature with some distinct band gap for the system arises from the corresponding DOS spectrum of (G...HCB)...TAB. The characteristic electronic charge transfer between the different layers of the DS is not observed either. The presence of electronic states just above and below the Fermi level for (G...HAT-CN)...TAB and (G...HAT-CN)...TTF facilitates the charge transfer in the DS formation process. The contributions of different types of atoms to the overall DOS can be understood from the projected density of states (PDOS) plots reported in Figure 10. As we have observed in the DOS spectrum, the peaks appear near the Fermi level mostly originated from the contribution of carbon atoms. However, the contributions of nitrogen atoms to the electronicenergy states located just below and above the Fermi level are also distinct in (G...HAT-CN)...TAB and (G...HCB)...TAB. It also appears from the PDOS plot that the highest occupied state of the


16

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(G...HAT-CN)...TTF is mostly populated by the electronic contributions of the sulfur atom. A figure containing the relevant band structure of the pristine graphene along with the three DS structures is included in the supporting information section (Figure S3). The locations of the highest occupied states and the lowest unoccupied states in the orbital energy spectrum of a particular system provide useful insights into the charge-transfer possibilities. The isosurface densities for the highest occupied and lowest unoccupied states of the three systems are reported in Figure 11. In (G...HAT-CN)...TAB, HOMO is spread across the graphene surface and LUMO is mainly concentrated on HAT-CN. Here, the electronic charge transfer mainly occurs between the graphene and HAT-CN. It is worth mentioning here that, within a simple approximation one can estimate charge transfer behavior from the relative shapes/positions of the HOMO and LUMO. If the two span different regions of the molecule, then excitation of a HOMO electron into the LUMO would constitute a spatial transfer of charge. Apparently, it is not the most rigorous way to account the intermolecular charge transfer, but it provides qualitative information for minimal effort. The larger HOMO-LUMO gap as well as the lower possibilities of charge delocalization over HCB has substantially reduced the electronic charge transfer process in the system. Subsequently, (G...HAT-CN)...TTF also exhibits some characteristics similar to the previous complex. Taking analogy from the previous discussions it can be argued that a favorable positioning of the HOMO-LUMO lobes accounts a favorable charge transfer in this particular system. It is to be mentioned here that all important results from the cluster-model section concerning the stability of selected DS structures have been fully confirmed by physically justified crystal model calculations.

III. Free-Energy Calculations


17


Page 18 of 41

Molecular mechanics interaction energies As mentioned above, the ultimate decision on the stability of DS complexes should be made on the basis of ΔG and not ΔE values. The respective free-energy calculations are based on the molecular mechanics, and the method’s ability to describe the energetics of DS complex formation should be tested first. Both static and dynamic calculations on the DS complexes have been performed using the MM force field. By means of a metadynamics simulation, the barrier heights, going from the DS state to the dissociated state on the graphene surface, have been thoroughly computed with the estimate of their difference in the binding free energy. The order of stability followed by MM interaction energies is similar to that of DFT interaction energies. A comparison of the three complexes has shown that C4…HAT-CN…TAB is the most stable DS complex. The interaction energy calculated from MM is -26.5 kcal/mol, while the interaction energy calculated from DFT is -26.7 kcal/mol. The interaction energies of the dimer complexes such as C4…TAB and HAT-CN…TAB is calculated to be -16.5 and -25.6 kcal/mol, respectively. Both of the dimer complexes, C4…HAT-CN and HAT-CN…TAB, are lower in energy (by 10 and 0.9 kcal/mol, respectively) relative to the DS complex. These results confirm that the DS with HAT-CN as an EA and TAB as an ED on the cluster model C4 is stable as predicted by the DFT method. In order to shed light on the stability of DS complexes, we calculated the difference in the free energy of both the DS (associated) state and collapsed (dissociated) state as well as the barrier height between the two states. Table 4 shows the free-energy landscape of all three DS complexes with the estimates of their barrier heights. The MetaD simulation efficiently describes the DS (associated) state and collapsed (dissociated) state of all the complexes using the coordination number62,63 (cf. Eq. 3 and 4 in SI) as the chosen collective variables (CVs). The addition of biased potential during the MetaD simulation has prompted a transition state between the states where both EA and ED molecules form a stacked


18

Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


complex and a dissociated state where each of them is stacked on the graphene surface separately. While sampling a DS complex such as G…EA...ED requires a simulation in three-dimensional (3D) space, the sampling space for the dissociated complexes on the graphene surface has been reduced to two dimensions (2D). Sampling in 2D space is much simpler than in 3D space. The coordination numbers are very high – as many as 80 arbitrary units – when the DS is formed, whereas they are close to zero when the DS collapses and forms a dissociated state. Figure 12 depicts the free-energy surface (FES) of the DS complex G…HAT-CN…TAB after 200 ns of a gas-phase WTMetaD simulation at 300 K as a function of Sstack and Sdiss as the two collective variables (CVs). Careful observation of the FES shows two successive minima, state A (DS state) and state B (dissociated state). State A occurs within the Sstack and Sdiss coordination numbers of 40–70 and 5–10 arbitrary units, respectively. For state B, both Sstack and Sdiss coordination numbers decrease to 0–10 and 0–5 arbitrary units, respectively, which indicates that when the DS configuration collapses, the ED molecule apparently loses the stacking interaction with the EA and finally reaches the dissociated state. For the better understanding of the stacking interaction in both states, a reweighted FES is provided in the supporting information (cf. Figure S1). A careful investigation on the trajectory confirms that formation of the dissociation state occurs through sliding of the ED molecule from top of the EA on to the graphene surface. However, the process of formation of the DS complex again needs a 180-degree flip of the ED molecule. In the dissociated state,TAB is found to be strongly engaged forming H-bonds either with the nitrogen atoms or nitrile groups present in the HAT-CN molecule. Moreover, a complete dissociation state is rarely found where all there is no formation of H-bonds. Both the DS and the dissociated state along with the TS is well sampled in the MetaD simulation, which suggests coordination numbers are already a suitable CV for this particular transition. At the end of the simulation, State A is found to be a stable minimum relative to the dissociated state B by


19


Page 20 of 41

2.18 kcal/mol with a barrier height of 7.61 kcal/mol going from state A to B (cf. Table 4). So, this simulation suggests that G…HAT-CN…TAB is a good candidate for DS complex, which seen previous cluster and crystal calculations.

Furthermore, we have studied MetaD simulation for G…HAT-CN…HMB DS complex, in where according to DFT results, HMB molecule does not form DS complex and prefer to stay over graphene surface. Figure 13 depicts the FES of the DS complex G…HAT-CN…HMB as a function of the Sstack and Sdiss coordination number after a 200-ns-long gas-phase WTMetaD simulation at room temperature (300 K). The FES again exhibits two highly pronounced minima: state A, corresponding to the DS state, and state B, corresponding to the dissociated state. State A occurs within the Sstack and Sdiss coordination numbers of 50–70 and 15–20 arbitrary units, respectively. For state B, both Sstack and Sdiss coordination numbers decrease to 0–20 and 0–10 arbitrary units, respectively, which again leads to the conclusion drawn in the previous DS complex. As the DS collapses, the coordination CV approaches zero and forms the dissociated state. In the dissociated state, HMB is stacked with graphene surface, but it interacts with HAT-CN via hydrophobic interactions. The main interactions occur between the methyl groups of HMB and the nitrile groups of HAT-CN. These contacts are relatively less stable than the previously observed HATCN…TAB H-bonding interactions of the dissociated state. Therefore, in this case, a complete dissociated state is formed. MetaD simulations have been successful in defining both the states along with TS. Finally, at the end, State B is the global minimum, where both ED and EA individually make stacking interactions with the periodic graphene sheet. State A is found to be 4.91 kcal/mol less stable than state B (cf. Table 4). The barrier height is calculated to be slightly over 5 kcal/mol going from state A to state B and it is lower by 2.57 kcal/mol when compared to the previous G…HAT-CN…TAB complex.


20

Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Based on MM interaction energy prediction, HMB is unstable to stay on the top layer of the DS complex and rather appears on the graphene surface. Therefore, the free energy calculated on the MM potential using the WTMetaD simulation is in good agreement with previous interactionenergy predictions. The FES of our final DS complex, G…HAT-CN…TMOB, after a 200-ns gas-phase WTMetaD simulation at 300 K is depicted in Figure 14 as a function of Sstack and Sdiss as the two CVs. The simulation again predicts two successive minima, areas A and B are the DS and dissociated states, respectively. The DS state A occurs within the Sstack and Sdiss coordination numbers of 5090 and 8– 10 arbitrary units, respectively. For state B, the dissociated state, both Sstack and Sdiss coordination numbers again decrease to 0–10 and 0–6 arbitrary units, respectively. Like in the previous cases, the sliding mechanism is again observed forming the dissociated state and mostly the hydrophobic contacts are formed between methoxy groups of TMOB and nitrile groups of HAT-CN. A close inspection also revealed that TMOB has sampled a wide range of stacking conformations due to its large size relative to TAB and HMB which results in wider range of coordination numbers in the Sstack state (see above). The positioning of the minima indicates that when the DS configuration collapses, the ED molecule gradually loses the stacking interaction with the EA, which leads to the formation of the dissociated state. The stacking interactions in states A and B can be envisaged in Figure S2. In this particular complex, the dissociated state is again the global minimum over DS state by mere 0.24 kcal/mol with a barrier height of 6.40 kcal/mol (cf. Table 4) going from state A to B. The complex barrier height has increased by 1.36 kcal/mol with respect to our previous DS complex, G…HAT-CN…HMB, with a significant decrease (4.67 kcal/mol) in the free-energy difference between the two states. Finally, the small difference between states A and B as well as the


21


Page 22 of 41

increased barrier heights may imply the formation of a stable DS complex at room temperature. This again follows the same trend as MM predicted interaction energies. For the very first time, the complex G…HAT-CN…TAB is shown to be more stable in the DS configuration state A than in the dissociated state B by 2.18 kcal/mol with a barrier height of 7.61 kcal/mol going from state A to B. The complex barrier height is increased by 1.21 kcal/mol with respect to the previous G…HAT-CN…TMOB complex with the stabilization of the DS state from 0.24 kcal/mol to 2.18 kcal/mol (cf. Table 4). Therefore, in this simulation study, large barrier heights stabilize the DS complexes. Overall, through the application of the WTMetaD simulation on the MM potential, the stability order has been successfully assigned for all the DS complexes. It has been predicted as: G…HAT-CN…TAB > G…HAT-CN…TMOB > G…HAT-CN…HMB. A comparison of the order of stability with the DFT results confirmed that they follow the same trend (cf. Table 2). So far, G…HAT-CN…TAB is found to be a good combination for the formation of stable DS complex. However, an increase in the concentration of TMOB could also result in a stable G…HAT-CN...TMOB DS complex, since the free energy difference between both the DS and dissociated state is very tiny. On the other hand, with HMB as electron donor, there is no possibility for the formation of G…HAT-CN....HMB DS complex. Conclusions In the present study, we have used sophisticated quantum mechanical techniques to understand the stabilization criterion of the ED-EA complexes at the 2D graphene surface. Our theoretical findings show potential implications to the design and fabricate nanostructures based on the noncovalent interaction. We have used cluster as well as crystal models to evaluate the stability of various DS complexes using DFT calculations. The stability of various double-sandwich structures has been verified. DFT calculations including long-range dispersion effects have demonstrated that dispersion energy and charge transfer play dominant roles in the formation of stable DS complexes. The strongest complexes are formed between subsystems with the highest donor and


22

Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


acceptor abilities. The probability of the formation of the DS complex is high for (C4...HATCN)...TAB and (C4...HAT-CN)...TTF. Furthermore, the results obtained from cluster-model calculations have been fully confirmed by crystal-model calculations. Finally, using the WTMetaD simulation on the MM potential, the stability of the DS complexes has been corroborated; they follow the same trend as the DFT results. It has been proven that the DS complex G…HATCN…TAB is more stable in DS configuration than the dissociated state B by 2.18 kcal/mol with a barrier height of 7.61 kcal/mol going from the associated to the dissociated form. This extensive study has offered the possibility of an unlimited variety of molecules with predictable functionalities. Associated Content Supporting Information: simulation protocol; well-tempered metadynamics; reweighted free energy surfaces (FES); one-dimensional potential of mean force (PMF); Band structures of Pristine Graphene alone and after depositions of the acceptor and donor molecules. Author Information Corresponding Author *E-mail: [email protected]. (P.H.). ORCID Pavel Hobza: 0000-0001-5292-6719 Vijay Madhav Miriyala: 0000-0003-3926-9691 Rabindranath Lo: 0000-0002-4436-3618 Susanta Haldar: 0000-0001-8642-109


23


Page 24 of 41

Amrit Sarmah: 0000-0001-6429-665 Acknowledgements This work was part of Research Project RVO: 61388963 of the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic. This work was supported by the Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development, and Innovations project “IT4 Innovations National Supercomputing Center LM2015070,” as well as from project LO1305 (P.H.). Notes The authors declare no competing financial interest References (1)

Novoselov, K. S.; Geim, A. K.; Morozov, S. V; Jiang, D.; Zhang, Y.; Dubonos, S. V; Grigorieva, I. V; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666–669.

(2)

Morozov, S. V.; Novoselov, K. S.; Katsnelson, M. I.; Schedin, F.; Elias, D. C.; Jaszczak, J. A.; Geim, A. K. Giant Intrinsic Carrier Mobilities in Graphene and Its Bilayer. Phys. Rev. Lett. 2008, 100, 016602.

(3)

Schedin, F.; Geim, A.; Morozov, S.; Hill, E.; Blake, P.; Katsnelson, M.; Novoselov, K. Detection of Individual Gas Molecules Adsorbed on Graphene. Nat. Mater. 2007, 6, 652– 655.

(4)

Stankovich, S.; Dikin, D. A.; Dommett, G. H. B.; Kohlhaas, K. M.; Zimney, E. J.; Stach, E. A.; Piner, R. D.; Nguyen, S. B. T.; Ruoff, R. S. Graphene-Based Composite Materials. Nature 2006, 442, 282–286.

(5)

Gao, W.; Alemany, L. B.; Ci, L.; Ajayan, P. M. New Insights into the Structure and Reduction of Graphite Oxide. Nat. Chem. 2009, 1, 403–408.

(6)

Williams, G.; Seger, B.; Kamt, P. V. TiO2-Graphene Nanocomposites. UV-Assisted Photocatalytic Reduction of Graphene Oxide. ACS Nano 2008, 2, 1487–1491.

(7)

Weiss, N. O.; Zhou, H.; Liao, L.; Liu, Y.; Jiang, S.; Huang, Y.; Duan, X. Graphene: An Emerging Electronic Material. Advanced Materials. 2012, 24, 5782–5825.

(8)

Bonaccorso, F.; Sun, Z.; Hasan, T.; Ferrari, A. C. Graphene Photonics and Optoelectronics. Nat. Photonics 2010, 4, 611–622.


24

Page 25 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9)

Yavari, F.; Koratkar, N. Graphene-Based Chemical Sensors. Journal of Physical Chemistry Letters. 2012, 3, 1746–1753.

(10)

Huang, X.; Zeng, Z.; Fan, Z.; Liu, J.; Zhang, H. Graphene-Based Electrodes. Advanced Materials. 2012, 24, 5979–6004.

(11)

Xiang, Q.; Yu, J.; Jaroniec, M. Graphene-Based Semiconductor Photocatalysts. Chemical Society Reviews. 2012, 41, 782–796.

(12)

Meyer, E. A.; Castellano, R. K.; Diederich, F. Interactions with Aromatic Rings in Chemical and Biological Recognition. Angewandte Chemie – International Edition. 2003, 42, 1210–1250.

(13)

Burley, S. K.; Petsko, G. A. Aromatic-Aromatic Interaction: A Mechanism of Protein Structure Stabilization. Science 1985, 229, 23–28.

(14)

Hong, B. H.; Lee, J. Y.; Lee, C. W.; Kim, J. C.; Bae, S. C.; Kim, K. S. Self-Assembled Arrays of Organic Nanotubes with Infinitely Long One-Dimensional H-Bond Chains [6]. Journal of the American Chemical Society. 2001, 123, 10748–10749.

(15)

Singh, N. J.; Lee, H. M.; Hwang, I. C.; Kim, K. S. Designing Ionophores and Molecular Nanotubes Based on Molecular Recognition. Supramolecular Chemistry. 2007, 19, 321– 332.

(16)

Georgakilas, V.; Otyepka, M.; Bourlinos, A. B.; Chandra, V.; Kim, N.; Kemp, K. C.; Hobza, P.; Zboril, R.; Kim, K. S. Functionalization of Graphene: Covalent and NonCovalent Approaches, Derivatives and Applications. Chemical Reviews. 2012, 112, 6156– 6214.

(17)

Haldar, S.; Kolář, M.; Sedlák, R.; Hobza, P. Adsorption of Organic Electron Acceptors onGraphene-like Molecules: Quantum Chemical and Molecular Mechanical Study. J. Phys. Chem. C, 2012, 116, 25328–25336.

(18)

Calborean, A.; Morari, C.; Maldivi, P. Combinded Molecular and Periodic DFT Aanalysisof the Adsorption of Co Macrocycles on Graphene. Journal of computational chemistry, 2018, 39, 130-138.

(19)

Voggu, R.; Rout, C. S.; Franklin, A. D.; Fisher, T. S.; Rao, C. N. R. Extraordinary Sensitivity of the Electronic Structure and Properties of Single-Walled Carbon Nanotubes to Molecular Charge-Transfer. J. Phys. Chem. C 2008, 112, 13053–13056.

(20)

Das, B.; Voggu, R.; Rout, C. S.; Rao, C. N. R. Changes in the Electronic Structure and Properties of Graphene Induced by Molecular Charge-Transfer. Chem. Commun. (Camb). 2008, 5155–5157.

(21)

Voggu, R.; Das, B.; Rout, C. S.; Rao, C. N. R. Effects of Charge Transfer Interaction of Graphene with Electron Donor and Acceptor Molecules Examined Using Raman Spectroscopy and Cognate Techniques. J. Phys. Condens. Matter 2008, 20, 472204.


25


Page 26 of 41

(22)

Subrahmanyam, K. S.; Voggu, R.; Govindaraj, A.; Rao, C. N. R. A Comparative Raman Study of the Interaction of Electron Donor and Acceptor Molecules with Graphene Prepared by Different Methods. Chem. Phys. Lett. 2009, 472, 96–98.

(23)

Varghese, N.; Ghosh, A.; Voggu, R.; Ghosh, S.; Rao, C. N. R. Selectivity in the Interaction of Electron Donor and Acceptor Molecules with Graphene and Single-Walled Carbon Nanotubes. J. Phys. Chem. C 2009, 113, 16855–16859.

(24)

Manna, A. K.; Pati, S. K. Tuning the Electronic Structure of Graphene by Molecular Charge Transfer: A Computational Study. Chem. – An Asian J. 2009, 4, 855–860.

(25)

Chen, L.; Wang, L.; Shuai, Z.; Beljonne, D. Energy Level Alignment and Charge Carrier Mobility in Noncovalently Functionalized Graphene. J. Phys. Chem. Lett. 2013, 4, 2158– 2165.

(26)

Das, A.; Sood, A. K.; Govindaraj, A.; Saitta, A. M.; Lazzeri, M.; Mauri, F.; Rao, C. N. R. Doping in Carbon Nanotubes Probed by Raman and Transport Measurements. Phys. Rev. Lett. 2007, 99, 136803.

(27)

Das, A.; Pisana, S.; Chakraborty, B.; Piscanec, S.; Saha, S. K.; Waghmare, U. V.; Novoselov, K. S.; Krishnamurthy, H. R.; Geim, A. K.; Ferrari, A. C.; et al. Monitoring Dopants by Raman Scattering in an Electrochemically Top-Gated Graphene Transistor. Nat. Nanotechnol. 2008, 3, 210–215.

(28)

Liu, H.; Liu, Y.; Zhu, D. Chemical Doping of Graphene. J. Mater. Chem. 2011, 21, 3335– 3345.

(29)

Antony, J.; Grimme, S. Structures and Interaction Energies of Stacked GrapheneNucleobase Complexes. Phys. Chem. Chem. Phys. 2008, 10, 2722–2729.

(30)

Szalay, P. S.; Galán-Mascarós, J. R.; Schottel, B. L.; Bacsa, J.; Pérez, L. M.; Ichimura, A. S.; Chouai, A.; Dunbar, K. R. Experimental and Computational Studies of ChargeTransfer and Reduction Products of 1,4,5,8,9,11-Hexaazatriphenylene-Hexacarbonitrile: HAT(CN)6. J. Clust. Sci. 2004, 15, 503–530.

(31)

M. Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus T. L.; de Jong, W. A. NWChem: A comprehensive and scalable open-source solution for large scale molecular simulationsComput. Phys. Commun. 2010, 181, 1477–1489.

(32)

TURBOMOLE V6.6 2014, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com

(33)

a) Grimme, S. Semiempirical GGA- type density functional constructed with a longrangedispersion correction. J. Comp. Chem., 2006, 27, 1787–1799. b) Mollenhauer, D.; Brieger, C.; Voloshina, E.; Paulus, B. Performance of Dispersion-Corrected DFT for the Weak Interaction between Aromatic Molecules and Extended Carbon-Based Systems. J. Phys. Chem. C 2015, 119, 1898−1904.

(34)

Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements HPu. J. Chem. Phys., 2010, 132, 154104.


26

Page 27 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(35)

NBO Version 3.1, E. D. Glendening, A. E. Reed, Eds.

(36)

Gaussian 09, Revision E.01, 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian, Inc., Wallingford CT, 2009.

(37)

Rezac, J. Cuby: An integrative framework for computational chemistry. J. Comput. Chem., 2016, 37, 1230.

(38)

Kathleen, L. The Structure of the Benzene Ring in Hexamethylbenzene. Proc. R. Soc. A. 1929, 123, 494–515.

(39)

Melander, M.; Jónsson, E. Ö.; Mortensen, J. J.; Vegge, T.; García, J. M.; Lastra Implementation of Constrained DFT for Computing Charge Transfer Rates within the Projector Augmented Wave Method. J. Chem. Theory Comput. 2016, 12 , 5367-5378.

(40)

Wu, Q.; Voorhis, T. V. Direct optimization method to study constrained systems within density-functional theory. Phys. Rev. A. 2005, 72, 024502.

(41)

Wu, Q.; Voorhis, T. V. Extracting electron transfer coupling elements from constrained density functional theory. J. Chem. Phys. 2006, 125, 164105.

(42)

Kaduk, B.; Kowalczyk, T.; Voorhis. T.V. Constrained Density Functional Theory. Chem. Rev. 2012, 112 , 321-370.

(43)

Henkelman, G.; Arnaldsson, A.; and Jónsson, H.; A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci. 2006, 36, 354-360.

(44)

Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source softwareproject for quantum simulations of materials, J. Phys. Condens. Matter. 2009, 21, 395502.


27


(45)

Page 28 of 41

Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations, Phys.

Rev. B. 1991, 43, 1993.

(46)

Becke, A. D. Density-functional exchange-energy approximation with correct asymptoticbehavior. Phys. Rev. A. 1988, 38, 3098–3100.

(47)

Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction, J. Comp. Chem. 2006, 27, 1787–1799.

(48)

Nguyen, C.V.; Iiyasov, V. V.; Nguyen, H. N. Tuning the electronic properties of armchair graphene nanoribbons by strain engineering. Phys. Scr. 2015, 90, 015802.

(49)

Ilyasov, V. V.; Meshi, B. C.; Nguyen, V. C.; Ershov, I. V.; Nguyen, D. C. Magnetism and transport properties of zigzag graphene nanoribbons/hexagonal boron nitride heterostructures. J. Appl. Phys. 2014, 115, 053708 .

(50)

Ilyasov, V. V.; Nguyen, C. V.; Ershov, I. V.; Nguyen, C. D.; Hieu, N. N. Modulation of the band structure in bilayer zigzag graphene nanoribbons on hexagonal boron nitride using the force and electric fields. Mater. Chem. Phys. 2015, 154, 78–83.

(51)

Ilyasov, V. V.; Meshi, B. C.; Nguyen, V. C.; Ershov, I. V.; Nguyen, D. C. Tuning the band structure, magnetic and transport properties of the zigzag graphene nanoribbons/hexagonal boron nitride heterostructures by transverse electric field. J. Chem. Phys. 2014, 141, 014708.

(52)

Ilyasov, V. V.; Nguyen, C. V.; Ershov, I. V.; Hieu, N. N. Electric field and substrateinduced modulation of spin-polarized transport in graphene nanoribbons on a3b5 semiconductors. J. Appl. Phys. 2015, 117, 174309.

(53)

Ilyasov, V. V.; Nguyen, C. V.; Ershov, I. V.; Hieu, N. N.; Effect of electric field on the electronic and magnetic properties of a graphene nanoribbon/aluminium nitride bilayer system. RSC Adv. 2015, 5, 49308–49316.

(54)

Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development andTesting of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157–1174.

(55)

Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graph. Model. 2006, 25, 247– 260.


28

Page 29 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(56)

Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269–10280.

(57)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 2001, 926 (May 2016).

(58)

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 1–4.

(59)

Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604–613.

(60)

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447.

(61)

Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119, 736–742.

(62)

Haldar, S.; Spiwok, V.; Hobza, P. On the Association of the Base Pairs on the SilicaSurface Based on Free Energy Biased Molecular Dynamics Simulation and Quantum Mechanical Calculations. J. Phys. Chem. C, 2013, 117, 11066– 11075.

(63)

Haldar, S.; Kührová, P.; Banáš, P.; Spiwok, V.; Šponer, J.; Hobza, P. Otyepka. M, Insightinto Stability and Folding of GNRA and UNCG Tetraloops Revealed by Microsecond Molecular Dynamics and Well Tempered Metadynamics. J. Chem. Theory Comput. 2015, 11, 3866–3877.


29


Page 30 of 41

Tables and Figures Table 1: The ionization potentials (IP) and electron affinities (EA) of electron donors and acceptors in eV, and Vs,max from the ESP of electron acceptors.

Monomer Ionization Potential (IP) Electron Affinity (EA) Acceptors (EA) HAT-CN

Donors (ED)

Vs,max (kcal/mol)

9.9

3.6

50.7

TACN

11.3

2.5

49.3

TANO

10.3

2.9

58.6

HCB

10.7

3.3

47.5

HMB

7.5

-1.7

TMOB

7.4

-2.0

TAB

6.5

-2.9

TTF

6.2

-0.8

Table 2: The DFT-D3 interaction energies (ΔE) of the complexes in kcal/mol, DFT-D3 dispersion energies Edisp in kcal/ mol, the distances between the monomers of the complexes and their centroids, charge transfer in electrons from NBO analysis, HOMO-LUMO gaps (H-L GAP) in eV, and the difference in the stabilization energies of (C4...EA)...ED and C4...ED complexes ∆E’ in kcal/mol.

Complex

ΔE

Edisp

Distance (Å)

NBO CT

H-L GAP

C4...HAT-CN

-46.5

-68.2

3.360

-0.22

0.25

C4...HCB

-30.6

-43.8

3.536

-0.19

0.27

C4...TACN

-17.7

-22.6

3.373

-0.09

0.34

C4...TANO

-22.9

-27.1

3.319

-0.11

0.25

C4...TMOB

-20.4

-34.2

3.574

-0.01

1.36

C4...HMB

-23.1

-40.2

3.493

-0.00

1.36

C4...TAB

-22.8

-31.9

3.135

-0.03

0.80

C4...TTF

-21.5

-35.0

4.147

-0.03

0.57

(C4...HAT-CN)...TMOB

-22.9

-30.6

3.115

-0.06

0.31


ΔE’

-2.5

30

Page 31 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(C4...HAT-CN)...HMB

-21.6

-32.2

3.326

-0.05

0.34

1.5

(C4...HAT-CN)...TAB

-26.7

-25.2

3.360

-0.14

0.40

-3.9

(C4...HAT-CN)...TTF

-26.0

-32.2

3.507

-0.32

0.36

-4.5

(C4...HCB)...TMOB

-17.3

-23.2

3.421

-0.07

0.36

3.1

(C4...HCB)...HMB

-19.0

-27.5

3.597

-0.08

0.34

4.1

(C4...HCB)...TAB

-25.9

-21.9

3.276

-0.20

0.42

-3.1

(C4...HCB)...TTF

-20.8

-24.7

3.276

-0.23

0.35

0.7

Table 3: Charge transfer in electrons between the donor top layer and the rest of the complex calculated from DFT Lowdin population analysis. Bader charges calculated at the same level theory are also presented here for comparison. ΔE’’ in kcal/mol is the difference between the interaction energies of DFT and cDFT of the DS complexes.

System

ΔE’’

Lowdin Charges

Bader Charges

(C4...HAT-CN)...TAB

5.3

-0.16

-0.18

(C4...HAT-CN)...HMB

0.6

-0.06

-0.10

(C4...HAT-CN)...TMOB

0.8

-0.07

-0.09

(C4...HAT-CN)...TTF

6.4

-0.31

-0.35

Table 4: The stability of the DS complexes found in the metadynamics simulation. Free energies are shown in kcal/ mol. Interaction energies in kcal/mol calculated using AMBER are also presented here.

System

DS

Collapsed

TS#b

AMBER IE

G…HAT-CN…HMB

4.91

0a

5.04

-22.0

G…HAT-CN…TMOB

0.24

0

6.40

-21.0

0

2.18

7.61

-26.5

G…HAT-CN…TAB a

A state with the energy of 0 kcal/mol is considered the global minimum. bThe height of the TS going from DS to the collapsed state.


31


Page 32 of 41

Figure 1. The molecular representation of various molecules used in this work; the model graphene system: a) C4; electron acceptors: b) HAT-CN, c) HCB, d) TACN, e) TANO; and electron donors: f) TMOB, g) HMB, h) TAB, i) TTF.


32

Page 33 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Vs, max = 50.7 kcal/mol

c)


d)



Figure 2. Molecular Electrostatic potential (ESP) and the Vsmax (in kcal/mol) values of selected electron acceptors, a) 1,4,5,8,9,11-hexaazatriphenylenehexacarbonitrile (HAT-CN), b) hexacyanobenzene (HCB), c) 1,3,5-triazine2,4,6tricarbonitrile (TACN) and d) 2,4,6-trinitro-1,3,5-triazine (TANO).


33


Page 34 of 41

Figure 3. Plot of ionization potentials (IP) and electron affinities (EA) in eV of various donor and acceptor subsystems used in this work.


34

Page 35 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.360

3.536

3.373

3.319

Figure 4. The optimized geometries of circumcircumcoronene C4 complexed with a) 1,4,5,8,9,11hexaazatriphenylenehexacarbonitrile (HAT-CN) b) hexacyanobenzene (HCB), c) 1,3,5-triazine2,4,6tricarbonitrile (TACN) and d) 2,4,6-trinitro-1,3,5-triazine (TANO). The distances in Å between the centroids are given next to the dotted line. [C: grey, N: blue, O: red, H: white]

a)

Figure 5. The optimized geometries of circumcircumcoronene (C4) complexed with (a) trimethoxybenzene (TMOB), (b) hexamethylbenzene (HMB), (c) triaminobenzene (TAB) and (d) tetrathiafulvalene (TTF). The distances in Å between the centroids are given next to the dotted line. [C: grey, N: blue, O: red, H: white, S: yellow]

3.326 3.234

3.580



Page 36 of 41

3.507

3.360

3.592

3.178

Figure 6. The optimized DS complex geometries of (a) C4…HAT-CN…TMOB, (b) C4…HAT-CN… HMB, (c) C4…HAT-CN…TAB and (d) C4…HAT-CN…TTF. The distances in Å between the centroids are given next to the dotted line. [C: grey, N: blue, O: red, H: white, S: yellow]

3.597

3.421 3.503

3.276

3.276

3.388

3.580

Figure 7. The optimized DS complex geometries of (a) C4…HCB…TMOB, (b) C4…HCB…HMB, (c) C4…HCB…TAB and (d) C4…HCB…TTF. The distances in Å between the centroids are given next to the dotted line. [C: grey, N: blue, O: red, H: white, S: yellow]


Page 37 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. Plot of interaction energies in kcal/mol of stable sandwich and double sandwich complexes.

Figure 9. The total density of states (DOS) spectrum of the three different DS systems (G...HATCN)...TAB, (G...HCB)...TAB and (G...HAT-CN)...TTF. The Fermi level is set at zero.



Figure 10. The projected density of states (PDOS) spectrum of the three different DS systems (G...HATCN)...TAB, (G...HCB)...TAB and (G...HAT-CN)...TTF. The various contributions from the different types of atoms are represented in different shades. The Fermi level is set at zero.

Figure 11. The pictorial representations of the molecular orbital densities for highest occupied and the lowest unoccupied states for the three systems (G...HAT-CN)...TAB, (G...HCB)...TAB and (G...HATCN)...TTF along with their HOMO-LUMO energy gaps (in eV).


Page 38 of 41

Page 39 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12. The FES constructed on the gas-phase WTMetaD simulation for the DS complex G…HATCN…TAB as a function of the Sstack and Sdiss, the coordinations number as CVs (a). The one-dimensional PMF projected on the Sstack coordination number for the calculation of the barrier heights (b). The simulation time is 200 ns long, which results in the DS state being the global minimum over the dissociated state, with the difference in the free energy being 2.18 kcal/mol.



Figure 13. The FES constructed on the gas-phase WTMetaD simulation for the DS complex G…HATCN…HMB as a function of the Sstack and Sdiss, the coordination numbers as CVs (a). The onedimensional potential of mean force (PMF) projected on the Sstack coordination number (b). The simulation time is 200 ns, which results in the dissociated state as the global minimum over the DS complex.

Figure 14. The FES constructed on the gas-phase WTMetaD simulation for the DS complex G…HATCN…TMOB as a function of the Sstack and Sdiss, the coordination numbers as CVs (a). The PMF projected on the Sstack coordination number, which shows the barrier height between the two states (b). The simulation time here is again 200 ns longer, which results in the dissociated state being the global minimum, with the difference in the free energy being only 0.24 kcal/mol over DS state.


Page 40 of 41

Page 41 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Structure and Properties of Double-Sandwich Complexes at Graphene

Recommend Documents