# Problem with circuit elements and reverse biased pn-diode

Hi everyone,
I have encoutered a problem with the build-in circuit elements and a reverse biased pn-diode. When the diode is put in reverse with the set_parameter() strategy everything works fine. But, if you are doing the same with a circuit element voltage source, the current through the source explodes. I have build a MWE based on the diode3d devsim example to reproduce the problem, executed with devsim 2.6.0 and python 3.10:

#
from devsim import *
from devsim.python_packages.simple_physics import *
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))

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", n2="GND", value=0.0)
v_source_bot = circuit_element(name="VSourceBot", n1="nbot", n2="GND", value=0.0)

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

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

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

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

solve(type="transient_dc", absolute_error=1e10, relative_error=1e-8, maximum_iterations=50)

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

When you now look at the last solution which is printed to console, you can see that the current through the voltage source “VSourceTop” is extremly high in the 1e8 regime, while the current from the device behaves as expected:

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+4.008e-10	+1.568e-20	+4.008e-10
top       -1.000e+00	-3.273e-20	-3.992e-10	-3.992e-10
Circuit Solution
VSourceBot.I	-4.180081894749284e-10
VSourceTop.I	-1.304888058639660e+08
nbot	0.000000000000000e+00
ntop	-9.999999999999999e-01

Any thoughts about this? My goal at the end would be to create a diode in reverse as a photodiode and attach a capacitor which integrates the generated current, therefore it would be really nice to use the circuit elements.

Thanks a lot!

Hi @GLuek

Welcome, and thanks for providing such detailed information. I am able to reproduced the issue.

VSourceBot.I	-1.645353898689660e-10
VSourceTop.I	-5.238027693518443e+07
nbot	-6.310887241768094e-30
ntop	-9.999999999999999e-01

If sweep the “bot” bias positive, I get a similarly bad result at the other contact.

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +1.000e+00	+1.574e-10	+6.155e-21	+1.574e-10
top       +0.000e+00	-1.275e-20	-1.566e-10	-1.566e-10
Circuit Solution
VSourceBot.I	5.299975324150556e+07
VSourceTop.I	1.645353898689661e-10
nbot	9.999999999999999e-01
ntop	0.000000000000000e+00

If I change the solve type from “transient_bdf1” to “dc”, the issue goes away.

Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       -9.861e-32	+3.383e-20	+6.155e-21	+3.998e-20
top       -1.000e+00	-1.275e-20	-2.468e-19	-2.596e-19
Circuit Solution
VSourceBot.I	-3.872222930951233e-20
VSourceTop.I	3.872222930951466e-20
nbot	-9.860761315262648e-32
ntop	-9.999999999999999e-01

but the current is not too precise between the terminal current reported by get_contact_current and the circuit solution.

The diode 3d mesh quality is likely very bad, and may make it difficult to get a good solution across biases.

The transient solver may need some fine tuning for the problem that you are trying to solve.

I am happy to help with these issues as they come up.

I suggest that you start with a 2D device to see if the simulator is behaving in a manner expected for your problem. From there, it would be able to tackle other issues as they come up.

Hi Juan,
thanks for the quick reply and your support. I really enjoy working with the devsim package!

I tested some more options, so here comes a quick update:
I first encountered the problem in a “handcrafted” 3d gmsh mesh with some mesh optimization and varying density of mesh points with the same result, so I started to simplify the problem which lead to the MWE example above. A quick test with the 2d diode example shows the same behaviour.
What I noticed is, that resolving the same bias point leads to better results. So I split the tdelta of 1e-8 for one bias change into four steps with a tdelta of 0.25e-9 per bias point. Additionally, I activated the extended precision:

set_parameter(name="debug_level", value="info")
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"
#[...]

v = -0.1
while v > -1.01:
circuit_alter(name="VSourceTop", value=v)
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-09,
maximum_iterations=30, tdelta=0.25e-9, charge_error=1.0)
print("First: -----------------------------------------------------------")
print_all_currents()
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-09,
maximum_iterations=30, tdelta=0.25e-9, charge_error=1.0)
print("Second: ----------------------------------------------------------")
print_all_currents()
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-09,
maximum_iterations=30, tdelta=0.25e-9, charge_error=1.0)
print("Third: -----------------------------------------------------------")
print_all_currents()
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-09,
maximum_iterations=30, tdelta=0.25e-9, charge_error=1.0)
print("Fourth: ----------------------------------------------------------")
print_all_currents()
v -= 0.1

solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-06,
maximum_iterations=30, tdelta=1e-8, charge_error=1.0)
print_all_currents()
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-06,
maximum_iterations=30, tdelta=1e-8, charge_error=1.0)
print_all_currents()
solve(type="transient_bdf1", absolute_error=1e10, relative_error=1e-06,
maximum_iterations=30, tdelta=1e-8, charge_error=1.0)
print_all_currents()

which leads to the following output:

First: -----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+1.603e-08	+1.568e-20	+1.603e-08
top       -1.000e+00	-3.273e-20	-1.596e-08	-1.596e-08
Circuit Solution
VSourceBot.I	-1.671947125336367e-08
VSourceTop.I	-5.219552234558595e+09
nbot	0.000000000000000e+00
ntop	-9.999999999999999e-01

Second: ----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+1.850e-12	+1.568e-20	+1.850e-12
top       -1.000e+00	-3.273e-20	-2.869e-12	-2.869e-12
Circuit Solution
VSourceBot.I	-8.538716838052146e-13
VSourceTop.I	8.538716837982201e-13
nbot	0.000000000000000e+00
ntop	-9.999999999999999e-01

Third: -----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+2.607e-16	+1.568e-20	+2.608e-16
top       -1.000e+00	-3.273e-20	-5.161e-16	-5.161e-16
Circuit Solution
VSourceBot.I	-1.463354518653647e-16
VSourceTop.I	1.463354480513783e-16
nbot	0.000000000000000e+00
ntop	-9.999999999999999e-01

Fourth: ----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+1.260e-19	+1.568e-20	+1.417e-19
top       -1.000e+00	-3.273e-20	-1.606e-19	-1.933e-19
Circuit Solution
VSourceBot.I	-1.257382689847078e-19
VSourceTop.I	1.257317060540657e-19
nbot	0.000000000000000e+00
ntop	-9.999999999999999e-01

The first solve still shows the problem with the high current. The second solve leads to better results, since now the current from top and bottom source match, albeit there is still a difference between device and circuit currents. When I do the “steady state” solve at the end, everything starts to match:

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +0.000e+00	+8.366e-20	+1.568e-20	+9.934e-20
top       -1.000e+00	-3.273e-20	-6.661e-20	-9.934e-20
Circuit Solution
VSourceBot.I	-9.933941995934235e-20
VSourceTop.I	9.933952479401672e-20

Furthermore, when I ease the timing from nanosecond range to microseconds, I get matching results within one loop of the four solves.

Next I tried changing the transient solver type, “transient_bdf2” simply terminates with an error: “There was a floating point exception of type “Invalid” during LU Factorization Matrix factorization failed”. Using “transient_tr” has an interesing behaviour:

First: -----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       -2.758e-22	+3.774e-07	+1.568e-20	+3.774e-07
top       -1.000e+00	-3.277e-20	-3.749e-07	-3.749e-07
Circuit Solution
VSourceBot.I	-3.863950992429586e-07
VSourceTop.I	-1.066310100459579e+10
nbot	-2.758305545076182e-22
ntop	-9.999999999999999e-01

Second: ----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +1.420e-20	-3.772e-07	+1.568e-20	-3.772e-07
top       -1.000e+00	-3.269e-20	+3.746e-07	+3.746e-07
Circuit Solution
VSourceBot.I	3.862886235699354e-07
VSourceTop.I	1.066310100459579e+10
nbot	1.420024325896346e-20
ntop	-9.999999999999999e-01

Third: -----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       -1.420e-20	+3.770e-07	+1.568e-20	+3.770e-07
top       -1.000e+00	-3.277e-20	-3.743e-07	-3.743e-07
Circuit Solution
VSourceBot.I	-3.861822645090371e-07
VSourceTop.I	-1.066310100459579e+10
nbot	-1.420009706869676e-20
ntop	-9.999999999999999e-01

Fourth: ----------------------------------------------------------

Solution:
Contact   Voltage        Electron    Hole        Total
Current     Current     Current
bot       +1.420e-20	-3.768e-07	+1.568e-20	-3.768e-07
top       -1.000e+00	-3.269e-20	+3.740e-07	+3.740e-07
Circuit Solution
VSourceBot.I	3.860759811020647e-07
VSourceTop.I	1.066310100459579e+10
nbot	1.420009706869676e-20
ntop	-9.999999999999999e-01

The solution is relatively stable with the high current across the four steps, but the sign of the current flips with each step.

I will check some more simple structures with different mesh densities next.

Hi @GLuek

I believe I have found the issue. I put a bad model for charge at an ohmic contact, which doesn’t even make sense.

Can you please confirm this by editing:

\$CONDA_PREFIX/lib/python3.11/site-packages/devsim/python_packages/simple_physics.py

where python3.11 may be a different version in your installation. The patch is to edit this method to
CreateSiliconPotentialOnlyContact

so that the node_charge_model is removed from each invocation of the contact_equation command.

So this:

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

becomes this:

if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
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="",
edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")

I have this recorded here and will get the fix out as r2.6.1 once we have confirmed this is the issue.