An efficient new method is presented to calculate the quantum transports
using periodic boundary conditions. This method allows the use of conventional
ground state ab initio programs without big changes. The computational effort
is only a few times of a normal ground state calculation, thus it makes
accurate quantum transport calculations for large systems possible.Comment: 9 pages, 6 figure

K. Hirose

Lin-Wang Wang

English

arXiv

An efficient new method is presented to calculate the quantum transports using periodic boundary conditions. This method allows the use of conventional ground state ab initio programs without big changes. The computational effort is only a few times of a normal groundstate calculations, thus is makes accurate quantum transport calculations for large systems possible

Wang, Lin-Wang

UNT Digital Library

Quantum Transport Calulations Using Periodi Boundary ConditionsLin-Wang WangComputational Researh Division, Lawrene Berkeley National Laboratory, Berkeley, CA 94720(Dated: June 14, 2004)An eÆient new method is presented to alulate the quantum transports using periodi boundaryonditions. This method allows the use of onventional ground state ab initio programs without bighanges. The omputational eort is only a few times of a normal ground state alulation, thus itmakes aurate quantum transport alulations for large systems possible.PACS numbers: 71.15.-m, 73.63.-b, 73.22.-fQuantum transport for moleules, nanowires, and nan-odevies is a fast growing researh area in both experi-ment and theory, with the potential of replaing the ur-rent Si based tehnology after the Moor's law reahes itslimit in about 15 years. In the theoretial ballisti trans-port alulations, a key step is to alulate the urrentvia the Landauer formula:I = 2eh Z RL Xn Tn(E)dE; (1)where L and R are left and right eletrode Fermi en-ergies (assuming the urrent ows from right to left in zdiretion), and Tn(E) is the transmission oeÆient forthe nth right hand eletrode hannel (band) at energyE. There are two major ways to alulate Tn(E). Oneis to use the Green's funtion G(r; r0; E) of the system.However sine G is a double variable funtion, ompu-tationally this approah an be quite expensive, thus itis mostly used for loalized basis set methods [1℄. Theother way to alulate Tn(E) is to solve the followingsattering states: H s(r) = E s(r) (2)and for z !1( 1): s(r) =Xn [AR(L)n R(L)n (r) +BR(L)n R(L)n (r)℄; (3)with onditions: ARn = 0 exept ARm = 1 for one m,and BLn = 0 (assuming dER(L)n (k)=dk > 0). In above,H is the single partile Hamiltonian, and R(L)n (r) =un;kn(r)exp(ikR(L)n z) are right going running waves inthe the right(R) and left(L) eletrodes, and R(L)nare the left going running waves. ER(L)n (kR(L)n ) =E are the eletrode band struture. The transmis-sion oeÆient for hannel m an be alulated asTm(E) = [Pn jALn j2(dELn (k)=dk)℄=(dERm(k)=dk), and thereetion oeÆient an be alulated as Rm(E) =[Pn jBRn j2(dERn (k)=dk)℄=(dERm(k)=dk). Transfer matrixmethod [2, 3℄ and the Lippmann-Shwinger equation [5℄have been used to solve Eqs(2),(3). Unfortunately, thetransfer matrix method is rather ompliated and om-putationally expensive to deal with the nonloal pseu-dopotentials [3℄ and it is often plagued by the numerialinstability due to the evanesent states in a multi-hanneleletrode [4℄. On the other hand, the use of Lippmann-Shwinger equation [5℄ requires the solution of a linearequation of the dimension of the full system, and it alsoneeds the Green's funtion of the two eletrode systemunder a potential bias. As a result, urrently this ap-proah is only used for jellium eletrode model and rela-tively small systems. Overall, ompared to the more ma-tured ground state alulations, all the urrent methodsfor transport alulations are ompliated and omputa-tionally expensive, and they an only be used to alulaterelatively small systems although there is a strong needto study the transports of large moleules and nanostru-tures. Here, we present a new and simple approah whihmakes the transport alulation similar to the groundstate alulation. In this approah, onventional periodisuperell methods and a speially designed perturbativeapproah are used to solve Eqs(2),(3). This allows usto use modern ab initio total energy programs withoutmuh hange. The omputational eort is similar to anormal ground state loal density approximation (LDA)alulation, hene it opens the door for quantum trans-port studies for large systems. W   (r)mV  (r)vLeft electrode Right electrodemoleculeBBzΩ ΩFIG. 1: A shemati view of the alulated system.To demonstrate our method, we have hosen a benzenemoleule onneted by two Cu quantum wires as shownin Fig.1. This system is hosen sine the ondutivityof the benzene moleule is well studied [5, 6℄ and similar2quantum wires have been used as eletrodes in previousquantum transport alulations [7℄. Two hydrogen atomsat the two ends of the benzene moleule are replaedby two sulfur atoms, whih are bonded to two entralCu atoms at the eletrode. The atomi positions of themoleule are relaxed at the zero voltage bias under LDAalulation. We have used norm onserving pseudopo-tentials and 30 Ryd planewave uto with a standardplanewave LDA program [8℄. We have inluded 5 and 6unit ells in the left and right eletrodes respetively, andthey are onneted at boundary B in Fig.1 by periodiboundary ondition. The x, y dimensions of the super-ell are 3 times the width of the Cu wire to avoid possibleneighbore-neighbore interations. After the Kohn-Shamsingle partile potential V0(r) is obtained from a LDAselfonsistent alulation at the zero bias, we have addeda potential V=2sin(z=L0) in the entral region of themoleule and shifted the rest of the right (left) eletrodeby V=2 ( V=2) to get the potential VV (r) for a bias V sys-tem. Although a selfonsistent treatment an be ahievedstraight forwardly under the urrent approah [sine thesattering states of Eq(2) will be alulated℄, the urrentnonselfonsistent treatment for nite bias V is suÆientin illustrating the new methodology. Note that, there isa jump of VV (r) at the boundary B, but that is not aproblem in our numerial alulations.E(eV)n0.40.20 0.6 0.8 1−4−6−220462134567,891011Γ X’kz (   /a)πFIG. 2: The band struture of the eletrode. Eah ontinuousline from   to X 0 is denoted as one band. The zero is theeletrode Fermi energy. The rosses are the kRn points, seetext for details.Figure 2 shows the band strutures En(kz) of the quan-tum wire eletrode. There are 9 Cu atoms in eah unitell of the eletrode. To simplify our alulation andanalysis, we did not inlude the 3d eletrons in our pseu-dopotential. Although this will introdue a signianterror in the eletrode total energy, the eletroni stru-ture near the Fermi surfae and the related transportproperties are intat. However, for bias larger than 2 V,our eletrode should be onsidered as a model eletrodedue to the lak of Cu 3d eletrons.To solve the sattering states of Eqs(2),(3), we rstalulate the eigenstates f i(r), Eig of the periodi su-perell under VV (r) with a Kz point (e.g, Kz = =2Lz,where Lz is the length of the superell and this Kzis not the kz of the eletrode as in Fig.2) using ourstandard LDA program [8℄. Let's rst assume that,using some methods, we an generate l degeneratedstates  i;(l)(r) they all satisfy the Shrodinger's equa-tion H i;(l) = Ei i;(l) (however, for our purpose, thisequation needs only to be satised within the interior ofthe superell, not near the boundary B of Fig.1). Theidea is to use a linear ombination of these states toonstrut the sattering states of Eq(2). First, withinthe R(L) eletrode,  i;(l)(r) an be deomposed into theeletrode states R(L)n (r) just as in Eq(3). From a givenEi, the available n and kR(L)n an be found from Fig.2,or say: En(kR(L)n ) + R(L) = Ei. The orrespondingR(L)n is then generated by numerial interpolations frompre-alulated eletrode states. The expansion oeÆ-ients AR(L)n (i; l), BR(L)n (i; l) for wavefuntion  i;(l) anbe easily alulated from the funtional produts like:R  i;(l)nd3r and R mnd3r, where  is one eletrodeunit ell at the middle of the eletrodes as shown in Fig.1.We found that, after inluding the possible evanesentstates [9℄, this expansion typially aptures more than99:99% of the weight of the original  i;(l). The next stepis to make a linear ombination of  i;(l) to get the sat-tering state  i;s of Eq(2): i;s = Pl Cl i;(l) (4)= PnPl[ClAR(L)n (i; l)n(r) + ClBR(L)n (i; l)n(r)℄;here the seond equation and the R and L are for r withinthe right and left eletrodes respetively. Note that, dueto the use of superell Kz point and both  i and  i areused as  i;(l),  i;s is no longer periodi at the superellboundary B. To make  i;s in Eq(4) as one satteringstate of Eq(3), we need to make it satisfy the boundaryonditions of Eq(3):ARn6=m = 0, BLn = 0, by seleting Cl.In order to have a solution for Cl, we need N indepen-dent  i;(l) (l = 1; N), if there are N nonzero ontributingeletrode states n in Eq(4) [ounting both the left andright eletrodes, but n and n are ounted as one℄.Notie that, sine both  i(r) and  i (r) satisfy theShrodinger's equation H = Ei , for systems withonly a single hannel (N=2), the eigen state  i(r) itselfis enough to onstrut the sattering states  i;s fromEq(4). So, the main task here is for multi-hannel ases.In our example, from Fig.2 we see that for a given energyEi, we ould have 4-5 hannels. To get the N degener-ated  i;(l), we will use a perturbative approah. First,we will onstrut Wm(r) = um; (r) (the k = 0 mth ele-trode state) when z is within the last unit ell of the righteletrode near boundary B of Fig.1, and Wm(r) = 0 forall the other z (see Fig.1). We will add jWm >< Wmjas a perturbation in the original H , and solve the follow-ing eigenstate equation using our standard LDA program(e.g, using onjugate gradient method):H 0i;m +  < Wmj 0i;m > Wm = (Ei +Ei;m) 0i;m: (5)3Here  is a very small number, hene Ei;m and  i;m  0i;m    i are both small. Suppose we have solved theabove equations for two dierent m's: m1 and m2. Thenwe an onstrut  i;(l) = F1 i;m1 + F2 i;m2 , withF1Ei;m1 + F2Ei;m2 = 0. After dropping the seondorder terms Ei;m1;2 i;m1;2 we have:H i;(l) +  Xj=1;2Fj < Wmj j 0i;mj > Wmj = Ei i;(l) (6)Notie that the Wmj terms are nonzero only nearthe boundary B, so for all the other plaes, we haveH i;(l) = Ei i;(l). Thus,  i;(l) are the wavefuntionswe needed. In the simple ases, when there are N totaleletrode states R(L)n with nonzero omponents in theexpansion of  i(r) in Eq(4), there will be N/2 right ele-trode states Rn . Then the perturbations by the relatedN/2 Wm states (whih have the same harateristis andross setion symmetries as Rn ) will introdue N/2 inde-pedent perturbative wavefuntion hanges  i;m. These i;m will generate (N=2   1) independent  i;(l) states(besides the original  i). Thus, the total number of  i;(l)states (ounting also  i;(l)) is just N, the exat number weneed to onstrut the sattering state  i;s from Eq(4).This argument remains true when there are evanesentstates or the number of eletrode states in the left andright eletrodes are not the same. Thus, using this pro-edure, we are guaranteed that there will be enough  i;(l)states for a given  i to generate a few orresponding sat-tering states  i;s.From the superell eigenstates f i(r); Eig, we an gen-erate a set of fkRn g from En(kRn )+R = Ei. These fkRn gare shown in Fig.2 as the rosses for a 1V bias ase (usingall the  i(r) with Ei between the two horizontal arrowsin Fig.2). As an be desribed by a phase aumulationmodel[10℄, on eah band, these kRn have roughly equal dis-tanes and their total number roughly equals the numberof eletrode unit ells. We have typially used 6   pointeletrode states as Wm in Eq(5) starting from the low-est band as annotated in Fig.2. This means we have tosolve f 0i;mg of Eq(5) 6 times using f ig as the initialwavefuntions (Notie that, this number 6 is roughly thenumber of hannels in the problem. The same prefatoris needed in the alulations of other methods like thetransfer matrix or Lippmann-Shwinger equation). Af-ter f 0i;mg are alulated, using Eq(4), we an onstruta sattering state from eah of these kRn shown in Fig.2.Two of these onstruted sattering states are shown inFig.3. Notie that the dash lines arePl Cl i;(l) of Eq(4),while the solid lines are the eletrode state deomposi-tions [the seond line of the Eq(4)℄. Within the eletrode,the eletrode state deomposition gives a very auratedesription of the total wavefuntion. From these sat-tering states, the transmission oeÆients Tm(Ei) an bealulated, and are shown in Fig.4 as the symbols. Thealulated Tm(Ei) +Rm(Ei) is typially very lose to 1,indiating the numerial stability of the urrent method.100.110.010.010.1101Z direction positionPlanar averaged wavefunction square(arbitrary unit)T=0.36T=0.008E=−3.264 eVE= −1.043 eVRight electrodeband=5band=5,5’Left electrodeRight electrodeband=3Left electrodeband=3FIG. 3: The onstruted sattering states from Eq(4). Thedashed and solid lines orrespond to the rst and seond linesof Eq(4) respetively. T is the transmission oeÆient, E isthe eigen energy, and band number indiate the n of n inEq(4). The bias of the system is 1V.ln(T  )n0 0.4 0.6 0.8 10.2−16−12−8−40n=2n=5n=3n=4kz (   /a)πFIG. 4: The alulated transmission oeÆients Tn(kRn ) (sym-bols) and tted smooth urves fn(k) (lines) for dierent bandsof the eletrode. The bias of the system is 1V.Notie that, unlike the other approahes disussed be-fore where the sattering states of arbitrary energy Ean be solved, here only the sattering states of energyfEig are alulated. This translates into nite numberof fkRn g points as shown in Fig.2 and Fig.4. In Eq(1),we need all the energies between L and R. This isan omplete analogy with the k-point integration prob-lem in onventional ground state bulk alulation [theenergy integral in Eq(1) an also be hanged into a k-point integral of the eletrode band struture of Fig.2℄.Thus, similar to the onventional bulk alulations, here4we will use an interpolation sheme to arry out the inte-gral in Eq(1). First, if the number of fkRn g is not enoughin Fig.4, we an hoose a dierent superellKz in our su-perell alulations (or hange the potential VV (r) nearboundary B), and repeat the above proedure. That willgive us more fEig and fkRn g points. In our system, wend one Kz alulation is suÆient. We have used asmooth urve fn(k) to interpolate the ln(Tn(kRn )) pointsshown in Fig.4. More speially, we have minimized:PkRn jln(Tn(kRn )) fn(kRn )j2+ R jd2fn(k)=dk2j2dk, with  1. Numerially, this orresponds to a simple linearequation with disretized k points. The resulting urve isshown in Fig.4 for the 1V bias ase. Using these fn(k), wean alulate the total transmission T (E) = Pn Tn(E)of the system. The results are shown in Fig.5 for dier-ent biases. We see that T (E) is inuened strongly bytwo fators. One is the relative energy levels between theeletrode states and the moleule states. When the biasinreases, the moleular levels drop relative to the righteletrode state levels. As a result, the magnitude of thetransmission dereases near the region of -3 eV. Anotherfator is the band struture of the eletrode. There is awell shape of T (E) near -0.3 eV. This is aused by bandgaps of the 2,3 bands at the X 0 point in Fig.2. The T (E)also shows a big drop at -3.3 eV. This is due to the endof 2,3 bands at the   point.1V4V3V2V0VE −EF (right electrode) (eV)−2−3 −1 0 10.20.00.40.60.81.01.2−4Total transmission coefficient T(E) FIG. 5: The alulated total transmission oeÆients T (E)(whih an be larger than 1) for dierent biases. The zero isthe right eletrode Fermi energy. For a given bias V, there arenet right to left urrent ow only within the [ V; 0℄ energywindow.After T (E)'s of Fig.5 are obtained, a simple energyintegration between -V to 0 will give us the total ur-rent I . The resulting I(V ), and the ondutane dI=dVare shown in Fig.6. We do see the well known peak anddip for this system in the ondutane around 2V. Ouralulated peak and dip positions of 1.8V and 2.3V or-responds well with the experimental results of 1.4V and2.4V [6℄, and this agreement is better than previously al-ulated results [5℄. We also see the marks of the eletrodeeletroni strutures. Near 0.3V, the ondutane showsa well shape, again due to the band gaps of 2,3 bandsat the X 0 point. Above 3.3V, we see a big drop, thenthe negative ondutane. The drop is due to the end ofthe 2,3 bands at the   points, and the negative ondu-tane is beause the onduting eletrode levels (energywindow) are moving away from the onduting moleu-lar levels. Here we see that, the eletroni struture ofthe eletrode is extremely important in determining theoverall ondutane of the system.0.100.2−0.1 −0.100.10.20 31 2 4(G  V)0Current (G  ) 0Conductance Voltage (V)FIG. 6: The alulated I-V urve and the orresponding on-dutane. G0 = 2e=h = 77s.In summary, we have presented a simple and numeri-ally stable sheme to alulate the quantum transport.Within this sheme, the onventional periodi superellsand ground state ab initio programs an be used with-out muh hange. The omputational eort is only a fewtimes of a onventional ground state alulation. Thispromised quantum transport alulations for muh largersystems whih annot be takled by othe methods. Theimplementation of this method is simple and straight for-ward based on any onventional ground state ab initioprograms.This work was supported by U.S. Department of En-ergy under Contrat No. DE-AC03-76SF00098. Thisresearh used the resoures of the National Energy Re-searh Sienti Computing Center.[1℄ J. Taylor, H. Guo, J. Wang, Phys. Rev. B 63, 245407(2001); P. Damle, A. Ghosh, S. Datta, Phys. Rev. B 64,201403 (2001); M. B. Nardelli, J.-L. Fattebert, J. Bern-hol, Phys. Rev. B 64, 245423 (2001); M. Brandbyge,et.al, Phys. Rev. B 65, 165401 (2002); E. Louis, et.al,Phys. Rev. B 67, 155321 (2003). [2℄ K. Hiros, M. Tsukada, Phys. Rev. B 51, 5278 (1994);[3℄ H.J. Choi, J. Ihm, Phys. Rev. B 59, 2267 (1999).[4℄ D. Z.-Y. Ting, E. T. Yu, T.C. MGillha, Phys. Rev. B45, 3583 (1992).[5℄ M. Di Ventra, S.T. Pentelides, N. D. Lang, Phys. Rev.Lett. 84, 979 (2000); M. Di Ventra, N. D. Lang, Phys.5Rev. B 65, 45402 (2001).[6℄ M.A. Reed, et.al, Siene 278, 252 (1997).[7℄ C. Roland, V. Meunier, B. Larade, H. Guo, Phys. Rev.B 66, 35332 (2002).[8℄ http://rd.lbl.gov/ linwang/PEtot/PEtot.html[9℄ When the energy Ei is lose to an extremum of a bandenergy (e.g, at the   and X 0 points, or inside an anti-rossing band gap of Fig.2), an evanesent state mightexist at the middle eletrode ell  where the deompo-sition is alulated. These evanesent states ome fromthe omplex k-point band struture of the eletrode [see:Y.C. Chang, Phys. Rev. B 25, 605 (1982)℄. If the al-ulated eletrode is very long, these evanesent statesshould have deayed to zero in the middle of the ele-trode. But in pratie, and espeially when the energy islose to the extrema, these evanesent states an have sig-niant ontributions. We have used un(km)exp(ikmz) atthose extrema km to approximate the evanesent states,and inlude them as n(r) in the deomposition of i;(l).Note that, these evanesent states arry no urrents sinedEn(k)=dkjk=km = 0.[10℄ N.V. Smith, N.B. Brookes, Y. Chang, and P.D. Johnson,Phys. Rev. B 49, 332 (1994).

Quantum transport calculations using periodic boundaryconditions

Crossref

Elastic quantum transport calculations using auxiliary periodic boundary conditions

An efficient new method is presented to calculate the quantum transports using periodic  boundary conditions. This method allows the use of conventional ground state ab initio programs  without big changes. The computational effort is only a few times of a normal ground state  calculations, thus is makes accurate quantum transport calculations for large systems  possible

Quantum Transport Calculations Using Periodic Boundary Conditions

Abstract

Similar works

Full text

Available Versions

UNT Digital Library

Crossref

eScholarship - University of California