wavefunc_max_change.c

Go to the documentation of this file.

/*
** wavefunc_max_change.c
** 
** Made by (Johnny Q. Hacker)
** Login   <solarion@borkborkbork>
** 
    Copyright (C) 2010 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/>.

** Started on  Wed Mar  3 16:24:09 2010 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdio.h>
#include <math.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include <libpostproc/lp_wavefunction.h>
#include <libpostproc/unit_conversion.h>

int main(int argc, char** argv) {
  g_type_init();
  if(argc != 3) {
    fprintf(stderr, "ERROR: only got %d arguments, but must have two: the names of two wavefunctions to compare!\n", argc);
    fprintf(stdout, "USAGE: %s a b\n", argv[0]);
    fprintf(stdout, "Returns the maximum difference in wave function probability densities between a and b as an absolute number (|b|^2 - |a|^2), and as a percentage of the average of a.\n");
    return 1;
  }
  
  GError *e = NULL;

  char* wffa=argv[1];
  char* wffb=argv[2];
  
  struct wavefunction a, b;
  a.file = g_file_new_for_commandline_arg(wffa);
  b.file = g_file_new_for_commandline_arg(wffb);
  b.charge = a.charge = +1.0;

  /*This is kind of cumbersome, but it works right now. Perhaps I'll write a
    one-wf version at some point.*/   
  GList* la = g_list_append(NULL, &a);
  GList* lb = g_list_append(NULL, &b);

  load_wavefunctions(la, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR: failed to load first wavefunction (%s): %s\n", wffa, e->message);
    return 4;
  }
  load_wavefunctions(lb, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR: failed to load first wavefunction (%s): %s\n", wffa, e->message);
    return 5;
  }

  /*Verify dimensions.*/
  if(a.N[0] != b.N[0]) {
    fprintf(stderr, "ERROR: wavefunction x dimensions are not identical! a.N[0]=%d b.N[0]=%d\n", a.N[0], b.N[0]);
    return 6;
  }
  if(a.N[1] != b.N[1]) {
    fprintf(stderr, "ERROR: wavefunction y dimensions are not identical! a.N[0]=%d b.N[0]=%d\n", a.N[1], b.N[1]);
    return 7;
  }
  if(a.N[2] != b.N[2]) {
    fprintf(stderr, "ERROR: wavefunction z dimensions are not identical! a.N[0]=%d b.N[0]=%d\n", a.N[2], b.N[2]);
    return 8;
  }

  struct density *density_a = get_density_for_list_of_wavefunctions(la, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR: failed to process first wavefunction (%s): %s\n", wffa, e->message);
    return 2;
  }
  struct density* density_b = get_density_for_list_of_wavefunctions(lb, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR: failed to process second wavefunction (%s): %s\n", wffb, e->message);
    return 3;
  }

  double delta;
  double max_delta = nan("Nan");
  double avg_a = 0;
  double avg_b = 0;
  double avg_delta = 0;
  /*Values of a, b at the max delta point.*/
  double bmax_del=0, amax_del=0;
  double aye=0, bee=0;
  int xmax=0, ymax=0, zmax=0;
  for(int k=0; k<density_a->N[2]; k++) {
    for(int j=0; j<density_a->N[1]; j++) {
      for(int i=0; i<density_a->N[0]; i++) {
  bee = density_b->storage[get_density_storage_index_(i, j, k, density_a->N[0], density_a->N[1])];
  aye = density_a->storage[get_density_storage_index_(i, j, k, density_a->N[0], density_a->N[1])];
  delta= bee - aye;
  //fprintf(stderr, "%d/%d/%d: %g\t%g\t%g\t%g\n", i, j, k, aye, bee, delta, max_delta);
  if(isnan(max_delta) || (fabs(delta) > fabs(max_delta))) {
    max_delta = delta;
    bmax_del = bee;
    amax_del = aye;
    xmax=i;
    ymax=j;
    zmax=k;
  }
  avg_delta += fabs(delta);
  avg_a += density_a->storage[get_density_storage_index_(i, j, k, density_a->N[0], density_a->N[1])];
  avg_b += density_b->storage[get_density_storage_index_(i, j, k, density_b->N[0], density_b->N[1])];
      }
    }
  }
  avg_a /= (density_a->N[0]*density_a->N[1]*density_a->N[2]);
  avg_b /= (density_b->N[0]*density_b->N[1]*density_b->N[2]);
  avg_delta /= (density_a->N[0]*density_a->N[1]*density_a->N[2]);
  printf("#max_delta\tavg_delta\ta@max\tb@max\t%%a@max\t%%b@max\t%% avg_a\t%%avg_b\t%%avg_delta\tix@max\tiy@max\tiz@max\tx@max(nm)\ty@max(nm\tz@max(nm)\n");
  printf("%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%d\t%d\t%d\t%g\t%g\t%g\n", max_delta, avg_delta, amax_del, bmax_del, max_delta/amax_del, max_delta/bmax_del, max_delta/avg_a, max_delta/avg_b, max_delta/avg_delta, xmax, ymax, zmax, bohr_to_nm(xmax*density_a->dx[0]), bohr_to_nm(xmax*density_a->dx[1]), bohr_to_nm(xmax*density_a->dx[2]));
  return 0;
}