dc_sweep_n_slice.c

Go to the documentation of this file.

/**Gets sweeps and slices of a dotcode wavefucntion.


    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/>.

*/

#include <stdlib.h>
#include <stdio.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <libpostproc/unit_conversion.h>
#include <libpostproc/lp_wavefunction.h>
#include <libdotcode/dotcode_wavefunction.h>
#include "dc_sweep_n_slice_args.h"

void write_slices_n_sweeps(const struct wavefunction* wf, const struct density* d, int* planes, int* sweeps, int* axes, int* gridsite, GFile* slices_base, GFile* sweeps_base, GFile* axes_base);

int main(int argc, char* argv[]) {
  g_type_init();
  struct args args;
  GError *e = NULL;
  get_args(argc, argv, &args, &e);
  if(e != NULL) {
    if(e->code != 1) {
      fprintf(stderr, "ERROR processing arguments: %s\n", e->message);
      g_error_free(e);
      return 1;
    }
  }
  load_wavefunctions(args.wavefunctions, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR loading a requested wavefunction: %s\n", e->message);
    return 2;
  }
  for(GList *l = args.wavefunctions; l!=NULL; l=l->next) {
    struct wavefunction *wf = (struct wavefunction*)(l->data);
    int gridsite[3];
    struct density* d = NULL;
    if(args.use_max_probability_point) {
      GList *wfl = g_list_append(NULL, wf);
      wf->charge=1.0;
      d = get_density_for_list_of_wavefunctions(wfl, &e);
      if(e != NULL) {
  gchar* s =  g_file_get_uri(wf->file);
  fprintf(stderr, "ERROR getting density for wavefunction (%s): %s\n", s, e->message);
  g_free(s);
  return 5;
      }
      g_list_free(wfl);
      double max = density_get_max_and_loc(d, 1, gridsite);
      printf("maxium probability density (*dx^3) is %g\n", max);
    }else{
      for(int i=0; i<3; ++i) {
  if(args.use_gridsite) {
    gridsite[i] = args.gridsite[i];
  }else{
    gridsite[i] = wf->N[i]/2 + args.point[i]/wf->dx[i];
  }
      }
    }
    printf("gridsite is set to <%d, %d, %d> (<%gnm, %gnm, %gnm>)\n", gridsite[0], gridsite[1], gridsite[2], bohr_to_nm((gridsite[0]-wf->N[0]/2)*wf->dx[0]), bohr_to_nm((gridsite[1]-wf->N[1]/2)*wf->dx[1]), bohr_to_nm((gridsite[2]-wf->N[2]/2)*wf->dx[2]));
    if(args.bands == NULL) {
      char* filename = g_file_get_uri(wf->file);
      GString *slices_basedir = g_string_new(filename);
      slices_basedir = g_string_append(slices_basedir, "_dcsns/slices");
      GString *axes_basedir = g_string_new(filename);
      axes_basedir = g_string_append(axes_basedir, "_dcsns/axes");
      GString *sweeps_basedir = g_string_new(filename);
      sweeps_basedir = g_string_append(sweeps_basedir, "_dcsns/sweeps");
      g_free(filename);
      GFile* slices_base = g_file_new_for_uri(slices_basedir->str);
      GFile* axes_base = g_file_new_for_uri(axes_basedir->str);
      GFile* sweeps_base = g_file_new_for_uri(sweeps_basedir->str);
      g_string_free(slices_basedir, TRUE);
      g_string_free(axes_basedir, TRUE);
      g_string_free(sweeps_basedir, TRUE);
      GList *wfl = g_list_append(NULL, wf);
      wf->charge=1.0;
      /*If we fetched the density to get the point of maximum density, we'll reuse it.*/
      if(d == NULL) {
  d = get_density_for_list_of_wavefunctions(wfl, &e);
  if(e != NULL) {
    gchar* s =  g_file_get_uri(wf->file);
    fprintf(stderr, "ERROR getting density for wavefunction (%s): %s\n", s, e->message);
    g_free(s);
    return 5;
  }
      }
      /*Convert the coordinates to nm*/
      for(int i=0; i<3; ++i) d->dx[i]=bohr_to_nm(d->dx[i]);
      write_slices_n_sweeps(wf, d, args.planes, args.sweeps, args.axes, gridsite, slices_base, sweeps_base, axes_base);
      free_density(d);
      g_object_unref(slices_base);
      g_object_unref(axes_base);
      g_object_unref(sweeps_base);
    }else{
      /*If we fetched the density to get the point of maximum density, free it.*/
      if(d != NULL) free_density(d);
      for(GList *l = args.bands; l!=NULL; l=l->next) {
  d = dotcode_wavefunction_get_density_in_band(wf, *((int*)(l->data)));
  char* filename = g_file_get_uri(wf->file);
  GString *slices_basedir = g_string_new(filename);
  g_string_append_printf(slices_basedir, "_dcsns/band_%s_slices", dotcode_band_short_names[*(int*)(l->data)]);
  GString *axes_basedir = g_string_new(filename);
  g_string_append_printf(axes_basedir, "_dcsns/band_%s_axes", dotcode_band_short_names[*(int*)(l->data)]);
  GString *sweeps_basedir = g_string_new(filename);
  g_string_append_printf(sweeps_basedir, "_dcsns/band_%s_sweeps", dotcode_band_short_names[*(int*)(l->data)]);
  g_free(filename);
  GFile* slices_base = g_file_new_for_uri(slices_basedir->str);
  GFile* axes_base = g_file_new_for_uri(axes_basedir->str);
  GFile* sweeps_base = g_file_new_for_uri(sweeps_basedir->str);
  /*Convert the coordinates to nm*/
  for(int i=0; i<3; ++i) d->dx[i]=bohr_to_nm(d->dx[i]);
  write_slices_n_sweeps(wf, d, args.planes, args.sweeps, args.axes, gridsite, slices_base, sweeps_base, axes_base);
  free_density(d);
  g_string_free(slices_basedir, TRUE);
  g_string_free(axes_basedir, TRUE);
  g_string_free(sweeps_basedir, TRUE);
  g_object_unref(slices_base);
  g_object_unref(axes_base);
  g_object_unref(sweeps_base);
      }
    }
  }
  return 0;
}


void write_slices_n_sweeps(const struct wavefunction* wf, const struct density* d, int* planes, int* sweeps, int* axes, int* gridsite, GFile* slices_base, GFile* sweeps_base, GFile* axes_base) {
  GError *e = NULL;
  for(int i=0; i<3; ++i) {
    char* filename;
    if(sweeps[i]) {
      g_file_make_directory_with_parents(sweeps_base, NULL, &e);
      if(e != NULL) {
  if(e->code != G_IO_ERROR_EXISTS) {
    gchar* s = g_file_get_uri(wf->file);
    fprintf(stderr, "ERROR making directory for sweep along %d(%c) (directory %s): %s", i, (i==0)?'X':(i==1)?'Y':(i==2)?'Z':'?', s, e->message);
    g_free(s);
    exit(6);
  }
  g_error_free(e);
  e=NULL;
      }
      write_density_sweep(i, d, &filename, &e, sweeps_base, FALSE);
      if(e != NULL) {
  gchar* s = g_file_get_uri(wf->file);
  fprintf(stderr, "ERROR writing sweep along axis %d(%c) for wavefunction %s: %s", i, (i==0)?'X':(i==1)?'Y':(i==2)?'Z':'?', s, e->message);
  g_free(s);
  exit(6);
      }
      g_free(filename);
    }
    if(axes[i]) {
      g_file_make_directory_with_parents(axes_base, NULL, &e);
      if(e != NULL) {
  if(e->code != G_IO_ERROR_EXISTS) {
    gchar* s = g_file_get_uri(wf->file);
    fprintf(stderr, "ERROR making directory for axis along %d(%c) (directory %s): %s", i, (i==0)?'X':(i==1)?'Y':(i==2)?'Z':'?', s, e->message);
    g_free(s);
    exit(6);
  }
  g_error_free(e);
  e=NULL;
      }
      write_density_line(i, gridsite, d, &filename, &e, axes_base, '\t', "nm", "");
      if(e != NULL) {
  gchar* s = g_file_get_uri(wf->file);
  fprintf(stderr, "ERROR writing line along axis %d(%c) for wavefunction %s: %s", i, (i==0)?'X':(i==1)?'Y':(i==2)?'Z':'?', s, e->message);
  g_free(s);
  exit(7);
      }
      g_free(filename);
    }
    if(planes[i]) {
      int coord=0;
      /*XY->Z coordinate is to be taken*/
      if(planes[0]) coord=2;
      /*XZ->Y*/
      else if(planes[1]) coord=1;
      /*YZ->X*/
      else coord=0;
      g_file_make_directory_with_parents(slices_base, NULL, &e);
      if(e != NULL) {
  if(e->code != G_IO_ERROR_EXISTS) {
    gchar* s = g_file_get_uri(wf->file);
    fprintf(stderr, "ERROR making directory for %d(%s) plane (directory %s): %s", i, (i==0)?"XY":(i==1)?"XZ":(i==2)?"YZ":"?", s, e->message);
    g_free(s);
    exit(6);
  }
  g_error_free(e);
  e=NULL;
      }
      write_density_slice(coord, &gridsite[coord], d, &filename, &e, slices_base, '\t', "nm", "");
      if(e != NULL) {
  gchar* s = g_file_get_uri(wf->file);
  fprintf(stderr, "ERROR writing %d(%s) plane: %s", i, (i==0)?"XY":(i==1)?"XZ":(i==2)?"YZ":"?", e->message);
  g_free(s);
  exit(8);
      }
      g_free(filename);
    }
  }
}