calculation_self_consistent_potential.h++

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 <math.h>
#include <glib.h>
#include <libkdotp/wavefunction_cartesian_effective_mass.h++>
#include <libpostproc/density_1d.h>
#include <libpostproc/potential_1d.h>
#include <libpostproc/fftw_wavefunc.h>
#include <libkdotp/hamiltonians.h++>
#include <libmodelxx/cartesian_grid.h++>
#include <libmodelxx/generic_structure.h++>
#include <libkdotp/physical_constants.h++>

#ifndef ___CALCULATION_SELF_CONSISTENT_POTENTIAL___
#define ___CALCULATION_SELF_CONSISTENT_POTENTIAL___

#define CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN 8438


namespace kdotp {
  namespace libmetacalc {
    namespace self_consistent {

      /*This is a function pointer that needs passed to get the effective mass*/
      typedef double (*ooms_averager)(double* block, unsigned int eigvec, unsigned int L, double* func);
      typedef double (*wavefunction_max_delta_callback)(double* block1, unsigned vec1, double* block2, unsigned vec2, unsigned L);

      //typedef void* (*potential_calculator_callback)(void* charge_density, void* privdat);

      /**Gets the Fermi energy for a set of eigenvalues (2D DoS)
       *\param Neigs number of eigenvalues in the array (extras are OK; just don't have fewer eigs than this
       *\param eigs the array of eigenvalues
       *\param last_eig pointer to an unsigned int to store the last eig.  If this is the last one, we can't be sure that Ef is the correct Ef; there may be a higher-energy subband that would fill!
       *\return the Fermi energy
       */
      double get_Ef_2D(unsigned int Neigs, double* eigs, double* avg_ooms, int* last_eig, double total_density);

      /**callback for mixing the potentials (search for "Mix the potentials" below)*/
      void potential_mixing_callback(int X, unsigned N, double dx, double *restrict to, double from, void* mixing_pct);

      /*The point of N_eigs_to_check_for_convergence is to prevent cycles by looking back N eigs.*/
      template<class hamiltonian, unsigned int D, unsigned N_eigs_to_check_for_convergence=2>
      class self_consistent_potential {
      public:
  ~self_consistent_potential() {
    /*We do a bunch of these frees.  Cluster them.  If they're not NULL, free them.
     *Note that these are pointers to pointers; we keep going until the array entry is NULL (if it's a real member of this class, it won't be NULL
     */
    /*Note that we don't care about the generic structure; we were handed it by someone else*/
    void** things_to_free[] = {
      (void**)&dopant_density,
      NULL
    };
    for(int i=0; things_to_free[i] != NULL; i++) {
      if(*things_to_free[i] != NULL) {
        free(*things_to_free[i]);
        *things_to_free[i] = NULL;
      }
    }
    delete H;
  }
  self_consistent_potential(
          struct ::kdotp::libmodelxx::structure::generic_structure *s,
          enum ::kdotp::libkdotp::hamiltonian::boundary_condition bc,
          ooms_averager averager_,
          wavefunction_max_delta_callback deltar_, double abserr_, double relerr_)
    : generic_structure(s),
      H(NULL),
      averager(averager_),
      deltar(deltar_),
      dopant_density(NULL),
      abserr(abserr_),
      relerr(relerr_) {
    H = new hamiltonian(s, bc, false);
    if(H == NULL) return;
    dopant_density = new_density_1d(((::kdotp::libmodelxx::grid::cartesian_grid<double, D>*)(s->doping_grid))->L[0], ((::kdotp::libmodelxx::grid::cartesian_grid<double, D>*)(s->doping_grid))->a[0], 0);
    for(unsigned int i=0; i<((::kdotp::libmodelxx::grid::cartesian_grid<double, D>*)(s->doping_grid))->L[0]; i++) {
      /*Dopant density is in cm^-3; convert to m^-3*/
      dopant_density->storage[i] = (((::kdotp::libmodelxx::grid::cartesian_grid<double, D>*)(s->doping_grid))->data)[i]*100*100*100;
      //fprintf(stderr, "density[%d] = %g\n", i, dopant_density->storage[i]);
    }
  }
  /**Performs the actual calculation.
   *\param Neigs number of eigenvalues to use to do the calculation
   *\param max_iterations quit after max_iterations has been completed (avoids infinite loops; this can usually be a small value)
   *\param convergence the convergence criterion; the max change in probability density must be less than this value to be converged.
   *\param max_converged_eigenstate make sure all of the eigenstates up through this one are converged.
   *\param eigenvectors pointer to hold the address of the array of eigenvectors
   *\param eigenvalues pointer to hold the address of the array of eigenvalues
   *\param ooms_parameter parameter to pass to the Hamiltonian to get the static dielectric constant
   *\param density_integration_error estimate of the error from integrating the density.
   *\param percent_current percent of the previous potential to mix with 1-percent_current of the second-to-last to determine the potential of the current calculation.
   *\return a GError containing error information.  NULL if successful.
   */
  inline GError* perform_calculation(unsigned int Neigs, unsigned int max_iterations, double convergence, unsigned int max_converged_eigenstate, double** eigenvectors, double**eigenvalues, int ooms_parameter, int epsilon_parameter, double* Ef, double *density_integration_error, double *potential_integration_error, double* efield_integration_error, GFile* output_directory, double percent_current=1) {
    return perform_calculation(false, 0, 0, Neigs, max_iterations, convergence, max_converged_eigenstate, eigenvectors, eigenvalues, ooms_parameter, epsilon_parameter, Ef, density_integration_error, potential_integration_error, efield_integration_error, output_directory, percent_current);
  }
  inline GError* perform_calculation(double e_min, double e_max, unsigned int max_iterations, double convergence, unsigned int max_converged_eigenstate, double** eigenvectors, double**eigenvalues, int ooms_parameter, int epsilon_parameter, double* Ef, double *density_integration_error, double *potential_integration_error, double* efield_integration_error, GFile* output_directory, double percent_current=1) {
    return perform_calculation(true, e_min, e_max, 0, max_iterations, convergence, max_converged_eigenstate, eigenvectors, eigenvalues, ooms_parameter, epsilon_parameter, Ef, density_integration_error, potential_integration_error, efield_integration_error, output_directory, percent_current);
  }
  GError* perform_calculation(bool use_energies, double e_min, double e_max, unsigned int Neigs, unsigned int max_iterations, double convergence, unsigned int max_converged_eigenstate, double** eigenvectors, double**eigenvalues, int ooms_parameter, int epsilon_parameter, double* Ef, double *density_integration_error, double *potential_integration_error, double* efield_integration_error, GFile* output_directory, double percent_current=1) {
    double* average_ooms;
    if(use_energies) {
      /*For safety's sake, take many eigenvalues*/
      Neigs = 50;
    }
    fprintf(stderr, "Neigs=%u\n", Neigs);
    average_ooms = new double[Neigs];
    double* mstar = new double[H->L];
    double* epsilon = new double[H->L];
    struct density_1d epsilon_density;
    double *evecs;
    double *evals;
    int Neigs_obtained;
    GError *err = NULL;
    struct density_1d *subband_density;
    struct density_1d *total_density;
    struct density_1d *total_potential=NULL;
    struct density_1d *previous_potential=NULL;
    char fn[4096];
    GFile* gf;
    GOutputStream* gos;
    GString* gs = g_string_new("");
    gsize wlen;
    /*dx in the total doping was in units of nm; we need m here.*/
    double total_doping = density_1d_integrate_density(dopant_density, abserr, relerr, density_integration_error, &err, 0)*1e-9;
    //double total_doping = density_1d_sum(dopant_density, CONSTANT)*(dopant_density->dx)*1e-9;
    fprintf(stderr, "total doping=%g\n", total_doping);
    if(err != NULL) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Failed to integrate dopant charge density (before iterating)");
      return e;
    }
    if(output_directory != NULL) {
      gf = g_file_get_child(output_directory, "doping.pdata");
      gos = (GOutputStream*)g_file_replace(gf, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &err);
      if(err != NULL) {
        GError *e = NULL;
        g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR opening bands file (%s)", g_file_get_uri(gf));
        return e;
      }
      g_string_printf(gs, "#gridsite\tdensity(*1e18)\n");
      if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
        GError *e = NULL;
        g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
        return e;
      }
      for(unsigned int i=0; i<H->L; i++) {
        g_string_printf(gs, "%u\t%g\n", i, dopant_density->storage[i]);
        if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
    return e;
        }
      }
      g_object_unref(gos);
    }

    int last_filled_eig;
    double* last_evecs[N_eigs_to_check_for_convergence];
    for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
      last_evecs[eig]=NULL;
    }
    bool converged = false;
    /*The effective mass across the structure won't change between iterations*/
    H->get_parameter_function(ooms_parameter, mstar);
    H->get_parameter_function(epsilon_parameter, epsilon);
    epsilon_density.storage = epsilon;
    epsilon_density.N = H->L;
    /*H has units in m; we need nm*/
    epsilon_density.dx = H->dx;
    total_density = new_density_1d(H->L, H->dx, 0);
    if(unlikely(total_density == NULL)) {
      delete[] average_ooms;
      delete[] mstar;
      delete[] epsilon;
      return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 5, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: failed to allocate storage for the total subband density");
    }
    subband_density = new_density_1d(H->L, H->dx, 0);
    if(unlikely(subband_density == NULL)) {
      free_density_1d(total_density);
      delete[] average_ooms;
      delete[] mstar;
      delete[] epsilon;
      return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 5, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: failed to allocate storage for the total subband density");
    }
    /*The internal potential is currently 0.*/
    for(unsigned int iter=0; iter<max_iterations; iter++) {
      fprintf(stderr, "ITER %u\n\n", iter);
      if(use_energies) {
        evecs = H->get_eigenvectors_double(e_min, e_max, &evals, &Neigs_obtained, &err);
      }else{
        evecs = H->get_eigenvectors_int(1, (int)Neigs, &evals, &Neigs_obtained, &err);
      }
      if(unlikely(err != NULL)) {
        *eigenvectors = evecs;
        *eigenvalues = evals;
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        GError *e = NULL;
        g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Failed to get eigenvalues in iteration %u/%u.", iter, max_iterations);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return e;
      }
      if(unlikely(Neigs_obtained != (int)Neigs)) {
        *eigenvectors = evecs;
        *eigenvalues = evals;
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 2, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Got incorrect number of eigenvalues from LAPACK: %u requested, %d obtained @iteration %u", Neigs, Neigs_obtained, iter);
      }
      for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
        if(last_evecs[eig] != NULL) {
    for(unsigned int i=0; i<max_converged_eigenstate; i++) {
      double delta;
      delta = deltar(evecs, i, last_evecs[eig], i, H->L);
      if(delta < 0) {
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
          if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 7, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Got error in deltar: %g\n", -delta);
      }
      if(convergence < delta) {
        fprintf(stderr, "%u: delta=%g UNCONVERGED (convergence=%g)\n", i, delta, convergence);
        converged = false;
        break;
      }else{
        fprintf(stderr, "%u: delta=%g CONVERGED (convergence=%g)\n", i, delta, convergence);
      }
      converged = true;
    }
    /*if just one converges, break.*/
    if(converged) break;
        }
      }
      fprintf(stderr, "Converged? %s\n", converged?"Yes":"No");
      /*Verify our potential*/
      double* potcheck = (double*)malloc(sizeof(double)*(H->L));
      if(potcheck == NULL) {
        return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 15, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Failed to allocate potential check.");
      }
      H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
      //for(unsigned i=0; i<H->L; i++) fprintf(stderr, "pc@%u: %g\n", i, potcheck[i]);
      free(potcheck);
      if(converged) {
        *eigenvectors = evecs;
        *eigenvalues = evals;
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return NULL;
      }
      for(unsigned int i=0; i<Neigs; i++)
        average_ooms[i] = averager(evecs, i, H->L, mstar);
      *Ef = get_Ef_2D(Neigs, evals, average_ooms, &last_filled_eig, total_doping);
      if(unlikely(last_filled_eig < 0)) {
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 3, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: @iter %u: last filled eigenstate (%u) is invalid; no levels are filled!", iter, last_filled_eig);
      }
      if(unlikely((unsigned)last_filled_eig >= Neigs)) {
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 3, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: @iter %u: last filled eigenstate (%u) is the same as the highest calculated eigenstate (%u); we can't be sure there aren't higher filled levels!", iter, Neigs, last_filled_eig);
      }
      init_density_1d(total_density, 0);
      for(unsigned int subband=0; subband <= (unsigned)last_filled_eig; subband++) {
        if(unlikely((err = init_density_1d_from_cartesian_effective_mass_block(evecs, subband, H->L, H->dx, subband_density)) != NULL)) {
    for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
      if(last_evecs[eig] != NULL) free(last_evecs[eig]);
    }
    free_density_1d(subband_density);
    if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
    free_density_1d(total_density);
    delete[] average_ooms;
    delete[] mstar;
    delete[] epsilon;
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: @iter %u: failed to set density from wavefunction for subband %u: ", iter, subband);
    return e;
        }
        GError *eerr=NULL;
        double dier;
        double sum = 0;
        for(unsigned i=0; i<subband_density->N; i++) sum += subband_density->storage[i];
        fprintf(stderr, "subband %u: density of this subband is: %g (sum is %g)\n", subband, density_1d_integrate_density(subband_density, abserr, relerr, &dier, &eerr, 1), sum);
        /*Multiply the single-electron charge density from the eigenstate by the number of electrons in this state.*/
        density_1d_mul_by_const(subband_density, 2*M_PI*kdp::constants::m0/(kdp::constants::hbar*kdp::constants::hbar*average_ooms[subband])*((*Ef)-evals[subband]));
        fprintf(stderr, "Subband %u: after multpliying by total areal charge density (%g), have total areal cahrge density of %g\n", subband, 2*M_PI*kdp::constants::m0/(kdp::constants::hbar*kdp::constants::hbar*average_ooms[subband])*((*Ef)-evals[subband]), density_1d_integrate_density(subband_density, abserr, relerr, &dier, &eerr, 2));
        op_density_1d_to_density_1d(total_density, ADD_EQUALS, subband_density);
      }
      /*We have |Psi|^2*dx (so the sum is 1) in the total density right now.  We need to convert
       *  it to a charge density.
       * |Psi|^2*dx / dx = |Psi|^2 -> 1/m
       *Note that we also multiply by a minus sign (negative charge; the magnitude will come in a few lines, when we add in the dopants
       */
      density_1d_mul_by_const(total_density, -1.0/(total_density->dx*1e-9));
      double edensity_ierr;
      GError *eee=NULL;
      double total_electron_density = density_1d_integrate_density(total_density, abserr, relerr, &edensity_ierr, &eee, 3)*1e-9;
      /*Renormalize such that we multiply by the percentage we actually want*/
      density_1d_mul_by_const(total_density, fabs(1.0-(fabs(total_electron_density)-fabs(total_doping))/fabs(total_doping)));
      double edensity2_ierr;
      double total_electron_density_2 = density_1d_integrate_density(total_density, abserr, relerr, &edensity2_ierr, &eee, 6)*1e-9;
      //fprintf(stderr, "total_density->dx = %g\n", total_density->dx);
      op_density_1d_to_density_1d(total_density, ADD_EQUALS, dopant_density);
      double tcdensity_ierr;
      double total_calculated_density = density_1d_integrate_density(total_density, abserr, relerr, &edensity_ierr, &eee, 4)*1e-9;
      fprintf(stderr, "total dopant density: %.15g total electron density: %.15g sum: %.15g -> %.15g (%g): total calculated density: %.15g \n", total_doping, total_electron_density, total_doping + total_electron_density, total_electron_density_2,  total_doping + total_electron_density_2, total_calculated_density);
      /*Now multiply by the magnitude of the charge of an electron*/
      /*Multiply the total potential by the grid spacing, since we're taking the second derivative.*/
      density_1d_mul_by_const(total_density, kdp::constants::e);
      //total_potential = get_potential_for_density_1d(total_density, &epsilon_density, NULL);
      total_potential = potential_1d_get_potential(0.0, 0.0, total_density, &epsilon_density, INTEGRATION_METHOD_NP, potential_integration_error, efield_integration_error, relerr, abserr, &err);
      if((total_potential == NULL) || (err != NULL)) {
        if(err == NULL) {
    err = g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 4, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: @iter %u/%u: error getting potential from total electron density distribution (got NULL)", iter, max_iterations);
        }
        GError *e = NULL;
        g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: @iter %u/%u: failed to get potential for charge density: ", iter, max_iterations);
        for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
    if(last_evecs[eig] != NULL) free(last_evecs[eig]);
        }
        free_density_1d(subband_density);
        if(total_potential != NULL) free_density_1d(total_potential);
        if(previous_potential != NULL) free_density_1d(previous_potential);
        free_density_1d(total_density);
        delete[] average_ooms;
        delete[] mstar;
        delete[] epsilon;
        return err;
      }
      /*Potential was integrated twice with x in the implicit-nm basis; convert to m*/
      density_1d_mul_by_const(total_potential, 1e-9*1e-9);
      if(percent_current < 1) {
        /*Mix the potentials*/
        if(previous_potential != NULL) {
    op_density_1d_to_density_1d_cb(total_potential, previous_potential, potential_mixing_callback, (void*)&percent_current);
    free(previous_potential);
        }
        /*Now that we're done mixing, update.*/
        previous_potential = new_density_1d_copy(total_potential);
      }
      /**Write out the iteration*/
      if(output_directory != NULL) {
        double sum, val;
        for(unsigned int j=0; j<Neigs; j++) {
    g_string_printf(gs, "i%u_eigvec_%d.pdata", iter, j);
    gf = g_file_get_child(output_directory, gs->str);
    gos = (GOutputStream*)g_file_replace(gf, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &err);
    if(err != NULL) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR opening eigvec file (%s)", g_file_get_uri(gf));
      return e;
    }
    g_string_printf(gs, "#gridsite\tpsi\n");
    if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
      return e;
    }
    for(unsigned int i=0; i<H->L; i++) {
      g_string_printf(gs, "%u\t%g\n", i, evecs[i + j*H->L]);
      if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
        GError *e = NULL;
        g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
        return e;
      }
    }
    g_object_unref(gos);
        }
        g_string_printf(gs, "i%u_pot.pdata", iter);
        gf = g_file_get_child(output_directory, gs->str);
        gos = (GOutputStream*)g_file_replace(gf, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &err);
        if(err != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR opening eigvec file (%s)", g_file_get_uri(gf));
    return e;
        }
        g_string_printf(gs, "#gridsite\tpot(eV)\n");
        if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
    return e;
        }
        for(unsigned int i=0; i<total_potential->N; i++) {
    g_string_printf(gs, "%u\t%g\n", i, total_potential->storage[i]);
    if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
      return e;
    }
        }
        g_object_unref(gos);

        g_string_printf(gs, "i%u_total_density.pdata", iter);
        gf = g_file_get_child(output_directory, gs->str);
        gos = (GOutputStream*)g_file_replace(gf, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &err);
        if(err != NULL) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR opening density file (%s)", g_file_get_uri(gf));
    return e;
        }
        g_string_printf(gs, "#gridsite\tdensity(m^-3)*a^2\n");
        if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
    GError *e = NULL;
    g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
    return e;
        }
        for(unsigned int i=0; i<total_density->N; i++) {
    g_string_printf(gs, "%u\t%g\n", i, total_density->storage[i]);
    if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &err)) {
      GError *e = NULL;
      g_propagate_prefixed_error(&e, err, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: ERROR writing to bands file (%s)", g_file_get_uri(gf));
      return e;
    }
        }
        g_object_unref(gos);
      }
      /**Final step: copy the new potential through.
       *Make sure we copy off the values into a new density so that we don't delete the
       *  potential array when we exit this function and free the host potential_1d structure
       */
      density_1d* tp = new_density_1d_copy(total_potential);
      H->replace_internal_potential(tp->storage);
      /*Free the struct itself; not the array inside.*/
      free(tp);
      free(evals);
      if(last_evecs[N_eigs_to_check_for_convergence-1] != NULL) free(last_evecs[N_eigs_to_check_for_convergence-1]);
      for(unsigned eig=0; eig<N_eigs_to_check_for_convergence-1; eig++) {
        last_evecs[eig+1]=last_evecs[eig];
      }
      last_evecs[0] = evecs;
    }
    /*If we've hit this point, we hit the max number of iterations without success.*/
    for(unsigned eig=0; eig<N_eigs_to_check_for_convergence; eig++) {
      if(last_evecs[eig] != NULL) free(last_evecs[eig]);
    }
    free_density_1d(subband_density);
    if(total_potential != NULL) free_density_1d(total_potential);
    if(previous_potential != NULL) free_density_1d(previous_potential);
    free_density_1d(total_density);
    delete[] average_ooms;
    delete[] mstar;
    delete[] epsilon;
    return g_error_new(CALCULATION_SELF_CONSISTENT_GERROR_DOMAIN, 1, "kdotp::libkdotp::calculation_self_consistent_potential::perform_calculation: Maximum number of iterations (%u) reached without convergence.", max_iterations);
  }
  struct ::kdotp::libmodelxx::structure::generic_structure *generic_structure;
  hamiltonian* H;
  ooms_averager averager;
  wavefunction_max_delta_callback deltar;
  struct density_1d *dopant_density;
  /** Tolerances for GSL-based integration.*/
  double abserr;
  double relerr;
      };

    }
  }
}

#endif