Resistor (1D) - Scharfetter Gummel Method

   #--------------------
   # Create 1D structure
   #--------------------
   #Grid
   line x loc=0.0 spac=0.1 tag=Top
   line x loc=0.5 spac=0.01
   line x loc=1.0 spac=0.1 tag=Bottom
   
   mater add name=Silicon
   region Silicon xlo=Top xhi=Bottom
   init
   
   #Contacts
   contact name=GND Silicon xlo=-0.1 xhi=0.1 add
   contact name=VSS Silicon xlo=0.9  xhi=1.1 add
   
   #-----------------------------
   # Declare solution variables
   #-----------------------------
   DevicePackage
   solution add name=DevPsi pde solve negative damp
   solution add name=Elec   pde solve !negative
   solution add name=Hole   pde solve !negative
   
   #-----------------------------
   # Define constants
   #-----------------------------
   set T 300.0
   set k 1.38066e-23
   set q 1.619e-19
   set Vt [expr {$k*$T/$q}]
   set ni 1.1e10
   set esi [expr 11.8 * 8.85418e-14]
   set eps [expr $esi / $q]
   set Emob 350.0
   set Hmob 150.0
   set small 1.0e-10
   
   #-----------------------------
   # Ionized dopant profile
   #-----------------------------
   sel z=1.0e20 name=Nd
   sel z=1.0e17 name=Na
   sel z=(Nd-Na) name=Doping
   
   #-----------------------------
   # Bulk Equations
   #-----------------------------
   set eqnP "$eps * grad(DevPsi) + Doping - Elec + Hole"
   set eqnE "ddt(Elec) - ($Emob) * $Vt * sgrad(Elec, DevPsi/$Vt)"
   set eqnH "ddt(Hole) - ($Hmob) * $Vt * sgrad(Hole, -DevPsi/$Vt)"
   pdbSetString Silicon DevPsi Equation $eqnP
   pdbSetString Silicon Elec   Equation $eqnE
   pdbSetString Silicon Hole   Equation $eqnH
   
   pdbSetDouble Silicon DevPsi DampValue $Vt
   pdbSetDouble Silicon DevPsi Abs.Error 1.0e-9
   pdbSetDouble Silicon Elec   Abs.Error 1.0e-5
   pdbSetDouble Silicon Hole   Abs.Error 1.0e-5
   
   #-----------------------------
   # Contact Equations
   #-----------------------------
   proc OhmicContact {Contact} {
       global Vt ni
       pdbSetBoolean $Contact Elec Flux 1
       pdbSetBoolean $Contact Hole Flux 1
       pdbSetBoolean $Contact DevPsi Flux 1
       pdbSetBoolean $Contact Elec Fixed 1
       pdbSetBoolean $Contact Hole Fixed 1
       pdbSetBoolean $Contact DevPsi Fixed 1
       pdbSetDouble $Contact Elec Flux.Scale 1.619e-19
       pdbSetDouble $Contact Hole Flux.Scale 1.619e-19
       pdbSetString $Contact DevPsi Equation "Nd - Na - Elec + Hole"
       pdbSetString $Contact Elec Equation "DevPsi - $Vt*log((Elec)/$ni) -$Contact"
       pdbSetString $Contact Hole Equation "DevPsi + $Vt*log((Hole)/$ni) -$Contact"
   }
   OhmicContact VSS
   OhmicContact GND
   
   #-----------------------------
   # Initial Conditions
   #-----------------------------
   #Bias contacts
   contact name=VSS voltage supply=0.0
   contact name=GND voltage supply=0.0
   
   #Initial Guess
   sel z= {
       (Doping>0.0)
       ?  (0.025*log( (Doping + 1e-10) / 1.1e10))
       :  (-0.025*log(-(Doping + 1e-10) / 1.1e10))} name = DevPsi
   sel z=$ni*exp(DevPsi/$Vt)  name=Elec
   sel z=$ni*exp(-DevPsi/$Vt) name=Hole
   
   #-----------------------------------
   # 1st DC Solve at Equilibrium (0V)
   #-----------------------------------
   device
   puts "Electron Flux [contact name=VSS sol=Elec flux]"
   puts "Hole Flux [contact name=VSS sol=Hole flux]"
   
   window row=1 col=3
   #Plot the equilibrium concentration profiles
   foreach var {Doping Elec Hole} {
       sel z=log10(abs($var+1.0))
       plot1d graph=Doping
   }
   
   #--------------------------------------
   # Ramp the DC Bias to make I-V plot
   #--------------------------------------
   #initialize an array, bias (Vb) and counter (i), and create and new graph window
   array set curr {}
   set i 1
   
   #Ramp from 0.0-1.2V and Plot
   for {set bias 0.0} {$bias < 1.2} {set bias [expr $bias+0.05]} {
       set Vbias($i) $bias
       incr i
       
       contact name=VSS supply=$bias
       device init
       
       set curr($i) [expr ([contact name=VSS sol=Elec flux] - [contact name=VSS sol=Hole flux])]
       
       sel z=log10(abs(Hole)+1.0)
       plot1d graph=HoleConc
       
       chart graph=IV curve=IV xval=$bias yval=$curr($i)
   }
   
   set winRefresh 1
   if {$winRefresh} {
       plot1d graph=Doping clear
       plot1d graph=HoleConc clear
       chart graph=IV curve=IV clear
   }