lp_wavefunction.h

Go to the documentation of this file.

/*
** wavefunction.h
** 

    Copyright (C) 2009 Joseph Pingenot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU Affero General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

** Made by Johnny Q. Hacker
** Login   <solarion@nathan>
** 
** Started on  Tue Nov 17 22:28:05 2009 Johnny Q. Hacker
** Last update Tue Nov 17 22:28:05 2009 Johnny Q. Hacker
*/

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

#ifndef     WAVEFUNCTION_H_
# define    WAVEFUNCTION_H_

#ifdef __cplusplus
extern "C" {
#endif


/*For clarity, X, Y, and Z specify respective array indices.*/
enum coord {X, Y, Z};

struct wavefunction {
  GFile* file;
  double charge;
  int N[3];
  double dx[3];
  int bands;
  struct wavefunc_header header;
  double* storage;
};

struct density {
  int N[3];
  double dx[3];
  double* storage;
};


/**Gets a new density with a double-resolution grid, using an average for interpolation
 *\param original original density
 *\param auto_scale flags on which dimension to use to scale: 1->Y 2->Y 4->Z (combine for multiple directions)
 *\return pointer to a density, or NULL if memory allocation failed.
 *\note the new grid doesn't add an outside row of new sites; it only places a new grid site between two existing ones
 */
  struct density* get_double_resolution_density(const struct density* original, int auto_scale);

/**Gets a new density with a half-resolution grid, using an average for assignment
 *\param original original density
 *\param e pointer to pointer to Gerror in which to store the error information
 *\param auto_scale flags on which dimension to use to scale: 1->Y 2->Y 4->Z (combine for multiple directions)
 *\return pointer to a density, or NULL if memory allocation failed.
 *\note the new grid doesn't add an outside row of new sites; it only places a new grid site between two existing ones
 */
  struct density* get_half_resolution_density(const struct density* original, GError **e, int auto_scale);

/**Gets a density from a list of *loaded* wavefunctions
 *\param list the list of *loaded* wavefunction structures.
 *\param 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!
 *\note The charge is used! If you just want pre-normalized |psi|^2,
 *set charge to 1.
 */
struct density* get_density_for_list_of_wavefunctions(GList *list, GError **err);
/**Loads the wavefunctions from the respective GFiles.
 *\param wf_list List of wavefunctions
 *\param err pointer to a pointer in which to store the GError (if
 *there's an error)
*/
void load_wavefunctions(GList* wf_list, GError** err);
void init_density(struct density* d, double initial_value);
void add_to_density(struct density* d, struct wavefunction* wf);
/**Puts a density into *one section* of a metadensity.
 *\param md the target meta_density
 *\param d the source density
 *\param Nn an array of ints representing the number of neighbors in
 *each direction.
 *\param x location of this block to be inserted.
 *\param y location of this block to be inserted.
 *\param z location of this block to be inserted.
 *\param invert_x if true, the x-coordinates are read backwards from
 *the source.
 *\param invert_y if true, the y-coordinates are read backwards from
 *the source.
 *\param invert_z if true, the z-coordinates are read backwards from
 *the source.
 *\param offsets offset in grid sites between neighborsa
 */
  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);
/**Writes out a single, ascii density slice or the entire box.
 *\param c the coordinate to output
 *\param index pointer to an int containing the index which should be
 *written out
 *\param d the density to be written
 *\param filename pointer to a pointer to a char, to provide the
 *filename of the written file.
 *\param err pointer to a pointer to a GError, for returning any error
 *info
 *\param location base directory to output, or NULL if cwd. If index
 *is NULL, this is the file to write the dump into.
 *\param sepchar character to use to separate the columns in the
 *output
 *\return nothing
 */
void write_density_slice(enum coord c, const int* index, const struct density* d, char** filename, GError **err, GFile *location, char sepchar, const char* coord_units, const char* density_units);
void write_density_sweep(enum coord c, const struct density* d, char** filename, GError **err, GFile* location, gboolean sweep_log);
void write_density_line(enum coord axis, const int* gridsite, const struct density* d, char** filename, GError **err, GFile *location, char sepchar, const char* coord_units, const char* density_units);

void density_mul_by_const(struct density* d, double c);
/**Gets a block of the meta_density back out.
 *\param md the meta denstiy
 *\param d an initialized density the sizeo f one block
 *\param buffzone the buffer zone size around the meta_density
 *\param Mn the number of neighbors in each direction
 *\param x the x-coordinate of the block
 *\param y the y-coordinate of the block
 *\param z the z-coordinate of the block
 *\param offsets offset in grid sites between neighborsa
 */
  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);

/**Writes the density into a povray visualization.
 *\param d the density to be written out
 *\param towrite the file into which the povray visualization should be written
 *\paarm err pointer to a pointer to a GError, used to "return" a GError pointer if an error ocurred.
 */
void write_density_povray(const struct density* d, GFile *towrite, int pixel_strategy, GError **err);

/**Expands a density box using material mat
 *\param od the original density
 *\param I a 3-int array indicating how many grid sites to expand outward.
 *\param value a double containing the value of the density throughout the new region
 *\returns a pointer to a new density.
 */
struct density* expand_density_box_with_value(const struct density* od, const int* I, double value, GError **err);
/**Expands a density box using the value from the outermost grid site.
 *\param od the original density
 *\param I a 3-int array containing the number of grid sites to expand
 *outwards
 */
struct density* expand_density_box(const struct density* od, const int* I, GError **err);

struct density* new_density(const int* N, const double* dx, double intial_value);
struct density* new_density_copy(const struct density* orig);
void free_density(struct density* d);
  int op_density_to_density(struct density* to, enum density_op op, const struct density* from);
  /**Executes the specified callback
   *\param d density struct
   *\param callback pointer to function to call
   *\param privdat private data for the callback
   *\return nothing
   */
  void density_execute_callback(const struct density* d, void (*callback)(int x, int y, int z, const struct density*d, void* privdat), void* privdat);
  /**same, but takes two density functions*/
  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);
  /**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.
   *\param d density struct
   *\param callback pointer to function to call
   *\param privdat array of private data pointers for the callback, one per thread.
   *\param x_report array of ints, one per thread.  Set to 1 if a thread is executed, and is left untouched otherwise.
   *\return nothing
   */
  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);


  /**returns the maximum percent (at the point) difference between two densities, or NaN if the dimensions differ
   *\param d1 density to compare
   *\param d2 density to compare
   *\param 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.
   *\param scaling value to use for scaling if scaleby==3, else ignored
   */
  double density_get_max_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
   *\param d1 density to compare
   *\param d2 density to compare
   *\param 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.
   *\param scaling value to use for scaling if scaleby==3, else ignored
   */
  double density_get_min_abs_delta_pct(const struct density* d1, const struct density* d2, int scaleby, double scaling);
  /**Gets the maximum value in the density
   *\param d density to inspect
   *\param abs if true, compare absolute value instead of raw value.
   *\return the maximum value in d
   */
  double density_get_max(const struct density* d, int abs);
  double density_get_max_and_loc(const struct density* d, int abs, int* location);
  /**Gets the minimum value in the density
   *\param d density to inspect
   *\param abs if true, compare absolute value instead of raw value.
   *\return the minimum value in d
   */
  double density_get_min(const struct density* d, int abs);
  double density_get_min_and_loc(const struct density* d, int abs, int* location);

struct dotcode_grid* dotcode_grid_from_density(const struct density *d, GError **err);
long unsigned int get_density_storage_index_(unsigned int x, unsigned int y, unsigned int z, unsigned int Nx, unsigned int Ny);

  /**gets pointer to the cell at this location
   *\param x  coordinate
   *\param y  coordinate
   *\param z  coordinate
   *\param d  density structure
   *\return pointer to the density value at the specified point
   *\note no error checking is performed
   */
  double* get_density_storage_at_(unsigned x, unsigned y, unsigned z,struct density* d);
  /**gets const pointer to the cell at this location
   *\param x  coordinate
   *\param y  coordinate
   *\param z  coordinate
   *\param d  density structure
   *\return const pointer to the density value at the specified point
   *\note no error checking is performed
   */
  const double* get_density_storage_at_c_(unsigned x, unsigned y, unsigned z, const struct density* d);

  /**Evaluates an integral of the density using GSL and a spline fit.
   *\param z_error_estimate: pointer to a double in which to store the estimated integration error from the final integration along z, or NULL if the error isn't to be saved.
   *\param y_error_estimate: pointer to an array (of size wf->N[2], one for each z point) of doubles in which to store the estimated integration error for each y integral, or NULL if the error isn't to be saved.
   *\param x_error_estimate: pointer to an array (of size wf->N[1]*wf->N[2], one for each (y, z) pairing) of doubles in which to store the estimated error for each x integral, or NULL if the error isn't to be saved.
   *\param abserr the requested absolute error (see GSL documentation)
   *\param relerr the requested relative error (see GSL documentation)
   *\param d the density over which to integrate
   *\param f a function to call at each point (if NULL, does nothing).
   *\param privdat private data to pass to f.
   *\return the integral
   *\note d will be fit to a cubic spline before integration.
   *\note THIS IS BROKEN
   */
  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);

  /**Evaluates an integral of the density using GSL and a spline fit, parallelized with OpenMP on the z coordinate
   *\param d the density over which to integrate
   *\param f a function to call at each point (if NULL, does nothing).
   *\param privdat array of private data to pass to the callback (one for each thread, indexed by omp_get_thread_num()).
   *\return the integral
   *\note d will be fit to a cubic spline before integration.
   *\note the same f can be used for both the single- and multi-threaded functions.
   *\note if privdat is NULL, NULL will be passed to f (this saves you having to allocate an array of NULLs
   *\note THIS IS SUPER-BROKEN
   */
  double density_integrate_spline_parallel(const struct density* d, double (*f)(double x, double y, double z, void* privdat), 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.
   *\param d the density over which to integrate
   *\param f a function to call at each point (if NULL, does nothing).
   *\param privdat array of private data to pass to the callback (one for each thread, indexed by omp_get_thread_num()).
   *\return the integral
   *\note d will be fit to a cubic spline before integration.
   *\note the same f can be used for both the single- and multi-threaded functions.
   *\note if privdat is NULL, NULL will be passed to f (this saves you having to allocate an array of NULLs
   */
  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);


#ifdef __cplusplus
}
#endif

#endif      /* !WAVEFUNCTION_H_ */