C ******************** SEMITIP2 ************************ C C E-FIELD COMPUTATIONS FOR TIP IN PROXIMITY TO SEMICONDUCTOR, C WITH CIRCULAR SYMMETRY C VERSION 1.0 - WRITTEN BY R. M. FEENSTRA, DEC 2002 C 2.0 - OCT/04, USE THIRD ORDER MATCHING AT INTERFACE C USE VARIABLE GRID RADIALLY AND INTO SEMICOND. C INCLUDE SURFACE STATES AND USE TIP MASKING ARRAY C CHECK LAST TWO ITERATIONS FOR STOPPING CRITERION C REVISE CHECKING OF ENERGY VALUES IN TABLE C 2.1 - DEC/04, REVISE STOPPING CRITERION IN GOLDEN SECTION SEARCHES C 2.2 - MAR/05, FURTHER REVISION (FOR BIAS=0 CASE) C 3.0 - NOV/07, ADD IINV PARAMETER, CHANGE COORDINATE SYSTEM TO EXACTLY C MATCH TIP (INTRODUCE C,Z0; MODIFY DELV,A,ETAT) 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 DEGENERATE DOPING, 0 OTHERWISE C IINV=1 OR 2 TO SUPPRESS VB OR CB OCCUPATION, 0 OTHERWISE C C INPUT AND OUTPUT PARAMETERS: C SEP=SEPARATION BETWEEN SEMICONDUCTOR AND END OF TIP (NM) C RAD=RADIUS OF HYPERBOLOID WHICH FORMS TIP (NM) C SLOPE=SLOPE OF TIP SHANK C ETAT=ETA PARAMETER FOR THE HYPERBOLIC COORDINATE SYSTEM C A=LOCATION OF FOCUS FOR THE HYPERBOLIC COORDINATE SYSTEM C Z0=LOCATION OF ORIGIN FOR THE HYPERBOLIC COORDINATE SYSTEM C C=ORIGIN PARAMETER FOR FOR THE HYPERBOLIC COORDINATE SYSTEM C VAC=ARRAY OF POTENTIAL VALUES IN VACUUM C SEM=ARRAY OF POTENTIAL VALUES IN SEMICONDUCTOR C TIP=ARRAY WITH VALUES OF .TRUE. IF INSIDE TIP, .FALSE. OTHERWISE C VSINT=ARRAY OF POTENTIAL VALUES AT SEMICONDUCTOR SURFACE C R=ARRAY OF R VALUES C S=ARRAY OF Z VALUES IN SEMICONDUCTOR C DELV=ARRAY OF Z-SPACING VALUES IN VACUUM C DELR=R-SPACING (FOR SMALL R) IN VACUUM AND SEMICONDUCTOR C DELS=Z-SPACING (FOR SMALL Z) IN SEMICONDUCTOR C NRDIM=R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS C NVDIM=Z-DIMENSION FOR VAC ARRAY C NSDIM=Z-DIMENSION FOR SEM ARRAY C NR=CURRENT VALUE OF R-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS C NV=CURRENT VALUE OF Z-DIMENSION FOR VAC ARRAY C NS=CURRENT VALUE OF Z-DIMENSION FOR SEM ARRAY C BIAS=BIAS OF TIP RELATIVE TO A POINT FAR INSIDE THE SEMICONDUCTOR C IWRIT >= 1 IF WANT WRITTEN OUTPUT, 0 OTHERWISE C NE=NUMBER OF POINTS IN CHARGE DENSITY TABLE C ESTART=STARTING ENERGY IN TABLE C DELE=ENERGY STEP SIZE IN TABLE C RHOC=CHARGE DENSITY TABLE C EF=FERMI-LEVEL POSITION C ITMAX=ARRAY OF MAXIMUM ITERATION SPECIFICATIONS C EP=ARRAY OF STOPPING CRITERIONS C IPMAX=MAXIMUM NUMBER OF SCALING LOOPS FOR GRID DOUBLING C PHI0=MAX SURFACE BAND BENDING C IERR=ERROR INDICATION (=1 IF GRID SIZE TOO SMALL, =2 IF ENERGY OUTSIDE TABLE) C C C SEMMIN USED TO FIND OPTIMAL UPDATE EQUATION FOR SOLN IN SEMICONDUCTOR C FUNCTION SEMMIN(PHI,STEMP,DENOM,RHOC,NE,IERR) DIMENSION RHOC(NE) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF DATA EEP/1.80943E-20/ ENER=EF-PHI IENER=NINT((ENER-ESTART)/DELE)+1 IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 100 RHO=RHOC(IENER) GO TO 110 100 RHO=RHOTOT(ENER,0.) IERR=2 110 CONTINUE TEMP=STEMP-RHO*EEP/EPSIL SEMMIN=ABS(PHI-TEMP/DENOM) RETURN END C C GOLDEN SECTION SEARCH ROUTINE (USED WITH SEMMIN) C SUBROUTINE GSECT2(XMIN,XMAX,EP,STEMP,DENOM,RHOC,NE,IERR) DIMENSION RHOC(NE) DATA GS/0.3819660/ IF (XMAX.EQ.XMIN) RETURN IF (EP.EQ.0.) RETURN IF (XMAX.LT.XMIN) THEN TEMP=XMAX XMAX=XMIN XMIN=TEMP END IF DELX=XMAX-XMIN XA=XMIN+DELX*GS FA=SEMMIN(XA,STEMP,DENOM,RHOC,NE,IERR) XB=XMAX-DELX*GS FB=SEMMIN(XB,STEMP,DENOM,RHOC,NE,IERR) 100 DELXSAV=DELX IF (DELX.LT.EP) RETURN IF (FB.LT.FA) GO TO 200 XMAX=XB DELX=XMAX-XMIN IF (DELX.EQ.DELXSAV) RETURN XB=XA FB=FA XA=XMIN+DELX*GS FA=SEMMIN(XA,STEMP,DENOM,RHOC,NE,IERR) GO TO 100 200 XMIN=XA DELX=XMAX-XMIN IF (DELX.EQ.DELXSAV) RETURN XA=XB FA=FB XB=XMAX-DELX*GS FB=SEMMIN(XB,STEMP,DENOM,RHOC,NE,IERR) GO TO 100 END C C SURFMIN USED TO FIND OPTIMAL UPDATE EQUATION FOR SOLN ON SURFACE C FUNCTION SURFMIN(PHI,RI,STEMP,DENOM,RHOS,NE,IERR) DIMENSION RHOS(NE) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF DATA EEP/1.80943E-20/ ENER=EF-PHI DELE1=DELE IENER=NINT((ENER-ESTART)/DELE)+1 IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 100 RHO=RHOS(IENER) GO TO 110 100 RHO=SRHO(ENER,DELE1) IERR=2 110 CONTINUE C ADD STATEMENT BELOW TO RESTRICT REGION OF SURFACE CHARGE C IF (RI.LT.0.) RHO=0. TEMP=STEMP-RHO*EEP*1.E7 SURFMIN=ABS(PHI-TEMP/DENOM) RETURN END C C GOLDEN SECTION SEARCH ROUTINE (USED WITH SURFMIN) C SUBROUTINE GSECT3(XMIN,XMAX,EP,RI,STEMP,DENOM,RHOS,NE,IERR) DIMENSION RHOS(NE) DATA GS/0.3819660/ IF (XMAX.EQ.XMIN) RETURN IF (EP.EQ.0.) RETURN IF (XMAX.LT.XMIN) THEN TEMP=XMAX XMAX=XMIN XMIN=TEMP END IF DELX=XMAX-XMIN XA=XMIN+DELX*GS FA=SURFMIN(XA,RI,STEMP,DENOM,RHOS,NE,IERR) XB=XMAX-DELX*GS FB=SURFMIN(XB,RI,STEMP,DENOM,RHOS,NE,IERR) 100 DELXSAV=DELX IF (DELX.LT.EP) RETURN IF (FB.LT.FA) GO TO 200 XMAX=XB DELX=XMAX-XMIN IF (DELX.EQ.DELXSAV) RETURN XB=XA FB=FA XA=XMIN+DELX*GS FA=SURFMIN(XA,RI,STEMP,DENOM,RHOS,NE,IERR) GO TO 100 200 XMIN=XA DELX=XMAX-XMIN IF (DELX.EQ.DELXSAV) RETURN XA=XB FA=FB XB=XMAX-DELX*GS FB=SURFMIN(XB,RI,STEMP,DENOM,RHOS,NE,IERR) GO TO 100 END C C MAIN PROGRAM C SUBROUTINE SEMITIP2(SEP,RAD,SLOPE,ETAT,A,Z0,C,VAC,TIP,SEM, &VSINT,R,S,DELV,DELR1,DELS1,NRDIM,NVDIM,NSDIM,NR1,NV1,NS1,BIAS1, &IWRIT,NE,ESTART1,DELE1,RHOC,RHOS,EF1,ITMAX,EP,IPMAX,PHI0,IERR) C DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM), &R(NRDIM),DELR(NRDIM),DELV(NRDIM),S(NSDIM),DELS(NSDIM),ITMAX(10), &EP(10),RHOC(NE),RHOS(NE) LOGICAL TIP(NRDIM,NVDIM) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG,IINV COMMON/SIZE/NR,NV,NS COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A1,ESTART,DELE,EF DATA PI/3.141592654/ BIAS=BIAS1 ESTART=ESTART1 DELE=DELE1 EF=EF1 NR=NR1 NV=NV1 NS=NS1 DELR0=DELR1 DELS0=DELS1 IERR=0 C C CONSTRUCT THE TIP AND VACUUM GRID AND INITIALIZE POTENTIAL C ETAT=1./SQRT(1.+1./SLOPE**2) A=RAD*SLOPE**2/ETAT A1=A SPRIME=A*ETAT Z0=SEP-SPRIME C=Z0/SPRIME DELETA=ETAT/FLOAT(NV) CETAT=ALOG((1.+ETAT)/(1.-ETAT)) IF (IWRIT.EQ.0) GO TO 100 WRITE(6,*) 'ETAT, A, Z0, C =',ETAT,A,Z0,C WRITE(16,*) 'ETAT, A, Z0, C =',ETAT,A,Z0,C WRITE(6,*) 'NR,NS,NV =',NR,NS,NV WRITE(16,*) 'NR,NS,NV =',NR,NS,NV WRITE(6,*) 'DELR,DELS,DELV =',DELR0,DELS0,(1.+C)*A*DELETA WRITE(16,*) 'DELR,DELS,DELV =',DELR0,DELS0,(1.+C)*A*DELETA 100 CONTINUE IF (NR.GT.NRDIM.OR.NV.GT.NVDIM.OR.NS.GT.NSDIM) GO TO 600 DO 150 I=1,NR R(I)=(2*NR*DELR0/PI)*TAN(PI*(I-0.5)/(2.*NR)) DELR(I)=DELR0/(COS(PI*(I-0.5)/(2.*NR)))**2 DELV(I)=(SQRT(A**2+R(I)**2)+C*A)*ETAT/FLOAT(NV) DO 120 J=1,NV-1 X2M1=(R(I)/A)**2 XSI=DSQRT(1.D0+X2M1) ETA=J*ETAT/FLOAT(NV) Z=A*ETA*(XSI+C) ZP=Z*(J+0.5)/FLOAT(J) RP=A*SQRT(X2M1*(1.-ETA**2)) IF (ZP.GT.(A*ETAT*(SQRT(1.+RP**2/((1.-ETAT**2)*A**2))+C) & -P(RP))) GO TO 110 C ADD STATEMENT BELOW FOR FLAT TIP APEX C IF (Z.GT.(SEP-P(RP)).AND.RP.LT.10.) GO TO 110 VAC(1,I,J)=0. VAC(2,I,J)=0. TIP(I,J)=.FALSE. GO TO 120 110 VAC(1,I,J)=BIAS VAC(2,I,J)=BIAS TIP(I,J)=.TRUE. 120 CONTINUE VAC(1,I,NV)=BIAS VAC(2,I,NV)=BIAS TIP(I,NV)=.TRUE. 150 CONTINUE DO 160 I=1,NR VSINT(1,I)=0. VSINT(2,I)=0. 160 CONTINUE DO 180 J=1,NS IF (J.EQ.NS) GO TO 165 S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS)) DELS(J)=DELS0/(COS(PI*J/(2.*NS)))**2 165 DO 170 I=1,NR SEM(1,I,J)=0. SEM(2,I,J)=0. 170 CONTINUE 180 CONTINUE IF (.NOT.TIP(1,1)) GO TO 190 IERR=1 WRITE(6,*) '*** ERROR - VACUUM GRID SPACING TOO LARGE' WRITE(16,*) '*** ERROR - VACUUM GRID SPACING TOO LARGE' 190 CONTINUE WRITE(6,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1) WRITE(16,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1) IF (IERR.EQ.1) GO TO 600 C C INITIAL GUESS C DO 250 I=1,NR DO 210 J=NV,1,-1 IF (.NOT.TIP(I,J)) GO TO 220 210 CONTINUE 220 JSAV=J DO 230 J=JSAV,1,-1 ETA=J*DELETA*FLOAT(NV)/FLOAT(JSAV+1) VAC(1,I,J)=BIAS*ALOG((1.+ETA)/(1.-ETA))/CETAT VAC(2,I,J)=VAC(1,I,J) 230 CONTINUE 250 CONTINUE PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8. C C ITERATE POISSON'S EQN C DO 590 IP=1,IPMAX C IF (IWRIT.NE.0) WRITE(6,*) 'SOLUTION #',IP IF (IWRIT.NE.0) WRITE(16,*) 'SOLUTION #',IP ITM=ITMAX(IP) EPI=EP(IP) CALL ITER2(VAC,TIP,SEM,VSINT,R,DELR,DELV,DELS,NRDIM,NVDIM, & NSDIM,EPI,ITM,PHI0,IWRIT,ETAT,C,RHOC,RHOS,NE,IERR) IF (IP.EQ.IPMAX) GO TO 600 C IF (NR*2.GT.NRDIM.OR.NV*2.GT.NVDIM.OR.NS*2.GT.NSDIM) GO TO 500 C C DOUBLE ALL ARRAYS C NR=NR*2 NS=NS*2 NV=NV*2 DELR0=DELR0/2. DELS0=DELS0/2. DELETA=DELETA/2. IF (IWRIT.EQ.0) GO TO 300 WRITE(6,*) 'NR,NS,NV =',NR,NS,NV WRITE(6,*) 'DELR,DELS,DELV =',DELR0,DELS0,(1.+C)*A*DELETA WRITE(16,*) 'NR,NS,NV =',NR,NS,NV WRITE(16,*) 'DELR,DELS,DELV =',DELR0,DELS0,(1.+C)*A*DELETA 300 CONTINUE DO 450 I=1,NR R(I)=(2*NR*DELR0/PI)*TAN(PI*(I-0.5)/(2.*NR)) DELR(I)=DELR0/(COS(PI*(I-0.5)/(2.*NR)))**2 DELV(I)=(SQRT(A**2+R(I)**2)+C*A)*ETAT/FLOAT(NV) DO 430 J=1,NV-1 X2M1=(R(I)/A)**2 XSI=DSQRT(1.D0+X2M1) ETA=J*ETAT/FLOAT(NV) Z=A*ETA*(XSI+C) ZP=Z*(J+0.5)/FLOAT(J) RP=A*SQRT(X2M1*(1.-ETA**2)) IF (ZP.GT.(A*ETAT*(SQRT(1.+RP**2/((1.-ETAT**2)*A**2)) & +C)-P(RP))) GO TO 410 C ADD STATEMENT BELOW FOR FLAT TIP APEX C IF (Z.GT.(SEP-P(RP)).AND.RP.LT.10.) GO TO 410 IF (J.EQ.1) GO TO 390 IF (MOD(J,2).NE.0) GO TO 380 TEMP=VAC(2,(I-1)/2+1,J/2) GO TO 400 380 TEMP=(VAC(2,(I-1)/2+1,J/2)+VAC(2,(I-1)/2+1,(J+1)/2))/2. GO TO 400 390 TEMP=(VAC(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2. 400 VAC(1,I,J)=TEMP TIP(I,J)=.FALSE. GO TO 430 410 VAC(1,I,J)=BIAS TIP(I,J)=.TRUE. 430 CONTINUE VAC(1,I,NV)=BIAS TIP(I,NV)=.TRUE. 450 CONTINUE DO 460 I=1,NR VSINT(1,I)=VSINT(2,(I-1)/2+1) 460 CONTINUE DO 463 J=1,NS-1 S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS)) DELS(J)=DELS0/(COS(PI*J/(2.*NS)))**2 463 CONTINUE DO 465 I=1,NR SEM(1,I,1)=(SEM(2,(I-1)/2+1,1)+VSINT(2,(I-1)/2+1))/2. 465 CONTINUE DO 475 J=2,NS,2 DO 470 I=1,NR SEM(1,I,J)=SEM(2,(I-1)/2+1,J/2) 470 CONTINUE 475 CONTINUE DO 485 J=3,NS,2 DO 480 I=1,NR SEM(1,I,J)=(SEM(2,(I-1)/2+1,J/2)+ & SEM(2,(I-1)/2+1,(J+1)/2))/2. 480 CONTINUE 485 CONTINUE DO 492 J=1,NV DO 490 I=1,NR VAC(2,I,J)=VAC(1,I,J) 490 CONTINUE 492 CONTINUE DO 494 I=1,NR VSINT(2,I)=VSINT(1,I) 494 CONTINUE DO 498 J=1,NS DO 496 I=1,NR SEM(2,I,J)=SEM(1,I,J) 496 CONTINUE 498 CONTINUE WRITE(6,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1) WRITE(16,*) 'LARGEST RADIUS, DEPTH =',R(NR),S(NS-1) GO TO 590 C 500 IF (NS*2.GT.NSDIM) GO TO 600 C C DOUBLE ONLY SEM ARRAY ELEMENTS C NS=NS*2 DELS0=DELS0/2. IF (IWRIT.EQ.0) GO TO 505 WRITE(6,*) 'NS =',NS WRITE(6,*) 'DELS =',DELS0 WRITE(16,*) 'NS =',NS WRITE(16,*) 'DELS =',DELS0 505 CONTINUE DO 508 J=1,NS-1 S(J)=(2*NS*DELS0/PI)*TAN(PI*J/(2.*NS)) DELS(J)=DELS0/COS(PI*J/(2.*NS)) 508 CONTINUE DO 510 I=1,NR SEM(1,I,1)=(SEM(2,I,1)+VSINT(2,I))/2. 510 CONTINUE DO 520 J=2,NS,2 DO 515 I=1,NR SEM(1,I,J)=SEM(2,I,J/2) 515 CONTINUE 520 CONTINUE DO 530 J=3,NS,2 DO 525 I=1,NR SEM(1,I,J)=(SEM(2,I,J/2)+ & SEM(2,I,(J+1)/2))/2. 525 CONTINUE 530 CONTINUE DO 550 J=1,NS DO 540 I=1,NR SEM(2,I,J)=SEM(1,I,J) 540 CONTINUE 550 CONTINUE WRITE(6,*) 'LARGEST DEPTH =',S(NS-1) WRITE(16,*) 'LARGEST DEPTH =',S(NS-1) 590 CONTINUE C 600 CONTINUE IF (IWRIT.NE.0) WRITE(6,*) 'NUMBER OF ITERATIONS = ',ITM IF (IWRIT.NE.0) WRITE(6,*) 'BAND BENDING AT MIDPOINT = ',PHI0 IF (IWRIT.NE.0) WRITE(16,*) 'NUMBER OF ITERATIONS = ',ITM IF (IWRIT.NE.0) WRITE(16,*) 'BAND BENDING AT MIDPOINT = ',PHI0 C NR1=NR NV1=NV NS1=NS DELR1=DELR0 DELS1=DELS0 RETURN END C INTEGER FUNCTION IADS(I) COMMON/SIZE/NR,NV,NS IF (I.EQ.0) GO TO 10 IF (I.GT.NR) GO TO 20 IADS=I RETURN 10 IADS=1 RETURN 20 IADS=NR RETURN END C C 3-D SOLUTION FOR POTENTIAL C SUBROUTINE ITER2(VAC,TIP,SEM,VSINT,R,DELR,DELV,DELS,NRDIM,NVDIM, &NSDIM,EP,ITMAX,PHI0,IWRIT,ETAT,C,RHOC,RHOS,NE,IERR) C DIMENSION VAC(2,NRDIM,NVDIM),SEM(2,NRDIM,NSDIM),VSINT(2,NRDIM), &R(NRDIM),DELR(NRDIM),DELV(NRDIM),DELS(NSDIM),RHOC(NE), &RHOS(NE) LOGICAL TIP(NRDIM,NVDIM) COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,DEGEN COMMON/SIZE/NR,NV,NS COMMON/STIP/BIAS,DELR0,DELS0,DELETA,A,ESTART,DELE,EF DATA EEP/1.80943E-20/ C PHISAV=0. C2=C*C C3=C2*C C2M1=C2-1. C2P1=C2+1. C2P2=C2+2. C2P6=C2+6. TC2P2=3.*C2+2. DO 500 ITER=1,ITMAX DO 360 I=1,NR X2M1=(R(I)/A)**2 XSI=DSQRT(1.D0+X2M1) XSI2=XSI*XSI XSI3=XSI2*XSI XSI4=XSI3*XSI XSI5=XSI4*XSI DELXSI=R(I)*DELR(I)/(XSI*A**2) DO 350 J=1,NV-1 IF (TIP(I,J)) GO TO 360 ETA=J*DELETA ETA2=ETA*ETA ETA3=ETA2*ETA ETA4=ETA3*ETA ETA5=ETA4*ETA OME2=1.-ETA**2 X2ME2=XSI**2-ETA**2 X2ME2C=XSI*(XSI+C)-ETA**2*(C*XSI+1.) X2ME2C2=X2ME2C*X2ME2C T1=X2M1*((XSI+C)**2-ETA**2*(XSI*C*2.+C2+1.))/X2ME2C T2=OME2*X2ME2/X2ME2C T4=-C*ETA*X2M1*OME2/X2ME2C T5=(C3+3.*C2*XSI+C*C2P2*XSI2+3.*C2*XSI3+4.*C*XSI4+ &2.*XSI5+ETA4*(C3+TC2P2*XSI+C*C2P6*XSI2+3.*C2*XSI3)- &2.*ETA2*(C*C2M1+3.*C2*XSI+C*C2P6*XSI2+TC2P2*XSI3+C*XSI4))/X2ME2C2 T6=-ETA*(C2+4.*C*XSI+C2*XSI2+2.*XSI4+ &ETA4*(2.+C2+4.*C*XSI+C2*XSI2)- &2.*ETA2*(C2+4.*C*XSI+C2P2*XSI2))/X2ME2C2 C IF (J.EQ.1) GO TO 335 TEMP=T1* &(VAC(1,IADS(I+1),J)+VAC(1,IADS(I-1),J))/DELXSI**2+ &T2*(VAC(1,I,J+1)+VAC(1,I,J-1))/DELETA**2+ &T4*(VAC(1,IADS(I+1),J+1)-VAC(1,IADS(I-1),J+1)- &VAC(1,IADS(I+1),J-1)+VAC(1,IADS(I-1),J-1))/(2.*DELXSI*DELETA) &+T5*(VAC(1,IADS(I+1),J)-VAC(1,IADS(I-1),J))/(2.*DELXSI) &+T6*(VAC(1,I,J+1)-VAC(1,I,J-1))/(2.*DELETA) GO TO 340 335 TEMP=T1* &(VAC(1,IADS(I+1),1)+VAC(1,IADS(I-1),1))/DELXSI**2+ &T2*(VAC(1,I,2)+VSINT(1,I))/DELETA**2+ &T4*(VAC(1,IADS(I+1),2)-VAC(1,IADS(I-1),2)- &VSINT(1,IADS(I+1))+VSINT(1,IADS(I-1)))/(2.*DELXSI*DELETA) &+T5*(VAC(1,IADS(I+1),1)-VAC(1,IADS(I-1),1))/(2.*DELXSI) &+T6*(VAC(1,I,2)-VSINT(1,I))/(2.*DELETA) 340 CONTINUE VAC(2,I,J)=TEMP/(2.*T1/DELXSI**2+2.*T2/DELETA**2) 350 CONTINUE 360 CONTINUE DO 370 J=1,NV DO 365 I=1,NR VAC(1,I,J)=VAC(2,I,J) 365 CONTINUE 370 CONTINUE DO 400 I=1,NR SURFOLD=VSINT(1,I) ENER=EF-VSINT(1,I) IF (TIP(I,3)) GO TO 375 C 3RD ORDER SOLUTION IN SEMICONDUCTOR & VACUUM STEMP= & ((3.2*VAC(1,I,1)-1.7*VAC(1,I,2)+0.4*VAC(1,I,3))/DELV(I)+ & EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0) DENOM=(1.9/DELV(I)+1.9*EPSIL/DELS0) GO TO 385 C 3RD ORDER SOLUTION IN SEMICONDUCTOR, 2ND ORDER IN VACUUM 375 IF (TIP(I,2)) GO TO 380 STEMP= & ((2.0*VAC(1,I,1)-0.5*VAC(1,I,2))/DELV(I)+ & EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0) DENOM=(1.5/DELV(I)+1.9*EPSIL/DELS0) GO TO 385 C 3RD ORDER SOLUTION IN SEMICONDUCTOR, 1ST ORDER IN VACUUM 380 STEMP= & (VAC(1,I,1)/DELV(I)+ & EPSIL*(3.2*SEM(1,I,1)-1.7*SEM(1,I,2)+0.4*SEM(1,I,3))/DELS0) DENOM=(1./DELV(I)+1.9*EPSIL/DELS0) 385 CONTINUE IENER=NINT((ENER-ESTART)/DELE)+1 IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 388 RHO=RHOS(IENER) GO TO 390 388 RHO=SRHO(ENER,DELE) IERR=2 390 CONTINUE C ADD STATEMENT BELOW TO RESTRICT REGION OF SURFACE CHARGE C IF (R(I).LT.0.) RHO=0. TEMP=STEMP-RHO*EEP*1.E7 SURFNEW=TEMP/DENOM DELSURF=ABS(BIAS)/1.E6 CALL GSECT3(SURFOLD,SURFNEW,DELSURF,R(I),STEMP,DENOM,RHOS,NE & ,IERR) VSINT(2,I)=(SURFOLD+SURFNEW)/2. 400 CONTINUE DO 410 I=1,NR VSINT(1,I)=VSINT(2,I) 410 CONTINUE IF (MOD(ITER,100).EQ.0) THEN PHISAV2=PHISAV PHISAV=PHI0 PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8. IF (IWRIT.NE.0) WRITE(6,*) 'ITER,PHI0 =',ITER,PHI0 IF (IWRIT.NE.0) WRITE(16,*) 'ITER,PHI0 =',ITER,PHI0 END IF DO 450 J=1,NS-1 DO 440 I=1,NR SEMOLD=SEM(1,I,J) ENER=EF-SEM(1,I,J) IF (J.EQ.1) GO TO 420 STEMP=(SEM(1,IADS(I+1),J)+SEM(1,IADS(I-1),J))/DELR(I)**2+ & (SEM(1,I,J+1)+SEM(1,I,J-1))/DELS(J)**2+ & (SEM(1,IADS(I+1),J)-SEM(1,IADS(I-1),J))/(2.*R(I)*DELR(I)) GO TO 430 420 STEMP=(SEM(1,IADS(I+1),1)+SEM(1,IADS(I-1),1))/DELR(I)**2+ & (SEM(1,I,2)+VSINT(1,I))/DELS(J)**2+ & (SEM(1,IADS(I+1),1)-SEM(1,IADS(I-1),1))/(2.*R(I)*DELR(I)) 430 CONTINUE IENER=NINT((ENER-ESTART)/DELE)+1 IF (IENER.LT.1.OR.IENER.GT.NE) GO TO 433 RHO=RHOC(IENER) GO TO 436 433 RHO=RHOTOT(ENER,0.) IERR=2 436 CONTINUE TEMP=STEMP-RHO*EEP/EPSIL DENOM=(2./DELR(I)**2+2./DELS(J)**2) SEMNEW=TEMP/DENOM DELSEM=AMAX1(1.E-6,ABS(BIAS)/1.E6) CALL GSECT2(SEMOLD,SEMNEW,DELSEM,STEMP,DENOM,RHOC,NE, & IERR) SEM(2,I,J)=(SEMOLD+SEMNEW)/2. 440 CONTINUE 450 CONTINUE DO 480 J=1,NS DO 470 I=1,NR SEM(1,I,J)=SEM(2,I,J) 470 CONTINUE 480 CONTINUE IF ((MOD(ITER,100).EQ.0.AND.ABS(PHI0-PHISAV).LT.EP.AND. & ABS(PHISAV-PHISAV2).LT.2.*EP)) GO TO 510 500 CONTINUE 510 ITMAX=ITER PHI0=(9.*VSINT(1,1)-VSINT(1,2))/8. RETURN END