func_sum.c

Go to the documentation of this file.

/*
 *func_sum
 *Outputs the sum of a wavefunction's probability.
 
    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/>.

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <glib-2.0/glib.h>
#include <libdotcode/function_parser.h>
#include "config.h"
#include "common.h"

#define EV_PER_HARTREE 27.2113834
#define NM_PER_BOHR 0.05291772083

/*This many gridsites in*/
#define INSET 10

int main(int argv, char* argc[]) {
  g_type_init();
  int i;
  struct wavefunc_header header;
  double *matrix;
  int err;
  int retval;
  /*probability of being within the specified region*/
  double prob;
  int x, y, z, band;
  double re, im;
  /*TOTAL sum*/
  double tsum;
  if(argv == 1) {
    fprintf(stderr, "Umm, you probably want to pass me a file or space-separated list of files to read.\n");
    exit(1);
  }
  for(i=1; i<argv; i++) {
    retval = read_wavefunction(argc[i], &header, &matrix, &err);
    if(retval != 0) {
      fprintf(stderr, "ERROR: read_wavefunction returned %d, err=%d (%s).\n", retval, err, strerror(err));
    }
    tsum=prob=0;
    for(z=0; z<header.Nz; ++z) {
      for(y=0; y<header.Ny; ++y) {
  for(x=0; x<header.Nx; ++x) {
    for(band=0; band<header.Nc; ++band) {
      re = matrix[band*2 + x*header.Nc*2 + y*header.Nx*header.Nc*2 + z*header.Ny*header.Nx*header.Nc*2];
      im = matrix[1 + band*2 + x*header.Nc*2 + y*header.Nx*header.Nc*2 + z*header.Ny*header.Nx*header.Nc*2];
      #ifdef DEBUG
      fprintf(stderr, "x=%d,y=%d,z=%d,b=%d: %g + i%g\n", x, y, z, band, re, im);
      #endif
      tsum += re*re + (im*im);
      if(
         (x < INSET) || ((header.Nx-x-1) < INSET)
         || (y < INSET) || ((header.Ny-y-1) < INSET)
         || (z < INSET) || ((header.Nz-z-1) < INSET)
         ){ 
        prob += re*re + (im*im);
      }
    }
  }
      }
    }
    printf("%s\t%g\n", argc[i], prob);
  }
  exit(0);
}