[ Home ] [ News ] [ Contact ] [ Search ]

 ==> Download Software
 nextnano³ documentation

 Copyright notice
 About us
 Useful Links
 Publications
 
 * password protected

 

 
 

Linear current equation

In this chapter it is described how to set up the current matrix for the linear current equation:

The following topics are addressed:

  1. How current-clusters are defined and set up.
  2. How the current matrix is assembled.
  3. How boundary conditiond are implemented.
  4. How the current matrix is stored.
  5. Special treatment of current in 1D.
  6. How to calculate the matrix-vector product.

 

Current clusters

Before the current matrix can be created it must be defined

  • where the current equation is solved.
  • which boundary conditions on domain/inner boundaries to take.

The relevant information is written to MODULE current: [first half]

Current clusters are - as Poisson contact clusters - the Poisson contact regions defined via the input system. Poisson contact regions are a set of points on the material grid. All adjoining points on the physical grid are assigned the corresponding boundary condition. The Poisson contact region is NOT the area where the Poisson equation has to be solved, the Poisson contact region defines the contacts specified via the keyword $poisson-boundary-conditions.

MODULE current

num_current_clusters ...
Number of current clusters


current_clusterM(1..num_current_clusters)
has TYPE


 TYPE     :: current_cluster
  INTEGER :: num_points    Number of points in that cluster
  INTEGER,DIMENSION(1:num_points),POINTER :: num_guenther_points
                         Number of points on material grid which belong to that cluster.
  INTEGER :: kind_point = 1 --> Dirichlet condition for Fermi level
                                Fermi level set to value [V]
                        = 2 --> set to Dirichlet value,
Fermi equilibrium value in [V]
                        = 3 --> set to Neumann value, [A/m2] for current
  REAL(8),DIMENSION(3) :: valueV ...  see kind_point
  Note: The field is a vector  E = (Ex,Ey,Ez) = Ex * x  + Ey * y + Ez * z.


MODULE current


The Poisson assembler takes all 8 (3D) octants for each point
.
Wherever there is a transition from an octant on a Poisson cluster to an octant outside a Poisson cluster and a field is specified for that cluster, the corresponding component of the field is taken.
If any of the 8 octants has the option kind_point=2,3,4 (i.e. any kind of Dirichlet) the assembler ignores all fields and sets the potential on that point to the specified Dirichlet value.
For kind_point=4 --> the band edge of the adjoining octant is taken.
Note: This option is used for Schottky contacts, i.e. at a contact/semiconductor interface. So one has to take the band edge in the semiconductor, not in the contact.

If colliding potentials are specified the cluster with the lowest value
is valid.


The same procedure for boundary points on boundary of simulation domain:

 currentbnd_x3DM(1..2,j,k,4): currentbnd_x3DM(1,j,k) --> i=1
                              currentbnd_x3DM(2,j,k) --> i=n_1
                                                1..4 -->
interfaces 9,10,11,12
 currentbnd_y3DM(1..2,i,k,4): currentbnd_x3DM(1,i,k) --> j=1
                              currentbnd_x3DM(2,i,k) --> j=n_2
                                                1..4 -->
interfaces 5,6,7,8
 currentbnd_z3DM(1..2,i,j,4): currentbnd_x3DM(1,i,j) --> k=1
                              currentbnd_x3DM(2,i,j) --> k=n_1
                                                1..4 -->
interfaces 1,2,3,4


Note: Numbering of interfaces (1..12)
  interface   1 <--> between octant 1 and 5
  interface   2 <--> between octant 2 and 6
  interface   3 <--> between octant 3 and 7
  interface   4 <--> between octant 4 and 8
  interface   5 <--> between octant 1 and 3
  interface   6 <--> between octant 2 and 4
  interface   7 <--> between octant 5 and 7
  interface   8 <--> between octant 6 and 8
  interface   9 <--> between octant 1 and 2
  interface 10 <--> between octant 3 and 4
  interface 11 <--> between octant 5 and 6
  interface 12 <--> between octant 7 and 8

In 2D:
Note: Numbering of interfaces (1..4)
  interface 1 <--> between octant 1 and 2
  interface 2 <--> between octant 1 and 3                     3
  interface 3 <--> between octant 3 and 4                 2  +  4
  interface 4 <--> between octant 2 and 4                     1

 


           = 1 --> Dirichlet condition for Fermi level, Fermi level set to value [V]
           = 2 --> set to Dirichlet value, Fermi equilibrium value in [V]
           = 3 --> set to Neumann value, [A/m2] for current
         
 bndval_x3DM(1..2,j,k,4)
 bndval_y3DM(1..2,i,k,4)
 bndval_z3DM(1..2,i,j,4)
 
    ...  field, if    currentbnd_.3DM=1   (directed orthonormal to the surface)
         potential currentbnd_.3DM=2   ("." stands for x, y or z)

Note: Up to now, on simulation domain boundary only field zero is implemented.


Must be allocated and defined before call of currentmatrix3D and is deallocated in deallocate_current.


The current clusters and contact models are provided in the modules:
   mod_num_of_ohmics
   mod_num_of_schottkys
   mod_num_of_dirichlets
   mod_num_of_neumanns
   mod_num_of_mixed
   mod_apply_dirichlet_to_cluster
   mod_dirichlet_grid_points
   mod_potential_for_dirichlet
   mod_apply_ohmic_to_cluster
   mod_ohmic_grid_points
   mod_voltage_for_ohmic
   mod_apply_schottky_to_cluster
   mod_schottky_grid_points
   mod_voltage_for_schottky
   mod_barrier_for_schottky

See there for more details.

 

With that information subroutine define_current_clusters defines the current problem in MODULE current:

ohmic,     voltage controlled  --> When voltage is applied, Dirichlet condition (for potential and Fermi levels) with
pot   = pot_built_in      + V_applied
fermi = fermi_equilibrium - V_applied.
ohmic,     current controlled  --> not yet implemented
Schottky, voltage controlled  --> When voltage is applied, Dirichlet condition (for potential and Fermi levels) with
pot   = pot_built_in      + V_applied
fermi = fermi_equilibrium - V_applied
.
Schottky, current controlled  --> not yet implemented
Dirichlet               --> not yet implemented
Neumann              --> not yet implemented
mixed, i.e. relation between Dirichlet
and Neumann (Cauchy)     -->
not yet implemented

Note: The grid points inside a current cluster (contact) are exempted from simulation (i.e. set to Dirichlet which is of no importance).

 

 

Setup of current matrix

After the Poisson clusters have been defined by subroutine define_current_clusters the current matrix is set up in SUBROUTINE currentmatrix3D (currentmatrix2D,currentmatrix1D).
The current solver allows one to exempt parts of the simulation domain from simulation (contacts). On these internal surfaces Neumann boundary conditions can be specified (zero by default).
The corresponding information is provided by SUBROUTINE irregular_boundary('Current',i,j,k,irr_bcC,which_octantsVL,fieldM):

Note: Up to now only Dirichlet boundary conditions on contacts are used.

  • SUBROUTINE irregular_boundary(i,j,k,irr_bcC,which_octantsVL,fieldM)

    (i,j,k) point on physical grid    [Input]
    which_octantsVL(1..8) .TRUE.  ---> octant within simulation region   [Output]
                          .FALSE. ---> octant outside simulation region, i.e. inside contact

    irr_bcC = 'in'  If all octants within simulation region.
              'ot'  If all octants outside simulation region.
              'on'  If at least one octant within simulation region, but not all.

    curM(1..8,1..3) : curM(n,1) field component in x-direction for box n
                      curM(n,2) field component in y-direction for box n
                      curM(n,3) field component in z-direction for box n
                                [units in A/m2]

NOTE: Up to now only Dirichlet conditions on contacts are used. So curM is zero.


Finally, the program in first SETUP always sets Neumann boundary conditions which are set to zero if none are defined. Later on, Dirichlet boundary conditions are set (i.e. Neumann are overwritten by Dirichlet).
The fields on the domain boundaries can be obtained via boundary_value:

  • FUNCTION boundary_value(ch,i,j,k,whichinterface) RESULT(bv)

    This may just be called for (i,j,k) on the domain boundary, otherwise stop.

    Input:
    (i,j,k) on physical grid
    ch = 'x' or 'y' or 'z'  (direction of differential)
    called at any setup of current matrix
    whichinterface=1..4 --> interfaces  9,10,11,12 for 'x'
                                    5, 6, 7, 8 for 'y'
                                    1, 2, 3, 4 for 'z'
    For whichinterface the same convention as for currentbnd_x3DM... is .TRUE.

    Output:
    If Neumann boundary condidtions are valid, the current is provided.
    Up to now: not used --> set to zero
    Otherwise:                       set to zero
    Note: In the first setup always Neumann conditions are implemented.
    They are overwritten with Dirichlet values (if specified) later on.

The mobility is provided by subroutine get_mobility3D(get_mobility2D). The mobility is needed on a point on the material grid. The specification is done as described above.

  • SUBROUTINE get_mobility3D(i1,j1,k1,oct1,i2,j2,k2,oct2,chargeC,mobility)

    Input:

      (i1,j1,k1) point on physical grid, oct1 octant(1..8)
      (i2,j2,k2) point on physical grid, oct2 octant(1..8)


                3  |  4                3  |  4
                   |                      |
             ---(i1,j1)---          ---(i2,j2)---
                   |                      |
                1  |  2                1  |  2


      chargeC = 'el' or 'hl'

    Output:

    Mobility between the specified two octants in [cm2/Vs]

    Note: They must be adjoining points and octant numbering must be consistent.

The recombination is provided by subroutine get_net_recombination. It refers to a point of the material grid.

  • SUBROUTINE get_net_recombination(i,j,k,oct,chargeC,net_recomb)

    Input:

         (i,j,k) point on Physical grid
         oct ... 1..8 number of octant

                   3   |   4
                       |
                 ----(i,j)----
                       |
                   1   |   2

       chargeC = 'el' or 'hl'

    Output:

    net_recomb ... net recombination rate in [e/m3/s]
    Comprises recombination and generation terms.
    - is positive for effective generation
    - is negative for effective recombination

    Note:
    This routine is used for option 'DriftDiffAll'.
    This means that one quasi-Fermi level is assumed for all carriers (which is the only option up to now): One Fermi level for all electrons and another Fermi level for all holes.

The computational molecule is assembled by first calculating the computational molecules for all eight octants and adding them later on. Subroutine set_cur_differential3D (set_cur_differential2D, set_cur_differential1D, FUNCTION differential (1D only)) provides the computational molecule for a specified octant.

  • SUBROUTINE set_cur_differential3D(i,j,k,oct,chargeC,modeC,starV)

    Defines differential for current equation for octant oct of point (i,j,k).

    Input:
     i,j,k point on physical grid
     oct 1...8 number of octant
     modeC

    NOTE: Numbering of octants

             3   |   4
                 |
         -----(i,j,k-)-----
                 |
             1   |   2


             7   |   8
                 |
         -----(i,j,k+)-----
                 |
             5   |   6


       chargeC = 'el' or 'hl'


    Output:
    The computational molecule is provided in starV.
    At first, each octant has its own numbering of points as specified below. For each octant there are four neighboring points that are needed for finite differences. Their numbering is as follows:

     octant1 (i  ,j  ,k-1) =  1  (i-,j-,k-)
             (i  ,j-1,k  ) =  2
             (i-1,j  ,k  ) =  3
             (i  ,j  ,k  ) =  4

     octant2 (i  ,j  ,k-1) =  5  (i+,j-,k-)
             (i  ,j-1,k  ) =  6
             (i+1,j  ,k  ) =  7
             (i  ,j  ,k  ) =  8

     octant3 (i  ,j  ,k-1) =  9  (i-,j+,k-)
             (i  ,j+1,k  ) = 10
             (i-1,j  ,k  ) = 11
             (i  ,j  ,k  ) = 12

     octant4 (i  ,j  ,k-1) = 13  (i+,j+,k-)
             (i  ,j+1,k  ) = 14
             (i+1,j  ,k  ) = 15
             (i  ,j  ,k  ) = 16

     octant5 (i  ,j  ,k+1) = 17  (i-,j-,k+)
             (i  ,j-1,k  ) = 18
             (i-1,j  ,  k) = 19
             (i  ,j  ,k  ) = 20

     octant6 (i  ,j  ,k+1) = 21  (i+,j-,k+)
             (i  ,j-1,k  ) = 22
             (i+1,j  ,k  ) = 23
             (i  ,j  ,k  ) = 24

     octant7 (i  ,j  ,k+1) = 25  (i-,j+,k+)
             (i  ,j+1  ,k) = 26
             (i-1,j  ,k  ) = 27
             (i  ,j  ,k  ) = 28

     octant8 (i  ,j  ,k+1) = 29  (i+,j+,k+)
             (i  ,j+1,k  ) = 30
             (i+1,j  ,k  ) = 31
             (i  ,j  ,k  ) = 32


    right handside: octant1 = 33
               octant2 = 34
               octant3 = 35
               octant4 = 36
               octant5 = 37
               octant6 = 38
               octant7 = 39
               octant8 = 40


     starV(1:40)

     if oct = 1 --> starV( 1: 4) and starV(33) is defined, rest=0
           = 2 --> starV( 5: 8)
    and starV(34) is defined, rest=0
           = 3 --> starV( 9:12)
    and starV(35) is defined, rest=0
           = 4 --> starV(13:16)
    and starV(36) is defined, rest=0
           = 5 --> starV(17:20)
    and starV(37) is defined, rest=0
           = 6 --> starV(21:24)
    and starV(38) is defined, rest=0
           = 7 --> starV(25:28)
    and starV(39) is defined, rest=0
           = 8 --> starV(29:32)
    and starV(40) is defined, rest=0

    The subroutine is called for each octant and the corresponding eight parts of the computational molecule are defined. Later on they are added to build the final computational molecule.


    Discretizes the differential operator
     - int [div j]  =  - surface_integral j dA  ==  rhs
    where   j = mobility * dens * grad(EF)         [A/m2]
         rhs = int recombination

    Note:  In this first version we simply use the density for dens.
         dens > 0, C/m3

    Units:
             mobility  --> m2/Vs
             dens    --> C/m3 = As/m3
             grad E_F --> V/m --> j = A/m2
    Note: EF is given in [V].

The matrix which is created at the setup is stored in MODULE current:

 

diagV             .. DIMENSION(1:lengvecV)
contains diagonal elements of current matrix


matrixindicesM    .. DIMENSION(1:nummaelemente,1:2)
contains for nonvanishing offdiagonal matrix elements row and column index
Because matrix is symmetric, only upper (triangle) half is stored.


matrixwerteV      .. DIMENSION(1:nummaelemente)
contains corresponding values



matrixposV(i)     .. DIMENSION(1..lengvecV)
matrix_pos(n)
points to the position, where the nonvanishing elements of row n (upper triangle) are placed in the arrays.
matrixindicesM and matrixwerteV

If the line doesn't have any offdiagonal elements it points to the next entry unequal to zero.


nummaelemente     .. number of offdiagonal matrix elements


bV                .. DIMENSION(1:lengvecV)
contains the complete right hand side of the current matrix equation
 


The setup is done in SUBROUTINE currentmatrix3D (currentmatrix2D,currentmatrix1D):

The program loops through all points:

  1)  CALL WhichKindOfBoundaryPoint3D(ch,i,j,k) [internal subroutine, in 2D/1D this distinction is done explicitly]

       -->     ch ='in', if  (i,j,k) is within domain.
                otherwise --> ch provides information about the position of the boundary point (see source code).

       The variables  beteiligtL_V and randL_V provide information for box integration:
       beteiligtL_V(i)=.TRUE. <--> octant i is within simulation region, so box integration is performed over it.
       randL_V(1..3): says if flux through surface in x-, y-, z-direction has to be considered in box integration.

   2)  CALL irregular_boundary provides for (i,j,k) in irr_bcC the following information:
      irr_bcC = 'ot' <--> (i,j,k) lies in a region which is exempted from simulation (i.e. a Poisson contact region)
                          --> set all beteiligtL_V=0
              = 'in' <--> (i,j,k) lies within simulation region (i.e. not in a Poisson contact region)
                          --> do not alter beteiligtL_V
              = 'on' <--> (i,j,k) lies on internal interface separating
                              simulation / nonsimulation region
(Poisson contact region)
                          --> set the corresponding octants of beteiligtL_V to 0

  3) Now the setup of matrix elements for that point:

      a)   if outer point (irr_bcC=='ot') --> CALL outer_point, writes for (i,j,k) a Dirichlet condition (=0)

      b)   if inner point which is no multiple point and does not adjoin a multiple point and not to a Poisson contact region:
            the differential is defined explicitly: CALL standardfall3D
      c)   else: perform box integration over all boxes specified by beteiligtL_V(:)=.TRUE.:
      CALL CalculateCoefficient3D

  4) Now allocate and define ska_bnd_x3DM, ska_bnd_y3DM, ska_bnd_z3DM.

 

 

Technical details of computational molecule

First the computational molecule is extracted for each of the eight octants. See description of subroutine set_cur_differential3D for further information. Finally the computational molecules are added to get the final computational molecule:

        1             2               3           4           5             6            7        |           8
    ---------------------------------------------------------------------------------------------------------------------
    (i,j,k-1)    (i,j-1,k)     (i-1,j,k)     (i,j,k)    (i+1,j,k)    (i,j+1,k)    (i,j,k+1)   |   right hand side

The column l=1,..,8 corresponds to koeffsum_V(l):

!--------------
!
Addition of boxes
!--------------


koeffsum_V(1) = koeff_M( 1) + koeff_M( 5) + koeff_M( 9) + koeff_M(13)
koeffsum_V(2) = koeff_M( 2) + koeff_M( 6) + koeff_M(18) + koeff_M(22)
koeffsum_V(3) = koeff_M( 3) + koeff_M(11) + koeff_M(19) + koeff_M(27)
koeffsum_V(4) = koeff_M( 4) + koeff_M( 8) + koeff_M(12) + koeff_M(16) + &
                koeff_M(20) + koeff_M(24) + koeff_M(28) + koeff_M(32)
koeffsum_V(5) = koeff_M( 7) + koeff_M(15) + koeff_M(23) + koeff_M(31)
koeffsum_V(6) = koeff_M(10) + koeff_M(14) + koeff_M(26) + koeff_M(30)
koeffsum_V(7) = koeff_M(17) + koeff_M(21) + koeff_M(25) + koeff_M(29)
koeffsum_V(8) = koeff_M(33) + koeff_M(34) + koeff_M(35) + koeff_M(36) + &
                koeff_M(37) + koeff_M(38) + koeff_M(39) + koeff_M(40)

Now the area of internal interfaces are defined (only for multiple points):

The matrix elements are written (just for the upper part of the matrix (symmetric) in CreateMatrix3D:
Creates the elements for single points and the first line of multiple points.
Note: For multiple points the balance from integration is fully contained in the first line.
         In Poisson equation the potential is assumed to be continuous.
         The vector of the Poisson equation is defined on Stefan's grid (1..lengvecV) which contains
         multiple point; so additional lines are added to account for that.
         Since the matrix must be symmetric and positive definite, we use following scheme:
 

Generate multiple points
Double points:
The matrix is supplemented by an identity relation for the second point. The following structure is obtained:
2nd line: -a + a = 0
Matrix must be symmetric, thus add + a and -a to 1st line.

  -1   0  ...  0   -1  4+a   -a   -1   0   0  ...  0   -1  ...   = Q
          ...  0    0   -a    a    0   0   0  ...          ...   = 0
      -1  ...  0    0   -1    0    4  -1   0  ...  0   -1  ...   = Q

Eightfold point:

-1 ... -1  0  ...  0  -1  6+a  -a   0   0   0   0   0   0  -1 ... -1 ... -1 ...  = Q
                           -a  7a  -a  -a  -a  -a  -a  -a                        = 0
                            0  -a   a   0   0   0   0   0                        = 0
                            0  -a   0   a   0   0   0   0                        = 0
                            0  -a   0   0   a   0   0   0                        = 0
                            0  -a   0   0   0   a   0   0                        = 0
                            0  -a   0   0   0   0   a   0                        = 0
                            0  -a   0   0   0   0   0   a                        = 0
    -1  ... -1  0  ... 0   -1   0   0   0   0   0   0   0   6  -1 .. -1 .. -1 .. = Q

(An alternative would be, but it didn't converge with conjugate gradient, maybe later check again...)
-1 ... -1  0  ...  0  -1  6+7a -a  -a  -a  -a  -a  -a  -a  -1 ... -1 ... -1
                           -a   a
                           -a   0   a   0   0   0   0   0
                           -a   0   0   a   0   0   0   0
                           -a   0   0   0   a   0   0   0
                           -a   0   0   0   0   a   0   0
                           -a   0   0   0   0   0   a   0
                           -a   0   0   0   0   0   0   a
    -1  ... -1  0  ... 0   -1   0   0   0   0   0   0   0   6  -1 .. -1 .. -1

 

 

Boundary Conditions (Dirichlet points):

After the matrix has been completely set up, Neumann boundary conditions are overwritten with Dirichlet values:

There is an auxiliary subroutine DIRICHLET3D(i,j,k,dirval) (2D/1D) which defines Dirichlet conditions for a certain point.
In detail: Matrix elements in that line are overwritten with zero, diagonal element is set to a default value
             and the
right hand side  bV is defined accordingly.
             Note: Matrix must remain symmetric.
                      So if you overwrite one line to a Dirichlet condition, you also have
                      to change the line which shares the matrix element in order to restore symmetry.

SUBROUTINE set_dirichlet_current sets the boundary conditions as specified in MODULE current.

  • SUBROUTINE set_dirichlet_current

    Sets Dirichlet points as specified in MODULE current.
    Current matrix must be defined.
    Multiple points are all set to the same value.
    The points with higher cluster number are valid.

    Loops through all current clusters and loops through all points of the current cluster and sets the according Dirichlet value.

    NOTE: Sets Dirichlet value for Fermi level, i.e. only to be used in connection with
    CurrentProblemC = 'solve-for-fermi-1' (LAPACK), 'solve-for-fermi-2' (conjugate gradient)

 

 

Current in 1D

In 1D, there is an additional option to calculate the current by directly integrating the current equation (MODULE mod_integrate_current1D).
(Other option: current_setup1D. Same as current_setup2D and current_setup3D.)
The alternative method, which directly integrates the current equation is done in subroutine integrate_current1D.
 

  • SUBROUTINE integrate_current1D(chargeC,i1,i2,fermiV,E_FL,E_FR, &
                                   curL,curR,curV,bndC,optionC)

    Solves transport equation by using the full density.
      j = mobility * density * grad(EF)
    One has only two quasi-Fermi levels (one for electrons, one for holes).

    For calculation of densities Fermi levels in MODULE fermi_level are valid.
    Note: You must set all of them (for quantized states) to the value in fermi_nV.

    Potential in MODULE potentials phiV is valid.

    chargeC = 'el' or 'hl'
    i1 , i2  ... >=1 , <= n_1 , i2>i1

    fermiV (1..lengvecV)
    Note: Fermi supposed to be continuous at multiple points but this can be changed to discontinuous Fermi levels
                         ... Output: Fermi level between points overwritten,
                                        remains untouched in other regions.

    E_FL,E_FR  (Input/Output) ... Fermi level in i1, i2
    curV              (Input/Output) ... Current (particle current, i1 --> i2 :: positive)
               (1..n_1-1)     curV(n) = current between n and n+1 on physical grid
    curR,curL  (Input/Output) ... Current that flows in or out.

    bndC   ...  'E_l_E_r'  ... left/right Fermi level as boundary condition
                'E_l_cur'  ... left Fermi level and current curL as boundary condition
                'E_r_cur'  ... right fermilevel and current curR as boundary condition

    Note: This subroutine only takes one iteration.
    EFL, EFR, Fermi levels in [J]
    cur in A/m2

    NOTE:
    Up to now, initial value for current is set to 1.0d0.
    It would make sense to use curV as input to extract a starting value.

 

This method integrates the equation:

  d/dx  exp(E_F/kT) = J / (mu*k*T*kappa)   , kappa = n exp(-E_F/kT)

(which is obtained from j = mu * n * grad E_F)

The integration is carried out from i1 to i2 (i1,i2 on physical grid).

The current is given on material grid points (i.e. between physical points) and it is stored in the variable
curV (1..n_1-1).

Recombination:   div j = d/dx j = -R  -->   j_R - j_L = - int R dx

The current in position i can be split up:  curV(i) = cur_0 + j_rec
                                                               with  j_recR - j_recL = - int R dx
NOTE: cur_0 is the current that enters the device from the left contact.

Thus we get for points i and i+1:

exp(E_F/kT)|i+1  -  exp(E_F/kT)| i  =  int_i-i+1 J / (mu*k*T*kappa) dx  =
= int_i_i+1 cur_0 / (mu*k*T*kappa) dx  +  int_i_i+1 j_rec/ (mu*k*T*kappa) dx

The added values of  int_i-i+1 cur_0/ (mu*k*T*kappa) dx are stored in exp_fermiV. (exp_fermiV(i1) = exp(E_FL/kT)),
the added values of  int_i_i+1 j_rec/ (mu*k*T*kappa) dx  in exp_fermi_recV. (exp_fermi_recV(i1) = 0).

At first, cur_0 is set to 1.0d0. Then the integration is carried out.
Finally cur_0 is scaled in a way that  exp_fermiV reaches the value exp(E_FR/kT) at the right which is the boundary condition.
Noteexp_fermiV(i1) = exp(E_FL/kT)

--> If one begins at point i and wants to rescale cur_0 such that exp_fermiV(end_point) = exp(E_FR/kT):

  exp(E_FR/kT) - exp(E_Fi/kT)  =  scale * int_i_end cur_0 / (mu*k*T*kappa) dx
  +  int_i_end j_rec / (mu*k*T*kappa) dx

--> scale = [exp(E_FR/kT) - exp(E_Fi/kT) - int_i_end j_rec / (mu*k*T*kappa) dx ] /  [int_i_end cur_0 / (mu*k*T*kappa) dx ]  =
    [exp(E_FR/kT) - exp(E_Fi/kT) - (exp_fermi_recV(end_point)-exp_fermi_recV(i)) / exp_fermiV(end_point) - exp_fermiV(i))


Recombination/generation and mobilities are obtained via the subroutines in MODULE get_transport_on_guent. A positive recombination rate means generation of electrons/holes.

This subroutine is called in the internal subroutine calculate_new_fermi_int_current in subroutine solve_current_problem1D.

 

 

Matrix-vector product

The subroutines that calculate the matrix-vector products are:
subroutine current_matvec, subroutine psolve_jac_current contained in MODULE current_routines.

 

   
Last modified: 09-Jun-2011