Resistor (1D) - full deck

#########   1D Resistor Deck  ###############

#---------------------
# Create 1D struture 
#---------------------
#Grid
line x loc=0.0 spac=0.1 tag=Top
line x loc=1.0 spac=0.1 tag=Bottom
mater add name=Silicon
region Silicon xlo=Top xhi=Bottom
init

#Contacts
contact name=VSS Silicon xlo=-0.1 xhi=0.1 add
contact name=GND 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
#-------------------------
set Nd 1.0e19
set Na 1.0e15
set Doping [expr {$Nd-$Na}]

#---------------------
# 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)"
pdbSetDouble Silicon DevPsi DampValue $Vt
pdbSetString Silicon DevPsi Equation $eqnP
pdbSetString Silicon Elec   Equation $eqnE
pdbSetString Silicon Hole   Equation $eqnH 

#---------------------
# Contact Equations
#---------------------
proc OhmicContact {Contact} {
    global Vt ni Nd Na
    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 Voltage on the Contacts
contact name=VSS voltage supply=0.0
contact name=GND voltage supply=0.0

#Initial Guess at Zero Bias
sel z=(-$Vt*log(-($Doping+$small)/$ni)) name=DevPsi
sel z=$ni*exp(DevPsi/0.025) name=Elec
sel z=$ni*exp(-DevPsi/0.025) name=Hole

#DC Solve at Zero Bias
device

#---------------------------------
# DC solve / plot I-V as output
#---------------------------------
set Win [CreateGraphWindow]
set bias 0.0
for {set bias 0.0} {$bias < 1.01} {set bias [expr $bias+0.1]} {
	contact name=VSS supply = $bias
	device
	set cur [expr abs([contact name=VSS sol=Elec flux] - [contact name=VSS sol=Hole flux])]
	AddtoLine $Win I $bias $cur
}