Semi-Empirical Model for Nano-Scale Device Simulations
Kurt Stokbro, Dan Erik Petersen, Søren Smidstrup, Anders Blom, Mads Ipsen, Kristen Kaasbjerg
SSemi-Empirical Model for Nano-Scale Device Simulations
Kurt Stokbro, ∗ Dan Erik Petersen, Søren Smidstrup, Anders Blom, and Mads Ipsen
QuantumWise A/S,Nørre Søgade 27A, 1. th,DK-1370 Copenhagen K, Denmark
Kristen Kaasbjerg
Center for Atomic-scale Materials Design (CAMd),Department of Physics, Technical University of Denmark,DK-2800 Kgs. Lyngby, Denmark (Dated: October 27, 2018)We present a new semi-empirical model for calculating electron transport in atomic-scale devices.The model is an extension of the Extended H¨uckel method with a self-consistent Hartree potential.This potential models the effect of an external bias and corresponding charge re-arrangementsin the device. It is also possible to include the effect of external gate potentials and continuumdielectric regions in the device. The model is used to study the electron transport through an organicmolecule between gold surfaces, and it is demonstrated that the results are in closer agreement withexperiments than ab initio approaches provide. In another example, we study the transition fromtunneling to thermionic emission in a transistor structure based on graphene nanoribbons.
PACS numbers: 73.40.-c, 73.63.-b, 72.10.-d, 72.80.Vp
I. INTRODUCTION
As the minimum feature sizes of electronic devices areapproaching the atomic scale, it becomes increasinglyimportant to include the effects of single atoms in de-vice simulations. In recent years, there have been sev-eral developments of atomic-scale electron transport sim-ulation models based on the Non-Equilibrium Green’sFunction (NEGF) formalism . The approaches canroughly be divided into two catagories: ab initio ap-proaches, where the electronic structure of the systemis calculated from first principles, typically with Den-sity Functional Theory (DFT) , and semi-empirical ap-proaches, where the electronic structure is calculated us-ing a model with adjustable parameters fitted to experi-ments or first-principles calculations. Examples of semi-empirical transport models are methods based on Slater-Koster tight-binding parameters and Extended H¨uckelparameters .The ab initio models have the advantage of predictivepower, and can often give reasonable results for systemswhere there is no prior experimental data. However, theuse of the Kohn-Sham one-particle states as quasiparti-cles is questionable, and it is well known that for manysystems the energies of the unoccupied levels are ratherpoorly described within DFT. Furthermore, solving theKohn-Sham equations can be computationally demand-ing, and solving for device structures with thousands ofatoms is only feasible on large parallel computers.The semi-empirical models have less predictive power,but when used within their application domain they cangive very accurate results. The models may also be fit-ted to experimental data, and can thus in some cases givemore accurate results than DFT-based methods. How-ever, the main advantage of the semi-empirical methods are their lower computational cost.In this paper we will present the formalism behind anew semi-empirical transport model based on the Ex-tended H¨uckel (EH) method. The model can be viewed asan extension to the work by Zahid et al. , with the maindifference being the treatment of the electrostatic interac-tions. Zahid et al. only describe part of the electrostaticinteractions in the device; most importantly, they usethe Fermi level of the electrodes as a fitting parameterand do not account for the charge transfer from the elec-trodes to the device. In the current work, the Fermi levelof the electrodes is determined self-consistently by usingthe methodology introduced by Brandbyge et al. In thisway, we include the charge transfer from the electrodesto the device region and describe all electrostatic interac-tions self-consistently. This is accomplished by defininga real-space electron density and numerically solving forthe Hartree potential on a real-space grid. Through amulti-grid Poisson solver, we include the self-consistentfield from an applied bias, and allow for including con-tinuum dielectric regions and electrostatic gates withinthe scattering region.The organization of the paper is the following: In sec-tion II we introduce the self-consistent Extended H¨uckel(EH-SCF) model, and in section III we present the for-malism for modelling nano-scale devices. In section IVwe apply the model to a molecular device, and in sec-tion V we consider a graphene nano-transistor where anelectrostatic gate is controlling the electron transport inthe device. Finally, in section VI, we conclude the paper. a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r II. THE SELF-CONSISTENTEXTENDED-H ¨UCKEL METHOD
In this section we describe the EH-SCF framework. InExtended H¨uckel theory, the electronic structure of thesystem is expanded in a basis set of local atomic orbitals(LCAOs) φ nlm ( r ) = R nl ( r ) Y lm (ˆ r ) , (1)where Y lm is a spherical harmonic and R nl is a superpo-sition of Slater orbitals R nl ( r ) = r n − − l (2 n )! (cid:2) C (2 η ) n +1 e − η r + C (2 η ) n +1 e − η r (cid:3) . (2)The LCAOs are described by the adjustable parameters η , η , C , and C , and these parameters must be definedfor the valence orbitals of each element.The central object in EH theory is the overlap matrix, S ij = (cid:26) δ ij if R i = R j (cid:82) V φ i ( r − R i ) φ j ( r − R j ) d r if R i (cid:54) = R j (3)where i is a composite index for nlm and R i is the posi-tion of the center of orbital i .From the overlap matrix, the one-electron Hamiltonianis defined by H ij = (cid:26) E i + δV H ( R i ) if i = j ( β i + β j )( E i + E j ) S ij + ( δV H ( R i ) + δV H ( R j )) S ij if i (cid:54) = j (4)where E i is an orbital energy, and β i is an adjustable pa-rameter (often chosen to be 1.75). δV H ( R i ) is the Hartreepotential corresponding to the induced electron densityon the atoms, i.e. the change in electron density com-pared to a superposition of neutral atomic-like densities.This term must be determined self-consistently, and isnot included in standard EH models . In the followingsection we describe how this term is calculated. A. Solving the Poisson Equation to Obtain theHartree Potential
To calculate the induced Hartree potential we need todetermine the spatial distribution of the electron den-sity. To this end, we introduce the Mulliken populationof atom number µm µ = (cid:88) i ∈ µ (cid:88) j D ij S ij , (5)where D ij is the density matrix. The total number ofelectrons can now be written as a sum of atomic contri-butions, N = (cid:80) µ m µ .We will represent the spatial distribution of eachatomic contribution by a Gaussian function, and use the following approximation for the spatial distribution of theinduced electron density: δn ( r ) = (cid:88) µ δm µ (cid:114) α µ π e − α µ | r − R µ | , (6)where the weight δm µ = m µ − Z µ of each Gaussian isthe excess charge of atom µ as obtained from the Mul-liken population m µ and the ion valence charge Z µ . Sub-sequently, the Hartree potential is calculated from thePoisson equation − ∇ · [ (cid:15) ( r ) ∇ δV H ( r )] = δn ( r ) , (7)which is solved with the appropriate boundary conditionson the leads and gate electrodes imposed by the appliedvoltages. Here, (cid:15) ( r ) is the spatially dependent dielectricconstant of the device constituents, and allows for theinclusion of dielectric screening regions.To see the significance of the width α µ of the Gaussianorbital, let us calculate the electrostatic potential froma single Gaussian electron density at position R µ . Theresult is δV H ( r ) = e ( m µ − Z µ ) Erf( √ α µ | r − R µ | ) | r − R µ | , (8)and from this equation we see that the on-site value ofthe Hartree potential is δV H ( R µ ) = ( m µ − Z µ ) γ µ , wherethe parameter γ µ = 2 e (cid:114) α µ π (9)is the on-site Hartree shift. The parameter γ µ is a well-known quantity in CNDO theory , and values of γ µ are listed for many elements in the periodic table. Thus,we fix the value of γ µ using CNDO theory, and then useEq. (9) to calculate the value of α µ for each element. III. EH-SCF METHOD FOR A NANO-SCALEDEVICE
Fig. 1 illustrates the setup of a molecular device sys-tem. The system consists of three regions: the centralregion, and the left and right electrode regions. The cen-tral region includes the active parts of the device andsufficient parts of the contacts, such that the propertiesof the electrode regions can be described as bulk materi-als. For metallic contacts, this will typically be achievedby extending the central region 5–10 ˚A into the contacts.The calculation of the electron transport properties ofthe system is divided into two parts. The first part isa self-consistent calculation for the electrodes, with peri-odic boundary conditions in the transport direction. Inthe directions perpendicular to the transport direction,we apply the same boundary conditions for the two elec-trodes and the central region, and these boundary con-ditions are described below.
FIG. 1. (Color online) Geometry of a nano-device consist-ing of a dithiol-triethynylene-phenylene molecule attached totwo (3x3) (111) gold electrodes. The left and right electroderegions are illustrated with wire boxes, and the propertiesof these regions are obtained from a calculation with periodicboundary conditions in all directions. The region between thetwo electrodes is the central device region, which is describedwith open boundary conditions in the transport direction, andperiodic boundary conditions in the directions perpendicularto the transport direction.
In the second part of the calculation, the electrodesdefine the boundary conditions for a self-consistent openboundary calculation of the properties of the central re-gion. The main steps in the open boundary calculation isthe determination of the density matrix, the evaluationof the real-space density, and, finally, the calculation ofthe Hartree potential. These steps will be described inmore detail in the following section.
A. Calculating the Self-Consistent Density Matrixof the Central Region
In this section we will describe the calculation of thedensity matrix of the central region. We assume that theself-consistent properties of the left and right electrodeshave already been obtained, and thus we also know theirrespective Fermi levels, ε FL and ε FR . We allow for an exter-nal bias V b to be applied between the two electrodes, anddefine the left and right chemical potentials µ L = ε FL − eV b and µ R = ε FL . The applied bias thus shifts all energiesin the left electrode, and a positive bias gives rise to anelectrical current from left to right.The density matrix for this non-equilibrium system,with two different chemical potentials, is found by fillingup the left and right originating states according to theirrespective chemical potentials ,ˆ D = (cid:90) ∞−∞ (cid:2) ˆ ρ L ( ε ) n F ( ε − µ L ) + ˆ ρ R ( ε ) n F ( ε − µ R ) (cid:3) dε, (10)where ˆ ρ L (ˆ ρ R ) is the contribution to the spectral den-sity of states from scattering states originating in the left(right) reservoir.The calculation of the spectral densities is performedusing NEGF theory, and we write the partial spectraldensities as ˆ ρ L,R ( ε ) = 12 π ˆ G ( ε ) ˆΓ L,R ( ε ) ˆ G † ( ε ) , (11)where ˆ G is the retarded Green’s function of the centralregion, and the broadening function ˆΓ = i [ ˆΣ − ˆΣ † ] is given by the self energies ˆΣ L and ˆΣ R , which arise due to thecoupling of the central region with the semi-infinite leftand right electrodes, respectively.Further details of the NEGF formalism can be foundin Refs. 1,4. Here we just note that to improve the nu-merical efficiency, the integral in Eq. (10) is divided intoan equilibrium and non-equilibrium part. The equilib-rium part is calculated on a complex contour far from thereal-axis poles of the Green’s function, and only the non-equilibrium part is performed along the real axis. Theequilibrium and non-equilibrium parts are then joinedusing the double-contour technique introduced by Brand-byge et al. From the density matrix we may now evaluate the real-space density in the central region using Eq. (6). It isimportant to note that near the left and right faces ofthe central region there will be contributions from theelectrode regions, and this “spill in” must be properlyaccounted for.Once the real-space density is known, the Hartree po-tential is calculated by solving the Poisson equation inEq. (7) using a real-space multi-grid method. On the leftand right faces of the central region the Hartree potentialis fixed by the electrode Hartree potentials, appropriatelyshifted according to the applied bias. In the directionsperpendicular to the transport directions, we apply theappropriate boundary conditions, fixed or periodic, asdemanded by e.g. the presence of gate electrodes.The so-obtained Hartree potential defines a newHamiltonian, via Eq. 4, and the steps in section III Amust be repeated until a self-consistent solution is ob-tained.
B. Transmission and Current
Once the self-consistent one-electron Hamiltonian hasbeen obtained, we can finally evaluate the transmissioncoefficients T ( ε ) = Tr[ˆΓ L ( ε ) ˆ G † ( (cid:15) )ˆΓ R ( ε ) ˆ G ( ε )] (12)and the current I = 2 eh (cid:90) ∞−∞ T ( ε )[ n F ( ε − µ L ) − n F ( ε − µ R )] dε. (13)In the following sections, we apply this formalism tothe calculation of the electrical properties of a moleculebetween gold electrodes, as well as a graphene nano-transistor. IV. TOUR WIRE BETWEEN GOLDELECTRODES
In this section we will investigate the electrical prop-erties of a phenylene ethynylene oligomer, also popularlycalled a Tour wire. We will compare the electrical proper-ties of the molecule when it is symmetrically and asym-metrically coupled with two Au(111) surfaces. In thesymmetric system, as illustrated in Fig. 1, the moleculeis connected with both gold surfaces through thiol bonds,whereas the asymmetric system only has a thiol bond toone of them.The system has previously been investigated exper-imentally by Kushmerick et al. and theoretically byTaylor et al. , and it has been found that the asymmet-rically coupled system shows strongly asymmetric I–Vcharacteristics .The calculations by Taylor et al. were based on DFT-LDA, and the asymmetric behaviour could be related tothe voltage drop in the system. This system is thereforean excellent testing ground for our semi-empirical model,since a correct description of the electrical properties re-quires not only a good model for the zero-bias electronicstructure, but also a good description of the bias-inducedeffects. A. Transmission Spectrum of the Symmetric TourWire Junction
To setup the symmetric system we first relaxed theisolated Tour wire using DFT-LDA . During the re-laxation, passivating hydrogen atoms were kept on thesulfur atoms. Afterwards, these hydrogen atoms were re-moved and the two sulfur atoms placed at the FCC sitesof two Au(111)-(3x3) surfaces. The height of the S atomabove the surface was 1.9 ˚A (corresponding to an Au–Sdistance of 2.53 ˚A).We next set up the EH model with Hoffmannparameters and perform a self-consistent calculationto obtain transmission spectra for different k-point sam-pling grids. The results are shown in the upper plot ofFig. 2. In each case, the same k-point grid was used forboth the self-consistent and transmission calculation, andwe see from the figure that using (1x1) k-point is insuffi-cient while (2x2) and (4x4) k-points give almost identicalresults. Thus, we will use a (2x2) k-point sampling gridfor the remainder of this study.In the lower plot of Fig. 2 we compare the transmis-sion spectra calculated with DFT-LDA, EH-SCF, andEH without the Hartree term of Eq. (4). For the DFT-LDA model we use similar parameters as Taylor et al. ,except for the k-point sampling which is (2x2) in the cur-rent study. The calculations in Ref. 21 were performedwith a (1x1) k-point sampling, which is insufficient , andthus the DFT-LDA results in this study will differ fromthose by Taylor et al. For the EH calculation we see a peak in the transmis-sion spectrum just around the Fermi level of the goldelectrodes. This peak arises from transmission throughthe LUMO orbital of the Tour wire.In the self-consistent EH calculation there will be acharge transfer from the gold surface to the LUMO or-
FIG. 2. (Color online) The upper plot shows the transmis-sion spectrum of the symmetric Tour wire device, calculatedwith the EH-SCF model for three different k-point samplings.The lower plot shows transmission spectra calculated witha (2x2) k-point sampling using different models: EH-SCF(solid), EH without the Hartree term (dotted), and DFT-LDA (blue dashed). Energies are given relative to the Fermilevel of the gold electrodes. bital, and we see that this gives rise to a shift of theorbital by 1 eV, illustrated by the arrow in Fig. 2.For the DFT-LDA calculation we see that the LUMOpeak is shifted further away from the gold Fermi level,and the HOMO and LUMO peaks of the transmissionspectrum are placed almost symmetrically around thegold Fermi level. We also note that the transmission co-efficient at the Fermi level, corresponding to the zero-biasconductance, is almost one order of magnitude higherwithin the DFT-LDA model. We will discuss this fur-ther below.We also note that Taylor et al. find a LUMO level evenfurther away from the gold Fermi level; this is related tothe insufficient k-point sampling.
B. I–V Characteristics of the Symmetric andAsymmetric Tour Wire Systems
We will now study both the symmetric and asymmet-ric Tour wire system and compare their respective I–Vcharacteristics. The geometry of the asymmetric system
FIG. 3. (Color online) Geometry of the asymmetric system.The Tour wire is attached to the left gold electrode througha thiol bond, while the right end of the molecule is hydrogen-terminated and there is no chemical bond to the right goldelectrode.FIG. 4. (Color online) I–V characteristics of the symmetric(upper figure) and asymmetric (lower figure) Tour wire device.The positive current direction is from left to right. is illustrated in Fig. 3. The geometry is similar to thatof Fig. 1, except for the right-most sulfur atom whichhas been replaced by a hydrogen atom with a C–H bondlength of 1.1 ˚A. The distance between the hydrogen atomand the right gold surface is 1.5 ˚A.We perform self-consistent calculations for both thesymmetric and asymmetric systems with the EH-SCFand DFT-LDA methods, and vary the bias from –1 to+1 V in steps of 0.1 V. The results are shown in Fig. 4.For the symmetric device we obtain rather similar, sym-metric I–V characteristics for both the EH-SCF andDFT-LDA methods. The main difference is that the zero-bias conductance is significantly higher with DFT-LDA,reflecting the higher transmission coefficient at the Fermilevel, as shown in Fig. 2.For the asymmetric device we see that both the DFT-LDA and EH-SCF models give rise to rectification – how-
FIG. 5. (Color online) Voltage drop of the symmetric andasymmetric Tour wire systems along a line that goes throughthe two sulfur atoms in the symmetric system, for an ap-plied bias of +1 V. The inset shows the voltage drop in theasymmetric system subtracted from the voltage drop in thesymmetric system. Both plots show results calculated withthe EH-SCF (solid) and DFT-LDA (dashed) models. ever, in opposite directions. Taylor et al. demonstratedthat the rectification was related to the voltage drop inthe system, and we therefore in Fig. 5 compare the volt-age drops obtained with the two methods. The EH-SCFvoltage drop is smooth, since the charge density is com-posed of a superposition of single rather broad Gaussianson each atom. The DFT-LDA model shows atomic-scaledetails, however, as illustrated by the inset, the relativechange in the voltage drop between the asymmetric andsymmetric system is quite similar for the EH-SCF andDFT-LDA methods. Both methods reveal that in theasymmetric system there is an additional voltage drop atthe contact with the weak bond. This is also one of themain results of Taylor et al. .The additional voltage drop at the weak contact meansthat the molecular levels of the Tour wire mainly fol-low the electrochemical potential of the right electrode .Since the voltage drop is similar for the DFT-LDA andEH-SCF models, the difference in the I–V characteris-tics must be related to the different electronic structureat zero bias in the two models. Within the DFT-LDAmodel, the transport at the Fermi level is dominated bythe HOMO. At negative bias, the left electrode has ahigher electrochemical potential, and electrons from theoccupied HOMO level can propagate to empty states inthe left electrode. Thus, for the DFT-LDA model, thecurrent is highest for a negative bias at the left electrode.For the EH-SCF model, on the other hand, the transportat the Fermi level is dominated by the LUMO, and thecurrent in this case is highest for a positive bias at theleft electrode.Comparing with the experimental results of Kushmeric et al. , we find that the EH-SCF rectification directionagrees with the experimental rectification direction, while FIG. 6. (Color online) Graphene nano-transistor consisting oftwo metallic zigzag nanoribbons connected by a semiconduct-ing armchair ribbon. The nanoribbons are passivated withhydrogen, and the width of the ribbons are is 7 ˚A. The deviceis sitting on top of a dielectric and the transport is controlledby an electrostatic back-gate. The contour plot illustrates theHartree potential for a gate potential of –1 V. the DFT-LDA model predicts rectification in the oppo-site direction. We note that the rectification directionobtained with our DFT-LDA model is similar to the re-sults of Taylor et al. .Thus, for this system it seems that the EH-SCF modelis in better agreement with the experimental results,compared to the DFT-LDA model. The example showsthat the EH-SCF model gives a very good descriptionof both the electronic structure and the voltage drop inthe system. The comparisons between the two methodsalso illustrates how small variations in the positions ofthe HOMO and LUMO levels may change the electricalproperties of the Tour wire device. V. Z-SHAPED GRAPHENENANO-TRANSISTOR
In this section we will compare the electrical proper-ties of a short (34 ˚A) and long (86 ˚A) graphene nano-transistor. The system consists of two electrodes consist-ing of metallic, zigzag-edge graphene nanoribbons con-nected through a semiconducting armchair-edge centralribbon. The system is placed 1.4 ˚A above a dielectricmaterial with dielectric constant (cid:15) = 4 (cid:15) , correspondingto SiO . The dielectric is 3 ˚A thick, and below the di-electric there is an electrostatic gate. The geometry ofthe short system is illustrated in Fig. 6. A similar systemwas investigated by Yan et al. using DFT-LDA.For the calculation we use EH parameters from Ref. 13which were derived by fitting to a reference band struc-ture of a graphene sheet calculated with DFT-LDA. Withthese parameters, we find a band gap of the central rib-bon of 2.2 eV, in agreement with DFT-LDA calculations,which illustrates the transferability of the EH parametersfrom 2D graphene to a 1D graphene nanoribbon. A. Transmission Spectrum
Fig. 7 shows the transmission spectrum for both thelong and the short system when there is no applied bias
FIG. 7. Zero-bias transmission spectrum for the short(dashed) and long (solid) graphene device. Energies are rela-tive to the Fermi level of the electrodes. and zero gate potential. The shape of the transmissionspectrum is directly related to the electronic structure ofthe central semiconducting ribbon.The transmission is strongly reduced in the energy re-gion from –0.7 to 1.5 eV, corresponding to the band gapof the central armchair ribbon. Since there are no energylevels in this interval, the electrons must tunnel in orderto propagate across the junction. For the longer devicethe electrons must tunnel a longer distance, and thus thetransmission is more strongly reduced.Outside the band gap, the transmission is close to 1and shows a number of oscillations. Since the centralribbon has a finite length, it resembles a molecule with anumber of discrete energy levels. The levels give rise topeaks in the transmission spectrum, and since the longersystem has more energy levels, the peaks are more closelyspaced there.In the following section we will see how this differencein the transmission spectrum gives rise to qualitativelydifferent transport mechanisms in the two devices.
B. Transistor Characteristics
We now calculate the current for an applied source–drain voltage of 0.2 V as a function of the applied gatepotential. Fig. 8 shows the current for the long and shortdevices, respectively, for gate potentials in the range –1 to1 V, for different electrode temperatures. We see that forthe short device there is only a small effect of the gatepotential and electron temperature, while for the longdevice the conductance falls off exponentially, reaching aminimum in the range 0 to 0.5 V. Moreover, the currentis strongly temperature-dependent.The lack of temperature-dependence for the short de-vice shows that the transport is completely dominatedby electron tunneling. For the long device, on the
FIG. 8. (Color online) The tunneling current for a source–drain voltage of 0.2 V, as a function of the gate potentialfor the long (solid) and short (dashed) graphene device, re-spectively. Three different values of the electron temperaturein the electrodes were considered: 150 K, 300 K and 450 K.The dotted lines illustrate the 1 /k B T slope for the differenttemperatures. other hand, there is a strong temperature dependence,and in this case the electron transport is dominatedby thermionic emission. The dotted lines illustrate the1 /k B T slope expected for thermionic emission. We seethat in the gate voltage range from –0.25 to –0.75 V, theI–V characterics follow these slopes well.Fig. 6 also shows the electrostatic profile through thedevice. We see that the gate potential is almost per-fectly screened by the graphene ribbon, i.e. the gate po-tential does not penetrate through the central ribbon.This means that for a layered structure, only the firstlayer would be strongly affected by the back-gate. Thishas some implications also for gated nanotube devices.In such a device, only the atoms facing the gate elec-trode will be strongly influenced, and this explains whyin Ref. 26 we found that the transport in the device wasdominated by tunneling even though the nanotube was110 ˚A long, and thus longer than the graphene junctionsstudied in this paper. Thus, to obtain efficient gatingof a nanotube, the gate electrode must wrap around thetube. VI. CONCLUSIONS
In this paper we have introduced a new semi-empiricalmodel for electron transport in nano-devices. The modelis based on the Extended H¨uckel method that extendsthe work by Zahid et al. to give a more complete de-scription of the electrostatic interactions in the device. Inparticular, the position of the electrode Fermi level andthe charge transfer between the contacts and the deviceare calculated self-consistently.Compared to DFT-based transport methods, the mainadvantage of our new method is that it is computationallyless expensive, as well as having the option of adjustingparameters to reproduce experimental data or computa-tionally very demanding many-body electronic structuremethods.The model includes a self-consistent Hartree potentialwhich takes into account the effect of an external bias aswell as continuum dielectric regions and external electro-static gates.We used the model to study a Tour wire between goldelectrodes, and found that the voltage drop in the devicecompares well with ab initio results, while the calculatedcurrent–voltage characteristics qualitatively agree betterwith experimental findings than the corresponding DFT-LDA results do.We also considered a graphene nano-transistor, andour study illustrated how the transport mechanismchanges from tunnelling to thermionic emission as thedevice is made longer.These applications show that the new method can givean accurate description of a broad range of nano-scaledevices. With its favorable computational speed, it is agood complement to ab initio -based transport methods. ACKNOWLEDGMENTS
This work was supported by the Danish Council forStrategic Research ’NABIIT’ under Grant No. 2106-04-0017, “Parallel Algorithms for Computational Nano-Science”, and European Commission STREP project No.MODECOM “NMP-CT-2006-016434”, EU. ∗ H. Haug and A.-P. Jauho,
Quantum Kinetics in Trans-port and Optics of Semiconductors (Berlin Springer Verlag,1996). N. D. Lang, Phys. Rev. B, , 5335 (1995). Y. Xue, Chemical Physics, , 151 (2002). M. Brandbyge, J.-L. Mozos, P. Ordej´on, J. Taylor, andK. Stokbro, Phys. Rev. B, , 165401 (2002). J. Taylor, H. Guo, and J. Wang, Phys. Rev. B, , 245407 (2001). A. Di Carlo, Physica B, , 211 (2002). A. Pecchia and A. Di Carlo, Reports in Prog. in Phys., ,1497 (2004). M. Magoga and C. Joachim, Phys. Rev. B, , 4722 (1997). S. Corbel, J. Cerda, and P. Sautet, Phys. Rev. B, , 1989(1999). J. Cerd´a and F. Soria, Phys. Rev. B, , 7965 (2000). E. G. Emberly and G. Kirczenow, Phys. Rev. B, , 10451(2001). F. Zahid, M. Paulsson, E. Polizzi, A. W. Ghosh, L. Sid-diqui, and S. Datta, J. of Chem. Phys., , 064707(2005). D. Kienle, J. I. Cerda, and A. W. Ghosh, J. Appl. Phys., , 043714 (2006). D. Kienle, K. H. Bevan, G.-C. Liang, L. Siddiqui, J. I.Cerda, and A. W. Ghosh, J. Appl. Phys., , 043715(2006). M. H. Whangbo and R. Hoffmann, J. Chem. Phys., ,5498 (1978). J. A. Pople and G. A. Segal, J. Chem. Phys, , 3289(1966). J. N. Murell and A. J. Harget,
Semi-empirical SCF Theoryof Molecules (Wiley, 1972). T. N. Todorov, J. Hoekstra, and A. P. Sutton, Philos.Mag. B, , 421 (2000). S. Datta,
Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, UK, 1997). J. G. Kushmerick, D. B. Holt, J. C. Yang, J. Naciri, M. H.Moore, and R. Shashidhar, Phys. Rev. Lett., , 86802(2002). J. Taylor, M. Brandbyge, and K. Stokbro, Phys. Rev.Lett., , 138301 (2002). The transport calculations where per-formed with Atomistix ToolKit, version2008.10. The manual is available online at . J. H. Ammeter, H. B. Burgi, J. C. Thibeault, and R. Hoff-mann, J. Am. Chem. Soc., , 3686 (1978). K. S. Thygesen and K. W. Jacobsen, Phys. Rev. B, ,033401 (2005). Q. Yan, B. Huang, J. Yu, F. Zheng, J. Zang, J. Wu, B.-L.Gu, F. Liu, and W. Duan, Nano Letters, , 1489 (2007). H. H. B. Sørensen, P. C. Hansen, D. E. Petersen, S. Skel-boe, and K. Stokbro, Phys. Rev. B,79