C---------------------------------------------------------------------- C Ensemble Monte Carlo C Gallium Arsenide C---------------------------------------------------------------------- C Bulk - Three Valleys C Modified Constant Time Technique C Vectorized Version - Release 1.1 C---------------------------------------------------------------------- C Written by: C C Umberto Ravaioli and Ben Shapo C University of Illinois at Urbana-Champaign C 1987-1988 C C Revised by: C C Chris Lee C University of Illinois at Urbana-Champaign C Beckman Institite C April 1989 C--------------------------------------------------------------------- C C--------------------------------------------------------------------- C C This program contains a basic Ensemble Monte Carlo code for C study of transport in bulk III-V semiconductors. Parameters C for GaAs are used in this version. The program is designed C as a first step for a device simulation program. C C The non-parabolic approximation is used, with three valleys C in the conduction band. The L and X satellite valleys are C considered to have spherical constant energy surfaces, i.e. C the effective mass used is the conductivity mass, taken as C an average of longitudinal and transverse mass. C The tensorial properties of the effective mass in satellite C valleys will be included in the device simulations built on C this code. C C This program can be used to calculate observables of the C electron ensemble in steady state or transient conditions. C The calculation for the drift velocity as a function of C electric field is set up in this program. C C The random flight selection is accomplished using an C adaptation of the Constant Time Technique to calculate C flight times. The procedure has been vectorized to a C good extent, and has a structure suitable for parallel C computation. This code is efficient for supercomputer C applications. It runs well on scalar machine, but may C be memory limited, since most variables have to be C dimensioned for vectorization. C C The program has been written to be as much as possible C machine independent. Only trivial modifications should C be necessary to run the code on any computer supporting C standard FORTRAN 77 language. Therefore, performance C could be improved in supercomputing applications C by using special machine utilities (in particular, the C gather-scatter routines should be used when available. In C the present code, gather-scatter is implemented by indirect C addressing, for machine independence. C C Although the program has a well-defined structure, the usual C guidelines for "structured programming" have not been C followed to avoid supercomputing performance degradation, C due to excessive subroutines calls. The strategy adopted C should not detract from the understandability of the code. C C The instructions have been written in upper case. Lower case C has been reserved for comments and for "temporary" C instructions use in debugging or testing phase. This C approach makes the final "clean up" of the program much C easier. C C Unless otherwise specified, all calculations are performed in C M.K.S. units! Quantities entered in practical units are C converted to M.K.S. This choice is motivated by two reasons: C a) Normalization would give only marginal improvement in a C Monte Carlo approach, since round-off error is never an C issue (as it is for instance in models based on the C solution of partial differential equations). Careful C arrangement of the scattering rate calculations allows C one to keep the values at all times below the overflow C limit of single precision in most computers. C b) The number of variables and derived constants is very C large, and the use of M.K.S. units facilitates testing, C debugging, and modifications of the code. Also, it is C simpler to understand the code itself, because of the C perfect consistency of M.K.S. units. C C===================================================================== C ISIZE=size of the arrays. ISIZE .ge. NELEC C PARAMETER (ISIZE=20000) C C===================================================================== C C COMMON BLOCK DEFINITIONS C C===================================================================== C [1] Absolute Physical Parameters and constants C--------------------------------------------------------------------- C BOLTZ=Boltmann's Constant C EM0 =Vacuum Electron Mass C EP0 =Dielectric Permittivity C HBAR =Modified Planck's Constant C PI =pi (3.14....) C Q =Absolute Value of the Electronic Charge C TEMP =Lattice Temperature C TK =Thermal Voltage (kT) C TWOPI=2*pi (6.28...) C COMMON/PHYS1/BOLTZ,EM0,EP0,HBAR,PI,Q,TEMP,TK,TWOPI C C--------------------------------------------------------------------- C [2] Material Related Parameters C (parameters are dimensioned to easily include two different C materials in the device simulations developed on the basic C Monte Carlo code) C--------------------------------------------------------------------- C Non-parabolicity parameters C ALFA1=non-parabolicity for Gamma valley C ALFA2= " " " L " C ALFA3= " " " X " C COMMON/MAT1/ALFA1(2),ALFA2(2),ALFA3(2) C ------------------------------------- C Deformation potential for intervalley phonon. C DEFGL=Deformation potential for Gamma to L transition C DEFGX= " " " Gamma to X " C DEFLX= " " " L to X " C DEFLL= " " " L to L " C DEFXX= " " " X to X " C COMMON/MAT2/DEFGL(2),DEFGX(2),DEFLX(2),DEFLL(2),DEFXX(2) C ------------------------------------- C DEFAC=Acoustic phonon Deformation Potential C RHO =Density of the material C VS =Sound Velocity C COMMON/MAT3/DEFAC(2),RHO(2),VS(2) C ------------------------------------- C GAP =Energy gap (valence band to Gamma valley); C GAPGL=Energy gap Gamma to L valley C GAPGX=Energy gap Gamma to X valley. C COMMON/MAT4/GAP(2),GAPGL(2),GAPGX(2) C ------------------------------------- C EMSTAR1=Effective Mass coefficient in Gamma valley C EMSTAR2=Effective Mass coefficient in L valley C EMSTAR3=Effective Mass coefficient in X valley C COMMON/MAT5/EMSTAR1(2),EMSTAR2(2),EMSTAR3(2) C ------------------------------------- C EM1=Total Conductivity Effective Mass in Gamma valley. C EM2=Total Conductivity Effective Mass in L valley. C EM3=Total Conductivity Effective Mass in X valley. C COMMON/MAT6/EM1(2),EM2(2),EM3(2) C ------------------------------------- C The following are input as relative dielectric constants, C but are then converted to total dielectric constants. C EPZER=Dielectric constant at zero frequency. C EPINF=Dielectric constant at infinite frequency. C COMMON/MAT7/EPZER(2),EPINF(2) C ------------------------------------- C TPOP =Polar optical phonon temperature C TKPOP=Polar optical phonon energy (TPOP*BOLTZ) C COMMON/MAT8/TPOP(2),TKPOP(2) C ------------------------------------- C TINTGL=Intervalley phonon temperature (Gamma to L) C TINTGX=Intervalley phonon temperature (Gamma to X) C TINTLX=Intervalley phonon temperature ( L to X) C TEQINL=Intervalley phonon temperature ( L to L) C TEQINX=Intervalley phonon temperature ( X to X) C COMMON/MAT9/TINTGL(2),TINTGX(2),TINTLX(2),TEQINL(2),TEQINX(2) C ------------------------------------- C TKINT1=Intervalley phonon energy (Gamma to L) C TKINT2=Intervalley phonon energy (Gamma to X) C TKINT3=Intervalley phonon energy ( L to X) C TKEIN2=Equivalent Intervalley phonon energy ( L valley ) C TKEIN3=Equivalent Intervalley phonon energy ( X valley ) C COMMON/MAT10/TKINT1(2),TKINT2(2),TKINT3(2),TKEIN2(2),TKEIN3(2) C ------------------------------------- C ZL=Number of equivalent valleys in L C ZX=Number of equivalent valleys in X C COMMON/MAT11/ZL,ZX COMMON/MAT12/DOP(5) C ------------------------------------- C Often used constants C C CC(i,j)=2*mass/Hbar**2 [ i=valley; j=material ] C C [ In the following variables, the ending 1,2,3 indicates C the valley, while (j) indicates the material ] C C TWOA1(j), TWOA2(j), TWOA3(j) = 2*Alfa C FOURA1(j), FOURA2(j), FOURA3(j) = 4*Alfa C TAINV1(j), TAINV2(j), TAINV3(j) = 1/(2*Alfa) C COMMON/MAT13/CC(3,2) COMMON/MAT14/TWOA1(2),TWOA2(2),TWOA3(2) COMMON/MAT15/FOURA1(2),FOURA2(2),FOURA3(2) COMMON/MAT16/TAINV1(2),TAINV2(2),TAINV3(2) C--------------------------------------------------------------------- C [3] Algorithm Parameters C--------------------------------------------------------------------- C NELEC =Number of simulated electrons (in this version it C should be a multiple of 1000) C NESTEP=Number of energy steps in scattering tables C EMAX =Maximum energy in scattering tables C ESTEP =Energy step in scattering tables C COMMON/ALG1/NELEC,NESTEP,EMAX,ESTEP C--------------------------------------------------------------------- C [4] Scattering Tables C--------------------------------------------------------------------- C SUMRATE=Partial scattering rate sums, exclusive of impurity C scattering [ 1000=Energy steps; 3=# of valleys; C 15=maximum number of scattering events condidered ] C COMMON/TABLE1/SUMRATE(1000,3,15) C --------------------------------------------- C RATEMAX=Maximum rate without impurity scattering C RLAMBDA=Maximum rate with impurity scattering C COMMON/TABLE2/RATEMAX(1000,3),RLAMBDA(1000,3) C --------------------------------------------- C PIMPG=Impurity scattering tables for Gamma valley. C PIMPL=Impurity scattering tables for L valley. C PIMPX=Impurity scattering tables for X valley. C Used separately for development flexibility. C Up to 6 different doping regions included (second index). C COMMON/TABLE3/PIMPG(1000,6),PIMPL(1000,6),PIMPX(1000,6) C===================================================================== C C DIMENSION BLOCKS C C===================================================================== C [1] Electron parameters C--------------------------------------------------------------------- C K2=Square momentum C KX=x-component of the momentum C KY=y-component of the momentum C KZ=z-component of the momentum C REAL K2(ISIZE),KX(ISIZE),KY(ISIZE),KZ(ISIZE) C ---------------------------------------------- C Updated vectors for momentum components after time step C REAL KXNEW(ISIZE) C ---------------------------------------------- C K2NEW =Updated square momentum after time step C KPERP =Momentum perpendicular to field C KPERP2=Square of KPERP C REAL K2NEW(ISIZE),KPERP(ISIZE),K2PERP(ISIZE) C ------------------------------------------- C ENER=Electron kinetic energy C AE2 = coefficient function of alfa and energy C ENEW= updated energy after time step C DIMENSION ENER(ISIZE),AE2(ISIZE),ENEW(ISIZE) C ------------------------------------------- C GAM=Energy(1+Alfa*Energy) C AE1=1+Alfa*Energy C IVAL=Valley in which electrons reside C ELMASS=Mass of the individual electrons C DIMENSION GAM(ISIZE),AE1(ISIZE) DIMENSION IVAL(ISIZE) DIMENSION ELMASS(ISIZE) C ------------------------------------------- C Useful variables used for vectorization C VAR1=Hbar/mass C VAR2=q/mass C VAR3 and VAR4 are constants for flight calculations C VAR5=Hbar**2/(2*mass) C VAR6=Alfa; VAR7=4*Alfa; VAR8=1/(2*Alfa) C DIMENSION VAR1(ISIZE),VAR2(ISIZE),VAR3(ISIZE) DIMENSION VAR4(ISIZE),VAR5(ISIZE),VAR6(ISIZE) DIMENSION VAR7(ISIZE),VAR8(ISIZE) C ------------------------------------------- C ETEMP=Conventional electron temperature C VX =Drift velocity C DIMENSION ETEMP(ISIZE),VX(ISIZE) C--------------------------------------------------------------------- C [2] Flight variables C--------------------------------------------------------------------- C JFLY=Vector with indices of particles which did not scatter C after the constant time step. C TFLY=Vector with time of flight of particle. In this C approach it is set equal to the constant time step. C DIMENSION JFLY(ISIZE),TFLY(ISIZE) C ------------------------------------------ C JSCAT=Vector with indices of particles which scattered C during the constant time step. C TSCAT=Time at which the scattering occurs during the C constant time step. C T2 =Time left to complete the constant time step C after scattering. C R1 =Vector of random numbers for free flights C R2 =Modulus of total momentum used in final C state computations C DIMENSION JSCAT(ISIZE),TSCAT(ISIZE),T2(ISIZE) DIMENSION R1(ISIZE),R2(ISIZE) C--------------------------------------------------------------------- C [3] Variables to compute Maxwellian (initial equilibrium) C--------------------------------------------------------------------- DIMENSION E(0:4000),FMB(0:4000),IECOUNT(4000) DIMENSION TRAPEZ(4000),NPAR(0:4000) DIMENSION SUMMAXW(500) REAL KABS(ISIZE) C--------------------------------------------------------------------- C [4] Variables for scattering selection C--------------------------------------------------------------------- C RATEGAM=Maximum scattering rate for electrons, C during constant time step C DIMENSION RATEGAM(ISIZE) C ------------------------------------------- C RTEST=Selected (-ln r) for random flight C RNEW =Residual of RTEST during flight after time step C DIMENSION RTEST(ISIZE),RNEW(ISIZE) C ------------------------------------------- C IE= Energy histogram index for particles C IENEW= " " " after time step C DIMENSION IE(ISIZE),IENEW(ISIZE) C ------------------------------------------- C Random number for scattering type selection after free flight C (normalized to RATEGAM for every electron) C DIMENSION RSCAT(ISIZE) C ------------------------------------------- C JMECH=Vector containing the index of the selected scattering C mechanism for each electron which scatters during C time step. C JCOUNT=Vector with number of electrons which scattered C with different mechanisms during time step. C JARRAY=Maps JMECH to compute final states vith vectorization. C DIMENSION JMECH(ISIZE),JCOUNT(15),JARRAY(ISIZE,15) C ---------------------------------------------------Polar Optical C EIN=Energy at beginning of new free flight, right C after scattering event C DIMENSION EIN(1000),RAGAM(1000),AEAS(1000) DIMENSION GAMAS(1000),RAGAMAS(1000),PRAGAMA(1000) DIMENSION FACT1(1000),FACT2(1000) DIMENSION COBETA(1000),FI(1000),R(1000) DIMENSION COGAM(1000),SICO(1000) DIMENSION SISI(1000),COSETA(1000),CONFR(1000) DIMENSION SIETA(1000),SIBETA(1000) C -----------------------------------------------------Intervalley C Vectorization maps for final state C First index corresponds to valley before scattering C Second index gives the type of intervalley scat. selected C C IMAP=Map for valley after intervalley scattering C EJUMP=Energy Gap between selected valleys C EPHON=Energy of intervalley phonon C EMAP=Effective mass map, set up already for expansion to C device with various material regions. C AMAP2=Map for alfa C DIMENSION IMAP(3,9),EJUMP(3,9) DIMENSION EPHON(3,9) DIMENSION EMAP(3),AMAP2(3) C--------------------------------------------------------------------- C INITIALIZE RANDOM NUMBER GENERATOR C--------------------------------------------------------------------- double precision seed SEED=1.45412342882337 DO 100 I=1,25 100 RANINIT=RAN(SEED) C C===================================================================== C Open file statements for use on APOLLO workstation C User needs to change for specific machine C C Unit 11 Final distribution file C Unit 12 Written output file C Unit 13 Input data file C OPEN(UNIT=11,FILE='final') OPEN(UNIT=12,FILE='bulkout') OPEN(UNIT=13,FILE='bulkin') C===================================================================== C INPUTS FROM FILE C===================================================================== C Input electric field (in kV/cm) C READ(13,*) FINPUT C--------------------------------------------------------------------- C Input maximum number of iterations (time steps) C READ(13,*) ITERMAX C--------------------------------------------------------------------- C Select initial condition C C JREAD=0 ====> Start from equilibrium initial condition C JREAD=1 ====> Start from final distribution of previous run C C JYES=0 ====> Cumulative averages of previous simulations C are not taken into account C JYES=1 ====> Previous averages are taken into account C (If JREAD=0, we set automatically JYES=0) C READ(13,*) JREAD READ(13,*) JYES C IF (JREAD.EQ.0) JYES=0 C===================================================================== C [ DATA 1 ] DATA BLOCK FOR UNIVERSAL CONSTANTS C--------------------------------------------------------------------- DATA BOLTZ/1.3807E-23/ DATA EM0/9.1095E-31/ DATA EP0/8.85E-12/ DATA HBAR/1.05459E-34/ DATA Q/1.60219E-19/ PI=ACOS(-1.) TWOPI=2.*PI C--------------------------------------------------------------------- C [ DATA 2 ] ALGORITHM RELATED INPUTS C--------------------------------------------------------------------- C Number of particles in simulation DATA NELEC/5000/ C--------------------------------------------------------------------- C Histogram intervals in scattering table DATA NESTEP/1000/ C--------------------------------------------------------------------- C EMAX=Highest energy in table (input in Kelvin) C ESTEP=Energy increment in table (in Kelvin) C STEPJ=Energy increment (in Joule) C DTFL=Time step for Constant Time Technique C DATA EMAX/10000./ ESTEP=EMAX/FLOAT(NESTEP) STEPJ=ESTEP*BOLTZ DATA DTFL/1.E-15/ C--------------------------------------------------------------------- C [ DATA 3 ] DOPING AND TEMPERATURE C--------------------------------------------------------------------- C DOP(i)=Doping level in various regions (ready for device simul.) C TEMP=Lattice temperature C DATA DOP(1),DOP(2),DOP(3)/1.E18,1.E18,.5E18/ DATA DOP(4),DOP(5)/1.E10,1.E16/ DATA TEMP/300./ C--------------------------------------------------------------------- C [ DATA 4 ] MATERIAL PARAMETERS C--------------------------------------------------------------------- C--------------------------------------------------------------------- C Gallium Arsenide C--------------------------------------------------------------------- C DATA ZL,ZX/4.,3./ DATA EMSTAR1(1),EMSTAR2(1),EMSTAR3(1)/.063,.222,.58/ DATA GAP(1),GAPGL(1),GAPGX(1)/1.42,.33,.522/ ALFA1(1)=(1.-EMSTAR1(1))**2/GAP(1) ALFA2(1)=ALFA1(1)*461./610. ALFA3(1)=ALFA1(1)*204./610. DATA DEFGL(1),DEFGX(1),DEFLX(1)/1.E9,1.E9,5.E8/ DATA DEFLL(1),DEFXX(1)/1.E9,7.E8/ DATA DEFAC(1),RHO(1),VS(1)/7.,5.37,5.22E5/ DATA EPZER(1),EPINF(1)/12.9,10.92/ DATA TPOP(1)/410./ TINTGL(1)=323. TINTGX(1)=347. TINTLX(1)=340. TEQINL(1)=336. TEQINX(1)=347. C--------------------------------------------------------------------- C J=1 Select GaAs for the simulation C In device simulation J=2 will select AlGaAs C J=1 C--------------------------------------------------------------------- C Calculate Thermal Voltage and Intervalley Energies C--------------------------------------------------------------------- C TK=BOLTZ*TEMP TKINT1(J)=BOLTZ*TINTGL(J) TKINT2(J)=BOLTZ*TINTGX(J) TKINT3(J)=BOLTZ*TINTLX(J) TKEIN2(J)=BOLTZ*TEQINL(J) TKEIN3(J)=BOLTZ*TEQINX(J) TKPOP(J)=BOLTZ*TPOP(J) C--------------------------------------------------------------------- C Useful Constants for flight calculations C--------------------------------------------------------------------- CON1=Q/HBAR CON2=-CON1*DTFL C--------------------------------------------------------------------- C Convert units to mks C--------------------------------------------------------------------- C Doping in cubic meters (loop has been unrolled) C DOP(1)=DOP(1)*1.0E6 DOP(2)=DOP(2)*1.0E6 DOP(3)=DOP(3)*1.0E6 DOP(4)=DOP(4)*1.0E6 DOP(5)=DOP(5)*1.0E6 C ------------------------------ C Gallium Arsenide parameters J=1 ALFA1(J)=ALFA1(J)/Q ALFA2(J)=ALFA2(J)/Q ALFA3(J)=ALFA3(J)/Q VS(J)=VS(J)/100. RHO(J)=RHO(J)*1000. DEFGL(J)=DEFGL(J)*Q*100. DEFGX(J)=DEFGX(J)*Q*100. DEFLX(J)=DEFLX(J)*Q*100. DEFLL(J)=DEFLL(J)*Q*100. DEFXX(J)=DEFXX(J)*Q*100. DEFAC(J)=DEFAC(J)*Q FX=FINPUT*1000.*100. EPZER(J)=EPZER(J)*EP0 EPINF(J)=EPINF(J)*EP0 GAPGL(J)=GAPGL(J)*Q GAPGX(J)=GAPGX(J)*Q GAP(J)=GAP(J)*Q C--------------------------------------------------------------------- EM1(J)=EMSTAR1(J)*EM0 EM2(J)=EMSTAR2(J)*EM0 EM3(J)=EMSTAR3(J)*EM0 CC(1,J)=2.*(EM1(J)/HBAR)/HBAR CC(2,J)=2.*(EM2(J)/HBAR)/HBAR CC(3,J)=2.*(EM3(J)/HBAR)/HBAR TWOA1(J)=2.*ALFA1(J) TWOA2(J)=2.*ALFA2(J) TWOA3(J)=2.*ALFA3(J) FOURA1(J)=4.*ALFA1(J) FOURA2(J)=4.*ALFA2(J) FOURA3(J)=4.*ALFA3(J) TAINV1(J)=1./(2.*ALFA1(J)) TAINV2(J)=1./(2.*ALFA2(J)) TAINV3(J)=1./(2.*ALFA3(J)) C--------------------------------------------------------------------- C [ DATA 5 ] Create map for intervalley scattering C--------------------------------------------------------------------- C IMAP(i,j) gives the FINAL valley after an intervalley scattering C The intervalley scatterings are grouped according to the index j C C i = Valley before scattering C C Non-equivalent intervalley transitions C C i = 1 j = 4 Absorption I - Gamma to L C j = 5 Emission I - Gamma to L C j = 6 Absorption II - Gamma to X C j = 7 Emission II - Gamma to X C i = 2 j = 4 Absorption I - L to Gamma C j = 5 Emission I - L to Gamma C j = 6 Absorption II - L to X C j = 7 Emission II - L to X C i = 3 j = 4 Absorption I - X to Gamma C j = 5 Emission I - X to Gamma C j = 6 Absorption II - X to L C j = 7 Emission II - X to L C C Equivalent intervalley transitions C C i = 2 j = 8 Absorption - L to L C j = 9 Emission - L to L C i = 3 j = 8 Absorption - X to X C j = 9 Emission - X to X C DATA IMAP(1,4),IMAP(1,5),IMAP(1,6),IMAP(1,7)/2,2,3,3/ DATA IMAP(2,4),IMAP(2,5),IMAP(2,6),IMAP(2,7)/1,1,3,3/ DATA IMAP(2,8),IMAP(2,9)/2,2/ DATA IMAP(3,4),IMAP(3,5),IMAP(3,6),IMAP(3,7)/1,1,2,2/ DATA IMAP(3,8),IMAP(3,9)/3,3/ C--------------------------------------------------------------------- C Energy change in jump between valleys C--------------------------------------------------------------------- C EJUMP(i,j) maps the change in energy gap for non-equivalent C intervalley transitions. The indices (i,j) have the same meaning C as in the table given for IMAP(i,j). C EJUMP(1,4)=-GAPGL(1) EJUMP(1,5)=-GAPGL(1) EJUMP(1,6)=-GAPGX(1) EJUMP(1,7)=-GAPGX(1) EJUMP(2,4)=GAPGL(1) EJUMP(2,5)=GAPGL(1) EJUMP(2,6)=GAPGL(1)-GAPGX(1) EJUMP(2,7)=GAPGL(1)-GAPGX(1) EJUMP(3,4)=GAPGX(1) EJUMP(3,5)=GAPGX(1) EJUMP(3,6)=GAPGX(1)-GAPGL(1) EJUMP(3,7)=GAPGX(1)-GAPGL(1) C--------------------------------------------------------------------- C Intervalley phonon energy map C--------------------------------------------------------------------- C EPHON(i,j) maps the energy change due to phonon absorption and C emission for non-equivalent intervalley transitions. The indices C (i,j) have the same meaning as in the table given for IMAP(i,j). C EPHON(1,4)=TKINT1(1) EPHON(1,5)=TKINT1(1) EPHON(1,6)=TKINT2(1) EPHON(1,7)=TKINT2(1) EPHON(2,4)=TKINT1(1) EPHON(2,5)=TKINT1(1) EPHON(2,6)=TKINT3(1) EPHON(2,7)=TKINT3(1) EPHON(3,4)=TKINT2(1) EPHON(3,5)=TKINT2(1) EPHON(3,6)=TKINT3(1) EPHON(3,7)=TKINT3(1) EPHON(2,8)=TKEIN2(1) EPHON(2,9)=TKEIN2(1) EPHON(3,8)=TKEIN3(1) EPHON(3,9)=TKEIN3(1) C--------------------------------------------------------------------- C Non-parabolicity constant map C--------------------------------------------------------------------- AMAP2(1)=ALFA1(1) AMAP2(2)=ALFA2(1) AMAP2(3)=ALFA3(1) C--------------------------------------------------------------------- C Impurity scattering parameter C--------------------------------------------------------------------- SCMAP=4./((Q/BOLTZ)*(Q/EPZER(1))*DOP(5)/TEMP) C--------------------------------------------------------------------- C Time of flight for Modified Constant Time iteration C--------------------------------------------------------------------- DO 120 I=1,NELEC TFLY(I)=DTFL 120 CONTINUE C--------------------------------------------------------------------- C Mass by valley map C--------------------------------------------------------------------- EMAP(1)=EM1(1) EMAP(2)=EM2(1) EMAP(3)=EM3(1) C--------------------------------------------------------------------- C Initialize global scattering counters C--------------------------------------------------------------------- C C CREATE SCATTERING TABLES IN SUBROUTINE SCAT C C--------------------------------------------------------------------- C J=1 selects GaAs J=1 CALL SCAT(J) C--------------------------------------------------------------------- C INITIAL CONDITION C--------------------------------------------------------------------- C If JREAD = 1, read from file C If JREAD = 0, create Maxwell-Boltzmann C IF (JREAD.EQ.1) THEN DO 130 I=1,NELEC READ(11,*)IVAL(I),KX(I),KY(I),KZ(I) 130 CONTINUE C IF (JYES.EQ.1) THEN READ(11,*) AVELOOLD,AVENEOLD,ITERTOT ELSE AVELOOLD=0.0 AVENEOLD=0.0 ITERTOT=0 ENDIF C ---------------------------------------------------------------- C Create necessary dimensioned variables C DO 140 I=1,NELEC ELMASS(I)=EMAP(IVAL(I)) VAR1(I)=HBAR/ELMASS(I) VAR2(I)=Q/ELMASS(I) VAR3(I)=-.5*VAR2(I)*DTFL**2 VAR4(I)=VAR1(I)*DTFL VAR5(I)=VAR1(I)*HBAR*.5 VAR6(I)=AMAP2(IVAL(I)) VAR7(I)=4.*VAR6(I) VAR8(I)=.5/VAR6(I) K2PERP(I)=KY(I)**2+KZ(I)**2 KPERP(I)=SQRT(K2PERP(I)) K2NEW(I)=KX(I)**2+K2PERP(I) GAM(I)=VAR5(I)*K2NEW(I) ENER(I)=(SQRT(1.+VAR7(I)*GAM(I))-1.)*VAR8(I) IE(I)=IFIX(ENER(I)/STEPJ)+1 IF(IE(I).GE.NESTEP) IE(I)=NESTEP 140 CONTINUE C--------------------------------------------------------------------- ELSE C--------------------------------------------------------------------- C C SET UP INITIAL MAXWELLIAN DISTRIBUTION C C--------------------------------------------------------------------- C Generates the vector NPAR(J) which holds the number C of electrons at energy J meV's assuming an initial C distribution of 1000 (NSTAND) particles. To obtain C a Maxwell-Boltzmann distribution for NELEC particles C where NELEC is a multiple of 1000, simply multiply C NPAR(J) by NELEC/NSTAND. This subroutine uses the C trapezoidal rule in piecewise integration of the C Maxwell-Boltzmann distribution. C--------------------------------------------------------------------- C TKEV = kT in [eV] C NSTAND = Elementary ensemble of electrons generated = 1,000 C NEMX = Number of Energy steps considered in the integral C DE = Energy step in [eV] C ROUND = Rounding factor to adjust to integer number of electrons C NSCALE = Factor which adjusts to actual number of electrons C selected (should be multiple of NSTAND) C NPAR(j)= Histogram of electrons obtained for the distribution C (j indicates the energy interval) C TKEV=TK/Q NSTAND=1000 NEMX=145 DE=.001 ROUND=.5 NSCALE=NELEC/NSTAND DO 160 I=1,500 SUMMAXW(I)=0.0 160 CONTINUE C--------------------------------------------------------------------- C Calculate the discretized energies C and the unnormalized distribution C--------------------------------------------------------------------- DO 170 I=0,NEMX E(I)=FLOAT(I)*DE FMB(I)=SQRT(E(I))*EXP(-E(I)/TKEV) 170 CONTINUE C--------------------------------------------------------------------- C Compute the area of each discrete trapezoid C--------------------------------------------------------------------- DO 180 I=1,NEMX TRAPEZ(I)=(FMB(I)+FMB(I-1))*DE/2. 180 CONTINUE C--------------------------------------------------------------------- C Sum up all the piecewise areas to obtain the normalization C Note that the loop DO 190 is not vectorizable C--------------------------------------------------------------------- SUMMAXW(1)=TRAPEZ(1) DO 190 I=2,NEMX SUMMAXW(I)=SUMMAXW(I-1)+TRAPEZ(I) 190 CONTINUE RNORM=SUMMAXW(NEMX) C--------------------------------------------------------------------- C Obtain the number of particles at each energy C by dividing by the normalization factor and C multiplying by NSTAND electrons. An empirical C rounding constant (ROUND) is added to assure C that the total number of electrons is NSTAND. C Scale distribution up to NELEC particles C--------------------------------------------------------------------- CTEMP=FLOAT(NSTAND)/RNORM DO 200 I=1,NEMX NPAR(I)=IFIX(TRAPEZ(I)*CTEMP+ROUND)*NSCALE 200 CONTINUE C--------------------------------------------------------------------- C Compute energy of each particle. These C values are stored in the array ENER(J) C * * * NOT vectorized C--------------------------------------------------------------------- KOUNT=0 C--------------------------------------------------------------------- C The scaling factor SCAL yields the correct average temperature C SCAL=3.61 DO 220 I=1,NEMX DO 230 J=1,NPAR(I) KOUNT=KOUNT+1 ENER(KOUNT)=DE*(FLOAT(I-1)+RAN(SEED)*SCAL) 230 CONTINUE 220 CONTINUE C--------------------------------------------------------------------- C Convert the energy variable ENER to MKS units DO 240 I=1,NELEC ENER(I)=ENER(I)*Q C--------------------------------------------------------------------- C Create momentum distribution (3-D) C--------------------------------------------------------------------- IE(I)=IFIX(ENER(I)/STEPJ)+1 AE1(I)=1.+ALFA1(1)*ENER(I) AE2(I)=2.*AE1(I)-1. GAM(I)=ENER(I)*AE1(I) K2(I)=CC(1,1)*GAM(I) KABS(I)=SQRT(K2(I)) 240 CONTINUE DO 241 I=1,NELEC KX(I)=KABS(I)*(2.*RAN(SEED)-1.) 241 CONTINUE DO 242 I=1,NELEC KPERP(I)=SQRT(K2(I)-KX(I)**2) K2PERP(I)=KPERP(I)**2 242 CONTINUE DO 243 I=1,NELEC KZ(I)=KPERP(I)*(2.*RAN(SEED)-1.) 243 CONTINUE DO 244 I=1,NELEC KY(I)=SQRT(K2PERP(I)-KZ(I)**2) 244 CONTINUE DO 245 I=1,NELEC KY(I)=SIGN(KY(I),RAN(SEED)-.5) 245 CONTINUE C ---------------------------------------------------------------- C Assign each electron to Gamma valley C Originally ELMASS(I)=EMAP(IVAL(I)) DO 246 I=1,NELEC IVAL(I)=1 ELMASS(I)=EMAP(1) VAR1(I)=HBAR/ELMASS(I) VAR2(I)=Q/ELMASS(I) VAR3(I)=-.5*VAR2(I)*DTFL**2 VAR4(I)=VAR1(I)*DTFL VAR5(I)=VAR1(I)*HBAR*.5 VAR6(I)=ALFA1(1) VAR7(I)=4.*ALFA1(1) VAR8(I)=.5/ALFA1(1) 246 CONTINUE ENDIF ENERSUM=0. DO 269 I=1,NELEC ENERSUM=ENERSUM+ENER(I) 269 CONTINUE TEMP=.66666667*ENERSUM/NELEC/BOLTZ C--------------------------------------------------------------------- c write(6,*)'start mc loop' C===================================================================== C C ELECTRON MONTE CARLO LOOP C C===================================================================== C avelo is the variable used to evaluate the drift velocity C avene is the variable used to evaluate the average energy C AVELO=0.0 AVENE=0.0 C--------------------------------------------------------------------- C Initialize RTEST(I) for random flights in 1st iteration C--------------------------------------------------------------------- DO 300 I=1,NELEC RTEST(I)=-ALOG(RAN(SEED)) 300 CONTINUE C--------------------------------------------------------------------- C Main time loop C--------------------------------------------------------------------- DO 1000 ITER=1,ITERMAX C ---------------------------------------------------------------- C Free flights C ---------------------------------------------------------------- C Update momentum DO 2100 I=1,NELEC KXNEW(I)=KX(I)+CON2*FX C ---------------------------------------------------------------- C Update energy and index K2NEW(I)=KXNEW(I)**2+K2PERP(I) GAM(I)=VAR1(I)*K2NEW(I)*HBAR*.5 ENEW(I)=(SQRT(1.+VAR7(I)*GAM(I))-1.)*VAR8(I) IENEW(I)=IFIX(ENEW(I)/STEPJ)+1 IF (IENEW(I).GE.NESTEP) IENEW(I)=NESTEP RATEGAM(I)=AMAX1(RLAMBDA(IE(I),IVAL(I)), # RLAMBDA(IENEW(I),IVAL(I))) C--------------------------------------------------------------------- C Determine residual portion of random flight variable C RNEW(I)=RTEST(I)-RATEGAM(I)*DTFL 2100 CONTINUE C --------------------------------------------- NFLY=0 NSCAT=0 C --------------------------------------------- C If RNEW(i) is equal or less than 0, scattering must occurr C during the current subflight. C C NFLY = Total # of particles which did not scatter C JFLY = Vector with indices of particles which did not scatter C NSCAT = Total # of particles which scattered in subflight C JSCAT = Vector with indices of particles which scattered C DO 2180 I=1,NELEC IF (RNEW(I).GE.0.) THEN NFLY=NFLY+1 JFLY(NFLY)=I ELSE NSCAT=NSCAT+1 JSCAT(NSCAT)=I ENDIF 2180 CONTINUE C--------------------------------------------------------------------- C Treat particles that do not scatter C--------------------------------------------------------------------- C DO 4010 I=1,NFLY RTEST(JFLY(I))=RNEW(JFLY(I)) KX(JFLY(I))=KXNEW(JFLY(I)) ENER(JFLY(I))=ENEW(JFLY(I)) IE(JFLY(I))=IENEW(JFLY(I)) 4010 CONTINUE C--------------------------------------------------------------------- C TREAT PARTICLES THAT DO SCATTER C--------------------------------------------------------------------- C TSCAT(i) = Time of scattering during the subflight C T2(i) = Time left to fly during subflight C DO 5010 I=1,NSCAT TSCAT(I)=RTEST(JSCAT(I))/RATEGAM(JSCAT(I)) T2(I)=DTFL-TSCAT(I) 5010 CONTINUE C--------------------------------------------------------------------- C Recalculate variables at time of scattering C--------------------------------------------------------------------- DO 5030 I=1,NSCAT KXNEW(JSCAT(I))=KX(JSCAT(I))-CON1*FX*TSCAT(I) K2NEW(JSCAT(I))=KXNEW(JSCAT(I))**2+K2PERP(JSCAT(I)) GAM(JSCAT(I))=VAR1(JSCAT(I))*K2NEW(JSCAT(I))*HBAR*.5 ENEW(JSCAT(I))=(SQRT(1.+VAR7(JSCAT(I))*GAM(JSCAT(I)))-1.)* # VAR8(JSCAT(I)) IENEW(JSCAT(I))=IFIX(ENEW(JSCAT(I))/STEPJ)+1 C ---------------------------------------------------------------- C NOTE: If maximum energy index exceed, IENEW is set equal to C the maximum C IF (IENEW(JSCAT(I)).GE.NESTEP) IENEW(JSCAT(I))=NESTEP 5030 CONTINUE C===================================================================== C SCATTERING CHOICE C===================================================================== C Create vector of random numbers for scattering C DO 5090 I=1,NSCAT R1(I)=-ALOG(RAN(SEED)) 5090 CONTINUE C--------------------------------------------------------------------- C Generate random variable, scaled with the the total C scattering rate for the energy of the particle C DO 5110 I=1,NSCAT RSCAT(I)=RAN(SEED)*RATEGAM(JSCAT(I)) 5110 CONTINUE C--------------------------------------------------------------------- C This loop select the scattering mechanism for each C particle, by comparing the random variable RSCAT(i) C to the partial scattering rate sums SUMRATE(j,k,l) C C JMECH(I)=Vector containing the scattering mechanism C DO 5120 I=1,NSCAT JMECH(I)=1 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),1)) # JMECH(I)=2 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),2)) # JMECH(I)=3 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),3)) # JMECH(I)=4 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),4)) # JMECH(I)=5 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),5)) # JMECH(I)=6 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),6)) # JMECH(I)=7 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),7)) # JMECH(I)=8 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),8)) # JMECH(I)=9 IF (RSCAT(I).GT.SUMRATE(IENEW(JSCAT(I)),IVAL(JSCAT(I)),9)) # JMECH(I)=10 IF (RSCAT(I).GT.RLAMBDA(IENEW(JSCAT(I)),IVAL(JSCAT(I)))) # JMECH(I)=11 5120 CONTINUE C ---------------------------------------------------------------- C Reset counter DO 5130 I=1,15 JCOUNT(I)=0 5130 CONTINUE C ---------------------------------------------------------------- C Fill JARRAY(J,K) (Not vectorizable) C JARRAY(J,K)=Vector which groups electrons by scattering type C JMECH(I) =Vector containing the scattering mechanism C DO 5140 I=1,NSCAT JCOUNT(JMECH(I))=JCOUNT(JMECH(I))+1 JARRAY(JCOUNT(JMECH(I)),JMECH(I))=JSCAT(I) 5140 CONTINUE c do 5555 i=1,15 c5555 write(6,*) i,jcount(i) C===================================================================== C FINAL STATE C C The variables with the suffix NEW are the values just before the C scattering occurs. C The indirect addressing of the arrays may cause vectorizing C problems. Sometimes, it is possible to vecorize a loop by C moving some of the loops out C===================================================================== C [ 1 ] Polar optical absorption C--------------------------------------------------------------------- DO 5210 I=1,JCOUNT(1) EIN(I)=ENEW(JARRAY(I,1))+TKPOP(1) RAGAM(I)=SQRT(GAM(JARRAY(I,1))) AEAS(I)=VAR6(JARRAY(I,1))*EIN(I) GAMAS(I)=EIN(I)*(1.+AEAS(I)) RAGAMAS(I)=SQRT(GAMAS(I)) PRAGAMA(I)=RAGAM(I)*RAGAMAS(I) FACT1(I)=(RAGAM(I)+RAGAMAS(I))/(RAGAMAS(I)-RAGAM(I)) FACT2(I)=AEAS(I)*ENEW(JARRAY(I,1)) 5210 CONTINUE C --------------------------------------------- C Rejection for angle (Not vectorized) C --------------------------------------------- DO 5270 I=1,JCOUNT(1) c write (6,*) 'i ',i,pragama(i) 5260 CONTINUE COBETA(I)=(GAM(JARRAY(I,1))+GAMAS(I)- # ((RAGAMAS(I)-RAGAM(I))*FACT1(I)**RAN(SEED))**2)/(2.*PRAGAMA(I)) CONFR(I)=((PRAGAMA(I)+FACT2(I)*COBETA(I))/ # (PRAGAMA(I)+FACT2(I)))**2 IF(RAN(SEED).GT.CONFR(I)) GO TO 5260 5270 CONTINUE DO 5280 I=1,JCOUNT(1) FI(I)=TWOPI*RAN(SEED) 5280 CONTINUE DO 5290 I=1,JCOUNT(1) R(I)=SQRT(K2NEW(JARRAY(I,1))) COSETA(I)=KY(JARRAY(I,1))/KPERP(JARRAY(I,1)) SIETA(I)=KZ(JARRAY(I,1))/KPERP(JARRAY(I,1)) COGAM(I)=KXNEW(JARRAY(I,1))/R(I) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- K2NEW(JARRAY(I,1))=GAMAS(I)*CC(IVAL(JARRAY(I,1)),1) R2(I)=SQRT(K2NEW(JARRAY(I,1))) COBETA(I)=AMIN1(COBETA(I),1.) SICO(I)=SQRT(1.-COBETA(I)**2+1.E-37)*COS(FI(I)) SISI(I)=SQRT(1.-COBETA(I)**2+1.E-37)*SIN(FI(I)) 5290 CONTINUE DO 5291 I=1,JCOUNT(1) KZ(JARRAY(I,1))=R2(I)*(COBETA(I)*KZ(JARRAY(I,1))/R(I)- # SICO(I)*COGAM(I)*SIETA(I)+SISI(I)*COSETA(I)) KY(JARRAY(I,1))=R2(I)*(COBETA(I)*KY(JARRAY(I,1))/R(I)- # SICO(I)*COGAM(I)*COSETA(I)-SISI(I)*SIETA(I)) KX(JARRAY(I,1))=R2(I)*(COBETA(I)*COGAM(I)+SICO(I)* # SQRT(1.-COGAM(I)**2+1.E-37)) ENER(JARRAY(I,1))=EIN(I) 5291 CONTINUE C--------------------------------------------------------------------- C [ 2 ] Polar optical emission C--------------------------------------------------------------------- DO 5410 I=1,JCOUNT(2) EIN(I)=AMAX1(ENEW(JARRAY(I,2))-TKPOP(1),1.E-24) RAGAM(I)=SQRT(GAM(JARRAY(I,2))) AEAS(I)=VAR6(JARRAY(I,2))*EIN(I) GAMAS(I)=EIN(I)*(1.+AEAS(I)) RAGAMAS(I)=SQRT(GAMAS(I)) PRAGAMA(I)=RAGAM(I)*RAGAMAS(I) FACT1(I)=(RAGAM(I)+RAGAMAS(I))/(RAGAM(I)-RAGAMAS(I)) FACT2(I)=AEAS(I)*ENEW(JARRAY(I,2)) 5410 CONTINUE C --------------------------------------------- C Rejection for angle (Not vectorized) C --------------------------------------------- DO 5470 I=1,JCOUNT(2) c write (6,*) 'i , jcount(i)',i,pragama(i) 5460 CONTINUE COBETA(I)=(GAM(JARRAY(I,2))+GAMAS(I)- # ((RAGAM(I)-RAGAMAS(I))*FACT1(I)**RAN(SEED))**2) # /(2.*PRAGAMA(I)) CONFR(I)=((PRAGAMA(I)+FACT2(I)*COBETA(I))/ # (PRAGAMA(I)+FACT2(I)))**2 IF(RAN(SEED).GT.CONFR(I)) GO TO 5460 5470 CONTINUE DO 5480 I=1,JCOUNT(2) FI(I)=TWOPI*RAN(SEED) 5480 CONTINUE DO 5490 I=1,JCOUNT(2) R(I)=SQRT(K2NEW(JARRAY(I,2))) COSETA(I)=KY(JARRAY(I,2))/KPERP(JARRAY(I,2)) SIETA(I)=KZ(JARRAY(I,2))/KPERP(JARRAY(I,2)) COGAM(I)=KXNEW(JARRAY(I,2))/R(I) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- K2NEW(JARRAY(I,2))=GAMAS(I)*CC(IVAL(JARRAY(I,2)),1) R2(I)=SQRT(K2NEW(JARRAY(I,2))) COBETA(I)=AMIN1(COBETA(I),1.) SICO(I)=SQRT(1.-COBETA(I)**2+1.E-37)*COS(FI(I)) SISI(I)=SQRT(1.-COBETA(I)**2+1.E-37)*SIN(FI(I)) 5490 CONTINUE DO 5491 I=1,JCOUNT(2) KZ(JARRAY(I,2))=R2(I)*(COBETA(I)*KZ(JARRAY(I,2))/R(I)- # SICO(I)*COGAM(I)*SIETA(I)+SISI(I)*COSETA(I)) KY(JARRAY(I,2))=R2(I)*(COBETA(I)*KY(JARRAY(I,2))/R(I)- # SICO(I)*COGAM(I)*COSETA(I)-SISI(I)*SIETA(I)) KX(JARRAY(I,2))=R2(I)*(COBETA(I)*COGAM(I)+SICO(I)* # SQRT(1.-COGAM(I)**2+1.E-37)) ENER(JARRAY(I,2))=EIN(I) 5491 CONTINUE C--------------------------------------------------------------------- C [ 3 ] Acoustic Phonon Scattering C--------------------------------------------------------------------- C Rejection for angle (Not vectorized) C --------------------------------------------- DO 5610 I=1,JCOUNT(3) 5620 CONTINUE COBETA(I)=2.*RAN(SEED)-1. CONFR(I)=((1.+COBETA(I))*ENEW(JARRAY(I,3))* # VAR6(JARRAY(I,3))+1.)**2 IF(RAN(SEED).GT.CONFR(I)) GO TO 5620 5610 CONTINUE DO 5630 I=1,JCOUNT(3) FI(I)=TWOPI*RAN(SEED) 5630 CONTINUE DO 5631 I=1,JCOUNT(3) R(I)=SQRT(K2NEW(JARRAY(I,3))) COSETA(I)=KY(JARRAY(I,3))/KPERP(JARRAY(I,3)) SIETA(I)=KZ(JARRAY(I,3))/KPERP(JARRAY(I,3)) SIBETA(I)=SQRT(1.-COBETA(I)**2) SICO(I)=SIBETA(I)*COS(FI(I)) SISI(I)=SIBETA(I)*SIN(FI(I)) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,3))=KXNEW(JARRAY(I,3))*COBETA(I)+ # KPERP(JARRAY(I,3))*SICO(I) KZ(JARRAY(I,3))=KZ(JARRAY(I,3))*COBETA(I)-KXNEW(JARRAY(I,3))* # SICO(I)*SIETA(I)+R(I)*SISI(I)*COSETA(I) KY(JARRAY(I,3))=KY(JARRAY(I,3))*COBETA(I)-KXNEW(JARRAY(I,3))* # SICO(I)*COSETA(I)-R(I)*SISI(I)*SIETA(I) ENER(JARRAY(I,3))=ENEW(JARRAY(I,3)) 5631 CONTINUE C--------------------------------------------------------------------- C [ 4 ] Non-Equivalent Intervalley Absorption I C--------------------------------------------------------------------- C Intervalley Absorption (non-eq.) GL,LG,XG DO 5710 I=1,JCOUNT(4) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5710 CONTINUE DO 5711 I=1,JCOUNT(4) EIN(I)=ENEW(JARRAY(I,4))+EJUMP(IVAL(JARRAY(I,4)),4) # +EPHON(IVAL(JARRAY(I,4)),4) EIN(I)=AMAX1(EIN(I),1.E-24) IVAL(JARRAY(I,4))=IMAP(IVAL(JARRAY(I,4)),4) ELMASS(JARRAY(I,4))=EMAP(IVAL(JARRAY(I,4))) VAR1(JARRAY(I,4))=HBAR/ELMASS(JARRAY(I,4)) VAR2(JARRAY(I,4))=Q/ELMASS(JARRAY(I,4)) VAR3(JARRAY(I,4))=-.5*VAR2(JARRAY(I,4))*DTFL**2 VAR4(JARRAY(I,4))=VAR1(JARRAY(I,4))*DTFL VAR5(JARRAY(I,4))=VAR1(JARRAY(I,4))*HBAR*.5 VAR6(JARRAY(I,4))=AMAP2(IVAL(JARRAY(I,4))) VAR7(JARRAY(I,4))=4.*VAR6(JARRAY(I,4)) VAR8(JARRAY(I,4))=.5/VAR6(JARRAY(I,4)) AEAS(I)=VAR6(JARRAY(I,4))*EIN(I) GAM(JARRAY(I,4))=EIN(I)*(1.+AEAS(I)) R(I)=SQRT(GAM(JARRAY(I,4))*CC(IVAL(JARRAY(I,4)),1)) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,4))=R(I)*COBETA(I) KPERP(JARRAY(I,4))=SQRT(R(I)**2-KX(JARRAY(I,4))**2) KY(JARRAY(I,4))=KPERP(JARRAY(I,4))*COS(FI(I)) KZ(JARRAY(I,4))=KPERP(JARRAY(I,4))*SIN(FI(I)) ENER(JARRAY(I,4))=EIN(I) 5711 CONTINUE C--------------------------------------------------------------------- C [ 5 ] Non-Equivalent Intervalley Emission I C--------------------------------------------------------------------- C Intervalley Emission (non-eq.) GL,LG,XG DO 5720 I=1,JCOUNT(5) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5720 CONTINUE DO 5721 I=1,JCOUNT(5) EIN(I)=ENEW(JARRAY(I,5))+EJUMP(IVAL(JARRAY(I,5)),5) # -EPHON(IVAL(JARRAY(I,5)),5) EIN(I)=AMAX1(EIN(I),1.E-24) IVAL(JARRAY(I,5))=IMAP(IVAL(JARRAY(I,5)),5) ELMASS(JARRAY(I,5))=EMAP(IVAL(JARRAY(I,5))) VAR1(JARRAY(I,5))=HBAR/ELMASS(JARRAY(I,5)) VAR2(JARRAY(I,5))=Q/ELMASS(JARRAY(I,5)) VAR3(JARRAY(I,5))=-.5*VAR2(JARRAY(I,5))*DTFL**2 VAR4(JARRAY(I,5))=VAR1(JARRAY(I,5))*DTFL VAR5(JARRAY(I,5))=VAR1(JARRAY(I,5))*HBAR*.5 VAR6(JARRAY(I,5))=AMAP2(IVAL(JARRAY(I,5))) VAR7(JARRAY(I,5))=4.*VAR6(JARRAY(I,5)) VAR8(JARRAY(I,5))=.5/VAR6(JARRAY(I,5)) AEAS(I)=VAR6(JARRAY(I,5))*EIN(I) GAM(JARRAY(I,5))=EIN(I)*(1.+AEAS(I)) R(I)=SQRT(GAM(JARRAY(I,5))*CC(IVAL(JARRAY(I,5)),1)) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,5))=R(I)*COBETA(I) KPERP(JARRAY(I,5))=SQRT(R(I)**2-KX(JARRAY(I,5))**2) KY(JARRAY(I,5))=KPERP(JARRAY(I,5))*COS(FI(I)) KZ(JARRAY(I,5))=KPERP(JARRAY(I,5))*SIN(FI(I)) ENER(JARRAY(I,5))=EIN(I) 5721 CONTINUE C--------------------------------------------------------------------- C [ 6 ] Non-Equivalent Intervalley Absorption II C--------------------------------------------------------------------- C Intervalley Absorption (non-eq.) GX,LX,XL DO 5730 I=1,JCOUNT(6) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5730 CONTINUE DO 5731 I=1,JCOUNT(6) EIN(I)=ENEW(JARRAY(I,6))+EJUMP(IVAL(JARRAY(I,6)),6) # +EPHON(IVAL(JARRAY(I,6)),6) EIN(I)=AMAX1(EIN(I),1.E-24) IVAL(JARRAY(I,6))=IMAP(IVAL(JARRAY(I,6)),6) ELMASS(JARRAY(I,6))=EMAP(IVAL(JARRAY(I,6))) VAR1(JARRAY(I,6))=HBAR/ELMASS(JARRAY(I,6)) VAR2(JARRAY(I,6))=Q/ELMASS(JARRAY(I,6)) VAR3(JARRAY(I,6))=-.5*VAR2(JARRAY(I,6))*DTFL**2 VAR4(JARRAY(I,6))=VAR1(JARRAY(I,6))*DTFL VAR5(JARRAY(I,6))=VAR1(JARRAY(I,6))*HBAR*.5 VAR6(JARRAY(I,6))=AMAP2(IVAL(JARRAY(I,6))) VAR7(JARRAY(I,6))=4.*VAR6(JARRAY(I,6)) VAR8(JARRAY(I,6))=.5/VAR6(JARRAY(I,6)) AEAS(I)=VAR6(JARRAY(I,6))*EIN(I) GAM(JARRAY(I,6))=EIN(I)*(1.+AEAS(I)) R(I)=SQRT(GAM(JARRAY(I,6))*CC(IVAL(JARRAY(I,6)),1)) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,6))=R(I)*COBETA(I) KPERP(JARRAY(I,6))=SQRT(R(I)**2-KX(JARRAY(I,6))**2) KY(JARRAY(I,6))=KPERP(JARRAY(I,6))*COS(FI(I)) KZ(JARRAY(I,6))=KPERP(JARRAY(I,6))*SIN(FI(I)) ENER(JARRAY(I,6))=EIN(I) 5731 CONTINUE C--------------------------------------------------------------------- C [ 7 ] Non-Equivalent Intervalley Emission II C--------------------------------------------------------------------- C Intervalley Emission (non-eq.) GX,LX,XL DO 5740 I=1,JCOUNT(7) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5740 CONTINUE DO 5741 I=1,JCOUNT(7) EIN(I)=ENEW(JARRAY(I,7))+EJUMP(IVAL(JARRAY(I,7)),7) # -EPHON(IVAL(JARRAY(I,7)),7) EIN(I)=AMAX1(EIN(I),1.E-24) IVAL(JARRAY(I,7))=IMAP(IVAL(JARRAY(I,7)),7) ELMASS(JARRAY(I,7))=EMAP(IVAL(JARRAY(I,7))) VAR1(JARRAY(I,7))=HBAR/ELMASS(JARRAY(I,7)) VAR2(JARRAY(I,7))=Q/ELMASS(JARRAY(I,7)) VAR3(JARRAY(I,7))=-.5*VAR2(JARRAY(I,7))*DTFL**2 VAR4(JARRAY(I,7))=VAR1(JARRAY(I,7))*DTFL VAR5(JARRAY(I,7))=VAR1(JARRAY(I,7))*HBAR*.5 VAR6(JARRAY(I,7))=AMAP2(IVAL(JARRAY(I,7))) VAR7(JARRAY(I,7))=4.*VAR6(JARRAY(I,7)) VAR8(JARRAY(I,7))=.5/VAR6(JARRAY(I,7)) AEAS(I)=VAR6(JARRAY(I,7))*EIN(I) GAM(JARRAY(I,7))=EIN(I)*(1.+AEAS(I)) R(I)=SQRT(GAM(JARRAY(I,7))*CC(IVAL(JARRAY(I,7)),1)) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,7))=R(I)*COBETA(I) KPERP(JARRAY(I,7))=SQRT(R(I)**2-KX(JARRAY(I,7))**2) KY(JARRAY(I,7))=KPERP(JARRAY(I,7))*COS(FI(I)) KZ(JARRAY(I,7))=KPERP(JARRAY(I,7))*SIN(FI(I)) ENER(JARRAY(I,7))=EIN(I) 5741 CONTINUE C--------------------------------------------------------------------- C [ 8 ] Equivalent Intervalley Absorption C--------------------------------------------------------------------- DO 5750 I=1,JCOUNT(8) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5750 CONTINUE DO 5751 I=1,JCOUNT(8) EIN(I)=ENEW(JARRAY(I,8))+EPHON(IVAL(JARRAY(I,8)),8) EIN(I)=AMAX1(EIN(I),1.E-24) K2NEW(JARRAY(I,8))=EIN(I)*(1.+VAR6(JARRAY(I,8))*EIN(I)) # *CC(IMAP(IVAL(JARRAY(I,8)),8),1) R(I)=SQRT(K2NEW(JARRAY(I,8))) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,8))=R(I)*COBETA(I) KPERP(JARRAY(I,8))=SQRT(K2NEW(JARRAY(I,8))- # KX(JARRAY(I,8))**2) KY(JARRAY(I,8))=KPERP(JARRAY(I,8))*COS(FI(I)) KZ(JARRAY(I,8))=KPERP(JARRAY(I,8))*SIN(FI(I)) ENER(JARRAY(I,8))=EIN(I) 5751 CONTINUE C--------------------------------------------------------------------- C [ 9 ] Equivalent Intervalley Emission C--------------------------------------------------------------------- DO 5760 I=1,JCOUNT(9) COBETA(I)=1.-2.*RAN(SEED) FI(I)=TWOPI*RAN(SEED) 5760 CONTINUE DO 5761 I=1,JCOUNT(9) EIN(I)=ENEW(JARRAY(I,9))-EPHON(IVAL(JARRAY(I,9)),9) EIN(I)=AMAX1(EIN(I),1.E-24) K2NEW(JARRAY(I,9))=EIN(I)*(1.+VAR6(JARRAY(I,9))*EIN(I)) # *CC(IMAP(IVAL(JARRAY(I,9)),9),1) R(I)=SQRT(K2NEW(JARRAY(I,9))) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KX(JARRAY(I,9))=R(I)*COBETA(I) KPERP(JARRAY(I,9))=SQRT(K2NEW(JARRAY(I,9))- # KX(JARRAY(I,9))**2) KY(JARRAY(I,9))=KPERP(JARRAY(I,9))*COS(FI(I)) KZ(JARRAY(I,9))=KPERP(JARRAY(I,9))*SIN(FI(I)) ENER(JARRAY(I,9))=EIN(I) 5761 CONTINUE C--------------------------------------------------------------------- C [ 10 ] Impurity Scattering C--------------------------------------------------------------------- DO 5770 I=1,JCOUNT(10) FI(I)=TWOPI*RAN(SEED) R1(I)=RAN(SEED) 5770 CONTINUE DO 5771 I=1,JCOUNT(10) R(I)=SQRT(K2NEW(JARRAY(I,10))) COBETA(I)=1.-2.*(1.-R1(I))/(1.+R1(I)*K2NEW(JARRAY(I,10))* # SCMAP) COBETA(I)=AMIN1(COBETA(I),1.) SIBETA(I)=SQRT(1.-COBETA(I)**2) COSETA(I)=KY(JARRAY(I,10))/KPERP(JARRAY(I,10)) SIETA(I)=KZ(JARRAY(I,10))/KPERP(JARRAY(I,10)) SICO(I)=SIBETA(I)*COS(FI(I)) SISI(I)=SIBETA(I)*SQRT(1-COS(FI(I))**2) C --------------------------------------------- C Calculate new wavevector and update energy C --------------------------------------------- KZ(JARRAY(I,10))=KZ(JARRAY(I,10))*COBETA(I)-KXNEW(JARRAY(I,10))* # SICO(I)*SIETA(I)+R(I)*SISI(I)*COSETA(I) KY(JARRAY(I,10))=KY(JARRAY(I,10))*COBETA(I)-KX(JARRAY(I,10))* # SICO(I)*COSETA(I)-R(I)*SISI(I)*SIETA(I) KX(JARRAY(I,10))=KX(JARRAY(I,10))*COBETA(I)+ # KPERP(JARRAY(I,10))*SICO(I) ENER(JARRAY(I,10))=ENEW(JARRAY(I,10)) 5771 CONTINUE C--------------------------------------------------------------------- C [ 11 ] Self-scattering C--------------------------------------------------------------------- DO 5780 I=1,JCOUNT(11) KX(JARRAY(I,11))=KXNEW(JARRAY(I,11)) ENER(JARRAY(I,11))=ENEW(JARRAY(I,11)) 5780 CONTINUE C===================================================================== C END OF FINAL STATE C===================================================================== C C--------------------------------------------------------------------- C Update K2PERP for flight C--------------------------------------------------------------------- DO 5810 I=1,NSCAT IE(JSCAT(I))=IFIX(ENER(JSCAT(I))/STEPJ)+1 IF (IE(JSCAT(I)).GT.NESTEP) IE(JSCAT(I))=NESTEP K2PERP(JSCAT(I))=KY(JSCAT(I))**2+KZ(JSCAT(I))**2 KPERP(JSCAT(I))=SQRT(K2PERP(JSCAT(I))) C --------------------------------------------- C Free flights for electrons C --------------------------------------------- KX(JSCAT(I))=KX(JSCAT(I))-CON1*FX*T2(I) K2(JSCAT(I))=KX(JSCAT(I))**2+K2PERP(JSCAT(I)) GAM(JSCAT(I))=VAR1(JSCAT(I))*K2(JSCAT(I))*HBAR*.5 ENEW(JSCAT(I))=(SQRT(1.+VAR7(JSCAT(I))*GAM(JSCAT(I)))-1.)* # VAR8(JSCAT(I)) IENEW(JSCAT(I))=IFIX(ENEW(JSCAT(I))/STEPJ)+1 IF (IENEW(JSCAT(I)).GE.NESTEP) IENEW(JSCAT(I))=NESTEP C--------------------------------------------------------------------- C Pick Rategam for new flight C--------------------------------------------------------------------- RATEGAM(JSCAT(I))=AMAX1(RLAMBDA(IE(JSCAT(I)),IVAL(JSCAT(I))), # RLAMBDA(IENEW(JSCAT(I)),IVAL(JSCAT(I)))) 5810 CONTINUE C--------------------------------------------------------------------- C Pick random # for new flight C--------------------------------------------------------------------- DO 5910 I=1,NSCAT RTEST(JSCAT(I))=-ALOG(RAN(SEED)) 5910 CONTINUE C ---------------------------------------------------------------- C If the particle is supposed to scatter a second time before the C end of the timestep, we force it to scatter just at the very C beginning of the next timestep. With small enough timestep, C double scattering should almost never occur. C ---------------------------------------------------------------- DO 5920 I=1,NSCAT RNEW(JSCAT(I))=RTEST(JSCAT(I))-RATEGAM(JSCAT(I))*T2(I) RTEST(JSCAT(I))=AMAX1(RNEW(JSCAT(I)),1.E-20) C ------------------------------------------------- C Update energy, index for particles that scattered C ------------------------------------------------- ENER(JSCAT(I))=ENEW(JSCAT(I)) IE(JSCAT(I))=IENEW(JSCAT(I)) C --------------------------------------------- C Update GAM for next loop C --------------------------------------------- GAM(JSCAT(I))=ENER(JSCAT(I))* # (1.+VAR6(JSCAT(I))*ENER(JSCAT(I))) 5920 CONTINUE C --------------------------------------------- C Calculate drift velocity of each particle C --------------------------------------------- DO 6010 I=1,NELEC ETEMP(I)=(2./3.)*ENER(I)/BOLTZ VX(I)=HBAR*KX(I)/ELMASS(I)/(1.+2.*VAR6(I)*ENER(I)) 6010 CONTINUE C===================================================================== C AVERAGES C===================================================================== C Update the running sum C AVELIT=0. AVENIT=0. DO 6015 ISS=1,NELEC AVELIT=AVELIT+VX(ISS) AVENIT=AVENIT+ENER(ISS) 6015 CONTINUE AVELO=AVELO+AVELIT AVENE=AVENE+AVENIT C--------------------------------------------------------------------- 1000 CONTINUE write(6,*)'end 1000' C--------------------------------------------------------------------- C Memorize status of ensemble at the end of the simulation C--------------------------------------------------------------------- IF (JYES.EQ.1) THEN AVELO=((AVELO/FLOAT(NELEC)*100)+(AVELOOLD*FLOAT(ITERTOT)))/ # (ITERMAX+ITERTOT) AVENE=((AVENE/FLOAT(NELEC))+(AVENEOLD*FLOAT(ITERTOT)))/ # FLOAT(ITERMAX+ITERTOT) ELSE AVELO=AVELO/FLOAT(NELEC)/FLOAT(ITERMAX)*100. AVENE=AVENE/FLOAT(NELEC)/FLOAT(ITERMAX) ENDIF C ---------------------------------------------------------------- ATEMP=AVENE*2./3./BOLTZ ITERTOT=ITERTOT+ITERMAX C ---------------------------------------------------------------- C Write out averages C ---------------------------------------------------------------- WRITE(12,*) 'Drift velocity = ',AVELO,' cm/s' WRITE(12,*) 'Average energy = ',AVENE,' J ',AVENE/Q,' eV' WRITE(12,*) 'Electron gas equivalent temperature = ',ATEMP,' K' WRITE(12,*) 'Total number of iterations = ',ITERTOT C ---------------------------------------------------------------- C Write out final distribution C ---------------------------------------------------------------- REWIND(11) DO 6030 I=1,NELEC WRITE(11,*)IVAL(I),KX(I),KY(I),KZ(I) 6030 CONTINUE WRITE(11,*)AVELO,AVENE,ITERTOT C ---------------------------------------------------------------- C Istograms of particle distribution C ---------------------------------------------------------------- DO 6040 I=1,500 IECOUNT(I)=0 6040 CONTINUE C --------------------------------------------------- DO 6050 I=1,NELEC IECOUNT(IE(I))=IECOUNT(IE(I))+1 6050 CONTINUE C --------------------------------------------------- REWIND 11 REWIND 13 CLOSE(11) CLOSE(13) STOP END C===================================================================== C END OF MAIN PROGRAM C===================================================================== C===================================================================== C C Create the scattering rate table C C===================================================================== SUBROUTINE SCAT(J) C ---------------------------------------------------------------- COMMON/TABLE1/SUMRATE(1000,3,15) COMMON/TABLE2/RATEMAX(1000,3),RLAMBDA(1000,3) COMMON/TABLE3/PIMPG(1000,6),PIMPL(1000,6),PIMPX(1000,6) C ---------------------------------------------------------------- C Physical Parameters C ---------------------------------------------------------------- COMMON/PHYS1/BOLTZ,EM0,EP0,HBAR,PI,Q,TEMP,TK,TWOPI C ---------------------------------------------------------------- C Material Related Parameters C ---------------------------------------------------------------- COMMON/MAT1/ALFA1(2),ALFA2(2),ALFA3(2) COMMON/MAT2/DEFGL(2),DEFGX(2),DEFLX(2),DEFLL(2),DEFXX(2) COMMON/MAT3/DEFAC(2),RHO(2),VS(2) COMMON/MAT4/GAP(2),GAPGL(2),GAPGX(2) COMMON/MAT5/EMSTAR1(2),EMSTAR2(2),EMSTAR3(2) COMMON/MAT6/EM1(2),EM2(2),EM3(2) COMMON/MAT7/EPZER(2),EPINF(2) COMMON/MAT8/TPOP(2),TKPOP(2) COMMON/MAT9/TINTGL(2),TINTGX(2),TINTLX(2),TEQINL(2),TEQINX(2) COMMON/MAT10/TKINT1(2),TKINT2(2),TKINT3(2),TKEIN2(2),TKEIN3(2) COMMON/MAT11/ZL,ZX COMMON/MAT12/DOP(5) COMMON/MAT13/CC(3,2) COMMON/MAT14/TWOA1(2),TWOA2(2),TWOA3(2) COMMON/MAT15/FOURA1(2),FOURA2(2),FOURA3(2) COMMON/MAT16/TAINV1(2),TAINV2(2),TAINV3(2) C ---------------------------------------------------------------- C Algorithm Parameters C ---------------------------------------------------------------- COMMON/ALG1/NELEC,NESTEP,EMAX,ESTEP C ---------------------------------------------------------------- DIMENSION PTAC(1000),PTAINT(1000),PTEINT(1000), # PTAINT3(1000),PTEINT3(1000),PTAOP(1000),PTEOP(1000) C ---------------------------------------------------------------- DIMENSION PTACS(1000),PTAINS(1000),PTEINS(1000),PTAINS3(1000), # PTEINS3(1000),PTAOPS(1000),PTEOPS(1000), # PTEIE(1000),PTAIE(1000) C ---------------------------------------------------------------- DIMENSION PTAC3(1000),PTAIN3(1000),PTEIN3(1000),PINAL3(1000), # PINEL3(1000),PTAOP3(1000),PTEOP3(1000),PTEIE3(1000), # PTAIE3(1000) C ---------------------------------------------------------------- DIMENSION SCRM(5),SCRM4(5) REAL CC C ---------------------------------------------------------------- C STEP = Step in energy for table [joules] C OMPOP = Polar optical phonon frequency C ---------------------------------------------------------------- c write(6,*) 'inside scat' STEP=ESTEP*BOLTZ OMPOP=TKPOP(J)/HBAR C--------------------------------------------------------------------- C C Constants for scattering rate calculations C C--------------------------------------------------------------------- SCRM(5) =(Q/BOLTZ)*(Q/EPZER(1))*DOP(5)/TEMP SCRM4(5)=4.0/SCRM(5) C ---------------------------------------------------------------- DIFCAP=1./EPINF(J)-1./EPZER(J) C ---------------------------------------------------------------- RAPT =TPOP(J)/TEMP RAPINT =TINTGL(J)/TEMP RAPINT3 =TINTGX(J)/TEMP RAPINT33=TINTLX(J)/TEMP REQINL =TEQINL(J)/TEMP REQINX =TEQINX(J)/TEMP C ---------------------------------------------------------------- BOSEOP =1./(EXP(RAPT)-1.) BOSEINT =1./(EXP(RAPINT)-1.) BOSEINT3=1./(EXP(RAPINT3)-1.) BOSEIN33=1./(EXP(RAPINT33)-1.) BOSEQL =1./(EXP(REQINL)-1.) BOSEQX =1./(EXP(REQINX)-1.) C ---------------------------------------------------------------- DEFQ =DEFGL(J)*DEFGL(J) DEFQ3 =DEFGX(J)*DEFGX(J) DEFQ33=DEFLX(J)*DEFLX(J) C ---------------------------------------------------------------- CONST1=Q/HBAR*Q CONST2=OMPOP*SQRT(EM1(J)/2.) CONST3=CONST1*CONST2*DIFCAP/4./PI C ---------------------------------------------------------------- DONST2 =OMPOP*SQRT(EM2(J)/2.) DONST23=OMPOP*SQRT(EM3(J)/2.) DONST3 =CONST1*DONST2*DIFCAP/4./PI DONST33=CONST1*DONST23*DIFCAP/4./PI C ---------------------------------------------------------------- C Optical phonon absorption C ---------------------------------------------------------------- CONST4 =CONST3*BOSEOP DONST4 =DONST3*BOSEOP DONST43=DONST33*BOSEOP C ---------------------------------------------------------------- C Optical phonon emission C ---------------------------------------------------------------- CONST5 =CONST3*(BOSEOP+1.) DONST5 =DONST3*(BOSEOP+1.) DONST53=DONST33*(BOSEOP+1.) C ---------------------------------------------------------------- C Acoustical phonon C ---------------------------------------------------------------- CONST6 =(DEFAC(J)/HBAR)*SQRT(TK)/(VS(J)*HBAR) CONST9 =PI*RHO(J) CONST10 =CONST6*EM1(J)*SQRT(2.*EM1(J))/CONST9*CONST6 DONST10 =CONST6*EM2(J)*SQRT(2.*EM2(J))/CONST9*CONST6 DONST103=CONST6*EM3(J)*SQRT(2.*EM3(J))/CONST9*CONST6 C ---------------------------------------------------------------- C Non-equivalent intervalley C ---------------------------------------------------------------- BRODO =CONST9*HBAR CONST12 =ZL*EM2(J)/HBAR*SQRT(EM2(J)/2.)/TKINT1(J)*DEFQ/BRODO CONST123=ZX*EM3(J)/HBAR*SQRT(EM3(J)/2.)/TKINT2(J)*DEFQ3/BRODO C ---------------------------------------------------------------- DONST12 =EM1(J)/HBAR*SQRT(EM1(J)/2.)/TKINT1(J)*DEFQ/BRODO DONST123=ZX*EM3(J)/HBAR*SQRT(EM3(J)/2.)/TKINT3(J)*DEFQ33/BRODO C ---------------------------------------------------------------- DONST312=EM1(J)/HBAR*SQRT(EM1(J)/2.)/TKINT2(J)*DEFQ3/BRODO DONST323=ZL*EM2(J)/HBAR*SQRT(EM2(J)/2.)/TKINT3(J)*DEFQ33/BRODO C ---------------------------------------------------------------- C Non Eq. intervalley absorption C ---------------------------------------------------------------- CONST13 =CONST12*BOSEINT CONST133=CONST123*BOSEINT3 DONST13 =DONST12*BOSEINT DONST133=DONST123*BOSEIN33 DONST313=DONST312*BOSEINT3 DONST332=DONST323*BOSEIN33 C ---------------------------------------------------------------- C Non Eq. intervalley emission C ---------------------------------------------------------------- CONST14 =CONST12*(BOSEINT+1.) CONST143=CONST123*(BOSEINT3+1.) DONST14 =DONST12*(BOSEINT+1.) DONST143=DONST123*(BOSEIN33+1.) DONST314=DONST312*(BOSEINT3+1.) DONST333=DONST323*(BOSEIN33+1.) C ---------------------------------------------------------------- C Equivalent intervalley C ---------------------------------------------------------------- DEFEQL =DEFLL(J)*DEFLL(J) DEFEQX =DEFXX(J)*DEFXX(J) CONEQ12 =ZL*EM2(J)/HBAR*SQRT(EM2(J)/2.)/TKEIN2(J)*DEFEQL/BRODO CONEQ123=ZX*EM3(J)/HBAR*SQRT(EM3(J)/2.)/TKEIN3(J)*DEFEQX/BRODO CONEQ14 =CONEQ12*(BOSEQL+1.) CONEQ143=CONEQ123*(BOSEQX+1.) CONEQ13 =CONEQ12*BOSEQL CONEQ133=CONEQ123*BOSEQX CONST16 =CONEQ14*(ZL-1.)/ZL CONST316=CONEQ143*(ZX-1.)/ZX CONST15 =CONEQ13*(ZL-1.)/ZL CONST315=CONEQ133*(ZX-1.)/ZX C ---------------------------------------------------------------- C Impurity scattering C ---------------------------------------------------------------- CONSIM =(Q/HBAR*Q/EPZER(J))**2*EM1(J)/HBAR/4./PI CONSIM2=CONSIM*(EM2(J)/EM1(J)) CONSIM3=CONSIM*(EM3(J)/EM1(J)) C===================================================================== C C CALCULATION OF SCATTERING RATES C C===================================================================== C C Gamma Valley C C--------------------------------------------------------------------- NPO1 =INT(TKPOP(J)/STEP)+1 NINA =INT((GAPGL(J)-TKINT1(J))/STEP)+1 NINE =INT((GAPGL(J)+TKINT1(J))/STEP)+1 NINNA=INT((GAPGX(J)-TKINT2(J))/STEP)+1 NINNE=INT((GAPGX(J)+TKINT2(J))/STEP)+1 C ---------------------------------------------------------------- C DO 10 - MAIN LOOP FOR ENERGY C ---------------------------------------------------------------- DO 10 I=1,NESTEP E=FLOAT(I)*ESTEP*BOLTZ C ---------------------------------------------------------------- C Acoustical phonon C ---------------------------------------------------------------- AE =ALFA1(J)*E AEQ =AE*AE UPAE =1.+AE GAM =E*UPAE RAGAM =SQRT(GAM) UPDAE =UPAE+AE FA =1.+1.333333*AEQ/UPDAE PTAC(I)=CONST10*RAGAM*FA C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to L) C ---------------------------------------------------------------- AMA=E+TKINT1(J) IF(I.LE.NINA) THEN PTAINT(I)=0. ELSE EAINT =AMA-GAPGL(J) REAINT =SQRT(EAINT*(1.+ALFA2(J)*EAINT)) TIP =UPAE/UPDAE TAP2 =(1.+ALFA2(J)*EAINT) PTAINT(I)=CONST13*REAINT*TIP*TAP2 ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley emission (to L) C ---------------------------------------------------------------- BIBE=E-TKINT1(J) IF(I.LE.NINE) THEN PTEINT(I)=0. ELSE EEINT =BIBE-GAPGL(J) REEINT =SQRT(EEINT*(1.+ALFA2(J)*EEINT)) TAP2 =(1.+ALFA2(J)*EEINT) PTEINT(I)=CONST14*REEINT*TIP*TAP2 ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to X) C ---------------------------------------------------------------- AMA3=E+TKINT2(J) IF(I.LE.NINNA) THEN PTAINT3(I)=0. ELSE EAINT3 =AMA3-GAPGX(J) REAINT3 =SQRT(EAINT3*(1.+ALFA3(J)*EAINT3)) TAP3 =(1.+ALFA3(J)*EAINT3) PTAINT3(I)=CONST133*REAINT3*TIP*TAP3 ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley emission (to X) C ---------------------------------------------------------------- BIBE3=E-TKINT2(J) IF(I.LE.NINNE) THEN PTEINT3(I)=0. ELSE EEINT3 =BIBE3-GAPGX(J) REEINT3 =SQRT(EEINT3*(1.+ALFA3(J)*EEINT3)) TAP3 =(1.+ALFA3(J)*EEINT3) PTEINT3(I)=CONST143*REEINT3*TIP*TAP3 ENDIF C ---------------------------------------------------------------- C Polar optical absorption C ---------------------------------------------------------------- DUPDAE =2.*UPDAE EAS =E+TKPOP(J) AEAS =ALFA1(J)*EAS UPAEAS =1.+AEAS GAMAS =EAS*UPAEAS RAGAMAS =SQRT(GAMAS) SOMGAMA =GAM+GAMAS PRAGAMA =RAGAM*RAGAMAS SRAGAMA =RAGAM+RAGAMAS DRAGAMA =RAGAMAS-RAGAM FIOREA =SRAGAMA/DRAGAMA ALOGA =LOG(FIOREA) PLUTOA =2.*UPAE*UPAEAS PIPPOA =PLUTOA+ALFA1(J)*SOMGAMA BAMBIA =PIPPOA+PLUTOA QUIA =PIPPOA*PIPPOA QUOA =-TWOA1(J)*PRAGAMA*BAMBIA FOA =(QUIA*ALOGA+QUOA)/(PLUTOA*DUPDAE) PTAOP(I)=CONST4*FOA/RAGAM C ---------------------------------------------------------------- C Polar optical emission C ---------------------------------------------------------------- IF(I.LE.NPO1) THEN PTEOP(I)=0. ELSE EEM =E-TKPOP(J) AEEM =ALFA1(J)*EEM UPAEEM =1.+AEEM GAMEM =EEM*UPAEEM RAGAMEM =SQRT(GAMEM) SOMGAME =GAM+GAMEM PRAGAME =RAGAM*RAGAMEM SRAGAME =RAGAM+RAGAMEM DRAGAME =RAGAM-RAGAMEM FIOREE =SRAGAME/DRAGAME ALOGE =LOG(FIOREE) PLUTOE =2.*UPAE*UPAEEM PIPPOE =PLUTOE+ALFA1(J)*SOMGAME BAMBIE =PIPPOE+PLUTOE QUIE =PIPPOE*PIPPOE QUOE =-TWOA1(J)*PRAGAME*BAMBIE FOE =(QUIE*ALOGE+QUOE)/(PLUTOE*DUPDAE) PTEOP(I)=CONST5*FOE/RAGAM ENDIF 10 CONTINUE C ---------------------------------------------------------------- C C L Valley C C ---------------------------------------------------------------- NPO2 =INT(TKPOP(J)/STEP)+1 NINT2=INT(TKEIN2(J)/STEP)+1 C ---------------------------------------------------------------- C DO 20 - MAIN LOOP FOR ENERGY C ---------------------------------------------------------------- DO 20 I=1,NESTEP E=FLOAT(I)*ESTEP*BOLTZ C ---------------------------------------------------------------- C Acoustical phonon C ---------------------------------------------------------------- AE2 =ALFA2(J)*E UPAE2 =1.+AE2 UPDAE2 =UPAE2+AE2 AEQ2 =AE2*AE2 GAM2 =E*UPAE2 DRAGAM =SQRT(GAM2) FA2 =1.+1.333333*AEQ2/UPDAE2 PTACS(I)=DONST10*DRAGAM*FA2 C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to Gamma) C ---------------------------------------------------------------- EFS =E+GAPGL(J) EAINS =EFS+TKINT1(J) REAINS =SQRT(EAINS) AEAINS =ALFA1(J)*EAINS UPAEAS =1.+AEAINS PAPPA =REAINS*UPAEAS*SQRT(UPAEAS) TOPPETE =UPAE2/UPDAE2 PTAINS(I)=DONST13*PAPPA*TOPPETE C ---------------------------------------------------------------- C Non Eq. intervalley emission (to Gamma) C ---------------------------------------------------------------- EEINS=EFS-TKINT1(J) IF(EEINS.LE.0.) THEN PTEINS(I)=0. ELSE REEINS =SQRT(EEINS) AEEINS =ALFA1(J)*EEINS UPAEEM =1.+AEEINS FRUTTA =REEINS*UPAEEM*SQRT(UPAEEM) PTEINS(I)=DONST14*FRUTTA*TOPPETE ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to X) C ---------------------------------------------------------------- EFS3 =E-(GAPGX(J)-GAPGL(J)) EAINS3=EFS3+TKINT3(J) IF(EAINS3.LE.0.) THEN PTAINS3(I)=0. ELSE REAINS3 =SQRT(EAINS3) AEAINS3 =ALFA3(J)*EAINS3 UPAEAS3 =1.+AEAINS3 PAPPA3 =REAINS3*UPAEAS3*SQRT(UPAEAS3) PTAINS3(I)=DONST133*PAPPA3*TOPPETE ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley emission (to X) C ---------------------------------------------------------------- EEINS3=EFS3-TKINT3(J) IF(EEINS3.LE.0.) THEN PTEINS3(I)=0. ELSE REEINS3 =SQRT(EEINS3) AEEINS3 =ALFA3(J)*EEINS3 UPAEEM3 =1.+AEEINS3 FRU3 =REEINS3*UPAEEM3 PTEINS3(I)=DONST143*FRU3*TOPPETE ENDIF C ---------------------------------------------------------------- C Polar optical absorption C ---------------------------------------------------------------- DUPDAE2 =2.*UPDAE2 EAS2 =E+TKPOP(J) AEAS2 =ALFA2(J)*EAS2 UPAEAS2 =1.+AEAS2 GAMAS2 =EAS2*UPAEAS2 RAGAMAS2 =SQRT(GAMAS2) SOMGAMA2 =GAM2+GAMAS2 PRAGAMA2 =DRAGAM*RAGAMAS2 SRAGAMA2 =DRAGAM+RAGAMAS2 DRAGAMA2 =RAGAMAS2-DRAGAM FIOR2 =SRAGAMA2/DRAGAMA2 DALOGA =LOG(FIOR2) PLUTO2 =2.*UPAE2*UPAEAS2 PIPPO2 =PLUTO2+ALFA2(J)*SOMGAMA2 BAMBI2 =PIPPO2+PLUTO2 QUI2 =PIPPO2**2 QUO2 =-TWOA2(J)*PRAGAMA2*BAMBI2 FOA2 =(QUI2*DALOGA+QUO2)/(PLUTO2*DUPDAE2) PTAOPS(I)=DONST4*FOA2/DRAGAM C ---------------------------------------------------------------- C Polar optical emission C ---------------------------------------------------------------- IF(I.LE.NPO2) THEN PTEOPS(I)=0. ELSE EEM2 =E-TKPOP(J) AEEM2 =ALFA2(J)*EEM2 UPAEM2 =1.+AEEM2 GAMEM2 =EEM2*UPAEM2 RAGAMEM2 =SQRT(GAMEM2) PRAGAME2 =DRAGAM*RAGAMEM2 SOMGAME2 =GAM2+GAMEM2 SRAGAME2 =DRAGAM+RAGAMEM2 DRAGAME2 =DRAGAM-RAGAMEM2 FIOR2E =SRAGAME2/DRAGAME2 ALOGE2 =LOG(FIOR2E) PLUTOE2 =2.*UPAE2*UPAEM2 PIPPOE2 =PLUTOE2+ALFA2(J)*SOMGAME2 BAMBIE2 =PIPPOE2+PLUTOE2 QUIE2 =PIPPOE2*PIPPOE2 QUOE2 =-TWOA2(J)*PRAGAME2*BAMBIE2 FOE2 =(QUIE2*ALOGE2+QUOE2)/(PLUTOE2*DUPDAE2) PTEOPS(I)=DONST5*FOE2/DRAGAM ENDIF C ---------------------------------------------------------------- C Equivalent intervalley emission C ---------------------------------------------------------------- IF(I.LE.NINT2) THEN PTEIE(I)=0. ELSE EIN =E-TKEIN2(J) UPINT2 =1.+ALFA2(J)*EIN PTEIE(I)=CONST16*SQRT(EIN*UPINT2)*UPINT2*UPAE2/UPDAE2 ENDIF C ---------------------------------------------------------------- C Equivalent intervalley absorption C ---------------------------------------------------------------- EIN =E+TKEIN2(J) UPINT2 =1+ALFA2(J)*EIN PTAIE(I)=CONST15*SQRT(EIN*UPINT2)*UPINT2*UPAE2/UPDAE2 20 CONTINUE C ---------------------------------------------------------------- C C X Valley C C ---------------------------------------------------------------- NPO3 =INT(TKPOP(J)/STEP)+1 NINT3=INT(TKEIN3(J)/STEP)+1 C ---------------------------------------------------------------- C DO 30 - MAIN ENERGY LOOP C ---------------------------------------------------------------- DO 30 I=1,NESTEP E=FLOAT(I)*ESTEP*BOLTZ C ---------------------------------------------------------------- C Acoustical phonon C ---------------------------------------------------------------- AE3 =ALFA3(J)*E UPAE3 =1.+AE3 UPDAE3 =UPAE3+AE3 AEQ3 =AE3*AE3 GAM3 =E*UPAE3 DRAGAM3 =SQRT(GAM3) FA3 =1.+1.333333*AEQ3/UPDAE3 PTAC3(I)=DONST103*DRAGAM3*FA3 C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to Gamma) C ---------------------------------------------------------------- EEFS =E+GAPGX(J) EEAIN3 =EEFS+TKINT2(J) REAIN3 =SQRT(EEAIN3) AEAIN3 =ALFA1(J)*EEAIN3 UPAEA3 =1.+AEAIN3 PA3 =REAIN3*UPAEA3*SQRT(UPAEA3) TOPP3 =UPAE3/UPDAE3 PTAIN3(I)=DONST313*PA3*TOPP3 C ---------------------------------------------------------------- C Non Eq. intervalley emission (to Gamma) C ---------------------------------------------------------------- EEIN3=EEFS-TKINT2(J) IF(EEIN3.LE.0.) THEN PTEIN3(I)=0. ELSE REEIN3 =SQRT(EEIN3) AEEIN3 =ALFA1(J)*EEIN3 UPAEM3 =1.+AEEIN3 FRUFRU =REEIN3*UPAEM3*SQRT(UPAEM3) PTEIN3(I)=DONST314*FRUFRU*TOPP3 ENDIF C ---------------------------------------------------------------- C Non Eq. intervalley absorption (to L) C ---------------------------------------------------------------- EFF3 =EEFS-GAPGL(J) EAIN33 =EFF3+TKINT3(J) GGAM =1.+ALFA2(J)*EAIN33 GIGAM =EAIN33*GGAM PINAL3(I)=DONST332*SQRT(GIGAM)*GGAM*TOPP3 C ---------------------------------------------------------------- C Non Eq. intervalley emission (to L) C ---------------------------------------------------------------- EEIN33=EFF3-TKINT3(J) IF(EEIN33.LE.0.) THEN PINEL3(I)=0. ELSE GGAM =1.+ALFA2(J)*EEIN33 GIGAM =EAIN33*GGAM PINEL3(I)=DONST333*SQRT(GIGAM)*GGAM*TOPP3 ENDIF C ---------------------------------------------------------------- C Polar optical absorption C ---------------------------------------------------------------- EAS3 =E+TKPOP(J) DUPDAE3 =2.*UPDAE3 AEAS3 =ALFA3(J)*EAS3 UPAEAS3 =1.+AEAS3 GAMAS3 =EAS3*UPAEAS3 RAGAMAS3 =SQRT(GAMAS3) SOMGAMA3 =GAM3+GAMAS3 PRAGAMA3 =DRAGAM3*RAGAMAS3 SRAGAMA3 =DRAGAM3+RAGAMAS3 DRAGAMA3 =RAGAMAS3-DRAGAM3 FIOR3 =SRAGAMA3/DRAGAMA3 TALOGA =ALOG(FIOR3) PLUTO3 =2.*UPAE3*UPAEAS3 PIPPO3 =PLUTO3+ALFA3(J)*SOMGAMA3 BAMBI3 =PIPPO3+PLUTO3 QUI3 =PIPPO3*PIPPO3 QUO3 =-TWOA3(J)*PRAGAMA3*BAMBI3 FOA3 =(QUI3*TALOGA+QUO3)/(PLUTO3*DUPDAE3) PTAOP3(I)=DONST43*FOA3/DRAGAM3 C ---------------------------------------------------------------- C Polar optical emission C ---------------------------------------------------------------- IF(I.LE.NPO3) THEN PTEOP3(I)=0. ELSE EEM3 =E-TKPOP(J) AEEM3 =ALFA3(J)*EEM3 UPAEM3 =1.+AEEM3 GAMEM3 =EEM3*UPAEM3 RAGAMEM3 =SQRT(GAMEM3) PRAGAME3 =DRAGAM3*RAGAMEM3 SOMGAME3 =GAM3+GAMEM3 SRAGAME3 =DRAGAM3+RAGAMEM3 DRAGAME3 =DRAGAM3-RAGAMEM3 FIOR3E =SRAGAME3/DRAGAME3 ALOGE3 =ALOG(FIOR3E) PLUTOE3 =2.*UPAE3*UPAEM3 PIPPOE3 =PLUTOE3+ALFA3(J)*SOMGAME3 BAMBIE3 =PIPPOE3+PLUTOE3 QUOE3 =-TWOA3(J)*PRAGAME3*BAMBIE3 FOE3 =(QUI3*ALOGE3+QUOE3)/(PLUTOE3*DUPDAE3) PTEOP3(I)=DONST53*FOE3/DRAGAM3 ENDIF C ---------------------------------------------------------------- C Equivalent intervalley emission C ---------------------------------------------------------------- IF(I.LE.NINT3) THEN PTEIE3(I)=0. ELSE EIN3 =E-TKEIN3(J) UPINT3 =1.+ALFA3(J)*EIN3 GAG3 =EIN3*UPINT3 PTEIE3(I)=CONST316*SQRT(GAG3)*UPINT3*TOPP3 ENDIF C ---------------------------------------------------------------- C Equivalent intervalley absorption C ---------------------------------------------------------------- EIN3 =E+TKEIN3(J) UPINT3 =1.+ALFA3(J)*EIN3 GAGA3 =EIN3*UPINT3 PTAIE3(I)=CONST315*SQRT(GAGA3)*UPINT3*TOPP3 30 CONTINUE C ---------------------------------------------------------------- C Compute impurity scattering rates C ---------------------------------------------------------------- C Values are for GaAs bulk C ---------------------------------------------------------------- DO 2000 I=1,NESTEP E=FLOAT(I)*ESTEP*BOLTZ C ---------------------------------------------------------------- C Gamma Valley C ---------------------------------------------------------------- AE =ALFA1(J)*E AEQ =AE*AE UPAE =1.+AE UPDAE =UPAE+AE GAM =E*UPAE RQ =GAM*CC(1,J) AMU =1./(RQ*SCRM4(5)) R =SQRT(RQ) UMU =1.+AMU PAR =1.+2.*AE*UMU PIMPG(I,1)=CONSIM/(RQ*R*UPDAE)*DOP(5)*(AEQ+AE*PAR*LOG(AMU/UMU)+ # PAR*PAR/(4.*AMU*UMU)) PIMPG(I,1)=AMIN1(PIMPG(I,1),1.E14) C ---------------------------------------------------------------- C L Valley C ---------------------------------------------------------------- AE2 =ALFA2(J)*E UPAE2 =1.+AE2 UPDAE2 =UPAE2+AE2 GAM2 =E*UPAE2 RQ2 =GAM2*CC(2,J)*(EM2(J)/EM1(J)) R =SQRT(RQ2) AMU =1./(RQ2*SCRM4(5)) UMU =1.+AMU PAR =1.+2.*UMU*AE2 COFFA =CONSIM2/(RQ2*R*UPDAE2) PIMPL(I,1)=COFFA*DOP(5)*(AE2**2+AE2*PAR*LOG(AMU/UMU)+PAR**2 # /(4.*AMU*UMU)) PIMPL(I,1)=AMIN1(PIMPL(I,1),1.E14) C ---------------------------------------------------------------- C X Valley C ---------------------------------------------------------------- AE3 =ALFA3(J)*E UPAE3 =1.+AE3 UPDAE3 =UPAE3+AE3 GAM3 =E*UPAE3 RQ3 =GAM3*CC(3,J)*(EM3(J)/EM1(J)) AMU =1./(RQ3*SCRM4(5)) R =SQRT(RQ3) UMU =1.+AMU PAR =1.+2.*UMU*AE3 COFFA =CONSIM3/(RQ3*R*UPDAE3) PIMPX(I,1)=COFFA*DOP(5)*(AE3**2+AE3*PAR*LOG(AMU/UMU)+PAR**2/ # (4.*AMU*UMU)) PIMPX(I,1)=AMIN1(PIMPX(I,1),1.E14) 2000 CONTINUE C ---------------------------------------------------------------- C Compute partial sums of scattering events SUMRATE(I,J,K) C C I = Energy Index C J = Valley Index C K = Scattering Process C ---------------------------------------------------------------- C Polar optical abs. C ---------------------------------------------------------------- DO 800 I=1,NESTEP SUMRATE(I,1,1)=PTAOP(I) SUMRATE(I,2,1)=PTAOPS(I) SUMRATE(I,3,1)=PTAOP3(I) 800 CONTINUE C ---------------------------------------------------------------- C Polar opt. emission C ---------------------------------------------------------------- DO 801 I=1,NESTEP SUMRATE(I,1,2)=SUMRATE(I,1,1)+PTEOP(I) SUMRATE(I,2,2)=SUMRATE(I,2,1)+PTEOPS(I) SUMRATE(I,3,2)=SUMRATE(I,3,1)+PTEOP3(I) 801 CONTINUE C ---------------------------------------------------------------- C Acoustic C ---------------------------------------------------------------- DO 802 I=1,NESTEP SUMRATE(I,1,3)=SUMRATE(I,1,2)+PTAC(I) SUMRATE(I,2,3)=SUMRATE(I,2,2)+PTACS(I) SUMRATE(I,3,3)=SUMRATE(I,3,2)+PTAC3(I) 802 CONTINUE C ---------------------------------------------------------------- C Intervalley Absorption C In the order: G to L; L to G; X to G C ---------------------------------------------------------------- DO 803 I=1,NESTEP SUMRATE(I,1,4)=SUMRATE(I,1,3)+PTAINT(I) SUMRATE(I,2,4)=SUMRATE(I,2,3)+PTAINS(I) SUMRATE(I,3,4)=SUMRATE(I,3,3)+PTAIN3(I) 803 CONTINUE C ---------------------------------------------------------------- C Intervalley Emission C In the order: G to L; L to G; X to G C ---------------------------------------------------------------- DO 804 I=1,NESTEP SUMRATE(I,1,5)=SUMRATE(I,1,4)+PTEINT(I) SUMRATE(I,2,5)=SUMRATE(I,2,4)+PTEINS(I) SUMRATE(I,3,5)=SUMRATE(I,3,4)+PTEIN3(I) 804 CONTINUE C ---------------------------------------------------------------- C Intervalley Absorption C In the order: G to X; L to X; X to L C ---------------------------------------------------------------- DO 805 I=1,NESTEP SUMRATE(I,1,6)=SUMRATE(I,1,5)+PTAINT3(I) SUMRATE(I,2,6)=SUMRATE(I,2,5)+PTAINS3(I) SUMRATE(I,3,6)=SUMRATE(I,3,5)+PINAL3(I) 805 CONTINUE C ---------------------------------------------------------------- C Intervalley Emission C In the order: G to X; L to X; X to L C ---------------------------------------------------------------- DO 806 I=1,NESTEP SUMRATE(I,1,7)=SUMRATE(I,1,6)+PTEINT3(I) SUMRATE(I,2,7)=SUMRATE(I,2,6)+PTEINS3(I) SUMRATE(I,3,7)=SUMRATE(I,3,6)+PINEL3(I) 806 CONTINUE C ---------------------------------------------------------------- C Equivalent Intervalley Absorption C In the order: G (always 0); L; X C ---------------------------------------------------------------- DO 807 I=1,NESTEP SUMRATE(I,1,8)=SUMRATE(I,1,7)+0. SUMRATE(I,2,8)=SUMRATE(I,2,7)+PTAIE(I) SUMRATE(I,3,8)=SUMRATE(I,3,7)+PTAIE3(I) 807 CONTINUE C ---------------------------------------------------------------- C Equivalent Intervalley Emission C In the order: G (always 0); L; X C ---------------------------------------------------------------- DO 808 I=1,NESTEP SUMRATE(I,1,9)=SUMRATE(I,1,8)+0. SUMRATE(I,2,9)=SUMRATE(I,2,8)+PTEIE(I) SUMRATE(I,3,9)=SUMRATE(I,3,8)+PTEIE3(I) 808 CONTINUE C ---------------------------------------------------------------- DO 809 I=1,NESTEP RATEMAX(I,1)=SUMRATE(I,1,9) RATEMAX(I,2)=SUMRATE(I,2,9) RATEMAX(I,3)=SUMRATE(I,3,9) 809 CONTINUE C ---------------------------------------------------------------- C PIMP..(I,1) is the impurity scattering rate in bulk GaAs C ( region 1 selected ). The program is dimensioned to C include up to 5 regions with different doping for extension C to device simulation. C ---------------------------------------------------------------- DO 810 I=1,NESTEP RLAMBDA(I,1)=RATEMAX(I,1)+PIMPG(I,1) RLAMBDA(I,2)=RATEMAX(I,2)+PIMPL(I,1) RLAMBDA(I,3)=RATEMAX(I,3)+PIMPX(I,1) 810 CONTINUE C ---------------------------------------------------------------- RETURN END C--------------------------------------------------------------------- C--------------------------------------------------------------------- C This is the function that generates the random number. C--------------------------------------------------------------------- FUNCTION RAN(SEED) DOUBLE PRECISION SEED,A,B,DMOD A=SEED*7.**5 B=2.**31-1. SEED=DMOD(A,B) RAN=1./2.**31*SEED RETURN END