wavefunction_simple_operations.c

Go to the documentation of this file.

/**Performs very simple operations on dotcode wavefunctions.

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

*/

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

#include <math.h>
#include <libdotcode/wavefunction_simple_operations.h>

double dotcode_wavefunction_get_percentage_in_band(const struct wavefunction* wf, int band) {
  struct density* d = dotcode_wavefunction_get_density_in_band(wf, band);
  if(d == NULL) return nan("NaN");
  /*d contains an implicit dx*dy*dz; remove this.*/
  //density_mul_by_const(d, 1.0/(wf->dx[0]*wf->dx[1]*wf->dx[1]));
  //double val = density_integrate_spline(NULL, NULL, NULL, 1e-16, 1e-10, d, NULL, NULL);
  double val = density_integrate_sum_parallel(d, NULL, NULL);
  free_density(d);
  return val;
}

void dotcode_wavefunction_integrate_sum(const struct wavefunction* b, void (*f)(double x, double y, double z, double* re, double* im, void* privdat), const struct wavefunction* a, double* value_re, double* value_im, void* privdat) {
  double f_re, f_im, re, im;
  int x, y, z, band;
  double sum_re = 0;
  double sum_im = 0;
  for(z=0; z<a->N[2]; ++z) {
    for(y=0; y<a->N[1]; ++y) {
      for(x=0; x<a->N[0]; ++x) {
  for(band=0; band<a->bands; ++band) {
    if(f != NULL) {
      f((x-a->N[0]/2)*a->dx[0], (y-a->N[1]/2)*a->dx[1], (z-a->N[2]/2)*a->dx[2], &f_re, &f_im, privdat);
      re = f_re * *get_wavefunction_storage_c_(x, y, z, 0, band, a) 
        + f_im * *get_wavefunction_storage_c_(x, y, z, 1, band, a);
      im = f_re * *get_wavefunction_storage_c_(x, y, z, 1, band, a) 
        - f_im * *get_wavefunction_storage_c_(x, y, z, 0, band, a);
    }else{
      re = *get_wavefunction_storage_c_(x, y, z, 0, band, a);
      im = *get_wavefunction_storage_c_(x, y, z, 1, band, a);
    }
    sum_re += *get_wavefunction_storage_c_(x, y, z, 0, band, b) * re
      + *get_wavefunction_storage_c_(x, y, z, 1, band, b) * im;
    sum_im += *get_wavefunction_storage_c_(x, y, z, 0, band, b) * im
      - *get_wavefunction_storage_c_(x, y, z, 1, band, b) * re;
  }
      }
    }
  }
  *value_re = sum_re;
  *value_im = sum_im;
}
void dotcode_wavefunction_integrate_sum_1(const struct wavefunction* b, void (*f)(double x, double y, double z, double* re, double* im, void* privdat), const struct wavefunction* a, double* value, void* privdat) {
  dotcode_wavefunction_integrate_sum(b, f, a, &value[0], &value[1], privdat);
}

void dotcode_wavefunction_integrate_sum_parallel(const struct wavefunction* b, void (*f)(double x, double y, double z, double* re, double* im, void* privdat), const struct wavefunction* a, double* value_re, double* value_im, void** privdat) {
  double f_re, f_im, re, im;
  int x, y, z, band;
  double sum_re = 0;
  double sum_im = 0;
  #ifdef OPENMP
  #pragma omp parallel for private(x, y, z, band, f_re, f_im, re, im) reduction(+:sum_re,sum_im)
  #endif
  for(z=0; z<a->N[2]; ++z) {
    for(y=0; y<a->N[1]; ++y) {
      for(x=0; x<a->N[0]; ++x) {
  for(band=0; band<a->bands; ++band) {
    if(f != NULL) {
      if(privdat != NULL) {
      f((x-a->N[0]/2)*a->dx[0], (y-a->N[1]/2)*a->dx[1], (z-a->N[2]/2)*a->dx[2], &f_re, &f_im, privdat[
        #ifdef OPENMP
        omp_get_thread_num()
              #else
              0
              #endif
        ]);
      }else{
        f((x-a->N[0]/2)*a->dx[0], (y-a->N[1]/2)*a->dx[1], (z-a->N[2]/2)*a->dx[2], &f_re, &f_im, NULL);
      }
      re = f_re * *get_wavefunction_storage_c_(x, y, z, 0, band, a) 
        + f_im * *get_wavefunction_storage_c_(x, y, z, 1, band, a);
      im = f_re * *get_wavefunction_storage_c_(x, y, z, 1, band, a) 
        - f_im * *get_wavefunction_storage_c_(x, y, z, 0, band, a);
    }else{
      re = *get_wavefunction_storage_c_(x, y, z, 0, band, a);
      im = *get_wavefunction_storage_c_(x, y, z, 1, band, a); 
    }
    sum_re += *get_wavefunction_storage_c_(x, y, z, 0, band, b) * re
      + *get_wavefunction_storage_c_(x, y, z, 1, band, b) * im;
    sum_im += *get_wavefunction_storage_c_(x, y, z, 0, band, b) * im
      - *get_wavefunction_storage_c_(x, y, z, 1, band, b) * re;
  }
      }
    }
  }
  *value_re = sum_re;
  *value_im = sum_im;
}
void dotcode_wavefunction_integrate_sum_parallel_1(const struct wavefunction* b, void (*f)(double x, double y, double z, double* re, double* im, void* privdat), const struct wavefunction* a, double* value, void** privdat) {
  dotcode_wavefunction_integrate_sum_parallel(b, f, a, &value[0], &value[1], privdat);
}