oneoff_1d.c++

Go to the documentation of this file.

/*
 *A quick program to perform a 1D k-dot-p calc

 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 <string>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <glib.h>
#include <gio/gio.h>
#include <gsl/gsl_errno.h>
#include <libmodelxx/matdb.h++>
#include <libmodelxx/matdb-dotcode.h++>
#include <libmodelxx/continuous_structure.h++>
#include <libmodelxx/generic_structure.h++>
#include <libmodelxx/cartesian_grid.h++>
#include <libkdotp/hamiltonians.h++>
#include <libmodelxx/material.h++>
#include <libmodelxx/utilities.h++>
#include <libmetacalc/calculation_self_consistent_potential.h++>
#include <libkdotp/wavefunction_cartesian_effective_mass.h++>
#include <libpostproccalx/alpha_calc.h>
#include <libmodelxx/function_specification.h++>
#include "oneoff_1d_argprocess.h++"

using namespace kdotp::libmodelxx::continuous_structure;
using namespace kdotp::libmodelxx::structure;
using namespace kdotp::libkdotp::hamiltonian;
using namespace kdotp::libmodelxx::utilities;
using namespace kdotp::libmetacalc::self_consistent;
using namespace kdotp::libkdotp::wavefunction;
using namespace kdotp::libmodelxx::postprocessing;
using namespace std;

static void free_value(void* value, void* pridat) {free(value);}

void yyparse();
int yylex();
extern int yydebug;

extern struct continuous_structure* the_structure;
extern GInputStream* instream;

typedef kdotp::libmodelxx::grid::cartesian_grid<kdotp::libmodelxx::structure::material*,1,1> matgrid;
typedef kdotp::libmodelxx::grid::cartesian_grid<double,1,1> potgrid;

static void gsl_error_handler (const char * reason, const char * file, int line, int gsl_errno) {
  //fprintf(stderr, "GSL error occured in file %s, line %d (%d: %s).  This should be caught by the calling function, or report a bug.\n", file, line, gsl_errno, reason);
}


int main(int argc, char **argv) {
  g_type_init();
  gsl_set_error_handler(&gsl_error_handler);
  //yydebug = 1;
  struct args args;
  int err;
  err = process_args(argc, argv, &args);
  if(err != 0) exit(1);
  GString *mdb_file = g_string_new(args.matdb);
  struct matdb* mdb = read_matdb_dotcode(mdb_file, &err);
  if(mdb == NULL) {fprintf(stderr, "ERROR reading in the materials database %s: %d\n", args.matdb, err); return 1;}

  GList* l;
  GError* gerr = NULL;
  GString *gs = g_string_new("");
  gsize wlen;
  GHashTable *gst = new_symbol_table();
  GString *alpha_filename = g_string_new("");
  GFile* alpha_file;
  GFile* dump_file;
  for(l=(args.filestreams); l != NULL; l=l->next) {
    parameter_list_init(args.template_parameters);
    alpha_filename = g_string_assign(alpha_filename, ((struct file_info*)(l->data))->filename->str);
    alpha_filename = g_string_append(alpha_filename, "_alpha.pdata");
    alpha_file = g_file_new_for_path(alpha_filename->str);
    g_string_free(alpha_filename, TRUE);
    GOutputStream* alpha_gos = (GOutputStream*)g_file_replace(alpha_file, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerr);
    if(gerr != NULL) {
      fprintf(stderr, "ERROR opening alpha file (%s): %s\n", g_file_get_uri(alpha_file), gerr->message);
      exit(32);
    }
    g_string_printf(gs, "#");
    char* paramstring;
    g_string_append(gs, paramstring = get_header_for_template_state(args.template_parameters));
    g_free(paramstring);
    g_string_append_printf(gs, "\tEp\talpha_iface (eVm)\talpha_efield\talpha method 1(eVm)\talpha method2(eVm)\n");
    if(!g_output_stream_write_all(alpha_gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
      fprintf(stderr, "ERROR writing to alpha file (%s): %s\n", g_file_get_uri(alpha_file), gerr->message);
      exit(32);
    }
    /*Now do the loop*/
    do {
      gerr = NULL;
      add_template_parameters_to_symbol_table(gst, args.template_parameters);
      GString *gfn=g_string_new(((struct file_info*)(l->data))->filename->str);
      gfn = g_string_append(gfn, "_output/");
      GFile* job_directory = get_gfile_directory_for_template_state(gfn->str, args.template_parameters);
      g_string_free(gfn, TRUE);
      /*make the directory and cd into it*/
      g_file_make_directory_with_parents(job_directory, NULL, &gerr);
      /*Set the dump to go into this directory.*/
      dump_file = g_file_get_child(job_directory, "structure.griddump");
      if(gerr != NULL) {
  if(gerr->code != G_IO_ERROR_EXISTS) {
    char* fn = g_file_get_uri(job_directory);
    fprintf(stderr, "ERROR creating the target directory for templated job (%s): %s\n", fn, gerr->message);
    g_free(fn);
    exit(5);
  }else{
    /*If the directory already exists, that's OK.*/
    g_error_free(gerr);
    gerr = NULL;
  }
      }
      /*For now, we read in the file over and over.  It's inefficient, but not *that* inefficient.*/
      instream = (GInputStream*)((struct file_info*)(l->data))->instream;
      yyparse();
      if(the_structure == NULL) {
  fprintf(stderr, "the_structure was NULL!\n");
  exit(1);
      }
      struct generic_structure *s = postprocess_structure(the_structure, mdb, gst, &gerr);
      if(gerr != NULL) {
  fprintf(stderr, "ERROR while postprocessing structure: %s\n", gerr->message);
  return 5;
      }
      fprintf(stderr, "got structure!\n");

      gerr = generic_structure_dump(s, dump_file);
      if(gerr != NULL) {
  fprintf(stderr, "ERROR while dumping structure: %s\n", gerr->message);
  return 5;
      }
      fprintf(stderr, "got structure!\n");
      matgrid *mg = (matgrid*)(s->material_grid);
      potgrid *pg = (potgrid*)(s->potential_grid);
      int e;
      int L = mg->L[0];

      printf("Checking factorial: %u!->%lu %u^%u->%lu\n",
       13, compiletime_integer_factorial<13>(),
       4, 4, compiletime_integer_exponentiation<4>(4));


      GFile* file = g_file_get_child(job_directory, "bands.pdata");
      GOutputStream* gos = (GOutputStream*)g_file_replace(file, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerr);
      if(gerr != NULL) {
  fprintf(stderr, "ERROR opening bands file (%s): %s\n", g_file_get_uri(file), gerr->message);
  exit(32);
      }
      g_string_printf(gs, "#pos(nm)\tCB\tHH/LH\tSO\tPotential\n");
      if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
  fprintf(stderr, "ERROR writing to bands file (%s): %s\n", g_file_get_uri(file), gerr->message);
  exit(32);
      }
      double x;
      for(int i=0; i<L; i++) {
  e=0;
  double Eg = mg->data[i]->get_property("E_g_Gamma", &e);
  if(e != 0) fprintf(stderr, "ERROR: got error %d for property E_g_Gamma\n", e);
  e=0;
  double Ev = mg->data[i]->get_property("E_v", &e);
  if(e != 0) fprintf(stderr, "ERROR: got error %d for property E_v\n", e);
  e=0;
  double Delta = mg->data[i]->get_property("Delta", &e);
  if(e != 0) fprintf(stderr, "ERROR: got error %d for property Delta\n", e);
  unsigned X=i;
  mg->point_index_to_realspace(&X, &x);
  g_string_printf(gs, "%g\t%g\t%g\t%g\t%g\n", x, Eg+Ev, Ev, Ev-Delta, pg->data[i]);
  if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
    fprintf(stderr, "ERROR writing to bands file (%s): %s\n", g_file_get_uri(file), gerr->message);
    exit(32);
  }
      }
      g_object_unref(gos);

      double dopant_ierr, potential_ierr, efield_ierr;
      self_consistent_potential< ::kdotp::libkdotp::hamiltonian::cartesian_effective_mass, 1> scpot_calc(s, ::kdotp::libkdotp::hamiltonian::HARDWALL, get_average_with_cartesian_effective_mass_block, get_cartesian_effective_mass_block_max_delta, 1e6, 1e-14);

      double* eigvals;
      double* eigvecs;
      int neigs=0;
      double Ef;
      gerr = scpot_calc.perform_calculation(3, 1001, 1e-8, 1, &eigvecs, &eigvals, ::kdotp::libkdotp::hamiltonian::OOMSTAR, ::kdotp::libkdotp::hamiltonian::EPSILON, &Ef, &dopant_ierr, &potential_ierr, &efield_ierr, job_directory);
      printf("dopant integration error: %g\npotential integration error: %g\nefield integration error: %g\n", dopant_ierr, potential_ierr, efield_ierr);

      /*Verify our potential*/
      /*
      double* potcheck = new double[scpot_calc.H->L];
      scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
      for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "pc@%u: %g\n", i, potcheck[i]);
      delete[] potcheck;
      */


      /*
      potcheck = new double[scpot_calc.H->L];
      scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
      for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "pc@%u: %g\n", i, potcheck[i]);
      delete[] potcheck;
      */

      if(gerr != NULL) {
  fprintf(stderr, "ERROR getting eigenvectors of Hamiltonian: %s\n", gerr->message);
      }else{
  printf("Got %d eigenvalues: %g, %g, %g\n", neigs, eigvals[0], eigvals[1], eigvals[2]);
  char fn[1024];
  double sum, val;
  printf("Checking normalization and writing wavefunctions:\n");
  double x;
  for(int j=0; j<3; j++) {
    //sprintf(fn, "eigvec_%g.pdata", eigvals[j]);
    printf("\t%g: ", eigvals[j]);
    g_string_printf(gs, "eigvec_%d.pdata", j);
    file = g_file_get_child(job_directory, gs->str);
    GOutputStream* gos = (GOutputStream*)g_file_replace(file, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerr);
    if(gerr != NULL) {
      fprintf(stderr, "ERROR opening eigenvec %d file (%s): %s\n", j, g_file_get_uri(file), gerr->message);
      exit(32);
    }
    g_string_printf(gs, "#pos(nm)\tpsi\n");
    if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
      fprintf(stderr, "ERROR writing to eigvec %d file (%s): %s\n", j, g_file_get_uri(file), gerr->message);
      exit(32);
    }
    sum = 0;
    for(int i=0; i<L; i++) {
      val=eigvecs[i + j*L];
      unsigned X=i;
      scpot_calc.H->matgrid->point_index_to_realspace(&X, &x);
      g_string_printf(gs, "%g\t%g\n", x, val);
      if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
        fprintf(stderr, "ERROR writing to eigvec %d file (%s): %s\n", j, g_file_get_uri(file), gerr->message);
        exit(32);
      }
      sum += val*val;
    }
    if(fabs(1.0-sum) < 1e-8) {
      printf("Normalized (%g)\n", sum);
    }else{
      printf("NON-NORMALIZED (%g)\n", sum);
    }
    g_object_unref(gos);
  }

  /*
  potcheck = new double[scpot_calc.H->L];
  scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
  for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "pc@%u: %g\n", i, potcheck[i]);
  delete[] potcheck;
  */

  printf("################################################################\n");
  printf("################################################################\n");
  printf("## Final output\n");
  printf("################################################################\n");
  printf("################################################################\n");
  printf("Ef=%g\n", Ef);


  struct ::potential pot;
  pot.N = scpot_calc.H->L;
  pot.spacing=scpot_calc.H->dx;

  fprintf(stderr, "pot.x=%p pot.potential=%p H->internal_potential=%p\n", (void*)pot.x, (void*)pot.potential, (void*)scpot_calc.H->internal_potential);
  pot.x = (double*)malloc(sizeof(double)*pot.N);
  pot.potential = (double*)malloc(sizeof(double)*pot.N);
  if((pot.x == NULL) || (pot.potential == NULL)) {
    fprintf(stderr, "ERROR: unable to allocate memory for pot x or bandedge data!\n");
    exit(5);
  }
  /*
  potcheck = new double[scpot_calc.H->L];
  scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
  for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "apc@%u: %g\n", i, potcheck[i]);
  delete[] potcheck;
  */

  fprintf(stderr, "N=%d dx=%g\n", pot.N, pot.spacing);

  for(unsigned int i=0; i<(unsigned int)pot.N; i++) {
    scpot_calc.H->matgrid->point_index_to_realspace(&i, &(pot.x[i]));
    fprintf(stderr, "%u/%d: x=%g (&X=%p, &pot.x[%d]=%p(%p/%p))\n", X, i, pot.x[i], (void*)&i, i, (void*)&(pot.x[i]), (void*)pot.x, (void*)scpot_calc.H->internal_potential);
  }
  /*
  potcheck = new double[scpot_calc.H->L];
  scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
  for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "bpc@%u: %g\n", i, potcheck[i]);
  delete[] potcheck;
  */

  scpot_calc.H->get_parameter_function(TOTAL_POTENTIAL, pot.potential);

  /*
  potcheck = new double[scpot_calc.H->L];
  scpot_calc.H->get_parameter_function(::kdotp::libkdotp::hamiltonian::TOTAL_POTENTIAL, potcheck);
  for(unsigned i=0; i<scpot_calc.H->L; i++) fprintf(stderr, "cpc@%u: %g\n", i, potcheck[i]);
  delete[] potcheck;
  */

  struct mstar_wavefunction_1d *wf = mstar_wavefunction_1d_new(scpot_calc.H->L);

  if(wf == NULL) {
    fprintf(stderr, "ERROR: unable to allocate memory for wavefunction data!\n");
    exit(5);
  }
  for(int i=0; i<wf->N; i++) {
    unsigned int X=i;
    scpot_calc.H->matgrid->point_index_to_realspace(&X, &(wf->points[i].x));
    scpot_calc.H->matgrid->point_index_to_realspace(&X, &(wf->points[i].x));
    wf->points[i].im=0;
    /*We want the first wavefunction*/
    wf->points[i].re=eigvecs[i];
  }

  /**Bands: 0->CB, 1->HH, 2->SO*/
  struct band bands[3];
  int errnum;
  for(int i=0; i<3; i++) {
    bands[i].N = scpot_calc.H->L;
    bands[i].spacing=scpot_calc.H->dx;
    bands[i].x = (double*)malloc(sizeof(double)*bands[i].N);
    bands[i].bandedge = (double*)malloc(sizeof(double)*bands[i].N);
    if((bands[i].x == NULL) || (bands[i].bandedge == NULL)) {
      fprintf(stderr, "ERROR: unable to allocate memory for pot x or bandedge data!\n");
      exit(5);
    }
    for(unsigned int j=0; j<(unsigned int)bands[i].N; j++) {
      scpot_calc.H->matgrid->point_index_to_realspace(&j, &(bands[i].x[j]));
    }
  }
  /*For alpha_calch.h compat, and implying nextnano compat, we add in the potential*/
  /*\note these band edge energies are relative to Ef, for compatibility with alpha_calc (and therefore with nextnano*/
  scpot_calc.H->get_parameter_function(EC, bands[0].bandedge);
  for(unsigned j=0; j<scpot_calc.H->L; j++) {
    /*Add the potential to the CB*/
    bands[0].bandedge[j] += -(pot.potential[j] + Ef);
    /*Set the HH/LH band edge*/
    bands[1].bandedge[j] = scpot_calc.H->matgrid->data[j]->get_property("E_v", &errnum) - pot.potential[j] - Ef;
    if(errnum != 0) fprintf(stderr, "ERROR getting property Ev at gridsite %u: %d\n", j, errnum);
    /*Subtract the Delta from the HH/LH band edge (note: HH/LH includes potential and Ef already)*/
    bands[2].bandedge[j] = bands[1].bandedge[j] - scpot_calc.H->matgrid->data[j]->get_property("Delta", &errnum);
    if(errnum != 0) fprintf(stderr, "ERROR getting property Delta at gridsite %u: %d\n", j, errnum);
    /**\note that the bandedge of the HH/LH band includes the potential (previous loop)*/
  }
  /*Use method 1*/
  /*First, find the ifaces.*/
  GList* ifaces=NULL;
  /*We start at 1, since we compare against the *previous* point*/
  /*While we're going through it, get Ep for the material.*/
  double* Ep_in_structure = new double[scpot_calc.H->L];
  for(unsigned i=1; i<scpot_calc.H->L; i++) {
    if(scpot_calc.H->matgrid->data[i-1] != scpot_calc.H->matgrid->data[i]) {
      double *midpoint = new double;
      double x_prev, x_cur;
      unsigned prev=i-1;
      scpot_calc.H->matgrid->point_index_to_realspace(&prev, &x_prev);
      scpot_calc.H->matgrid->point_index_to_realspace(&i, &x_cur);
      /*cur and prev are fixed; the interface is halfway between.*/
      *midpoint = x_prev + (x_cur - x_prev)*0.5;
      ifaces = g_list_append(ifaces, midpoint);
      printf("Got interface: %g\n", *midpoint);
    }
    Ep_in_structure[i] = scpot_calc.H->matgrid->data[i]->get_property("E_p", &errnum);
    if(errnum != 0) fprintf(stderr, "ERROR getting property E_p at gridsite %u: %d\n", i, errnum);
  }
  delete[] Ep_in_structure;

  struct alpha_integrand_verificator vf1;
  struct alpha_integrand_verificator_2 vf2;

  double Ep = get_average_with_cartesian_effective_mass_block(eigvecs, 0, scpot_calc.H->L, Ep_in_structure);

  printf("Ep_avg=%g\n", Ep);

  double ierr_alpha, ierr_alpha_m2;
  double alpha_iface = alpha_interface_term_mstar_nostrain(wf, &bands[1], &bands[1], &bands[2], &bands[0], &pot, ifaces, Ep);
  g_list_foreach(ifaces, free_value, NULL);
  g_list_free(ifaces);
  double alpha_efield = alpha_electric_term_mstar_nostrain(wf, &bands[1], &bands[1], &bands[2], &bands[0], &pot, &vf1, &ierr_alpha, 1e-7, 1e-7, FALSE, Ep);
  double alpha_method2 = alpha_from_del_psi2_nostrain(wf, &bands[1], &bands[1], &bands[2], &bands[0], &pot, &vf2, &ierr_alpha_m2, 1e-7, 1e-7, FALSE, Ep);

  /*Write out the band ifnormation as it'll be passed into alpah_calc.*/
  file = g_file_get_child(job_directory, "bands2.pdata");
  GOutputStream* gos = (GOutputStream*)g_file_replace(file, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerr);
  if(gerr != NULL) {
    fprintf(stderr, "ERROR opening bands2 file (%s): %s\n", g_file_get_uri(file), gerr->message);
    exit(32);
  }
  g_string_printf(gs, "#x(nm)\tcb(eV)\tHH/LH(eV)\tSO(eV)\tpot(eV)\tpsi^2(nm^-1)\n");
  if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
    fprintf(stderr, "ERROR writing to bands2 file (%s): %s\n", g_file_get_uri(file), gerr->message);
    exit(32);
  }
  for(int i=0; i<bands[0].N; i++) {
    g_string_printf(gs, "%g\t", bands[0].x[i]);
    for(int j=0; j<3; j++) {
      g_string_append_printf(gs, "%g\t", bands[j].bandedge[i]);
    }
    g_string_append_printf(gs, "%g\n", pot.potential[i]);
    g_string_append_printf(gs, "%g\n", wf->points[i].re*wf->points[i].re);
    g_string_append_printf(gs, "#x(nm)\tcb(eV)\tHH/LH(eV)\tSO(eV)\tpot(eV)\tpsi^2(nm^-1)\n");
    if(!g_output_stream_write_all(gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
      fprintf(stderr, "ERROR writing to bands2 file (%s): %s\n", g_file_get_uri(file), gerr->message);
      exit(32);
    }
  }
  g_object_unref(gos);

  gs = g_string_assign(gs, paramstring=get_string_for_template_state(args.template_parameters));
  g_free(paramstring);
  g_string_append_printf(gs, "\t%g\t%g\t%g\t%g\t%g\n", Ep, alpha_iface, alpha_efield, alpha_iface + alpha_efield, alpha_method2);
  if(!g_output_stream_write_all(alpha_gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
    fprintf(stderr, "ERROR writing to alpha file (%s): %s\n", g_file_get_uri(alpha_file), gerr->message);
    exit(32);
  }

  GFile* vf1_file = g_file_get_child(job_directory, "alpha_method_1_efield_instrumentation.pdata");
  GFile* vf2_file = g_file_get_child(job_directory, "alpha_method_2_instrumentation.pdata");
  alpha_write_verificator_data(vf1_file, &vf1);
  alpha_write_verificator2_data(vf2_file, &vf2);
  g_object_unref(vf1_file);
  g_object_unref(vf2_file);
  for(int i=0; i<3; i++) {
    free(bands[i].x);
    free(bands[i].bandedge);
  }
  free(pot.potential);
  free(pot.x);
  mstar_wavefunction_1d_free(wf);
      }
    }while(!update_parameters_plus_check_if_done(args.template_parameters));
    g_object_unref(alpha_gos);
    g_object_unref(alpha_file);
    g_object_unref(dump_file);
  }
  return 0;
}