Source Jam - aeneas-1.1

00001 /*
00002    This file belongs to Aeneas. Aeneas is a GNU package released under GPL 3 license.
00003    This code is a simulator for Submicron 3D Semiconductor Devices. 
00004    It implements the Monte Carlo transport in 3D tetrahedra meshes
00005    for the simulation of the semiclassical Boltzmann equation for both electrons.
00006    It also includes all the relevant quantum effects for nanodevices.
00007 
00008    Copyright (C) 2007 Jean Michel Sellier <sellier@dmi.unict.it>
00009  
00010    This program is free software; you can redistribute it and/or modify
00011    it under the terms of the GNU General Public License as published by
00012    the Free Software Foundation; either version 3, or (at your option)
00013    any later version.
00014 
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019 
00020    You should have received a copy of the GNU General Public License
00021    along with this program. If not, see <http://www.gnu.org/licenses/>.
00022 */
00023 
00024 // Created on : 11 june 2007, Siracusa, Jean Michel Sellier
00025 // Last modified : 16 august 2007, Siracusa, Jean Michel Sellier
00026 
00027 // drift step of the EMC algorithm
00028 
00029 void drift(double tau)
00030 {
00031  register int j;
00032  int n=4;
00033  int p=0;
00034  int ifr[3]={0.,0.,0.};
00035  int indx[5];
00036  double d;
00037  double a[5][5],b[5];
00038  double dkx,dky,dkz;
00039  double hmt, x, y, z;
00040  double ksquared,gk,thesquareroot;
00041 
00042  if(IV==9) return;
00043 // if(EL==-1) return; // just in case...
00044  if(i_dom[EL]==SIO2) return;
00045 
00046 // Valid for every material
00047  dkx=-QH*0.25*(Ex[noeud_geo[0][EL]-1]+Ex[noeud_geo[1][EL]-1]
00048               +Ex[noeud_geo[2][EL]-1]+Ex[noeud_geo[3][EL]-1])*tau;
00049  dky=-QH*0.25*(Ey[noeud_geo[0][EL]-1]+Ey[noeud_geo[1][EL]-1]
00050               +Ey[noeud_geo[2][EL]-1]+Ey[noeud_geo[3][EL]-1])*tau;
00051  dkz=-QH*0.25*(Ez[noeud_geo[0][EL]-1]+Ez[noeud_geo[1][EL]-1]
00052               +Ez[noeud_geo[2][EL]-1]+Ez[noeud_geo[3][EL]-1])*tau;
00053 // printf("KX=%g KY=%g KZ=%g\n",KX,KY,KZ);
00054 // printf("dkx=%g dky=%g dkz=%g\n",dkx,dky,dkz);
00055  
00056  x=y=z=0.;
00057  for(j=0;j<4;j++){
00058   x+=L[j]*coord[0][noeud_geo[j][EL]-1];
00059   y+=L[j]*coord[1][noeud_geo[j][EL]-1];
00060   z+=L[j]*coord[2][noeud_geo[j][EL]-1];
00061  }
00062  x/=4.; y/=4.; z/=4.;
00063 //if(z<-1.){
00064 //   printf("%g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]); 
00065 //   for(j=0;j<4;j++) printf("%d %g %g %g\n",j,coord[0][noeud_geo[j][EL]-1],coord[1][noeud_geo[j][EL]-1],coord[2][noeud_geo[j][EL]-1]);
00066 //   printf("%d %d x = %g y = %g z = %g\n",EL,noeud_geo[3][EL]-1,x,y,z);
00067 //}
00068 
00069 // x*=1.e-6; y*=1.e-6; z*=1.e-6;
00070   
00071  if(i_dom[EL]!=SIO2){
00072 // Electron drift process
00073    hmt=HM[i_dom[EL]][IV]*tau;
00074    ksquared=KX*KX+KY*KY+KZ*KZ;
00075    gk=HHM[i_dom[EL]][IV]*ksquared;
00076    thesquareroot=sqrt(1.+4.*alphaK[i_dom[EL]][IV]*gk);
00077 //   printf("!!! %g\n",thesquareroot);
00078  }
00079 // updating the position and pseudo-wave vector
00080 // (second order Runge-Kutta method)
00081 //printf("hmt = %g\n",hmt);
00082 // printf("X=%g Y=%g Z=%g\n",x,y,z);
00083 // printf("KX=%g KY=%g KZ=%g\n",KX,KY,KZ);
00084 // printf("dkx=%g dky=%g dkz=%g\n",dkx,dky,dkz);
00085 // printf("hmt = %g sq= %g\n",hmt,thesquareroot);
00086 // printf("%g %g %g\n",hmt*(KX+0.5*dkx)/thesquareroot,hmt*(KY+0.5*dky)/thesquareroot,
00087 //                     hmt*(KZ+0.5*dkz)/thesquareroot);
00088    x+=hmt*(KX+0.5*dkx)/thesquareroot;
00089    y+=hmt*(KY+0.5*dky)/thesquareroot;
00090    z+=hmt*(KZ+0.5*dkz)/thesquareroot;
00091    KX+=dkx;
00092    KY+=dky;
00093    KZ+=dkz;
00094 // printf("X=%g Y=%g Z=%g\n\n",x,y,z);
00095 
00096 // updating the barycentric coordinates
00097 //   x*=1.e6; y*=1.e6; z*=1.e6;
00098 //if(z>1.5) printf("x = %g y = %g z = %g\n",x,y,z);
00099 //   memset(&a,0,sizeof(a));
00100 //   memset(&b,0,sizeof(b));
00101    {int k;
00102     for(k=0;k<4;k++){
00103       a[1][k+1]=coord[0][noeud_geo[k][EL]-1];
00104       a[2][k+1]=coord[1][noeud_geo[k][EL]-1];
00105       a[3][k+1]=coord[2][noeud_geo[k][EL]-1];
00106       a[4][k+1]=1.;
00107    }}
00108    b[1]=4.*x;
00109    b[2]=4.*y;
00110    b[3]=4.*z;
00111    b[4]=4.;
00112 //printf("%g %g %g\n",coord[0][noeud_geo[1][EL]-1],coord[1][noeud_geo[1][EL]-1],
00113 //                    coord[2][noeud_geo[1][EL]-1]);
00114 //printf("%g %g %g\n",x,y,z);
00115 //printf("EL = %d\n",EL);
00116 //     {int ii,jj;
00117 //     for(ii=1;ii<=4;ii++){
00118 //        for(jj=1;jj<=4;jj++) printf("a[%d][%d] = %g   ",ii,jj,a[ii][jj]);
00119 //        printf("\n");
00120 //     }}
00121 // LU decomposition solver A.x=b
00122 #include "../linear_solver/ludcmp.h"
00123 #include "../linear_solver/lubksb.h"
00124 //   printf("before %g %g %g %g\n",L[0],L[1],L[2],L[3]);
00125    for(j=0;j<4;j++) L[j]=b[j+1];
00126 //   printf("** %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00127    if(L[0]<0. || L[1]<0. || L[2]<0. || L[3]<0.) goto label1;
00128    else return;
00129    label1:
00130 //    printf("EL = %d %g %g %g %g\n",EL,L[0],L[1],L[2],L[3]);
00131 //   printf("star %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00132     for(j=0;j<4;j++) if(L[j]<0.) p=j;
00133 // on the boundaries
00134 // =================
00135     if(p==0){
00136       ifr[0]=i_front[noeud_geo[1][EL]-1];
00137       ifr[1]=i_front[noeud_geo[2][EL]-1];
00138       ifr[2]=i_front[noeud_geo[3][EL]-1];
00139     }
00140     if(p==1){
00141       ifr[0]=i_front[noeud_geo[0][EL]-1];
00142       ifr[1]=i_front[noeud_geo[2][EL]-1];
00143       ifr[2]=i_front[noeud_geo[3][EL]-1];
00144     }
00145     if(p==2){
00146       ifr[0]=i_front[noeud_geo[0][EL]-1];
00147       ifr[1]=i_front[noeud_geo[1][EL]-1];
00148       ifr[2]=i_front[noeud_geo[3][EL]-1];
00149     }
00150     if(p==3){
00151       ifr[0]=i_front[noeud_geo[0][EL]-1];
00152       ifr[1]=i_front[noeud_geo[1][EL]-1];
00153       ifr[2]=i_front[noeud_geo[2][EL]-1];
00154     }
00155 // Schottky contact
00156     if(vois[p][EL]==-1 &&
00157        (ifr[0]==SCHOTTKY &&
00158         ifr[1]==SCHOTTKY &&
00159         ifr[2]==SCHOTTKY)){
00160             IV=9;
00161             return;
00162         }
00163 // ohmic contacts
00164     if(vois[p][EL]==-1 &&
00165        (ifr[0]==OHMIC &&
00166         ifr[1]==OHMIC &&
00167         ifr[2]==OHMIC)){
00168             IV=9;
00169             return;
00170         }
00171 // insulator contacts
00172     if(vois[p][EL]==-1 || i_dom[vois[p][EL]]==SIO2/* &&
00173        (ifr[0]==INSULATOR ||
00174         ifr[1]==INSULATOR ||
00175         ifr[2]==INSULATOR)*/){
00176         reflection:;
00177         double max=-4.;
00178         double x[2],y[2],z[2];
00179 // initial position
00180         x[0]=(L[0]*coord[0][noeud_geo[0][EL]-1]+
00181               L[1]*coord[0][noeud_geo[1][EL]-1]+
00182               L[2]*coord[0][noeud_geo[2][EL]-1]+
00183               L[3]*coord[0][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00184         y[0]=(L[0]*coord[1][noeud_geo[0][EL]-1]+
00185               L[1]*coord[1][noeud_geo[1][EL]-1]+
00186               L[2]*coord[1][noeud_geo[2][EL]-1]+
00187               L[3]*coord[1][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00188         z[0]=(L[0]*coord[2][noeud_geo[0][EL]-1]+
00189               L[1]*coord[2][noeud_geo[1][EL]-1]+
00190               L[2]*coord[2][noeud_geo[2][EL]-1]+
00191               L[3]*coord[2][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00192         if(L[0]>L[1] && L[0]>L[2]) max=L[0];
00193         if(L[1]>L[0] && L[1]>L[2]) max=L[1];
00194         if(L[2]>L[1] && L[2]>L[0]) max=L[2];
00195         for(j=0;j<3;j++){
00196           L[j]/=max;
00197           if(L[j]<0.) L[j]=-L[j];
00198 //          L[j]=rnd();
00199         }
00200 // this is an elastic scattering
00201        L[3]=4.-(L[0]+L[1]+L[2]);
00202 // final position
00203         x[1]=(L[0]*coord[0][noeud_geo[0][EL]-1]+
00204               L[1]*coord[0][noeud_geo[1][EL]-1]+
00205               L[2]*coord[0][noeud_geo[2][EL]-1]+
00206               L[3]*coord[0][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00207         y[1]=(L[0]*coord[1][noeud_geo[0][EL]-1]+
00208               L[1]*coord[1][noeud_geo[1][EL]-1]+
00209               L[2]*coord[1][noeud_geo[2][EL]-1]+
00210               L[3]*coord[1][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00211         z[1]=(L[0]*coord[2][noeud_geo[0][EL]-1]+
00212               L[1]*coord[2][noeud_geo[1][EL]-1]+
00213               L[2]*coord[2][noeud_geo[2][EL]-1]+
00214               L[3]*coord[2][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00215 // update the pseudo-wave vector
00216        KX=(x[1]-x[0])/tau;
00217        KY=(y[1]-y[0])/tau;
00218        KZ=(z[1]-z[0])/tau;
00219        return;
00220     }
00221 // scattering on SiO2
00222 /*    if(i_dom[vois[p][EL]]==SIO2){
00223         double max;
00224         if(L[0]>L[1] && L[0]>L[2]) max=L[0];
00225         if(L[1]>L[0] && L[1]>L[2]) max=L[1];
00226         if(L[2]>L[1] && L[2]>L[0]) max=L[2];
00227         for(j=0;j<3;j++){
00228           L[j]/=max;
00229           if(L[j]<0.) L[j]=-L[j];
00230         }
00231 // this is an elastic scattering
00232        L[3]=4.-(L[0]+L[1]+L[2]);
00233        return;
00234     }*/
00235 // in case of abrupt heterojunction, one has
00236 // to calculate the reflection and transmission
00237 // probability in order to see jow the particle
00238 // behaves.
00239    if(i_dom[EL]!=i_dom[vois[p][EL]] && i_dom[vois[p][EL]]!=SIO2
00240                                     && i_dom[vois[p][EL]]!=-1){
00241 //     printf("############### %d %d\n",EL,vois[p][EL]);
00242 //     printf("*************** %d %d\n\n",i_dom[EL],i_dom[vois[p][EL]]);
00243      double k1,k2,T,R;
00244      double ksquared,thesquareroot,energy;
00245 // calculates the particle energy
00246      ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00247      thesquareroot=sqrt(1.+4.*alphaK[i_dom[EL]][IV]*HHM[i_dom[EL]][IV]*ksquared);
00248      energy=(thesquareroot-1.)/(2.*alphaK[i_dom[EL]][IV]);
00249      if(IV==2) energy+=EG[i_dom[EL]];
00250      energy*=Q; // <--- energy in Joule
00251 // calculates the Transmission probability
00252      k1=sqrt(2.*MSTAR[i_dom[EL]][IV]*M*energy)/HBAR;
00253      if(IV==2 && NOVALLEY[i_dom[vois[p][EL]]]==1){
00254       if((energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][1])*Q)<0.) goto reflection;
00255       k2=sqrt(2.*MSTAR[i_dom[vois[p][EL]]][1]*M*
00256              (energy-(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q))/HBAR;
00257      }
00258      else{
00259       if((energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q)<0.) goto reflection;
00260       k2=sqrt(2.*MSTAR[i_dom[vois[p][EL]]][IV]*M*
00261              (energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q))/HBAR;
00262      }
00263      if(IV==2 && NOVALLEY[i_dom[vois[p][EL]]]==1){
00264        T=pow(2.*(k1/(MSTAR[i_dom[EL]][IV]*M))
00265          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][1]*M)),2.);
00266        R=pow((k1/(MSTAR[i_dom[EL]][IV]*M)-k2/(MSTAR[i_dom[vois[p][EL]]][1]*M))
00267          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][1]*M)),2.);
00268      }
00269      else{
00270        T=pow(2.*(k1/(MSTAR[i_dom[EL]][IV]*M))
00271          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M)),2.);
00272        R=pow((k1/(MSTAR[i_dom[EL]][IV]*M)-k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M))
00273             /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M)),2.);
00274      }
00275 //     printf("k1 = %g     k2 = %g   R = %g  T = %g\n",k1,k2,R*R/(R*R+T*T),T*T/(R*R+T*T)); 
00276      if(rnd()<R*R/(R*R+T*T)) goto reflection;
00277    }
00278 
00279 //   transmission:;
00280 // inside the device (between two different tetrahedra of the same material)
00281    EL=vois[p][EL];
00282 //   printf("new EL=%d\n",EL);
00283 //   printf("%g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00284       memset(&a,0,sizeof(a));
00285       memset(&b,0,sizeof(b));
00286       {int k;
00287        for(k=0;k<4;k++){
00288          a[1][k+1]=coord[0][noeud_geo[k][EL]-1];
00289          a[2][k+1]=coord[1][noeud_geo[k][EL]-1];
00290          a[3][k+1]=coord[2][noeud_geo[k][EL]-1];
00291          a[4][k+1]=1.;
00292       }}
00293       b[1]=4.*x;
00294       b[2]=4.*y;
00295       b[3]=4.*z;
00296       b[4]=4.;
00297 // LU decomposition solver A.x=b
00298 #include "../linear_solver/ludcmp.h"
00299 #include "../linear_solver/lubksb.h"
00300       for(j=0;j<4;j++) L[j]=b[j+1];
00301       if(L[0]>=0. && L[1]>=0. && L[2]>=0. && L[3]>=0.){
00302 //         printf("final %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00303          return;
00304       }
00305       else goto label1;
00306 }
drift.h
