oneoff_optimize_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/>.

*/

/*The following line will be removed in the auto-generated MPI variant.*/
#undef HAVE_MPI

#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/grid.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 <libmetacalc/digital_alloy_optimizer_1d.h++>
#include <libmetacalc/continuous_alloy_optimizer_1d.h++>
#ifdef HAVE_MPI
#include <mpi.h>
#include <libmetacalc/digital_alloy_optimizer_1d_mpi.h++>
#endif
#include "oneoff_optimize_1d_argprocess.h++"

using namespace kdotp::libmodelxx::continuous_structure;
using namespace kdotp::libmodelxx::structure;
using namespace kdotp::libmodelxx::grid;
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;
using namespace kdotp::libmetacalc::optimizer;

static void free_value(void* value, void* pridat) {free(value);}
static void delete_double(void* value, void* pridat) {delete (double*)value;}
void set_point_callback(unsigned X, const material* mat, void* privdat);
double alpha_efield_change_calculation_callback(void* privdat, unsigned assignment, long unsigned iteration);
double calculate_alpha_for_structure(void* privdat, unsigned assignment, long unsigned iteration);
void write_raw_bands_file(GFile* file, struct generic_structure* genstruct);

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;

struct set_point_callback_privdat {
  struct generic_structure *genstruct;
  self_consistent_potential< ::kdotp::libkdotp::hamiltonian::cartesian_effective_mass, 1>* scpot_calc;
  GFile* base_job_directory;
  GOutputStream* alpha_by_iteration_gos;
  GFile* alpha_by_iteration_file;
  /*This stores the different alpha values from the current computation.*/
  double alpha_iface;
  double alpha_efield;
  double alpha_method2;
  double alpha_iface_0field;
  double alpha_efield_0field;
  double alpha_method2_0field;
  double Ep;
  gboolean do_nofield_comparison;
  #ifdef HAVE_MPI
  int N_nodes;
  int node;
  #endif
};


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();
  #ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
  #endif
  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);
  fprintf(stderr, "Got minx=%g maxx=%g stepsize=%g multiplier=%g min_stepsize=%g\n", args.minx, args.maxx, args.stepsize, args.multiplier, args.min_stepsize);
  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); 
    #ifdef HAVE_MPI
    MPI_Finalize();
    #endif
    return 1;
  }

  /*Insert the fake material into the mdb.*/
  ::kdotp::libmodelxx::structure::insert_external_material(mdb, args.optimizer_material);

  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;
  #ifdef HAVE_MPI
  int node;
  int N_nodes;
  MPI_Comm_rank(MPI_COMM_WORLD, &node);
  MPI_Comm_size(MPI_COMM_WORLD, &N_nodes);
  #endif
  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, "\tstartstruct\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);
    }
    gerr = NULL;
    GString *gfn=g_string_new(((struct file_info*)(l->data))->filename->str);
    gfn = g_string_append(gfn, "_output/");
    do {
      add_template_parameters_to_symbol_table(gst, args.template_parameters);
      GFile* file_directory = get_gfile_directory_for_template_state(gfn->str, args.template_parameters);
      g_file_make_directory_with_parents(file_directory, NULL, &gerr);
      if(gerr != NULL) {
  if(gerr->code != G_IO_ERROR_EXISTS) {
    char* fn = g_file_get_uri(file_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(unsigned sstruct=0; sstruct<args.N_starting_points; sstruct++) {
  /*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);
          #ifdef HAVE_MPI
    MPI_Finalize();
          #endif
    return 5;
  }
  fprintf(stderr, "got structure!\n");

  material* markermat = new material(mdb, args.optimizer_material, NULL, 1);
  material* mat1=NULL;
  material* mat0=NULL;
  if(!(args.continuous)) {
    mat1 = new material(mdb, args.mat1.mat1, args.mat1.mat2, args.mat1.x);
    mat0 = new material(mdb, args.mat0.mat1, args.mat0.mat2, args.mat0.x);
  }

  GError *e = NULL;
  /**We have the structure set up.  Now we need to set up the optimizer.
   *\note that this is used for both set_point and calculate_alpha_for_structure.
   */
  GString *iter_dirname = g_string_new("");
  g_string_append_printf(iter_dirname, "%s_optimizer_start_%u", (args.continuous)?"continuous":"digital", sstruct);
  GFile* job_directory = g_file_get_child(file_directory, iter_dirname->str);
  g_string_free(iter_dirname, TRUE);
  g_file_make_directory_with_parents(job_directory, NULL, &gerr);
  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;
    }
  }
  GString *alpha_by_iter_filename = g_string_new("alpha_sstruct");
  #ifdef HAVE_MPI
  g_string_append_printf(alpha_by_iter_filename, "_node%d", node);
  #endif
  g_string_append_printf(alpha_by_iter_filename, "=%u.pdata", sstruct);
  GFile* alpha_by_iter_file = g_file_get_child(job_directory, alpha_by_iter_filename->str);
  g_string_free(alpha_by_iter_filename, TRUE);
  GOutputStream* alpha_by_iter_gos = (GOutputStream*)g_file_replace(alpha_by_iter_file, NULL, FALSE, G_FILE_CREATE_NONE, NULL, &gerr);
  if(gerr != NULL) {
    fprintf(stderr, "ERROR opening alpha_by_iter file (%s): %s\n", g_file_get_uri(alpha_by_iter_file), gerr->message);
    exit(32);
  }
  g_string_printf(gs, "#");
  g_string_append(gs, paramstring = get_header_for_template_state(args.template_parameters));
  g_free(paramstring);
  g_string_append_printf(gs, "iteration\ttotal\tiface\tefield");
  if(args.do_nofield_comparison) {
    g_string_append_printf(gs, "total 0field\tiface 0field\tefield 0field\ttotal difference\tiface difference\tefield difference\n");
  }else{
    gs = g_string_append(gs, "\n");
  }
  if(!g_output_stream_write_all(alpha_by_iter_gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
    fprintf(stderr, "ERROR writing to alpha_by_iteration file (%s): %s\n", g_file_get_uri(alpha_by_iter_file), gerr->message);
    exit(32);
  }

  GFile* dump_file = g_file_get_child(job_directory, "structure.griddump");
  gerr = generic_structure_dump(s, dump_file);
  if(gerr != NULL) {
    fprintf(stderr, "ERROR while dumping structure: %s\n", gerr->message);
          #ifdef HAVE_MPI
    MPI_Finalize();
          #endif
    return 5;
  }
  fprintf(stderr, "got structure!\n");
  struct set_point_callback_privdat spcp = {
    s,
    NULL,
    job_directory,
    alpha_by_iter_gos,
    alpha_by_iter_file,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    args.do_nofield_comparison
    #ifdef HAVE_MPI
    N_nodes,
    node,
    #endif
  };
        #ifdef HAVE_MPI
  digital_alloy_optimizer_1d_mpi* digital_optimizer_mpi=NULL;
        #endif
  digital_alloy_optimizer_1d* optimizer=NULL;
  continuous_alloy_optimizer_1d* continuous_optimizer=NULL;
  if(args.continuous) {
    continuous_optimizer = new continuous_alloy_optimizer_1d(markermat, args.mat1.mat1, args.mat0.mat1, s, &set_point_callback, args.blocksize, continuous_alloy_optimizer_1d::MAXIMIZE_ABS, &e, (void*)&spcp, mdb, args.minx, args.maxx, args.seed_value, args.stepsize, args.N_blocks, args.assignments, args.assignment_expr);
  }else{
          #ifdef HAVE_MPI
    if(args.mpi) {
      digital_optimizer_mpi = new digital_alloy_optimizer_1d_mpi(markermat, mat1, mat0, s, &set_point_callback, args.blocksize, digital_alloy_optimizer_1d_mpi::MAXIMIZE_ABS, &e, (void*)&spcp, N_nodes, node, args.N_blocks, args.seed_value, args.assignments, args.assignment_expr, args.temperature);
    }else{
            #endif
      optimizer = new digital_alloy_optimizer_1d(markermat, mat1, mat0, s, &set_point_callback, args.blocksize, digital_alloy_optimizer_1d::MAXIMIZE_ABS, &e, (void*)&spcp, args.N_blocks, args.seed_value, args.assignments, args.assignment_expr);
      #ifdef HAVE_MPI
    }
            #endif
  }
  if(e != NULL) {
    fprintf(stderr, "ERROR setting up optimizer: %s", e->message);
          #ifdef HAVE_MPI
    MPI_Finalize();
          #endif
    return 102;
  }
  if(args.continuous) {
    e = continuous_optimizer->do_optimization(args.do_nofield_comparison?&alpha_efield_change_calculation_callback:&calculate_alpha_for_structure, (void*)&spcp, (void*)&spcp, args.multiplier, args.min_stepsize);
  }else{
    #ifdef HAVE_MPI
    if(args.mpi) {
      e = digital_optimizer_mpi->do_optimization(args.do_nofield_comparison?&alpha_efield_change_calculation_callback:&calculate_alpha_for_structure, (void*)&spcp, (void*)&spcp, args.anneal_start, args.anneal_end, args.num_anneal_steps, args.anneal_cycles, args.calculation_steps_per_annealing_step);
    }else{
      #endif
      e = optimizer->do_optimization(args.do_nofield_comparison?&alpha_efield_change_calculation_callback:&calculate_alpha_for_structure, (void*)&spcp, (void*)&spcp);
      #ifdef HAVE_MPI
    }
            #endif
  }
  if(e != NULL) {
    fprintf(stderr, "ERROR running optimizer: %s", e->message);
          #ifdef HAVE_MPI
    MPI_Finalize();
          #endif
    return 103;
  }
  gs = g_string_assign(gs, paramstring=get_string_for_template_state(args.template_parameters));
  g_free(paramstring);
  g_string_append_printf(gs, "\t%u\t%g\t%g\t%g\t%g\t%g", sstruct, spcp.Ep, spcp.alpha_iface, spcp.alpha_efield, spcp.alpha_iface + spcp.alpha_efield, spcp.alpha_method2);
  if(args.do_nofield_comparison) {
    g_string_append_printf(gs, "\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", spcp.alpha_iface_0field, spcp.alpha_efield_0field, spcp.alpha_iface_0field + spcp.alpha_efield_0field, spcp.alpha_method2, spcp.alpha_iface - spcp.alpha_iface_0field, spcp.alpha_efield - spcp.alpha_efield_0field, (spcp.alpha_iface - spcp.alpha_iface_0field) + (spcp.alpha_efield - spcp.alpha_efield_0field));
  }else{
    gs = g_string_append(gs, "\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);
  }
  dump_file = g_file_get_child(job_directory, "final-structure.griddump");
  gerr = generic_structure_dump(s, dump_file);
  if(gerr != NULL) {
    fprintf(stderr, "ERROR while dumping structure: %s\n", gerr->message);
          #ifdef HAVE_MPI
    MPI_Finalize();
          #endif
    return 5;
  }
  GFile* file = g_file_get_child(job_directory, "final-bands.pdata");
  write_raw_bands_file(file, spcp.genstruct);
  if(args.continuous) {
    delete continuous_optimizer;
  }else{
    delete optimizer;
  }
  g_object_unref(alpha_by_iter_gos);
  g_object_unref(alpha_by_iter_file);
  generic_structure_free(s);
      }
      g_object_unref(file_directory);
    }while(!update_parameters_plus_check_if_done(args.template_parameters));
    g_string_free(gfn, TRUE);
    g_object_unref(alpha_gos);
    g_object_unref(alpha_file);
  }
  #ifdef HAVE_MPI
  MPI_Finalize();
  #endif
  return 0;
}

void set_point_callback(unsigned X, const material* mat, void* privdat) {
  struct set_point_callback_privdat* spcp = (struct set_point_callback_privdat*)privdat;
  /*We don't need to do anything if we've not yet set up the calculation (and therefore also not the hamiltonian).*/
  if(spcp->scpot_calc == NULL) return;
  /**\note that it'll likely be faster to introduce a special routine to just change the nearest neighbors of this point.*/
}

double alpha_efield_change_calculation_callback(void* privdat, unsigned assignment, long unsigned iteration) {
  struct set_point_callback_privdat* spcp = (struct set_point_callback_privdat*)privdat;
  /*We set the job directory to "with_field/" and "without_field/" for each respective run
    Additionally, we save off the alpha values that are being overwritten.
    Technically, Ep should change (it's an average), but not a lot.
   */
  GFile* original_base_job_dir = spcp->base_job_directory;

  GError* e = NULL;
  void* original_applied_potential = spcp->genstruct->potential_grid;
  if(e != NULL) {
    fprintf(stderr, "ERROR making 0-potential grid: %s\n", e->message);
    g_error_free(e);
    exit(64);
  }
  spcp->genstruct->potential_grid = set_up_grid<double>(spcp->genstruct->g, 0, &e);
  spcp->base_job_directory = g_file_get_child(original_base_job_dir, "without_field");
  calculate_alpha_for_structure(privdat, assignment, iteration);

  generic_structure_free_potential(spcp->genstruct);
  spcp->genstruct->potential_grid = original_applied_potential;
  spcp->alpha_iface_0field = spcp->alpha_iface;
  spcp->alpha_efield_0field = spcp->alpha_efield;
  spcp->alpha_method2_0field = spcp->alpha_method2;
  g_object_unref(spcp->base_job_directory);
  spcp->base_job_directory = g_file_get_child(original_base_job_dir, "with_field");
  calculate_alpha_for_structure(privdat, assignment, iteration);

  g_object_unref(spcp->base_job_directory);
  spcp->base_job_directory = original_base_job_dir;
  /*Now write out the values we got*/
  GString* gs = g_string_new("");
  g_string_append_printf(gs, "%lu\t%g\t%g\t%g", iteration, spcp->alpha_iface + spcp->alpha_efield, spcp->alpha_iface, spcp->alpha_efield);
  g_string_append_printf(gs, "%g\t%g\t%g", spcp->alpha_iface_0field + spcp->alpha_efield_0field, spcp->alpha_iface_0field, spcp->alpha_efield_0field);
  g_string_append_printf(gs, "%g\t%g\t%g", (spcp->alpha_iface - spcp->alpha_iface_0field) + (spcp->alpha_efield - spcp->alpha_efield_0field), spcp->alpha_iface - spcp->alpha_iface_0field, spcp->alpha_efield - spcp->alpha_efield_0field);
  gsize wlen;
  GError *gerr = NULL;
  if(!g_output_stream_write_all(spcp->alpha_by_iteration_gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
    fprintf(stderr, "ERROR writing to alpha_by_iteration file (%s): %s\n", g_file_get_uri(spcp->alpha_by_iteration_file), gerr->message);
    exit(32);
  }
  /*Stopping memory leaks*/
  g_string_free(gs, TRUE);
  /*And return the *difference**/
  return (spcp->alpha_iface - spcp->alpha_iface_0field) + (spcp->alpha_efield - spcp->alpha_efield_0field);
}

double calculate_alpha_for_structure(void* privdat, unsigned assignment, long unsigned iteration) {
  struct set_point_callback_privdat* spcp = (struct set_point_callback_privdat*)privdat;
  GError *gerr = NULL;
  matgrid *mg = (matgrid*)(spcp->genstruct->material_grid);
  potgrid *pg = (potgrid*)(spcp->genstruct->potential_grid);
  int L = mg->L[0];
  GString *gs=g_string_new("");
  GString *iter_dirname = g_string_new("optimizer_iteration_");
  g_string_append_printf(iter_dirname, "%lu", iteration);
  #ifdef HAVE_MPI
  g_string_append_printf(iter_dirname, "_node_%d_assignment_%u", spcp->node, assignment);
  #endif
  GFile* job_directory = g_file_get_child(spcp->base_job_directory, iter_dirname->str);
  g_string_free(iter_dirname, TRUE);
  /*make the directory and cd into it*/
  g_file_make_directory_with_parents(job_directory, NULL, &gerr);
  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;
    }
  }
  /*Set the dump to go into this directory.*/
  GFile* dump_file = g_file_get_child(job_directory, "structure.griddump");
  GFile* file = g_file_get_child(job_directory, "bands.pdata");
  write_raw_bands_file(file, spcp->genstruct);
  g_object_unref(file);
  file = NULL;
  gsize wlen;

  double dopant_ierr, potential_ierr, efield_ierr;

  if(spcp->scpot_calc == NULL) {
    /*The goal is to make this not get regenerated on each loop.  The main time-eater, though, is in the potential calculation, so it's not a huge issue*/
    spcp->scpot_calc = new self_consistent_potential< ::kdotp::libkdotp::hamiltonian::cartesian_effective_mass, 1>(spcp->genstruct, ::kdotp::libkdotp::hamiltonian::HARDWALL, get_average_with_cartesian_effective_mass_block, get_cartesian_effective_mass_block_max_delta, 1e6, 1e-14);
  }

  double* eigvals=NULL;
  double* eigvecs=NULL;
  int neigs=0;
  double Ef=0;
  gerr = (spcp->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, 0.75);
  printf("dopant integration error: %g\npotential integration error: %g\nefield integration error: %g\n", dopant_ierr, potential_ierr, efield_ierr);

  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;
  (spcp->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;
      }
      g_object_unref(file);
      if(fabs(1.0-sum) < 1e-8) {
  printf("Normalized (%g)\n", sum);
      }else{
  printf("NON-NORMALIZED (%g)\n", sum);
      }
      g_object_unref(gos);
    }
    printf("################################################################\n");
    printf("################################################################\n");
    printf("## Final output\n");
    printf("################################################################\n");
    printf("################################################################\n");
    printf("Ef=%g\n", Ef);
    struct ::potential pot;
    pot.N = (spcp->scpot_calc)->H->L;
    pot.spacing=(spcp->scpot_calc)->H->dx;

    fprintf(stderr, "pot.x=%p pot.potential=%p H->internal_potential=%p\n", (void*)pot.x, (void*)pot.potential, (void*)(spcp->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);
    }
    fprintf(stderr, "N=%d dx=%g\n", pot.N, pot.spacing);
    for(unsigned int i=0; i<(unsigned int)pot.N; i++) {
      (spcp->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*)(spcp->scpot_calc)->H->internal_potential);
    }
    (spcp->scpot_calc)->H->get_parameter_function(TOTAL_POTENTIAL, pot.potential);
    struct mstar_wavefunction_1d *wf = mstar_wavefunction_1d_new((spcp->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;
      (spcp->scpot_calc)->H->matgrid->point_index_to_realspace(&X, &(wf->points[i].x));
      (spcp->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 = (spcp->scpot_calc)->H->L;
      bands[i].spacing=(spcp->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++) {
  (spcp->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*/
    (spcp->scpot_calc)->H->get_parameter_function(EC, bands[0].bandedge);
    for(unsigned j=0; j<(spcp->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] = (spcp->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] - (spcp->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[(spcp->scpot_calc)->H->L];
    for(unsigned i=1; i<(spcp->scpot_calc)->H->L; i++) {
      if((spcp->scpot_calc)->H->matgrid->data[i-1] != (spcp->scpot_calc)->H->matgrid->data[i]) {
  double *midpoint = new double;
  double x_prev, x_cur;
  unsigned prev=i-1;
  (spcp->scpot_calc)->H->matgrid->point_index_to_realspace(&prev, &x_prev);
  (spcp->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] = (spcp->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, (spcp->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, delete_double, 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(file);
    g_object_unref(gos);
    spcp->Ep = Ep;
    spcp->alpha_iface = alpha_iface;
    spcp->alpha_efield = alpha_efield;
    spcp->alpha_method2 = alpha_method2;

    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);
    alpha_free_verificator_data(&vf1);
    alpha_free_verificator2_data(&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);
    free(eigvecs);
    free(eigvals);
  }
  /*This should go away after I've vetted the spcot class for this use.*/
  delete spcp->scpot_calc;
  spcp->scpot_calc = NULL;
  g_object_unref(dump_file);
  /*If do_nofield_comparison is true, we're being called from another callback, which is writing the values out to the file.*/
  if(!spcp->do_nofield_comparison) {
    gs = g_string_assign(gs, "");
    g_string_append_printf(gs, "%lu\t%g\t%g\t%g\n", iteration, spcp->alpha_iface + spcp->alpha_efield, spcp->alpha_iface, spcp->alpha_efield);
    if(!g_output_stream_write_all(spcp->alpha_by_iteration_gos, gs->str, gs->len, &wlen, NULL, &gerr)) {
      fprintf(stderr, "ERROR writing to alpha_by_iteration file (%s): %s\n", g_file_get_uri(spcp->alpha_by_iteration_file), gerr->message);
      exit(32);
    }
    /*Stopping memory leaks*/
    g_string_free(gs, TRUE);
  }
  return spcp->alpha_iface + spcp->alpha_efield;
}

void write_raw_bands_file(GFile* file, struct generic_structure* genstruct) {
  matgrid *mg = (matgrid*)(genstruct->material_grid);
  potgrid *pg = (potgrid*)(genstruct->potential_grid);
  GString *gs = g_string_new("");
  gsize wlen;
  int L = mg->L[0];
  GError *gerr = NULL;
  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;
  int e;
  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);
}