lp_wavefunction.h File Reference

#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <libdotcode/function_parser.h>
#include <libdotcode/binary_data_common.h>
#include <libpostproc/density_common.h>
#include <libdotcode/function_parser.h>
#include <libdotcode/binary_data_common.h>
#include <libpostproc/density_common.h>

Include dependency graph for lp_wavefunction.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	wavefunction
struct	density
Enumerations
enum	coord { X, Y, Z }
Functions
struct density *	get_double_resolution_density (const struct density *original, int auto_scale)
struct density *	get_half_resolution_density (const struct density *original, GError **e, int auto_scale)
struct density *	get_density_for_list_of_wavefunctions (GList *list, GError **err)
void	load_wavefunctions (GList *wf_list, GError **err)
void	init_density (struct density *d, double initial_value)
void	add_to_density (struct density *d, struct wavefunction *wf)
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)
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)
void	write_density_sweep (enum coord c, const struct density *d, char **filename, GError **err, GFile *location, gboolean sweep_log)
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)
void	density_mul_by_const (struct density *d, double c)
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)
void	write_density_povray (const struct density *d, GFile *towrite, int pixel_strategy, GError **err)
struct density *	expand_density_box_with_value (const struct density *od, const int *I, double value, GError **err)
struct density *	expand_density_box (const struct density *od, const int *I, GError **err)
struct density *	new_density (const int *N, const double *dx, double intial_value)
struct density *	new_density_copy (const struct density *orig)
void	free_density (struct density *d)
int	op_density_to_density (struct density *to, enum density_op op, const struct density *from)
void	density_execute_callback (const struct density *d, void(*callback)(int x, int y, int z, const struct density *d, void *privdat), void *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)
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)
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)
double	density_get_max_abs_delta_pct (const struct density *d1, const struct density *d2, int scaleby, double scaling)
double	density_get_min_abs_delta_pct (const struct density *d1, const struct density *d2, int scaleby, double scaling)
double	density_get_max (const struct density *d, int abs)
double	density_get_max_and_loc (const struct density *d, int abs, int *location)
double	density_get_min (const struct density *d, int abs)
double	density_get_min_and_loc (const struct density *d, int abs, int *location)
struct dotcode_grid *	dotcode_grid_from_density (const struct density *d, GError **err)
long unsigned int	get_density_storage_index_ (unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny)
double *	get_density_storage_at_ (unsigned x, unsigned y, unsigned z, struct density *d)
const double *	get_density_storage_at_c_ (unsigned x, unsigned y, unsigned z, const struct density *d)
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)
double	density_integrate_spline_parallel (const struct density *d, double(*f)(double x, double y, double z, void *privdat), void **privdat)
double	density_integrate_sum (const struct density *d, double(*f)(double x, double y, double z, void *privdat), void *privdat)
double	density_integrate_sum_parallel (const struct density *d, double(*f)(double x, double y, double z, void *privdat), void **privdat)

---

Enumeration Type Documentation

enum coord

Enumerator:

X
Y
Z

Definition at line 42 of file lp_wavefunction.h.

{X, Y, Z};

---

Function Documentation

void add_to_density	(	struct density *	d,
		struct wavefunction *	wf
	)

Definition at line 148 of file lp_wavefunction.c.

References internal_add_to_density().

                                                                {
  internal_add_to_density(d, wf);
}

Here is the call graph for this function:

void density_2_execute_callback	(	const struct density *	d1,
		const struct density *	d2,
		void(*)(int x, int y, int z, const struct density *d1, const struct density *d2, void *privdat)	callback,
		void *	privdat
	)

same, but takes two density functions

Definition at line 1280 of file lp_wavefunction.c.

References density::N.

Referenced by density_get_max_abs_delta_pct(), and density_get_min_abs_delta_pct().

                                                                                                                                                                                                             {
  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);
      }
    }
  }
}

Here is the caller graph for this function:

void density_2_execute_callback_parallel	(	const struct density *	d1,
		const struct density *	d2,
		void(*)(int x, int y, int z, const struct density *d1, const struct density *d2, void *privdat)	callback,
		void **	privdat,
		int *	x_report
	)

Definition at line 1317 of file lp_wavefunction.c.

References density::N.

Referenced by density_get_max_abs_delta_pct(), and density_get_min_abs_delta_pct().

                                                                                                                                                                                                                                      {
  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
     );
      }
    }
  }
}

Here is the caller graph for this function:

void density_execute_callback	(	const struct density *	d,
		void(*)(int x, int y, int z, const struct density *d, void *privdat)	callback,
		void *	privdat
	)

Executes the specified callback

Parameters:

d	density struct
callback	pointer to function to call
privdat	private data for the callback

Returns:

nothing

Definition at line 1270 of file lp_wavefunction.c.

References density::N.

                                                                                                                                                    {
  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_execute_callback_parallel	(	const struct density *	d,
		void(*)(int x, int y, int z, const struct density *d, void *privdat)	callback,
		void **	privdat,
		int *	x_report
	)

OpenMP-parallelized callback executors. This is identical to the serial executors above, but parallelized with OpenMP. The most obvious effect is a slight change in the API: the executor itself is passed an *array* of pointers to void. The callback itself still accepts a pointer to void (not a pointer to a pointer to void). The specified callback is executed in an OpenMP parallelized for loop (on the z index, for whatever that matters). Each thread is assigned (by thread number) a pointer to void to use as its privdat argument to the callback. Each thread, then, executes the callback on its specified points (subject to OpenMP scheduling; no guarantees can be made with regard to order) with the same callback privdat pointer. See the density_get_max function for an example of how to use this.

Parameters:

d	density struct
callback	pointer to function to call
privdat	array of private data pointers for the callback, one per thread.
x_report	array of ints, one per thread. Set to 1 if a thread is executed, and is left untouched otherwise.

Returns:

nothing

Definition at line 1290 of file lp_wavefunction.c.

References density::N.

Referenced by density_get_max_and_loc(), and density_get_min_and_loc().

                                                                                                                                                                             {
  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
     );
      }
    }
  }
}

Here is the caller graph for this function:

double density_get_max	(	const struct density *	d,
		int	abs
	)

Gets the maximum value in the density

Parameters:

d	density to inspect
abs	if true, compare absolute value instead of raw value.

Returns:

the maximum value in d

Definition at line 1592 of file lp_wavefunction.c.

References density_get_max_and_loc(), and minmax_dat::location.

Referenced by density_get_max_abs_delta_pct(), and get_potential_for_density_inhomog_multigrid().

                                                         {
  int location[3];
  return density_get_max_and_loc(d, abs, location);
}

Here is the call graph for this function:

Here is the caller graph for this function:

double density_get_max_abs_delta_pct	(	const struct density *	d1,
		const struct density *	d2,
		int	scaleby,
		double	scaling
	)

returns the maximum percent (at the point) difference between two densities, or NaN if the dimensions differ

Parameters:

d1	density to compare
d2	density to compare
scaleby	if 3, use scaling to scale before comparing (the difference will be multiplied by this value) if 2, get diff as percent of global maximum of d2 if 1, get diff as percent of global maximum of d1 if 0, get diff as percent of local value in d1. if -1, get diff as percent of local value in d2.
scaling	value to use for scaling if scaleby==3, else ignored

Definition at line 1385 of file lp_wavefunction.c.

References density_2_execute_callback(), density_2_execute_callback_parallel(), density_get_max(), get_max_diff_d1(), get_max_diff_d2(), get_max_diff_scaled(), diff_data::minmax, density::N, and diff_data::scaling.

Referenced by get_potential_for_density_inhomog_simple().

                                                                                                                      {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double density_get_max_and_loc	(	const struct density *	d,
		int	abs,
		int *	location
	)

Note:

that max[i] is already one or zero.

Definition at line 1596 of file lp_wavefunction.c.

References density_execute_callback_parallel(), get_max(), get_max_abs(), and minmax_dat::minmax.

Referenced by density_get_max(), and main().

                                                                                {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double density_get_min	(	const struct density *	d,
		int	abs
	)

Gets the minimum value in the density

Parameters:

d	density to inspect
abs	if true, compare absolute value instead of raw value.

Returns:

the minimum value in d

Definition at line 1627 of file lp_wavefunction.c.

References density_get_min_and_loc(), and minmax_dat::location.

Referenced by density_get_min_abs_delta_pct().

                                                         {
  int location[3];
  return density_get_min_and_loc(d, abs, location);
}

Here is the call graph for this function:

Here is the caller graph for this function:

double density_get_min_abs_delta_pct	(	const struct density *	d1,
		const struct density *	d2,
		int	scaleby,
		double	scaling
	)

returns the minimum percent (at the point) difference between two densities, or NaN if the dimensions differ

Parameters:

d1	density to compare
d2	density to compare
scaleby	if 3, use scaling to scale before comparing (the difference will be multiplied by this value) if 2, get diff as percent of global maximum of d2 if 1, get diff as percent of global maximum of d1 if 0, get diff as percent of local value in d1. if -1, get diff as percent of local value in d2.
scaling	value to use for scaling if scaleby==3, else ignored

Definition at line 1467 of file lp_wavefunction.c.

References density_2_execute_callback(), density_2_execute_callback_parallel(), density_get_min(), get_min_diff_d1(), get_min_diff_d2(), get_min_diff_scaled(), diff_data::minmax, density::N, and diff_data::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;
}

Here is the call graph for this function:

double density_get_min_and_loc	(	const struct density *	d,
		int	abs,
		int *	location
	)

Note:

that min[i] is already one or zero.

Definition at line 1631 of file lp_wavefunction.c.

References density_execute_callback_parallel(), get_min(), get_min_abs(), and minmax_dat::minmax.

Referenced by density_get_min().

                                                                                {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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(*)(double x, double y, double z, void *privdat)	f,
		void *	privdat
	)

Evaluates an integral of the density using GSL and a spline fit.

Parameters:

z_error_estimate,:	pointer to a double in which to store the estimated integration error from the final integration along z, or NULL if the error isn't to be saved.
y_error_estimate,:	pointer to an array (of size wf->N[2], one for each z point) of doubles in which to store the estimated integration error for each y integral, or NULL if the error isn't to be saved.
x_error_estimate,:	pointer to an array (of size wf->N[1]*wf->N[2], one for each (y, z) pairing) of doubles in which to store the estimated error for each x integral, or NULL if the error isn't to be saved.
abserr	the requested absolute error (see GSL documentation)
relerr	the requested relative error (see GSL documentation)
d	the density over which to integrate
f	a function to call at each point (if NULL, does nothing).
privdat	private data to pass to f.

Returns:

the integral

Note:

d will be fit to a cubic spline before integration.

THIS IS BROKEN

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.

Integrating x results in a yz plane of values

Note:

that the y coords are contiguous in memory

Definition at line 1714 of file lp_wavefunction.c.

References gsl_integrand_x::accel, density::dx, gsl_integrand_x::f, get_density_storage_at_c(), gsl_integrand_x_func(), gsl_integrand_x_nofunc(), gsl_integrand_yz(), gsl_integrand_x::iter, MAX_GSL_INTEGRATION_INTERVALS, density::N, gsl_integrand_x::outfile, gsl_integrand_x::privdat, gsl_integrand_x::spline, gsl_integrand_x::y, and gsl_integrand_x::z.

                                                                                                                                                                                                                                              {
   /**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;
 }

Here is the call graph for this function:

double density_integrate_spline_parallel	(	const struct density *	d,
		double(*)(double x, double y, double z, void *privdat)	f,
		void **	privdat
	)

Evaluates an integral of the density using GSL and a spline fit, parallelized with OpenMP on the z coordinate

Parameters:

d	the density over which to integrate
f	a function to call at each point (if NULL, does nothing).
privdat	array of private data to pass to the callback (one for each thread, indexed by omp_get_thread_num()).

Returns:

the integral

Note:

d will be fit to a cubic spline before integration.

the same f can be used for both the single- and multi-threaded functions.

if privdat is NULL, NULL will be passed to f (this saves you having to allocate an array of NULLs

THIS IS SUPER-BROKEN

double density_integrate_sum	(	const struct density *	d,
		double(*)(double x, double y, double z, void *privdat)	f,
		void *	privdat
	)

Evaluates an integral of the density using GSL and a simple summing mechanism: sum goes inwards from both ends, assuming that each end is a minimum.

Parameters:

d	the density over which to integrate
f	a function to call at each point (if NULL, does nothing).
privdat	array of private data to pass to the callback (one for each thread, indexed by omp_get_thread_num()).

Returns:

the integral

Note:

d will be fit to a cubic spline before integration.

the same f can be used for both the single- and multi-threaded functions.

if privdat is NULL, NULL will be passed to f (this saves you having to allocate an array of NULLs

Definition at line 1825 of file lp_wavefunction.c.

References density::dx, gsl_integrand_x::f, get_density_storage_at_c_(), density::N, gsl_integrand_x::y, and gsl_integrand_x::z.

                                                                                                                               {
  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;
}

Here is the call graph for this function:

double density_integrate_sum_parallel	(	const struct density *	d,
		double(*)(double x, double y, double z, void *privdat)	f,
		void **	privdat
	)

Definition at line 1843 of file lp_wavefunction.c.

References density::dx, gsl_integrand_x::f, get_density_storage_at_c_(), density::N, gsl_integrand_x::y, and gsl_integrand_x::z.

Referenced by dotcode_wavefunction_get_percentage_in_band().

                                                                                                                                         {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void density_mul_by_const	(	struct density *	d,
		double	c
	)

Definition at line 276 of file lp_wavefunction.c.

References c, get_density_storage_index(), density::N, and density::storage.

Referenced by density_un_normalize(), get_double_resolution_density(), get_half_resolution_density(), get_potential_for_density(), and main().

                                                       {
  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;
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct dotcode_grid* dotcode_grid_from_density	(	const struct density *	d,
		GError **	err
	)		[read]

Definition at line 723 of file lp_wavefunction.c.

References dotcode_tensor::data, dotcode_grid::data, DOTCODE_ARRAY, DOTCODE_FLOAT8, dotcode_free_type(), density::dx, e, ev_to_hartree(), get_density_storage_index(), dotcode_grid::gridsites, dotcode_tensor::indices, dotcode_tensor::metadata, density::N, dotcode_tensor::Nd, dotcode_grid::Ng, density::storage, dotcode_tensor::t, unlikely, and WAVEFUNC_GERROR_DOMAIN.

Referenced by main().

                                                                                      {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* expand_density_box	(	const struct density *	od,
		const int *	I,
		GError **	err
	)		[read]

Expands a density box using the value from the outermost grid site.

Parameters:

od	the original density
I	a 3-int array containing the number of grid sites to expand outwards

Definition at line 903 of file lp_wavefunction.c.

References expand_density_box_with_value(), get_density_storage_index(), density::N, density::storage, and WAVEFUNC_GERROR_DOMAIN.

Referenced by main().

                                                                                         {
  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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* expand_density_box_with_value	(	const struct density *	od,
		const int *	I,
		double	value,
		GError **	err
	)		[read]

Expands a density box using material mat

Parameters:

od	the original density
I	a 3-int array indicating how many grid sites to expand outward.
value	a double containing the value of the density throughout the new region

Returns:

a pointer to a new density.

Definition at line 814 of file lp_wavefunction.c.

References density::dx, get_density_storage_index(), density::N, density::storage, and WAVEFUNC_GERROR_DOMAIN.

Referenced by expand_density_box(), and main().

                                                                                                                  {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void free_density

(

struct density *

d

)

Definition at line 717 of file lp_wavefunction.c.

References density::storage.

Referenced by dotcode_wavefunction_get_percentage_in_band(), get_potential_for_density(), get_potential_for_density_inhomog_multigrid(), get_potential_for_density_inhomog_simple(), and main().

                                     {
  if(d == NULL) return;
  if(d->storage != NULL) free(d->storage);
  free(d);
}

Here is the caller graph for this function:

struct density* get_density_for_list_of_wavefunctions	(	GList *	list,
		GError **	err
	)		[read]

Gets a density from a list of *loaded* wavefunctions

Parameters:

list	the list of *loaded* wavefunction structures.
err	pointer to a pointer in which to store the GError (if there's an error)

Returns:

a density that is the sum of the individual wavefunctions.

Note:

Wavefunctions must be *loaded* before calling this function!

The charge is used! If you just want pre-normalized |psi|^2, set charge to 1.

Definition at line 152 of file lp_wavefunction.c.

References wavefunction::dx, density::dx, wavefunction::file, init_density(), internal_add_to_density(), density::N, wavefunction::N, density::storage, and WAVEFUNC_GERROR_DOMAIN.

Referenced by main().

                                                                                 {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

Gets a block of the meta_density back out.

Parameters:

md	the meta denstiy
d	an initialized density the sizeo f one block
buffzone	the buffer zone size around the meta_density
Mn	the number of neighbors in each direction
x	the x-coordinate of the block
y	the y-coordinate of the block
z	the z-coordinate of the block
offsets	offset in grid sites between neighborsa

Gets a block of the meta_density back out.

Parameters:

md	the meta denstiy
d	an initialized density the sizeo f one block
buffzone	the buffer zone size around the meta_density
Mn	the number of neighbors in each direction
x	the x-coordinate of the block
y	the y-coordinate of the block
z	the z-coordinate of the block
offsets	offset between neighbors in each dimension

Definition at line 251 of file lp_wavefunction.c.

References get_density_storage_index(), density::N, and density::storage.

Referenced by main().

                                                                                                                                                                {
  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])];
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double* get_density_storage_at_	(	unsigned	x,
		unsigned	y,
		unsigned	z,
		struct density *	d
	)

gets pointer to the cell at this location

Parameters:

x	coordinate
y	coordinate
z	coordinate
d	density structure

Returns:

pointer to the density value at the specified point

Note:

no error checking is performed

Definition at line 103 of file lp_wavefunction.c.

References get_density_storage_at().

Referenced by dotcode_wavefunction_get_density_in_band(), and set_del_phi_term_simple().

                                                                                       {
  return get_density_storage_at(x, y, z, d);
}

Here is the call graph for this function:

Here is the caller graph for this function:

const double* get_density_storage_at_c_	(	unsigned	x,
		unsigned	y,
		unsigned	z,
		const struct density *	d
	)

gets const pointer to the cell at this location

Parameters:

x	coordinate
y	coordinate
z	coordinate
d	density structure

Returns:

const pointer to the density value at the specified point

Note:

no error checking is performed

Definition at line 111 of file lp_wavefunction.c.

References get_density_storage_at_c().

Referenced by density_integrate_sum(), density_integrate_sum_parallel(), get_double_resolution_density(), get_potential_for_density(), get_potential_for_density_inhomog_multigrid(), main(), and set_del_phi_term_simple().

                                                                                                     {
  return get_density_storage_at_c(x, y, z, d);
}

Here is the call graph for this function:

Here is the caller graph for this function:

long unsigned int get_density_storage_index_	(	unsigned int	x,
		unsigned int	y,
		unsigned int	z,
		unsigned int	Nx,
		unsigned int	Ny
	)

Definition at line 95 of file lp_wavefunction.c.

References get_density_storage_index().

Referenced by main(), set_del_phi_term(), and write_density_povray().

                                                                                                                               {
  return get_density_storage_index(x, y, z, Nx, Ny);
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* get_double_resolution_density	(	const struct density *	original,
		int	auto_scale
	)		[read]

Gets a new density with a double-resolution grid, using an average for interpolation

Parameters:

original	original density
auto_scale	flags on which dimension to use to scale: 1->Y 2->Y 4->Z (combine for multiple directions)

Returns:

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

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))

Definition at line 1006 of file lp_wavefunction.c.

References c, density_mul_by_const(), density::dx, e, get_density_storage_at_c_(), get_density_storage_index(), density::N, new_density(), and density::storage.

Referenced by get_potential_for_density_inhomog_multigrid(), and main().

                                                                                              {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* get_half_resolution_density	(	const struct density *	original,
		GError **	e,
		int	auto_scale
	)		[read]

Gets a new density with a half-resolution grid, using an average for assignment

Parameters:

original	original density
e	pointer to pointer to Gerror in which to store the error information
auto_scale	flags on which dimension to use to scale: 1->Y 2->Y 4->Z (combine for multiple directions)

Returns:

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

Gets a new density with a half-resolution grid, using an average for assignment

Parameters:

original

original density

Returns:

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

Definition at line 1221 of file lp_wavefunction.c.

References density_mul_by_const(), density::dx, get_density_storage_index(), density::N, new_density(), density::storage, and WAVEFUNC_GERROR_DOMAIN.

Referenced by get_potential_for_density_inhomog_multigrid(), and main().

                                                                                                        {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void init_density	(	struct density *	d,
		double	initial_value
	)

Definition at line 117 of file lp_wavefunction.c.

References density::dx, get_density_storage_index(), density::N, and density::storage.

Referenced by get_density_for_list_of_wavefunctions().

                                                   {
  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;
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void load_wavefunctions	(	GList *	wf_list,
		GError **	err
	)

Loads the wavefunctions from the respective GFiles.

Parameters:

wf_list	List of wavefunctions
err	pointer to a pointer in which to store the GError (if there's an error)

Definition at line 187 of file lp_wavefunction.c.

References wavefunction::bands, wavefunc_header::dx, wavefunction::dx, wavefunc_header::dy, wavefunc_header::dz, e, wavefunction::file, wavefunction::header, wavefunction::N, wavefunc_header::Nc, wavefunc_header::Nx, wavefunc_header::Ny, wavefunc_header::Nz, read_wavefunction_gfile(), wavefunction::storage, and WAVEFUNC_GERROR_DOMAIN.

Referenced by main().

                                                      {
  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;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* new_density	(	const int *	N,
		const double *	dx,
		double	intial_value
	)		[read]

Definition at line 616 of file lp_wavefunction.c.

References density::dx, get_density_storage_index(), density::N, and density::storage.

Referenced by dotcode_wavefunction_get_density_in_band(), get_double_resolution_density(), get_half_resolution_density(), get_potential_for_density(), get_potential_for_density_inhomog_multigrid(), get_potential_for_density_inhomog_simple(), and main().

                                                                                  {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

struct density* new_density_copy

(

const struct density *

orig

)

[read]

Definition at line 644 of file lp_wavefunction.c.

References density::dx, get_density_storage_index(), density::N, and density::storage.

Referenced by get_potential_for_density().

                                                             {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int op_density_to_density	(	struct density *	to,
		enum density_op	op,
		const struct density *	from
	)

Definition at line 673 of file lp_wavefunction.c.

References ADD_EQUALS, ASSIGN, DIV_EQUALS, density::dx, get_density_storage_at(), get_density_storage_at_c(), get_density_storage_index(), MINUS_EQUALS, MUL_EQUALS, density::N, and density::storage.

Referenced by get_potential_for_density().

                                                                                              {
  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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

Puts a density into *one section* of a metadensity.

Parameters:

md	the target meta_density
d	the source density
Nn	an array of ints representing the number of neighbors in each direction.
x	location of this block to be inserted.
y	location of this block to be inserted.
z	location of this block to be inserted.
invert_x	if true, the x-coordinates are read backwards from the source.
invert_y	if true, the y-coordinates are read backwards from the source.
invert_z	if true, the z-coordinates are read backwards from the source.
offsets	offset in grid sites between neighborsa

Definition at line 214 of file lp_wavefunction.c.

References get_density_storage_index(), density::N, and density::storage.

Referenced by main().

                                                                                                                                                                                                    {
  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])];
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

Definition at line 1875 of file lp_wavefunction.c.

References density::dx, e, gsl_integrand_x::f, get_density_storage_index(), density::N, density::storage, str, X, Y, and Z.

Referenced by write_slices_n_sweeps().

                                                                                                                                                                                                         {
  /*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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void write_density_povray	(	const struct density *	d,
		GFile *	towrite,
		int	pixel_strategy,
		GError **	err
	)

Writes the density into a povray visualization.

Parameters:

d	the density to be written out
towrite	the file into which the povray visualization should be written err pointer to a pointer to a GError, used to "return" a GError pointer if an error ocurred.

Definition at line 921 of file lp_wavefunction.c.

References density::dx, e, get_density_storage_index_(), get_max_density_value(), density::N, NUM_STEPS, density::storage, and THRESHOLD.

Referenced by main(), and write_slices().

                                                                                                     {
  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;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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
	)

Writes out a single, ascii density slice or the entire box.

Parameters:

c	the coordinate to output
index	pointer to an int containing the index which should be written out
d	the density to be written
filename	pointer to a pointer to a char, to provide the filename of the written file.
err	pointer to a pointer to a GError, for returning any error info
location	base directory to output, or NULL if cwd. If index is NULL, this is the file to write the dump into.
sepchar	character to use to separate the columns in the output

Returns:

nothing

Definition at line 301 of file lp_wavefunction.c.

References density::dx, e, get_density_storage_index(), density::N, density::storage, str, X, Y, and Z.

Referenced by get_potential_for_density_inhomog_simple(), main(), write_slices(), and write_slices_n_sweeps().

                                                                                                                                                                                                    {
  /*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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void write_density_sweep	(	enum coord	c,
		const struct density *	d,
		char **	filename,
		GError **	err,
		GFile *	location,
		gboolean	sweep_log
	)