PN diode (3D) - Quasi Fermi Method

   ###################################################################
   #                          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
   
   ###################################################################
   #                     Define the Grid                             #
   ###################################################################
   set MeshSpac 0.10
   set Side 2.0
   line x loc=0.0   spac=$MeshSpac   tag=Top
   line x loc=$Side spac=$MeshSpac   tag=Bottom
   
   line y loc=0.0   spac=$MeshSpac   tag=Lside
   line y loc=$Side spac=$MeshSpac   tag=Rside
   
   line z loc=0.0   spac=$MeshSpac   tag=Zlo
   line z loc=1.0   spac=$MeshSpac*2 tag=Zhi
   
   mater add name=Silicon
   region Silicon xlo=Top xhi=Bottom ylo=Lside yhi=Rside zlo=Zlo zhi=Zhi
   init
   
   ###################################################################
   #                  Contacts and Doping                            #
   ###################################################################
   # place the contacts on the corners of the square
   set Contact 0.2
   set buf [ expr $MeshSpac/2]
   contact name=GND Silicon xlo=-$buf xhi=$Contact+$buf ylo=-$buf yhi=$Contact+$buf zlo=-$buf zhi=$buf add
   contact name=VSS Silicon xlo=$Side-$Contact-$buf xhi=$Side+$buf ylo=$Side-$Contact-$buf yhi=$Side+$buf zlo=-$buf zhi=$buf add
   
   sel z=1.0e20*((x+y)<2.0) name=Nd
   sel z=1.0e17 name=Na
   sel z=(Nd-Na) name=Doping
   struct vtk=doping.vtk
   #Plot the doping concentration
   sel z=log10(abs(Doping+1.0))
   window xwin col=2 row=2
   #plot2d zvalue=0.2 levels=05 graph= Doping\ Concentration\ at\ z=0.2
   
   ###################################################################
   #                       Equations                                 #
   ###################################################################
   # Bulk Equations
   DevicePackage
   solution name=Potential nosolve
   solution add name=DevPsi solve negative damp
   solution add name=Qfn solve negative 
   solution add name=Qfp solve negative
   solution add name=Elec const solve val= "1.0e10*exp((DevPsi-Qfn)/0.025)"
   solution add name=Hole const solve val= "1.0e10*exp((Qfp-DevPsi)/0.025)"
   
   set eqnP "$eps * grad(DevPsi) + Doping - Elec + Hole"
   set eqnQfn "ddt(Elec) - $Emob * Elec * grad(Qfn)"
   set eqnQfp "ddt(Hole) - $Hmob * Hole * grad(Qfp)"
   
   pdbSetString Silicon DevPsi Equation $eqnP
   pdbSetString Silicon Qfn    Equation $eqnQfn
   pdbSetString Silicon Qfp    Equation $eqnQfp
   
   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 ohmic.contact {Contact} {
       set vt 0.02558357
       set ni 1.1e10
       pdbSetBoolean $Contact Qfn Flux 1
       pdbSetBoolean $Contact Qfp Flux 1
       pdbSetBoolean $Contact DevPsi Flux 1
       pdbSetBoolean $Contact Qfn Fixed 1
       pdbSetBoolean $Contact Qfp Fixed 1
       pdbSetBoolean $Contact DevPsi Fixed 1
       pdbSetDouble $Contact Qfn Flux.Scale 1.619e-19
       pdbSetDouble $Contact Qfp Flux.Scale 1.619e-19
       pdbSetString $Contact Qfn Equation "Qfn-$Contact"
       pdbSetString $Contact Qfp Equation "Qfp-$Contact"
       pdbSetString $Contact DevPsi Equation "Doping - $ni*exp((DevPsi-Qfn)/$vt) + $ni*exp((Qfp-DevPsi)/$vt)"
       pdbSetString $Contact Equation "1.619e-19 * (Flux_Hole - Flux_Elec)"
   }
   ohmic.contact VSS
   ohmic.contact GND
   
   #Bias contacts
   contact name=VSS voltage supply=0.0
   contact name=GND voltage supply=0.0
   
   ###################################################################
   #                         Solve                                   #
   ###################################################################
   #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
   
   #First DC solve at Equilibrium (0V)
   device
   
   puts "Electron Flux [contact name=VSS sol=Qfn flux]"
   puts "Hole Flux [contact name=VSS sol=Qfp flux]"
   
   #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 < 2.0} {set bias [expr $bias+0.05]} {
       contact name=VSS supply=$bias
       device init
       sel z=Elec
       plot2d zvalue=0.2 levels=100 graph= Elec\ Concentration\ at\ z=0.2
       sel z=Hole
       plot2d zvalue=0.2 levels=100 graph= Hole\ Concentration\ at\ z=0.2
       sel z=DevPsi
       plot2d zvalue=0.2 levels=200 graph= Electrostatic\ Potential\ at\ z=0.2
       set curr($i) [expr ([contact name=VSS sol=Qfn flux] - [contact name=VSS sol=Qfp flux])]
       chart graph=IV curve=IV xval=$bias yval=$curr($i)
       incr i
   }