View source for PN diode (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*(x<0.5) 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
    }