gimme-alpha.c

Go to the documentation of this file.

/*
** gimme-alpha.c
**
** Made by (Johnny Q. Hacker)
** Login   <solarion@borkborkbork>
**

Copyright (C) 2008 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/>.

** Started on  Fri Sep 26 16:27:45 2008 Johnny Q. Hacker
** Last update Wed Jun 17 15:07:42 2009 Johnny Q. Hacker
*/

#include <stdio.h>
#include <stdlib.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <math.h>
#include <libcommon/gimme_alpha_types.h>
#include "process-args.h"
#include <libnextnano/read_nextnano_inputfile.h>
#include <libnextnano/read_nextnano_effectivemass_output.h>
#include <libnextnano/read_nextnano_bandstructure.h>
#include <libnextnano/read_nextnano_8x8kp_output.h>
#include <libcommon/gimme-alpha_structs.h>
#include <libpostproccalx/alpha_calc.h>
#include <libmodeling/wavefunction.h>

___malloc double* new_double(double val) {
  double* v = malloc(sizeof(double));
  if(unlikely(v == NULL)) return v;
  *v = val;
  return v;
             }

void point_hash_print(GHashTable *t) {
  double *v;
  if(t == NULL) {
    fprintf(stderr, "\tPoint hash does not exist!\n");
    return;
  }
  GList *keys = g_hash_table_get_keys(t);
  if(keys == NULL) {
    fprintf(stderr, "\tPoint hash is empty.\n");
    return;
  }
  for(GList *k = keys; k != NULL; k=g_list_next(k)) {
    fprintf(stderr, "\t%s=", (char*)(k->data));
    v = g_hash_table_lookup(t, k);
    if(v == NULL) {
      fprintf(stderr, "%s\n", "[NULL]");
    }else{
      fprintf(stderr, "%g\n", *v);
    }
  }
  g_list_free(keys);
             }

void print_value_double(gpointer data, gpointer private) {
  printf((char*)private, *(double*)data);
}

int main(int argc, char* argv[]) {
  g_type_init();
  GFile* pwd = g_file_new_for_commandline_arg(".");
  struct args *arg = process_args(argc, argv);
  if(arg == NULL) exit(2);
  GList *bdl;
  char *interror_title = "";
  if(!(arg->no_print_integration_error)) {
    interror_title="\tintegration error(eV m)";
  }
  if(arg->separate_terms) {
    printf("#directory\ttotal(eV m)\tinterface(eV m)\tefield(eV m)%s\n", interror_title);
  }
  for(bdl = arg->base_dirs; bdl != NULL; bdl = g_list_next(bdl)) {
    GFile *bd = bdl->data;
    GError *err = NULL;
    //fprintf(stderr, "back in gimme-alpha; about to get the gfile\n");
    GFile *keywords = g_file_get_child(bd, "keywords.in");
    //fprintf(stderr, "back in gimme-alpha; got the gfile\n");
    /*
      GFileInputStream *kwds_istream = g_file_read(keywords, NULL, &err);
      if(err != NULL) {
      fprintf(stderr, "Got error opening keywords file in %s: %s\n", arg->dir, err->message);
      exit(3);
      }
      GHashTable *kwds = read_nextnano_inputfile(G_INPUT_STREAM(db_istream), &err);
    */
    GFile *database = g_file_get_child(bd, "database.in");
    GFileInputStream *db_istream = g_file_read(database, NULL, &err);
    if(err != NULL) {
      gchar *fn = g_file_get_uri(bd);
      fprintf(stderr, "Got error opening database file %s: %s\n", fn, err->message);
      g_free(fn);
      exit(1);
    }
    GHashTable *db = read_nextnano_dbfile(G_INPUT_STREAM(db_istream), &err);
    /*
      fprintf(stderr, "Database:\n");
      print_stringy_inputfile(db);
    */
    GFile *input_file = g_file_get_child(bd, "input_file1.in");
    GFileInputStream *input_istream = g_file_read(input_file, NULL, &err);
    if(err != NULL) {
      gchar *fn = g_file_get_uri(bd);
      fprintf(stderr, "Got error opening input file %s: %s\n", fn, err->message);
      g_free(fn);
      exit(1);
    }
    GHashTable *ifile = read_nextnano_inputfile_stringy(G_INPUT_STREAM(input_istream), &err);
    struct mstar_wavefunction_1d *wf = read_effectivemass_complex_wavefunction2(CB, 1, 1, 1, NEU, 1, 1, bd, ifile, &err);
    if(unlikely((wf == NULL) || (err != NULL))) {
      fprintf(stderr, "Got error opening nextnano input file: %s\n", (err == NULL)?"[internal wavefunction conversion failure; is likely a change in nextnano wavefunction output format]":err->message);
      exit(1);
    }
    double *pot_x, *pot_band;
    int pot_n;
    pot_band = read_nextnano_bandstructure_file(POTENTIAL_FILE, bd, ifile, &err, &pot_x, &pot_n);
    struct potential pot = {.x = pot_x, .N = pot_n, .potential=pot_band};
    int hh_n;
    double *hh_x, *hh_band;
    hh_band = read_nextnano_bandstructure_file(HH_FILE, bd, ifile, &err, &hh_x, &hh_n);
    struct band hh_b = {.x = hh_x, .N = hh_n, .bandedge=hh_band};
    int lh_n;
    double *lh_x, *lh_band;
    lh_band = read_nextnano_bandstructure_file(LH_FILE, bd, ifile, &err, &lh_x, &lh_n);
    struct band lh_b = {.x = lh_x, .N = lh_n, .bandedge=lh_band};
    if(!(arg->avg_hh_lh & 1)){
      /*don't average the heavy and light holes; use just the heavy hole
  ("don't use the light hole")*/
      for(int l=0; l<lh_b.N; l++) {
  lh_b.bandedge[l] = hh_b.bandedge[l];
      }
    }
    if(!(arg->avg_hh_lh & 2)) {
      /*don't average the heavy and light holes; use just the light hole
  ("don't use the heavy hole")
      */
      for(int l=0; l<hh_b.N; l++) {
  hh_b.bandedge[l] = lh_b.bandedge[l];
      }
    }
    int so_n;
    double *so_x, *so_band;
    so_band = read_nextnano_bandstructure_file(SO_FILE, bd, ifile, &err, &so_x, &so_n);
    struct band so_b = {.x = so_x, .N = so_n, .bandedge=so_band};
    int cb_n;
    double *cb_x, *cb_band;
    cb_band = read_nextnano_bandstructure_file(CB_FILE, bd, ifile, &err, &cb_x, &cb_n);
    struct band cb_b = {.x = cb_x, .N = cb_n, .bandedge=cb_band};

    /*for(int i=0; i<cb_n; i++) printf("%g\t%g\n", cb_x[i], cb_band[i]);
      fprintf(stderr, "Interfaces found at:");
      for(GList *l=g_list_first(arg->interfaces); l!= NULL; l=g_list_next(l)) fprintf(stderr, "%g ", *((double*)(l->data)));
      fprintf(stderr, "\n");
    */

    double alpha_efield=0.0;
    double integration_error=0;
    double alpha_interface=0.0;
    double alpha=0.0;
    struct alpha_integrand_verificator vf;
    struct alpha_integrand_verificator_2 vf2;
    struct alpha_integrand_verificator *vfp = NULL;
    struct alpha_integrand_verificator_2 *vfp2 = NULL;
    if(arg->del_psi2) {
      if(arg->instrument_efield_alpha) vfp2 = &vf2;
      alpha = alpha_from_del_psi2_nostrain(wf, &hh_b, &lh_b, &so_b, &cb_b, &pot, vfp2, &integration_error, arg->request_integration_abs_error, arg->request_integration_rel_error, arg->no_integrate, arg->Ep);
    }else{
      if(!(arg->no_interface_term)) {
  /*Add all of the x points in the CB list if they're within an interface sequene.*/
  GList *l;
  /*Note that this is a "shallow copy" i.e. data points to the exact same thing.  Is OK; we just want to access it, not copy the data.  We'll free the list again when we're done (but, not the data pointed to by l->data!!)*/
        GList *ifaces = g_list_copy(arg->interfaces);
  GList *new_ifaces = NULL;
  for(int i=0; i<(cb_b.N-1); i++) {
    for(l=arg->ifaces_sequences; l!=NULL; l=l->next) {
      if((cb_b.x[i] >= ((double*)(l->data))[0]) && (cb_b.x[i] <= ((double*)(l->data))[1])) {
        /*Just like before, we only care about accessing the data; we don't really want to make a separate copy of it.*/
        //fprintf(stderr, "iface: %g-> (%g:%g)\n", cb_b.x[i], ((double*)(l->data))[0], ((double*)(l->data))[1]);
        new_ifaces = g_list_prepend(new_ifaces, new_double((cb_b.x[i] + cb_b.x[i])/2));
        ifaces = g_list_append(ifaces, new_ifaces->data);
      }else{
        //fprintf(stderr, "!iface: %g (%g:%g)\n", cb_b.x[i], ((double*)(l->data))[0], ((double*)(l->data))[1]);
      }
    }
  }
      alpha_interface=alpha_interface_term_mstar_nostrain(wf, &hh_b, &lh_b, &so_b, &cb_b, &pot, ifaces, arg->Ep);
      /*Make sure to free the interstitial storage*/
      for(GList *ll = new_ifaces; ll!=NULL; ll=ll->next) free(ll->data);
      g_list_free(new_ifaces);
      g_list_free(ifaces);
    }
    if(!(arg->no_efield_term)) {
      if(arg->instrument_efield_alpha) vfp = &vf;
      alpha_efield = alpha_electric_term_mstar_nostrain(wf, &hh_b, &lh_b, &so_b, &cb_b, &pot, vfp, &integration_error, arg->request_integration_abs_error, arg->request_integration_rel_error, arg->no_integrate, arg->Ep);
    }
  }
  if((vfp != NULL) || (vfp2 != NULL)) {
    fprintf(stderr, "Writing integration debug output.\n");
    GFile *instfile;
    if(vfp != NULL) {
      instfile = g_file_new_for_commandline_arg("./alpha_efield_instrumentation_data.pdata");
      alpha_write_verificator_data(instfile, vfp);
    }
    if(vfp2 != NULL) {
      instfile = g_file_new_for_commandline_arg("./alpha_method2_instrumentation_data.pdata");
      alpha_write_verificator2_data(instfile, vfp2);
    }
  }
  char interror_string[512];
  interror_string[0] = '\0';
  if(!(arg->no_print_integration_error)) {
    interror_title="\tintegration error(eV m)";
    snprintf(interror_string, 512, "\t%g", integration_error);
    interror_string[511] = '\0';
  }
  gchar* fn = g_file_get_uri(bd);
  if(arg->del_psi2) {
    printf("%s\t%g", fn, alpha);
  }else{
    printf("%s\t%g", fn, alpha_interface + alpha_efield);
  }
  g_free(fn);
  if((!(arg->del_psi2)) && (arg->separate_terms)) {
    printf("\t%g\t%g%s\n", alpha_interface, alpha_efield, interror_string);
  }else{
    printf("%s\n", interror_string);
  }
}
return 0;
}