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 : 19 june 2007, Siracusa, Jean Michel Sellier
00025 // Last modified : 19 june 2007, Siracusa, Jean Michel Sellier
00026 
00027 // we calculate here the quantum Bohm potential.
00028 // this approach avoids one the use of the Scroedinger equation
00029 // which is computationally heavy. 
00030 // this method is directly relied to the pilot-wave theory
00031 // of quantum mechanics.
00032 
00033 // See famous Bohm's paper for more informations
00034 
00035 void Bohm_potential(void)
00036 {
00037  int i,j;
00038  double *v;
00039  double *nabla2R;
00040  double *R;
00041  double *ac[4],*bc[4],*cc[4],*dc[4];
00042  double *xi[3], l[4];
00043  double dx,dy,dz;
00044  
00045 // infinitesimal increments = 5 Angstrom
00046  dx=0.5e-9; dy=0.5e-9; dz=0.5e-9;
00047  
00048  v=malloc((Ng+1)*sizeof(double));
00049  nabla2R=malloc((Ng+1)*sizeof(double));
00050  R=malloc((Ng+1)*sizeof(double));
00051  for(i=0;i<3;i++) xi[i]=malloc((Ne+1)*sizeof(double));
00052  for(i=0;i<4;i++){
00053    ac[i]=malloc((Ne+1)*sizeof(double));
00054    bc[i]=malloc((Ne+1)*sizeof(double));
00055    cc[i]=malloc((Ne+1)*sizeof(double));
00056    dc[i]=malloc((Ne+1)*sizeof(double));
00057  }
00058 
00059 // barycentric coordinates of the Guass point
00060  l[0]=0.25;  l[1]=0.25;  l[2]=0.25;  l[3]=0.25;
00061 
00062 // calculation of the center of mass coordinates for every element
00063  for(i=0;i<Ne;i++){
00064       xi[0][i]=l[0]*coord[0][noeud_geo[0][i]-1]
00065               +l[1]*coord[0][noeud_geo[1][i]-1]
00066               +l[2]*coord[0][noeud_geo[2][i]-1]
00067               +l[3]*coord[0][noeud_geo[3][i]-1];
00068       xi[1][i]=l[0]*coord[1][noeud_geo[0][i]-1]
00069               +l[1]*coord[1][noeud_geo[1][i]-1]
00070               +l[2]*coord[1][noeud_geo[2][i]-1]
00071               +l[3]*coord[1][noeud_geo[3][i]-1];
00072       xi[2][i]=l[0]*coord[2][noeud_geo[0][i]-1]
00073               +l[1]*coord[2][noeud_geo[1][i]-1]
00074               +l[2]*coord[2][noeud_geo[2][i]-1]
00075               +l[3]*coord[2][noeud_geo[3][i]-1];
00076 //    printf("%g %g %g\n",xi[0][j][i],xi[1][j][i],xi[2][j][i]);
00077     }
00078 
00079 // calculation of the theta functions for every element
00080  for(i=0;i<Ne;i++){
00081     int l;
00082     int n=4;
00083     int indx[5];
00084     double d;
00085     double a[5][5];
00086     double b[5];
00087 
00088 // definition of the a matrix for the linear system
00089     {int k;
00090      for(k=0;k<4;k++){
00091        a[k+1][1]=coord[0][noeud_geo[k][i]-1];
00092        a[k+1][2]=coord[1][noeud_geo[k][i]-1];
00093        a[k+1][3]=coord[2][noeud_geo[k][i]-1];
00094        a[k+1][4]=1.;
00095      }}
00096 // LU decomposition of the A matrix
00097 #include "../linear_solver/ludcmp.h"
00098     for(l=0;l<4;l++){
00099       if(l==0){b[1]=1.;b[2]=0.;b[3]=0.;b[4]=0.;}
00100       if(l==1){b[1]=0.;b[2]=1.;b[3]=0.;b[4]=0.;}
00101       if(l==2){b[1]=0.;b[2]=0.;b[3]=1.;b[4]=0.;}
00102       if(l==3){b[1]=0.;b[2]=0.;b[3]=0.;b[4]=1.;}
00103 #include "../linear_solver/lubksb.h"
00104       ac[l][i]=b[1]; bc[l][i]=b[2]; cc[l][i]=b[3]; dc[l][i]=b[4];
00105 // printf("%g %g %g %g\n",b[1],b[2],b[3],b[4]);
00106     } // end of l-cycle
00107 
00108 // definition of R field on each node
00109     for(j=0;j<4;j++) R[noeud_geo[j][i]-1]=sqrt(NE[noeud_geo[j][i]-1]);
00110  } // end of i-cycle
00111 
00112 // just in case...
00113  for(i=0;i<Ng;i++) nabla2R[i]=0.;
00114 
00115 // definition of laplacian of R on each node
00116  for(i=0;i<Ne;i++){
00117    double h1,h2,h3;
00118 // second order finite difference scheme
00119    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*(xi[0][i]+dx)+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00120      +R[noeud_geo[1][i]-1]*(ac[1][i]*(xi[0][i]+dx)+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00121      +R[noeud_geo[2][i]-1]*(ac[2][i]*(xi[0][i]+dx)+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00122      +R[noeud_geo[3][i]-1]*(ac[3][i]*(xi[0][i]+dx)+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00123    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00124      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00125      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00126      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00127    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*(xi[0][i]-dx)+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00128      +R[noeud_geo[1][i]-1]*(ac[1][i]*(xi[0][i]-dx)+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00129      +R[noeud_geo[2][i]-1]*(ac[2][i]*(xi[0][i]-dx)+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00130      +R[noeud_geo[3][i]-1]*(ac[3][i]*(xi[0][i]-dx)+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00131    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dx*dx);
00132    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*(xi[1][i]+dy)+cc[0][i]*xi[2][i]+dc[0][i])
00133      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*(xi[1][i]+dy)+cc[1][i]*xi[2][i]+dc[1][i])
00134      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*(xi[1][i]+dy)+cc[2][i]*xi[2][i]+dc[2][i])
00135      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*(xi[1][i]+dy)+cc[3][i]*xi[2][i]+dc[3][i]);
00136    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00137      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00138      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00139      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00140    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*(xi[1][i]-dy)+cc[0][i]*xi[2][i]+dc[0][i])
00141      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*(xi[1][i]-dy)+cc[1][i]*xi[2][i]+dc[1][i])
00142      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*(xi[1][i]-dy)+cc[2][i]*xi[2][i]+dc[2][i])
00143      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*(xi[1][i]-dy)+cc[3][i]*xi[2][i]+dc[3][i]);
00144    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dy*dy);
00145    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*(xi[2][i]+dz)+dc[0][i])
00146      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*(xi[2][i]+dz)+dc[1][i])
00147      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*(xi[2][i]+dz)+dc[2][i])
00148      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*(xi[2][i]+dz)+dc[3][i]);
00149    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00150      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00151      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00152      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00153    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*(xi[2][i]-dz)+dc[0][i])
00154      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*(xi[2][i]-dz)+dc[1][i])
00155      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*(xi[2][i]-dz)+dc[2][i])
00156      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*(xi[2][i]-dz)+dc[3][i]);
00157    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dz*dz);
00158  }
00159 
00160  for(i=0;i<Ng;i++) v[i]=0.; // just in case...
00161  
00162 // updates the potential by means of Bohm potential
00163  for(i=0;i<Ne;i++){
00164   for(j=0;j<4;j++)
00165    if(R[noeud_geo[j][i]-1]!=0.)
00166      V[noeud_geo[j][i]-1]+=-(HBAR*HBAR/(2.*MSTAR[i_dom[i]][1]*M*Q))
00167                          *nabla2R[noeud_geo[j][i]-1]/R[noeud_geo[j][i]-1];
00168 // printf("QV = %g\n",-(HBAR*HBAR/(2.*MSTAR[i_dom[i]][1]*M*Q))*nabla2R[noeud_geo[0][i]-1]/R[noeud_geo[0][i]-1]);
00169 //printf("QV = %g\n",V[noeud_geo[0][i]-1]);
00170  }
00171  
00172 // free precedently allocated memory
00173  free(v);
00174  free(R);
00175  free(nabla2R);
00176  for(i=0;i<3;i++) free(xi[i]);
00177  for(i=0;i<4;i++){
00178    free(ac[i]);
00179    free(bc[i]);
00180    free(cc[i]);
00181    free(dc[i]);
00182  }
00183 }
Bohm_potential.h
