calculation_self_consistent_potential.c++

Go to the documentation of this file.

/*
 *Performs a self-consistent potential calculation

    Copyright (C) 2010 Joseph Pingenot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 of the License, 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 Affero General Public License for more details.

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

*/


#include <libmetacalc/calculation_self_consistent_potential.h++>

#include <math.h>
#include <libkdotp/physical_constants.h++>

using namespace kdp::constants;

namespace kdotp {
  namespace libmetacalc {
    namespace self_consistent {

      double get_Ef_2D(unsigned int Neigs, double* eigs, double* avg_ooms, int* last_eig, double total_density) {
  double sumA=0;
  double sumB=0;
  double Ef;
  for(unsigned int i=0; i<Neigs; i++) {
    /*Sum up from the individual subbands*/
    fprintf(stderr, "i=%u\n", i);
    sumA = sumB = 0;
    for(int j=i; j>=0; j--) {
      fprintf(stderr, "\tj=%d sumA=%g sumB=%g avg_ooms[%d]=%g eigs[%d]=%g\n", j, sumA, sumB, j, avg_ooms[j], j, eigs[j]);
      sumA += eigs[j]/avg_ooms[j];
      sumB += 1.0/avg_ooms[j];
    }
    Ef = (total_density*hbar*hbar/(2*M_PI*m0) + sumA)/sumB;
    fprintf(stderr, " Ef=%g: ", Ef);
    if((i==(Neigs-1)) || (Ef <= eigs[i+1])) {
      if(i==(Neigs-1)) {
        fprintf(stderr, " <eigs[%d]=%g\n", i+1, eigs[i+1]);
      }else{
        fprintf(stderr, " i==Neigs(%u)\n", Neigs);
      }
      *last_eig = i;
      return Ef;
    }else{
      fprintf(stderr, "\n");
    }
  }
  *last_eig=Neigs+1;
  return 0;
      }

      /**callback for mixing the potentials (search for "Mix the potentials" below)
       *\note simple: just add the two together!
       */
      void potential_mixing_callback(int X, unsigned N, double dx, double *restrict to, double from, void* mixing_pct) {
  *to = from + (*to - from)*(*(double*)mixing_pct);
      }

    }
  }
}