Discription:

The Fermi-Energy is calculated in the temperature range T_min <= T <= T_max at max temperatures. For the upper edge of the valence-band E_V = 7.00eV was chosen. Change the parameter to calculate the Fermi-Energy for a other configuration. The chosen value for E_G = E_C - E_V = 1.12eV correspond to the band-gap of silicon. For the effective mass m_n^* = m_n^* = 0.3·m_e is chosen where m_e is the electron mass.

#include<math.h>
#include<stdio.h>
const float EC  = 7.00 + 1.12 ; /* lower edge of the conduction-band in eV */
const float EV  = 7.00        ; /* upper edge of the valence   -band in eV */
const float EA  = 7.00  + 0.03; /* level of the acceptors            in eV */
const float ED  = 7.00  + 1.12 - 0.03; /* level of the donators      in eV */
const float NA  =     1.00E+12; /* acceptor-concentration         in cm^-3 */
const float ND  =     1.00E+15; /* donntor -concentration         in cm^-3 */
const float kB  =  1.38066E-23; /* Boltzman konstant                in J/K */
const float e   =  1.60230E-19; /* elemetal charge                    in C */
const float mstar = 0.3*9.10939E-31; /* effevtive mass               in kg */
const float hplanck =  6.626076E-34; /* Planck konstant              in Js */ 
const float pi      =   3.141592654; 
const float TempMin =          0.00; /* min temperature               in K */
const float TempMax =       1000.00; /* max temperature               in K */
const int   max     =           200; /* number of computations             */

FILE *out;
int i;
double temp,ct,ef;

 float SearchEF(float actfermi)
 {
  float exparg1,exparg2,exparg3,exparg4,Faktor;
  float exp1,exp2,exp3,exp4,dummy,ctemp;

  Faktor  = e/(kB*temp);
  exparg1 = (actfermi -       ED)*Faktor; /* EFn-Term */
  exparg2 = (EV       - actfermi)*Faktor; /*   p-Term */
  exparg3 = (actfermi -       EC)*Faktor; /*   n-Term */
  exparg4 = (EA       - actfermi)*Faktor; /* EFp-Term */
  if (exparg1 > - 70.00) exp1 = exp(exparg1); else exp1 = 0.00;
  if (exparg2 > - 70.00) exp2 = exp(exparg2); else exp2 = 0.00;
  if (exparg3 > - 70.00) exp3 = exp(exparg3); else exp3 = 0.00;
  if (exparg4 > - 70.00) exp4 = exp(exparg4); else exp4 = 0.00;
  ctemp = exp(3.00/2.00*log(temp));
  dummy = 1E+6*(ND/(1.00 + exp1) - NA/(1.00 + exp4)) + ct*ctemp*(exp2 - exp3);
  return dummy;
 }

 int IterationFinished(double IFW, double IFD)
 {
  int h1,h2;

  h1 = 0; h2 = 0;
  if (fabs(IFW) < 1E-6) h1 = 1;
  if (IFD < 1E-10) h2 = 1;
  h1 = !(h1 || h2);
  return h1;
 }

 double CalculateEF(void)
 {
    double exparg,dummy;
    double EFL,EFM,EFR;
    double W3,delta;

    ct = 2.00*pi*mstar*kB/(hplanck*hplanck);
    ct = 2.00*exp(3.00/2.00*log(ct)); /* use 2.00 or (float)2 */
    if (temp == 0.00)
    {
     return (EC + ED)/2.00;
    }
    else
    {
     exparg = (EC - ED)*e/(kB*temp);
     if (exparg > 70.00)
     {
      dummy  = exp(3.00/2.00*log(temp));
      dummy  = (EC + ED)/2.00 + kB*temp/(2.00*e)*log(1E+16*ND/(ct*dummy));
      return dummy;
     }
     else
     {
      EFL = EV;
      EFR = EC;
      do
      {
       EFM = (EFL + EFR)/2.00;
       W3  = SearchEF(EFM);
       if (W3 > 0.00)
       {
        EFL   = EFM;
        EFR   = EFR;
        delta = fabs(EFR - EFM);
       }
       else
       {
        EFL   = EFL;
        EFR   = EFM;
        delta = fabs(EFL - EFM);
       }
 /*    printf("D = %f",delta);
       printf("  Error = %f\n",W3); */
      } /* do */
      while (IterationFinished((EFR-EFL)/EFM,delta) == 1);
      return EFM;
     } /* if (exparg > 70) */
    } /* if (temp == 0) */
  }

void main(void)
{
  out = fopen("FERMI.DAT","w");
  printf("\n\n\n computation of the fermi-enery temperature dependence\n");
  printf(" by Jens Korallus, October 1999 \n\n");
  for (i = 1; i <= max; i++)
    {
      temp = ((TempMax - TempMin)/(max - 1)*(i - 1)) + TempMin;
      ef   = CalculateEF();
      printf(" T = %f",temp);
      printf(" EF = %f\n",ef);
      fprintf(out,"%f",temp);
      fprintf(out,"  %f\n",ef);
    }
  printf("finished!");
}