The visualization of the MOS_2d example

Hello, everyone! I am trying to learn the example of MOS2. But I can’t acquire a potential or carrier density distribution in the VISIT output file. I can only see the mesh in the VISIT. Below is the case file I used. Could anyone help me? Thank you.


# 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

import mos_2d_create
from mos_2d_physics import *

import devsim
oxide_regions   = ("oxide",)
silicon_regions = ("gate", "bulk")
all_regions     = ("gate", "bulk", "oxide")

for i in all_regions:
    createSolution(device, i, "Potential")

for i in silicon_regions:
    setSiliconParameters(device, i)
    createSiliconPotentialOnly(device, i)
for i in oxide_regions:
    setOxideParameters(device, i)
    createOxidePotentialOnly(device, "oxide")

createSiliconPotentialOnlyContact(device, "gate", "gate")
createSiliconPotentialOnlyContact(device, "bulk", "drain")
createSiliconPotentialOnlyContact(device, "bulk", "source")
createSiliconPotentialOnlyContact(device, "bulk", "body")

createSiliconOxideInterface(device, "bulk_oxide")
createSiliconOxideInterface(device, "gate_oxide")

devsim.solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
devsim.solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)

createSolution(device, "gate", "Electrons")
createSolution(device, "gate", "Holes")
devsim.set_node_values(device=device, region="gate", name="Electrons", init_from="IntrinsicElectrons")
devsim.set_node_values(device=device, region="gate", name="Holes",     init_from="IntrinsicHoles")
createSiliconDriftDiffusion(device, "gate")
createSiliconDriftDiffusionAtContact(device, "gate", "gate")

createSolution(device, "bulk", "Electrons")
createSolution(device, "bulk", "Holes")
devsim.set_node_values(device=device, region="bulk", name="Electrons", init_from="IntrinsicElectrons")
devsim.set_node_values(device=device, region="bulk", name="Holes",     init_from="IntrinsicHoles")
createSiliconDriftDiffusion(device, "bulk")

createSiliconDriftDiffusionAtContact(device, "bulk", "drain")
createSiliconDriftDiffusionAtContact(device, "bulk", "source")
createSiliconDriftDiffusionAtContact(device, "bulk", "body")

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

devsim.element_from_edge_model(edge_model="ElectricField", device=device, region="bulk")

devsim.write_devices(file="mos_2d_dd.msh", type="devsim")

with open("", "w", encoding="utf-8") as ofh:
    ofh.write('import devsim\n')
    for p in devsim.get_parameter_list():
        ofh.write('devsim.set_parameter(name="%s", value=%s)\n' % (p, v))
    for i in devsim.get_device_list():
        for p in devsim.get_parameter_list(device=i):
            v=repr(devsim.get_parameter(device=i, name=p))
            ofh.write('devsim.set_parameter(device="%s", name="%s", value=%s)\n' % (i, p, v))

    for i in devsim.get_device_list():
        for j in devsim.get_region_list(device=i):
            for p in devsim.get_parameter_list(device=i, region=j):
                v=repr(devsim.get_parameter(device=i, region=j, name=p))
                ofh.write('devsim.set_parameter(device="%s", region="%s", name="%s", value=%s)\n' % (i, j, p, v))

devsim.write_devices( file="mos_2dd", type="vtk")


# 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

import devsim

def printCurrents(device, contact, bias):
    ecurr=devsim.get_contact_current(contact=contact, equation="ElectronContinuityEquation", device=device)
    hcurr=devsim.get_contact_current(contact=contact, equation="HoleContinuityEquation", device=device)
    tcurr = ecurr + hcurr
    print("%s %g %g %g %g" % (contact, bias, ecurr, hcurr, tcurr))

#### Constants
def setOxideParameters(device, region):
    devsim.set_parameter(device=device, region=region, name="Permittivity",   value=3.9*eps)
    devsim.set_parameter(device=device, region=region, name="ElectronCharge", value=q)

def setSiliconParameters(device, region):
    for name, value in (
        ("Permittivity", 11.1*eps),
      ("ElectronCharge", q),
      ("n_i", 1.0e10),
      ("kT", eps * T),
      ("V_t", k*T/q),
      ("mu_n", 400),
      ("mu_p", 200),
        devsim.set_parameter(device=device, region=region, name=name, value=value)

def createSolution(device, region, name):
    devsim.node_solution(device=device, region=region, name=name)
    devsim.edge_from_node_model(device=device, region=region, node_model=name)

def createSiliconPotentialOnly(device, region):
    ie = devsim.node_model(device=device, region=region, name="IntrinsicElectrons", equation="n_i*exp(Potential/V_t)")
    res = devsim.node_model(device=device, region=region, name="IntrinsicElectrons:Potential", equation="diff(%s, Potential)" % ie)
    for name, equation in (
        ("IntrinsicHoles",                         "n_i^2/IntrinsicElectrons"),
      ("IntrinsicHoles:Potential",               "diff(n_i^2/IntrinsicElectrons, Potential)"),
      ("IntrinsicCharge",                        "IntrinsicHoles-IntrinsicElectrons + NetDoping"),
      ("IntrinsicCharge:Potential",              "diff(IntrinsicHoles-IntrinsicElectrons, Potential)"),
      ("PotentialIntrinsicNodeCharge",           "-ElectronCharge*IntrinsicCharge"),
      ("PotentialIntrinsicNodeCharge:Potential", "diff(-ElectronCharge*IntrinsicCharge, Potential)"),
        devsim.node_model(device=device, region=region, name=name, equation=equation)

    for name, equation in (
        ("ElectricField",                  "(Potential@n0-Potential@n1)*EdgeInverseLength"),
      ("ElectricField:Potential@n0",     "EdgeInverseLength"),
      ("ElectricField:Potential@n1",     "-ElectricField:Potential@n0"),
      ("PotentialEdgeFlux",              "Permittivity*ElectricField"),
      ("PotentialEdgeFlux:Potential@n0", "diff(Permittivity*ElectricField,Potential@n0)"),
      ("PotentialEdgeFlux:Potential@n1", "-PotentialEdgeFlux:Potential@n0"),
        devsim.edge_model(device=device, region=region, name=name, equation=equation)

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

def createSiliconPotentialOnlyContact(device, region, contact):
    bias_name="%sbias" % contact
    format_dict = { "contact" : contact}
    devsim.set_parameter(device=device, region=region, name=bias_name, value=0.0)
    for name, equation in (
        ("celec_%(contact)s", "1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5))"),
      ("chole_%(contact)s", "1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5))"),
      ("%(contact)snodemodel", '''
        ifelse(NetDoping > 0,
      ("%(contact)snodemodel:Potential", "1"),
        name_sub = name % format_dict
        equation_sub = equation % format_dict
        devsim.contact_node_model(device=device, contact=contact, name=name_sub, equation=equation_sub)

    devsim.contact_equation(device=device, contact=contact, name="PotentialEquation",
                            node_model="%snodemodel" % contact)

def createSiliconDriftDiffusion(device, region):
    for name, equation in (
        ("PotentialNodeCharge",           "-ElectronCharge*(Holes -Electrons + NetDoping)"),
      ("PotentialNodeCharge:Electrons", "+ElectronCharge"),
      ("PotentialNodeCharge:Holes",     "-ElectronCharge"),
        devsim.node_model(device=device, region=region, name=name, equation=equation)

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

    createBernoulli(device, region)
    createElectronCurrent(device, region)
    createHoleCurrent(device, region)

    NCharge="-ElectronCharge * Electrons"
    devsim.node_model(device=device, region=region, name="NCharge", equation=NCharge)
    devsim.node_model(device=device, region=region, name="NCharge:Electrons", equation=dNChargedn)

    PCharge="-ElectronCharge * Holes"
    devsim.node_model(device=device, region=region, name="PCharge", equation=PCharge)
    devsim.node_model(device=device, region=region, name="PCharge:Holes", equation=dPChargedp)

    ni=devsim.get_parameter(device=device, region=region, name="n_i")
    devsim.set_parameter(device=device, region=region, name="n1", value=ni)
    devsim.set_parameter(device=device, region=region, name="p1", value=ni)
    devsim.set_parameter(device=device, region=region, name="taun", value=1e-5)
    devsim.set_parameter(device=device, region=region, name="taup", value=1e-5)

    USRH="-ElectronCharge*(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
    dUSRHdn="simplify(diff(%s, Electrons))" % USRH
    dUSRHdp="simplify(diff(%s, Holes))" % USRH
    devsim.node_model(device=device, region=region , name="USRH", equation=USRH)
    devsim.node_model(device=device, region=region , name="USRH:Electrons", equation=dUSRHdn)
    devsim.node_model(device=device, region=region , name="USRH:Holes", equation=dUSRHdp)

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

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

def createSiliconDriftDiffusionAtContact(device, region, contact):
    format_dict = { "contact" : contact }
    for name, equation in (
        ("%(contact)snodeelectrons",           "ifelse(NetDoping > 0, Electrons - celec_%(contact)s, Electrons - n_i^2/chole_%(contact)s)"),
      ("%(contact)snodeholes",               "ifelse(NetDoping < 0, Holes - chole_%(contact)s, Holes - n_i^2/celec_%(contact)s)"),
      ("%(contact)snodeelectrons:Electrons", "1.0"),
      ("%(contact)snodeholes:Holes",         "1.0"),
        name_sub = name % format_dict
        equation_sub = equation % format_dict
        devsim.contact_node_model(device=device, contact=contact, name=name_sub, equation=equation_sub)

    devsim.contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
                            node_model="%snodeelectrons" % contact, edge_current_model="ElectronCurrent")

    devsim.contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
                            node_model="%snodeholes" % contact, edge_current_model="HoleCurrent")

def createOxidePotentialOnly(device, region):
    for name, equation in (
        ("ElectricField",                   "(Potential@n0 - Potential@n1)*EdgeInverseLength"),
      ("ElectricField:Potential@n0",      "EdgeInverseLength"),
      ("ElectricField:Potential@n1",      "-EdgeInverseLength"),
      ("PotentialEdgeFlux",               "Permittivity*ElectricField"),
      ("PotentialEdgeFlux:Potential@n0",  "diff(Permittivity*ElectricField, Potential@n0)"),
      ("PotentialEdgeFlux:Potential@n1",  "-PotentialEdgeFlux:Potential@n0"),
        devsim.edge_model(device=device, region=region, name=name, equation=equation)

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

def createSiliconOxideInterface(device, interface):
    for name, equation in (
        ("continuousPotential",              "Potential@r0-Potential@r1"),
      ("continuousPotential:Potential@r0", "1"),
      ("continuousPotential:Potential@r1", "-1"),
        devsim.interface_model(device=device, interface=interface, name=name, equation=equation)

    devsim.interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model="continuousPotential", type="continuous")

def createSiliconSiliconInterface(device, interface):
    for variable in ("Potential", "Electrons", "Holes"):
        format_dict = { "var", variable }
        for name, equation in (
            ("continuous%(var)s", "%(var)s@r0-%(var)s@r1"),
          ("continuous%(var)s:%(var)s@r0", "1"),
          ("continuous%(var)s:%(var)s@r1", "-1"),
            name_sub = name % format_dict
            equation_sub = equation % format_dict
            devsim.interface_model(device=device, interface=interface, name=name_sub, equation=equation_sub)
        eqname = "%sEquation" % variable
        ieqname = "continuous%s" % variable
        devsim.interface_equation(device=device, interface=interface, name=eqname,
                                  interface_model=ieqname, type="continuous")

def createBernoulli(device, region):
    #### test for requisite models here
    vdiffstr="(Potential@n0 - Potential@n1)/V_t"
    for name, equation in (
        ("vdiff", vdiffstr),
      ("vdiff:Potential@n0",  "V_t^(-1)"),
      ("vdiff:Potential@n1",  "-vdiff:Potential@n0"),
      ("Bern01",              "B(vdiff)"),
      ("Bern01:Potential@n0", "dBdx(vdiff) * vdiff:Potential@n0"),
      ("Bern01:Potential@n1", "-Bern01:Potential@n0"),
      ("Bern10",              "Bern01 + vdiff"),
      ("Bern10:Potential@n0", "Bern01:Potential@n0 + vdiff:Potential@n0"),
      ("Bern10:Potential@n1", "Bern01:Potential@n1 + vdiff:Potential@n1"),
        devsim.edge_model(device=device, region=region, name=name, equation=equation)

def createElectronCurrent(device, region):
    Jn="ElectronCharge*mu_n*EdgeInverseLength*V_t*(Electrons@n1*Bern10 - Electrons@n0*Bern01)"
    devsim.edge_model(device=device, region=region, name="ElectronCurrent", equation=Jn)
    for variable in ("Electrons", "Potential"):
        der = "simplify(diff(%s, %s))" % (Jn, variable)
        devsim.edge_model(device=device, region=region, name="ElectronCurrent", equation=der)

def createHoleCurrent (device, region):
    Jp="-ElectronCharge*mu_p*EdgeInverseLength*V_t*(Holes@n1*Bern01 - Holes@n0*Bern10)"
    devsim.edge_model(device=device, region=region, name="HoleCurrent", equation=Jp)
    for variable in ("Holes", "Potential"):
        der = "simplify(diff(%s, %s))" % (Jp, variable)
        devsim.edge_model(device=device, region=region, name="HoleCurrent", equation=der)


# 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

import devsim

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

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

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

bulk_doping=        -1e15
body_doping=        -1e19
drain_doping=       1e19
source_doping=      1e19
gate_doping=        1e20

y_channel_spacing=    1e-8
y_diffusion_spacing=  1e-6
y_gate_top_spacing=   1e-8
y_gate_mid_spacing=   (gate_thickness * 0.25)
y_gate_bot_spacing=   1e-8
y_oxide_mid_spacing=  (oxide_thickness * 0.25)
x_channel_spacing=  1e-6
max_y_spacing=      1e-4
max_x_spacing=      1e-2
y_bulk_mid_spacing=(device_thickness * 0.25)

x_bulk_left=  0.0
x_bulk_right= (x_bulk_left + device_width)
x_center=     (0.5 * (x_bulk_left + x_bulk_right))
x_gate_left=  (x_center - 0.5 * (gate_width))
x_gate_right= (x_center + 0.5 * (gate_width))
x_device_left=(x_bulk_left - air_thickness)
x_device_right=(x_bulk_right + air_thickness)

y_bulk_top=     0.0
y_oxide_top=    (y_bulk_top - oxide_thickness)
y_oxide_mid=    (0.5 * (y_oxide_top + y_bulk_top))
y_gate_top=     (y_oxide_top - gate_thickness)
y_gate_mid=     (0.5 * (y_gate_top + y_oxide_top))
y_device_top=   (y_gate_top - air_thickness)
y_bulk_bottom=  (y_bulk_top + device_thickness)
y_bulk_mid=     (0.5 * (y_bulk_top + y_bulk_bottom))
y_device_bottom=(y_bulk_bottom + air_thickness)
y_diffusion=    (y_bulk_top + diffusion_thickness)

devsim.create_2d_mesh(  mesh="mos")
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_device_top, ps=max_y_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_gate_top, ps=y_gate_top_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_gate_mid, ps=y_gate_mid_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_oxide_top, ps=y_gate_bot_spacing, ns=y_oxide_mid_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_oxide_mid, ps=y_oxide_mid_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_bulk_top, ns=y_oxide_mid_spacing, ps=y_channel_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_diffusion, ps=y_diffusion_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_bulk_mid, ps=y_bulk_mid_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_bulk_bottom, ns=y_bulk_bottom_spacing, ps=max_y_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="y", pos=y_device_bottom, ps=max_y_spacing)

devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_center, ps=x_channel_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_gate_left, ps=x_channel_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_bulk_left, ps=x_diffusion_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_gate_right, ps=x_channel_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_bulk_right, ps=x_diffusion_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_device_right, ps=max_x_spacing)
devsim.add_2d_mesh_line(mesh="mos", dir="x", pos=x_device_left, ps=max_x_spacing)

devsim.add_2d_region(mesh="mos", material="Air"    , region="air")
devsim.add_2d_region(mesh="mos", material="Silicon", region="bulk", xl=x_bulk_left, xh=x_bulk_right, yl=y_bulk_bottom, yh=y_bulk_top)
devsim.add_2d_region(mesh="mos", material="Silicon", region="gate", xl=x_gate_left, xh=x_gate_right, yl=y_oxide_top, yh=y_gate_top)
devsim.add_2d_region(mesh="mos", material="Oxide"  , region="oxide", xl=x_gate_left, xh=x_gate_right, yl=y_bulk_top, yh=y_oxide_top)

devsim.add_2d_contact(mesh="mos", name="gate", region="gate", yl=y_gate_top, yh=y_gate_top, material="metal")
devsim.add_2d_contact(mesh="mos", name="body", region="bulk", yl=y_bulk_bottom, yh=y_bulk_bottom, material="metal")
devsim.add_2d_contact(mesh="mos", name="source", region="bulk", yl=y_bulk_top, yh=y_bulk_top, xl=x_device_left, xh=x_gate_left, material="metal")
devsim.add_2d_contact(mesh="mos", name="drain" , region="bulk", yl=y_bulk_top, yh=y_bulk_top, xl=x_gate_right, xh=x_device_right, material="metal")

devsim.add_2d_interface(mesh="mos", name="gate_oxide", region0="gate", region1="oxide")
devsim.add_2d_interface(mesh="mos", name="bulk_oxide", region0="bulk", region1="oxide")
devsim.create_device(mesh="mos", device=device)

format_dict= {
    'gate_doping' : gate_doping,
  'source_doping' : source_doping,
  'drain_doping' : drain_doping,
  'body_doping' : body_doping,
  'bulk_doping' : bulk_doping,
  'x_gate_left' : x_gate_left,
  'x_gate_right' : x_gate_right,
  'x_diffusion_decay' : x_diffusion_decay,
  'y_diffusion' : y_diffusion,
  'y_diffusion_decay' : y_diffusion_decay,
  'y_bulk_bottom' : y_bulk_bottom,

devsim.node_model(name="NetDoping"   , device=device, region="gate", equation="%(gate_doping)s" % format_dict)
devsim.node_model(name="DrainDoping" , device=device, region="bulk", equation="0.25*%(drain_doping)s*erfc((x-%(x_gate_left)s)/%(x_diffusion_decay)s)*erfc((y-%(y_diffusion)s)/%(y_diffusion_decay)s)" % format_dict)
devsim.node_model(name="SourceDoping", device=device, region="bulk", equation="0.25*%(source_doping)s*erfc(-(x-%(x_gate_right)s)/%(x_diffusion_decay)s)*erfc((y-%(y_diffusion)s)/%(y_diffusion_decay)s)" % format_dict)
devsim.node_model(name="BodyDoping", device=device, region="bulk", equation="0.5*%(body_doping)s*erfc(-(y-%(y_bulk_bottom)s)/%(y_diffusion_decay)s)" % format_dict)
devsim.node_model(name="NetDoping"   , device=device, region="bulk", equation="DrainDoping + SourceDoping + %(bulk_doping)s + BodyDoping" % format_dict)

devsim.write_devices( file="mos_2d", type="vtk")
devsim.write_devices( file="mos_2d.flps", type="floops")

Hi @ghost

Please try the tecplot format in write_devices. The vtk format is problematic in that it is a multi-file format, and I may not have done a very good implementation.

Please note that VisIt has at least 2 different Tecplot readers, so please experiment and see which one works best for you.

Hello, Juan! Thanks for you generous advice. I can visualize the results now. But there is a new question I want to disturb you.
The question is that I am confused to the setting of bias. In the example of Mos_2d. The only setting about bias is in the “createSiliconPotentialOnlyContact” as below:

def createSiliconPotentialOnlyContact(device, region, contact):
    bias_name="%sbias" % contact
    format_dict = { "contact" : contact}
    devsim.set_parameter(device=device, region=region, name=bias_name, value=0)

The example seems to set all bias as 0. I try to change it as followed:

def createSiliconPotentialOnlyContact(device, region, contact):
    bias_name="%sbias" % contact
    format_dict = { "contact" : contact}
    # devsim.set_parameter(device=device, region=region, name=bias_name, value=0)
    devsim.set_parameter(device="mymos", region="bulk", name="sourcebias", value=0)
    devsim.set_parameter(device="mymos", region="bulk", name="bodybias", value=0)
    devsim.set_parameter(device="mymos", region="bulk", name="drainbias", value=0.5)
    devsim.set_parameter(device="mymos", region="gate", name="gatebias", value=1)

However, the potential distribution in the result is not as I set it. The potential in the drain is about 1.And potential in the gate is about 1.5. Why does this happen?

I also found that there is a function “ramp” that can be used to change the bias and get the current trend. But it also only gives the value of drain or source? And in the simulation, does each contact need a initial potential value? Could you tell me how to set the bias or give me some information about the setting of bias?

To answer the first part of your question, the potential is offset from the contacts with the built in potential. This comes from tying the fermi level at equilibrium to the contact potential.

I will have to follow up later concerning usage of the ramp command.

Thank you for your reply. I am trying to learn from other user’s topic.

Hi @ghost

You should definitely use the set_parameter command to set the appropriate bias on each contact. A value of 0.0 for each is appropriate for the equilibrium solution.

Please note when using the devsim.python_packages.ramp.rampbias command uses a parameter value that is set on the device, and not on the region. So please omit the region option from the set_parameter command.

devsim.set_parameter(device="mymos", name="sourcebias", value=0)

If you do not do this, the value will remain fixed, since a parameter set on a region will always override the value set on a device.

Using the python help command, you can see what the rampbias command options are.

rampbias(device, contact, end_bias, step_size, min_step, max_iter, rel_error, abs_error, callback)
    Ramps bias with assignable callback function

where the device and contact options are to name the specific contact being ramped. The specified contact then will read the current parameter setting for the '%sbias' % contact:

    start_bias = ds.get_parameter(device=device, name=GetContactBiasName(contact))

starting from the solution at this bias point, it will then ramp the contact bias value up to end_bias using step_size steps. If the simulation fails to converge at a bias, it will automatically reduce the step size to attempt to get convergence. The ramp will stop if the step_size is reduced below min_step due to too many failures.

The max_iter, rel_error, and abs_error are the convergence criteria for a successful simulation.

The callback function is run at the end of each successful step. For example, it can be used to print the currents at each of the contacts.

Hi Juan:
I really appreciate you taking the time to answer my simple question. I have learned Devsim for some weeks. But I still can’t fully understand it. There are some questions I want to disturb you again.
First, in the Devsim, the equations are assembled on nodal and edge models. For example, the Potential equation:

def CreateSiliconPotentialOnly(device, region):
      Creates the physical models for a Silicon region
    if not InNodeModelList(device, region, "Potential"):
        print("Creating Node Solution Potential")
        CreateSolution(device, region, "Potential")
    elec_i = "n_i*exp(Potential/V_t)"
    hole_i = "n_i^2/IntrinsicElectrons"
    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")

In my understanding, the nodel model gives the source term in each nodel and it is one side of the equation in the calculation. And the edge nodel decide the traffic around the node and it is another side of equation. So the equation it describes is:

I don’t know if i am wrong. Because I found some problem in the Electron Continuity Equation. The code in the example as below:

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)

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")

It seems that the equation in this part is:

However, it is different from this book(, page 85 shows:

Although the code use"Gn", but the description should be equation of SRH recombination rate. So is there an extra negative sign in the equation in this code? Or am I misunderstood?

Besides, I want to know why the Electron Continuity Equation need define “time_node_model” ?

I also noticed when calculate the current density Jn&Jp, the example use Bernoulli Function. Why don’t we directly calculate the following formula:

Can’t the equation be easily described if we ignore the effect of temperature in the Devsim?

And could you explain the equation of contact_electrons_model in the CreateSiliconDriftDiffusionAtContact?

def CreateSiliconDriftDiffusionAtContact(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)

where can I find information about it?

Hello @ghost

Regarding: CreateSIliconPotentialOnly, it solves for an equilibrium solution, as an initital guess for future drift diffusion simulation.

n = n_i \exp\left(\frac{\psi}{V_t}\right)\\ p = \frac{{n_i}^2}{n}

Once this equilibrium solution is available, the n and p are solved as independent solution variables. The values for n and p from the equilibrum solution are used

Regarding the time_node_model, the equations you present as 4.2 and 4.3 are steady state equations. Under time varying conditions, terms such as \frac{\partial n}{\partial t}, \frac{\partial p}{\partial t} are needed to account for accumulation of carriers. You may find a reference to this type of equation as “Fick’s Law”.

Regarding the G and R terms.

G = G_n = G_p

only when you have direct tunneling. SRH recombination assumes an intermediate trap level with a time dependent term set to zero.

\frac{\partial n_t}{\partial t} = 0 = G_n - R_n + R_p - G_p

The Bernoulli function is part of the Scharfetter-Gummel algorithm, which is part of most drift-diffusion finite volume simulators. It is known to be among the most numerically stable formulations for simulation. It is easily extended for the case of a gradient in temperature or bandgap.

It is only possible to solve for J_n or J_p only in the absence of generation or recombination and only a monopolar simulation.

For an ohmic contact, a charge-neutral solution is assumed.

\frac{1}{q}\nabla \cdot \nabla \psi = 0 = n - \frac{{n_i}^2}{n} + N_A - N_D


\frac{1}{q}\nabla \cdot \nabla \psi = 0 = \frac{{n_i}^2}{p} - p + N_A - N_D

So the electrons and holes are constrained to those values at the contact nodes.

I highly recommend a book like:
Analysis and Simulation of Semiconductor Devices by Selberherr (1984)

for classical methods in semiconductor device TCAD simulation.

Also see this document with some of my notes:

Hi @Juan
Thanks for your kindness. I will learn your notes and that book in detail

Actually, I am trying to calculate the heat generation in the Mosfet. Considering only Joule heat, the heat source term is this formula:

But the Jn&Jp&E are both edge model. Could you tell me how to get a nodel model from them?

Besides, I can’t understand that not all geometric boundaries have a boundary condition (the contact or interface) when construct the Poisson equation.But when calculating the thermal conductivity equation, it is necessary to set a boundary condition for all boundary of the area. How can I set up the adiabatic or isothermal boundaries?

HI @ghost,

I have not coded lattice temperature equations before. My thought is that eq 4.23 is a term that would go into a continuity equation. The total equation would involve gradient and time derivative of temperature. However, I do not know this for sure.

The simplest boundary conditions would be a contact setting a constant temperature.

T - T_\mathrm{contact} = 0

Alternatively, a thermal circuit could be added with a resistor and capacitor to model heating and cooling. Looking up the meaning of adiabatic, I think this is handled naturally on boundaries without a contact or interface. This is since no heat would be flowing out of the device.

So I am thinking of isothermal is perhaps a Dirichlet boundary condition, and adiabatic is a Von Neumann boundary condition.

Hi @Juan
Thanks for your advice. By reading the book you suggested, I understand the setting of nature boundary and artificial boundary. For the heat generation simulation, I create a new topic to discuss.