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

Include dependency graph for lp_wavefunction.c:

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

Go to the source code of this file.

Classes
struct	diff_data
struct	minmax_dat
struct	gsl_integrand_x
Defines
#define	NUM_STEPS 12
#define	THRESHOLD 0.5
#define	WAVEFUNC_GERROR_DOMAIN 4096
#define	MAX_PIXEL_VAL 255
#define	MAX_GSL_INTEGRATION_INTERVALS 4096
Functions
hotfunc ___const ___always_inline long unsigned int	get_wavefunction_storage_index (unsigned int x, unsigned int y, unsigned int z, unsigned int im, unsigned int b, unsigned int Nx, unsigned int Ny, int Nb)
hotfunc ___const ___always_inline long unsigned int	get_density_storage_index (unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny)
hotfunc ___const ___always_inline long unsigned int	get_density_slice_storage_index (enum coord c, unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny)
long unsigned int	get_density_storage_index_ (unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny)
hotfunc ___always_inline double *	get_density_storage_at (unsigned x, unsigned y, unsigned z, struct density *d)
double *	get_density_storage_at_ (unsigned x, unsigned y, unsigned z, struct density *d)
hotfunc ___always_inline const double *	get_density_storage_at_c (unsigned x, unsigned y, unsigned z, const struct density *d)
const double *	get_density_storage_at_c_ (unsigned x, unsigned y, unsigned z, const struct density *d)
void	init_density (struct density *d, double value)
void	internal_add_to_density (struct density d, struct wavefunction wf)
void	add_to_density (struct density d, struct wavefunction wf)
struct density *	get_density_for_list_of_wavefunctions (GList list, GError *err)
void	load_wavefunctions (GList wf_list, GError *err)
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	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	density_mul_by_const (struct density *d, double c)
void	density_un_normalize (struct density *d)
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_slice (enum coord c, const int index, const struct density d, char filename, GError err, GFile *location, double maxval, int pixel_strategy)
double	get_max_density_value (const struct density *d)
void	write_density_sweep (enum coord c, const struct density d, char filename, GError err, GFile location, int pixel_strategy)
struct density *	new_density (const int N, const double dx, double initial_value)
struct density *	new_density_copy (const struct density *orig)
int	op_density_to_density (struct density to, enum density_op op, const struct density from)
void	free_density (struct density *d)
struct dotcode_grid *	dotcode_grid_from_density (const struct density d, 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)
void	write_density_povray (const struct density d, GFile towrite, int pixel_strategy, GError **err)
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)
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)
static void	get_max_diff_scaled (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
static void	get_min_diff_scaled (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
static void	get_max_diff_d1 (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
static void	get_min_diff_d1 (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
static void	get_max_diff_d2 (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
static void	get_min_diff_d2 (int x, int y, int z, const struct density d1, const struct density d2, void *privdat)
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)
static void	get_max (int x, int y, int z, const struct density d, void privdat)
static void	get_min (int x, int y, int z, const struct density d, void privdat)
static void	get_max_abs (int x, int y, int z, const struct density d, void privdat)
static void	get_min_abs (int x, int y, int z, const struct density d, void privdat)
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)
hotfunc double	gsl_integrand_x_func (double x, void *privdat)
hotfunc double	gsl_integrand_x_nofunc (double x, void *privdat)
hotfunc double	gsl_integrand_yz (double coord, void *privdat)
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_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)
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)

Define Documentation

#define MAX_GSL_INTEGRATION_INTERVALS 4096

Definition at line 52 of file lp_wavefunction.c.

Referenced by density_integrate_spline().

#define MAX_PIXEL_VAL 255

Definition at line 50 of file lp_wavefunction.c.

Referenced by write_density_sweep_slice().

#define NUM_STEPS 12

Definition at line 45 of file lp_wavefunction.c.

Referenced by write_density_povray().

#define THRESHOLD 0.5

Definition at line 46 of file lp_wavefunction.c.

Referenced by write_density_povray().

#define WAVEFUNC_GERROR_DOMAIN 4096

Definition at line 48 of file lp_wavefunction.c.

Referenced by dotcode_grid_from_density(), expand_density_box(), expand_density_box_with_value(), get_density_for_list_of_wavefunctions(), get_half_resolution_density(), load_wavefunctions(), and write_density_sweep().

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
	)

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_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:

void density_un_normalize ( struct density * d )

The wavefunction itself is normalized such that adding all the values togetheer yields 1. We need the density at each point; hence must divide by dx*dy*dz

Definition at line 295 of file lp_wavefunction.c.

References density_mul_by_const(), and density::dx.

                                             {
  double volume = 1.0/((d->dx[0])*(d->dx[1])*(d->dx[2]));
  density_mul_by_const(d, volume);
}

Here is the call 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 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:

hotfunc ___const ___always_inline long unsigned int get_density_slice_storage_index	(	enum coord	c,
		unsigned int	x,
		unsigned int	y,
		unsigned int	z,
		unsigned int	Nx,
		unsigned int	Ny
	)

Definition at line 78 of file lp_wavefunction.c.

References X, Y, and Z.

                                                                                                                                                                                    {
  switch(c) {
  case X:
    x=0;
    Nx=1;
    break;
  case Y:
    y=0;
    Ny=1;
    break;
  case Z:
    z=0;
    break;
  }
  return x + Nx*(y + Ny*(z));
}

hotfunc ___always_inline double* get_density_storage_at	(	unsigned	x,
		unsigned	y,
		unsigned	z,
		struct density *	d
	)

Definition at line 99 of file lp_wavefunction.c.

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

Referenced by get_density_storage_at_(), and op_density_to_density().

                                                                                                               {
  return &(d->storage[get_density_storage_index(x, y, z, d->N[0], d->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:

hotfunc ___always_inline const double* get_density_storage_at_c	(	unsigned	x,
		unsigned	y,
		unsigned	z,
		const struct density *	d
	)

Definition at line 107 of file lp_wavefunction.c.

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

Referenced by density_integrate_spline(), get_density_storage_at_c_(), get_max(), get_max_abs(), get_min(), get_min_abs(), and op_density_to_density().

                                                                                                                             {
  return &(d->storage[get_density_storage_index(x, y, z, d->N[0], d->N[1])]);
}

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:

hotfunc ___const ___always_inline long unsigned int get_density_storage_index	(	unsigned int	x,
		unsigned int	y,
		unsigned int	z,
		unsigned int	Nx,
		unsigned int	Ny
	)

Definition at line 75 of file lp_wavefunction.c.

Referenced by density_mul_by_const(), dotcode_grid_from_density(), expand_density_box(), expand_density_box_with_value(), get_density_from_metadensity(), get_density_storage_at(), get_density_storage_at_c(), get_density_storage_index_(), get_double_resolution_density(), get_half_resolution_density(), get_max_density_value(), get_max_diff_d1(), get_max_diff_d2(), get_max_diff_scaled(), get_min_diff_d1(), get_min_diff_d2(), get_min_diff_scaled(), init_density(), internal_add_to_density(), new_density(), new_density_copy(), op_density_to_density(), put_density_into_metadensity(), write_density_line(), write_density_slice(), and write_density_sweep_slice().

                                                                                                                                                                {
  return x + Nx*(y + Ny*(z));
}

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

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:

static void get_max	(	int	x,
		int	y,
		int	z,
		const struct density *	d,
		void *	privdat
	)		`[static]`

Definition at line 1554 of file lp_wavefunction.c.

References get_density_storage_at_c(), minmax_dat::location, and minmax_dat::minmax.

Referenced by density_get_max_and_loc().

                                                                                 {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(mmd->minmax >= *get_density_storage_at_c(x, y, z, d))) {
    mmd->minmax = *get_density_storage_at_c(x, y, z, d);
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_max_abs	(	int	x,
		int	y,
		int	z,
		const struct density *	d,
		void *	privdat
	)		`[static]`

Definition at line 1573 of file lp_wavefunction.c.

References get_density_storage_at_c(), minmax_dat::location, and minmax_dat::minmax.

Referenced by density_get_max_and_loc().

                                                                                     {
  struct minmax_dat* mmd = (struct minmax_dat*)privdat;
  if(!(fabs(mmd->minmax) >= fabs(*get_density_storage_at_c(x, y, z, d)))) {
    mmd->minmax = fabs(*get_density_storage_at_c(x, y, z, d));
    mmd->location[0]=x;
    mmd->location[1]=y;
    mmd->location[2]=z;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double get_max_density_value ( const struct density * d )

Definition at line 544 of file lp_wavefunction.c.

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

Referenced by write_density_povray(), and write_density_sweep().

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

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_max_diff_d1	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Note:: that this comparison will fail if diff is a NaN

Definition at line 1362 of file lp_wavefunction.c.

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

Referenced by density_get_max_abs_delta_pct().

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

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_max_diff_d2	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Note:: that this comparison will fail if diff is a NaN

Definition at line 1373 of file lp_wavefunction.c.

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

Referenced by density_get_max_abs_delta_pct().

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

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_max_diff_scaled	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Note:: that this comparison will fail if diff is a NaN

Definition at line 1350 of file lp_wavefunction.c.

References get_density_storage_index(), diff_data::minmax, density::N, diff_data::scaling, and density::storage.

Referenced by density_get_max_abs_delta_pct().

                                                                                                                        {
  struct diff_data* dd = (struct diff_data*)privdat;
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])*dd->scaling);  
  /**\note that this comparison will fail if diff is a NaN */
  if(!(diff <= (*dd->minmax))) (*dd->minmax) = diff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_min	(	int	x,
		int	y,
		int	z,
		const struct density *	d,
		void *	privdat
	)		`[static]`

Definition at line 1563 of file lp_wavefunction.c.

References get_density_storage_at_c(), minmax_dat::location, and minmax_dat::minmax.

Referenced by density_get_min_and_loc().

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

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_min_abs	(	int	x,
		int	y,
		int	z,
		const struct density *	d,
		void *	privdat
	)		`[static]`

Definition at line 1582 of file lp_wavefunction.c.

References get_density_storage_at_c(), minmax_dat::location, and minmax_dat::minmax.

Referenced by density_get_min_and_loc().

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

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_min_diff_d1	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Definition at line 1368 of file lp_wavefunction.c.

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

Referenced by density_get_min_abs_delta_pct().

                                                                                                                    {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])]);  
  //fprintf(stderr, "MIN: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  if(!(diff >= *(double*)privdat)) *(double*)privdat = diff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_min_diff_d2	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Definition at line 1379 of file lp_wavefunction.c.

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

Referenced by density_get_min_abs_delta_pct().

                                                                                                                    {
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])/d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])]);  
  //fprintf(stderr, "MIN: %d %d %d: %g - %g = %g\n", x, y, z, d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])], d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])], diff);
  if(!(diff >= *(double*)privdat)) *(double*)privdat = diff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void get_min_diff_scaled	(	int	x,
		int	y,
		int	z,
		const struct density *	d1,
		const struct density *	d2,
		void *	privdat
	)		`[static]`

Definition at line 1356 of file lp_wavefunction.c.

References get_density_storage_index(), diff_data::minmax, density::N, diff_data::scaling, and density::storage.

Referenced by density_get_min_abs_delta_pct().

                                                                                                                        {
  struct diff_data* dd = (struct diff_data*)privdat;
  double diff = fabs((d1->storage[get_density_storage_index(x, y, z, d1->N[0], d1->N[1])] - d2->storage[get_density_storage_index(x, y, z, d2->N[0], d2->N[1])])*dd->scaling);  
  if(!(diff >= (*dd->minmax))) (*dd->minmax) = diff;
}

Here is the call graph for this function:

Here is the caller graph for this function:

hotfunc ___const ___always_inline long unsigned int get_wavefunction_storage_index	(	unsigned int	x,
		unsigned int	y,
		unsigned int	z,
		unsigned int	im,
		unsigned int	b,
		unsigned int	Nx,
		unsigned int	Ny,
		int	Nb
	)

Gets the index into the storage member for a wavefunction, or for a density

Parameters:

x	x index component
y	y index component
z	z index component
b	band index component
im	1 if imaginary component is wanted; 0 else.
Nx	Number of x sites
Ny	Number of y sites
Nb	Number of bands

Returns:: long unsigned int, which is the address in the storage array.

Note:: there is no error checking, for speed purposes. E.g. if im is 3, it'll just be blindly used as a 3 and you'll be addressing the imaginary component of the next index.; Nz is not needed.

Definition at line 71 of file lp_wavefunction.c.

                                                                                                                                                                                                              {
  return im + 2*(b + Nb*(x + Nx*(y + Ny*(z))));
             }

hotfunc double gsl_integrand_x_func	(	double	x,
		void *	privdat
	)

Definition at line 1675 of file lp_wavefunction.c.

References gsl_integrand_x::accel, gsl_integrand_x::f, gsl_integrand_x::iter, gsl_integrand_x::outfile, gsl_integrand_x::privdat, gsl_integrand_x::spline, gsl_integrand_x::y, and gsl_integrand_x::z.

Referenced by density_integrate_spline().

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

Here is the caller graph for this function:

hotfunc double gsl_integrand_x_nofunc	(	double	x,
		void *	privdat
	)

Definition at line 1683 of file lp_wavefunction.c.

References gsl_integrand_x::accel, gsl_integrand_x::iter, gsl_integrand_x::outfile, gsl_integrand_x::spline, gsl_integrand_x::y, and gsl_integrand_x::z.

Referenced by density_integrate_spline().

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

Here is the caller graph for this function:

hotfunc double gsl_integrand_yz	(	double	coord,
		void *	privdat
	)

Definition at line 1690 of file lp_wavefunction.c.

References gsl_integrand_x::accel, gsl_integrand_x::iter, gsl_integrand_x::outfile, and gsl_integrand_x::spline.

Referenced by density_integrate_spline().

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

Here is the caller graph for this function:

void init_density	(	struct density *	d,
		double	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 internal_add_to_density	(	struct density *	d,
		struct wavefunction *	wf
	)

Definition at line 133 of file lp_wavefunction.c.

References wavefunction::bands, wavefunction::charge, get_density_storage_index(), get_wavefunction_storage_index(), wavefunction::N, density::N, density::storage, and wavefunction::storage.

Referenced by add_to_density(), and get_density_for_list_of_wavefunctions().

                                                                         {
  register int i, j, k, b;
  double re, im;
  for(k=0; k<d->N[2]; k++) {
    for(j=0; j<d->N[1]; j++) {
      for(i=0; i<d->N[0]; i++) {
  for(b=0; b<wf->bands; b++){
    re = wf->storage[get_wavefunction_storage_index(i, j, k, b, 0, wf->N[0], wf->N[1], wf->bands)];
    im = wf->storage[get_wavefunction_storage_index(i, j, k, b, 0, wf->N[0], wf->N[1], wf->bands)];
    d->storage[get_density_storage_index(i, j, k, d->N[0], d->N[1])] += (re*re + im*im)*wf->charge;
  }
      }
    }
  }
}

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	initial_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,
		int	pixel_strategy
	)

Definition at line 561 of file lp_wavefunction.c.

References c, e, get_max_density_value(), density::N, WAVEFUNC_GERROR_DOMAIN, write_density_sweep_slice(), X, Y, and Z.

Referenced by get_potential_for_density_inhomog_multigrid(), get_potential_for_density_inhomog_simple(), write_slices_n_sweeps(), and write_sweeps().

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

Here is the call graph for this function:

Here is the caller graph for this function:

void write_density_sweep_slice	(	enum coord	c,
		const int *	index,
		const struct density *	d,
		char **	filename,
		GError **	err,
		GFile *	location,
		double	maxval,
		int	pixel_strategy
	)

Note:: that the indices in the dotcode go from bottom to top, while the columns in the ppm file go from top to bottom; we need to go backwards here in order for the image to correspond to what the user is expecting to see

Definition at line 412 of file lp_wavefunction.c.

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

Referenced by write_density_sweep().

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

Here is the call graph for this function:

Here is the caller graph for this function:

lp_wavefunction.c File Reference

Classes

Defines

Functions

Define Documentation

Function Documentation