dc_get_position.c

Go to the documentation of this file.

/**Gets the composition of a dotcode wavefunction
 

    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 <stdio.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <math.h>
#include <libpostproc/lp_wavefunction.h>
#include <libdotcode/dotcode_wavefunction.h>
#include <libdotcode/wavefunction_simple_operations.h>
#include <libdotcode/dotcode_data.h>
#include <libpostproc/unit_conversion.h>

static void get_x(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = x;
  *im = 0;
}
static void get_y(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = y;
  *im = 0;
}
static void get_z(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = z;
  *im = 0;
}

static void get_sigma_x(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = (x - *(double*)privdat)*(x - *(double*)privdat);
  *im = 0;
}
static void get_sigma_y(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = (y - *(double*)privdat)*(y - *(double*)privdat);
  *im = 0;
}
static void get_sigma_z(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = (z - *(double*)privdat)*(z - *(double*)privdat);
  *im = 0;
}

static void get_x_110(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = 1/sqrt(2)*(x+y);
  *im = 0;
}
static void get_y_110(double x, double y, double z, double* re, double* im, void* privdat) {
  *re = 1/sqrt(2)*(y-x);
  *im = 0;
}

static void get_sigma_x_110(double x, double y, double z, double* re, double* im, void* privdat) {
  double x_110 = 1/sqrt(2)*(x+y);
  *re = (x_110 - *(double*)privdat)*(x_110 - *(double*)privdat);
  *im = 0;
}
static void get_sigma_y_110(double x, double y, double z, double* re, double* im, void* privdat) {
  double y_110 = 1/sqrt(2)*(y-x);
  *re = (y_110 - *(double*)privdat)*(y_110 - *(double*)privdat);
  *im = 0;
}

int main(int argc, char** argv) {
  g_type_init();
  gboolean basis_110 = TRUE;
  printf("#file\t<x>(nm)\t<y>(nm)\t<z>(nm)\tsigma_x(nm)\tsigma_y(nm)\tsigma_z(nm)\n");
  for(int wfn=1; wfn<argc; ++wfn) {
    GList* l=NULL;
    GFile* f = g_file_new_for_commandline_arg(argv[wfn]);
    struct wavefunction wf = {.file=f, .charge=1, .N={0, 0, 0}, .dx={0, 0, 0}, .bands=0, .storage=NULL};
    l = g_list_append(l, &wf);
    GError *e = NULL;
    load_wavefunctions(l, &e);
    if(e == NULL) {
      printf("%s", argv[wfn]);
      double im = 0;
      double avg_x=0;
      double avg_y=0;
      double avg_z=0;
      double sigma_x=0;
      double sigma_y=0;
      double sigma_z=0;
      dotcode_wavefunction_integrate_sum_parallel(&wf, basis_110?get_x_110:get_x, &wf, &avg_x, &im, NULL);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting <x>: imaginary part was non-zero!\n");
  return 1;
      }
      dotcode_wavefunction_integrate_sum_parallel(&wf, basis_110?get_y_110:get_y, &wf, &avg_y, &im, NULL);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting <y>: imaginary part was non-zero!\n");
  return 1;
      }
      dotcode_wavefunction_integrate_sum_parallel(&wf, get_z, &wf, &avg_z, &im, NULL);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting <z>: imaginary part was non-zero!\n");
  return 1;
      }
      #ifdef OPENMP
      int N_threads=omp_get_max_threads();
      #else
      int N_threads=1;
      #endif
      void* callback_data[N_threads];
      /*Since the callback data is only read, we can have one copy for all threads*/
      for(int i=0; i<N_threads; ++i) callback_data[i] = &avg_x;
      dotcode_wavefunction_integrate_sum_parallel(&wf, basis_110?get_sigma_x_110:get_sigma_x, &wf, &sigma_x, &im, callback_data);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting sigma_x: imaginary part was non-zero!\n");
  return 1;
      }
      for(int i=0; i<N_threads; ++i) callback_data[i] = &avg_y;
      dotcode_wavefunction_integrate_sum_parallel(&wf, basis_110?get_sigma_y_110:get_sigma_y, &wf, &sigma_y, &im, callback_data);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting sigma_y: imaginary part was non-zero!\n");
  return 1;
      }
      for(int i=0; i<N_threads; ++i) callback_data[i] = &avg_z;
      dotcode_wavefunction_integrate_sum_parallel(&wf, get_sigma_z, &wf, &sigma_z, &im, callback_data);
      if(fabs(im) >= 1e-8) {
  fprintf(stderr, "ERROR getting sigma_z: imaginary part was non-zero!\n");
  return 1;
      }
      sigma_x = sqrt(sigma_x);
      sigma_y = sqrt(sigma_y);
      sigma_z = sqrt(sigma_z);
      printf("\t%g\t%g\t%g\t%g\t%g\t%g\n", bohr_to_nm(avg_x), bohr_to_nm(avg_y), bohr_to_nm(avg_z), bohr_to_nm(sigma_x), bohr_to_nm(sigma_y), bohr_to_nm(sigma_z));

      dotcode_wavefunction_free_resources(&wf);
    }else{
      fprintf(stderr, "ERROR loading wavefunction (%s): %s", argv[wfn], e->message);
    }
    g_list_free(l);
    if(e != NULL) g_error_free(e);
  }
  return 0;
}