fftw_wavefunc.c

Go to the documentation of this file.

/*
** fftw_wavefunc.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  Mon Nov 16 06:34:41 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdlib.h>
#include <stdio.h>
#include <fftw3.h>
#include <math.h>
#include <glib-2.0/glib.h>
#ifdef OPENMP
#include <omp.h>
#endif
#include <libpostproc/fftw_wavefunc.h>

#define FFTW_WAVEFUNC_GERROR_DOMAIN 32

#ifdef ACMLMATH
  #include <acml_mv.h>
  #define COS fastcos
#else
  #define COS cos
#endif


void init_karray(fftw_complex* karray, int *ranges, double re, double im) {
  int i, j, k;
  #ifdef OPENMP
  #pragma omp parallel for private(i, j, k)
  #endif
  for(k=0; k<ranges[2]; k++) {
    for(j=0; j<ranges[1]; j++) {
      for(i=0; i<ranges[0]; i++) {
  karray[i + ranges[0]*(j + ranges[1]*k)][0] = re;
  karray[i + ranges[0]*(j + ranges[1]*k)][1] = im;
      }
    }
  }
}

void init_karray_1d(fftw_complex* karray, int range, double re, double im) {
  int i;
  #ifdef OPENMP
  #pragma omp parallel for private(i)
  #endif
  for(i=0; i<range; i++) {
    karray[i][0] = re;
    karray[i][1] = im;
  }
}

/*NOTES:
**For index ordering, examine the figures at
****http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
*/
struct density* get_potential_for_density(const struct density* rho, const struct density* epsilon, struct density** rhoprime_out, const struct density* dielectric_term_k, const struct density* dielectric_term_r, GError **e) {
  for(int i=0; i<3; i++) {
    if(rho->N[i] != epsilon->N[i]) {
      *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density: charge density and dielectric constant grid sizes for dimension %d are not the same! (%d and %d, respectively)\n", i, rho->N[i], epsilon->N[i]);
      return NULL;
    }
  }
  /*If instructed to, intialize the FFTW threading*/
  #ifdef OPENMP
  int errval = fftw_init_threads();
  if(errval == 0) {
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density: Failed to initialize FFTW threading: error %d", errval);
    return NULL;
  }
  int N_threads = omp_get_max_threads();
  fftw_plan_with_nthreads(N_threads);
  #endif
  
  fprintf(stderr, "In get_potential_for_density: \n");
  /*Construct rhoprime = rho/epsilon @ each lattice site*/
  struct density* rhoprime = new_density_copy(rho);
  if(rhoprime == NULL) {
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for rhoprime.");
    return NULL;
  }
  fprintf(stderr, "\t allocated rhoprime\n");
  if(op_density_to_density(rhoprime, DIV_EQUALS, epsilon)) {
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 101, "fftw_wavefunc.c::get_potential_for_density: failed to divide rho by epsilon at each point");
    free_density(rhoprime);
    return NULL;
  }
  fprintf(stderr, "\tdivided by epsilon\n");
  /*Add in the realspace dielectric term if it's been specified.*/
  if(dielectric_term_r != NULL) {
    op_density_to_density(rhoprime, ADD_EQUALS, dielectric_term_r);
  }
  fprintf(stderr, "\tadded dielectric term\n");
  /*Save a copy for output, if requested.*/
  if(rhoprime_out != NULL) *rhoprime_out = rhoprime;

  /*Centralizes the size of k-space arrays*/
  int k_ranges[3] = {rhoprime->N[0]/2+1, rhoprime->N[1], rhoprime->N[2]};
  //if((rhoprime->N[0])%2 == 1) k_ranges[0]++;

  fftw_complex* rhoprime_k = (fftw_complex*)malloc(k_ranges[0]*k_ranges[1]*k_ranges[2]*sizeof(fftw_complex));
  if(rhoprime_k == NULL) {
    if(rhoprime_out == NULL) free_density(rhoprime);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for k-space rhoprime.");
    return NULL;
  }
  //fprintf(stderr, "\tinitializing rhoprime_k array with 0\n");
  //init_karray(rhoprime_k, k_ranges, 0, 0);
  /*Note that indices are in opposite order from how we've been using
    them; is to make sure fftw gets the mapping into memory right.*/
  fftw_plan p_f = fftw_plan_dft_r2c_3d(rhoprime->N[2], rhoprime->N[1], rhoprime->N[0], rhoprime->storage, rhoprime_k, FFTW_ESTIMATE);
  fprintf(stderr, "\tExecuting FFT\n");
  fftw_execute(p_f);
  /*Free the plan now that we're done with it.*/
  fftw_destroy_plan(p_f);

  /*Calculate the potential in k-space*/
  /*The "cosine term" is
   *cos(2*pi*i/(N[0]/2))+cos(2*pi*j/N[1])+cos(2*pi*k/N[2])-4
   */
  fftw_complex* phi_k = (fftw_complex*)malloc(k_ranges[0]*k_ranges[1]*k_ranges[2]*sizeof(fftw_complex));
  if(phi_k == NULL ){
    if(rhoprime_out == NULL) free_density(rhoprime);
    free(rhoprime_k);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for k-space potential.");
    return NULL;
  }
  //fprintf(stderr, "\tinitializing phi_k with 0\n");
  //init_karray(phi_k, k_ranges, 0, 0);
  /*Algo is from Numerical Recipes, 3rd ed., p 1055*/
  /*NOTE THAT DELTA^2 IS NO LONGER IN HERE; DIVIDE THE WAVEFUNCTION BY
    A *SINGLE* GRID SIDE LENGTH BEFORE CALLING HERE!*/
  int k, j, i;
  fprintf(stderr, "\t setting phi_k\n");
  #ifdef OPENMP
  #pragma omp parallel for private(k, j, i)
  #endif
  for(k=0; k<k_ranges[2]; k++) {
    for(j=0; j<k_ranges[1]; j++) {
      for(i=0; i<k_ranges[0]; i++) {
  phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0] 
    /*TODO: fix this to not depend on dx->[0] only*/
    = rhoprime_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0] 
    / (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N[0])))
      + COS(2.0*M_PI*((double)j)/((double)(rhoprime->N[1])))
      + COS(2.0*M_PI*((double)k)/((double)(rhoprime->N[2])))
      - 3.0));
  phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1] 
    /*TODO: fix this to not depend on dx->[0] only*/
    =rhoprime_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1] 
    / (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N[0])))
      + COS(2.0*M_PI*((double)j)/((double)(rhoprime->N[1])))
      + COS(2.0*M_PI*((double)k)/((double)(rhoprime->N[2])))
      - 3.0));
  if(dielectric_term_k != NULL)
    /**\note the dielectric term is real*/
    /**\note the d^2/dx^2 on the left-hand side introduces *dx^2, so we end up with one dx; the two dxes are implicitly cancelled in psi.*/
    phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0] += *get_density_storage_at_c_(i, j, k, dielectric_term_k);
#ifdef DEBUG3
  fprintf(stderr, "i=%d j=%d k=%d -> %g + i%g\n", i, j, k, phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][0], phi_k[i + k_ranges[0]*(j + k_ranges[1]*k)][1]);
#endif
      }
    }
  }
  /*Set the 0th-component to 0; there's an inf and a nan otherwise;
    rho is disconnected from phi_k at this point.*/
  phi_k[0][0] = phi_k[0][1] = 0.0;

  /*reverse plan*/
  /*The final potential (density structure works as well for both)
    will be saved in phi.*/
  struct density* phi = new_density(rho->N, rho->dx, 0.0);
  if(phi == NULL) {
    if(rhoprime_out == NULL) free_density(rhoprime);
    free(rhoprime_k);
    free(phi_k);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density: failed to allocate storage for realspace potential.");
    return NULL;
  }
  fftw_plan p_r = fftw_plan_dft_c2r_3d(phi->N[2], phi->N[1], phi->N[0], phi_k, phi->storage, FFTW_ESTIMATE);
  fprintf(stderr, "\treverse-transforming\n");
  fftw_execute(p_r);
  /*Free the plan now that we're done with it.*/
  fftw_destroy_plan(p_r);
  /*Free our structures*/
  if(rhoprime == NULL) free_density(rhoprime);
  free(rhoprime_k);
  free(phi_k);
  if(phi->storage == NULL) {
    fprintf(stderr, "After reverse FFT, phi was NULL!\n");
  }
  /*After reverse-transforming, we need to scale the data; the
    transforms are not normalized and must be scaled by 1/#real_elements
    see http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
  */
  fprintf(stderr, "\tmultiplying by size of array\n");
  density_mul_by_const(phi, -1.0/((phi->N[0])*(phi->N[1])*(phi->N[2])));
  return phi;
}

struct density_1d* get_potential_for_density_1d(const struct density_1d* rho, const struct density_1d* epsilon, struct density_1d** rhoprime_out, GError **e) {
  if(rho->N != epsilon->N) {
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 100, "fftw_wavefunc.c::get_potential_for_density_1d: charge density_1d and dielectric constant grid sizes for are not the same! (%d and %d, respectively)\n", rho->N, epsilon->N);
    return NULL;
  }

  /*Construct rhoprime = rho/epsilon @ each lattice site*/
  struct density_1d* rhoprime = new_density_1d_copy(rho);
  if(rhoprime == NULL) {
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for rhoprime");
    return NULL;
  }
  if(op_density_1d_to_density_1d(rhoprime, DIV_EQUALS, epsilon)) {
    free_density_1d(rhoprime);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 101, "fftw_wavefunc.c::get_potential_for_density_1d: failed to divide rho by epsilon at each point");
    return NULL;
  }
  /*Save a copy for output, if requested.*/
  if(rhoprime_out != NULL) *rhoprime_out = rhoprime;

  /*Centralizes the size of k-space arrays*/
  int k_range = rhoprime->N/2 + 1;
  if((rhoprime->N)%2 == 1) k_range++;

  fftw_complex* rhoprime_k = (fftw_complex*)malloc(k_range*sizeof(fftw_complex));
  if(rhoprime_k == NULL) {
    if(rhoprime_out == NULL) free_density_1d(rhoprime);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for k-space rhoprime");
    return NULL;
  }
  init_karray_1d(rhoprime_k, k_range, 0, 0);
  /*Note that indices are in opposite order from how we've been using
    them; is to make sure fftw gets the mapping into memory right.*/
  fprintf(stderr, "rhoprime->N=%u krange=%d\n", rhoprime->N, k_range);
  fftw_plan p_f = fftw_plan_dft_r2c_1d(rhoprime->N, rhoprime->storage, rhoprime_k, FFTW_ESTIMATE);
  fftw_execute(p_f);

  /*Calculate the potential in k-space*/
  /*The "cosine term" is
   *cos(2*pi*i/(N[0]/2))+cos(2*pi*j/N[1])+cos(2*pi*k/N[2])-4
   */
  fftw_complex* phi_k = (fftw_complex*)malloc(k_range*sizeof(fftw_complex));
  if(phi_k == NULL ){
    if(rhoprime_out == NULL) free_density_1d(rhoprime);
    free(rhoprime_k);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for k-space potential");
    return NULL;
  }
  init_karray_1d(phi_k, k_range, 0, 0);
  /*Algo is from Numerical Recipes, 3rd ed., p 1055*/
  /*NOTE THAT DELTA^2 IS NO LONGER IN HERE; DIVIDE THE WAVEFUNCTION BY
    A *SINGLE* GRID SIDE LENGTH BEFORE CALLING HERE!*/
  for(int i=0; i<k_range; i++) {
    phi_k[i][0] 
      /*TODO: fix this to not depend on dx->[0] only*/
      = rhoprime_k[i][0] 
      / (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N)))
        - 1.0))*rhoprime->dx*rhoprime->dx;
    phi_k[i][1] 
      /*TODO: fix this to not depend on dx->[0] only*/
      = rhoprime_k[i][1] 
      / (2.0*(COS(2.0*M_PI*((double)i)/((double)(rhoprime->N)))
        - 1.0)) * rhoprime->dx*rhoprime->dx;
#ifdef DEBUG3
    fprintf(stderr, "i=%d -> %g + i%g\n", i, phi_k[i][0], phi_k[i][1]);
#endif
  }
  /*Set the 0th-component to 0; there's an inf and a nan otherwise;
    rho is disconnected from phi_k at this point.*/
  //phi_k[0][0] = phi_k[0][1] = 0.0;

  /*reverse plan*/
  /*The final potential (density_1d structure works as well for both)
    will be saved in phi.*/
  struct density_1d* phi = new_density_1d(rho->N, rho->dx, 0.0);
  if(phi == NULL) {
    if(rhoprime_out == NULL) free_density_1d(rhoprime);
    free(rhoprime_k);
    free(phi_k);
    *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 1, "fftw_wavefunc.c::get_potential_for_density_1d: failed to allocate storage for realspace potential");
    return NULL;
  }
  fftw_plan p_r = fftw_plan_dft_c2r_1d(phi->N, phi_k, phi->storage, FFTW_ESTIMATE);
  fftw_execute(p_r);
  /*Free our structures*/
  if(rhoprime == NULL) free_density_1d(rhoprime);
  free(rhoprime_k);
  free(phi_k);
  if(phi->storage == NULL) {
    fprintf(stderr, "After reverse FFT, phi was NULL!\n");
  }
  /*After reverse-transforming, we need to scale the data; the
    transforms are not normalized and must be scaled by 1/#real_elements
    see http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
  */
  density_1d_mul_by_const(phi, -1.0/(phi->N));
  return phi;
}

/**sets the del_phi term
 *\param del_phi_term density struct to store the del_phi term; must be epsilon*N[i]*2-1 (i.e. doubled-grid size)
 *\param epsilon dielectric distribution; must be same dimensions as phi
 *\param phi electrostatic potential from the previous iteration
 *\return nothing
 */
void set_del_phi_term(struct density *del_phi_term, const struct density *epsilon, const struct density* phi) {
  /*First, set everything but the final x, y, and z faces*/
  long x, y, z, ex, ox, ey, oy, ez, oz;
  double phi_a, phi_b, phi_c, phi_d, phi_e, phi_f, phi_g, phi_h, epsilon_a, epsilon_b, epsilon_c, epsilon_d, epsilon_e, epsilon_f, epsilon_g, epsilon_h, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC;
  #ifdef OPENMP
#pragma omp parallel for private(y, x, ex, ox, ey, oy, ez, oz, phi_a, phi_b, phi_c, phi_d, phi_e, phi_f, phi_g, phi_h, epsilon_a, epsilon_b, epsilon_c, epsilon_d, epsilon_e, epsilon_f, epsilon_g, epsilon_h, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC)
  #endif
  for(z=0; z<(epsilon->N[2]-1); ++z) {
    ez=2*z;
    oz=ez+1;
    /**even-z loop*/
    for(y=0; y<(epsilon->N[1]-1); ++y){
      ey=2*y;
      oy=ey+1;
      /**even-y loop*/
      for(x=0; x<(epsilon->N[0]-1); ++x){
  /*Get phi and epsilon for the surrounding sites

    top level:      y+1   B    C

                            y   A    D   z
                                x   x+1



    bottom level:   y+1   F    G
  
                            y   E    H   z+1
                x   x+1
  */         
  phi_a=phi->storage[get_density_storage_index_(x, y, z, phi->N[0], phi->N[1])];
  phi_b=phi->storage[get_density_storage_index_(x, y+1, z, phi->N[0], phi->N[1])];
  phi_c=phi->storage[get_density_storage_index_(x+1, y+1, z, phi->N[0], phi->N[1])];
  phi_d=phi->storage[get_density_storage_index_(x+1, y, z, phi->N[0], phi->N[1])];
  phi_e=phi->storage[get_density_storage_index_(x, y, z+1, phi->N[0], phi->N[1])];
  phi_f=phi->storage[get_density_storage_index_(x, y+1, z+1, phi->N[0], phi->N[1])];
  phi_g=phi->storage[get_density_storage_index_(x+1, y+1, z+1, phi->N[0], phi->N[1])];
  phi_h=phi->storage[get_density_storage_index_(x+1, y, z+1, phi->N[0], phi->N[1])];
  epsilon_a=epsilon->storage[get_density_storage_index_(x, y, z, epsilon->N[0], epsilon->N[1])];
  epsilon_b=epsilon->storage[get_density_storage_index_(x, y+1, z, epsilon->N[0], epsilon->N[1])];
  epsilon_c=epsilon->storage[get_density_storage_index_(x+1, y+1, z, epsilon->N[0], epsilon->N[1])];
  epsilon_d=epsilon->storage[get_density_storage_index_(x+1, y, z, epsilon->N[0], epsilon->N[1])];
  epsilon_e=epsilon->storage[get_density_storage_index_(x, y, z+1, epsilon->N[0], epsilon->N[1])];
  epsilon_f=epsilon->storage[get_density_storage_index_(x, y+1, z+1, epsilon->N[0], epsilon->N[1])];
  epsilon_g=epsilon->storage[get_density_storage_index_(x+1, y+1, z+1, epsilon->N[0], epsilon->N[1])];
  epsilon_h=epsilon->storage[get_density_storage_index_(x+1, y, z+1, epsilon->N[0], epsilon->N[1])];
  /**The possible terms to average*/
  /*Top ones: clockwise*/
  AB = log(epsilon_b/epsilon_a)*(phi_b - phi_a)/epsilon->dx[1];
  CD = log(epsilon_d/epsilon_c)*(phi_d - phi_c)/epsilon->dx[1];
  DA = log(epsilon_a/epsilon_d)*(phi_a - phi_d)/epsilon->dx[0];
  BC = log(epsilon_c/epsilon_b)*(phi_c - phi_b)/epsilon->dx[0];
  /*Bottom ones: counter clockwise*/
  FE = log(epsilon_e/epsilon_f)*(phi_e - phi_f)/epsilon->dx[1];
  HG = log(epsilon_g/epsilon_h)*(phi_g - phi_h)/epsilon->dx[1];
  EH = log(epsilon_h/epsilon_e)*(phi_h - phi_e)/epsilon->dx[0];
  GF = log(epsilon_f/epsilon_g)*(phi_f - phi_g)/epsilon->dx[0];
  /*The z directions*/
  BF = log(epsilon_f/epsilon_b)*(phi_f - phi_b)/epsilon->dx[2];
  EA = log(epsilon_a/epsilon_e)*(phi_a - phi_e)/epsilon->dx[2];
  DH = log(epsilon_h/epsilon_d)*(phi_h - phi_d)/epsilon->dx[2];
  GC = log(epsilon_c/epsilon_g)*(phi_c - phi_g)/epsilon->dx[2];
  fprintf(stderr, "del_phi_term: %ld %ld %ld: AB=%g CD=%g DA=%g BC=%g FE=%g HG=%g EH=%g GF=%g BF=%g EA=%g DH=%g GC=%g\n", x, y, z, AB, CD, DA, BC, FE, HG, EH, GF, BF, EA, DH, GC);
  ex=2*x;
  ox=ex+1;
  /*We're setting the points from A*/
  /*use oz oy ox: central point*/
  del_phi_term->storage[get_density_storage_index_(ox, oy, oz, del_phi_term->N[0], del_phi_term->N[1])]
    = (AB+CD+DA+BC+FE+HG+EH+GF+BF+EA+DH+GC)/12.0;
  /*use oz oy ex: y-z facial*/
  del_phi_term->storage[get_density_storage_index_(ex, oy, oz, del_phi_term->N[0], del_phi_term->N[1])]
    = (AB + BF + FE + EA)/4.0;
  /*use oz ey ox: x-z facial*/
  del_phi_term->storage[get_density_storage_index_(ox, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
    = (DA + EA + EH + DH)/4.0;
  /*use oz ey ex: z axial */
  del_phi_term->storage[get_density_storage_index_(ex, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
    = EA;
  /*use ez oy ox: x-y facial*/
  del_phi_term->storage[get_density_storage_index_(ox, oy, ez, del_phi_term->N[0], del_phi_term->N[1])]
    = (AB + BC + CD + DA)/4.0;
  /*use ez oy ex: y axial*/
  del_phi_term->storage[get_density_storage_index_(ex, oy, ez, del_phi_term->N[0], del_phi_term->N[1])]
    = AB;
  /*use ez ey ox: x axial*/
  del_phi_term->storage[get_density_storage_index_(ex, ey, oz, del_phi_term->N[0], del_phi_term->N[1])]
      = DA;
  /*use ez ey ex: ignored*/
      }
    }
  }
}

/**Sets the effective potential for the simple method.
 *\param del_phi density to be asigned the del_phi term.
 *\param epsilon dielectric constant times epsilon_0
 *\param phi phi from the previous iteration
 *\param periodic true if this is to be periodic boundary conditions.
 */
static void set_del_phi_term_simple(struct density* del_phi, const struct density* epsilon, const struct density* phi, int periodic) {
  /*Set the interior points with symmetric differencing*/
  int x, y, z;
  /**\note the dx terms cancel out (if all dx[0]==dx[1]==dx[2]) when multiplied by the dx^2 from the del^2 phi term, as is the case with rho!
  double inverse_4_delta2_x = 1.0/(4.0*del_phi->dx[0]*del_phi->dx[0]);
  double inverse_4_delta2_y = 1.0/(4.0*del_phi->dx[1]*del_phi->dx[1]);
  double inverse_4_delta2_z = 1.0/(4.0*del_phi->dx[2]*del_phi->dx[2]);
  */
  double inverse_4_delta2_x = 1.0/4.0;
  double inverse_4_delta2_y = 1.0/4.0;
  double inverse_4_delta2_z = 1.0/4.0;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y, z)
  #endif
  for(z=1; z<del_phi->N[2]-1; z++) {
    for(y=1; y<del_phi->N[1]-1; y++) {
      for(x=1; x<del_phi->N[0]-1; x++) {
  *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
  //fprintf(stderr, "%d %d %d: %g", x, y, z, *get_density_storage_at_c_(x, y, z, del_phi));
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
  //fprintf(stderr, " %g", *get_density_storage_at_c_(x, y, z, del_phi));
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
  //fprintf(stderr, " %g\n", *get_density_storage_at_c_(x, y, z, del_phi));
      }
    }
  }
  /*Set the z faces*/
  z=0;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y)
  #endif
  for(y=1; y<del_phi->N[1]-1; y++) {
    for(x=1; x<del_phi->N[0]-1; x++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, del_phi->N[2]-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, del_phi->N[2]-1, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z, phi));
      }
    }
  }
  z=del_phi->N[2]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y)
  #endif
  for(y=1; y<del_phi->N[1]-1; y++) {
    for(x=1; x<del_phi->N[0]-1; x++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, 0, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, 0, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      }
    }
  }
  /*Set the y faces*/
  y=0;
  #ifdef OPENMP
  #pragma omp parallel for private(x, z)
  #endif
  for(z=1; z<del_phi->N[2]-1; z++) {
    for(x=1; x<del_phi->N[0]-1; x++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, del_phi->N[1]-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, del_phi->N[1]-1, z, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y, z, phi));
      }
    }
  }
  y=del_phi->N[1]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(x, z)
  #endif
  for(z=1; z<del_phi->N[2]-1; z++) {
    for(x=1; x<del_phi->N[0]-1; x++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_y*log((*get_density_storage_at_c_(x, 0, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, 0, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += 4.0*inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      }
    }
  }
  /*Finally, the x faces*/
  x=0;
  #ifdef OPENMP
  #pragma omp parallel for private(y, z)
  #endif
  for(z=1; z<del_phi->N[2]-1; z++) {
    for(y=1; y<del_phi->N[1]-1; y++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(del_phi->N[0]-1, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(del_phi->N[0]-1, y, z, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x+1, y, z, epsilon))/(*get_density_storage_at_c_(x, y, z, epsilon)))*(*get_density_storage_at_c_(x+1, y, z, phi) - *get_density_storage_at_c_(x, y, z, phi));
      }
    }
  }
  x=del_phi->N[0]-1;
  #ifdef OPENMP
  #pragma omp parallel for private(y, z)
  #endif
  for(z=1; z<del_phi->N[2]-1; z++) {
    for(y=1; y<del_phi->N[1]-1; y++) {
      *get_density_storage_at_(x, y, z, del_phi) = inverse_4_delta2_y*log((*get_density_storage_at_c_(x, y+1, z, epsilon))/(*get_density_storage_at_c_(x, y-1, z, epsilon)))*(*get_density_storage_at_c_(x, y+1, z, phi) - *get_density_storage_at_c_(x, y-1, z, phi));
      *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_z*log((*get_density_storage_at_c_(x, y, z+1, epsilon))/(*get_density_storage_at_c_(x, y, z-1, epsilon)))*(*get_density_storage_at_c_(x, y, z+1, phi) - *get_density_storage_at_c_(x, y, z-1, phi));
      if(periodic) {
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(0, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(0, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      }else{
  *get_density_storage_at_(x, y, z, del_phi) += inverse_4_delta2_x*log((*get_density_storage_at_c_(x, y, z, epsilon))/(*get_density_storage_at_c_(x-1, y, z, epsilon)))*(*get_density_storage_at_c_(x, y, z, phi) - *get_density_storage_at_c_(x-1, y, z, phi));
      }
    }
  }
}


struct density* get_potential_for_density_inhomog_simple(const struct density* rho, const struct density* epsilon, double close_enough, int global, int save_iterations, GError **e) {
  fprintf(stderr, "get_potential_for_density_inhomog_simple: starting iteration.\n");
  /*Get our initial phi*/
  GError *err = NULL;
  struct density* previous_phi = get_potential_for_density(rho, epsilon, NULL, NULL, NULL, &err);
  if((previous_phi == NULL) || (err != NULL)) {
    if(err != NULL) {
      g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential: ");
    }else{
      *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential: no further data because error was NULL; this is definitely a bug in the code! ");
    }
    if(previous_phi != NULL) free_density(previous_phi);
    return NULL;
  } 
  /*Dump the output*/
  int level=0;
  GString *iter_dir = g_string_new("del_phi_iter=0");
  GFile* base_dir=g_file_new_for_path(iter_dir->str);
  g_string_free(iter_dir, TRUE);
  char* filename;
  if(save_iterations) {
    write_density_sweep(Z, previous_phi, &filename, &err, base_dir, FALSE);
    g_free(filename);
  }
  err = NULL;
  g_object_unref(base_dir);
  struct density* phi = NULL;
  struct density* del_phi = new_density(epsilon->N, epsilon->dx, 0);
  double diff;
  do {
    ++level;
    set_del_phi_term_simple(del_phi, epsilon, previous_phi, 1);
    phi = get_potential_for_density(rho, epsilon, NULL, NULL, del_phi, &err);
    if((phi == NULL) || (err != NULL)) {
      if(err != NULL) {
  g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting potential for iteration %d: ", level);
      }else{
  *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_simple: error getting initial potential for iteration %d: no further data because error was NULL; this is definitely a bug in the code! ", level);
      }
      if(phi != NULL) free_density(phi);
      free_density(del_phi);
      free_density(previous_phi);
      return NULL;
    } 
    iter_dir = g_string_new("del_phi_simple_iter=");
    g_string_append_printf(iter_dir, "%d/phi", level);
    base_dir=g_file_new_for_path(iter_dir->str);
    g_string_free(iter_dir, TRUE);
    err = NULL;
    if(save_iterations) {
      write_density_sweep(Z, phi, &filename, &err, base_dir, FALSE);
      g_free(filename);
    }
    err = NULL;
    g_object_unref(base_dir);
    iter_dir = g_string_new("del_phi_simple_iter=");
    g_string_append_printf(iter_dir, "%d/del_phi", level);
    base_dir=g_file_new_for_path(iter_dir->str);
    g_string_free(iter_dir, TRUE);
    err = NULL;
    if(save_iterations) {
      write_density_sweep(Z, del_phi, &filename, &err, base_dir, FALSE);
      g_free(filename);
    }
    err = NULL;
    GFile* dumpfile = g_file_get_child(base_dir, "del_phi.dump");
    if(save_iterations) write_density_slice(Z, NULL, del_phi, &filename, &err, dumpfile, ' ', "nm", "V");
    err = NULL;
    g_object_unref(base_dir);
    fprintf(stderr, "@%d:\n", level);
    if(global) {
      diff = density_get_max_abs_delta_pct(previous_phi, phi, 1, 0);
    }else{
      diff = density_get_max_abs_delta_pct(previous_phi, phi, 0, 0);
    }
    fprintf(stderr, "\tdiff=%g\n", diff);
    free_density(previous_phi);
    previous_phi=phi;
  } while(diff >= close_enough);
  free_density(del_phi);
  /**\note if we converged, then phi is contained in previous_phi from the end of the previous loop!!*/
  return previous_phi;
}

struct density* get_potential_for_density_inhomog_multigrid(const struct density* rho, const struct density* epsilon, double close_enough, GError **e) {
  /*Get the higher-res grid versions; note that the density is LAPACK-based so the grid spacing is implicitly included in the value of the density (f*dx*dy*dz).*/
  /*Get our initial phi*/
  GError *err = NULL;
  struct density* previous_phi = get_potential_for_density(rho, epsilon, NULL, NULL, NULL, &err);
  if((previous_phi == NULL) || (err != NULL)) {
    if(err != NULL) {
      g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential: ");
    }else{
      *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential: no further data because error was NULL; this is definitely a bug in the code! ");
    }
    if(previous_phi != NULL) free_density(previous_phi);
    return NULL;
  } 
  int level=0;
  /*Dump the output*/
  GString *iter_dir = g_string_new("del_phi_multigrid_iter=");
  g_string_append_printf(iter_dir, "%d", level);
  GFile* base_dir=g_file_new_for_path(iter_dir->str);
  g_string_free(iter_dir, TRUE);
  char* filename;
  write_density_sweep(Z, previous_phi, &filename, &err, base_dir, FALSE);
  err = NULL;
  g_free(filename);
  g_object_unref(base_dir);
  struct density* hires_rho = NULL;
  struct density* hires_epsilon = NULL;
  struct density* del_phi_term = NULL;
  struct density* phi = NULL;
  double final_max_pct_diff;
  do {
    ++level;
    if(phi != NULL) {
      /*To save space, start by removing the old del_phi term*/
      free_density(del_phi_term);
      /*Also, update previous_phi now*/
      free_density(previous_phi);
      previous_phi = phi;
      /*And also to save space, only do one density at a time*/
      struct density* hires_rho_ = get_double_resolution_density(rho, 3);
      free_density(hires_rho);
      hires_rho = hires_rho_;
      struct density* hires_epsilon_ = get_double_resolution_density(epsilon, 0);
      free_density(hires_epsilon);
      hires_epsilon = hires_epsilon_;
    }else{
      hires_rho = get_double_resolution_density(hires_rho, 3);
      hires_epsilon = get_double_resolution_density(hires_epsilon, 0);
    }
    /*Allocate and then set the del_phi term*/
    del_phi_term = new_density(hires_rho->N, hires_rho->dx, 0);
    set_del_phi_term(del_phi_term, epsilon, previous_phi);
    /*Get the updated phi*/
    struct density* phi = get_potential_for_density(hires_rho, hires_epsilon, NULL, del_phi_term, NULL, &err);
    if((phi == NULL) || (err != NULL)) {
      if(err != NULL) {
  g_propagate_prefixed_error(e, err, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting potential for iteration %d: ", level);
      }else{
  *e = g_error_new(FFTW_WAVEFUNC_GERROR_DOMAIN, 102, "fftw_wavefunc.c::get_potential_for_density_inhomog_multigrid: error getting initial potential for iteration %d: no further data because error was NULL; this is definitely a bug in the code! ", level);
      }
      if(phi != NULL) free_density(phi);
      free_density(del_phi_term);
      free_density(hires_rho);
      free_density(hires_epsilon);
      free_density(previous_phi);
      return NULL;
    } 
    iter_dir = g_string_new("del_phi_multigrid_iter=");
    g_string_append_printf(iter_dir, "%d/phi", level);
    base_dir=g_file_new_for_path(iter_dir->str);
    g_string_free(iter_dir, TRUE);
    err = NULL;
    write_density_sweep(Z, phi, &filename, &err, base_dir, FALSE);
    err = NULL;
    g_free(filename);
    g_object_unref(base_dir);
    iter_dir = g_string_new("del_phi_multigrid_iter=");
    g_string_append_printf(iter_dir, "%d/del_phi", level);
    base_dir=g_file_new_for_path(iter_dir->str);
    g_string_free(iter_dir, TRUE);
    err = NULL;
    write_density_sweep(Z, del_phi_term, &filename, &err, base_dir, FALSE);
    err = NULL;
    g_free(filename);
    g_object_unref(base_dir);
    /*parallelized finding maximum deviation from previous phi*/
    int N_threads = 0;
    #ifdef OPENMP
    N_threads = omp_get_max_threads();
    #else
    N_threads = 1;
    #endif
    double max_pct_diff[N_threads];
    for(int i=0; i<N_threads; i++) max_pct_diff[i] = nan("NaN");
    long x, y, z;
    double diff;
    int thread;
    double max_previous_phi = density_get_max(previous_phi, 1);
    #ifdef OPENMP
    #pragma omp parallel for private(y, x, diff, thread)
    #endif
    for(z=0; z<previous_phi->N[2]; z++) {
      thread = 0;
      #ifdef OPENMP
      thread = omp_get_thread_num();
      #endif
      for(y=0; y<previous_phi->N[1]; y++) {
  for(x=0; x<previous_phi->N[0]; x++) {
    diff = abs(((*get_density_storage_at_c_(2*x, 2*y, 2*z, phi))-(*get_density_storage_at_c_(x, y, z, previous_phi)))/max_previous_phi);
    //we use !<= since NaN will evaluate to false on any comparison
    if(!(diff <= max_pct_diff[thread])) {
      max_pct_diff[thread] = diff;
    }
    //fprintf(stderr, "fastmax@%d %ld %ld %ld: (%g - %g)/%g = %g (/->%g)\n", level, x, y, z, *get_density_storage_at_c_(2*x, 2*y, 2*z, phi),  *get_density_storage_at_c_(x, y, z, previous_phi), max_previous_phi, diff, (*get_density_storage_at_c_(2*x, 2*y, 2*z, phi))/(*get_density_storage_at_c_(x, y, z, previous_phi)));
  }
      }
    }
    final_max_pct_diff = max_pct_diff[0];
    for(int i=1; i<N_threads; i++) {
      if(!(max_pct_diff[i] <= final_max_pct_diff)) final_max_pct_diff = max_pct_diff[i];
    }
    fprintf(stderr, "@%d: Max difference is %g:\n", level, final_max_pct_diff);
  } while(final_max_pct_diff >= close_enough);
  /*free the hires grids*/
  free_density(hires_rho);
  free_density(hires_epsilon);
  free_density(del_phi_term);
  /*Halve the grid resolution of phi until it's the original res*/
  err = NULL;
  for(int i=0; i<level; i++) {
    struct density* phi_ = get_half_resolution_density(phi, &err, 0);
    if(err != NULL) {
      free_density(phi);
      g_propagate_prefixed_error(e, err, "fftw_wavefunc::get_potential_for_density_inhomog@level %d: error occurred halving the resolution of a phi we got from a recursive call:\n\t", level);
      return NULL;
    }
    free_density(phi);
    phi = phi_;
  }
  return phi;
}