NMOS IV (2D) - Scharfetter Gummel Method

   ###################################################################
   #                          Constants                              #
   ###################################################################
   set T    300
   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 eox  [expr 3.9 * 8.85418e-14]
   set eps  [expr $esi / $q]
   set epo  [expr $eox / $q]
   set Emob 350.0
   set Hmob 150.0
   set small 1.0e-10
   set Nc 3.2e+19
   set Nv 1.8e+19
   ###################################################################
   #                          Equations                              #
   ###################################################################
   # math 
   math device dim=1 row bcgs ilu tol=1e-30
   math diffuse dim=1 umf none col !scale
   pdbSetDouble Math iterLimit 100 ;# set itteration limit
   
   #add solutions
   DevicePackage ;# Call for Device Simulations
   solution add name=DevPsi pde solve negative damp continuous
   solution add Silicon name=Elec   pde solve !negative
   solution add Silicon name=Hole   pde solve !negative
   
   solution add Oxide   name=Elec const solve val=($small)
   solution add Oxide   name=Hole const solve val=($small)
   
   #set solution equations
   set eqnP "$eps * grad(DevPsi) + Doping - Elec + Hole"
   set eqnE "ddt(Elec) - 400.0 * 0.025 * sgrad(Elec, DevPsi/0.025)"
   set eqnH "ddt(Hole) - 200.0 * 0.025 * sgrad(Hole, -DevPsi/0.025)"
   pdbSetDouble Silicon DevPsi DampValue 0.025
   pdbSetDouble Silicon DevPsi Abs.Error 1.0e-10
   pdbSetString Silicon DevPsi Equation $eqnP
   pdbSetDouble Silicon Elec Abs.Error 1.0e5
   pdbSetString Silicon Elec Equation $eqnE
   pdbSetDouble Silicon Hole Abs.Error 1.0e5
   pdbSetString Silicon Hole Equation $eqnH
   
   set eqnP "$epo * grad(DevPsi) - (1.0e-10)"
   #set eqnP "$epo * grad(DevPsi)"
   pdbSetDouble Oxide DevPsi DampValue $Vt
   pdbSetDouble Oxide DevPsi Abs.Error 1.0e-9
   pdbSetString Oxide DevPsi Equation $eqnP
   
   # Solve For Energy Band Structure 
   pdbSetDouble Silicon Eg 1.1
   pdbSetDouble Silicon Affinity 4.05
   solution add name=Econd solve Silicon const val = "((-[pdbGetDouble Silicon Affinity])-(DevPsi))"
   solution add name=Eval solve Silicon const val = "((-[pdbGetDouble Silicon Affinity])-([pdbGetDouble Silicon Eg])-(DevPsi))"
   set e " Econd - ($Vt) * log( $Nc / Elec ) "
   solution add name=Qfn solve Silicon const val = "($e)"
   set e " Eval + ($Vt) * log( $Nv / Hole ) "
   solution add name=Qfp solve Silicon const val = "($e)"
   
   pdbSetDouble Oxide Eg 9
   pdbSetDouble Oxide Affinity 0.95
   solution add name=Econd solve Oxide const val = "((-[pdbGetDouble Oxide Affinity])-(DevPsi))"
   solution add name=Eval solve Oxide const val = "((-[pdbGetDouble Oxide Affinity])-([pdbGetDouble Oxide Eg])-(DevPsi))"
   
   # Contact Equations
   pdbSetBoolean sub Elec   Fixed 1
   pdbSetBoolean sub Hole   Fixed 1
   pdbSetBoolean sub DevPsi Fixed 1
   pdbSetString sub Hole Equation {Doping - Elec + Hole}
   pdbSetString sub DevPsi Equation {DevPsi + 0.025*log((Hole+1.0e-10)/1.0e10) - sub}
   pdbSetString sub Elec Equation {DevPsi - 0.025*log((Elec+1.0e-10)/1.0e10) - sub}
   pdbSetDouble  sub Elec   Flux.Scale 1.619e-19
   pdbSetDouble  sub Hole   Flux.Scale 1.619e-19
   
   pdbSetBoolean drain Elec   Fixed 1
   pdbSetBoolean drain Hole   Fixed 1
   pdbSetBoolean drain DevPsi Fixed 1
   pdbSetString drain Elec Equation {Doping - Elec + Hole}
   pdbSetString drain Hole Equation {DevPsi + 0.025*log((Hole+1.0e-10)/1.0e10) - drain}
   pdbSetString drain DevPsi Equation {DevPsi - 0.025*log((Elec+1.0e-10)/1.0e10) - drain}
   pdbSetDouble  drain Elec   Flux.Scale 1.619e-19
   pdbSetDouble  drain Hole   Flux.Scale 1.619e-19
   
   pdbSetBoolean src Elec   Fixed 1
   pdbSetBoolean src Hole   Fixed 1
   pdbSetBoolean src DevPsi Fixed 1
   pdbSetString src  Elec Equation {Doping - Elec + Hole}
   pdbSetString src Hole Equation {DevPsi + 0.025*log((Hole+1.0e-10)/1.0e10) - src}
   pdbSetString src DevPsi Equation {DevPsi - 0.025*log((Elec+1.0e-10)/1.0e10) - src}
   pdbSetDouble  src Elec   Flux.Scale 1.619e-19
   pdbSetDouble  src Hole   Flux.Scale 1.619e-19
   
   set WFN -0.0 ;#FB=0.40295, Al=-0.5
   pdbSetBoolean front DevPsi Flux  1
   pdbSetBoolean front DevPsi Fixed 1
   pdbSetString  front DevPsi Equation "DevPsi+$WFN-front"
   
   ###################################################################
   #                          Define Grid                            #
   ###################################################################
   #line x loc=-0.02 spac=0.01 tag=M
   line x loc=-0.01 spac=0.001 tag=TopOx
   line x loc=0.0 spac=0.001 tag=TopSi
   line x loc=0.025 spac=0.005
   line x loc=0.04 spac=0.01 tag=Jn
   line x loc=0.6 spac=0.1 tag=Bottom
   
   line y loc=0.0 spac=0.1 tag=Left
   line y loc=0.1 spac=0.01 tag=GL
   line y loc=0.2 spac=0.01
   line y loc=0.3 spac=0.01 tag=GR
   line y loc=0.4 spac=0.1 tag=Right
   
   #add material and region
   mater add name=Oxide
   mater add name=Silicon
   
   region Oxide   xlo=TopOx xhi=TopSi  ylo=GL yhi=GR
   region Silicon xlo=TopSi xhi=Bottom ylo=Left yhi=Right
   init 
   
   contact name=drain silicon xlo=-0.005 xhi=0.005 ylo=0.32 yhi=0.51 add supply=0.0
   contact name=front oxide xlo=-0.015 xhi=-0.005 ylo=0.09 yhi=4.41 add supply=0.0
   contact name=sub silicon xlo=0.595 xhi=0.605 ylo=-0.05 yhi=0.55 add supply=0.0
   contact name=src silicon xlo=-0.005 xhi=0.005 ylo=-0.01 yhi=0.08 add supply=0.0
   
   ###################################################################
   #                       Doping Profile                            #
   ###################################################################
   set buff 1.0e10
   sel z=(1.0e21*((y>0.32)?1.0:(exp(-(y-0.32)*(y-0.32)/(0.0002)))*(exp(-(x*x)/(0.0002))))) name=DrainDoping
   sel z=(1.0e21*((y<0.08)?1.0:(exp(-(y-0.08)*(y-0.08)/(0.0002)))*(exp(-(x*x)/(0.0002))))) name=SourceDoping
   sel z=(5e16+$buff) name=NA
   
   sel z=DrainDoping+SourceDoping-NA name=Doping
   
   
   ###################################################################
   #                       Solve IV Curve                            #
   ###################################################################
   set xval 0.0001
   set yval 0.2
   #setup window to prepare to plot
   window xwin col=3
   
   for {set Vgs 0} {$Vgs<=8} {set Vgs [expr $Vgs+1]} {
       # set initial conditions at each Vgs
       sel z=0.5*(Doping+sqrt(Doping*Doping+4.0e20))/1.0e10 name=arg
       sel z=0.025*log(abs(arg)+1) name=DevPsi
       sel z=1.0e10*exp(DevPsi/0.025)  name=Elec
       sel z=1.0e10*exp(-DevPsi/0.025) name=Hole
       # Solve for proper electrostatics at each Vgs
       contact name=front supply=$Vgs
       device  ;# Call device 3 times to ensure convergence
       device
       device
       for {set Vds 0.0} {$Vds<10.0} {set Vds [expr $Vds+0.2]} {
       # Sweep the Drain Voltage
           contact name=drain supply=$Vds
           device
           device
           
           sel z=Econd
       plot2d graph=shades clear
       plot2d bound levels=100 graph=shades xmin=0.001 xmax=0.25 ;#line=20 
           plot3d graph=persp azimuth=140 altitude=50 xmin=0.001 xmax=0.15
           set cur [expr abs([contact name=drain sol=Elec flux] - [contact name=drain sol=Hole flux])]
           chart graph=IV curve=Vgs$Vgs xval=$Vds yval=$cur
       }
   }
   
   puts "Press Enter to Continue"
   gets stdin