##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ##PNJunctionAbrupt.py ## ##PURPOSE: ## To serve as a starting point for calculating the field and potential ## in an abrupt PN junction, i.e one in which the donor and acceptor ## concentrations are described by step functions ## ##KEYWORDS: ## Mode: int ## 0 - use an explicit method ## 1 - use an implicit method ## TvFlag: int ## 0 - Don't plot ## 1 - Plot according to TvPPT ## TvPPT: int ## 0 - Plot every step ## 1 - Plot every 200th step ## NSteps: int ## The number of steps to execute ## SaveFig: int ## 0 - Don't save any figures ## 1 - Save the figures at the end ## ##EXAMPLES ## run PNJunctionAbrupt.py ## ################################################################### import matplotlib.pylab as mplot from numpy import * from pylab import * from LU_decomposition import * import pdb execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py') PlotDir = '/nfs/slac/g/ki/ki03/lances/PNSimulations/' ##KEYWORDS keywords = ['Mode=', 'TvFlag=', 'NSteps=', 'SaveFig=','SaveMultFigs=',\ 'TvPPT='] Mode = 0 TvFlag = 0 NSteps = 10 SaveMultFigs=0 SaveFig = 0 TvPPT = 0 opts, extraparams = getopt.getopt(sys.argv[1:],'',keywords) for o,p in opts: if o in ['--Mode']: Mode = int(p) elif o in ['--TvFlag']: TvFlag = int(p) elif o in ['--NSteps']: NSteps = int(p) elif o in ['--SaveFig']: SaveFig= int(p) elif o in ['--SaveMultFigs']: SaveMultFigs=int(p) elif o in ['--TvPPT']: TvPPT = int(p) ##MODE STRING if Mode == 0: ModeStr = 'Explicit' else: ModeStr = 'Implicit' ##CONSTANTS: e_0 = 8.85419e-12 # F * m^-1 - 8.85419 e-14 F*cm^-1 e_si = 11.8 # Relative Permittivity for Si q = 1.6022e-19 # Coulombs per electron k = 1.38066e-23 # Joules per Kelvin T = double(300) # Kelvin ##Effective Masses m_e = 9.1e-31 # Kg m_ee = 0.26 # effective mass of electrons/mass of electron m_eh = 0.386 # effective mass of holes/mass of electron ##Mobilities mu_n = 1000./1e4 # m^2/V s mu_p = 400./1e4 # m^2/V s ##Diffusion constants D_n = k*T*mu_n/q D_p = k*T*mu_p/q ##Parameters x_0 = 0 # In meters x_1 = 1.e-6 # In meters dx = (x_1-x_0)/5000. # In meters dx = double(dx) X = frange(x_0,x_1,dx) # In Microns N = len(X) NN = len(X)-1 MN = N/2 ##Intrinsic Concentration n_i = 1.e16 ##Donor Concentration N_d_max = 1.e22 N_d = zeros(N, dtype=double) N_d[0:N/2] = 0 N_d[N/2:N] = double(N_d_max) ##Acceptor Concentration N_a_max = 1.e22 N_a = zeros(N, dtype=double) N_a[0:N/2] = double(N_a_max) N_a[N/2:N] = 0 ##Difference DC = N_d-N_a ##Assume Phi at x=-inf is 0, calculate Phi_max according to ##Ashcroft Mermin p. 575 E_g = 1.12*q # Joules = eV*q N_c = 2.5*(m_ee)**1.5*(T/300)**1.5*10.**18 P_v = 2.5*(m_eh)**1.5*(T/300)**1.5*10.**18 ##Built-in-potential calculated in Ashcroft-Mermin Phi_Max = (E_g+k*T*log(N_d_max*N_a_max/(N_c*P_v)))/q ##Built-in potential calculated in Neudeck V_bi = k*T*log(N_d_max*N_a_max/n_i**2)/q ##Calculation of depletion region width x_p = ((2*e_si*e_0*V_bi/q)*(N_d_max/(N_a_max*(N_a_max+N_d_max))))**.5 x_n = ((2*e_si*e_0*V_bi/q)*(N_a_max/(N_d_max*(N_a_max+N_d_max))))**.5 E_Max = (q*N_d_max/(e_si*e_0))*x_n Phi_Max = (E_g+k*T*log(N_d_max*N_a_max/(n_i**2)))/q ##Plot the concentrations if TvFlag == 1: mplot.figure(0) mplot.hold(True) mplot.plot(X, N_a, 'r') mplot.plot(X, N_d, 'b') ##Now setup and integrate PhiPINF = double(Phi_Max) #Potential at infinity PhiMINF = double(0.) #Potential at minus infinity Phi = zeros(N, dtype=double) E = zeros(N, dtype=double) E_bar = zeros(N, dtype = double) E_del = zeros(N, dtype = double) ##Set up concentrations n = zeros(N, dtype=double) n[0:N/2] = n_i**2/N_a_max n[N/2:N] = N_d_max p = zeros(N, dtype=double) p[0:N/2] = N_a_max p[N/2:N] = n_i**2/N_d_max ##Boundary conditions according to Kurata p 37 n[NN] = DC[NN]/2*(1+(1+(2*n_i/DC[NN])**2)**.5) p[0] = -DC[0]/2*(1+(1+(2*n_i/DC[0])**2)**.5) n[0] = n_i**2/p[0] p[NN] = n_i**2/n[NN] Phi[0] = PhiMINF Phi[1] = PhiMINF+1.e-8 dt = 1.e-22 cfl = dt/dx ##Set up the matrices that will be used to calculate the FD's ##******************************************************SPATIAL DER. MATRICES ##The equation we look to solve is M*E = N1*rho ##M and N are difference operator matrices and f=N1*rho ##This is the matrix N1_m = -ones(N); N1_m[0]=0; N1_m[N-1]=0; N1_l = ones(N-1); N1_l[N-2]=0; N1 = diag(N1_m)+diag(N1_l,-1) ##These are vectors that can be used without wasting 0's N1_mv= -ones(N-2) N1_lv= ones(N-2) ##Now we can try to solve for M*E = f by LU decomposition ##Boundary conditions have E=0 at edges a1 = ones(N)*2; a1[0]=1; a1[N-1]=1; b1 = -ones(N-1) c1 = -ones(N-1); c1[0]=0; b1[N-2]=0 L, U, alpha1, gamma1 = LU_decomposition(a1,b1,c1) ############################################################################# for Sn in arange(NSteps): ##Calculate the field with an implicit method rho_i+1/2, rho_i-1/2, etc. rho = q*(p - n - N_a + N_d) ##Do the difference on the charge density f = zeros(N, dtype=double) f[1:NN-1] = dx*(rho[0:NN-2]-rho[1:NN-1])/(e_si*e_0) f[0] = 0; f[NN]=0 #f = dx*dot(N1, rho)/(e_si*e_0) ##Now we can try to solve for M*E = f by LU decomposition E_new = LU_solve(alpha1, b1, gamma1, f) E = E_new.copy() E_bar[0:NN] = (E[0:NN]+E[1:N])/2 E_del[0:NN] = (E[1:N] -E[0:NN])/dx del E_new ##Plot the E Field every time if TvFlag == 1 and TvPPT==0: mplot.figure(0) mplot.clf() mplot.plot(X,E) ##Calculate the electron density n_new = zeros(N,dtype=double) p_new = zeros(N,dtype=double) if Mode == 0: ##Do full range for i in arange(N-2)+1: n_new[i] = n[i] + (cfl)*(\ n[i]*mu_n*E_del[i] + \ mu_n*E_bar[i]*(n[i+1]-n[i-1])/(2*dx) + \ D_n*(n[i+1]-2*n[i]+n[i-1])/dx**2) p_new[i] = p[i] + (cfl)*(\ -p[i]*mu_p*E_del[i] - \ mu_p*E_bar[i]*(p[i+1]-p[i-1])/(2*dx) + \ D_p*(p[i+1]-2*p[i]+p[i-1])/dx**2) ##Handle boundaries n_new[0] = n[0] ; p_new[0]=p[0] n_new[NN] = n[NN] ; p_new[NN]=p[NN] elif Mode == 1: #Form the difference and average of E, keeping it the same length as rho a2 = 1+cfl*( mu_p*E_del[0:N]+2*D_p/dx**2); a2[0]=1; a2[NN]=1; b2 = cfl*( -mu_p*E_bar[1:N]/(2*dx) -D_p/dx**2); b2[NN-1]=0 c2 = cfl*( mu_p*E_bar[0:NN]/(2*dx) -D_p/dx**2); c2[0]=0 L2, U2, alpha2, gamma2 = LU_decomposition(a2,b2,c2) p_new = LU_solve(alpha2, b2, gamma2, p) del L2, U2, alpha2, gamma2, a2, b2, c2 a3 = 1-cfl*( mu_n*E_del[0:N ]-2*D_n/dx**2); a3[0]=1; a3[NN]=1; b3 = cfl*( mu_n*E_bar[1:N]/(2*dx) -D_n/dx**2); b3[NN-1]=0 c3 = cfl*(-mu_n*E_bar[0:NN]/(2*dx) -D_n/dx**2); c3[0]=0 L3, U3, alpha3, gamma3 = LU_decomposition(a3,b3,c3) n_new = LU_solve(alpha3, b3, gamma3, n) del L3, U3, alpha3, gamma3, a3, b3, c3 ##Compute an error sum Error = max(((p_new**2-p**2))**.5) Verbosity=1 if Verbosity == 1: print str(Sn)+', E min: ' +str(min(E))+' , E max: ' + \ str(max(E)) + ' , Error: '+str(Error) ##Plot Every Step if TvFlag == 1 and TvPPT == 0: mplot.figure(1) mplot.clf() mplot.plot(X,n_new) mplot.plot(X,p_new) ##Only Plot every 200th Step if TvFlag == 1 and TvPPT == 1 and Sn%200 == 0: if Sn == 0: mplot.close('all') mplot.hold(True) else: mplot.figure(0) mplot.plot(X,E) mplot.figure(1) mplot.semilogy(X,p) mplot.semilogy(X,n) mplot.figure(2) mplot.plot(X,rho) ##Clear Memory if Sn < NSteps-1: del n,p p = p_new.copy() n = n_new.copy() del n_new, p_new if Sn % 500 == 0 and Sn > 0 and SaveMultFigs == 1: print 'Saving Figure' PlotLast = 1 ##After loop is complete, plot if PlotLast == 1: mplot.figure(0) mplot.clf() mplot.plot(X,E) if SaveFig == 1: mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_EField_'+str(NSteps)+'Steps.png') mplot.figure(1) mplot.clf() mplot.semilogy(X,n_new) mplot.semilogy(X,p_new) if SaveFig == 1: mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_pv_nc_'+str(NSteps)+'Steps.png') rho=p - n - N_a + N_d mplot.figure(2) mplot.clf() mplot.plot(X,rho) mplot.ylim([-1.5*N_d_max,1.5*N_d_max]) if SaveFig == 1: mplot.savefig(PlotDir+'PNJunction_'+ModeStr+'_Rho_'+str(NSteps)+'Steps.png')