aeneas.c

Go to the documentation of this file.

/*
   This file belongs to Aeneas. Aeneas is a GNU package released under GPL 3.
   This code is a simulator for Submicron 3D Semiconductor Devices. 
   It implements the Monte Carlo transport in 3D tetrahedra meshes
   for the simulation of the semiclassical Boltzmann equation for both electrons.
   It also includes all the relevant quantum effects for nanodevices.

   Copyright (C) 2007 Jean Michel Sellier <sellier@dmi.unict.it>
                                          <aeneas@southnovel.eu>
 
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 3, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/*
 ************************************************************************

  Program Name              : Aeneas (GNU package)
  Version                   : 1.1
  File Name                 : aeneas.c
  Programmer                : Jean Michel Sellier
  Copyright                 : 2007, Jean Michel Sellier
  Input File                : specified by the user
  Output file               : various mesh, BB, vtk, 3D output files
  Created on                : 06 june 2007 Siracusa, Jean Michel Sellier
  Last modified on          : 16 september 2007 Siracusa, Jean Michel Sellier
  Description               : Monte Carlo simulation of the Wigner equation
                              on a 3D unstructured mesh.
  Dependencies              : None

  Comments                  : WARNING!!! The input mesh file is described in micron!!!

  Important bibliography    : Elements finis: theorie, applications, mise en oeuvre ,
                              Alexandre Ern, Jean-Luc Guermond, SMAI, Springer (french book).

                              Numerical Simulation of Submicron Semiconductor Devices,
                              Kazutaka Tomizawa, Artech House (english book).
  ************************************************************************
*/

#include<getopt.h>
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<memory.h>
#include<string.h>
#include<time.h>
//#include<gl/gl.h>
//#include<gl/glu.h>
//#include<gl/glut.h>


// the following will be changed once the dynamic allocation will be implemented
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#define MAXNUMELEMENTS 100000  // maximum number of elements per mesh
#define MAXNUMNODES 10000 // maximum number of nodes per mesh
#define NPMAX 5000000        // maximum number of particles
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

// read and write
#define READ  0
#define WRITE 1

// yes and no
#define NO  0
#define YES 1

// definition of various output format flags
#define POINT3D 0
#define VTK     1
#define MESH    2
#define ALL     100

// definition of the material reference table
#define NOAMTIA 13   // Number Of All Material Taken Into Account (excluding SiO2)
// ********************
#define SIO2 -1      // reference number to Silicon Oxide
#define SILICON 0    // reference number to Silicon
#define GAAS 1       // reference number to GaAs
#define GERMANIUM 2  // reference number to Germanium
#define INSB 3       // reference number to InSb
#define ALSB 4       // reference number to AlSb
#define ALXINXSB 5   // reference number to Al_x In_x Sb
#define ALXIN1XSB 6  // reference number to Al_x In_(1-x) Sb
#define ALAS 7       // reference number to AlAs
#define ALP 8        // reference number to AlP
#define GAP 9        // reference number to GaP
#define GASB 10      // reference number to GaSb
#define INAS 11      // reference number to InAs
#define INP 12       // reference number to InP
// ********************

// definition of the various possible contacts
#define NOCONTACT -1
#define INSULATOR 0
#define OHMIC     1
#define SCHOTTKY  2
#define CONTACT   3  // applied voltage -- with no carrier reservoir 
                     // and no carrier can go outside this contact

// maximum number of steps in the discretised energy
#define DIME 2003

// Definition of all the physical constants
// The system used is the following :
//   [L] = Angstrom
//   [M] = 1.e-31 Kg
//   [T] = picosecond
//   [C] = 1.e-19 Coulomb

double EPSR[NOAMTIA+1];
double EPF[NOAMTIA+1];
double alphaK[NOAMTIA+1][4];
double EMIN[NOAMTIA+1][4];
double HWO[NOAMTIA+1][6];
double DTK[NOAMTIA+1][6];
double ZF[NOAMTIA+1][6];
double RHO[NOAMTIA+1];
double DA[NOAMTIA+1];
double UL[NOAMTIA+1];
double EG[NOAMTIA+1];
double PII[NOAMTIA+1];
double ETH[NOAMTIA+1];

// permittivity of vacuum [C]*V^-1*Angstrom^-1 (this definition is used
// only in the assemble_matrix_A routine)
const double epslon0=8.854187817e-3;
// Relative Permittivity of Silicon Oxide (SiO2)
const double EPSR_SiO2=3.9;

double VIC[MAXNUMNODES+1]; // Volt on the I-th vertex (Contact)

const double thresh=1.e-16; // threshold value for linear sparse solver 

const double EPS=1.e-16;

const double TINY=1.e-16;

// \hat{\psi}_n(\hat{\xi}_l) see pag.322 of the french book
double base_geo[4][4]; // i=1..ng, j=1..lg

// \hat{\theta}_n(\hat{\xi}_l) see pag.322 of the french book
double base_ref[4][4]; // i=1..nf, j=1..lg

// \frac{\partial \hat{\psi}_n}{\partial \hat{x}_k}(\hat{\xi}_l) see pag.323 of the french book
double d_base_geo[3][4]; // k=1..dim, i=1..ng

// \frac{\partial \hat{\theta}_n}{\partial \hat{x}_k}(\hat{\xi}_l) see pag.323 of the french book
double d_base_ref[3][4]; // k=1..dim, i=1..nf

// The following declarations are related to the mesh
// *****
int Ng; // Number of vertices
int Ne; // Number of elements
int ng = 4; // Number of vertices per element
int nf = 4; // Number of degree of freedom
int dim = 3; // dimension of the space
int i_front[MAXNUMNODES+1];
int i_dom[MAXNUMELEMENTS+1];
int noeud_geo[4][MAXNUMELEMENTS+1]; // i=1..ng
double coord[3][MAXNUMNODES+1]; // i=1..dim
// *****

// Global Gaussian quadrature weights
double poids_K[4][MAXNUMELEMENTS+1]; // i=1..lg, j=1..Ne

// elements of the inverse jacobian for the various transformations
double inv_jac_K[3][3][MAXNUMELEMENTS+1]; // i=1..dim, j=1..dim, k=1..Ne

// derivee d'une fonction de forme globale
double d_base_K[3][4][MAXNUMELEMENTS+1]; // i=1..dim, j=1..nf, m=1..Ne

// neighbourhood table
int vois[4][MAXNUMELEMENTS+1];

// the right hand side vector for the 3D Poisson equation
double b[MAXNUMNODES+1];

// the electrostatic potential calculated by Poisson 3D
double V[MAXNUMNODES+1];
double Vi[MAXNUMNODES+1]; // initial applied potential

// the electric field components
double Ex[MAXNUMNODES+1];
double Ey[MAXNUMNODES+1];
double Ez[MAXNUMNODES+1];

// for the linear solver part...
int ija[100*MAXNUMELEMENTS+1];
double sa[100*MAXNUMELEMENTS+1];

// Gauss coordinates and weights
// Gauss integration with power p=4 in the polynomial in xi[0], xi[1], xi[2].
int lg = 4; // number of points for the standard Gauss integration (formule de quadrature a' lg points)
double xi[3][4]; // i=1..dim, j=1..lg
double omega[4]; // i=1..lg

// some variables for the linear biconjugate gradient method
int itol;
int itmax;
int Quantum_Flag;
double tol;

double p[MAXNUMNODES+2];
double pp[MAXNUMNODES+2];
double r[MAXNUMNODES+2];
double rr[MAXNUMNODES+2];
double z[MAXNUMNODES+2];
double zz[MAXNUMNODES+2];

// time structures
time_t binarytime;
struct tm *nowtm;

time_t fbinarytime;
struct tm *fnowtm;

// MC variables
int EL;
int number_of_steps, it;
int NOVALLEY[NOAMTIA+1]; // see montecarlo_setup.h for more informations
int ISEED, NP1, INUM, IV;
int NOELEC[MAXNUMELEMENTS+1];
int MEAN;
//double **resr;
double  *da;  // k = 1./epsr*eps0
double *NA;
double *ND, DT, TS, *NE;
double KX, KY, KZ, TF;
double DDMAX, *EPP, L[4];
double P[NPMAX+1][10];
double BKTQ, QH, TL, CIMP, QD2;
double MSTAR[NOAMTIA+1][6];
double *VOLUME, TEMPO;
double SWK[NOAMTIA+1][3][15][DIME+1];
double AF[NOAMTIA+1][3],EC[NOAMTIA+1][3];
double HM[NOAMTIA+1][3],GM[NOAMTIA+1];
double SMH[NOAMTIA+1][3],HHM[NOAMTIA+1][3];

double *XVEL; // X-component of Velocity of ELectrons
double *YVEL; // Y-component of Velocity of ELectrons
double *ZVEL; // Z-component of Velocity of ELectrons
double *ENEL; // ENergy of ELectrons

int err;
int every_num_steps;
int NHTFLAG;
int POISSONUPDATEFLAG;
int SAVEPARTICLESFLAG;
int SAVEPARTICLESFORMAT;
int SAVEFIELDSFLAG;
int SAVEFIELDSFORMAT;

int NPMAX0;
int WINX,WINY;
double up,down;

// Impact Ionization flag
int IIFLAG;

char MESHNAME[256];
double XVAL[NOAMTIA+1];

//int RELATIVEDISTANCESFLAG;

// All files here...
FILE *fp;

struct option longopts[] =
{
  { "version", no_argument, NULL, 'v' },
  { "help", no_argument, NULL, 'h' }
};

// All strings here...
static char *progname;

//int the_main_single_step(void);
#include "headers/included_files.h"

// ********************************************************************
//                         M    A    I    N
// ********************************************************************

int main(int argc, char *argv[])
{
 int index;
 int optc;
 int h = 0, v = 0, lose = 0, z = 0;

  progname = argv[0];

  while ((optc = getopt_long (argc, argv, "hv", longopts, (int *) 0))
         != EOF)
    switch (optc)
      {
      case 'v':
        v = 1;
        break;
      case 'h':
        h = 1;
        break;
      default:
        lose = 1;
        break;
      }
  
  if (optind == argc - 1)
    z = 1;
  else if (lose || optind < argc)
    {
      /* Print error message and exit.  */
      if (optind < argc)
        printf("Too many arguments\n");
        printf("Try `%s --help' for more information.\n",progname);
      exit(1);
    }

  /* `help' should come first.  If `help' is requested, ignore the other
     options. */
  if (h)
    {
      /* Print help info and exit.  */
      /* TRANSLATORS: --help output 1
         no-wrap */
      printf("\
Aeneas, the GNU 3D simulator for III-V semiconductor devices.\n");
      printf ("\n");
      /* TRANSLATORS: --help output 2
         no-wrap */
      printf ("\
Usage: %s [OPTION] file...\n",progname);

      printf ("\n");
      /* TRANSLATORS: --help output 3 : options 1/2
         no-wrap */
      printf("\
  -h, --help          display this help and exit\n\
  -v, --version       display version information and exit\n");

      printf ("\n");
      /* TRANSLATORS: --help output 5 (end)
         TRANSLATORS, please don't forget to add the contact address for
         your translation!
         no-wrap */
      printf ("\
Report bugs to sellier@dmi.unict.it or aeneas@southnovel.eu\n");

      exit (0);
    }

  if (v)
    {
      /* Print version number.  */
      printf("aeneas - GNU aeneas 1.1\n");
      /* xgettext: no-wrap */
      printf("\n");
      printf("\
Copyright (C) %s Sellier Jean Michel.\n\
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A\n\
PARTICULAR PURPOSE.\n\
You may redistribute copies of GNU %s under the terms\n\
of the GNU General Public License.\n\
For more information about these matters, see the file named COPYING.\n",
              "2007","Aeneas");
      exit (0);
    }
  else if (z){
// In case of filename specified
// =============================
     fp=fopen(argv[1],"r"); // here we open th input file...
// File Control, just in case the file does not exist...
     if(fp==NULL){
      printf("%s: fatal error in opening the input file %s\n",
             progname,argv[1]);
      exit(EXIT_FAILURE);
     }

// =========================

// creation of the output directories
  printf("\ntrying to create the output directories.\n");
  err = system("mkdir e_density");
  err = system("mkdir e_density/BB_mesh_format");
  err = system("mkdir e_density/VTK_format");
  err = system("mkdir e_velocity");
  err = system("mkdir e_velocity/xvel_BB_mesh_format");
  err = system("mkdir e_velocity/yvel_BB_mesh_format");
  err = system("mkdir e_velocity/zvel_BB_mesh_format");
  err = system("mkdir e_velocity/VTK_format");
  err = system("mkdir e_energy");
  err = system("mkdir e_energy/BB_mesh_format");
  err = system("mkdir e_energy/VTK_format");
  err = system("mkdir E");
  err = system("mkdir E/Ex_BB_mesh_format");
  err = system("mkdir E/Ey_BB_mesh_format");
  err = system("mkdir E/Ez_BB_mesh_format");
  err = system("mkdir E/VTK_format");
  err = system("mkdir e_current");
  err = system("mkdir e_current/e_current_x_BB_mesh_format");
  err = system("mkdir e_current/e_current_y_BB_mesh_format");
  err = system("mkdir e_current/e_current_z_BB_mesh_format");
  err = system("mkdir e_current/VTK_format");
  err = system("mkdir potential");
  err = system("mkdir potential/BB_mesh_format");
  err = system("mkdir potential/VTK_format");
  err = system("mkdir particles");
  err = system("mkdir particles/position_mesh_format");
  err = system("mkdir particles/position_Point3D_format");
  err = system("mkdir particles/energy_Point3D_format");
  err = system("mkdir particles/potential_Point3D_format");

// Printing of the initial runtime
// ===============================
  binarytime=time(NULL);
  nowtm=localtime(&binarytime);
  printf("\n\nSimulation started : %s\n",asctime(nowtm));


//  initial value for random number generator
  ISEED = 38467.;

// Open the input file and check if it exists
// If the input file does not exist then the code stops.
//  printf("\n\nLoading the input mesh file...\n");
//  if(1) load_mesh(); // this row has to be putted before everything calculation!!!
//  if(0) load_mshV1(); // this subroutine is recommended only for GMSH experts
// warning : the user has to be very experienced in GMSH

// ***************************
// HERE WE READ THE INPUT FILE
// ***************************
  read_input_file();

// closure of the user specified input file
  fclose(fp);

// If no error has occured then print a nice message :)
  printf("Mesh input file correctly processed...\n\n");

// ==============================
// ==============================
// ========================
// Material constants here!
// ========================
// Silicon electrons in the X-valley
// Germanium electrons in the X-valley
// All the other semiconductor compounds have electrons in the Gamma- and L-valley
// definition of the number of valley for every material
 NOVALLEY[SILICON]=1;
 NOVALLEY[GAAS]=2;
 NOVALLEY[GERMANIUM]=1;
 NOVALLEY[INSB]=2;
 NOVALLEY[ALSB]=2;
 NOVALLEY[ALAS]=2;
 NOVALLEY[ALP]=2;
 NOVALLEY[GAP]=2;
 NOVALLEY[GASB]=2;
 NOVALLEY[INAS]=2;
 NOVALLEY[INP]=2;
 NOVALLEY[ALXINXSB]=2;
 NOVALLEY[ALXIN1XSB]=2;
// Relative dielectric constant for Silicon
 EPSR[SILICON]=11.7;
// Relative dielectric constant for Germanium
 EPSR[GERMANIUM]=16.2;
// Relative dielectric constant for InSb
 EPSR[INSB]=17.65;
// Relative dieletric constant for GaAs
 EPSR[GAAS]=12.90;
// Relative dieletric constant for AlSb
 EPSR[ALSB]=12.04;
// Relative dieletric constant for AlAs
 EPSR[ALAS]=10.06;
// Relative dieletric constant for AlP
 EPSR[ALP]=9.80;
// Relative dieletric constant for GaP
 EPSR[GAP]=11.10;
// Relative dieletric constant for GaSb
 EPSR[GASB]=15.69;
// Relative dieletric constant for InAs
 EPSR[INAS]=15.15;
// Relative dieletric constant for InP
 EPSR[INP]=12.61;
// Relative dieletric constant for Al_x In_x Sb
 EPSR[ALXINXSB]=XVAL[ALXINXSB]*EPSR[ALSB]+XVAL[ALXINXSB]*EPSR[INSB];
// Relative dieletric constant for Al_x In_(1-x) Sb
 EPSR[ALXIN1XSB]=XVAL[ALXIN1XSB]*EPSR[ALSB]+(1.-XVAL[ALXIN1XSB])*EPSR[INSB];
// InSb high frequency dieletric constant
 EPF[INSB]=15.68;
// GaAs high frequency dieletric constant
 EPF[GAAS]=10.92;
// AlSb high frequency dieletric constant
 EPF[ALSB]=9.88;
// AlAs high frequency dieletric constant
 EPF[ALAS]=8.16;
// AlP high frequency dieletric constant
 EPF[ALP]=7.54;
// GaP high frequency dieletric constant
 EPF[GAP]=9.08;
// GaSb high frequency dieletric constant
 EPF[GASB]=14.44;
// InAs high frequency dieletric constant
 EPF[INAS]=12.75;
// InP high frequency dieletric constant
 EPF[INP]=9.61;
// Al_x In_x Sb high frequency dieletric constant
 EPF[ALXINXSB]=XVAL[ALXINXSB]*(EPF[ALSB]+EPF[INSB]);
// Al_x In_(1-x) Sb high frequency dieletric constant
 EPF[ALXIN1XSB]=XVAL[ALXIN1XSB]*EPF[ALSB]+(1.-XVAL[ALXIN1XSB])*EPF[INSB];
 for(index=0;index<NOAMTIA;index++){
   int i;
   for(i=0;i<6;i++){
     HWO[index][i]=0.;
     DTK[index][i]=0.;
     ZF[index][i]=0.;
   }
 }
// Silicon optical phonon scattering energy (eV)
 HWO[SILICON][0]=0.0120;
 HWO[SILICON][1]=0.0185;
 HWO[SILICON][2]=0.0190;
 HWO[SILICON][3]=0.0474;
 HWO[SILICON][4]=0.0612;
 HWO[SILICON][5]=0.0590;
// Germanium optical phonon scattering energy (eV)
 HWO[GERMANIUM][0]=0.09;
// InSb optical phonon scattering energy (eV)
 HWO[INSB][0]=0.0282;
// GaAs optical phonon scattering energy (eV)
 HWO[GAAS][0]=0.03536;
// AlSb optical phonon scattering energy (eV)
 HWO[ALSB][0]=0.0360;
// AlAs optical phonon scattering energy (eV)
 HWO[ALAS][0]=0.05009;
// AlP optical phonon scattering energy (eV)
 HWO[ALP][0]=0.06212;
// GaP optical phonon scattering energy (eV)
 HWO[GAP][0]=0.04523;
// GaSb optical phonon scattering energy (eV)
 HWO[GASB][0]=0.02529;
// InAs optical phonon scattering energy (eV)
 HWO[INAS][0]=0.03008;
// InP optical phonon scattering energy (eV)
 HWO[INP][0]=0.04240;
// Al_x In_x Sb optical phonon scattering energy (eV)
 HWO[ALXINXSB][0]=XVAL[ALXINXSB]*(HWO[ALSB][0]+HWO[INSB][0]);
// Al_x In_(1-x) Sb optical phonon scattering energy (eV)
 HWO[ALXIN1XSB][0]=XVAL[ALXIN1XSB]*HWO[ALSB][0]+(1.-XVAL[ALXIN1XSB])*HWO[INSB][0];
// Silicon optical coupling constants (eV/m)
 DTK[SILICON][0]=0.05e11;
 DTK[SILICON][1]=0.08e11;
 DTK[SILICON][2]=0.03e11;
 DTK[SILICON][3]=0.20e11;
 DTK[SILICON][4]=1.14e11;
 DTK[SILICON][5]=0.20e11;
// Germanium optical coupling constants (eV/m)
 DTK[GERMANIUM][0]=1.889e11;
// InSb optical coupling constants (eV/m)
 DTK[INSB][0]=0.47e11;
// GaAs optical coupling constant (eV/m)
 DTK[GAAS][0]=1.11e11;
// AlSb optical coupling constants (eV/m)
 DTK[ALSB][0]=0.55e11;
// AlAs optical coupling constants (eV/m)
 DTK[ALAS][0]=3.0e11;
// AlP optical coupling constants (eV/m)
 DTK[ALP][0]=0.95e11;
// GaP optical coupling constants (eV/m)
 DTK[GAP][0]=5.33e11;
// GaSb optical coupling constants (eV/m)
 DTK[GASB][0]=0.94e11;
// InAs optical coupling constants (eV/m)
 DTK[INAS][0]=3.59e11;
// InP optical coupling constants (eV/m)
 DTK[INP][0]=2.46e11;
// Al_x In_x Sb optical coupling constants (eV/m)
 DTK[ALXINXSB][0]=XVAL[ALXINXSB]*(DTK[ALSB][0]+DTK[INSB][0]);
// Al_x In_(1-x) Sb optical coupling constants (eV/m)
 DTK[ALXIN1XSB][0]=XVAL[ALXIN1XSB]*DTK[ALSB][0]+(1.-XVAL[ALXIN1XSB])*DTK[INSB][0];
// Silicon optical phonon Z-factor
 ZF[SILICON][0]=1.;
 ZF[SILICON][1]=1.;
 ZF[SILICON][2]=4.;
 ZF[SILICON][3]=4.;
 ZF[SILICON][4]=1.;
 ZF[SILICON][5]=4.;
// Germanium optical phonon Z-factor
 ZF[GERMANIUM][0]=1.;
// InSb optical phonon Z-factor
 ZF[INSB][0]=1.;
// AlSb optical phonon Z-factor
 ZF[ALSB][0]=1.;
// optical phonon Z-factor
 ZF[ALAS][0]=1.;
// optical phonon Z-factor
 ZF[ALP][0]=1.;
// optical phonon Z-factor
 ZF[GAP][0]=1.;
// optical phonon Z-factor
 ZF[GASB][0]=1.;
// optical phonon Z-factor
 ZF[INAS][0]=1.;
// optical phonon Z-factor
 ZF[INP][0]=1.;
// Al_x In_x Sb optical phonon Z-factor
 ZF[ALXINXSB][0]=XVAL[ALXINXSB]*(ZF[ALSB][0]+ZF[INSB][0]);
// Al_x In_(1-x) Sb optical phonon Z-factor
 ZF[ALXIN1XSB][0]=XVAL[ALXIN1XSB]*ZF[ALSB][0]+(1.-XVAL[ALXIN1XSB])*ZF[INSB][0];
// Silicon Crystal Density (Kg/m^3)
 RHO[SILICON]=2330.;
// Germanium Crystal Density (Kg/m^3)
 RHO[GERMANIUM]=5320.;
// InSb Crystal Density (Kg/m^3)
 RHO[INSB]=5790.;
// GaAs Crystal Density (Kg/m^3)
 RHO[GAAS]=5360.;
// AlSb Crystal Density (Kg/m^3)
 RHO[ALSB]=4288.;
// AlAs Crystal Density (Kg/m^3)
 RHO[ALAS]=5650.;
// AlP Crystal Density (Kg/m^3)
 RHO[ALP]=7410.;
// GaP Crystal Density (Kg/m^3)
 RHO[GAP]=5850.;
// GaSb Crystal Density (Kg/m^3)
 RHO[GASB]=3970.;
// InAs Crystal Density (Kg/m^3)
 RHO[INAS]=4280.;
// InP Crystal Density (Kg/m^3)
 RHO[INP]=5130.;
// Al_x In_x Sb Crystal Density (Kg/m^3)
 RHO[ALXINXSB]=XVAL[ALXINXSB]*(RHO[ALSB]+RHO[INSB]); 
// Al_x In_(1-x) Sb Crystal Density (Kg/m^3)
 RHO[ALXIN1XSB]=XVAL[ALXIN1XSB]*RHO[ALSB]+(1.-XVAL[ALXIN1XSB])*RHO[INSB];
// Silicon acoustic deformation potential (Joule)
 DA[SILICON]=9.*1.60217733e-19;
// Germanium acoustic deformation potential (Joule)
 DA[GERMANIUM]=9.*1.60217733e-19;
// InSb acoustic deformation potential (Joule)
 DA[INSB]=5.96*1.60217733e-19;
// GaAs acoustic deformation potential (Joule)
 DA[GAAS]=7.*1.60217733e-19;
// AlSb acoustic deformation potential (Joule)
 DA[ALSB]=4.6*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[ALAS]=9.3*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[ALP]=9.3*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[GAP]=7.4*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[GASB]=9.*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[INAS]=8.2*1.60217733e-19;
// acoustic deformation potential (Joule)
 DA[INP]=6.2*1.60217733e-19;
// Al_x In_x Sb acoustic deformation potential (Joule)
 DA[ALXINXSB]=XVAL[ALXINXSB]*(DA[ALSB]+DA[INSB]);
// Al_x In_(1-x) Sb acoustic deformation potential (Joule)
 DA[ALXIN1XSB]=XVAL[ALXIN1XSB]*DA[ALSB]+(1.-XVAL[ALXIN1XSB])*DA[INSB];
// Silicon longitudinal sound velocity (m/sec)
 UL[SILICON]=9040.;
// Germanium longitudinal sound velocity (m/sec)
 UL[GERMANIUM]=5400.;
// InSb longitudinal sound velocity (m/sec)
 UL[INSB]=4060.;
// GaAs longitudinal sound velocity (m/sec)
 UL[GAAS]=5240.;
// AlSb longitudinal sound velocity (m/sec)
 UL[ALSB]=4600.;
// longitudinal sound velocity (m/sec)
 UL[ALAS]=5650.;
// longitudinal sound velocity (m/sec)
 UL[ALP]=7410.;
// longitudinal sound velocity (m/sec)
 UL[GAP]=5850.;
// longitudinal sound velocity (m/sec)
 UL[GASB]=3970.;
// longitudinal sound velocity (m/sec)
 UL[INAS]=4280.;
// longitudinal sound velocity (m/sec)
 UL[INP]=5130.;
// Al_x In_x Sb longitudinal sound velocity (m/sec)
 UL[ALXINXSB]=XVAL[ALXINXSB]*(UL[ALSB]+UL[INSB]);
// Al_x In_(1-x) Sb longitudinal sound velocity (m/sec)
 UL[ALXIN1XSB]=XVAL[ALXIN1XSB]*UL[ALSB]+(1.-XVAL[ALXIN1XSB])*UL[INSB];
// Silicon energy gap
 EG[SILICON]=1.21-3.333e-4*TL;
 printf("EG[SILICON]      = %g\n",EG[SILICON]);
// Germanium energy gap
 EG[GERMANIUM]=0.747-3.587e-4*TL;
 printf("EG[GERMANIUM]    = %g\n",EG[GERMANIUM]);
// InSb energy gap
 EG[INSB]=0.2446-2.153e-4*TL;//0.174;???
 printf("EG[INSB]         = %g\n",EG[INSB]);
// GaAs energy gap
 EG[GAAS]=1.54-4.036e-4*TL;//0.32;???
 printf("EG[GAAS]         = %g\n",EG[GAAS]);
// AlSb energy gap
 EG[ALSB]=1.696-2.20e-4*TL;
 printf("EG[ALSB]         = %g\n",EG[ALSB]);
// energy gap
 EG[ALAS]=2.314-3.0e-4*TL;
// energy gap
 EG[ALP]=2.51-3.333e-4*TL;
// energy gap
 EG[GAP]=2.35-2.667e-4*TL;
// energy gap
 EG[GASB]=0.81-3.667e-4*TL;
// energy gap
 EG[INAS]=0.434-2.601e-4*TL;
// energy gap
 EG[INP]=1.445-3.296e-4*TL;
// Al_x In_x Sb energy gap
 EG[ALXINXSB]=XVAL[ALXINXSB]*(EG[ALSB]+EG[INSB]);
 printf("EG[ALxINxSB]     = %g\n",EG[ALXINXSB]);
// Al_x In_(1-x) Sb energy gap
 EG[ALXIN1XSB]=XVAL[ALXIN1XSB]*EG[ALSB]+(1.-XVAL[ALXIN1XSB])*EG[INSB];
 printf("EG[ALxIN(1-x)SB] = %g\n",EG[ALXIN1XSB]);
 printf("\n");
// Silicon energy minimum
 EMIN[SILICON][1]=0.0;
// Germanium energy minimum
 EMIN[GERMANIUM][1]=0.173;
// InSb energy minimum of GAMMA-valley
 EMIN[INSB][1]=0.0;
// InSb energy minimum 0f L-valley
 EMIN[INSB][2]=1.038;
// GaAs energy minimum of GAMMA-valley
 EMIN[GAAS][1]=0.0;
// GaAs energy minimum of L-valley (eV)
 EMIN[GAAS][2]=0.323;
// AlSb energy minimum of GAMMA-valley
 EMIN[ALSB][1]=0.507;
// AlSb energy minimum 0f L-valley
 EMIN[ALSB][2]=0.292;
// energy minimum of GAMMA-valley
 EMIN[ALAS][1]=0.767;
// energy minimum 0f L-valley
 EMIN[ALAS][2]=0.332;
// energy minimum of GAMMA-valley
 EMIN[ALP][1]=1.237;
// energy minimum 0f L-valley
 EMIN[ALP][2]=0.824;
// energy minimum of GAMMA-valley
 EMIN[GAP][1]=0.496;
// energy minimum 0f L-valley
 EMIN[GAP][2]=0.415;
// energy minimum of GAMMA-valley
 EMIN[GASB][1]=0.0;
// energy minimum 0f L-valley
 EMIN[GASB][2]=0.217;
// energy minimum of GAMMA-valley
 EMIN[INAS][1]=0.0;
// energy minimum 0f L-valley
 EMIN[INAS][2]=1.078;
// energy minimum of GAMMA-valley
 EMIN[INP][1]=0.0;
// energy minimum 0f L-valley
 EMIN[INP][2]=0.832;
// Al_x In_x Sb energy minimum of GAMMA-valley
 EMIN[ALXINXSB][1]=XVAL[ALXINXSB]*(EMIN[ALSB][1]+EMIN[INSB][1]);
// Al_x In_x Sb energy minimum 0f L-valley
 EMIN[ALXINXSB][2]=XVAL[ALXINXSB]*(EMIN[ALSB][2]+EMIN[INSB][2]);
// Al_x In_(1-x) Sb energy minimum of GAMMA-valley
 EMIN[ALXIN1XSB][1]=XVAL[ALXIN1XSB]*EMIN[ALSB][1]+(1.-XVAL[ALXIN1XSB])*EMIN[INSB][1];
// Al_x In_(1-x) Sb energy minimum 0f L-valley
 EMIN[ALXIN1XSB][2]=XVAL[ALXIN1XSB]*EMIN[ALSB][2]+(1.-XVAL[ALXIN1XSB])*EMIN[INSB][2];
// Definition of effective mass for all materials in all valleys
 MSTAR[SILICON][1]=0.32;
 MSTAR[GAAS][1]=0.067; // Gamma-valley
 MSTAR[GAAS][2]=0.350; // L-valley
 MSTAR[GERMANIUM][1]=0.22;
 MSTAR[INSB][1]=0.013;
 MSTAR[INSB][2]=0.442;
 MSTAR[ALSB][1]=0.14;
 MSTAR[ALSB][2]=0.442;
 MSTAR[ALAS][1]=0.149;
 MSTAR[ALAS][2]=pow(1.456*pow(0.165,2.),1./3.);
 MSTAR[ALP][1]=0.166;
 MSTAR[ALP][2]=pow(1.556*pow(0.162,2.),1./3.);
 MSTAR[GAP][1]=0.122;
 MSTAR[GAP][2]=pow(1.498*pow(0.358,2.),1./3.);
 MSTAR[GASB][1]=0.048;
 MSTAR[GASB][2]=pow(1.4*pow(0.139,2.),1./3.);
 MSTAR[INAS][1]=0.031;
 MSTAR[INAS][2]=pow(1.590*pow(0.326,2.),1./3.);
 MSTAR[INP][1]=0.083;
 MSTAR[INP][2]=pow(1.893*pow(0.425,2.),1./3.);
 MSTAR[ALXINXSB][1]=XVAL[ALXINXSB]*(MSTAR[ALSB][1]+MSTAR[INSB][1]);
 MSTAR[ALXINXSB][2]=XVAL[ALXINXSB]*(MSTAR[ALSB][2]+MSTAR[INSB][2]);
 MSTAR[ALXIN1XSB][1]=XVAL[ALXIN1XSB]*MSTAR[ALSB][1]+(1.-XVAL[ALXIN1XSB])*MSTAR[INSB][1];
 MSTAR[ALXIN1XSB][2]=XVAL[ALXIN1XSB]*MSTAR[ALSB][2]+(1.-XVAL[ALXIN1XSB])*MSTAR[INSB][2];
// non-parabolicity coefficient for Silicon
 alphaK[SILICON][1]=0.5;
// non-parabolicity coefficient for Germanium
 alphaK[GERMANIUM][1]=0.65;
// non-parabolicity coefficient for InSb in the GAMMA-valley
 alphaK[INSB][1]=pow(1.-MSTAR[INSB][1],2.)/EG[INSB];//5.59;
 printf("alphaK_gamma[InSb] = %g\n",alphaK[INSB][1]); 
// non-parabolicity coefficient for InSb in the L-valley
 alphaK[INSB][2]=pow(1.-MSTAR[INSB][2],2.)/EG[INSB];//1.789;
 printf("alphaK_L[InSb]     = %g\n",alphaK[INSB][2]);
// non-parabolicity coefficient for GaAs in the GAMMA-valley
 alphaK[GAAS][1]=pow(1.-MSTAR[GAAS][1],2.)/EG[GAAS];//0.611;
 printf("alphaK_gamma[GaAs] = %g\n",alphaK[GAAS][1]);
// non-parabolicity coefficient for GaAs in the L-valley
 alphaK[GAAS][2]=pow(1.-MSTAR[GAAS][2],2.)/EG[GAAS];//0.242;
 printf("alphaK_L[GaAs]     = %g\n",alphaK[GAAS][2]);
// non-parabolicity coefficient for AlSb in the GAMMA-valley
 alphaK[ALSB][1]=pow(1.-MSTAR[ALSB][1],2.)/EG[ALSB];//0.321;
 printf("alphaK_gamma[AlSb] = %g\n",alphaK[ALSB][1]);
// non-parabolicity coefficient for AlSb in the L-valley
 alphaK[ALSB][2]=pow(1.-MSTAR[ALSB][2],2.)/EG[ALSB];//0.135;
 printf("alphaK_L[AlSb]     = %g\n",alphaK[ALSB][2]);
 for(index=ALAS;index<=INP;index++){
// non-parabolicity coefficient in the GAMMA-valley
 alphaK[index][1]=pow(1.-MSTAR[index][1],2.)/EG[index];
 printf("alphaK_gamma[%d] = %g\n",index,alphaK[index][1]);
// non-parabolicity coefficient in the L-valley
 alphaK[index][2]=pow(1.-MSTAR[index][2],2.)/EG[index];
 printf("alphaK_L[%d]     = %g\n",index,alphaK[index][2]);
 }
// non-parabolicity coefficient for Al_x In_x Sb in the GAMMA-valley
 alphaK[ALXINXSB][1]=XVAL[ALXINXSB]*(alphaK[ALSB][1]+alphaK[INSB][1]);
// non-parabolicity coefficient for Al_x In_x Sb in the L-valley
 alphaK[ALXINXSB][2]=XVAL[ALXINXSB]*(alphaK[ALSB][2]+alphaK[INSB][2]);
// non-parabolicity coefficient for Al_x In_(1-x) Sb in the GAMMA-valley
 alphaK[ALXIN1XSB][1]=XVAL[ALXIN1XSB]*alphaK[ALSB][1]+(1.-XVAL[ALXIN1XSB])*alphaK[INSB][1];
// non-parabolicity coefficient for Al_x In_(1-x) Sb in the L-valley
 alphaK[ALXIN1XSB][2]=XVAL[ALXIN1XSB]*alphaK[ALSB][2]+(1.-XVAL[ALXIN1XSB])*alphaK[INSB][2];
 printf("\n");
// constant for the Impact-Ionization Keldysh model (1/sec)
 PII[SILICON]=6.80e14;
 PII[GERMANIUM]=2.0e12;
 PII[INSB]=2.0e12;//7.5e14;
 PII[GAAS]=2.5e15;
 PII[ALSB]=1.5e15;
 PII[ALAS]=1.5e15;
 PII[ALP]=1.5e15;
 PII[GAP]=2.5e14;
 PII[GASB]=1.5e15;
 PII[INAS]=7.5e14;
 PII[INP]=7.5e13;
 PII[ALXINXSB]=XVAL[ALXINXSB]*(PII[ALSB]+PII[INSB]);
 PII[ALXIN1XSB]=XVAL[ALXIN1XSB]*PII[ALSB]+(1.-XVAL[ALXIN1XSB])*PII[INSB];
// Ionization threshold
 ETH[SILICON]=3.45;
 ETH[GERMANIUM]=EG[GERMANIUM];
 ETH[INSB]=0.036+EG[INSB];
 ETH[GAAS]=0.30+EG[GAAS];
 ETH[ALSB]=0.30+EG[ALSB];
 ETH[ALAS]=0.30+EG[ALAS];
 ETH[ALP]=0.30+EG[ALP];
 ETH[GAP]=0.10+EG[GAP];
 ETH[GASB]=0.30+EG[GASB];
 ETH[INAS]=0.03+EG[INAS];
 ETH[INP]=0.50+EG[INP];
 ETH[ALXINXSB]=XVAL[ALXINXSB]*(ETH[ALSB]+ETH[INSB]);
 ETH[ALXIN1XSB]=XVAL[ALXIN1XSB]*ETH[ALSB]+(1.-XVAL[ALXIN1XSB])*ETH[INSB];
// ==============================
// ==============================

// definition of the neighbourhood table
  if(NHTFLAG==WRITE){
   printf("Computing the neighbourhood table...\n");
   neighbourhood_table();
   printf("Neighbourhood table calculated and saved.\n\n");
  }

// reading the neighbourhood table
  if(NHTFLAG==READ){
   printf("Reading the neighbourhood table...\n");
   read_neighbourhood_table();
   printf("Neighbourhood table read.\n\n");
  }

// various dynamical allocations
// #include "headers/services/dynamical_allocations.h"

// various parameters setup for montecarlo
  printf("\n\nSetting up the materials table...\n");
  for(index=0;index<NOAMTIA;index++) montecarlo_setup(index);
  printf("Materials table setted.\n");

// simulated device configuration
  printf("\n\nSetting up the device...\n");
  device_setup();
  printf("Device setted.\n");

// save_particles_mesh_format(0);
// exit(0);

// Initialise the arrays xi[][] and omega[]
  initialize_GAUSS();

// Definition of all the **base_geo and **base_ref elements.
  def_base_all();

// Definition of all the **d_base_geo and **d_base_ref elements.
  def_d_base_all();

// Definition of all the **poids_K elements.
  def_global_weights_all();

// definition of all components of ***inv_jac_K
  define_inverse_jacobian();

// definition of the arrays ****d_base_K
  define_d_base_K();

// every boundary which has not been specified
// is considered to be of insulator type as default
{
 int i;
 int ij[4];
 int j,iflag;
 for(i=0;i<Ng;i++)
  if(i_front[i]!=OHMIC && i_front[i]!=SCHOTTKY) i_front[i]=4;
 for(i=0;i<Ne;i++){
   iflag=0;
//   printf("%d of %d\n",i,Ne);
   if(vois[0][i]==-1){ij[1]=1; ij[2]=2; ij[3]=3; iflag=1;}
   if(vois[1][i]==-1){ij[1]=0; ij[2]=2; ij[3]=3; iflag=1;}
   if(vois[2][i]==-1){ij[1]=1; ij[2]=0; ij[3]=3; iflag=1;}
   if(vois[3][i]==-1){ij[1]=1; ij[2]=2; ij[3]=0; iflag=1;}
   if(iflag==1) for(j=1;j<=3;j++){
//    printf("j=%d %d\n",j,iij[j]);
    if(i_front[noeud_geo[ij[j]][i]-1]!=OHMIC &&
       i_front[noeud_geo[ij[j]][i]-1]!=SCHOTTKY)
       i_front[noeud_geo[ij[j]][i]-1]=INSULATOR;
   }
 }
}

// definition of the A matrix of the 3D Poisson related system
// in a row-indexed sparse storage mode
  assemble_matrix_A();

// definition of the b rhs of the 3D Poisson related system
  assemble_rhs_B();
// {int i; for(i=1;i<=Ng;i++) printf("main b[%d]=%g\n",i,b[i]);}

// resolution of 3D Poisson equation (finite element method)
// (only for the applied electrostatic potential!)
  poisson3D(Vi,itol,tol,itmax);
  {int j; <