Multiscale Device Simulator

In this first 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:

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 (http://geuz.org/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.

Modeling of the device with GMSH

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

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.

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

TiberCAD input file

• Device Section »
• Model section »
• Solver section »
• Physics section »
• Simulation section »
• Output »

1. Device structure (UP)

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

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.

2. Simulation models (UP)

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 simulation models.

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) , Macrostrain (to calculate macroscopical strain) and others.

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

Here, we declare that the simulation to be created will belong to the model class driftdiffusion (model driftdiffusion)
In the next Options block, we define the name of the TiberCAD simulation (simulation_name = driftdiffusion_1) and the TiberCAD regions to which it will be applied (physical_regions = all, where all means the whole device)

In this example, just one simulation for the driftdiffusion model is created. In general, several TiberCAD simulations belonging to the same model can be created; each of them must have a different name. This name must also be different from the model name. As we will see in the following, it is possible to specify physical and solver parameters for a single TiberCAD simulation, by referring to its name. Anyway, it is also possible to specify parameters common to all the TiberCAD simulations belonging to the same model, by referring to the name of the TiberCAD model instead (in this example driftdiffusion)

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

Next, some physical_model blocks follow.
The physical_model blocks are used in general to define the physical model to be used to describe a physical property or quantity associated to the present TiberCAD model.

Here, for example, we choose a Schottky-Read-Hall model ( model  =  srh) to describe recombination in our driftdiffusion calculations.

Also, we choose a field dependent model to describe electron and hole mobility (model = field_dependent), with the low field mobility determined by a doping-dependent model (Masetti model, see reference manual) (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 BC_Regions section.

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

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

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, specified in the sweep block (see after). [0.0] means that the default voltage value is 0.0 V.

3. Solver parameters (UP)

Now, the parameters for the solver are defined. The block name driftdiffusion means that all which is written in the block is valid for all the TiberCAD simulations belonging to the TiberCAD model driftdiffusion (see above).
Alternatively, the block name can be the name of one of the defined TiberCAD simulations.

These values are tuned for the specific simulation and usually don't need to be changed; however you can refer to the manual for further details.

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

4. Physics parameters (UP)

In the Physics section usually some physical parameters are defined.
In this case we state that the Fermi-Dirac statistic will be applied (statistics =  FD).

5. Run simulations (UP)

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.

Dimension of the problem is defined by dimension = 1.
Results of the calculation will be in the format of xmgr data visualization program (output_format = grace), suited to 1D case, since it consists of ascii text files with column data.
In plot we define the output variables to be calculated.

Note that the output variable ContactCurrents is needed to get in output the IV characteristic.

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.

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:

DriftDiffusion (name: driftdiffusion)
Solving equilibrium
......
iterations: 6, residual = 1.95054e-16
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:

Sweep value Vb = 0.05

contact name: contact voltage: contact current:
cathode 0 6.31825e-05
anode 0.05 -6.3187e-05

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

Output (UP)

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

Thus, in driftdiffusion_materials.dat we have the information about the material regions of the device (not very interesting in 1D, but useful in 2 and 3D).

In driftdiffusion_nodal_Vb_0.000.dat and all the other files for each bias step, we have the output for the nodal 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.

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

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