PMOS capacitor example (1D)

#--------------------
# Create 1D structure
#--------------------
#Grid
line x loc=-0.05 spac=0.01 tag=OxTop
line x loc=0.0 spac=0.01 tag=Top
line x loc=0.5 spac=0.01
line x loc=1.0 spac=0.1 tag=Bottom

mater add name=Silicon
mater add name=Oxide
region Oxide   xlo=OxTop xhi=Top
region Silicon xlo=Top xhi=Bottom
init  

#Contacts
set buf 0.001
contact name=Gate Oxide   xlo=[expr {-0.05-$buf}] xhi=[expr {-0.05+$buf}] add
contact name=GND  Silicon xlo=0.9  xhi=1.1 add 

#-----------------------------
# 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 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 

#-----------------------------
# Declare solution variables
#-----------------------------
DevicePackage
solution add name=DevPsi pde solve negative damp
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)

#solution add name=Ec    const solve val=()
#solution add name=Efn	  const solve val=()
#solution add name=Efp    const solve val=()
#solution add name=Ev    const solve val=() 

#-----------------------------
# Ionized dopant profile
#-----------------------------
sel z=1.0e17 name=Nd
sel z=1.0e20 name=Na
sel z=(Nd-Na) name=Doping 

#plot doping
#sel z=log10(abs(Doping)+1.0)
#plot.1d symb=1 

#-----------------------------
# Bulk Equations
#-----------------------------
#Silicon
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 

#Oxide
set eqnP "$eps * grad(DevPsi)"
pdbSetString Silicon DevPsi Equation $eqnP
pdbSetDouble Silicon DevPsi DampValue $Vt
pdbSetDouble Silicon DevPsi Abs.Error 1.0e-9
 
#-----------------------------
# Interface Equations
#-----------------------------
solution name=DevPsi continuous  

#-----------------------------
# 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 GND

proc MetalContact {Contact} {
       set WF "0.0"
       pdbSetBoolean $Contact DevPsi Flux 1
       pdbSetBoolean $Contact DevPsi Fixed 1
       pdbSetString $Contact DevPsi Equation "DevPsi - $Contact + $WF"
}
MetalContact Gate

#-----------------------------
# Initial Conditions
#-----------------------------
#Bias contacts
contact name=Gate voltage supply=0.0
contact name=GND  voltage supply=0.0

#Initial Guess
sel z= {(Doping>0.0) 
            ?  ( 0.025*log( (Doping+$small) / $ni))
            :  (-0.025*log(-(Doping+$small) / $ni))} 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=Gate sol=Elec flux]"
puts "Hole Flux [contact name=Gate sol=Hole flux]"

#--------------------------------------------------------------
#Plot the equilibrium concentration profiles and Band Diagram
#--------------------------------------------------------------
if {1} {;#concentration
    foreach var {Doping Elec Hole} {
    	sel z=log10(abs($var)+1.0)
    	plot.1d label=$var !cle
    }
} else {;#band diagram
    sel z=DevPsi
}

#note - you can't ramp up a moscap bc the depletion region messes
#       up the Elec and Hole boundary conditions in that region