###################################################################
# 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=$MeshSpac-$buf xhi=$Contact+$buf ylo=$MeshSpac-$buf yhi=$Contact+$buf zlo=-$buf zhi=$buf add
contact name=VSS Silicon xlo=$Side-$Contact-$buf xhi=$Side+$buf-$MeshSpac ylo=$Side-$Contact-$buf yhi=$Side+$buf-$MeshSpac zlo=-$buf zhi=$buf add
sel z=1.0e20 name=Nd
sel z=1.0e17 name=Na
sel z=(Nd-Na) name=Doping
#Plot the doping concentration
sel z=log10(abs(Doping+1.0))
window xwin col=2
plot2d zvalue=0.2 levels=05 graph= Doping\ Concentration\ at\ z=0.2
###################################################################
# Equations #
###################################################################
solution add name=DevPsi pde solve negative damp
solution add name=Elec pde solve !negative
solution add name=Hole pde solve !negative
solution add name=Circuit solve negative damp circuit const val=0.0
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
#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=Elec flux]"
puts "Hole Flux [contact name=VSS sol=Hole 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 < 1.2} {set bias [expr $bias+0.05]} {
contact name=VSS supply=$bias
device init
set curr($i) [expr ([contact name=VSS sol=Elec flux] - [contact name=VSS sol=Hole flux])]
chart graph=IV curve=IV xval=$bias yval=$curr($i)
incr i
}
gets stdin