In this tutorial we will show a simple 1D Si pn diode example.
We will calculate the IV characteristic of the pn diode, by solving the Poisson equation and calculating the current with a drift-diffusion scheme.

The semi-classical transport simulation of electrons and holes is based on the driftdiffusion approximation.
Beside the electric potential, the electro-chemical potentials are used as variables such that the system of PDEs to be solved reads as follows:

equation 1

equation 2

equation 3

where P is the electric polarization due to e.g. piezoelectric effects and R is the net recombination rate, i.e. recombination rate minus generation rate.

In order to execute correctly the example you should have the following files in the same working directory:
diode_1D.tib : input file for TiberCAD
diode_1D.msh : mesh file produced by GMSH
GMSH is an automatic three-dimensional finite element mesh generator, please refer to GMSH manual for any details.

The mesh has been obtained from the GMSH geo script file diode_1D.geo.



In this script, first the geometrical entities Points and Lines are defined. Then, two physical regions are defined, each associated to one of the two geometrical entities: Physical Line("p_side") and Physical Line("n_side").

Physical Line("p_side") = {1};
Physical Line("n _side")= {2};

In this way, these physical regions are made available for tiberCAD, and will be associated to a tiberCAD Region (see in the following).

Then, two Physical Points are defined, anode and cathode, and associated to the first and to the last point of our 1D device.

Physical Point("anode") = {1};
Physical Point("cathode") = {3};

These points are needed to impose some boundary conditions and in this way they are made available for TiberCAD, and will be associated each to a boundary condition region (see in the following).

A non-interactive mode is also available in GMSH, without graphical user interface. For example, to mesh this first 1D tutorial in non-interactive mode, just type:

gmsh diode.geo -1 -o diode.msh

where diode.geo is the geometrical description of the device with GMSH syntax.

Some command line options are

  • -1 , -2, -3 to perform 1D, 2D or 3D mesh generation
  • -o mesh_file.msh to specify the name of the mesh file to be generated



Let's have a look to the input file; for further details you can refer to the reference manual.


  meshfile = diode_1D.msh
  material = Si
  Region n_side
      Nd = 1e18
      type = donor

  Region p_side
      Nd = 1e18
      type = acceptor

First we define the device structure: it is composed by two Regions, both constituted of the material Silicon, respectively n and p doped with a concentration of 1e18 cm-3.
With Region n_side and Region p_side the physical regions previously defined in GMSH are associated to the respective tiberCAD Region.



Now, we have to define the simulations which we are going to execute in our calculations.

One simulation is a particular set of equations, boundary conditions, physical parameters, solver parameters, which describes a physical problem to be solved by TiberCAD.
A valid tiberCAD simulation must belong to one of the predefined tiberCAD Modules.

TiberCAD, being a multiphysics software tool, includes several simulation models, each one describing a physical problem to be solved, e.g driftdiffusion (to solve Poisson and DriftDiffusion equations), efaschroedinger (to solve Schroedinger equation in effective mass approximation), Elasticity (to calculate macroscopical strain and ..) and others.

To create a TiberCAD simulation, we first have to declare the tiberCAD Module to which our simulation belongs:

Module driftdiffusion
  name = dd

Here, we declare that the simulation to be created will belong to the Module driftdiffusion (Module driftdiffusion)
Then we define the name of the TiberCAD simulation (name = dd) In this example, just one simulation for the driftdiffusion Module is created. In general, several tiberCAD simulations belonging to the same Module can be created; each of them must have a different name.

Module driftdiffusion
  name = dd
  plot = (Ec, Ev, Eg, ElField, ElPotential, eQFermi, hQFermi, eDensity,
  hDensity, NetRecombination, eCurrentDensity, hCurrentDensity,
  CurrentDensity, ContactCurrent, eMobility, hMobility, IonizedDonors,
      statistics = boltzmann
# use a field dependent mobility model with the
# gradient of the electrochem. potential as field
# parameter
    mobility field_dependent
      low_field_model = doping_dependent
# use SRH and Auger recombination
    recombination (srh, auger) {}

Thus, in this simulation we are going to calculate Poisson and Drift-diffusion (Module driftdiffusion) for all the device.

Next, the Physics section follows.
Physics is used in general to define the physical model to be used to describe a physical property or quantity associated to the simulation.

Here, for example, we choose a Schottky-Read-Hall (SRH) and Auger models to describe recombination in our driftdiffusion calculations.

recombination (srh, auger)

Also, we choose a field dependent model to describe electron and hole mobility (keyword field_dependent), with the low field mobility determined by a doping-dependent model (Masetti model, see reference manual) (low_field_model = doping_dependent).

mobility field_dependent
  low_field_model = doping_dependent
Now, we have to specify the Boundary conditions for the PDE problem represented by the simulation dd: they are defined in Contact blocks.

Here, two boundary conditions are defined, corresponding to the two contacts of our pn diode, anode and cathode.

Contact anode
  type = ohmic
  voltage = $Vb[0.0]
Contact cathode
  type = ohmic
  voltage = 0.0

Here the anode and cathode Contacts are associated to the corresponding anode and cathode regions defined in GMSH as Boundary condition regions (which in the present case of 1D simulation are actually just a point).

Besides, the type of Boundary condition is assigned (in this case ohmic) with its (initial) value, given by voltage.

Note that the anode voltage is expressed by the notation $Vb[0.0].
This means that the anode voltage will be given the value of the variable Vb, which will be specified in the sweep block (see below). [0.0] means that the default voltage value is 0.0 V.



Module sweep
  variable = $Vb
  solve = dd
  start = 0.0
  stop = 1.2
  steps = 60
  plot_data = true # to write data for each step

The sweep block allows us to perform the calculation of the IV characteristic of the diode: it uses the variable $Vb, defined in driftdiffusion, which will assume a range of values between 0 and 1.2 V. These values will be assigned in turn to the anode contact, in order to execute a series of complete dd simulations for each bias condition.



Finally, in Simulation section, we define the actual simulation to be performed, specified by

solve = sweep

this determines the execution of dd inside the sweep cycle;
this means that we are going to perform a calculation of Poisson and drift-diffusion for all the biases defined in the sweep section.

Results of the calculation will be in ascii column data format (output_format = grace), suited to 1D case, since it consists of ascii text files with column data.

  verbose = 2
  solve = sweep
  resultpath = diode1D_results
  output_format = grace


Now we can run TiberCAD....

tibercad diode.tib

During the execution, the screen output shows some information about how the simulation is going:

First, about the creation of the Device regions; if there is an inconsistency between the mesh file and input file device description, here an error message is displayed and the execution terminated.

Then, the defined Modules are initialised and the selected options and available output variables are shown on screen

If all is ok, the calculation is started:
We solve: sweep
In this example only one simulation is performed, sweep, which is a set of drift-diffusion calculations.

First, the equilibrium solution is found, by solving Poisson equation:

Ramping to sweep value $Vb = 0
INFO: I will first solve a currently unsolved system (dd):

>>>>---------------------------------- (dd) ---------------
Setting max_iterations for nonlinear solver to 150
Solving equilibrium
it 1, |du| = 11.0349, |r| = 1.88523 alpha = 1
it 2, |du| = 5.30845, |r| = 0.199306 alpha = 1
it 3, |du| = 1.29314, |r| = 0.0264031 alpha = 1
it 4, |du| = 0.152896, |r| = 0.000721821 alpha = 1
it 5, |du| = 0.0034889, |r| = 5.61159e-007 alpha = 1

iterations: 6, |du| = 2.57687e-006, |r| = 3.61564e-013
Equilibrium done

Then, current continuity equations are solved together with Poisson eq., for each step of the sweep.
The voltage and current found for each contact is displayed:

Ramping to sweep value $Vb = 0.02
Trying $Vb = 0.02

>>>>---------------------------------- (dd) ------------
it 1, |du| = 0.773634, |r| = 0.000488238 alpha = 1
it 2, |du| = 0.0769738, |r| = 2.33393e-005 alpha = 1
it 3, |du| = 0.00280745, |r| = 5.16137e-008 alpha = 1

iterations: 4, |du| = 4.23387e-006, |r| = 5.19017e-012

contact name: voltage: current:
anode 0.02 -2.84824e-009
cathode 0 2.69305e-009


The simulation ends correctly (hopefully...) with the Sweep value Vb = 1.2.

Ramping to sweep value $Vb = 1.2
Trying $Vb = 1.2

>>>>---------------------------------- (dd) --------------
it 1, |du| = 0.773634, |r| = 0.00968341 alpha = 1
it 2, |du| = 0.00292795, |r| = 2.16939e-005 alpha = 1

iterations: 3, |du| = 2.76748e-005, |r| = 7.49219e-007

contact name: voltage: current:
anode 1.2 -2526
cathode 0 2526

Solve time: 0h 0m 0.3s

Solve time: 0h 0m 38s

Total solve time: 0h 0m 38s
Simulation finished...



Now, the output directory contain the simulation results, as defined in output_format and plot.


In dd_Vb_0_msh.dat and all the other files for each bias step, we have the output for the quantities which have been calculated, e.g. conduction and valence bands, (quasi)fermi levels, electron and hole density and mobility.

Here is, for example, the band profile obtained at the last step of calculation, for V = 2V.

OUTPUT :: band profile for V = 2V


The IV characteristic of the diode is contained in the file sweep_dd.dat, where the currents at both the contacts are reported.

OUTPUT :: IV characteristic


In the file dd_Vb_0_msh.dat etc., you can find also the other quantities which have been calculated, in this case the electron, hole and total current density (Jn_x, Jp_x , J_x).