1 #ifndef ELECTRON_DENSITY_H
2 #define ELECTRON_DENSITY_H
3 
4 #include <map>
5 #include <string>
6 #include <fstream>
7 #include "ATC_TypeDefs.h"
8 #include "Function.h"
9 
10 const double tol = 1.e-8;
11 
12 namespace ATC {
13 
14   /**
15    *  @class  ElectronChargeDensity
16    *  @brief  Base class for models of extrinsic electric charges
17    */
18 
19   class ElectronChargeDensity
20   {
21     public:
ElectronChargeDensity()22       ElectronChargeDensity()   {};
~ElectronChargeDensity()23       virtual ~ElectronChargeDensity() {};
electron_charge_density(const FIELD_MATS &,DENS_MAT &)24       virtual bool electron_charge_density(const FIELD_MATS & /* fields */,
25                                            DENS_MAT & /* flux */) const { return false; };
26 
27 
D_electron_charge_density(const FieldName,const FIELD_MATS &,DENS_MAT &)28       virtual void D_electron_charge_density(const FieldName /* fieldName */,
29                                              const FIELD_MATS & /* fields */,
30                                              DENS_MAT & /* flux */) const
31         { throw ATC_Error("Charge density D_electron_charge_density unimplemented function");}
32 
band_edge_potential(const FIELD_MATS &,DENS_MAT &)33       virtual void band_edge_potential(const FIELD_MATS & /* fields */,
34                                        DENS_MAT & /* density */) const
35         { throw ATC_Error("Charge density band_edge_potential unimplemented function");}
36   };
37   //-----------------------------------------------------------------------
38   /**
39    *  @class  ElectronChargeDensityInterpolation
40    *  @brief  Class for models of electron charge density as a tabular function of electric potential
41    */
42 
43   class ElectronChargeDensityInterpolation : public ElectronChargeDensity
44   {
45     public:
46       ElectronChargeDensityInterpolation(std::fstream &matfile,std::map<std::string,double> & parameters);
~ElectronChargeDensityInterpolation()47       virtual ~ElectronChargeDensityInterpolation() {};
electron_charge_density(const FIELD_MATS & fields,DENS_MAT & flux)48       virtual bool electron_charge_density(const FIELD_MATS &fields,
49                                            DENS_MAT &flux) const
50       {
51         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
52         const DENS_MAT & phi = phi_field->second;
53         int nNodes  = phi.nRows();
54         flux.reset(nNodes,1,false);
55         for (int i = 0; i < nNodes; i++) {  // a mapping of a vector
56           flux(i,0) = n_.f(phi(i,0));
57         }
58         flux *= -1.;
59         return true;
60       };
D_electron_charge_density(const FieldName,const FIELD_MATS & fields,DENS_MAT & coef)61       virtual void D_electron_charge_density(const FieldName /* field */,
62                                              const FIELD_MATS &fields,
63                                              DENS_MAT &coef) const
64       {
65         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
66         const DENS_MAT & phi = phi_field->second;
67         int nNodes  = phi.nRows();
68         coef.reset(nNodes,1,false);
69         for (int i = 0; i < nNodes; i++) {
70           coef(i,0) = n_.dfdt(phi(i,0));
71           coef(i,0) = n_.dfdt(phi(i,0));
72         }
73         coef *= -1.;
74       }
75     private:
76       InterpolationFunction n_;
77   };
78   //-----------------------------------------------------------------------
79   /**
80    *  @class  ElectronChargeDensityLinear
81    *  @brief  Class for models of electron charge density proportional to electric potential
82    */
83 
84   class ElectronChargeDensityLinear : public ElectronChargeDensity
85   {
86     public:
87       ElectronChargeDensityLinear(std::fstream &matfile,std::map<std::string,double> & parameters);
~ElectronChargeDensityLinear()88       virtual ~ElectronChargeDensityLinear() {};
electron_charge_density(const FIELD_MATS & fields,DENS_MAT & flux)89       virtual bool electron_charge_density(const FIELD_MATS &fields,
90                                            DENS_MAT &flux) const
91       {
92         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
93         flux = phi_field->second;
94         flux *= -C_;
95         return true;
96       };
D_electron_charge_density(const FieldName,const FIELD_MATS & fields,DENS_MAT & coef)97       virtual void D_electron_charge_density(const FieldName /* field */,
98                                              const FIELD_MATS &fields,
99                                              DENS_MAT &coef) const
100       {
101         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
102         const DENS_MAT & phi = phi_field->second;
103         int nNodes  = phi.nRows();
104         coef.reset(nNodes,1,false);
105         coef = -C_;
106       }
107     private:
108       double C_;
109   };
110   //-----------------------------------------------------------------------
111   /**
112    *  @class  ElectronChargeDensityExponential
113    *  @brief  Class for models of electron charge density dependent on difference between electric potential and the Fermi level n = n_i exp ( (phi-E_i) / kB T)
114    */
115 
116   class ElectronChargeDensityExponential : public ElectronChargeDensity
117   {
118     public:
119       ElectronChargeDensityExponential(std::fstream &matfile,std::map<std::string,double> & parameters);
~ElectronChargeDensityExponential()120       virtual ~ElectronChargeDensityExponential() {};
121 
n(const double phi,double T)122       double n(const double phi, double T)  const
123       {
124         return -intrinsicConcentration_*exp((phi-intrinsicEnergy_)/(kBeV_*T));
125       }
dndphi(const double phi,double T)126       double dndphi(const double phi, double T)  const
127       {
128         return n(phi,T)/(kBeV_*T);
129       }
electron_charge_density(const FIELD_MATS & fields,DENS_MAT & density)130       virtual bool electron_charge_density(const FIELD_MATS &fields,
131                      DENS_MAT &density) const
132       {
133         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
134         FIELD_MATS::const_iterator T_field   = fields.find(TEMPERATURE);
135         double T = 300;
136         bool hasTref = (referenceTemperature_ > 0 );
137         const DENS_MAT & phi = phi_field->second;
138         int nNodes = phi.nRows();
139         density.resize(nNodes,1);
140         if (hasTref) {
141           T = referenceTemperature_;
142           for (int i = 0; i < nNodes; i++) {
143             density(i,0) = n(phi(i,0),T); }
144         }
145         else {
146           const DENS_MAT & temp = T_field->second;
147           for (int i = 0; i < nNodes; i++) {
148             density(i,0) = n(phi(i,0),temp(i,0)); }
149         }
150         density *= -1.;
151         return true;
152       };
D_electron_charge_density(const FieldName,const FIELD_MATS & fields,DENS_MAT & coef)153       virtual void D_electron_charge_density(const FieldName /* field */,
154                                              const FIELD_MATS &fields,
155                                              DENS_MAT &coef) const
156       {
157         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
158         FIELD_MATS::const_iterator T_field   = fields.find(TEMPERATURE);
159         double T = 300;
160         bool hasTref = (referenceTemperature_ > 0 );
161         const DENS_MAT & phi = phi_field->second;
162         int nNodes = phi.nRows();
163         coef.resize(nNodes,1);
164         if (hasTref) {
165           T = referenceTemperature_;
166           for (int i = 0; i < nNodes; i++) {
167             coef(i,0) = dndphi(phi(i,0),T); }
168         }
169         else {
170           const DENS_MAT & temp = T_field->second;
171           for (int i = 0; i < nNodes; i++) {
172             coef(i,0) = dndphi(phi(i,0),temp(i,0)); }
173         }
174         coef *= -1.;
175       };
176     protected:
177       double intrinsicConcentration_,intrinsicEnergy_;
178       double referenceTemperature_;
179   };
180 
181   //-----------------------------------------------------------------------
182   /**
183    *  @class  ElectronChargeDensityFermiDirac
184    *  @brief  Class for models of electron charge density based on Fermi-Dirac statistics
185    */
186 
187   class ElectronChargeDensityFermiDirac : public ElectronChargeDensity
188   {
189     public:
190       ElectronChargeDensityFermiDirac(std::fstream &matfile,std::map<std::string,double> & parameters);
~ElectronChargeDensityFermiDirac()191       virtual ~ElectronChargeDensityFermiDirac() {};
fermi_dirac(const double E,const double T)192       double fermi_dirac(const double E, const double T) const
193       {
194         double f = 1.0;
195         if      (T > 0)   f = 1.0 / ( exp((E-Ef_)/kBeV_/T)+1.0 );
196         else if (E > Ef_) f = 0;
197         return f;
198       };
electron_charge_density(const FIELD_MATS & fields,DENS_MAT & density)199       virtual bool electron_charge_density(const FIELD_MATS &fields,
200                                            DENS_MAT &density) const
201       {
202         // psi : the inhomogeneous solution
203         FIELD_MATS::const_iterator psi_field = fields.find(ELECTRON_WAVEFUNCTION);
204 
205         const DENS_MAT & psi = psi_field->second;
206         FIELD_MATS::const_iterator psis_field = fields.find(ELECTRON_WAVEFUNCTIONS);
207         // if (psis_field==fields.end())
208         //throw ATC_Error("Wavefunctions not defined");
209         const DENS_MAT & psis = psis_field->second;
210         FIELD_MATS::const_iterator E_field = fields.find(ELECTRON_WAVEFUNCTION_ENERGIES);
211         const DENS_MAT & Es = E_field->second;
212         FIELD_MATS::const_iterator T_field   = fields.find(ELECTRON_TEMPERATURE);
213         const DENS_MAT & Ts = T_field->second;
214         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
215         const DENS_MAT & phi = phi_field->second;
216 
217         int nNodes = psi.nRows();
218         density.reset(nNodes,1);
219         double T  = referenceTemperature_;
220         int count = 0;
221         for (int i = 0; i < nNodes; i++) {
222           if (!hasReferenceTemperature_) { T = Ts(i,0); }
223           int j = 0;
224           for (j = 0; j < psis.nCols(); j++) {
225             double E = Es(j,0); // Eigenvalue
226             double f = fermi_dirac(E,T);
227             if (f < tol) break;
228             else  count++;
229             density(i,0) -= psis(i,j)*psis(i,j)*f;  // < 0
230           }
231           if (donorIonization_) {
232             double E = -1.0* phi(i,0);// units [eV] E = - |e| phi
233             if ( E + Eb_ > Ef_+Ed_) density(i,0) += Nd_;  // > 0
234           }
235         }
236         return true;
237       };
D_electron_charge_density(const FieldName,const FIELD_MATS & fields,DENS_MAT & coef)238       virtual void D_electron_charge_density(const FieldName /* fieldName */,
239                                              const FIELD_MATS &fields,
240                                              DENS_MAT &coef) const
241       {
242         FIELD_MATS::const_iterator phi_field = fields.find(ELECTRIC_POTENTIAL);
243         const DENS_MAT & phi = phi_field->second;
244         int nNodes  = phi.nRows();
245         coef.reset(nNodes,1,false);
246       }
247 
band_edge_potential(const FIELD_MATS & fields,DENS_MAT & density)248       virtual void band_edge_potential(const FIELD_MATS &fields,
249                                        DENS_MAT &density) const
250       {
251         FIELD_MATS::const_iterator p_field   = fields.find(ELECTRIC_POTENTIAL);
252         const DENS_MAT & phi = p_field->second;
253         int nNodes  = phi.nRows();
254         density.reset(nNodes,1,false);
255         density = Eb_;
256       };
257 
258     protected:
259       double Ef_;
260       double referenceTemperature_;
261       double Ed_, Nd_;
262       double Eb_;
263       bool hasReferenceTemperature_, donorIonization_;
264   };
265 }
266 
267 #endif
268 
269 
270