lp_wavefunction.c

Go to the documentation of this file.

/*
** wavefunction.c
** 

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

** Made by (Johnny Q. Hacker)
** Login   <solarion@nathan>
** 
** Started on  Tue Nov 17 22:28:22 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <stdlib.h>
#include <stdio.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <math.h>
#include <libpostproc/lp_wavefunction.h>
#include <libdotcode/binary_data_common.h>
#include <libpostproc/unit_conversion.h>
#include "gcc_hints.h"

#ifdef OPENMP
#include <omp.h>
#endif

#define NUM_STEPS 12
#define THRESHOLD 0.5

#define WAVEFUNC_GERROR_DOMAIN 4096

#define MAX_PIXEL_VAL 255

#define MAX_GSL_INTEGRATION_INTERVALS 4096

/**
 *Gets the index into the storage member for a wavefunction, or for a density
 *\param x x index component
 *\param y y index component
 *\param z z index component
 *\param b band index component
 *\param im 1 if imaginary component is wanted; 0 else.
 *\param Nx Number of x sites
 *\param Ny Number of y sites
 *\param Nb Number of bands
 *\returns long unsigned int, which is the address in the storage
 *array.
 *\note there is no error checking, for speed purposes. E.g. if im is
 *3, it'll just be blindly used as a 3 and you'll be addressing the
 *imaginary component of the next index.
 *\note Nz is not needed.
 */
hotfunc ___const ___always_inline long unsigned int get_wavefunction_storage_index(unsigned int x, unsigned int y, unsigned int z, unsigned int im, unsigned int b, unsigned int Nx, unsigned int Ny, int Nb) {
  return im + 2*(b + Nb*(x + Nx*(y + Ny*(z))));
             }

hotfunc ___const ___always_inline long unsigned int get_density_storage_index(unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny) {
  return x + Nx*(y + Ny*(z));
}
hotfunc ___const ___always_inline long unsigned int get_density_slice_storage_index(enum coord c, unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny) {
  switch(c) {
  case X:
    x=0;
    Nx=1;
    break;
  case Y:
    y=0;
    Ny=1;
    break;
  case Z:
    z=0;
    break;
  }
  return x + Nx*(y + Ny*(z));
}

long unsigned int get_density_storage_index_(unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny) {
  return get_density_storage_index(x, y, z, Nx, Ny);
}

hotfunc ___always_inline double* get_density_storage_at(unsigned x, unsigned y, unsigned z, struct density* d) {
  return &(d->storage[get_density_storage_index(x, y, z, d->N[0], d->N[1])]);
}

double* get_density_storage_at_(unsigned x, unsigned y, unsigned z, struct density* d) {
  return get_density_storage_at(x, y, z, d);
}

hotfunc ___always_inline const double* get_density_storage_at_c(unsigned x, unsigned y, unsigned z, const struct density* d) {
  return &(d->storage[get_density_storage_index(x, y, z, d->N[0], d->N[1])]);
}

const double* get_density_storage_at_c_(unsigned x, unsigned y, unsigned z, const struct density* d) {
  return get_density_storage_at_c(x, y, z, d);
}



void init_density(struct density* d, double value) {
  int i, j, k;
#ifdef DEBUG2
  fprintf(stderr, "wavefunc::init_density: d->N={%d,%d,%d} d->dx={%g,%g,%g}\n", d->N[0], d->N[1], d->N[2], d->dx[0], d->dx[1], d->dx[2]);
#endif
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] = value;
      }
    }
  }
}

/*This is where the complex wavefunction becomes real (and
  pre-normalized like LAPACK) density.*/
void internal_add_to_density(struct density* d, struct wavefunction* wf) {
  register int i, j, k, b;
  double re, im;
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  for(b=0; b<wf->bands; b++){
    re = wf->storage[get_wavefunction_storage_index(i, j, k, b, 0, wf->N[0], wf->N[1], wf->bands)];
    im = wf->storage[get_wavefunction_storage_index(i, j, k, b, 0, wf->N[0], wf->N[1], wf->bands)];
    d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] += (re*re + im*im)*wf->charge;
  }
      }
    }
  }
}
void add_to_density(struct density* d, struct wavefunction* wf) {
  internal_add_to_density(d, wf);
}

struct density* get_density_for_list_of_wavefunctions(GList *list, GError** err) {
  if(list==NULL) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 1, "lp_wavefunction.c::get_density_for_list_of_wavefunctions: ERROR: list was empty!\n");
    return NULL;
  }
  struct density *d;
  struct wavefunction *wf = (struct wavefunction*)(list->data);
  if((d = (struct density*)malloc(sizeof(struct density))) == NULL) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 2, "lp_wavefunction.c::get_density_for_list_of_wavefunctions: ERROR: failed to allocate memory for total density structure.");
    return NULL;
  }
  if((d->storage = (double*)malloc(sizeof(double)*(wf->N[0])*(wf->N[1])*(wf->N[2]))) == NULL) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 2, "lp_wavefunction.c::get_density_for_list_of_wavefunctions: ERROR: failed to allocate memory for total density's storage.");
    free(d);
    return NULL;
  }
  for(int i=0; i<3; i++) d->N[i] = wf->N[i];
  for(int i=0; i<3; i++) d->dx[i] = wf->dx[i];
  init_density(d, 0);
  GList *l;
  for(l=list; l != NULL; l=l->next) {
    wf = (struct wavefunction*)l->data;
    for(int i=0; i<3; i++) {
      if((wf->N[i] != d->N[i]) || (wf->dx[i] != d->dx[i])) {
  *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 2, "lp_wavefunction.c::get_density_for_list_of_wavefunctions: ERROR: wavefunction (%s) has an inconsistent size N=(%d, %d, %d) expected (%d, %d, %d); dx=(%g, %g, %g) expected (%g %g %g)", g_file_get_uri(wf->file), wf->N[0], wf->N[1], wf->N[2], d->N[0], d->N[1], d->N[2], wf->dx[0], wf->dx[1], wf->dx[2], d->dx[0], d->dx[1], d->dx[2]);
  free(d->storage);
  free(d);
  return NULL;
      }
      internal_add_to_density(d, wf);
    }
  }
  return d;
}

void load_wavefunctions(GList* wf_list, GError** err) {
  GList *l;
  struct wavefunction* wf;
  int ie;
  GError *e = NULL;
  for(l=wf_list; l!=NULL; l=l->next) {
    wf = (struct wavefunction *)(l->data);
    e=NULL;
    ie = read_wavefunction_gfile(wf->file, &(wf->header), &(wf->storage), &e);
    if((ie != 0) || (e != NULL)) {
      if(e == NULL) {
  *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, ie, "wavefunction::load_wavefunctions: Got a non-zero return value from function_parser::read_wavefunction_gfile: %d\n", ie);
      }else{
  *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, ie, "wavefunction::load_wavefunctions: Got following GError from function_parser::read_wavefunction_gfile: %d/%s\n", e->code, e->message);
      }
      return;
    }
    wf->N[0] = wf->header.Nx;
    wf->N[1] = wf->header.Ny;
    wf->N[2] = wf->header.Nz;
    wf->bands = wf->header.Nc;
    wf->dx[0] = wf->header.dx;
    wf->dx[1] = wf->header.dy;
    wf->dx[2] = wf->header.dz;
  }
}

void put_density_into_metadensity(struct density* md, const struct density* d, const int *Nn, int x, int y, int z, gboolean invert_x, gboolean invert_y, gboolean invert_z, const unsigned* offsets){
  if(d == NULL) return;
  if(d->storage == NULL) return;
  int xoffset = x*(d->N[0]+offsets[0]);
  int yoffset = y*(d->N[1]+offsets[1]);
  int zoffset = z*(d->N[2]+offsets[2]);
  int source_x, source_y, source_z;
  int k, j, i;
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i, source_z, source_y, source_x)
  #endif
  for(k=0; k<d->N[2]; k++) {
    /*The -1 is to ensure that we start at the last *valid* location
     *in the source (and end at the beginning of the source)
     */
    source_z = invert_z?(d->N[2]-k-1):k;
    for(j=0; j<d->N[1]; j++) {
      source_y = invert_y?(d->N[1]-j-1):j;
      for(i=0; i<d->N[0]; i++) {
  source_x = invert_x?(d->N[0]-i-1):i;
  md->storage[get_density_storage_index(i+xoffset, j+yoffset, k+zoffset, md->N[0], md->N[1])]
    = d->storage[get_density_storage_index(source_x, source_y, source_z, d->N[0], d->N[1])];
      }
    }
  }
}

/**Gets a block of the meta_density back out.
 *\param md the meta denstiy
 *\param d an initialized density the sizeo f one block
 *\param buffzone the buffer zone size around the meta_density
 *\param Mn the number of neighbors in each direction
 *\param x the x-coordinate of the block
 *\param y the y-coordinate of the block
 *\param z the z-coordinate of the block
 *\param offsets offset between neighbors in each dimension
 */
void get_density_from_metadensity(const struct density* md, struct density* d, const int* buffzone, const int *Nn, int x, int y, int z, const unsigned* offsets){
  if(d == NULL) return;
  if(d->storage == NULL) return;
  int xoffset = x*(d->N[0]+offsets[0]);
  int yoffset = y*(d->N[1]+offsets[1]);
  int zoffset = z*(d->N[2]+offsets[2]);
  if(buffzone != NULL) {
    xoffset += buffzone[0];
    yoffset += buffzone[1];
    zoffset += buffzone[2];
  }
  int k, j, i;
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i)
  #endif
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]
    = md->storage[get_density_storage_index(i+xoffset, j+yoffset, k+zoffset, md->N[0], md->N[1])];
      }
    }
  }
}

void density_mul_by_const(struct density* d, double c) {
  int i, j, k;
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k)
  #endif
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] *= c;
      }
    }
  }
}

/**
 *The wavefunction itself is normalized such that adding all the values
 *togetheer yields 1. We need the density at each point; hence must divide
 * by dx*dy*dz
 */
void density_un_normalize(struct density* d) {
  double volume = 1.0/((d->dx[0])*(d->dx[1])*(d->dx[2]));
  density_mul_by_const(d, volume);
}


void write_density_slice(enum coord c, const int* index, const struct density* d, char** filename, GError **err, GFile* location, char sepchar, const char* coord_units, const char* density_units) {
  /*If we're iterating on an index, go to the next loop unless the index is right*/
  GError *e = NULL;
  GFile* f;
  GString *s = g_string_new("slice_");
  if(index != NULL) {
    switch(c) {
    case X:
      g_string_append_printf(s, "X=%d.pdata", *index);
      break;
    case Y:
      g_string_append_printf(s, "Y=%d.pdata", *index);
      break;
    case Z:
      g_string_append_printf(s, "Z=%d.pdata", *index);
      break;
    }
    if(location == NULL) {
      f = g_file_new_for_path(s->str);
    }else{
      f = g_file_resolve_relative_path(location, s->str);
    }
  }else{
    f = location;
  }
  *filename = g_file_get_uri(f);
  //GFileIOStream *gfos = g_file_create_readwrite(f, G_FILE_CREATE_REPLACE_DESTINATION, NULL, &e);
  GFileOutputStream *fos = g_file_replace(f, NULL, FALSE, 0, NULL, &e);
  if(e != NULL) {
    *err = e;
    return;
  }
  //GOutputStream *os =
  //g_io_stream_get_output_stream((GIOStream*)gfos);
  GDataOutputStream *dos = g_data_output_stream_new((GOutputStream*)fos);
  if(index != NULL) {
    switch(c) {
    case X:
      g_string_printf(s, "#Y(%s)%cZ(%s)%c(%s)\n", coord_units, sepchar, coord_units, sepchar, density_units);
      break;
    case Y:
      g_string_printf(s, "#X(%s)%cZ(%s)%c(%s)\n", coord_units, sepchar, coord_units, sepchar, density_units);
      break;
    case Z:
      g_string_printf(s, "#X(%s)%cY(%s)%c(%s)\n", coord_units, sepchar, coord_units, sepchar, density_units);
      break;
    }
  }else{
    g_string_printf(s, "#X(%s)%cY(%s)%cZ(%s)%c(%s)\n", coord_units, sepchar, coord_units, sepchar, coord_units, sepchar, density_units);
  }
  e=NULL;
  char *str;
  g_data_output_stream_put_string(dos, s->str, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    //g_io_stream_close((GIOStream*)gfos, NULL, &e);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_slice: had an error (%s) while writing string (%s), but had another error (%s) while closing the file (%s)!\n", (*err)->message, s->str, e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
  for(int k=0; k<d->N[2]; k++) {
    if((index != NULL) && (c == Z) && (k != *index)) continue;
    for(int j=0; j<d->N[1]; j++) {
      if((index != NULL) && (c == Y) && (j != *index)) continue;
      for(int i=0; i<d->N[0]; i++) {
  if((index != NULL) && (c == X) && (i != *index)) continue;
  e=NULL;
  if(index != NULL) {
    switch(c) {
    case X:
      g_string_printf(s, "%g%c%g%c%g\n", j*(d->dx[1]), sepchar, k*(d->dx[2]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
      break;
    case Y:
      g_string_printf(s, "%g%c%g%c%g\n", i*(d->dx[0]), sepchar, k*(d->dx[2]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
      break;
    case Z:
      g_string_printf(s, "%g%c%g%c%g\n", i*(d->dx[0]), sepchar, j*(d->dx[1]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
      break;
    }
  }else{
    g_string_printf(s, "%g%c%g%c%g%c%g\n", i*(d->dx[0]), sepchar, j*(d->dx[1]), sepchar, k*(d->dx[2]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
  }
  g_data_output_stream_put_string(dos, s->str, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    //g_io_stream_close(os, NULL, &e);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_slice: had an error (%s) while writing string (%s), but had another error (%s) while closing the file (%s)!\n", (*err)->message, s->str, e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
      }
    }
  }
  //g_io_stream_close((GIOStream*)gfos, NULL, &e);
  g_output_stream_close((GOutputStream*)fos, NULL, &e);
  if(e != NULL) {
    *err = e;
  }
  g_string_free(s, TRUE);
}

void write_density_sweep_slice(enum coord c, const int* index, const struct density* d, char** filename, GError **err, GFile *location, double maxval, int pixel_strategy) {
  /*If we're iterating on an index, go to the next loop unless the index is right*/
  GError *e;
  GString *s = g_string_new("");
  switch(c) {
  case X:
    g_string_append_printf(s, "X=%d.pbm", *index);
    break;
  case Y:
    g_string_append_printf(s, "Y=%d.pbm", *index);
    break;
  case Z:
    g_string_append_printf(s, "Z=%d.pbm", *index);
    break;
  }
  GFile* f;
  if(location == NULL) {
    f = g_file_new_for_path(s->str);
  }else{
    f = g_file_resolve_relative_path(location, s->str);
  }
  *filename = g_file_get_uri(f);
  e = NULL;
  //GFileIOStream *gfos = g_file_create_readwrite(f, G_FILE_CREATE_REPLACE_DESTINATION, NULL, &e);
  GFileOutputStream *fos = g_file_replace(f, NULL, FALSE, 0, NULL, &e);
  if(e != NULL) {
    *err = e;
    return;
  }
  //GOutputStream *os = g_io_stream_get_output_stream((GIOStream*)gfos);
  char *str;
  /*Calculate the width and height of the resulting image*/
  s = g_string_new("P6\n");
  /*Height and width of the image. Format is width, whitespace, height*/
  int invert_coord = Z;
  switch(c) {
  case X:
    /*In x sweep, width is y height is z*/
    g_string_append_printf(s, "%d\t%d\n", d->N[1], d->N[2]);
    invert_coord=Z;
    break;
  case Y:
    /*In y sweep, width is x height is z*/
    g_string_append_printf(s, "%d\t%d\n", d->N[0], d->N[2]);
    invert_coord=Z;
    break;
  case Z:
    /*In z sweep, width is x height is y*/
    g_string_append_printf(s, "%d\t%d\n", d->N[0], d->N[1]);
    invert_coord=Y;
    break;
  }
  /*Notes from the above:
  ** x is always a width (Y, Z sweeps)
  ** z is always a height (X, Y sweeps)
  ** y is a width if Z sweep, height if X
  ** The below loop should be right: rows are x, then y if not x
  */
  g_string_append_printf(s, "%d\n", MAX_PIXEL_VAL);
#ifdef DEBUG3
  fprintf(stderr, "wavefunc::write_density_sweep_slice: write header %s\n", s->str);
#endif
  e=NULL;
  g_output_stream_write((GOutputStream*)fos, s->str, s->len, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    g_io_stream_close((GIOStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_slice: had an error (%s) while writing string (%s), but had another error (%s) while closing the file (%s)!\n", (*err)->message, s->str, e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
  double val;
  double v;
  double oo_maxval = 1.0/maxval;
  char set[3];
  /**\note that the indices in the dotcode go from bottom to top, while the columns in the ppm file go from top to bottom; we need to go backwards here in order for the image to correspond to what the user is expecting to see*/
  for(int k=0; k<(d->N[2]); k++) {
    if((c == Z) && (k != *index)) continue;
    for(int j=0; j<(d->N[1]); j++) {
      if((c == Y) && (j != *index)) continue;
      for(int i=0; i<(d->N[0]); i++) {
  if((c == X) && (i != *index)) continue;
  e=NULL;
  v = (d->storage[get_density_storage_index((invert_coord==X)?(d->N[0]-i-1):i, (invert_coord==Y)?(d->N[1]-j-1):j, (invert_coord==Z)?(d->N[2]-k-1):k, d->N[0], d->N[1])]);
  /*oo_maxval is already absolute value*/
  if(pixel_strategy == 1) {
    val = log(oo_maxval*fabs(v)+1)/log(2);
  }else if(pixel_strategy == 2) {
    val = (1.0-exp(oo_maxval*fabs(v)))/(1.0-M_E);
  }else{
    val = fabs(oo_maxval*v);
  }
  /*use the following color scheme: blue is positive, red is neg.*/
  if(v>0) {
    set[0]=set[1]=0;
    set[2] = (char)(MAX_PIXEL_VAL*val);
  }else{
    set[2]=set[1]=0;
    set[0] = (char)(MAX_PIXEL_VAL*val);
  }
#ifdef DEBUG3
  fprintf(stderr, "wavefunc::write_density_sweep_slice: i=%d j=%d k=%d val=%g->%d \n", i, j, k, val, set[0]);
#endif
  g_output_stream_write((GOutputStream*)fos, set, 3, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    //g_io_stream_close((GIOStream*)gfos, NULL, &e);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_slice: had an error (%s) while writing (%d,%d,%d->%g->%d), but had another error (%s) while closing the file (%s)!\n", (*err)->message, i, j, k, val, set[0], e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
      }
    }
  }
  //g_io_stream_close((GIOStream*)gfos, NULL, &e);
  g_output_stream_close((GOutputStream*)fos, NULL, &e);
  if(e != NULL) {
    *err = e;
  }
  g_string_free(s, TRUE);
}

/*Gets the *magnitude* of the largest value.*/
double get_max_density_value(const struct density*d) {
  /*Find the max value.*/
  double maxval = -1;
  double v;
  for(int k=0; k<d->N[2]; k++) {
    for(int j=0; j<d->N[1]; j++) {
      for(int i=0; i<d->N[0]; i++) {
  v = d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])];
  if(fabs(v) > maxval) {
    maxval = fabs(v);
  }
      }
    }
  }
  return maxval;
}

void write_density_sweep(enum coord c, const struct density* d, char** filename, GError **err, GFile* location, int pixel_strategy) {
  /*If we're iterating on an index, go to the next loop unless the index is right*/
  GError *e;
  GString *s = g_string_new("sweep_");
  switch(c) {
  case X:
    g_string_append_printf(s, "x");
    break;
  case Y:
    g_string_append_printf(s, "y");
    break;
  case Z:
    g_string_append_printf(s, "z");
    break;
  }
  GFile* f;
  if(location == NULL) {
    f = g_file_new_for_path(s->str);
  }else{
    f = g_file_resolve_relative_path(location, s->str);
  }
  e=NULL;
  g_file_make_directory_with_parents(f, NULL, &e);
  if((e != NULL) && (e->code != G_IO_ERROR_EXISTS)) {
    char* fn = g_file_get_uri(f);
    fprintf(stderr, "E1\n");
    if(e == NULL) {
      *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 32, "wavefunction::write_density_sweep: unable to make directory for writing density sweep (%s)", fn);
    }else{
      g_propagate_prefixed_error(err, e, "wavefunction::write_density_sweep: unable to make directory for writing density sweep (%s): %s", fn, e->message);
    }
    free(fn);
    g_string_free(s, TRUE);
    return;
  }
  double maxval = get_max_density_value(d);
  /*Iterate over the sweep var*/
  char *fn;
  for(int i=0; i<(d->N[c]); i++) {
    e = NULL;
    write_density_sweep_slice(c, &i, d, &fn, &e, f, maxval, pixel_strategy);
    if(e != NULL) {
      *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 32, "Error writing slice into %s: %s", fn, e->message);
      g_error_free(e);
      free(fn);
      return;
    }
    free(fn);
  }
  *filename = s->str;
  /*We don't want to return false, since we're not getting the filename.*/
  g_string_free(s, FALSE);
}


struct density* new_density(const int* N, const double* dx, double initial_value) {
  struct density* d = (struct density*)malloc(sizeof(struct density));
  if(d == NULL) {
    return NULL;
  }
  for(int i=0; i<3; i++) {
    d->N[i] = N[i];
    d->dx[i] = dx[i];
  }
  d->storage = (double*)malloc(N[0]*N[1]*N[2]*sizeof(double));
  if(d->storage == NULL) {
    free(d);
    return NULL;
  }
  int k, j, i;
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i)
  #endif
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] 
    = initial_value;
      }
    }
  }
  return d;
}
struct density* new_density_copy(const struct density* orig) {
  struct density* d = (struct density*)malloc(sizeof(struct density));
  if(d == NULL) {
    return NULL;
  }
  for(int i=0; i<3; i++) {
    d->N[i] = orig->N[i];
    d->dx[i] = orig->dx[i];
  }
  d->storage = (double*)malloc((d->N[0])*(d->N[1])*(d->N[2])*sizeof(double));
  if(d->storage == NULL) {
    free(d);
    return NULL;
  }
  int k, j, i;
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i)
  #endif
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] 
    = orig->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])];
      }
    }
  }
  return d;
}

int op_density_to_density(struct density* to, enum density_op op, const struct density* from) {
  for(int i=0; i<3; i++) {
    if((to->N[i] != from->N[i]) || (to->dx[i] != from->dx[i])) {
      return 1;
    }
  }
  /*The switch can probably be lifted out of the loops by the compiler.*/
  int k, j, i;
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i)
  #endif
  for(k=0; k<to->N[2]; k++) {
    for(j=0; j<to->N[1]; j++) {
      for(i=0; i<to->N[0]; i++) {
  switch(op) {
  case ASSIGN:
    *get_density_storage_at(i, j, k, to) = *get_density_storage_at_c(i, j, k, from);
  case ADD_EQUALS:
    to->storage[get_density_storage_index(i, j, k, to->N[0], to->N[1])] 
      += from->storage[get_density_storage_index(i, j, k, from->N[0], from->N[1])]
      ;
    break;
  case MINUS_EQUALS:
    to->storage[get_density_storage_index(i, j, k, to->N[0], to->N[1])] 
      -= from->storage[get_density_storage_index(i, j, k, from->N[0], from->N[1])]
      ;
    break;
  case MUL_EQUALS:
    to->storage[get_density_storage_index(i, j, k, to->N[0], to->N[1])] 
      *= from->storage[get_density_storage_index(i, j, k, from->N[0], from->N[1])]
      ;
    break;
  case DIV_EQUALS:
    to->storage[get_density_storage_index(i, j, k, to->N[0], to->N[1])] 
      /= from->storage[get_density_storage_index(i, j, k, from->N[0], from->N[1])]
      ;
    break;
  }
      }
    }
  }
  return 0;
}

void free_density(struct density* d) {
  if(d == NULL) return;
  if(d->storage != NULL) free(d->storage);
  free(d);
}

struct dotcode_grid* dotcode_grid_from_density(const struct density *d, GError **err) {
  struct dotcode_grid* g = (struct dotcode_grid*)malloc(sizeof(struct dotcode_grid));
  if(unlikely(g == NULL)) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 100, "wavefunction::dotcode_grid_from_density: Error allocating memory for new dotcode_grid.");
    return NULL;
  }
  g->data = (struct dotcode_tensor*)malloc(sizeof(struct dotcode_tensor));
  if(unlikely(g->data == NULL)) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 101, "wavefunction::dotcode_grid_from_density: Error allocating memory for tensor for the new grid.");
    free(g);
    return NULL;
  }
  g->data->Nd = 3;
  g->data->metadata = g_strdup("2.0");
  g->data->indices = (int*)malloc(2*(g->data->Nd)*sizeof(int));
  if(unlikely(g->data->indices == NULL)) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 102, "wavefunction::dotcode_grid_from_density: Error allocating memory for indicies in the tensor for the new grid.");
    free(g->data);
    free(g);
    return NULL;
  }
  for(int i=0; i<3; i++) {
    (g->gridsites)[i] = (struct dotcode_array*)malloc(sizeof(struct dotcode_array));
    if(unlikely((g->gridsites)[i] == NULL)) {
      *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 103, "wavefunction::dotcode_grid_from_density: Error allocating memory for the %d-th array for the new grid.", i);
      GError *e = NULL;
      for(int j=i-1; j>=0; j++) {
  dotcode_free_type(DOTCODE_ARRAY, (g->gridsites)[j], 1, &e, TRUE);
  if(e != NULL) g_error_free(e);
      }
      free(g->data->indices);
      free(g->data);
      free(g);
      return NULL;
    }
    (g->gridsites)[i]->t = DOTCODE_FLOAT8;
    (g->gridsites)[i]->imin = g->data->indices[2*i] = 0;
    (g->gridsites)[i]->imax = g->data->indices[2*i + 1] = d->N[i]-1;
    (g->gridsites)[i]->storage = (double*)malloc((d->N[i])*sizeof(double));
    if(unlikely((g->gridsites)[i]->storage == NULL)) {
      *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 104, "wavefunction::dotcode_grid_from_density: Error allocating memory for storage in the %d-th array for the new grid.", i);
      GError *e = NULL;
      for(int j=i-1; j>=0; j++) {
  dotcode_free_type(DOTCODE_ARRAY, (g->gridsites)[j], 1, &e, TRUE);
  if(e != NULL) g_error_free(e);
      }
      free(g->data->indices);
      free((g->gridsites)[i]);
      free(g->data);
      free(g);
      return NULL;
    }
    for(int j=0; j<d->N[i]; j++) ((double*)((g->gridsites)[i]->storage))[j] = j*(d->dx[i]);
    /*Put 0 in the center*/
    //g->Ng[i] = d->N[i]-1;
    g->Ng[i] = 0;
  }
  g->data->t = DOTCODE_FLOAT8;
  g->data->metadata = g_strdup("2.0");
  g->data->data = (double*)malloc((d->N[0])*(d->N[1])*(d->N[2])*sizeof(double));
  if(unlikely(g->data->data == NULL)) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 105, "wavefunction::dotcode_grid_from_density: Error allocating memory for storage in the tensor for the new grid.");
    GError *e = NULL;
    for(int j=2; j>=0; j++) {
      dotcode_free_type(DOTCODE_ARRAY, (g->gridsites)[j], 1, &e, TRUE);
      if(e != NULL) g_error_free(e);
    }
    free(g->data->indices);
    free(g->data);
    free(g);
    return NULL;
  }
  /*i is z
   *j is y
   *k is x
   */
  int i, j, k;
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k)
  #endif
  for(i=0; i<d->N[2]; i++) {
    for(j=0; j<d->N[1]; j++) {
      for(k=0; k<d->N[0]; k++) {
  ((double*)(g->data->data))[k + (d->N[0])*(j + (d->N[1])*i)] = ev_to_hartree((d->storage)[get_density_storage_index(k, j, i, d->N[0], d->N[1])]);
      }
    }
  }
  
  return g;
}

struct density* expand_density_box_with_value(const struct density* od, const int* I, double value, GError **err) {
  struct density* d = (struct density*)malloc(sizeof(struct density));
  if(d == NULL) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 106, "libpostproc::wavefunction::expand_density_box_with_material: failed to allocate the new density structure.");
    return NULL;
  }
  for(int i=0; i<3; i++) {
    d->dx[i] = od->dx[i];
    d->N[i] = od->N[i]+2*I[i];
    fprintf(stderr, "dx[%d](%g)->%g d->N[%d](%d)->%d (%d+2*%d)\n", i, od->dx[i], d->dx[i], i, od->N[i], d->N[i], od->N[i], I[i]);
  }
  d->storage = (double*)malloc(sizeof(double)*(d->N[0])*(d->N[1])*(d->N[2]));
  if(d->storage == NULL) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 106, "libpostproc::wavefunction::expand_density_box_with_material: failed to allocate value storage for the new density (size is %dx%dx%d=%ld", d->N[0], d->N[1], d->N[2], ((long int)d->N[0])*((long int)d->N[1])*((long int)d->N[2]));
    free(d);
    return NULL;
  }
  /*Explicitly breaking out the regions saves 1 *branch* per loop*/
  /*First, the top of the box*/
  int i, j, k;
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k)
  #endif
  for(k=0; k<I[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  (d->storage)[get_density_storage_index(i, j, k, d->N[0], d->N[1])] = value;
      }
    }
  }
  /*Now the z-interior of the box*/
  int yoffset, zoffset;
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k, yoffset, zoffset)
  #endif
  for(k=0; k<od->N[2]; k++) {
    zoffset = I[2];

    /*The lower-y outer region*/
    for(j=0; j<I[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  (d->storage)[get_density_storage_index(i, j, k+zoffset, d->N[0], d->N[1])] = value;
      }
    }

    /*Interior of the box in y (and z)*/
    for(j=0; j<od->N[1]; j++) {
      yoffset = I[1];

      for(i=0; i<I[0]; i++) {
  (d->storage)[get_density_storage_index(i, j+yoffset, k+zoffset, d->N[0], d->N[1])] = value;
      }

      /*We're in the middle region. Copy over from the source.*/
      for(i=0; i<od->N[0]; i++) {
  (d->storage)[get_density_storage_index(i+I[0], j+yoffset, k+zoffset, d->N[0], d->N[1])] = (od->storage)[get_density_storage_index(i, j, k, od->N[0], od->N[1])];
      }

      /*we're out the back region.*/
      for(i=0; i<I[0]; i++) {
  (d->storage)[get_density_storage_index(i+(od->N[0])+I[0], j+yoffset, k+zoffset, d->N[0], d->N[1])] = value;
      }

    }

    /*The upper-y outer region*/
    for(j=0; j<I[1]; j++) {
      yoffset = I[1] + od->N[1];
      for(i=0; i<d->N[0]; i++) {
  (d->storage)[get_density_storage_index(i, j+yoffset, k+zoffset, d->N[0], d->N[1])] = value;
      }
    }

  }
  /*Finally, the bottom of the box*/
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k, zoffset)
  #endif
  for(k=0; k<I[2]; k++) {
    zoffset = I[2] + od->N[2];
    for(j=0; j<d->N[2]; j++) {
      for(i=0; i<d->N[0]; i++) {
  (d->storage)[get_density_storage_index(i, j, k+zoffset, d->N[0], d->N[1])] = value;
      }
    }
  }
  return d;
}

struct density* expand_density_box(const struct density* od, const int* I, GError **err) {
  double val = od->storage[0];
  for(int i=0; i<od->N[2]; i++) {
    if((i != 0) && (i != (od->N[2]-1))) continue;
    for(int j=0; j<od->N[1]; j++) {
      if((j != 0) && (j != (od->N[1]-1))) continue;
      for(int k=0; k<od->N[0]; k++) {
  if((k != 0) && (k != (od->N[0]-1))) continue;
  if((od->storage)[get_density_storage_index(k, j, i, od->N[0], od->N[1])] != val) {
    *err = g_error_new(WAVEFUNC_GERROR_DOMAIN, 106, "libpostproc::wavefunction::expand_density_box: value is not identical across the entire outside of the box! (was, %g is %g at %d/%d/%d)", val, (od->storage)[get_density_storage_index(k, j, i, od->N[0], od->N[1])], k, j, i);
    return NULL;
  }
      }
    }
  }
  return expand_density_box_with_value(od, I, val, err);
}

void write_density_povray(const struct density* d, GFile *towrite, int pixel_strategy, GError **err) {
  GError *e = NULL;
  GFileOutputStream *fos = g_file_replace(towrite, NULL, FALSE, 0, NULL, &e);
  if(e != NULL) {
    g_propagate_prefixed_error(err, e, "libpostproc::lp_wavefunction.c::write_density_povray: Error opening output file.");
    return;
  }

  double oo_maxval = 1.0/get_max_density_value(d);

  /*First, write out the initial chunks.*/
  GString* s = g_string_new("//Density vidualziation generated by libpostproc\n");
  s = g_string_append(s, "//\n//Important switches for the end-user are located here.\n");
  s = g_string_append(s, "//Below the ====BLOCKS==== label are the individual data blocks.\n\n");
  s = g_string_append(s, "#include \"colors.inc\"\n");
  s = g_string_append(s, "global_settings { ambient_light rgb<1, 1, 1> }\n");
  s = g_string_append(s, "//background {rgb<1,1,1>}\n\n");
  s = g_string_append(s, "#declare pos_block_color = rgb<0, 0, 1>;\n\n");
  s = g_string_append(s, "#declare neg_block_color = rgb<1, 0, 0>;\n\n");
  g_string_append_printf(s, "camera {\nlocation <%g, %g, %g>\nlook_at <%g, %g, %g>\n}\n\n",
   d->N[0]*d->dx[0]*1.5, d->N[1]*d->dx[1]*1.5, -d->N[2]*d->dx[2]*1.5,
   d->N[0]*d->dx[0]*0.5, d->N[1]*d->dx[1]*0.5, -d->N[2]*d->dx[2]*0.5);
  g_string_append_printf(s, "light_source {<%g, %g, %g> color rgb<1, 1, 1> shadowless}\n",
       d->N[0]*d->dx[0]*1.5, d->N[1]*d->dx[1]*1.5, -d->N[2]*d->dx[2]*1.5);
  s = g_string_append(s, "//====BLOCKS=====\n");
  s = g_string_append(s, "merge {\n");
  g_output_stream_write((GOutputStream*)fos, s->str, s->len, NULL, &e);
  if(e != NULL) {
    g_propagate_prefixed_error(err, e, "libpostproc::wavefunction::write_density_povray: Error while writing header to output stream: %s", e->message);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    return;
  }

  double v, val;
  for(int k=0; k<d->N[2]; k++) {
    for(int j=0; j<d->N[1]; j++) {
      for(int i=0; i<d->N[0]; i++) {
  v = d->storage[get_density_storage_index_(i, j, k, d->N[0], d->N[1])];
  /*oo_maxval is already absolute value*/
  if(pixel_strategy == 1) {
    val = log(oo_maxval*fabs(v)+1)/log(2);
  }else if(pixel_strategy == 2) {
    val = (1.0-exp(oo_maxval*fabs(v)))/(1.0-M_E);
  }else{
    val = fabs(oo_maxval*v);
  }

  if(val < THRESHOLD) continue;
  /*Break the floating-point chunk into pieces. 1 is max, 0 is min, and we have N steps*/
  val *= NUM_STEPS;
  val = floor(val);
  /*Don't waste the processing time if it's totally clear!*/
  if(val == 0) continue;
  val /= NUM_STEPS;


  g_string_printf(s, "box {\n\t<%g, %g, %g>\n\t<%g, %g, %g>\n\ttexture {\n", i*d->dx[0], j*d->dx[1], -k*d->dx[2], (i+1)*d->dx[0], (j+1)*d->dx[1], -(k+1)*d->dx[2]);
  g_string_append_printf(s, "\t\tpigment{\n\t\t\tcolor %s_block_color\n\t\t\ttransmit %g\n\t\t}\n",  (v>=0)?"pos":"neg", 1.0-val);
  s = g_string_append(s, "\t}\n}\n\n");
  g_output_stream_write((GOutputStream*)fos, s->str, s->len, NULL, &e);
  if(e != NULL) {
    g_propagate_prefixed_error(err, e, "libpostproc::wavefunction::write_density_povray: Error while writing block <%d,%d,%d> to output stream: %s", i, j, k, e->message);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    return;
  }
      }
    }
  }

  /*File footer*/
  g_string_printf(s, "}\n");
  g_output_stream_write((GOutputStream*)fos, s->str, s->len, NULL, &e);
  if(e != NULL) {
    g_propagate_prefixed_error(err, e, "libpostproc::wavefunction::write_density_povray: Error while writing file footer: %s\n", e->message);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    return;
  }

  g_output_stream_close((GOutputStream*)fos, NULL, &e);
  if(e != NULL) {
    *err = e;
  }
}


struct density* get_double_resolution_density(const struct density* original, int auto_scale) {
  int N_new[3];
  double dx[3];
  for(int i=0; i<3; i++) {
    N_new[i] = (original->N[i])*2-1;
    dx[i] = original->dx[i]/2.0;
  }
  struct density* nd = new_density(N_new, dx, 0);
  if(nd == NULL) return NULL;
  /*Now do the interpolation.  We loop over the old points, setting the new ones in the inner loop (avoids branch prediction misses if we iterated over the new grid.*/
  long x, y, z;
  double a, b, c, d, e, f, g, h;
  #ifdef OPENMP
  #pragma omp parallel for private(y, x, a, b, c, d, e, f, g, h)
  #endif
  for(z=0; z<original->N[2]-1; z++) {
    for(y=0; y<original->N[1]-1; y++) {
      for(x=0; x<original->N[0]-1; x++) {
  /*For notational simplicity, we'll define the corners of the cube here.
   *\note Starting from the top (x=0 y=0 z=0) and going clockwise (y=1 x=0 z=0), and then the same for the upper 4 points*/
  a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
  b=original->storage[get_density_storage_index(x, y+1, z, original->N[0], original->N[1])];
  c=original->storage[get_density_storage_index(x+1, y+1, z, original->N[0], original->N[1])];
  d=original->storage[get_density_storage_index(x+1, y, z, original->N[0], original->N[1])];
  e=original->storage[get_density_storage_index(x, y, z+1, original->N[0], original->N[1])];
  f=original->storage[get_density_storage_index(x, y+1, z+1, original->N[0], original->N[1])];
  g=original->storage[get_density_storage_index(x+1, y+1, z+1, original->N[0], original->N[1])];
  h=original->storage[get_density_storage_index(x+1, y, z+1, original->N[0], original->N[1])];
  /*Set the original 4 sites (the other sites in the new cell will be set by later iterations*/
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
  nd->storage[get_density_storage_index(2*(x+1), 2*y, 2*z, N_new[0], N_new[1])] = d;
  nd->storage[get_density_storage_index(2*x, 2*(y+1), 2*z, N_new[0], N_new[1])] = b;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*(z+1), N_new[0], N_new[1])] = e;
  /*Set the axial sites*/
  nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z, N_new[0], N_new[1])] = (a+d)/2.0;
  nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b)/2.0;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+e)/2.0;
  /*Set the face sites*/
  nd->storage[get_density_storage_index(2*x+1, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b+c+d)/4.0;
  nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+d+e+h)/4.0;
  nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z+1, N_new[0], N_new[1])] = (a+b+e+f)/4.0;
  /*Set the interior point*/
  nd->storage[get_density_storage_index(2*x+1, 2*y+1, 2*z+1, N_new[0], N_new[1])] = (a+b+c+d+e+f+g+h)/8.0;
  #ifdef DEBUG3
  fprintf(stderr, "density doubling %ld %ld %ld: a=%g b=%g c=%g d=%g e=%g f=%g g=%g h=%g\n", x, y, z, a, b, c, d, e, f, g, h);
  fprintf(stderr, "\toriginal: %ld %ld %ld = %g, %ld %ld %ld = %g, %ld %ld %ld = %g, %ld %ld %ld = %g\n", 
    2*x, 2*y, 2*z, *get_density_storage_at_c_(2*x, 2*y, 2*z, nd),
    2*(x+1), 2*y, 2*z, *get_density_storage_at_c_(2*(x+1), 2*y, 2*z, nd),
    2*x, 2*(y+1), 2*z, *get_density_storage_at_c_(2*x, 2*(y+1), 2*z, nd),
    2*x, 2*y, 2*(z+1), *get_density_storage_at_c_(2*x, 2*y, 2*(z+1), nd)
    );
  fprintf(stderr, "\taxial: %ld %ld %ld = %g, %ld %ld %ld = %g, %ld %ld %ld = %g\n",
    2*x+1, 2*y, 2*z, *get_density_storage_at_c_(2*x+1, 2*y, 2*z, nd),
    2*x, 2*y+1, 2*z, *get_density_storage_at_c_(2*x, 2*y+1, 2*z, nd),
    2*x, 2*y, 2*z+1, *get_density_storage_at_c_(2*x, 2*y, 2*z+1, nd)
    );
  fprintf(stderr, "\tface: %ld %ld %ld = %g, %ld %ld %ld = %g, %ld %ld %ld = %g\n",
    2*x+1, 2*y+1, 2*z, *get_density_storage_at_c_(2*x+1, 2*y+1, 2*z, nd),
    2*x+1, 2*y, 2*z+1, *get_density_storage_at_c_(2*x+1, 2*y, 2*z+1, nd),
    2*x, 2*y+1, 2*z+1, *get_density_storage_at_c_(2*x, 2*y+1, 2*z+1, nd)
    );
  fprintf(stderr, "\tinterior: %ld %ld %ld = %g\n",
    2*x+1, 2*y+1, 2*z+1, *get_density_storage_at_c_(2*x+1, 2*y+1, 2*z+1, nd)
    );
  #endif
      }
    }
  }
  /*All of the points but the final x y and z planes have been set.  There are still axial and facial sites to be set on these.*/
  /*x face*/
  x=original->N[0]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(y, a, b, e, f)
  #endif
  for(z=0; z<original->N[2]-1; ++z) {
    for(y=0; y<original->N[1]-1; ++y) {
  a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
  b=original->storage[get_density_storage_index(x, y+1, z, original->N[0], original->N[1])];
  e=original->storage[get_density_storage_index(x, y, z+1, original->N[0], original->N[1])];
  f=original->storage[get_density_storage_index(x, y+1, z+1, original->N[0], original->N[1])];
  /*Set the original sites (the other sites in the new cell will be set by later iterations*/
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
  nd->storage[get_density_storage_index(2*x, 2*(y+1), 2*z, N_new[0], N_new[1])] = b;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*(z+1), N_new[0], N_new[1])] = e;
  /*Set the axial sites*/
  nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b)/2.0;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+e)/2.0;
  /*Set the face sites*/
  nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z+1, N_new[0], N_new[1])] = (a+b+e+f)/4.0;
    }
  }
  /*y face*/
  y=original->N[1]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(x, a, d, e, h)
  #endif
  for(z=0; z<original->N[2]-1; ++z) {
    for(x=0; x<original->N[0]-1; ++x) {
  a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
  d=original->storage[get_density_storage_index(x+1, y, z, original->N[0], original->N[1])];
  e=original->storage[get_density_storage_index(x, y, z+1, original->N[0], original->N[1])];
  h=original->storage[get_density_storage_index(x+1, y, z+1, original->N[0], original->N[1])];
  /*Set the original sites (the other sites in the new cell will be set by later iterations*/
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
  nd->storage[get_density_storage_index(2*(x+1), 2*y, 2*z, N_new[0], N_new[1])] = d;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*(z+1), N_new[0], N_new[1])] = e;
  /*Set the axial sites*/
  nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z, N_new[0], N_new[1])] = (a+d)/2.0;
  nd->storage[get_density_storage_index(2*x, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+e)/2.0;
  /*Set the face sites*/
  nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+d+e+h)/4.0;
    }
  }
  z=original->N[2]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(x, a, b, c, d)
  #endif
  for(y=0; y<original->N[1]-1; ++y) {
    for(x=0; x<original->N[0]-1; ++x) {
      /*we only get the one set of z sites since the next z site is not extant*/
      a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
      b=original->storage[get_density_storage_index(x, y+1, z, original->N[0], original->N[1])];
      c=original->storage[get_density_storage_index(x+1, y+1, z, original->N[0], original->N[1])];
      d=original->storage[get_density_storage_index(x+1, y, z, original->N[0], original->N[1])];
      /*Set the original sites (the other sites in the new cell will be set by later iterations*/
      nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
      nd->storage[get_density_storage_index(2*(x+1), 2*y, 2*z, N_new[0], N_new[1])] = d;
      nd->storage[get_density_storage_index(2*x, 2*(y+1), 2*z, N_new[0], N_new[1])] = b;
      /*Set the axial sites*/
      nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z, N_new[0], N_new[1])] = (a+d)/2.0;
      nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b)/2.0;
      /*Set the face sites*/
      nd->storage[get_density_storage_index(2*x+1, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b+c+d)/4.0;
    }
  }
  /*Finally, we have some axial sites that still haven't been set along the very edges:*/
  /*x axis*/
  y=original->N[1]-1;
  z=original->N[2]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(a, d)
  #endif
  for(x=0; x<original->N[0]-1; ++x) {
    a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
    d=original->storage[get_density_storage_index(x+1, y, z, original->N[0], original->N[1])];
    /*Set the original sites (the other sites in the new cell will be set by later iterations*/
    nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
    nd->storage[get_density_storage_index(2*(x+1), 2*y, 2*z, N_new[0], N_new[1])] = d;
    /*Set the axial sites*/
    nd->storage[get_density_storage_index(2*x+1, 2*y, 2*z, N_new[0], N_new[1])] = (a+d)/2.0;
  }
  /*y axis*/
  x=original->N[0]-1;
  z=original->N[2]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(a, b)
  #endif
  for(y=0; y<original->N[1]-1; ++y) {
    a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
    b=original->storage[get_density_storage_index(x, y+1, z, original->N[0], original->N[1])];
    /*Set the original sites (the other sites in the new cell will be set by later iterations*/
    nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
    nd->storage[get_density_storage_index(2*x, 2*(y+1), 2*z, N_new[0], N_new[1])] = b;
    /*Set the axial sites*/
    nd->storage[get_density_storage_index(2*x, 2*y+1, 2*z, N_new[0], N_new[1])] = (a+b)/2.0;
  }
  /*z axis*/
  x=original->N[0]-1;
  y=original->N[1]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(a, e)
  #endif
  for(z=0; z<original->N[2]-1; ++z) {
    a=original->storage[get_density_storage_index(x, y, z, original->N[0], original->N[1])];
    e=original->storage[get_density_storage_index(x, y, z+1, original->N[0], original->N[1])];
    /*Set the original sites (the other sites in the new cell will be set by later iterations*/
    nd->storage[get_density_storage_index(2*x, 2*y, 2*z, N_new[0], N_new[1])] = a;
    nd->storage[get_density_storage_index(2*x, 2*y, 2*(z+1), N_new[0], N_new[1])] = e;
    /*Set the axial sites*/
    nd->storage[get_density_storage_index(2*x, 2*y, 2*z+1, N_new[0], N_new[1])] = (a+e)/2.0;
  }
  /*Set the original very corner site*/
  nd->storage[get_density_storage_index(2*(original->N[0]-1), 2*(original->N[1]-1), 2*(original->N[2]-1), N_new[0], N_new[1])] = original->storage[get_density_storage_index(original->N[0]-1, original->N[1]-1, original->N[2]-1, original->N[0], original->N[1])];
  if(auto_scale != 0) {
    double scaling=1;
    int invert=0;
    if(auto_scale < 0) {
      /*Units are grid size, not inverse grid size, so invert*/
      invert=1;
      /*Set it so we get the bits back out*/
      auto_scale=-auto_scale;
    }
    /**since units are f/(dx*dy*dz) -> f/(dx/2*dy/2*dz/2) we need to *divide* f by 2*2*2 to maintain the proper density (f/(dx*dy*dz)->(f/(2*2*2*dx/2*dy/2*dz/2))*/
    if(auto_scale & 1) {
      scaling /= 2;
    }
    if(auto_scale & 2) {
      scaling /= 2;
    }
    if(auto_scale & 4) {
      scaling /= 2;
    }
    if(invert) scaling = 1.0/scaling;
    fprintf(stderr, "DOUBLE-RES SCALING BY %g\n", scaling);
    density_mul_by_const(nd, scaling);
  }
  return nd;
}


/**Gets a new density with a half-resolution grid, using an average for assignment
 *\param original original density
 *\return pointer to a density, or NULL if memory allocation failed.
 *\note the new grid doesn't add an outside row of new sites; it only places a new grid site between two existing ones
 */
struct density* get_half_resolution_density(const struct density* original, GError **e, int auto_scale) {
  int N_new[3];
  double dx[3];
  for(int i=0; i<3; i++) {
    if(original->N[i]%2 != 1) {
      *e = g_error_new(WAVEFUNC_GERROR_DOMAIN, 1, "lp_wavefunction.c::get_half_resolution_density: ERROR: %d-th coordinate is not odd, so there's no clean half-resolution grid!  (This feature is possible and may be added at some point in the future; at the moment, this function is solely for bringining the density back to the original resolution after doubling", i);
      return NULL;
    }
    N_new[i] = (original->N[i]+1)/2;
    dx[i] = original->dx[i]*2.0;
  }
  struct density* nd = new_density(N_new, dx, 0);
  if(nd == NULL) {
    *e = g_error_new(WAVEFUNC_GERROR_DOMAIN, 1, "lp_wavefunction.c::get_half_resolution_density: ERROR: failed to allocate space for half-resolution density.");
    return NULL;
  }
  for(long z=0; z<N_new[2]; ++z) {
    for(long y=0; y<N_new[1]; ++y) {
      for(long x=0; x<N_new[0]; ++x) {
  /*Simple version: just grab the point and neglect the double-res points*/
  nd->storage[get_density_storage_index(x, y, z, N_new[0], N_new[1])] = original->storage[get_density_storage_index(2*x, 2*y, 2*z, original->N[0], original->N[1])];
      }
    }
  }
  if(auto_scale != 0) {
    double scaling=1;
    int invert=1;
    if(auto_scale < 0) {
      /*Units are grid size, not inverse grid size, so don't (this is inverse operation of doubling the resolution)*/
      invert=0;
      /*Set it so we get the bits back out*/
      auto_scale=-auto_scale;
    }
    if(auto_scale & 1) {
      scaling /= 2;
    }
    if(auto_scale & 2) {
      scaling /= 2;
    }
    if(auto_scale & 4) {
      scaling /= 2;
    }
    if(invert) scaling = 1.0/scaling;
    fprintf(stderr, "HALF-RES SCALING BY %g\n", scaling);
    density_mul_by_const(nd, scaling);
  }
  return nd;
}

void density_execute_callback(const struct density* d, void (*callback)(int x, int y, int z, const struct density*d, void* privdat), void* privdat) {
  for(int z=0; z<d->N[2]; ++z) {
    for(int y=0; y<d->N[1]; ++y) {
      for(int x=0; x<d->N[0]; ++x) {
  callback(x, y, z, d, privdat);
      }
    }
  }
}

void density_2_execute_callback(const struct density* d1, const struct density* d2, void (*callback)(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat), void* privdat) {
  for(int z=0; z<d1->N[2]; ++z) {
    for(int y=0; y<d1->N[1]; ++y) {
      for(int x=0; x<d1->N[0]; ++x) {
  callback(x, y, z, d1, d2, privdat);
      }
    }
  }
}

void density_execute_callback_parallel(const struct density* d, void (*callback)(int x, int y, int z, const struct density*d, void* privdat), void** privdat, int* x_report) {
  int x, y, z;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y, z)
  #endif
  for(z=0; z<d->N[2]; ++z) {
    x_report[
             #ifdef OPENMP
       omp_get_thread_num()
             #else
       0
             #endif
       ] = 1;
    for(y=0; y<d->N[1]; ++y) {
      for(x=0; x<d->N[0]; ++x) {
  callback(x, y, z, d, 
     #ifdef OPENMP
     privdat[omp_get_thread_num()]
     #else
     privdat[0]
     #endif
     );
      }
    }
  }
}

void density_2_execute_callback_parallel(const struct density* d1, const struct density* d2, void (*callback)(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat), void** privdat, int* x_report) {
  int x, y, z;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y, z)
  #endif
  for(z=0; z<d1->N[2]; ++z) {
    x_report[
             #ifdef OPENMP
       omp_get_thread_num()
             #else
       0
             #endif
       ] = 1;
    for(y=0; y<d1->N[1]; ++y) {
      for(x=0; x<d1->N[0]; ++x) {
  callback(x, y, z, d1, d2, 
     #ifdef OPENMP
     privdat[omp_get_thread_num()]
     #else
     privdat[0]
     #endif
     );
      }
    }
  }
}


struct diff_data {
  double* minmax;
  /*Value to scale by, ready to *multiply* (i.e. should be inverse of e.g. the global maximum*/
  double scaling;
};
static void get_max_diff_scaled(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  struct diff_data* dd = (struct diff_data*)privdat;
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])*dd->scaling);  
  /**\note that this comparison will fail if diff is a NaN */
  if(!(diff <= (*dd->minmax))) (*dd->minmax) = diff;
}
static void get_min_diff_scaled(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  struct diff_data* dd = (struct diff_data*)privdat;
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])*dd->scaling);  
  if(!(diff >= (*dd->minmax))) (*dd->minmax) = diff;
}

static void get_max_diff_d1(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])]);  
  //fprintf(stderr, "MAX: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  /**\note that this comparison will fail if diff is a NaN */
  if(!(diff <= *(double*)privdat)) *(double*)privdat = diff;
}
static void get_min_diff_d1(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])]);  
  //fprintf(stderr, "MIN: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  if(!(diff >= *(double*)privdat)) *(double*)privdat = diff;
}
static void get_max_diff_d2(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])]);  
  //fprintf(stderr, "MAX: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  /**\note that this comparison will fail if diff is a NaN */
  if(!(diff <= *(double*)privdat)) *(double*)privdat = diff;
}
static void get_min_diff_d2(int x, int y, int z, const struct density* d1, const struct density* d2, void* privdat) {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])]);  
  //fprintf(stderr, "MIN: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  if(!(diff >= *(double*)privdat)) *(double*)privdat = diff;
}

double density_get_max_abs_delta_pct(const struct density* d1, const struct density* d2, int scaleby, double scaling) {
  double max = nan("NaN");
  for(int i=0; i<3; i++) if(d1->N[i] != d2->N[i]) return max;
  #ifdef OPENMP
  int N_threads = omp_get_max_threads();
  fprintf(stderr, "Got %d threads.\n", N_threads);
  struct diff_data dd_[N_threads];
  double max_[N_threads];
  void* maxset[N_threads];
  void* ddset[N_threads];
  int x_report[N_threads];
  int i;
  for(i=0; i<N_threads; ++i) {
    max_[i] = max;
    maxset[i] = (void*)&(max_[i]);
    dd_[i].minmax = &(max_[i]);
    ddset[i] = (void*)&(dd_[i]);
    x_report[i] = 0;
  }
  switch(scaleby) {
  case -1:
    density_2_execute_callback_parallel(d1, d2, get_max_diff_d2, maxset, x_report);
    break;
  case 0:
    density_2_execute_callback_parallel(d1, d2, get_max_diff_d1, maxset, x_report);
    break;
  case 1:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = 1.0/density_get_max(d1, 1);
    }
    density_2_execute_callback_parallel(d1, d2, get_max_diff_scaled, ddset, x_report);
    break;
  case 2:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = 1.0/density_get_max(d2, 1);
    }
    density_2_execute_callback_parallel(d1, d2, get_max_diff_scaled, ddset, x_report);
    break;
  case 3:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = scaling;
    }
    density_2_execute_callback_parallel(d1, d2, get_max_diff_scaled, ddset, x_report);
    break;
  default:
    return max;
  }
  max = max_[0];
  for(int i=1; i<N_threads; ++i) {
    /**\note that max[i] is already one or zero.*/
    fprintf(stderr, "\tmax_[%d]=%g max=%g\n", i, max_[i], max);
    if(x_report[i] && (!(max_[i] <= max))) max = max_[i];
  }
  #else
  struct diff_data dd;
  dd.minmax = &max;
  switch(scaleby) {
  case -1:
    density_2_execute_callback(d1, d2, get_max_diff_d2, &max);
    break;
  case 0:
    density_2_execute_callback(d1, d2, get_max_diff_d1, &max);
    break;
  case 1:
    dd.scaling = 1.0/density_get_max(d1, 1);
    density_2_execute_callback(d1, d2, get_max_diff_scaled, &dd);
    break;
  case 2:
    dd.scaling = 1.0/density_get_max(d2, 1);
    density_2_execute_callback(d1, d2, get_max_diff_scaled, &dd);
    break;
  case 3:
    dd.scaling = scaling;
    density_2_execute_callback(d1, d2, get_max_diff_scaled, &dd);
    break;
  default:
    return max;
  }
  #endif
  return max;
}

double density_get_min_abs_delta_pct(const struct density* d1, const struct density* d2, int scaleby, double scaling) {
  double min = nan("NaN");
  for(int i=0; i<3; i++) if(d1->N[i] != d2->N[i]) return min;
  #ifdef OPENMP
  int N_threads = omp_get_max_threads();
  fprintf(stderr, "Got %d threads.\n", N_threads);
  struct diff_data dd_[N_threads];
  double min_[N_threads];
  void* minset[N_threads];
  void* ddset[N_threads];
  int x_report[N_threads];
  int i;
  for(i=0; i<N_threads; ++i) {
    min_[i] = min;
    minset[i] = &(min_[i]);
    dd_[i].minmax = &(min_[i]);
    ddset[i] = &(dd_[i]);
    x_report[i] = 0;
  }
  switch(scaleby) {
  case -1:
    density_2_execute_callback_parallel(d1, d2, get_min_diff_d2, minset, x_report);
    break;
  case 0:
    density_2_execute_callback_parallel(d1, d2, get_min_diff_d1, minset, x_report);
    break;
  case 1:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = 1.0/density_get_min(d1, 1);
    }
    density_2_execute_callback_parallel(d1, d2, get_min_diff_scaled, ddset, x_report);
    break;
  case 2:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = 1.0/density_get_min(d2, 1);
    }
    density_2_execute_callback_parallel(d1, d2, get_min_diff_scaled, ddset, x_report);
    break;
  case 3:
    for(i=0; i<N_threads; ++i) {
      dd_[i].scaling = scaling;
    }
    density_2_execute_callback_parallel(d1, d2, get_min_diff_scaled, ddset, x_report);
    break;
  default:
    return min;
  }
  min = min_[0];
  for(int i=1; i<N_threads; ++i) {
    /**\note that min[i] is already one or zero.*/
    fprintf(stderr, "\tmin_[%d]=%g min=%g\n", i, min_[i], min);
    if((x_report[i]) && (!(min_[i] >= min))) min = min_[i];
  }
  #else
  struct diff_data dd;
  dd.minmax = &min;
  switch(scaleby) {
  case -1:
    density_2_execute_callback(d1, d2, get_min_diff_d2, &min);
    break;
  case 0:
    density_2_execute_callback(d1, d2, get_min_diff_d1, &min);
    break;
  case 1:
    dd.scaling = 1.0/density_get_min(d1, 1);
    density_2_execute_callback(d1, d2, get_min_diff_scaled, &dd);
    break;
  case 2:
    dd.scaling = 1.0/density_get_min(d2, 1);
    density_2_execute_callback(d1, d2, get_min_diff_scaled, &dd);
    break;
  case 3:
    dd.scaling = scaling;
    density_2_execute_callback(d1, d2, get_min_diff_scaled, &dd);
    break;
  default:
    return min;
  }
  #endif
  return min;
}

struct minmax_dat {
  double minmax;
  int location[3];
};

static void get_max(int x, int y, int z, const struct density* d, void* privdat) {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(mmd->minmax >= *get_density_storage_at_c(x, y, z, d))) {
    mmd->minmax = *get_density_storage_at_c(x, y, z, d);
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}
static void get_min(int x, int y, int z, const struct density* d, void* privdat) {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(mmd->minmax <= *get_density_storage_at_c(x, y, z, d))) {
    mmd->minmax = *get_density_storage_at_c(x, y, z, d);
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}

static void get_max_abs(int x, int y, int z, const struct density* d, void* privdat) {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(fabs(mmd->minmax) >= fabs(*get_density_storage_at_c(x, y, z, d)))) {
    mmd->minmax = fabs(*get_density_storage_at_c(x, y, z, d));
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}
static void get_min_abs(int x, int y, int z, const struct density* d, void* privdat) {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(fabs(mmd->minmax) <= fabs(*get_density_storage_at_c(x, y, z, d)))) {
    mmd->minmax = fabs(*get_density_storage_at_c(x, y, z, d));
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}

double density_get_max(const struct density* d, int abs) {
  int location[3];
  return density_get_max_and_loc(d, abs, location);
}
double density_get_max_and_loc(const struct density* d, int abs, int* location) {
  double max = nan("NaN");
  int N_threads = 0;
  #ifdef OPENMP
  N_threads = omp_get_max_threads();
  #endif
  struct minmax_dat max_[N_threads];
  /*maxset is the array of private data pointers*/
  void* maxset[N_threads];
  int x_report[N_threads];
  for(int i=0; i<N_threads; ++i) {
    max_[i].minmax = nan("NaN");
    maxset[i] = &max_[i];
    x_report[i] = 0;
  }
  if(abs) {
    density_execute_callback_parallel(d, get_max_abs, maxset, x_report);
  }else{
    density_execute_callback_parallel(d, get_max, maxset, x_report);
  }
  max = max_[0].minmax;
  for(int j=0; j<3; ++j) location[j] = max_[0].location[j];
  for(int i=1; i<N_threads; ++i) {
    /**\note that max[i] is already one or zero.*/
    if(x_report[i] && (!(max_[i].minmax <= max))) {
      max = max_[i].minmax;
      for(int j=0; j<3; ++j) location[j] = max_[0].location[j];
    }
  }
  return max;
}
double density_get_min(const struct density* d, int abs) {
  int location[3];
  return density_get_min_and_loc(d, abs, location);
}
double density_get_min_and_loc(const struct density* d, int abs, int* location) {
  double min = nan("NaN");
  int N_threads = 0;
  #ifdef OPENMP
  N_threads = omp_get_max_threads();
  #endif
  struct minmax_dat min_[N_threads];
  /*minset is the array of private data pointers*/
  void* minset[N_threads];
  int x_report[N_threads];
  for(int i=0; i<N_threads; ++i) {
    min_[i].minmax = nan("NaN");
    minset[i] = &min_[i];
    x_report[i] = 0;
  }
  if(abs) {
    density_execute_callback_parallel(d, get_min_abs, minset, x_report);
  }else{
    density_execute_callback_parallel(d, get_min, minset, x_report);
  }
  min = min_[0].minmax;
  for(int j=0; j<3; ++j) location[j] = min_[0].location[j];
  for(int i=1; i<N_threads; ++i) {
    /**\note that min[i] is already one or zero.*/
    if(x_report[i] && (!(min_[i].minmax <= min))) {
      min = min_[i].minmax;
      for(int j=0; j<3; ++j) location[j] = min_[0].location[j];
    }
  }
  return min;
}


 struct gsl_integrand_x {
   gsl_interp_accel* accel;
   gsl_spline* spline;
   double(*f)(double x, double y, double z, void* privdat);
   void* privdat;
   double y;
   double z;
   FILE* outfile;
   int iter;
 };

 hotfunc double gsl_integrand_x_func(double x, void* privdat) {
   struct gsl_integrand_x* gix = (struct gsl_integrand_x*)privdat;
   double val = gix->f(x, gix->y, gix->z, gix->privdat);
   double splineval = gsl_spline_eval(gix->spline, x, gix->accel);
   fprintf(gix->outfile, "%d\t%g\t%g\t%g\t%g\t%g\t%g\n", ++gix->iter, x, gix->y, gix->z, val, splineval, val*splineval);
   return val*splineval;
 }

 hotfunc double gsl_integrand_x_nofunc(double x, void* privdat) {
   struct gsl_integrand_x* gix = (struct gsl_integrand_x*)privdat;
   double splineval = gsl_spline_eval(gix->spline, x, gix->accel);
   fprintf(gix->outfile, "%d\t%g\t%g\t%g\t%g\n", ++gix->iter, x, gix->y, gix->z, splineval);
   return splineval;
 }
 
 hotfunc double gsl_integrand_yz(double coord, void* privdat) {
   struct gsl_integrand_x* gix = (struct gsl_integrand_x*)privdat;
   double splineval = gsl_spline_eval(gix->spline, coord, gix->accel);
   fprintf(gix->outfile, "%d\t%g\t%g\n", ++gix->iter, coord, splineval);
   return splineval;
 }

/**Note to self: this will be much cleaner for a templated version.  What would happen is this:
 1) Allocate array of doubles
 2) Integrate each point in this dimension (e.g. if we're in 3D, do a 2D integral at each z point)
 3) Integrate the array of integrals
 4) We're done!
 
 What's more, one could set an extra argument (of dimension D-1, if template is for dimension D) which specifies the index to integrate at each step, e.g.
 int integrate_first[] = {2, 1};
 would specify that first the Z coordinate is to be integrated over (3D) then in the 2D integral, integrate over Y.  When we're at 1D, we have no further coordinate to specify.

 Of course, 0 has to be a specialized template.

 if this isnt' a class member, remember that we must fully specialize the template to provide the termination condition.

 Of course, this would obscure the fact that the strid is critically important to performance.  
 And we'd have to fetch out full on classes of the D-1 dim.
*/
double density_integrate_spline(double* z_error_estimate, double* y_error_estimate, double* x_error_estimate, double relerr, double abserr, const struct density* d, double (*f)(double x, double y, double z, void* privdat), void* privdat) {
   /**Integrating x results in a yz plane of values
    *\note that the y coords are contiguous in memory
    */
   double* yzvals = (double*)malloc(sizeof(double)*d->N[1]*d->N[2]);
   /*This array has the coordinate we're integrating over*/
   double* coords = (double*)malloc(sizeof(double)*d->N[0]);
   if((yzvals == NULL) || (coords == NULL)) {
     if(yzvals != NULL) free(yzvals);
     if(coords != NULL) free(coords);
     return nan("NaN");
   }
   for(int x=0; x<d->N[0]; ++x) coords[x] = (x - d->N[0]/2)*d->dx[0];
   fprintf(stderr, "coords[0]=%g\tcoords[%d]=%g\n", coords[0], d->N[0]-1, coords[d->N[0]-1]);
   int index, z, y, status;
   gsl_interp_accel* accel = gsl_interp_accel_alloc();
   gsl_spline* spline =  gsl_spline_alloc(gsl_interp_cspline, d->N[0]);
   gsl_integration_workspace *w = gsl_integration_workspace_alloc(MAX_GSL_INTEGRATION_INTERVALS);
   struct gsl_integrand_x gix = {.accel=accel, .spline=spline, .f=f, .privdat=privdat, 0, 0};
   gsl_function gsl_func = {.function = gsl_integrand_x_nofunc, .params=&gix};
   if(f != NULL) gsl_func.function = gsl_integrand_x_func;
   double err_estimate;
   GString* filename = g_string_new("");
   status = 0;
   for(z=0; z<d->N[0]; ++z) {
     gix.z = (z - d->N[2]/2)*d->dx[2];
     for(y=0; y<d->N[1]; ++y) {
       gix.y = (y - d->N[1]/2)*d->dx[1];
       index = z*d->N[1] + y;
       g_string_printf(filename, "x_integration_y_%d_z_%d.pdata", y, z);
       gix.outfile = fopen(filename->str, "w+");
       gix.iter = -1;
       /*\note The x coordinates are contiguous; if we get the address for x=0, just move to the next double to get x=1 and so on.  That's why this next line works*/
       gsl_spline_init(spline, coords, get_density_storage_at_c(0, y, z, d), d->N[0]);
       status |= gsl_integration_qag(&gsl_func, coords[0], coords[d->N[0]-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, &(yzvals[index]), (x_error_estimate==NULL)?(&err_estimate):(&x_error_estimate[index]));
       fclose(gix.outfile);
     }
   }
   free(coords);
   gsl_spline_free(spline);
   if(status != 0) {
     free(yzvals);
     gsl_interp_accel_free(accel);
     gsl_spline_free(spline);
     free(coords);
     return nan("NaN");
   }

   /*Now do the y integration; we end up with a z array of values.*/
   coords = (double*)malloc(sizeof(double)*d->N[1]);
   double* zvals = (double*)malloc(sizeof(double)*d->N[1]);
   for(int y=0; y<d->N[1]; ++y) coords[y] = (y - d->N[1]/2)*d->dx[1];
   fprintf(stderr, "coords[0]=%g\tcoords[%d]=%g\n", coords[0], d->N[1]-1, coords[d->N[1]-1]);
   if((zvals == NULL) || (coords == NULL)){
     free(yzvals);
     gsl_interp_accel_free(accel);
     gsl_spline_free(spline);
     if(coords != NULL) free(coords);
     if(zvals != NULL) free(zvals);
     return nan("NaN");
   }
   spline =  gsl_spline_alloc(gsl_interp_cspline, d->N[1]);
   gsl_func.function = gsl_integrand_yz;
   for(int z=0; z<d->N[0]; ++z) {
     g_string_printf(filename, "y_integration_z_%d.pdata", z);
     gix.outfile = fopen(filename->str, "w+");
     gix.iter = -1;
     gsl_spline_init(spline, coords, &yzvals[z*d->N[1]], d->N[1]);
     status |= gsl_integration_qag(&gsl_func, coords[0], coords[d->N[0]-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, zvals, (y_error_estimate==NULL)?(&err_estimate):(&y_error_estimate[z]));
     fclose(gix.outfile);
   }
   free(coords);
   free(yzvals);
   gsl_spline_free(spline);
   if(status != 0) {
     gsl_interp_accel_free(accel);
     gsl_spline_free(spline);
     free(coords);
     free(zvals);
     return nan("NaN");
   }

   /*Finally, do the z integration.*/
   coords = (double*)malloc(sizeof(double)*d->N[2]);
   for(int z=0; z<d->N[2]; ++z) coords[z] = (z - d->N[2]/2)*d->dx[2];
   fprintf(stderr, "coords[0]=%g\tcoords[%d]=%g\n", coords[0], d->N[2]-1, coords[d->N[2]-1]);
   if(coords == NULL){
     gsl_interp_accel_free(accel);
     gsl_spline_free(spline);
     free(zvals);
     if(coords != NULL) free(coords);
     return nan("NaN");
   }
   spline =  gsl_spline_alloc(gsl_interp_cspline, d->N[2]);
   gsl_spline_init(spline, coords, zvals, d->N[2]);
   double final_value = nan("NaN");
   g_string_printf(filename, "y_integration_z_%d.pdata", z);
   gix.outfile = fopen(filename->str, "w+");
   gix.iter = -1;
   status = gsl_integration_qag(&gsl_func, coords[0], coords[d->N[0]-1], abserr, relerr, MAX_GSL_INTEGRATION_INTERVALS, GSL_INTEG_GAUSS61, w, &final_value, (z_error_estimate==NULL)?(&err_estimate):(z_error_estimate));
   fclose(gix.outfile);

   gsl_interp_accel_free(accel);
   gsl_spline_free(spline);
   free(coords);
   free(zvals);

   if(status != 0) return nan("NaN");
   return final_value;
 }

double density_integrate_sum(const struct density* d, double (*f)(double x, double y, double z, void* privdat), void* privdat) {
  double val;
  int x, y, z;
  double sum = 0;
  for(z=0; z<d->N[2]; ++z) {
    for(y=0; y<d->N[1]; ++y) {
      for(x=0; x<d->N[0]; ++x) {
  val = *get_density_storage_at_c_(x, y, z, d);
  if(f != NULL) {
    val *= f((x-d->N[0])*d->dx[0], (y-d->N[1])*d->dx[0], (z-d->N[2])*d->dx[0], privdat);
  }
  sum += val;
      }
    }
  }
  return sum;
}

double density_integrate_sum_parallel(const struct density* d, double (*f)(double x, double y, double z, void* privdat), void** privdat) {
  double val;
  int x, y, z;
  double sum = 0;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y, z, val) reduction(+:sum)
  #endif
  for(z=0; z<d->N[2]; ++z) {
    for(y=0; y<d->N[1]; ++y) {
      for(x=0; x<d->N[0]; ++x) {
  val = *get_density_storage_at_c_(x, y, z, d);
  if(f != NULL) {
    if(privdat != NULL) {
      val *= f((x-d->N[0])*d->dx[0], (y-d->N[1])*d->dx[0], (z-d->N[2])*d->dx[0], privdat[
                             #ifdef OPENMP
                             omp_get_thread_num()
                             #else
                             0
                             #endif
]);

    }else{
      val *= f((x-d->N[0])*d->dx[0], (y-d->N[1])*d->dx[0], (z-d->N[2])*d->dx[0], NULL);
    }
  }
  sum += val;
      }
    }
  }
  return sum;
}

void write_density_line(enum coord axis, const int* gridsite, const struct density* d, char** filename, GError **err, GFile* location, char sepchar, const char* coord_units, const char* density_units) {
  /*If we're iterating on an index, go to the next loop unless the index is right*/
  GError *e = NULL;
  GFile* f;
  GString *s = g_string_new("line_");
  switch(axis) {
  case X:
    g_string_append_printf(s, "Y=%d_Z=%d.pdata", gridsite[1], gridsite[2]);
    break;
  case Y:
    g_string_append_printf(s, "X=%d_Z=%d.pdata", gridsite[0], gridsite[2]);
    break;
  case Z:
    g_string_append_printf(s, "X=%d_Y=%d.pdata", gridsite[0], gridsite[1]);
    break;
  }
  if(location == NULL) {
    f = g_file_new_for_path(s->str);
  }else{
    f = g_file_resolve_relative_path(location, s->str);
  }
  *filename = g_file_get_uri(f);
  //GFileIOStream *gfos = g_file_create_readwrite(f, G_FILE_CREATE_REPLACE_DESTINATION, NULL, &e);
  GFileOutputStream *fos = g_file_replace(f, NULL, FALSE, 0, NULL, &e);
  if(e != NULL) {
    *err = e;
    return;
  }
  //GOutputStream *os =
  //g_io_stream_get_output_stream((GIOStream*)gfos);
  GDataOutputStream *dos = g_data_output_stream_new((GOutputStream*)fos);
  switch(axis) {
  case X:
    g_string_printf(s, "#X(%s)%c(%s)\n", coord_units, sepchar, density_units);
    break;
  case Y:
    g_string_printf(s, "#Y(%s)%c(%s)\n", coord_units, sepchar, density_units);
    break;
  case Z:
    g_string_printf(s, "#Z(%s)%c(%s)\n", coord_units, sepchar, density_units);
    break;
  }
  e=NULL;
  char *str;
  g_data_output_stream_put_string(dos, s->str, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    //g_io_stream_close((GIOStream*)gfos, NULL, &e);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_line: had an error (%s) while writing string (%s), but had another error (%s) while closing the file (%s)!\n", (*err)->message, s->str, e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
  for(int k=0; k<d->N[2]; k++) {
    if(((axis == X) || (axis == Y)) && (k != gridsite[2])) continue;
    for(int j=0; j<d->N[1]; j++) {
      if(((axis == X) || (axis == Z)) && (j != gridsite[1])) continue;
      for(int i=0; i<d->N[0]; i++) {
  if(((axis == Z) || (axis == Y)) && (i != gridsite[0])) continue;
  e=NULL;
  switch(axis) {
  case X:
    g_string_printf(s, "%g%c%g\n", i*(d->dx[0]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
    break;
  case Y:
    g_string_printf(s, "%g%c%g\n", j*(d->dx[1]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
    break;
  case Z:
    g_string_printf(s, "%g%c%g\n", k*(d->dx[2]), sepchar, d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])]);
    break;
  }
  g_data_output_stream_put_string(dos, s->str, NULL, &e);
  if(e != NULL) {
    *err = e;
    g_string_free(s, TRUE);
    e=NULL;
    //g_io_stream_close(os, NULL, &e);
    g_output_stream_close((GOutputStream*)fos, NULL, &e);
    if(e != NULL) {
      fprintf(stderr, "wavefunction::write_density_line: had an error (%s) while writing string (%s), but had another error (%s) while closing the file (%s)!\n", (*err)->message, s->str, e->message, str=g_file_get_uri(f));
      g_free(str);
    }
    return;
  }
      }
    }
  }
  //g_io_stream_close((GIOStream*)gfos, NULL, &e);
  g_output_stream_close((GOutputStream*)fos, NULL, &e);
  if(e != NULL) {
    *err = e;
  }
  g_string_free(s, TRUE);
}