Some questions about 2D-MOS example

1

Hello everyone, I am study the example of 2D MOS example. Could any help me in understanding the following scripts:

node_model(name="Donors",    device=device, region="gate", equation="%(gate_doping)1.15e + 1" % (mydict))
node_model(name="Acceptors", device=device, region="gate", equation="1")
node_model(name="NetDoping",    device=device, region="gate", equation="Donors - Acceptors")

node_model(name="DrainDoping",  device=device, region="bulk", equation="0.25*%(drain_doping)1.15e*erfc((x-%(x_gate_left)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
node_model(name="SourceDoping", device=device, region="bulk", equation="0.25*%(source_doping)1.15e*erfc(-(x-%(x_gate_right)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
node_model(name="BodyDoping",   device=device, region="bulk", equation="0.5*%(body_doping)1.15e*erfc(-(y-%(y_bulk_bottom)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
node_model(name="Donors",       device=device, region="bulk", equation="DrainDoping + SourceDoping + 1")
node_model(name="Acceptors",    device=device, region="bulk", equation="%(bulk_doping)1.15e + BodyDoping" % mydict)
node_model(name="NetDoping",    device=device, region="bulk", equation="Donors - Acceptors")

My questions are:
1、why equation are expressed as

equation="%(gate_doping)1.15e + 1" % (mydict)

here it is different from “%s ” %string and why 1.15e is here and whta does it mean？

2、why “*%” is in the equation and what does it mean?
3、why equation for donors is “1”
4、where we get some information about how to give the eqution for doping

node_model(name="DrainDoping",  device=device, region="bulk", equation="0.25*%(drain_doping)1.15e*erfc((x-%(x_gate_left)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)

Sorry i can not understand it

2

%(gate_doping)1.15e is a string interpolation method ensuring that number in the script show the full floating point representation of the value for mydict['gate_doping']. Python has many ways of doing this, and this format is an older method.

It is one way of specifying a minimum doping density. Sometimes if the density is 0, then there can be floating point errors or poor convergence with the simulator.

For the physics scripts, the most important quantity is the NetDoping, which goes into the Potential equation. The Acceptors and Donors may be required individually for some of the advanced mobility models.

The erfc function can be used to create an analytical doping profile that replicates doping from a semiconductor process.

1 Like
3

@Juan，Thanks.

But I searched erfc function in the Devsim mannual and only a very simple description is listed as below

erfc(exp1) complementary error function

Do you have any other more reference for this function？ I am wondering how to decide the input parameter of this function when creating an analytical doping profile.

4

Hi @caihng

Doing a quick search for:
erfc semiconductor processing

I found this book:
FUNDAMENTALS OF SEMICONDUCTOR MANUFACTURING AND PROCESS CONTROL

https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.465.9904&rep=rep1&type=pdf

which uses the erfc in a diffusion model. Please note the erfc function may have slightly different definitions in the literature.

You may choose to use other functions, such as the step function. It is also possible to use a numerical doping profile, if you would like.

5

Thanks @Juan ，I will have a look at it :smiley:

6

@Juan Hello Juan，i want to import my Gmsh file into Devsim, But it seems that MeshFormat 4.1 0 8 is not supported. Do you know what format i need to use ?

error: line: 3: ERROR: MeshFormat 4.1 0 8 not supported
7

@Juan, I tried MeshFormat 2.0，It seems ok

1 Like
8

@Juan Hello, Juan, I am trying to do a 2D MOS simulation based on the example and some error occured as below:

number of equations 10535
while evaluating node model Drain_CTnodemodel on Device: mos2d on Region: Drain
There was a Divide-by-zero floating point exception while evaluating log(((5.000000000000000e-01 * abs((-NetDoping + pow((pow(NetDoping,2) + (4 * pow(n_i,2))),5.000000000000000e-01)))) * pow(n_i,(-1))))
Invalid argument while evaluating function ifelse
while evaluating model Drain_CTnodemodel: expression (-Drain_CT_bias + Potential + (ifelse ((NetDoping > 0), (-V_t * log(((5.000000000000000e-01 * abs((NetDoping + pow((pow(NetDoping,2) + (4 * pow(n_i,2))),5.000000000000000e-01)))) * pow(n_i,(-1))))), (V_t * log(((5.000000000000000e-01 * abs((-NetDoping + pow((pow(NetDoping,2) + (4 * pow(n_i,2))),5.000000000000000e-01)))) * pow(n_i,(-1)))))))) evaluates to invalid
Traceback (most recent call last):

  File "D:\06 Devsim\projects\SiC 2dMOSFET\gmsh_mos2d.py", line 57, in <module>
    solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)

error: DEVSIM FATAL

My scripts are as below:

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create
device = "mos2d"
device_regions=("Drain","Ndrift","Pwell","Nwell","Pbase","Pshield")
oxide_regions=("Gox",)
gate_regions = ("Poly",)
regions = ("Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Gox","Poly")
interfaces = ("Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface")

###### Create Solutions for Potential and set region materials ####
for i in regions:
    CreateSolution(device, i, "Potential")

for i in device_regions:
    SetSiCParameters(device, i, 300)
    CreateSiliconPotentialOnly(device, i)

for i in oxide_regions:
    SetOxideParameters(device, i, 300)
    CreateOxidePotentialOnly(device, i, "log_damp")

for i in gate_regions:
    SetPolyParameters(device, i, 300)
    CreatePolyPotentialOnly(device, i)

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
    tmp = get_region_list(device=device, contact=i)
    r = tmp[0]
    print("%s %s" % (r, i))  #print region and contact name
    CreateSiliconPotentialOnlyContact(device, r, i)
    set_parameter(device=device, name=GetContactBiasName(i), value=0.0)

for i in interfaces:
    CreateSiliconOxideInterface(device, i)  #continuous potential at interface

solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
print("123456789")
solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
##write_devices -file gmsh_mos2d_potentialonly.flps -type floops
write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

for i in device_regions:
    CreateSolution(device, i, "Electrons")
    CreateSolution(device, i, "Holes")
    set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
    set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
    CreateSiliconDriftDiffusion(device, i, "mu_n", "mu_p")

for c in contacts:
    tmp = get_region_list(device=device, contact=c)
    r = tmp[0]
    CreateSiliconDriftDiffusionAtContact(device, r, c)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)

for r in device_regions:
    node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data
#


for r in device_regions:
    element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
    element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
    element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)
#
rampbias(device, "gate",  0.5, 0.5, 0.001, 100, 1e-10, 1e30, printAllCurrents)
rampbias(device, "drain", 0.5, 0.1, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

for scripts such as “simple physics”、“simple dd” I nearly do no change except some additional equations

and My script to load the Gmsh file is as below:

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from devsim import *
import pandas as pd
device ="mos2d"

device_width    =1.0e-4
gate_width      =1.0e-5
diffusion_width =0.4

air_thickness       =1e-7
oxide_thickness     =1e-5
gate_thickness      =1e-5
device_thickness    =1e-4
diffusion_thickness =1e-5

x_diffusion_decay  =1e-20
y_diffusion_decay  =1e-10

########## Doping Concentration Inputs #########
#p doping
Pwell_doping       =3e17
Gate_doping        =1e20
Pbase_doping       =3e16
Pshield_doping     =5e17
#n doping
Nwell_doping       =1e18
Ndrift_doping      =1e16
Nsub_doping        =1e19
########## Doping Concentration Inputs #########


# region_list = ["Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Poly"]
interface_list = ["Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface"]

material = ["SiC", "Poly", "Oxide"]
SiC_region = ["Drain", "Ndrift", "Pwell", "Nwell", "Pbase", "Pshield"]
Poly_region = ["Poly"]
Oxide_region = ["Gox"]
# Ma_region = pd.DataFrame({"SiC":SiC_region,"Poly":Poly_region, "Oxide":Oxide_region})
create_gmsh_mesh(file="SiC_MOSFET_2d_Gmsh.msh", mesh="mos2d")

######### Create regions from gmsh, region names in gmsh file and devsim keep the same ##
for i in material:
    for j in eval(i+"_region"):
        add_gmsh_region(mesh="mos2d", gmsh_name= j, region= j, material= i)

# add_gmsh_region(mesh="mos2d", gmsh_name="Drain", region="Drain", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Ndrift", region="Ndrift", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Pwell",  region="Pwell", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Nwell",  region="Nwell", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Pbase",  region="Pbase", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Pshield",  region="Pshield", material="SiC")
# add_gmsh_region(mesh="mos2d", gmsh_name="Poly",  region="Poly", material="Poly")
# add_gmsh_region(mesh="mos2d", gmsh_name="Gox",  region="Gox", material="Oxide")

######### Create contacts from gmsh ##
add_gmsh_contact(mesh="mos2d", gmsh_name="Gate_Contact",   region="Poly", name="Gate_CT", material="metal")
add_gmsh_contact(mesh="mos2d", gmsh_name="Drain_Contact",  region="Drain", name="Drain_CT", material="metal")
# add_gmsh_contact(mesh="mos2d", gmsh_name="Source_Contact", region="Pwell", name="Source_CT", material="metal")

##### add interface ############
interfaces = ["Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface"]
for i in interfaces:
    add_gmsh_interface (mesh="mos2d", gmsh_name= i, region0= i.split("_")[0], region1= i.split("_")[1], name= i)
# add_gmsh_interface (mesh="mos2d", gmsh_name="gate_oxide_interface", region0="gate", region1="oxide", name="gate_oxide")




#interface not decalred in Gmsh file
# add_gmsh_interface (mesh="mos2d", gmsh_name="gate_oxide_interface", region0="gate", region1="oxide", name="gate_oxide")
# add_gmsh_interface (mesh="mos2d", gmsh_name="bulk_oxide_interface", region0="bulk", region1="oxide", name="bulk_oxide")
finalize_mesh(mesh="mos2d")
# create_device(mesh="mos2d", device="mos2d")
create_device(mesh="mos2d", device = device)


#### all variable substitutions are immediate, since they are locked into the mesh
mydict = {}
mydict["Pwell"] = Pwell_doping
mydict["Poly"] = Gate_doping
mydict["Pbase"] = Pbase_doping
mydict["Pshield"] = Pshield_doping

mydict["Nwell"] = Nwell_doping
mydict["Ndrift"] = Ndrift_doping
mydict["Drain"] = Nsub_doping

region_list = ["Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Poly"]
# mydict["x_gate_left"] = x_gate_left
# mydict["x_gate_right"] = x_gate_right
# mydict["x_diffusion_decay"] = x_diffusion_decay
# mydict["y_diffusion"] = y_diffusion
# mydict["y_bulk_bottom"] = y_bulk_bottom
# mydict["y_diffusion_decay"] = y_diffusion_decay

######## for different area, doping profile and Net doping is required #######
for i in region_list:
    print("doping concentration in region %s is %.2e." %(i, mydict[i]))
    node_model(name="Donors", device=device, region=i, equation="%s + 1" %mydict[i])
    node_model(name="Acceptors", device=device, region=i, equation="1")
    node_model(name="NetDoping",    device=device, region=i, equation="Donors - Acceptors")
######## for different area, doping profile and Net doping is required #######



# node_model(name="Donors",    device=device, region="gate", equation="%(gate_doping)1.15e + 1" % (mydict))
# #node_model(name="Donors",    device=device, region="gate", equation="%s + 1" % mydict["gate_doping"])
# node_model(name="Acceptors", device=device, region="gate", equation="1")
# node_model(name="NetDoping",    device=device, region="gate", equation="Donors - Acceptors")

# node_model(name="DrainDoping",  device=device, region="bulk", equation=
#            "0.25*%(drain_doping)1.15e*erfc((x-%(x_gate_left)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="SourceDoping", device=device, region="bulk", equation=
#            "0.25*%(source_doping)1.15e*erfc(-(x-%(x_gate_right)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="BodyDoping",   device=device, region="bulk", equation=
#            "0.5*%(body_doping)1.15e*erfc(-(y-%(y_bulk_bottom)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="Donors",       device=device, region="bulk", equation="DrainDoping + SourceDoping + 1")
# node_model(name="Acceptors",    device=device, region="bulk", equation="%(bulk_doping)1.15e + BodyDoping" % mydict)
# node_model(name="NetDoping",    device=device, region="bulk", equation="Donors - Acceptors")
write_devices(file="gmsh_mos2d_out.msh")

Could you please help me to find the cause of this error and give some suggestions?

9

@Juan Hello Juan, i added the interface condition, some other errors occur as below:

Gox_Poly_interface
Gox_Nwell_interface
Gox_Pbase_interface
Gox_Ndrift_interface
Gox_Pshield_interface
Ndrift_Pshield_interface
Pwell_Ndrift_interface
Pwell_Nwell_interface
Pwell_Pbase_interface
Ndrift_Drain_interface
Nwell_Pshield_interface
123456789
number of equations 10568
On interface Ndrift_Drain_interface
Could not find equation ElectronContinuityEquation on region0 Ndrift for interface equation ElectronContinuityEquation
Could not find equation ElectronContinuityEquation on region1 Drain for interface equation ElectronContinuityEquation
Traceback (most recent call last):

  File "D:\06 Devsim\projects\SiC 2dMOSFET\gmsh_mos2d.py", line 74, in <module>
    solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)

error: DEVSIM FATAL

To be more clear, I paste my script here:
gmsh_2d.py

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from devsim.python_packages.simple_physics_SiC import *
from devsim.python_packages.simple_physics_Poly import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create
device = "mos2d"
device_regions=("Drain","Ndrift","Pwell","Nwell","Pbase","Pshield")
oxide_regions=("Gox",)
gate_regions = ("Poly",)
regions = ("Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Gox","Poly")

# SiC_interface = ("Ndrift_Pshield_interface","Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
#                  "Ndrift_Drain_interface","Nwell_Pshield_interface")
# Oxide_interface = ("Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
#                    "Gox_Ndrift_interface","Gox_Pshield_interface")
interfaces = ("Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface")

###### Create Solutions for Potential and set region materials ####
for i in regions:
    CreateSolution(device, i, "Potential")

for i in device_regions:
    SetSiCParameters(device, i, 300)
    CreateSiCPotentialOnly(device, i)

for i in oxide_regions:
    SetOxideParameters(device, i, 300)
    CreateOxidePotentialOnly(device, i, "log_damp")

for i in gate_regions:
    SetPolyParameters(device, i, 300)
    CreatePolyPotentialOnly(device, i)

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
    tmp = get_region_list(device=device, contact=i)
    r = tmp[0]
    print("%s %s" % (r, i))  #print region and contact name
    if "Poly" in r:
        CreatePolyPotentialOnlyContact(device, r, i)
        set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
    elif "Poly" not in r:
        CreateSiCPotentialOnlyContact(device, r, i)
        set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
    # set_parameter(device=device, name=GetContactBiasName(i), value=0.0)

#### Set up interface ####
for i in interfaces:
    print(i)
    if "Gox" in i:
        CreatePolyOxideInterface(device, i)  #continuous potential at interface
    if "Gox" not in i:
        CreateSiCSiCInterface(device, i)  #continuous potential at interface
    
print("123456789")
solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
print("123456789")
# DriftDiffusionInitialSolution(device, region)
###
### Drift diffusion simulation at equilibrium
###
solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
##write_devices -file gmsh_mos2d_potentialonly.flps -type floops
write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

for i in device_regions:
    CreateSolution(device, i, "Electrons")
    CreateSolution(device, i, "Holes")
    set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
    set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
    CreatePolyDriftDiffusion(device, i, "mu_n", "mu_p")

for c in contacts:
    tmp = get_region_list(device=device, contact=c)
    r = tmp[0]
    CreatePolyDriftDiffusionAtContact(device, r, c)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)

for r in device_regions:
    node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data
#

for r in device_regions:
    element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
    element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
    element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)
#
rampbias(device, "gate",  0.5, 0.5, 0.001, 100, 1e-10, 1e30, printAllCurrents)
rampbias(device, "drain", 0.5, 0.1, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

gmsh create:

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from devsim import *
import pandas as pd
device ="mos2d"

# device_width    =1.0e-4
# gate_width      =1.0e-5
# diffusion_width =0.4

# air_thickness       =1e-7
# oxide_thickness     =1e-5
# gate_thickness      =1e-5
# device_thickness    =1e-4
# diffusion_thickness =1e-5

# x_diffusion_decay  =1e-20
# y_diffusion_decay  =1e-10

########## Doping Concentration Inputs #########
#p doping
Pwell_doping       =3e17
Gate_doping        =1e20
Pbase_doping       =3e16
Pshield_doping     =5e17
#n doping
Nwell_doping       =1e18
Ndrift_doping      =1e16
Nsub_doping        =1e19
########## Doping Concentration Inputs #########


# region_list = ["Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Poly"]
interface_list = ["Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface"]

material = ["SiC", "Poly", "Oxide"]
SiC_region = ["Drain", "Ndrift", "Pwell", "Nwell", "Pbase", "Pshield"]
Poly_region = ["Poly"]
Oxide_region = ["Gox"]
# Ma_region = pd.DataFrame({"SiC":SiC_region,"Poly":Poly_region, "Oxide":Oxide_region})
create_gmsh_mesh(file="SiC_MOSFET_2d_Gmsh.msh", mesh="mos2d")

######### Create regions from gmsh, region names in gmsh file and devsim keep the same ##
for i in material:
    for j in eval(i+"_region"):
        add_gmsh_region(mesh="mos2d", gmsh_name= j, region= j, material= i)


######### Create contacts from gmsh ##
add_gmsh_contact(mesh="mos2d", gmsh_name="Gate_Contact",   region="Poly", name="Gate_CT", material="metal")
add_gmsh_contact(mesh="mos2d", gmsh_name="Drain_Contact",  region="Drain", name="Drain_CT", material="metal")
add_gmsh_contact(mesh="mos2d", gmsh_name="Source_Contact", region="Pwell", name="Source_CT", material="metal")

##### add interface ############
interfaces = ["Gox_Poly_interface", "Gox_Nwell_interface","Gox_Pbase_interface",
              "Gox_Ndrift_interface","Gox_Pshield_interface","Ndrift_Pshield_interface",
              "Pwell_Ndrift_interface","Pwell_Nwell_interface","Pwell_Pbase_interface",
              "Ndrift_Drain_interface","Nwell_Pshield_interface"]
for i in interfaces:
    add_gmsh_interface (mesh="mos2d", gmsh_name= i, region0= i.split("_")[0], region1= i.split("_")[1], name= i)


finalize_mesh(mesh="mos2d")
# create_device(mesh="mos2d", device="mos2d")
create_device(mesh="mos2d", device = device)


#### all variable substitutions are immediate, since they are locked into the mesh
mydict = {}
mydict["Pwell"] = Pwell_doping
mydict["Poly"] = Gate_doping
mydict["Pbase"] = Pbase_doping
mydict["Pshield"] = Pshield_doping

mydict["Nwell"] = Nwell_doping
mydict["Ndrift"] = Ndrift_doping
mydict["Drain"] = Nsub_doping

region_list = ["Drain","Ndrift","Pwell","Nwell","Pbase","Pshield","Poly"]
# mydict["x_gate_left"] = x_gate_left
# mydict["x_gate_right"] = x_gate_right
# mydict["x_diffusion_decay"] = x_diffusion_decay
# mydict["y_diffusion"] = y_diffusion
# mydict["y_bulk_bottom"] = y_bulk_bottom
# mydict["y_diffusion_decay"] = y_diffusion_decay

######## for different area, doping profile and Net doping is required #######
for i in region_list:
    print("doping concentration in region %s is %.2e." %(i, mydict[i]))
    node_model(name="Donors", device=device, region=i, equation="%s + 1" %mydict[i])
    node_model(name="Acceptors", device=device, region=i, equation="1")
    node_model(name="NetDoping",    device=device, region=i, equation="Donors - Acceptors")
######## for different area, doping profile and Net doping is required #######



# node_model(name="Donors",    device=device, region="gate", equation="%(gate_doping)1.15e + 1" % (mydict))
# #node_model(name="Donors",    device=device, region="gate", equation="%s + 1" % mydict["gate_doping"])
# node_model(name="Acceptors", device=device, region="gate", equation="1")
# node_model(name="NetDoping",    device=device, region="gate", equation="Donors - Acceptors")

# node_model(name="DrainDoping",  device=device, region="bulk", equation=
#            "0.25*%(drain_doping)1.15e*erfc((x-%(x_gate_left)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="SourceDoping", device=device, region="bulk", equation=
#            "0.25*%(source_doping)1.15e*erfc(-(x-%(x_gate_right)1.15e)/%(x_diffusion_decay)1.15e)*erfc((y-%(y_diffusion)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="BodyDoping",   device=device, region="bulk", equation=
#            "0.5*%(body_doping)1.15e*erfc(-(y-%(y_bulk_bottom)1.15e)/%(y_diffusion_decay)1.15e)" % mydict)
# node_model(name="Donors",       device=device, region="bulk", equation="DrainDoping + SourceDoping + 1")
# node_model(name="Acceptors",    device=device, region="bulk", equation="%(bulk_doping)1.15e + BodyDoping" % mydict)
# node_model(name="NetDoping",    device=device, region="bulk", equation="Donors - Acceptors")
write_devices(file="gmsh_mos2d_out.msh")

simple physics of SiC and Poly gate

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
# celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
# chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
celec_model = "(0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"   #revised
chole_model = "(0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"  #revised

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
T      = 300 # K
# eps_si = 11.1
eps_ox = 3.9
eps_sic = 9.7
# eps_Poly = 11.1
# eps = {"SiC": 9.7,"Poly": 11.1,"Gox":3.9}
# TODO: make this temperature dependent
# n_i    = 1.0e10 # #/cm^3  # for silicon
n_i    = 1e-8 # #/cm^3  #for SiC
# n_i = {"SiC":1e-8,"Poly":1e10}
# constant in our approximation
# mu_n   = 400 #for Si
# mu_p   = 200 #for Si
mu_n   = 950 # for SiC
mu_p   = 125 # for SiC
# mu_n = {"SiC":950,"Poly":400}  #注意浓度、电场的调制等效应
# mu_p = {"SiC":125,"Poly":200}
# taun = {"SiC":2.5e-6,"Poly":2.5e-6}
# taup = {"SiC":0.5e-6,"Poly":0.5e-6}
taun = 2.5e-6 # for SiC
taup = 0.5e-6 # for SiC


def GetContactBiasName(contact):
    return "{0}_bias".format(contact)

def GetContactNodeModelName(contact):
    return "{0}nodemodel".format(contact)

def PrintCurrents(device, contact):
    '''
       print out contact currents
    '''
    # TODO add charge
    contact_bias_name = GetContactBiasName(contact)
    electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
    hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
    total_current   = electron_current + hole_current                                        
    voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
    print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
    conteq="Permittivity*ElectricField"
    contact_bias_name  = GetContactBiasName(contact)
    contact_model_name = GetContactNodeModelName(contact)
    eq = "Potential - {0}".format(contact_bias_name)
    CreateContactNodeModel(device, contact, contact_model_name, eq)
    CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

    #TODO: make everyone use dfield
    if not InEdgeModelList(device, region, contactcharge_edge):
        CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
        CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

    contact_equation(device=device, contact=contact, name="PotentialEquation",
                     node_model=contact_model_name, edge_charge_model= contactcharge_edge)




#####
##### Constants
#####
def SetOxideParameters(device, region, T):
    '''
      Sets physical parameters
    '''
    set_parameter(device=device, region=region, name="Permittivity",   value=eps_ox * eps_0)
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)

def SetSiCParameters(device, region, T):
    '''
      Sets physical parameters assuming constants
    '''
    #### TODO: make T a free parameter and T dependent parameters as models
    set_parameter(device=device, region=region, name="Permittivity",   value=eps_sic * eps_0) #revise for sic
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)
    set_parameter(device=device, region=region, name="n_i",            value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="T",              value=T)
    set_parameter(device=device, region=region, name="kT",             value=k * T)
    set_parameter(device=device, region=region, name="V_t",            value=k*T/q)
    set_parameter(device=device, region=region, name="mu_n",           value=mu_n) #revise for sic
    set_parameter(device=device, region=region, name="mu_p",           value=mu_p) #revise for sic
    #default SRH parameters
    set_parameter(device=device, region=region, name="n1", value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="p1", value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="taun", value=taun) #revise for sic
    set_parameter(device=device, region=region, name="taup", value=taup) #revise for sic


def CreateSiCPotentialOnly(device, region):
    '''
      Creates the physical models for a Silicon region
    '''
    if not InNodeModelList(device, region, "Potential"):
        print("Creating extra Node Solution Potential")
        CreateSolution(device, region, "Potential")
    # elec_i = "n_i*exp(Potential/V_t)"
    # hole_i = "n_i^2/IntrinsicElectrons"
    elec_i = "n_i*exp(Potential/V_t)"  #revised
    hole_i = "n_i*exp(-Potential/V_t)" #revised
    #hole_i = "n_i^2/IntrinsicElectrons"#revised
    charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
    pcharge_i = "-ElectronCharge * IntrinsicCharge"
    
    # require NetDoping
    for i in (
        ("IntrinsicElectrons", elec_i),
         ("IntrinsicHoles", hole_i),
         ("IntrinsicCharge", charge_i),
         ("PotentialIntrinsicCharge", pcharge_i)
    ):
        n = i[0]
        e = i[1]
        CreateNodeModel(device, region, n, e)
        CreateNodeModelDerivative(device, region, n, e, "Potential")

    ### TODO: Edge Average Model
    for i in (
        ("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
      ("PotentialEdgeFlux", "Permittivity * ElectricField")
    ):
        n = i[0]
        e = i[1]
        CreateEdgeModel(device, region, n, e)
        CreateEdgeModelDerivatives(device, region, n, e, "Potential")

    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreateSiCPotentialOnlyContact(device, region, contact, is_circuit=False):
    '''
      Creates the potential equation at the contact
      if is_circuit is true, than use node given by GetContactBiasName
    '''
    # Means of determining contact charge
    # Same for all contacts
    # if not InNodeModelList(device, region, "contactcharge_node"):
    #     CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
    # #### TODO: This is the same as D-Field
    # if not InEdgeModelList(device, region, "contactcharge_edge"):
    #     CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
    #     CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")
    # contact_model_name = GetContactNodeModelName(contact) 
    # if is_circuit:
    #     contact_equation(device=device, contact=contact, name="PotentialEquation",
    #                      node_model=contact_model_name, edge_model="",
    #                      node_charge_model="", edge_charge_model="contactcharge_edge",
    #                      #node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
    #                      node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
    # else:
    #     contact_equation(device=device, contact=contact, name="PotentialEquation",
    #                      node_model=contact_model_name, edge_model="",
    #                      node_charge_model="", edge_charge_model="contactcharge_edge",
    #                      #node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
    #                      node_current_model="", edge_current_model="")
    # set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)
    contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
    -V_t*log({1}/n_i), \
    V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)
    contact_model_name = GetContactNodeModelName(contact)  #挪到前面
    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
    # Simplify it too complicated
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
    if is_circuit:
        CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

    if is_circuit:
        contact_equation(device=device, contact=contact, name="PotentialEquation",
                          node_model=contact_model_name, edge_model="",
                          node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
                          node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
    else:
        contact_equation(device=device, contact=contact, name="PotentialEquation",
                          node_model=contact_model_name, edge_model="",
                          node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
                          node_current_model="", edge_current_model="")

def CreateSRH(device, region):
    USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
    Gn = "-ElectronCharge * USRH"
    Gp = "+ElectronCharge * USRH"
    CreateNodeModel(device, region, "USRH", USRH)
    CreateNodeModel(device, region, "ElectronGeneration", Gn)
    CreateNodeModel(device, region, "HoleGeneration", Gp)
    for i in ("Electrons", "Holes"):
        CreateNodeModelDerivative(device, region, "USRH", USRH, i)
        CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
        CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)
        print("125df")

def CreateECE(device, region, mu_n):
    CreateElectronCurrent(device, region, mu_n)

    NCharge = "-ElectronCharge * Electrons"
    CreateNodeModel(device, region, "NCharge", NCharge)
    CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

    equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
             time_node_model = "NCharge",
             edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region, mu_p):
    CreateHoleCurrent(device, region, mu_p)
    PCharge = "ElectronCharge * Holes"
    CreateNodeModel(device, region, "PCharge", PCharge)
    CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

    equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
             time_node_model = "PCharge",
             edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
    pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
    CreateNodeModel(device, region, "PotentialNodeCharge", pne)
    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

    # equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
    #          node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
    #          time_node_model="", variable_update="log_damp")

    evsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
             time_node_model="", variable_update="default")  #revised
             #time_node_model="", variable_update="log_damp")
def CreateSiCDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p"):
    print("xxd")
    CreatePE(device, region)
    CreateBernoulli(device, region)
    CreateSRH(device, region)
    CreateECE(device, region, mu_n)
    CreateHCE(device, region, mu_p)


def CreateSiCDriftDiffusionAtContact(device, region, contact, is_circuit=False): 
    '''
      Restrict electrons and holes to their equilibrium values
      Integrates current into circuit
    '''
    contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
    contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
    contact_electrons_name = "{0}nodeelectrons".format(contact)
    contact_holes_name = "{0}nodeholes".format(contact)

    CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
    #TODO: The simplification of the ifelse statement is time consuming
#  CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

    CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

    #TODO: keyword args
    if is_circuit:
        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
                         node_model=contact_electrons_name,
                         edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
                         node_model=contact_holes_name,
                         edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

    else:
        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
                         node_model=contact_electrons_name,
                         edge_current_model="ElectronCurrent")

        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
                         node_model=contact_holes_name,
                         edge_current_model="HoleCurrent")


def CreateOxidePotentialOnly(device, region, update_type="default"):
    '''
      Create electric field model in oxide
      Creates Potential solution variable if not available
    '''
    if not InNodeModelList(device, region, "Potential"):
        print("Creating Node Solution Potential")
        CreateSolution(device, region, "Potential")

    efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
    # this needs to remove derivatives w.r.t. independents
    CreateEdgeModel(device, region, "ElectricField", efield)
    CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
    dfield="Permittivity*ElectricField"
    CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
    CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             edge_model="PotentialEdgeFlux", variable_update=update_type)


def CreateSiCOxideInterface(device, interface):
    '''
      continuous potential at interface
    '''
    model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
    interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiCSiCInterface(device, interface):
    '''
      Enforces potential, electron, and hole continuity across the interface
    '''
    CreateSiCOxideInterface(device, interface)
    ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
    interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
    hname = CreateContinuousInterfaceModel(device, interface, "Holes")
    interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")

other scripts like mode create.py is not changed

Do you know why it cannot be solved as below:

Could not find equation ElectronContinuityEquation on region0 Ndrift for interface equation ElectronContinuityEquation
10

here simple physics for Poly gate:

# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
# celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
# chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
celec_model = "(0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"   #revised
chole_model = "(0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"  #revised

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
T      = 300 # K
# eps_si = 11.1
# eps_ox = 3.9
# eps_sic = 9.7
eps_Poly = 11.1
# eps = {"SiC": 9.7,"Poly": 11.1,"Gox":3.9}
# TODO: make this temperature dependent
n_i    = 1.0e10 # #/cm^3  # for silicon
# n_i    = 1e-8 # #/cm^3  #for SiC
# n_i = {"SiC":1e-8,"Poly":1e10}
# constant in our approximation
mu_n   = 400 #for Si
mu_p   = 200 #for Si
# mu_n   = 950 # for SiC
# mu_p   = 125 # for SiC
# mu_n = {"SiC":950,"Poly":400}  #注意浓度、电场的调制等效应
# mu_p = {"SiC":125,"Poly":200}
# taun = {"SiC":2.5e-6,"Poly":2.5e-6}
# taup = {"SiC":0.5e-6,"Poly":0.5e-6}
taun = 2.5e-6 # for Poly
taup = 0.5e-6 # for Poly


def GetContactBiasName(contact):
    return "{0}_bias".format(contact)

def GetContactNodeModelName(contact):
    return "{0}nodemodel".format(contact)

def PrintCurrents(device, contact):
    '''
       print out contact currents
    '''
    # TODO add charge
    contact_bias_name = GetContactBiasName(contact)
    electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
    hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
    total_current   = electron_current + hole_current                                        
    voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
    print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
    conteq="Permittivity*ElectricField"
    contact_bias_name  = GetContactBiasName(contact)
    contact_model_name = GetContactNodeModelName(contact)
    eq = "Potential - {0}".format(contact_bias_name)
    CreateContactNodeModel(device, contact, contact_model_name, eq)
    CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

    #TODO: make everyone use dfield
    if not InEdgeModelList(device, region, contactcharge_edge):
        CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
        CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

    contact_equation(device=device, contact=contact, name="PotentialEquation",
                     node_model=contact_model_name, edge_charge_model= contactcharge_edge)




#####
##### Constants
#####
def SetOxideParameters(device, region, T):
    '''
      Sets physical parameters
    '''
    set_parameter(device=device, region=region, name="Permittivity",   value=eps_si * eps_0)
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)

def SetPolyParameters(device, region, T):
    '''
      Sets physical parameters assuming constants
    '''
    #### TODO: make T a free parameter and T dependent parameters as models
    set_parameter(device=device, region=region, name="Permittivity",   value=eps_si * eps_0) #revise for sic
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)
    set_parameter(device=device, region=region, name="n_i",            value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="T",              value=T)
    set_parameter(device=device, region=region, name="kT",             value=k * T)
    set_parameter(device=device, region=region, name="V_t",            value=k*T/q)
    set_parameter(device=device, region=region, name="mu_n",           value=mu_n) #revise for sic
    set_parameter(device=device, region=region, name="mu_p",           value=mu_p) #revise for sic
    #default SRH parameters
    set_parameter(device=device, region=region, name="n1", value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="p1", value=n_i) #revise for sic
    set_parameter(device=device, region=region, name="taun", value=taun) #revise for sic
    set_parameter(device=device, region=region, name="taup", value=taup) #revise for sic


def CreatePolyPotentialOnly(device, region):
    '''
      Creates the physical models for a Silicon region
    '''
    if not InNodeModelList(device, region, "Potential"):
        print("Creating extra Node Solution Potential")
        CreateSolution(device, region, "Potential")
    # elec_i = "n_i*exp(Potential/V_t)"
    # hole_i = "n_i^2/IntrinsicElectrons"
    elec_i = "n_i*exp(Potential/V_t)"  #revised
    hole_i = "n_i*exp(-Potential/V_t)" #revised
    #hole_i = "n_i^2/IntrinsicElectrons"#revised
    charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
    pcharge_i = "-ElectronCharge * IntrinsicCharge"
    

    
    # require NetDoping
    for i in (
        ("IntrinsicElectrons", elec_i),
         ("IntrinsicHoles", hole_i),
         ("IntrinsicCharge", charge_i),
         ("PotentialIntrinsicCharge", pcharge_i)
    ):
        n = i[0]
        e = i[1]
        CreateNodeModel(device, region, n, e)
        CreateNodeModelDerivative(device, region, n, e, "Potential")



    # log(((5.000000000000000e-01 * abs((-NetDoping + pow((pow(NetDoping,2) + (4 * pow(n_i,2))),5.000000000000000e-01)))) * pow(n_i,(-1))))


    ### TODO: Edge Average Model
    for i in (
        ("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
      ("PotentialEdgeFlux", "Permittivity * ElectricField")
    ):
        n = i[0]
        e = i[1]
        CreateEdgeModel(device, region, n, e)
        CreateEdgeModelDerivatives(device, region, n, e, "Potential")

    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreatePolyPotentialOnly(device, region):
    '''
      Creates the physical models for a Silicon region
    '''
    if not InNodeModelList(device, region, "Potential"):
        print("Creating extra Node Solution Potential")
        CreateSolution(device, region, "Potential")
    # elec_i = "n_i*exp(Potential/V_t)"
    # hole_i = "n_i^2/IntrinsicElectrons"
    elec_i = "%s*exp(Potential/V_t)"  %n_i["SiC"] #revised
    hole_i = "%s*exp(-Potential/V_t)" %n_i["SiC"] #revised
    #hole_i = "n_i^2/IntrinsicElectrons"#revised
    charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
    pcharge_i = "-ElectronCharge * IntrinsicCharge"
    
    # require NetDoping
    for i in (
        ("IntrinsicElectrons", elec_i),
         ("IntrinsicHoles", hole_i),
         ("IntrinsicCharge", charge_i),
         ("PotentialIntrinsicCharge", pcharge_i)
    ):
        n = i[0]
        e = i[1]
        CreateNodeModel(device, region, n, e)
        CreateNodeModelDerivative(device, region, n, e, "Potential")

    ### TODO: Edge Average Model
    for i in (
        ("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
      ("PotentialEdgeFlux", "Permittivity * ElectricField")
    ):
        n = i[0]
        e = i[1]
        CreateEdgeModel(device, region, n, e)
        CreateEdgeModelDerivatives(device, region, n, e, "Potential")

    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreatePolyPotentialOnlyContact(device, region, contact, is_circuit=False):
    '''
      Creates the potential equation at the contact
      if is_circuit is true, than use node given by GetContactBiasName
    '''
    # Means of determining contact charge
    # Same for all contacts
    # if not InNodeModelList(device, region, "contactcharge_node"):
    #     CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
    # #### TODO: This is the same as D-Field
    # if not InEdgeModelList(device, region, "contactcharge_edge"):
    #     CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
    #     CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")
    # contact_model_name = GetContactNodeModelName(contact) 
    # if is_circuit:
    #     contact_equation(device=device, contact=contact, name="PotentialEquation",
    #                      node_model=contact_model_name, edge_model="",
    #                      node_charge_model="", edge_charge_model="contactcharge_edge",
    #                      #node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
    #                      node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
    # else:
    #     contact_equation(device=device, contact=contact, name="PotentialEquation",
    #                      node_model=contact_model_name, edge_model="",
    #                      node_charge_model="", edge_charge_model="contactcharge_edge",
    #                      #node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
    #                      node_current_model="", edge_current_model="")
    # set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)
    contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
    -V_t*log({1}/n_i), \
    V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)
    contact_model_name = GetContactNodeModelName(contact)  #挪到前面
    CreateContactNodeModel(device, contact, contact_model_name, contact_model)
    # Simplify it too complicated
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
    if is_circuit:
        CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

    if is_circuit:
        contact_equation(device=device, contact=contact, name="PotentialEquation",
                         node_model=contact_model_name, edge_model="",
                         node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
                         node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
    else:
        contact_equation(device=device, contact=contact, name="PotentialEquation",
                         node_model=contact_model_name, edge_model="",
                         node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
                         node_current_model="", edge_current_model="")

def CreateSRH(device, region):
    USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
    Gn = "-ElectronCharge * USRH"
    Gp = "+ElectronCharge * USRH"
    CreateNodeModel(device, region, "USRH", USRH)
    CreateNodeModel(device, region, "ElectronGeneration", Gn)
    CreateNodeModel(device, region, "HoleGeneration", Gp)
    for i in ("Electrons", "Holes"):
        CreateNodeModelDerivative(device, region, "USRH", USRH, i)
        CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
        CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)
        print("125df")

def CreateECE(device, region, mu_n):
    CreateElectronCurrent(device, region, mu_n)

    NCharge = "-ElectronCharge * Electrons"
    CreateNodeModel(device, region, "NCharge", NCharge)
    CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

    equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
             time_node_model = "NCharge",
             edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region, mu_p):
    CreateHoleCurrent(device, region, mu_p)
    PCharge = "ElectronCharge * Holes"
    CreateNodeModel(device, region, "PCharge", PCharge)
    CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

    equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
             time_node_model = "PCharge",
             edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
    pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
    CreateNodeModel(device, region, "PotentialNodeCharge", pne)
    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
    CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

    # equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
    #          node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
    #          time_node_model="", variable_update="log_damp")

    evsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
             time_node_model="", variable_update="default")  #revised
             #time_node_model="", variable_update="log_damp")
def CreatePolyDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p"):
    CreatePE(device, region)
    CreateBernoulli(device, region)
    CreateSRH(device, region)
    CreateECE(device, region, mu_n)
    CreateHCE(device, region, mu_p)


def CreatePolyDriftDiffusionAtContact(device, region, contact, is_circuit=False): 
    '''
      Restrict electrons and holes to their equilibrium values
      Integrates current into circuit
    '''
    contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
    contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
    contact_electrons_name = "{0}nodeelectrons".format(contact)
    contact_holes_name = "{0}nodeholes".format(contact)

    CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
    #TODO: The simplification of the ifelse statement is time consuming
#  CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

    CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
    CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

    #TODO: keyword args
    if is_circuit:
        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
                         node_model=contact_electrons_name,
                         edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
                         node_model=contact_holes_name,
                         edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

    else:
        contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
                         node_model=contact_electrons_name,
                         edge_current_model="ElectronCurrent")

        contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
                         node_model=contact_holes_name,
                         edge_current_model="HoleCurrent")


def CreateOxidePotentialOnly(device, region, update_type="default"):
    '''
      Create electric field model in oxide
      Creates Potential solution variable if not available
    '''
    if not InNodeModelList(device, region, "Potential"):
        print("Creating Node Solution Potential")
        CreateSolution(device, region, "Potential")

    efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
    # this needs to remove derivatives w.r.t. independents
    CreateEdgeModel(device, region, "ElectricField", efield)
    CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
    dfield="Permittivity*ElectricField"
    CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
    CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
    equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
             edge_model="PotentialEdgeFlux", variable_update=update_type)


def CreatePolyOxideInterface(device, interface):
    '''
      continuous potential at interface
    '''
    model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
    interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreatePolyPolyInterface(device, interface):
    '''
      Enforces potential, electron, and hole continuity across the interface
    '''
    CreatePolyOxideInterface(device, interface)
    ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
    interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
    hname = CreateContinuousInterfaceModel(device, interface, "Holes")
    interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")
11

Hi @caihng

I am unable to run your scripts at the moment.

The error regarding ElectronContinuityEquation for interface equation means that you should not add ElectronContinityEquation with the interface_equation command, unless the equation command has also been created for the same equation on both regions of the interface.

12

Invalid or Divide-by-zero floating point exceptions are caused if your simulation evaluates a function that results in these problems. For your case, it is possible that either the NetDoping or n_i is 0.0.

In an ifelse expression, you must prevent floating point errors from occurring in the if model expression or the else model expression. In your testing, you can create additional test models, and use devsim.get_node_model_values results in an FPE, without needing to do a solve.

13

hello @Juan ，

I want to plot IV and CV curve of a mos or diode and i looked into the example of ssac.diode

i don’t understand why the capacitor is defined as ssac_imag

cap=get_circuit_node_value(node="V1.I", solution="ssac_imag")/ (-2*math.pi)
    # cur=get_contact_current(device=device, contact="top", equation="PotentialEquation")
    print("capacitance {0} {1}".format(v, cap))

Also i want to know what is the unit of frequency here.

    solve(type="ac", frequency=1)

if frequency is 1MHz, do i need to write it as 1e6 ？

14

Hi @caihng

You need to look at the small signal circuit model of the device you are characterizing. A simple diode model has an admittance of:

y = g + j 2 \pi f c

and using

i_{src} + y \cdot v_{src} = 0

where v_{src}, is the applied small-signal voltage, and i_{src} is the resulting current, you can derive that

c = -\frac{\mathrm{Im}\left(i_{src}\right)}{2 \pi f}

The frequency, f is in Hz. A low frequency will reveal the depletion capacitance. A high frequency will be a combination of the diffusion and depletion capacitance.

A similar small-signal circuit can be derived for a MOS gate capacitance.

15

@Juan

cap=get_circuit_node_value(node="V1.I", solution="ssac_imag")/ (-2*math.pi)

But why frequency is not mentioned here ?

Do we need to modify it to be

cap=get_circuit_node_value(node="V1.I", solution="ssac_imag")/ (-2*f*math.pi)

I tried to put “f” in the equation, but the C didn’t change…

16

Hi @caihng

The simulation frequency is specified here:

    solve(type="ac", frequency=1)

in Hz. You should put the same value in the capacitance equation as you put on the “ac” solve command.

17

@Juan

yes, i changed the f value in the solve eqution, but when i change the cap eqution at the same time，the result of cap does not change with f which seems not right ? i think different f should mean different C

f=1e6
solve(type="ac", frequency=f)
cap=get_circuit_node_value(node="V1.I", solution="ssac_imag")/ (-2*f*math.pi)

Do you have any other examples to solve the ac signal for a diode or capacitor ?

18

The capacitance will be dependent on bias, recombination lifetime. To see the effect in some devices, you may have to go to frequencies in the GHz or THz region.

Please see this project where frequency dependent results are shown:
https://github.com/devsim/devsim_bjt_example
Also see the pdf files in the doc directory.