C ***************** INTCURR ****************** C C ROUTINES FOR COMPUTING TUNNEL CURRENT FOR EFFECTIVE MASS C BULK BANDS, USING BARDEEN FORMALISM AND T&H APPROXIMATION C C VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, 2007-2009 C 1.1 - DEC/10, ELIMINATE USAGE OF 2-BAND MODEL, AND MODIFY C OUTPUT OF CURRENT VALUES C C OUTPUT INTO LOGICAL UNIT NUMBER LUN+I: C I=0, QUANTUM NO. AND ENERGY OF LOCALIZED STATES, VS BIAS C I=1, WAVEFUNCTIONS OF LOCALIZED STATES (AT FINAL BIAS) C I=3, EXTENDED STATE WAVEFCN AT KPARR=0 (AT FINAL BIAS) C SUBROUTINE CURR2(ICUT,VAC,TIP,SEM,VSINT,NV,NR,NS,NVDIM,NRDIM, &NSDIM,S,SEP,DELV,BIAS,EF,CHI,EFTIP,DELPHI,AVBH,AVBL,AVBSO, &ESO,E2HH,E2SO,nee,nwk,PHI0,NVDIM1,NVDIM2,NSDIM1,NSDIM2, &EXPANI,FRAC,NLOC,CURRV,CURRV0,CURRC,CURRC0,CURR,CURR0,IWRIT) C DIMENSION VAC(2,NRDIM,NVDIM),VSINT(2,NRDIM),SEM(2,NRDIM,NSDIM), &S(NSDIM),DELV(NRDIM),BARR(NVDIM1),PROF(NSDIM),NLOC(4), &CDESEM(NSDIM),CDLSEM(NSDIM),CDEVAC(NVDIM1),CDLVAC(NVDIM1) LOGICAL TIP(NRDIM,NVDIM) real kappa,lambda DOUBLE PRECISION SUM1,SUM2,SUM1S,SUM2S,SUM1P,SUM2P COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV DATA RQUANT/12900./ PI=4.*ATAN(1.) tk1=tk tk2=tk CDESURF=0. CDLSURF=0. DO 65 J=1,NS CDESEM(J)=0. CDLSEM(J)=0. 65 CONTINUE DO 75 J=1,NV+1 CDEVAC(J)=0. CDLVAC(J)=0. 75 CONTINUE C C GET A CUT FROM THE POTENTIAL C CALL POTCUT1(ICUT,VAC,TIP,SEM,VSINT,NRDIM,NVDIM,NSDIM, &NV,NS,SEP,S,R,DELV,PHI0,BIAS,CHI,DELPHI,BARR,PROF,NBARR1, &NVDIM1,NVDIM2,PMIN,PMAX,0) c c valence band c EV=0. LUN=30 CALL VBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMAX,AVBL,TK1,TK2, &BIAS,EF,NEE,NWK,NL,EV,EFTIP,EGAP,LUN,CURRVL,CURRV0L,IWRIT) IF (IWRIT.NE.0) THEN write(6,*) 'number of VB light-hole localized states =',nl write(16,*) 'number of VB light-hole localized states =',nl END IF NLOC(1)=NL IF (IWRIT.GE.2) CLOSE(LUN+1) IF (IWRIT.GE.2) CLOSE(LUN+2) IF (IWRIT.GE.4) CLOSE(LUN+3) LUN=40 CALL VBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMAX,AVBH,TK1,TK2, &BIAS,EF,NEE,NWK,NL,EV,EFTIP,E2HH,LUN,CURRVH,CURRV0H,IWRIT) IF (IWRIT.NE.0) THEN write(6,*) 'number of VB heavy-hole localized states =',nl write(16,*) 'number of VB heavy-hole localized states =',nl END IF NLOC(2)=NL IF (IWRIT.GE.2) CLOSE(LUN+1) IF (IWRIT.GE.2) CLOSE(LUN+2) IF (IWRIT.GE.4) CLOSE(LUN+3) EVSO=EV-ESO LUN=50 CALL VBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMAX,AVBSO,TK1,TK2, &BIAS,EF,NEE,NWK,NL,EVSO,EFTIP,E2SO,LUN,CURRVSO,CURRV0SO,IWRIT) IF (IWRIT.NE.0) THEN write(6,*) 'number of VB split-off localized states =',nl write(16,*) 'number of VB split-off localized states =',nl END IF NLOC(3)=NL IF (IWRIT.GE.2) CLOSE(LUN+1) IF (IWRIT.GE.2) CLOSE(LUN+2) IF (IWRIT.GE.4) CLOSE(LUN+3) currv=currvL+currvH+currvSO currv0=currv0L+currv0H+currv0SO IF (IWRIT.NE.0) THEN IF ((IINV.EQ.1.OR.IINV.EQ.3).AND.CURRV.GT.0.) THEN CURRV=0. WRITE(6,*) 'VB EXTENDED INVERSION CURRENT SET TO ZERO' WRITE(16,*) 'VB EXTENDED INVERSION CURRENT SET TO ZERO' END IF IF ((IINV.EQ.1.OR.IINV.EQ.3).AND.CURRV0.GT.0.) THEN CURRV0=0. WRITE(6,*) 'VB LOCALIZED INVERSION CURRENT SET TO ZERO' WRITE(16,*) 'VB LOCALIZED INVERSION CURRENT SET TO ZERO' END IF IF (NWK.NE.1) THEN write(6,*) 'valence band current ext,loc =',currv,currv0 write(16,*) 'valence band current ext,loc =',currv,currv0 END IF END IF c c conduction band c EC=EGAP LUN=60 CALL CBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMIN,ACB,TK1,TK2, &BIAS,EF,NEE,NWK,NL,EC,EFTIP,EGAP,LUN,CURRC,CURRC0,IWRIT) IF (IWRIT.NE.0) THEN write(6,*) 'number of CB localized states =',nl write(16,*) 'number of CB localized states =',nl END IF NLOC(4)=NL IF (IWRIT.GE.2) CLOSE(LUN+1) IF (IWRIT.GE.2) CLOSE(LUN+2) IF (IWRIT.GE.4) CLOSE(LUN+3) IF (IWRIT.NE.0) THEN IF ((IINV.EQ.2.OR.IINV.EQ.3).AND.CURRC.LT.0.) THEN CURRC=0. WRITE(6,*) 'CB EXTENDED INVERSION CURRENT SET TO ZERO' WRITE(16,*) 'CB EXTENDED INVERSION CURRENT SET TO ZERO' END IF IF ((IINV.EQ.2.OR.IINV.EQ.3).AND.CURRC0.LT.0.) THEN CURRC0=0. WRITE(6,*) 'CB LOCALIZED INVERSION CURRENT SET TO ZERO' WRITE(16,*) 'CB LOCALIZED INVERSION CURRENT SET TO ZERO' END IF IF (NWK.NE.1) THEN write(6,*) 'conduction band current ext,loc =',currc,currc0 write(16,*) 'conduction band current ext,loc =',currc,currc0 END IF END IF curr=currv+currc curr0=currv0+currc0 return end C C VB TUNNEL CURRENT C SUBROUTINE VBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMAX,EFFM,TK1,TK2, &BIAS,EF,NE,NWK,NLOC,EV,EFTIP,E2BAND,LUN,CURRV,CURRV0,IWRIT) C REAL KAPPA DIMENSION BARR(NVDIM1),PROF(NSDIM),S(NSDIM),BARR2(NVDIM2), &PROF2(NSDIM2),S2(NSDIM2),JSEM(NSDIM2),NEXSEM(NSDIM1) DOUBLE PRECISION SUM,SUM2,SUM3,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/RQUANT/12900./ PI=4.*ATAN(1.) NLOC=0 CURRV=0. CURRV0=0. WKFTIP=sqrt(C*EFTIP) C C EXTENDED STATES C emax=EV IF (NWK.NE.1) THEN emin=amin1(ef-10.*tk1,ef+bias-10.*tk2) ELSE emin=ef-10.*tk1 END IF dele=(emax-emin)/ne sum=0. if (dele.le.0.) go to 200 wkmax=sqrt(C*EFFM*(emax-emin)) SEMSTEP=(2.*pi/wkmax)/EXPANI KAPPA=SQRT(C*(amax1(BARR(NBARR1),BARR(1))-EMIN)) VACSTEP=(2.*pi/KAPPA)/EXPANI CALL POTCUT2(SEP,NV,PHI0,S,NS,NSDIM,BARR,NBARR1,BARR2, &NBARR2,NVDIM1,NVDIM2,PROF,PROF2,NSDIM2,S2,NS2,VACSTEP,SEMSTEP, &FRAC,JSEM,NEXSEM,NEXVAC,0,NSDIM1) DELVAC=sep/float(NBARR2-1) delwk=wkmax/nwk do 120 iwky=0,nwk-1 wky=iwky*delwk do 115 iwkx=0,nwk-1 wkx=iwkx*delwk wkparr=sqrt(wkx**2+wky**2) eparr=wkparr**2/(EFFM*C) nwkdeg=8 if (iwkx.eq.0) nwkdeg=nwkdeg/2 if (iwky.eq.0) nwkdeg=nwkdeg/2 if (iwkx.eq.iwky) nwkdeg=nwkdeg/2 if (iwky.gt.iwkx) go to 115 do 110 ie=1,ne ener=emax-(ie-0.5)*dele if (eparr.ge.(emax-ener)) go to 110 occtip=fd(ener-bias,ef,tk2) occsem=fd(ener,ef,tk1) occdiff=occtip-occsem call VBwf(wf,wfderiv,WKSEM,ener,wkparr,sep,bias,BARR2, & NVDIM2,NBARR2,PROF2,NS2,NSDIM2,S2,EFFM,EV,E2BAND,PSIVAC, & PSISEM) IF (IWKX.EQ.0.AND.IWKY.EQ.0) THEN C C WRITE OUTPUT IF DESIRED C IF (IWRIT.GE.4) THEN do 50 J=NBARR2,1,-1 write(LUN+3,*) -(J-1)*DELVAC,PSIVAC(J) 50 CONTINUE do 60 J=1,NS2 write(LUN+3,*) S2(J),PSISEM(J) 60 CONTINUE END IF END IF c eperp, kappa in vacuum EPERP=ener-WKPARR**2/C KAPPA=SQRT(C*(BARR2(NBARR2)-EPERP)) trans=2.*nwkdeg*(2.*wf)**2*WKFTIP/(WKSEM/EFFM) sum=sum+trans*occdiff 110 continue 115 continue 120 continue currv=sum*dele*delwk**2/(4.*PI**2*RQUANT) C C LOCALIZED STATES C 200 IF (PMAX.LE.0.) RETURN emax=pmax+EV c could delete statements below! IF (NWK.NE.1) THEN emin=amax1(EV,amin1(ef-10.*tk1,ef+bias-10.*tk2)) ELSE emin=amax1(EV,ef-10.*tk1) END IF dele=(emax-emin)/ne emin2=ef-10.*tk1 dele2=(emax-emin2)/ne sum=0. if (dele.le.0.) go to 460 wkmax=sqrt(C*EFFM*(emax-emin)) SEMSTEP=(2.*pi/wkmax)/EXPANI KAPPA=SQRT(C*(amax1(BARR(NBARR1),BARR(1))-EMIN)) VACSTEP=(2.*pi/KAPPA)/EXPANI CALL POTCUT2(SEP,NV,PHI0,S,NS,NSDIM,BARR,NBARR1,BARR2, &NBARR2,NVDIM1,NVDIM2,PROF,PROF2,NSDIM2,S2,NS2,VACSTEP,SEMSTEP, &FRAC,JSEM,NEXSEM,NEXVAC,0,NSDIM1) DELVAC=sep/float(NBARR2-1) delwk=wkmax/nwk do 320 iwky=0,nwk-1 wky=iwky*delwk do 315 iwkx=0,nwk-1 wkx=iwkx*delwk wkparr=sqrt(wkx**2+wky**2) eparr=wkparr**2/(EFFM*C) n=0 nsav=0 nwkdeg=8 if (iwkx.eq.0) nwkdeg=nwkdeg/2 if (iwky.eq.0) nwkdeg=nwkdeg/2 if (iwkx.eq.iwky) nwkdeg=nwkdeg/2 if (iwky.gt.iwkx) go to 315 do 310 ie=1,ne ener=emax-(ie-0.5)*dele if (eparr.ge.(emax-ener)) go to 310 occtip=fd(ener-bias,ef,tk2) occsem=fd(ener,ef,tk1) occdiff=occtip-occsem call VBloc(n,wf,wfderiv,ener,wkparr,sep,bias,BARR2,NVDIM2, & NBARR2,PROF2,NS2,NSDIM2,S2,EFFM,EV,E2BAND,PSIVAC,PSISEM) if (n.eq.nsav) go to 310 IF (IWKX.EQ.0.AND.IWKY.EQ.0) THEN IF (PSISEM(1).NE.0) THEN NLOC=NLOC+1 C C FOUND LOCALIZED STATE; WRITE OUTPUT IF DESIRED C IF (IWRIT.GT.1) THEN write(6,*) 'VB localized state at energy ',ener write(16,*)'VB localized state at energy ',ener END IF IF (IWRIT.GE.1) THEN WRITE(LUN,*) BIAS,N,ENER,OCCDIFF END IF IF (IWRIT.GE.2) THEN do 250 J=NBARR2,1,-1 write(LUN+1,*) -(J-1)*DELVAC,PSIVAC(J) 250 CONTINUE do 260 J=1,NS2 write(LUN+1,*) S2(J),PSISEM(J) 260 CONTINUE END IF END IF END IF c eperp, kappa in vacuum EPERP=ener-WKPARR**2/C KAPPA=SQRT(C*(BARR2(NBARR2)-EPERP)) trans=nwkdeg*(n-nsav)*2.*(2.*wf)**2*WKFTIP sum=sum+trans*occdiff nsav=n 310 continue 315 continue 320 continue currv0=sum*delwk**2/(C*2.*PI*RQUANT) 460 RETURN END C C CB TUNNEL CURRENT C SUBROUTINE CBCURR1(BARR,NBARR1,S,PROF,NS,NSDIM,NVDIM,NV,NVDIM1, &NSDIM1,NSDIM2,NVDIM2,SEP,PHI0,EXPANI,FRAC,PMIN,EFFM,TK1,TK2, &BIAS,EF,NE,NWK,NLOC,EC,EFTIP,E2BAND,LUN,CURRC,CURRC0,IWRIT) C REAL KAPPA DIMENSION BARR(NVDIM1),PROF(NSDIM),S(NSDIM),BARR2(NVDIM2), &PROF2(NSDIM2),S2(NSDIM2),JSEM(NSDIM2),NEXSEM(NSDIM1) DOUBLE PRECISION SUM,SUM2,SUM3,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/RQUANT/12900./ PI=4.*ATAN(1.) NLOC=0 CURRC=0. CURRC0=0. WKFTIP=sqrt(C*EFTIP) C C EXTENDED STATES C emin=EC IF (NWK.NE.1) THEN emax=amax1(ef+10.*tk1,ef+bias+10.*tk2) ELSE emax=ef+10.*tk1 END IF dele=(emax-emin)/ne sum=0. if (dele.le.0.) go to 200 wkmax=sqrt(C*EFFM*(emax-emin)) SEMSTEP=(2.*pi/wkmax)/EXPANI KAPPA=SQRT(C*(amax1(BARR(NBARR1),BARR(1))-EMIN)) VACSTEP=(2.*pi/KAPPA)/EXPANI CALL POTCUT2(SEP,NV,PHI0,S,NS,NSDIM,BARR,NBARR1,BARR2, &NBARR2,NVDIM1,NVDIM2,PROF,PROF2,NSDIM2,S2,NS2,VACSTEP,SEMSTEP, &FRAC,JSEM,NEXSEM,NEXVAC,0,NSDIM1) DELVAC=sep/float(NBARR2-1) delwk=wkmax/nwk do 120 iwky=0,nwk-1 wky=iwky*delwk do 115 iwkx=0,nwk-1 wkx=iwkx*delwk wkparr=sqrt(wkx**2+wky**2) eparr=wkparr**2/(EFFM*C) nwkdeg=8 if (iwkx.eq.0) nwkdeg=nwkdeg/2 if (iwky.eq.0) nwkdeg=nwkdeg/2 if (iwkx.eq.iwky) nwkdeg=nwkdeg/2 if (iwky.gt.iwkx) go to 115 do 110 ie=1,ne ener=emin+(ie-0.5)*dele if (eparr.ge.(ener-emin)) go to 110 occtip=fd(ener-bias,ef,tk2) occsem=fd(ener,ef,tk1) occdiff=occtip-occsem call CBwf(wf,wfderiv,WKSEM,ener,wkparr,sep,bias,BARR2, & NVDIM2,NBARR2,PROF2,NS2,NSDIM2,S2,EFFM,EC,E2BAND,PSIVAC, & PSISEM) IF (IWKX.EQ.0.AND.IWKY.EQ.0) THEN C C WRITE OUTPUT IF DESIRED C IF (IWRIT.GE.4) THEN do 50 J=NBARR2,1,-1 write(LUN+3,*) -(J-1)*DELVAC,PSIVAC(J) 50 CONTINUE do 60 J=1,NS2 write(LUN+3,*) S2(J),PSISEM(J) 60 CONTINUE END IF END IF c eperp, kappa in vacuum EPERP=ener-WKPARR**2/C KAPPA=SQRT(C*(BARR2(NBARR2)-EPERP)) trans=2.*nwkdeg*(2.*wf)**2*WKFTIP/(WKSEM/EFFM) sum=sum+trans*occdiff 110 continue 115 continue 120 continue currc=sum*dele*delwk**2/(4.*PI**2*RQUANT) C C LOCALIZED STATES C 200 IF (PMIN.GE.0.) RETURN emin=pmin+EC c could delete statements below! IF (NWK.NE.1) THEN emax=amin1(EC,amax1(ef+10.*tk1,ef+bias+10.*tk2)) ELSE emax=amin1(EC,ef+10.*tk1) END IF dele=(emax-emin)/ne emax2=ef+10.*tk1 dele2=(emax2-emin)/ne sum=0. if (dele.le.0.) go to 460 wkmax=sqrt(C*EFFM*(emax-emin)) SEMSTEP=(2.*pi/wkmax)/EXPANI KAPPA=SQRT(C*(amax1(BARR(NBARR1),BARR(1))-EMIN)) VACSTEP=(2.*pi/KAPPA)/EXPANI CALL POTCUT2(SEP,NV,PHI0,S,NS,NSDIM,BARR,NBARR1,BARR2, &NBARR2,NVDIM1,NVDIM2,PROF,PROF2,NSDIM2,S2,NS2,VACSTEP,SEMSTEP, &FRAC,JSEM,NEXSEM,NEXVAC,0,NSDIM1) DELVAC=sep/float(NBARR2-1) delwk=wkmax/nwk do 320 iwky=0,nwk-1 wky=iwky*delwk do 315 iwkx=0,nwk-1 wkx=iwkx*delwk wkparr=sqrt(wkx**2+wky**2) eparr=wkparr**2/(EFFM*C) n=0 nsav=0 nwkdeg=8 if (iwkx.eq.0) nwkdeg=nwkdeg/2 if (iwky.eq.0) nwkdeg=nwkdeg/2 if (iwkx.eq.iwky) nwkdeg=nwkdeg/2 if (iwky.gt.iwkx) go to 315 do 310 ie=1,ne ener=emin+(ie-0.5)*dele if (eparr.ge.(ener-emin)) go to 310 occtip=fd(ener-bias,ef,tk2) occsem=fd(ener,ef,tk1) occdiff=occtip-occsem call CBloc(n,wf,wfderiv,ener,wkparr,sep,bias,BARR2,NVDIM2, & NBARR2,PROF2,NS2,NSDIM2,S2,EFFM,EC,E2BAND,PSIVAC,PSISEM) if (n.eq.nsav) go to 310 IF (IWKX.EQ.0.AND.IWKY.EQ.0) THEN IF (PSISEM(1).NE.0) THEN NLOC=NLOC+1 C C FOUND LOCALIZED STATE; WRITE OUTPUT IF DESIRED C IF (IWRIT.GT.1) THEN write(6,*) 'CB localized state at energy ',ener write(16,*)'CB localized state at energy ',ener END IF IF (IWRIT.GE.1) THEN WRITE(LUN,*) BIAS,N,ENER,OCCDIFF END IF IF (IWRIT.GE.2) THEN do 250 J=NBARR2,1,-1 write(LUN+1,*) -(J-1)*DELVAC,PSIVAC(J) 250 CONTINUE do 260 J=1,NS2 write(LUN+1,*) S2(J),PSISEM(J) 260 CONTINUE END IF END IF END IF c eperp, kappa in vacuum EPERP=ener-WKPARR**2/C KAPPA=SQRT(C*(BARR2(NBARR2)-EPERP)) trans=nwkdeg*(n-nsav)*2.*(2.*wf)**2*WKFTIP sum=sum+trans*occdiff nsav=n 310 continue 315 continue 320 continue currc0=sum*delwk**2/(C*2.*PI*RQUANT) 460 RETURN END C C C real function fd(e,ef,tk) if (tk.eq.0.) go to 100 if ((e-ef)/tk.gt.-40.and.(e-ef)/tk.lt.40) go to 200 if ((e-ef)/tk.gt.40.) fd=0. if ((e-ef)/tk.lt.-40.) fd=1. return 100 fd=0. if (e.eq.ef) fd=0.5 if (e.lt.ef) fd=1. return 200 fd=1./(1.+exp((e-ef)/tk)) end c c integrate VB wavefunction from tip across to sample, return c values of wf and derivative of wf at tip surface c subroutine VBwf(wf,wfderiv,WKSEM,e,wkparr,sep,bias,BARR2,NVDIM2, &NBARR2,PROF2,NS2,NSDIM2,S2,effm,EV,E2BAND,PSIVAC,PSISEM) dimension PROF2(NSDIM2),BARR2(NVDIM2),S2(NSDIM2) DOUBLE PRECISION PSI,DPSI,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/ PI=4.*ATAN(1.) c wf=0. wfderiv=0. WKSEM=1. eperp=e+wkparr**2/(C*effm) if (eperp.ge.EV) return c c determine initial conditions for wavefunction c eperp=e-wkparr**2/(C) PSI=1.D0 PSIVAC(NBARR2)=PSI imin=1 imax=NBARR2 if (eperp.lt.BARR2(1).and.eperp.lt.BARR2(NBARR2)) go to 80 c eperp above vacuum barrier; find crossing points do 30 i=1,NBARR2 if (eperp.lt.BARR2(i)) go to 40 30 continue 40 imin=i do 50 i=NBARR2,1,-1 PSIVAC(I)=PSI if (eperp.lt.BARR2(i)) go to 60 50 continue 60 imax=i if (imax.gt.imin) go to 80 write(6,*) '*** error - eperp above vacuum barrier' write(16,*) '*** error - eperp above vacuum barrier' return 80 DPSI=PSI*sqrt(c*(BARR2(imax)-eperp)) wf=psi wfderiv=dpsi c c integrate through vacuum c DELVAC=sep/float(NBARR2-1) do 100 i=imax-1,1,-1 c if ((BARR2(i)-eperp).le.0.) go to 100 PSI=PSI+DPSI*DELVAC PSIVAC(I)=PSI DPSI=DPSI+C*(BARR2(i)-eperp)*PSI*DELVAC 100 CONTINUE c c match across vacuum-semiconductor interface c PSI=PSI DPSI=DPSI*effm c c integrate through semiconductor c eperp=e+wkparr**2/(C*effm) PSI=PSI+DPSI*S2(1) PSISEM(1)=PSI ebarr=eperp-(EV+PROF2(1)) C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*S2(1) do 200 i=2,NS2 dels=(S2(I)-S2(I-1)) PSI=PSI+DPSI*dels if (.not.(dabs(psi).lt.1.d100)) go to 500 PSISEM(I)=PSI ebarr=eperp-(EV+PROF2(i)) C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*dels 200 CONTINUE c c determine amplitude, and evaluate across barrier c WKSEM=sqrt(C*effm*(EV-eperp)) phase=atan(psi*wksem/dpsi) amp=sqrt(2.)*sin(phase)/psi wf=wf*amp wfderiv=wfderiv*amp C C NORMALIZE WAVEFUNCTION C DO 250 I=NBARR2,1,-1 PSIVAC(I)=PSIVAC(I)*AMP 250 CONTINUE DO 300 I=1,NS2 PSISEM(I)=PSISEM(I)*AMP 300 CONTINUE delsmax=S2(NS2)-S2(NS2-1) if (delsmax/(2.*pi/wksem).GT.0.25) then write(6,*) '*** CAUTION *** RATIO OF SEMICOND. STEP SIZE ', & 'TO WAVELENGTH = ',delsmax/(2.*pi/wksem) write(16,*) '*** CAUTION *** RATIO OF SEMICOND. STEP SIZE ', & 'TO WAVELENGTH = ',delsmax/(2.*pi/wksem) end if return C C OVERFLOW OF WAVEFUNCTION C 500 wf=0. wfderiv=0. return end c c integrate CB wavefunction from tip across to sample, return c values of wf and derivative of wf at tip surface c subroutine CBwf(wf,wfderiv,WKSEM,e,wkparr,sep,bias,BARR2,NVDIM2, &NBARR2,PROF2,NS2,NSDIM2,S2,effm,EC,E2BAND,PSIVAC,PSISEM) dimension PROF2(NSDIM2),BARR2(NVDIM2),S2(NSDIM2) DOUBLE PRECISION PSI,DPSI,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/ PI=4.*ATAN(1.) c wf=0. wfderiv=0. WKSEM=1. eperp=e-wkparr**2/(C*effm) if (eperp.le.EC) return c c determine initial conditions for wavefunction c eperp=e-wkparr**2/(C) PSI=1.D0 PSIVAC(NBARR2)=PSI imin=1 imax=NBARR2 if (eperp.lt.BARR2(1).and.eperp.lt.BARR2(NBARR2)) go to 80 c eperp above vacuum barrier; find crossing points do 30 i=1,NBARR2 if (eperp.lt.BARR2(i)) go to 40 30 continue 40 imin=i do 50 i=NBARR2,1,-1 PSIVAC(I)=PSI if (eperp.lt.BARR2(i)) go to 60 50 continue 60 imax=i if (imax.gt.imin) go to 80 write(6,*) '*** error - eperp above vacuum barrier' write(16,*) '*** error - eperp above vacuum barrier' return 80 DPSI=PSI*sqrt(c*(BARR2(imax)-eperp)) wf=psi wfderiv=dpsi c c integrate through vacuum c DELVAC=sep/float(NBARR2-1) do 100 i=imax-1,1,-1 c if ((BARR2(i)-eperp).le.0.) go to 100 PSI=PSI+DPSI*DELVAC PSIVAC(I)=PSI DPSI=DPSI+C*(BARR2(i)-eperp)*PSI*DELVAC 100 CONTINUE c c match across vacuum-semiconductor interface c PSI=PSI DPSI=DPSI*effm c c integrate through semiconductor c eperp=e-wkparr**2/(C*effm) PSI=PSI+DPSI*S2(1) PSISEM(1)=PSI ebarr=(EC+PROF2(1))-eperp C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*S2(1) do 200 i=2,NS2 dels=(S2(I)-S2(I-1)) PSI=PSI+DPSI*dels if (.not.(dabs(psi).lt.1.d100)) go to 500 PSISEM(I)=PSI ebarr=(EC+PROF2(i))-eperp C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*dels 200 CONTINUE c c determine amplitude, and evaluate across barrier c WKSEM=sqrt(C*effm*(eperp-EC)) phase=atan(psi*wksem/dpsi) amp=sqrt(2.)*sin(phase)/psi wf=wf*amp wfderiv=wfderiv*amp C C NORMALIZE WAVEFUNCTION C DO 250 I=NBARR2,1,-1 PSIVAC(I)=PSIVAC(I)*AMP 250 CONTINUE DO 300 I=1,NS2 PSISEM(I)=PSISEM(I)*AMP 300 CONTINUE delsmax=S2(NS2)-S2(NS2-1) if (delsmax/(2.*pi/wksem).GT.0.25) then write(6,*) '*** CAUTION *** RATIO OF SEMICOND. STEP SIZE ', & 'TO WAVELENGTH = ',delsmax/(2.*pi/wksem) write(16,*) '*** CAUTION *** RATIO OF SEMICOND. STEP SIZE ', & 'TO WAVELENGTH = ',delsmax/(2.*pi/wksem) end if return C C OVERFLOW OF WAVEFUNCTION C 500 wf=0. wfderiv=0. return end c c integrate localized VB wavefunction from tip across to sample, c looking for switches in sign of wf to enumerate localized states c subroutine VBloc(nsign,wf,wfderiv,e,wkparr,sep,bias,BARR2, &NVDIM2,NBARR2,PROF2,NS2,NSDIM2,S2,effm,EV,E2BAND,PSIVAC,PSISEM) dimension PROF2(NSDIM2),BARR2(NVDIM2),S2(NSDIM2) DOUBLE PRECISION PSI,DPSI,SUM1,SUM2,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/ c wf=0. wfderiv=0. isav=NS2 eperp=e+wkparr**2/(C*effm) if (eperp.le.EV) return nsign=0 c c determine initial conditions for wavefunction c eperp=e-wkparr**2/(C) SUM1=0.D0 SUM2=0.D0 PSI=1.D0 PSIVAC(NBARR2)=PSI imin=1 imax=NBARR2 if (eperp.lt.BARR2(1).and.eperp.lt.BARR2(NBARR2)) go to 80 c eperp above vacuum barrier; find crossing points do 30 i=1,NBARR2 if (eperp.lt.BARR2(i)) go to 40 30 continue 40 imin=i do 50 i=NBARR2,1,-1 PSIVAC(I)=PSI if (eperp.lt.BARR2(i)) go to 60 50 continue 60 imax=i if (imax.gt.imin) go to 80 write(6,*) '*** error - eperp above vacuum barrier' write(16,*) '*** error - eperp above vacuum barrier' return 80 DPSI=PSI*sqrt(c*(BARR2(imax)-eperp)) wf=psi wfderiv=dpsi c c integrate through vacuum c DELVAC=sep/float(NBARR2-1) do 100 i=imax-1,1,-1 c if ((BARR2(i)-eperp).le.0.) go to 100 PSI=PSI+DPSI*DELVAC PSIVAC(I)=PSI SUM1=SUM1+PSI**2*DELVAC DPSI=DPSI+C*(BARR2(i)-eperp)*PSI*DELVAC 100 CONTINUE c c match across vacuum-semiconductor interface c PSI=PSI DPSI=DPSI*effm c c integrate through semiconductor c eperp=e+wkparr**2/(C*effm) PSI=PSI+DPSI*S2(1) PSISEM(1)=PSI SUM1=SUM1+PSI**2*S2(1) ebarr=eperp-(EV+PROF2(1)) C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*S2(1) do 200 i=2,NS2-1 dels=(S2(I)-S2(I-1)) PSISAV=PSI PSI=PSI+DPSI*dels PSISEM(I)=PSI if (PSISAV*PSI.lt.0.) then nsign=nsign+1 isav=i SUM2=SUM2+SUM1 SUM1=0.D0 end if SUM1=SUM1+PSI**2*dels ebarr=eperp-(EV+PROF2(i)) C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*dels 200 CONTINUE C C NORMALIZE WAVEFUNCTION C IF (SUM2.NE.0.) THEN AMP=1./SQRT(SUM2) ELSE AMP=1./SQRT(SUM1) END IF wf=wf*AMP wfderiv=wfderiv*AMP DO 250 I=NBARR2,1,-1 PSIVAC(I)=PSIVAC(I)*AMP 250 CONTINUE DO 300 I=1,ISAV-1 PSISEM(I)=PSISEM(I)*AMP 300 CONTINUE do 350 i=isav,NS2 psisem(i)=0. 350 continue c return end c c integrate localized CB wavefunction from tip across to sample, c looking for switches in sign of wf to enumerate localized states c subroutine CBloc(nsign,wf,wfderiv,e,wkparr,sep,bias,BARR2, &NVDIM2,NBARR2,PROF2,NS2,NSDIM2,S2,effm,EC,E2BAND,PSIVAC,PSISEM) dimension PROF2(NSDIM2),BARR2(NVDIM2),S2(NSDIM2) DOUBLE PRECISION PSI,DPSI,SUM1,SUM2,PSIVAC(NVDIM2),PSISEM(NSDIM2) c value of C below is 2m/hbar^2 in units of 1/(eV nm^2) DATA C/26.254/ c wf=0. wfderiv=0. isav=NS2 eperp=e-wkparr**2/(C*effm) if (eperp.ge.EC) return nsign=0 c c determine initial conditions for wavefunction c eperp=e-wkparr**2/(C) SUM1=0.D0 SUM2=0.D0 PSI=1.D0 PSIVAC(NBARR2)=PSI imin=1 imax=NBARR2 if (eperp.lt.BARR2(1).and.eperp.lt.BARR2(NBARR2)) go to 80 c eperp above vacuum barrier; find crossing points do 30 i=1,NBARR2 if (eperp.lt.BARR2(i)) go to 40 30 continue 40 imin=i do 50 i=NBARR2,1,-1 PSIVAC(I)=PSI if (eperp.lt.BARR2(i)) go to 60 50 continue 60 imax=i if (imax.gt.imin) go to 80 write(6,*) '*** error - eperp above vacuum barrier' write(16,*) '*** error - eperp above vacuum barrier' return 80 DPSI=PSI*sqrt(c*(BARR2(imax)-eperp)) wf=psi wfderiv=dpsi c c integrate through vacuum c DELVAC=sep/float(NBARR2-1) do 100 i=imax-1,1,-1 c if ((BARR2(i)-eperp).le.0.) go to 100 PSI=PSI+DPSI*DELVAC PSIVAC(I)=PSI SUM1=SUM1+PSI**2*DELVAC DPSI=DPSI+C*(BARR2(i)-eperp)*PSI*DELVAC 100 CONTINUE c c match across vacuum-semiconductor interface c PSI=PSI DPSI=DPSI*effm c c integrate through semiconductor c eperp=e-wkparr**2/(C*effm) PSI=PSI+DPSI*S2(1) PSISEM(1)=PSI SUM1=SUM1+PSI**2*S2(1) ebarr=(EC+PROF2(1))-eperp C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*S2(1) do 200 i=2,NS2-1 dels=(S2(I)-S2(I-1)) PSISAV=PSI PSI=PSI+DPSI*dels PSISEM(I)=PSI if (PSISAV*PSI.lt.0.) then nsign=nsign+1 isav=i SUM2=SUM2+SUM1 SUM1=0.D0 end if SUM1=SUM1+PSI**2*dels ebarr=(EC+PROF2(i))-eperp C if (ebarr.gt.0.) ebarr=ebarr-ebarr**2/E2BAND DPSI=DPSI+C*effm*(ebarr)*PSI*dels 200 CONTINUE C C NORMALIZE WAVEFUNCTION C IF (SUM2.NE.0.) THEN AMP=1./SQRT(SUM2) ELSE AMP=1./SQRT(SUM1) END IF wf=wf*AMP wfderiv=wfderiv*AMP DO 250 I=NBARR2,1,-1 PSIVAC(I)=PSIVAC(I)*AMP 250 CONTINUE DO 300 I=1,ISAV-1 PSISEM(I)=PSISEM(I)*AMP 300 CONTINUE do 350 i=isav,NS2 psisem(i)=0. 350 continue c return end