C ************** SEMIRHO.F **************** C C CHARGE DENSITIES FOR A SEMICONDUCTOR C VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, DEC 2002 C 1.1 - JUL/03, CHECK EXPONENTS IN RHOD AND RHOA C 1.2 - AUG/03, ADD T=0 CASES TO EFFIND C 1.3 - APR/04, LIMIT ETA<-8 CASE IN FJINT C 1.4 - NOV/07, ADD IINV PARAMETER C 1.5 - FEB/09, ADD CASE FOR CD=CA IN EFFIND C VERSION 2.0 - MAR/09, ADD ADDITIONAL TABLES TO SEMIRHO (FOR SELF-CONSISTENT C COMPUTATIONS WITH QUANTUM CHARGE DENSITIES) C C SEMI COMMON BLOCK: C EGAP=BAND GAP (EV) C ED=DONOR BINDING ENERGY (EV) C EA=ACCEPTOR BINDING ENERGY (EV) C ACB=CONDUCTION BAND EFFECTIVE MASS C AVB=VALENCE BAND EFFECTIVE MASS C CD=DONOR DENSITY (CM^-3) C CA=ACCEPTOR DENSITY (CM^-3) C EPSIL=DIELECTRIC CONSTANT OF SEMICONDUCTOR C TK=TEMPERATURE (EV) C IDEG=1 IF DONORS & ACCEPTORS ARE DEGENERATE, 0 OTHERWISE C IINV=1 OR 2 TO SUPPRESS VB OR CB OCCUPATION, 3 FOR BOTH, 0 OTHERWISE C C INPUT PARAMETERS: C DELE=SPACING FOR CHARGE DENSITY TABLE (EV) C ESTART=STARTING ENERGY FOR CHARGE DENSITY TABLE (EV) C NE=NUMBER OF POINTS IN CHARGE DENSITY TABLE C C OUTPUT: C EF=FERMI-LEVEL POSITION RELATIVE TO VB MAX C RHOC=TABLE OF CHARGE DENSITIES (CM^-3) C C C FERMI-DIRAC INTEGRAL, J=3 OR 1 FOR 3/2 OR 1/2 INTEGRAL C REAL FUNCTION FJINT(J,ETA) EXTERNAL FJ COMMON/FJCOM/JF,ETAF DATA PI/3.141592654/ IF (ETA.GT.40) GO TO 300 IF (ETA.LT.-8) GO TO 200 JF=J ETAF=ETA FJINT=TRAP(FJ,0.,20.+ETA,1000) RETURN 200 FJINT=SQRT(PI)*EXP(AMAX1(-40.,ETA))/2. IF (J.EQ.3) FJINT=FJINT*3./2. RETURN 300 FJINT=SQRT(ETA**(J+2))/(J/2.+1) RETURN END C C TRAPEZOIDAL INTEGRATION ROUTINE C FUNCTION TRAP(F,XMIN,XMAX,NSTEP) DOUBLE PRECISION SUM DELX=(XMAX-XMIN)/FLOAT(NSTEP) SUM=(F(XMIN)+F(XMAX))/2. IF (NSTEP.LT.2) GO TO 200 DO 100 I=1,NSTEP-1 X=XMIN+I*DELX SUM=SUM+F(X) 100 CONTINUE 200 SUM=SUM*DELX TRAP=SUM RETURN END C FUNCTION FJ(X) COMMON/FJCOM/J,ETA FJ=SQRT(X**J)/(1.+EXP(X-ETA)) RETURN END C C ELECTRON DENSITY IN CONDUCTION BAND C FUNCTION RHOC(EF,PHI) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV DATA C/6.815E21/ RHOC=0. IF (IINV.EQ.2.OR.IINV.EQ.3) RETURN IF (TK.NE.0.) GO TO 200 IF ((EF-EGAP-PHI).LE.0.) GO TO 150 RHOC=(2.*C/3.)*SQRT((ACB*(EF-EGAP-PHI))**3) RETURN 150 RHOC=0. RETURN 200 RHOC=C*SQRT((ACB*TK)**3)*FJINT(1,(EF-EGAP-PHI)/TK) RETURN END C C HOLE DENSITY IN VALENCE BAND C FUNCTION RHOV(EF,PHI) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV DATA C/6.815E21/ RHOV=0. IF (IINV.EQ.1.OR.IINV.EQ.3) RETURN IF (TK.NE.0.) GO TO 200 IF ((-EF+PHI).LE.0.) GO TO 150 RHOV=(2.*C/3.)*SQRT((AVB*(-EF+PHI))**3) RETURN 150 RHOV=0. RETURN 200 RHOV=C*SQRT((AVB*TK)**3)*FJINT(1,(-EF+PHI)/TK) RETURN END C C IONIZED DONOR DENSITY C FUNCTION RHOD(EF,PHI) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV RHOD=CD IF (IDEG.EQ.1) RETURN EXPO=EF-EGAP+ED-PHI IF (TK.NE.0.) GO TO 50 IF (EXPO.LE.0.) GO TO 100 IF (EXPO.GT.0.) GO TO 200 50 EXPO=EXPO/TK IF (EXPO.LT.-40) GO TO 100 IF (EXPO.GT.40) GO TO 200 RHOD=CD/(1.+2.*EXP(EXPO)) RETURN 100 RHOD=CD RETURN 200 RHOD=0. RETURN END C C IONIZED ACCEPTOR DENSITY C FUNCTION RHOA(EF,PHI) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV RHOA=CA IF (IDEG.EQ.1) RETURN EXPO=EA-EF+PHI IF (TK.NE.0.) GO TO 50 IF (EXPO.LE.0.) GO TO 100 IF (EXPO.GT.0.) GO TO 200 50 EXPO=EXPO/TK IF (EXPO.LT.-40) GO TO 100 IF (EXPO.GT.40) GO TO 200 RHOA=CA/(1.+4.*EXP(EXPO)) RETURN 100 RHOA=CA RETURN 200 RHOA=0. RETURN END C C TOTAL DENSITY OF POSITIVE CHARGES C FUNCTION RHOTOT(EF,PHI) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV RHOE=RHOC(EF,PHI)+RHOA(EF,PHI) RHOH=RHOV(EF,PHI)+RHOD(EF,PHI) RHOTOT=-RHOE+RHOH RETURN END C FUNCTION ARHO(EF1) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV EF=EF1 TEMP=0. ARHO=ABS(RHOTOT(EF,TEMP)) RETURN END C C GOLDEN SECTION SEARCH ROUTINE (USED FOR FINDING FERMI-LEVEL) C SUBROUTINE GSECT1(F,XMIN,XMAX,EP) DATA GS/0.3819660/ DELX=XMAX-XMIN F1=F(XMIN) F2=F(XMAX) XA=XMIN+DELX*GS FA=F(XA) XB=XMAX-DELX*GS FB=F(XB) 100 IF (DELX.LT.EP) RETURN IF (FB.LT.FA) GO TO 200 XMAX=XB DELX=XMAX-XMIN XB=XA FB=FA XA=XMIN+DELX*GS FA=F(XA) GO TO 100 200 XMIN=XA DELX=XMAX-XMIN XA=XB FA=FB XB=XMAX-DELX*GS FB=F(XB) GO TO 100 END C C FIND FERMI-LEVEL C SUBROUTINE EFFIND(EF) C COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV EXTERNAL ARHO C IINVSAV=IINV IINV=0 IF (CD.EQ.0.AND.CA.EQ.0.) GO TO 300 IF (TK.EQ.0) GO TO 200 ESTART=-0.1 NE=1000 DELE=(EGAP+0.2)/FLOAT(NE) RMIN=ARHO(ESTART) IESAV=1 DO 100 IE=2,NE ENER=ESTART+(IE-1)*DELE RTMP=ARHO(ENER) IF (RTMP.GT.RMIN) GO TO 100 IESAV=IE RMIN=RTMP 100 CONTINUE EFMIN=ESTART+(IESAV-2)*DELE EFMAX=ESTART+(IESAV)*DELE CALL GSECT1(ARHO,EFMIN,EFMAX,1.E-6) EF=(EFMIN+EFMAX)/2. GO TO 400 200 IF (CA.GT.CD) EF=EA/2. IF (CA.LT.CD) EF=EGAP-(ED/2.) IF (CA.EQ.CD) EF=(EGAP-ED+EA)/2. GO TO 400 300 EF=EGAP/2.+0.75*TK*ALOG(AVB/ACB) 400 IINV=IINVSAV RETURN END C C OUTPUT TABLE OF CHARGE DENSITY VALUES C SUBROUTINE SEMIRHO(DELE,ESTART,NE,NEDIM,RHOB) C DIMENSION RHOB(NEDIM) C IF (NE.GT.NEDIM) THEN WRITE(6,*) '*** ERROR - NE > NEDIM; PROGRAM HALTED' WRITE(6,*) 'TYPE ENTER TO CONTINUE' READ(5,*) STOP END IF DO 300 I=1,NE EF1=(I-1)*DELE+ESTART RHOB(I)=RHOTOT(EF1,0.) 300 CONTINUE C RETURN END