NMOSFET (2D) - full deck

# Define solution variables
DevicePackage
solution name=Potential nosolve
solution add name=DevPsi solve  negative 
solution add name=Elec   solve !negative
solution add name=Hole   solve !negative

# Create 2D structure mesh

# Horizontal lines
line x loc=0.00 spac=0.01   tag=TopOx
line x loc=0.03 spac=0.001  tag=TopSi
line x loc=0.04 spac=0.0025
line x loc=0.06 spac=0.01
line x loc=0.09 spac=0.005
line x loc=0.10 spac=0.01
line x loc=0.12 spac=0.025
line x loc=0.19 spac=0.05
line x loc=2.00 spac=0.1    tag=Bottom

# Vertical lines
line y loc=0.00 spac=0.05   tag=Left
line y loc=0.10 spac=0.01
line y loc=0.15 spac=0.01
line y loc=0.25 spac=0.007  tag=GLE
line y loc=0.35 spac=0.05
line y loc=1.25 spac=0.1
line y loc=2.15 spac=0.05
line y loc=2.25 spac=0.007  tag=GRE
line y loc=2.35 spac=0.01
line y loc=2.40 spac=0.01
line y loc=2.50 spac=0.05   tag=Right

# Regions
region oxide   xlo=TopOx xhi=TopSi  ylo=GLE yhi=GRE
region silicon xlo=TopSi xhi=Bottom ylo=Left yhi=Right
init

# Create contacts

# Define contacts place and size
contact name=gate      oxide   xlo=0.00 xhi=0.00 ylo=0.25 yhi=2.25 add depth=1.0 width=1.0
contact name=src       silicon xlo=0.03 xhi=0.03 ylo=0.00 yhi=0.25 add depth=1.0 width=1.0
contact name=drain     silicon xlo=0.03 xhi=0.03 ylo=2.25 yhi=2.50 add depth=1.0 width=1.0
contact name=substrate silicon xlo=2.00 xhi=2.00 ylo=0.00 yhi=2.50 add depth=1.0 width=1.0

# Define initial bias
contact name=gate      voltage supply=0.0
contact name=src       voltage supply=0.0
contact name=drain     voltage supply=0.0
contact name=substrate voltage supply=0.0

# Save structure
struct outf=nMosfet.str

# Define doping
set buff 1.0e10
sel z=((1.0e20)*(x<0.11)*((y<0.25)+(y>2.25))) name=ND
sel z=(1.0e15*(x>0.03)+$buff) name=NA
sel z=ND-NA name=Doping

# Plot structure and contacts. The !cle allows features to be added to the current graphic.
#plot.2d grid
plot.2d bound fill        
#plot.2d contact=gate      !cle
#plot.2d contact=src       !cle
#plot.2d contact=drain     !cle
#plot.2d contact=substrate !cle
sel z=Doping
contour val=0.0
sel z=log10(abs(Doping))
for {set i 15} {$i <= 20.0} {incr i} {
    contour val=$i
}

# Define constants
set T    300
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

# Define Poisson's solution and continuity equations (silicon)
# Poisson's equation
set eqnP "$eps * grad(DevPsi) + Doping - Elec + Hole"
pdbSetDouble Silicon DevPsi DampValue $Vt
pdbSetDouble Silicon DevPsi Abs.Error 1.0e-9
pdbSetString Silicon DevPsi Equation  $eqnP
# Electron continuity
set eqnE "ddt(Elec) - ($Emob) * $Vt * sgrad(Elec,  DevPsi/$Vt)"
pdbSetDouble Silicon DevPsi Abs.Error 1.0e9
pdbSetString Silicon Elec   Equation  $eqnE
# Hole continuity
set eqnH "ddt(Hole) - ($Hmob) * $Vt * sgrad(Hole, -DevPsi/$Vt)"
pdbSetDouble Silicon DevPsi Abs.Error 1.0e9
pdbSetString Silicon Hole   Equation  $eqnH 
# 
pdbSetDouble Silicon DevPsi DampValue $Vt

# Define Poisson's solution and continuity equations (oxide)
# Poisson's equation
set eqnP "$epo * grad(DevPsi) - (1.0e-10)"
pdbSetDouble Oxide DevPsi DampValue $Vt
pdbSetDouble Oxide DevPsi Abs.Error 1.0e-9
pdbSetString Oxide DevPsi Equation $eqnP

# Define ohmic contact equations
proc ohmic.contact {Contact} {
	set vt 0.02558357
	set ni 1.1e10
	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"
}
ohmic.contact src
ohmic.contact drain
ohmic.contact substrate

# Define metal contact and Work-Function Difference
pdbSetDouble Aluminum 0.0
proc MetalContact {Contact Mat} {
	set WFN [pdbGetDouble $Mat WFN] 
	pdbSetBoolean $Contact DevPsi Flux  1
	pdbSetBoolean $Contact DevPsi Fixed 1
	pdbSetString  $Contact DevPsi Equation "DevPsi + $WFN - $Contact"; # Efm+WFN=El=-DevPsi, Efm=-$Contact
}
MetalContact gate Aluminum

# Define continuous interface equtions between different materials
solution name=DevPsi continuous

# Initial guess procedure, assumes charge neutrality
proc InitialGuess {Doping} {
	sel z= {(Doping>0.0)  
		?  ( 0.025*log(abs(Doping+1.0e-10)/1.0e10))
		:  (-0.025*log(abs(Doping+1.0e-10)/1.0e10))} name=DevPsi
	sel z=1.0e10*exp(DevPsi/0.025)  name=Elec
	sel z=1.0e10*exp(-DevPsi/0.025) name=Hole
} 
InitialGuess Doping

# Ramp drain voltage up first
set delta 0.01
for {set bias 0.0} {$bias <= 0.050} {set bias [expr {$bias + $delta}]} {
	puts \n
	contact name=drain voltage supply=$bias
	puts "Ramping drain to $bias"
	device
}

# Ramp gate up
set Win_IdVg [CreateGraphWindow]
set delta 0.01
for {set bias 0.0} {$bias <= 1.0} {set bias [expr {$bias + $delta}]} {
	puts \n
	contact name=gate voltage supply=$bias
	puts "Ramping gate to $bias"
	device

	# get current from the drain (choose regular or log)
	set cur [expr abs([contact name=drain sol=Elec flux] - [contact name=drain sol=Hole flux])]
	# set cur log10($cur+1e-10)
	AddtoLine $Win_IdVg IdVg $bias $cur; #channel, name, x-val, y-val
}