C ************** SEMIRHOMULT.F **************** C C CHARGE DENSITIES FOR A SEMICONDUCTOR CONTAINING MULTIPLE REGIONS C C VERSION 6.0 - FEB/11, DERIVED FROM semirho-6.0 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 FERMI-DIRAC OCCUPATION FACTOR 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 ELECTRON DENSITY IN CONDUCTION BAND C FUNCTION RHOCB(IREG,EF,Pot) PARAMETER(NREGDIM=2) COMMON/SEMI/TK,EGAP(NREGDIM),ED(NREGDIM),EA(NREGDIM),ACB(NREGDIM), &AVB(NREGDIM),CD(NREGDIM),CA(NREGDIM),IDEG(NREGDIM),IINV(NREGDIM), &DELVB(NREGDIM) c constant below is (2/rt(pi))*2*(m/(2*pi*hbar^2))^1.5 in eV^-1.5 cm^-3 DATA C/6.815E21/ RHOCB=0. IF (IINV(IREG).EQ.2.OR.IINV(IREG).EQ.3) RETURN IF (TK.NE.0.) GO TO 200 IF ((EF-EGAP(IREG)-DELVB(IREG)-Pot).LE.0.) GO TO 150 RHOCB=(2.*C/3.)* & SQRT((ACB(IREG)*(EF-EGAP(IREG)-DELVB(IREG)-Pot))**3) RETURN 150 RHOCB=0. RETURN 200 RHOCB=C*SQRT((ACB(IREG)*TK)**3)* & FJINT(1,(EF-EGAP(IREG)-DELVB(IREG)-Pot)/TK) RETURN END C C HOLE DENSITY IN VALENCE BAND C FUNCTION RHOVB(IREG,EF,Pot) PARAMETER(NREGDIM=2) COMMON/SEMI/TK,EGAP(NREGDIM),ED(NREGDIM),EA(NREGDIM),ACB(NREGDIM), &AVB(NREGDIM),CD(NREGDIM),CA(NREGDIM),IDEG(NREGDIM),IINV(NREGDIM), &DELVB(NREGDIM) c constant below is (2/rt(pi))*2*(m/(2*pi*hbar^2))^1.5 in eV^-1.5 cm^-3 DATA C/6.815E21/ RHOVB=0. IF (IINV(IREG).EQ.1.OR.IINV(IREG).EQ.3) RETURN IF (TK.NE.0.) GO TO 200 IF ((-EF+DELVB(IREG)+Pot).LE.0.) GO TO 150 RHOVB=(2.*C/3.)*SQRT((AVB(IREG)*(-EF+DELVB(IREG)+Pot))**3) RETURN 150 RHOVB=0. RETURN 200 RHOVB=C*SQRT((AVB(IREG)*TK)**3)*FJINT(1,(-EF+DELVB(IREG)+Pot)/TK) RETURN END C C IONIZED DONOR DENSITY C FUNCTION RHOD(IREG,EF,Pot) PARAMETER(NREGDIM=2) COMMON/SEMI/TK,EGAP(NREGDIM),ED(NREGDIM),EA(NREGDIM),ACB(NREGDIM), &AVB(NREGDIM),CD(NREGDIM),CA(NREGDIM),IDEG(NREGDIM),IINV(NREGDIM), &DELVB(NREGDIM) RHOD=CD(IREG) IF (IDEG(IREG).EQ.1) RETURN EXPO=EF-EGAP(IREG)-DELVB(IREG)+ED(IREG)-Pot 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(IREG)/(1.+2.*EXP(EXPO)) RETURN 100 RHOD=CD(IREG) RETURN 200 RHOD=0. RETURN END C C IONIZED ACCEPTOR DENSITY C FUNCTION RHOA(IREG,EF,Pot) PARAMETER(NREGDIM=2) COMMON/SEMI/TK,EGAP(NREGDIM),ED(NREGDIM),EA(NREGDIM),ACB(NREGDIM), &AVB(NREGDIM),CD(NREGDIM),CA(NREGDIM),IDEG(NREGDIM),IINV(NREGDIM), &DELVB(NREGDIM) RHOA=CA(IREG) IF (IDEG(IREG).EQ.1) RETURN EXPO=EA(IREG)-EF+DELVB(IREG)+Pot 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(IREG)/(1.+4.*EXP(EXPO)) RETURN 100 RHOA=CA(IREG) RETURN 200 RHOA=0. RETURN END C C TOTAL DENSITY OF BULK CHARGES C FUNCTION RHOB(IREG,EF,Pot) PARAMETER(NREGDIM=2) RHOE=RHOCB(IREG,EF,Pot)+RHOA(IREG,EF,Pot) RHOH=RHOVB(IREG,EF,Pot)+RHOD(IREG,EF,Pot) RHOB=-RHOE+RHOH RETURN END C FUNCTION ARHO(EF) COMMON/ARHO1/IREG1 IREG=IREG1 ARHO=ABS(RHOB(IREG,EF,0.)) RETURN END C C FIND FERMI-LEVEL C SUBROUTINE EFFIND(IREG,EF) C PARAMETER(NREGDIM=2) COMMON/SEMI/TK,EGAP(NREGDIM),ED(NREGDIM),EA(NREGDIM),ACB(NREGDIM), &AVB(NREGDIM),CD(NREGDIM),CA(NREGDIM),IDEG(NREGDIM),IINV(NREGDIM), &DELVB(NREGDIM) COMMON/ARHO1/IREG1 EXTERNAL ARHO IREG1=IREG C IINVSAV=IINV(IREG) IINV(IREG)=0 IF (CD(IREG).EQ.0.AND.CA(IREG).EQ.0.) GO TO 300 IF (TK.EQ.0) GO TO 200 ESTART=-0.1+DELVB(IREG) NE=1000 DELE=(EGAP(IREG)+DELVB(IREG)+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 GSECT(ARHO,EFMIN,EFMAX,1.E-6) EF=(EFMIN+EFMAX)/2. GO TO 400 200 IF (CA(IREG).GT.CD(IREG)) EF=EA(IREG)/2. IF (CA(IREG).LT.CD(IREG)) EF=EGAP(IREG)-(ED(IREG)/2.) IF (CA(IREG).EQ.CD(IREG)) EF=(EGAP(IREG)-ED(IREG)+EA(IREG))/2. GO TO 400 300 EF=EGAP(IREG)/2.+0.75*TK*ALOG(AVB(IREG)/ACB(IREG)) 400 IINV(IREG)=IINVSAV RETURN END C C OUTPUT TABLE OF CHARGE DENSITY VALUES FOR VB, CB, AND TOTAL C SUBROUTINE SEMIRHO(IREG,DELE,ESTART,NE,NEDIM,RHOBTAB,ICOMP, &RHOCBTAB,RHOVBTAB) C PARAMETER(NREGDIM=2) DIMENSION RHOBTAB(NREGDIM,NEDIM),RHOCBTAB(NREGDIM,NEDIM), &RHOVBTAB(NREGDIM,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 RHOCBSAV=RHOCB(IREG,EF1,0.) RHOVBSAV=RHOVB(IREG,EF1,0.) RHOBTAB(IREG,I)=-RHOCBSAV-RHOA(IREG,EF1,0.)+RHOVBSAV+ & RHOD(IREG,EF1,0.) IF (ICOMP.EQ.1) THEN RHOCBTAB(IREG,I)=RHOCBSAV RHOVBTAB(IREG,I)=RHOVBSAV END IF 300 CONTINUE C RETURN END