Transient and DC simulation of an illuminated photodiode

[Frangipani]

Hello everyone,

I recently started using Devsim. I have some trouble with this library.

Platform :

Windows with winpython

Goal :

Simulating the DC and transitional behavior of an illuminated PIN diode.

For the transient simulation, I want to create this circuit :

First problem :

When I connect my diode directly to the voltage source (with no resistor and no illumination). It seems to work in DC and tr ( up to approx. 0.6 V). I tried several things to connect the diode cathode to the voltage source and the anode to the resistor. However, I’m having trouble getting the right control combinations to make the connection properly. Could you tell me how to make the connection to get the given diagram, please?

Second problem :
I want to create a lighting profile that varies over time. To do this, I thought I’d change the parameters p1 and n1 via the set_parameter command during simulation. This should simulate the creation of electron-hole pairs. Do you think this method will work?

import os
os.environ["DEVSIM_MATH_LIBS"]="C:../scripts/Devsim/Library/bin/mkl_rt.2.dll"


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

import test_common
import diode_common
import numpy             as np
import matplotlib.pyplot as plt


reset_devsim()


#Physical constants
q      = 1.6e-19         #Coulomb
k      = 1.380650e-23    #J/K 
eps_0  = 8.85e-14        #F/cm^2
eps_Si = 11.1            #
eps_Ox = 3.9             #
n_i    = 1e10            # /cm^3
mu_n   = 400
mu_p   = 200
tau_n  = 1e-8            #s
tau_p  = tau_n/3         #s
T      = 300             #K

#Doping
Na_P   = 1e18            # /cm^3
Na_I   = 1e15            # /cm^3
Nd_N   = 1e18            # /cm^3
Nd_I   = 1e15            # /cm^3

#Photodiode schematic
"""
┌───────────────────────┐top
│          P+           │top_prime
│                       │PI_prime
│───────────────────────│PI
│          I            │PI_second
│                       │
│───────────────────────│Itop
│                       │
│                       │
│                       │
│                       │
│                       │
│                       │
│          I            │Imiddile
│                       │
│                       │
│                       │
│                       │
│                       │
│                       │
│                       │
│───────────────────────│Ibot
│          I            │
│                       │NI_second
│───────────────────────│NI
│          N+           │NI_prime
│                       │bot_prime
└───────────────────────┘bot
"""

top     =    0      #m
top_prime = 0.5e-6
PI_prime = 4e-6
PI      =    5e-6   #m
PI_second = 6e-6
Itop   =    50e-6  #m
Imiddle =    107.5E-6
Ibot    =    165e-6 #m
NI_second = 214e-6
NI      =    215e-6 #m
NI_prime = 216e-6
bot_prime = 219.5e-6
bot     =    220e-6 #m
 



name = "PIN_Diode_1D" 

def CreateMesh(device, region):
    ps=100e-9
    create_1d_mesh(mesh="dio")
    add_1d_mesh_line(mesh="dio", pos=top        , ps=ps/1000, tag="top"  )
    add_1d_mesh_line(mesh="dio", pos=top_prime  , ps=ps/100, tag="top_prime"  )
    add_1d_mesh_line(mesh="dio", pos=PI_prime   , ps=ps/100, tag="PI_prime"   )
    add_1d_mesh_line(mesh="dio", pos=PI         , ps=ps/1000, tag="PI"   )
    add_1d_mesh_line(mesh="dio", pos=PI_second , ps=ps/100, tag="PI_second")
    add_1d_mesh_line(mesh="dio", pos=Itop      , ps=ps/10, tag="Itop")
    add_1d_mesh_line(mesh="dio", pos=Imiddle    , ps=ps/1, tag="Imiddle")
    add_1d_mesh_line(mesh="dio", pos=Ibot       , ps=ps/10, tag="Ibot" )
    add_1d_mesh_line(mesh="dio", pos=NI_second , ps=ps/100, tag="NI_second")
    add_1d_mesh_line(mesh="dio", pos=NI         , ps=ps/1000, tag="NI"   )
    add_1d_mesh_line(mesh="dio", pos=NI_prime   , ps=ps/100, tag="NI_prime"   )
    add_1d_mesh_line(mesh="dio", pos=bot_prime  , ps=ps/100, tag="bot_prime"  )
    add_1d_mesh_line(mesh="dio", pos=bot        , ps=ps/1000, tag="bot"  )
    
    add_1d_contact(mesh="dio", name="top",  tag="top", material="metal")
    add_1d_contact(mesh="dio", name="bot",  tag="bot", material="metal")
    add_1d_region(mesh="dio", material="Si", region=region, tag1="top", tag2="bot")
    finalize_mesh(mesh="dio")
    create_device(mesh="dio", device=device)   

def SetNetDoping(device,region):
    #Step profile
    # CreateNodeModel(device, region, "Acceptors", str(Na_P)+"*step("+str(PI)+"-x)+tanh(0)+"+str(Na_I))
    # CreateNodeModel(device, region, "Donors"   , str(Nd_N)+"*step(x-"+str(NI)+")+"+str(Nd_I))

    #Soft step profile
    # CreateNodeModel(device, region, "Acceptors", "0.5*"+str(Na_P)+"*erf(1000000*(-x+"+str(PI)+"))+0.5*"+str(Na_P)+"+"+str(Na_I))
    # CreateNodeModel(device, region, "Donors", "0.5*"+str(Nd_N)+"*erf(1000000*(x-"+str(NI)+"))+0.5*"+str(Nd_N)+"+"+str(Nd_I))
     
    #Exp profile
    CreateNodeModel(device, region, "Acceptors", str(Na_P)+"*exp(-10*x/"+str(PI)+")+"+str(Na_I))
    CreateNodeModel(device, region, "Donors", str(Nd_N)+"*exp(10/"+str(PI)+"*(1*x-"+str(bot)+"))+"+str(Nd_I))
    
    CreateNodeModel(device, region, "NetDoping", "Donors-Acceptors")

def InitialSolution(device, region, circuit_contacts=None):
    # Create Potential, Potential@n0, Potential@n1
    CreateSolution(device, region, "Potential")
    # Create potential only physical models
    CreateSiliconPotentialOnly(device, region)

    # Set up the contacts applying a bias
    for i in get_contact_list(device=device):
        if circuit_contacts and i in circuit_contacts:
            CreateSiliconPotentialOnlyContact(device, region, i, True)
        else:
            ###print "FIX THIS"
            ### it is more correct for the bias to be 0, and it looks like there is side effects
            set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
            CreateSiliconPotentialOnlyContact(device, region, i)


def DriftDiffusionInitialSolution(device, region, circuit_contacts=None):
    ####
    #### drift diffusion solution variables
    ####
    CreateSolution(device, region, "Electrons")
    CreateSolution(device, region, "Holes")
    ####
    #### create initial guess from dc only solution
    ####
    set_node_values(
        device=device, region=region, name="Electrons", init_from="IntrinsicElectrons"
    )
    set_node_values(
        device=device, region=region, name="Holes", init_from="IntrinsicHoles"
    )

    ###
    ### Set up equations
    ###
    CreateSiliconDriftDiffusion(device, region)
    for i in get_contact_list(device=device):
        if circuit_contacts and i in circuit_contacts:
            CreateSiliconDriftDiffusionAtContact(device, region, i, True)
        else:
            CreateSiliconDriftDiffusionAtContact(device, region, i)

def SetSiliconParameters(device,region) :
    set_parameter(device=device, region=region, name="Permittivity", value=eps_Si*eps_0)    # Permittivity
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)             # ElectronCharge
    set_parameter(device=device, region=region, name="taun", value=tau_n)                   # taun
    set_parameter(device=device, region=region, name="taup", value=tau_p)                   # taup
    set_parameter(device=device, region=region, name="mu_n", value=mu_n)                    # mu_n
    set_parameter(device=device, region=region, name="mu_p", value=mu_p)                    # mu_p
    set_parameter(device=device, region=region, name="n_i", value=n_i)                      # n_i
    set_parameter(device=device, region=region, name="n1", value=n_i)                       # n1
    set_parameter(device=device, region=region, name="p1", value=n_i)                       # p1
    set_parameter(device=device, region=region, name="V_t", value=k*T/q)                    # V_t
    set_parameter(device=device, region=region, name="T", value=T)                          # T
    

    set_parameter(name="extended_solver", value=True)
    set_parameter(name="extended_model",  value=True)
    set_parameter(name="extended_equation",value=True)



"""
Simulation
"""
device = name
region = "Region"
CreateMesh(device,region)
SetSiliconParameters(device,region)
SetNetDoping(device,region)



InitialSolution(device, region)
solve(type="dc", absolute_error=1, relative_error=1e-20, maximum_iterations=100)



DriftDiffusionInitialSolution(device,region)
solve(type="dc", absolute_error=0.01, relative_error=1e-12, maximum_iterations=100)

#DC simulation
v = 0.59
set_parameter(device=device, name=GetContactBiasName("top"), value=v)
set_parameter(device=device, name=GetContactBiasName("bot"), value=0)
solve(type="dc", absolute_error=1, relative_error=1e-10, maximum_iterations=100)
PrintCurrents(device, "top")
PrintCurrents(device, "bot")

#Tran simulation
#time_step = 0.05e-3
#total_time = 1e-2
#current_time = 0
#
#while current<total_time:
#	solve(type="transient_dc", absolute_error=1, relative_error=1e-10, maximum_iterations=100)
#	solve(type=transient_bdf1",absolute_error=1, relative_error=1e-10, maximum_iterations=100)
#	current_time += time_step

Information :

I used the following example: Reverse recovery of pin diode - #3 by Saidon70

I have read the following forum articles :

Transient diode example

How to add optical generation (photovoltaics)

Ssac and transient in generation term

The calculated precision of Current

Reverse recovery of pin diode

I looked at the following example codes :

diode_1d.py, diode_commonpy, trand_diode .py, circ1.py, circ2.py, circ3.py, circ4.py, transient_circ.py, transient_circ2.py, transient_circ3.py.

Maybe I missed something obvious. I apologize in advance, but this is the first time I’ve used TCAD.

I’d be grateful for any help you could give me.

Have a nice day

[Juan]

Hi @Frangipani

Thanks for trying the software. devsim is not that user friendly and your efforts are impressive.

Please note that extended precision does not apply to the circuit boundary conditions, so you will need to be using a relative_error in your solve command at least 1e-15 or higher to get convergence.

Also setting the explicit path to the MKL should not be necessary if you are using a Python virtual environment or Anaconda conda environment.

Please let me experiment with your script and I’ll see if it works with no illumination. Unfortunately, I will not be able to help with the illumination part, but maybe some of the post authors, such as @simbilod or @GLuek ,may have some suggestions. I also suggest looking at this Thesis and project:

https://doi.org/10.15368/theses.2019.60
https://github.com/Bob95132/EE599-Thesis

[Juan]

also check out the “Projects” tab here:
https://docs.google.com/spreadsheets/d/11TpoCrNzKwWDDmjKtP1d5bCiwn8WYNJkm1ZleJI0_k4/edit?usp=sharing

[Juan]

it looks like your example is simulating at a bias without a circuit. Please see these examples for setting up a circuit bias:

examples/diode/ssac_diode.py
examples/diode/tran_diode.py

[GLuek]

Hi @Frangipani,

I have modified the example gmsh_diode3d.py which shows how to formulate circuit elements and connect them to the device (see below). In combination with the optical carrier generation from Ssac and transient in generation term, it should be all you need. I have succesfully simulated an image sensor with this approach.

# Copyright 2013 DEVSIM LLC
#
# SPDX-License-Identifier: Apache-2.0
from devsim import *
from devsim.python_packages.simple_physics import *
import devsim.python_packages.Klaassen as kla
import diode_common


def print_currents(device, contact):
    """
    Print out contact currents
    """
    ece_name = "ElectronContinuityEquation"
    hce_name = "HoleContinuityEquation"
    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
    if contact in ["top", "bot"]:
        voltage = get_circuit_node_value(solution='dcop', node=f"n{contact}")
    else:
        voltage = get_parameter(device=device, name=GetContactBiasName(contact))
    data = (voltage, electron_current, hole_current, total_current)
    print(f"{contact:10}{voltage:+.3e}\t{electron_current:+.3e}\t{hole_current:+.3e}\t{total_current:+.3e}")
    return data


def print_all_currents():
    """
    Prints all currents to console and returns outmap of values. !Currently only functional when all contacts are
    silicon contacts!
    :return:
    """
    print("\nSolution:")
    print("{0:10}{1:15}{2:12}{3:12}{4:10}".format("Contact", "Voltage", "Electron", "Hole", "Total"))
    print("                         Current     Current     Current")
    device_list = get_device_list()
    for device in device_list:
        contact_list = get_contact_list(device=device)
        outmap = {}
        for contact in contact_list:
                outmap[contact] = print_currents(device, contact)
    print_circuit_solution()
    return outmap


def print_circuit_solution():
    print('Circuit Solution')
    nodes = get_circuit_node_list()
    for node in get_circuit_node_list():
        r = get_circuit_node_value(solution='dcop', node=node)
        print("%s\t%1.15e" % (node, r))


set_parameter(name="debug_level", value="info")
#set_parameter(name="threads_available", value=8)
#set_parameter(name="threads_task_size", value=1024)
set_parameter(name="extended_solver", value="True")
set_parameter(name="extended_model", value="True")
set_parameter(name="extended_equation", value="True")

device="diode3d"
region="Bulk"

diode_common.Create3DGmshMesh(device, region)

# this is the devsim format
write_devices (file="gmsh_diode3d_out.msh")

diode_common.SetParameters(device=device, region=region)

v_source_top = circuit_element(name="VSourceTop", n1="ntop_res", n2="GND", value=0.0)
v_source_bot = circuit_element(name="VSourceBot", n1="nbot_res", n2="GND", value=0.0)

r_source_top = circuit_element(name="RSourceTop", n1="ntop_res", n2="ntop", value=1.0)
r_source_bot = circuit_element(name="RSourceBot", n1="nbot_res", n2="nbot", value=1.0)

c_source_top = circuit_element(name="CSourceTop", n1="ntop", n2="GND", value=1e-12)


# Connecting circuit to device ntop->top, nbot->bot, top and bot are defined in the mesh
circuit_node_alias(node="ntop", alias="top_bias")
circuit_node_alias(node="nbot", alias="bot_bias")

####
#### NetDoping
####
node_model(device=device, region=region, name="Acceptors", equation="1.0e18*step(0.5e-5-z);")
node_model(device=device, region=region, name="Donors",    equation="1.0e18*step(z-0.5e-5);")
node_model(device=device, region=region, name="NetDoping", equation="Donors-Acceptors;")
node_model(device=device, region=region, name="AbsDoping", equation="abs(NetDoping)")
node_model(device=device, region=region, name="logDoping", equation="log(AbsDoping)/log(10)")

diode_common.InitialSolution(device, region, circuit_contacts=["top", "bot"])


####
#### Initial DC solution
####
solve(type="transient_dc", absolute_error=1e-10, relative_error=1e-12, maximum_iterations=30)

###
### Drift diffusion simulation at equilibrium
###
diode_common.DriftDiffusionInitialSolution(device, region, circuit_contacts=["top", "bot"])
kla.Set_Mobility_Parameters(device, region)
kla.Klaassen_Mobility(device, region)

solve(type="transient_dc", absolute_error=1e-10, relative_error=1e-9, maximum_iterations=50)

v = -0.1
while v > -3.01:
    circuit_alter(name="VSourceTop", value=v)
    solve(type="transient_bdf1", absolute_error=1.0, relative_error=1e-09,
          maximum_iterations=100, tdelta=1e-8, charge_error=1.0)
    print_all_currents()
    v -= 0.1

print("\n\n\nSteady State: -----------------------------------------------------------")
for i in range(50):
    solve(type="transient_bdf1", absolute_error=1.0, relative_error=1e-09,
              maximum_iterations=30, tdelta=1e-8, charge_error=1.0)
    print_all_currents()

[Frangipani]

Hi @Juan and @GLuek,
thank you for your relevant answers.
I’ll take a look and try to implement it in my code.
Kind regards