getChainPotential.c

Go to the documentation of this file.

/*
** getChainPotential.c
** 


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

** Made by (Johnny Q. Hacker)
** Login   <solarion@nathan>
** 
** Started on  Sat Nov 14 17:15:07 2009 Johnny Q. Hacker
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdlib.h>
#include <stdio.h>
#include <glib-2.0/glib.h>
#include <gio/gio.h>
#include "getChainPotential_args.h"
#include <libpostproc/lp_wavefunction.h>
#include <libpostproc/fftw_wavefunc.h>
#include <libpostproc/unit_conversion.h>
#include <libdotcode/read_binary_data.h>
#include <libdotcode/write_binary_data.h>
#include <libdotcode/matdb-dotcode.h>
#include <libdotcode/grid1.h>

void write_sweeps(struct args* args, struct density* d, char* subdir);
void write_slices(struct args* args, struct density* d, char* subdir, const char* coord_units, const char* density_units);

int main(int argc, char** argv) {
  g_type_init();
  GError *e = NULL;
  struct args *args = process_args(argc, argv, &e);
  if(args == NULL) {
    fprintf(stderr, "Failed to process the args: ");
    if(e != NULL) {
      fprintf(stderr, "%s\n", e->message);
    }else{
      fprintf(stderr, "(got no GError, for some reason. This is a problem in the program one way or the other!)\n");
    }
    return 1;
  }
  /*At this point, the wfs haven't been loaded into the
   *wavefunction struct. Do that now.*/
  e=NULL;
  load_wavefunctions(args->wf_list, &e);
  if(e!=NULL){
    fprintf(stderr, "Got error loading the wavefunctions: %s\n", e->message);
    return 2;
  }
  e=NULL;
  load_wavefunctions(args->center_wf_list, &e);
  if(e!=NULL){
    fprintf(stderr, "Got error loading the center wavefunctions: %s\n", e->message);
    return 2;
  }
  e=NULL;
  /*Sum up all of the neighbors in a block*/
  struct density *d = NULL;
  if(args->wf_list != NULL) {
    d = get_density_for_list_of_wavefunctions(args->wf_list, &e);
    if(e != NULL) {
      fprintf(stderr, "ERROR: error while getting list of wavefunctions: %s\n", e->message);
      return 7;
    }
  }
  struct density *center_d = NULL;
  if(args->center_wf_list != NULL) {
    center_d = get_density_for_list_of_wavefunctions(args->center_wf_list, &e);
    if(e != NULL) {
      fprintf(stderr, "ERROR: error while getting list of center wavefunctions: %s\n", e->message);
      return 8;
    }
  }

  /*For now, just make a rho=epsilon_0*/
  /*epsilon_0 is 8.854 187 817... x 10-12 F m-1
   *Units are in Bohr, and 1 bohr = 0.529 177 208 59 x 10-10 m
   *(source: NIST physical constants database
   */
  struct density* rhoprime;
  /*Read in the density*/
  e=NULL;
  struct wavefunction* inwavefunc = (struct wavefunction*)(args->wf_list->data);
  if(inwavefunc == NULL) {
    inwavefunc = (struct wavefunction*)(args->center_wf_list->data);
    if(inwavefunc == NULL) {
      fprintf(stderr, "ERROR: Unable to find the first wavefunction!\n");
      return 23;
    }
  }
  GFile* grid1_file = g_file_get_child(inwavefunc->file, ".trun.displ/grid");
  struct dotcode_grid1* g1 = dotcode_grid1_get_from_gfile(grid1_file, &e);
  if(e != NULL) {
    char* pathstr = g_file_get_uri(grid1_file);
    fprintf(stderr, "ERROR: unable to open and process the grid (%s): %s\n", pathstr, e->message);
    g_free(pathstr);
    g_error_free(e);
    return 20;
  }
  GFile* mdbfile = g_file_get_child(inwavefunc->file, ".matdb");
  struct matdb* mdb = read_matdb_dotcode_gio(mdbfile, &e);
  if(e != NULL) {
    char* pathstr = g_file_get_uri(mdbfile);
    fprintf(stderr, "ERROR: unable to open and process the matdb (%s): %s\n", pathstr, e->message);
    g_free(pathstr);
    g_error_free(e);
    return 21;
  }
  struct density* epsilon = (struct density*)malloc(sizeof(struct density));
  if(epsilon == NULL) {
    fprintf(stderr, "ERROR: unable to allocate memory for epsilon.\n");
    return 23;
  }
  for(int i=0; i<3; i++) {
    epsilon->N[i] = g1->N[i];
    epsilon->dx[i] = g1->dx[i];
  }
  epsilon->storage = dotcode_grid1_get_property_grid(g1, mdb, "epsilon_0", &e);
  if(e != NULL) {
    char* pathstr1 = g_file_get_uri(grid1_file);
    char* pathstr2 = g_file_get_uri(mdbfile);
    fprintf(stderr, "ERROR: unable to process the grid (%s) and matdb (%s) into  a epsilon_0 grid: %s\n", pathstr1, pathstr2, e->message);
    g_free(pathstr1);
    g_free(pathstr2);
    g_error_free(e);
    return 22;
  }
  /*need to multiply epsilon by epsilon_0*/
  density_mul_by_const(epsilon, 8.854187817e-12);

  struct density* meta_density;
  struct density* meta_epsilon;
  /*this algo almost certainly blows our cache and is slow.*/
  if((args->N_neighbors[0] == 0) && (args->N_neighbors[1] == 0) && (args->N_neighbors[2] == 0)) {
    meta_density = d;
    meta_epsilon = epsilon;
  }else{
    int N[3];
    long unsigned int t = 1;
    int extra_x = 0;
    int extra_y = 0;
    int extra_z = 0;
    /*If a direction has a non-zero number of neighbors, we need to add an extra neighbor in that direction to get an infinite chain with periodic BCs.*/
    if(args->extra_neighbor[0]) extra_x = 1;
    if(args->extra_neighbor[1]) extra_y = 1;
    if(args->extra_neighbor[2]) extra_z = 1;
    for(int i=0; i<3; i++) {
      /*There are 2*N_neighbors[i] blocks (one on each side), plus the center block.*/
      N[i] = (2*(args->N_neighbors[i])+1+args->extra_neighbor[i])*(d->N[i] + args->neighbor_offset[i]);
      t *= N[i];
    }
    meta_density = new_density(N, d->dx, 0.0);
    if(meta_density == NULL) {
      fprintf(stderr, "Error: unable to allocate storage for an %dx%dx%d box (%lu total doubles (meta_density))\n", (2*(args->N_neighbors[0])+1)*(d->N[0]), (2*(args->N_neighbors[1])+1)*(d->N[1]), (2*(args->N_neighbors[2])+1)*(d->N[2]), t);
      return 3;
    }
    //Set the meta_epsilon default value from the 0, 0, 0 corner.
    /**\todo set this from user-specified value*/
    meta_epsilon = new_density(N, d->dx, *get_density_storage_at_c_(0, 0, 0, epsilon));
    if(meta_epsilon == NULL) {
      fprintf(stderr, "Error: unable to allocate storage for an %dx%dx%d box (%lu total doubles (meta epsilon))\n", (2*(args->N_neighbors[0])+1)*(d->N[0]), (2*(args->N_neighbors[1])+1)*(d->N[1]), (2*(args->N_neighbors[2])+1)*(d->N[2]), t);
      return 3;
    }
    int x_counter = 0, y_counter = 0, z_counter=0;
    gboolean invert_x, invert_y, invert_z;
    /*If we have an even number of neighbors, our first coordinate
      is an even neighbor and therefore the even coordinates are even
      neighbors and therefore shouldn't be flipped.  If we have an odd
      number of neighbors, the first coordinate is an odd neighbor,
      and therefore the even coordinates are odd neighbors and should
      be flipped.
      TO CHECK THE LOGIC: the center coordinate is *always* an
      even coordinate!
      NOTE that the counter will flip before it's checked.
      THE ALGORITHM being used here is that 1*-1 = 1; -1*-1 = 1; 0*-1
      = 0;  So we start at 0 if we're not counting (and we'll never be
      flip; we'll flip when the counter is +1, and we'll not flip when
      the counter is -1.
      UPDATE: the previous algo was crap. We need to flip the
      appropriate axes when we're on an "odd" block.
    */
    if((args->N_neighbors[0] % 2) == 0) {
      x_counter = 1;
    }else{
      x_counter = -1;
    }
    if((args->N_neighbors[1] % 2) == 0) {
      y_counter = 1;
    }else{
      y_counter = -1;
    }
    if((args->N_neighbors[2] % 2) == 0) {
      z_counter = 1;
    }else{
      z_counter = -1;
    }
    for(int z=0; z<(2*(args->N_neighbors[2])+1+extra_z); z++) {
      z_counter *= -1;
      for(int y=0; y<(2*(args->N_neighbors[1])+1+extra_y); y++) {
  y_counter *= -1;
  for(int x=0; x<(2*(args->N_neighbors[0])+1+extra_x); x++) {
    x_counter *= -1;
    if((x_counter*y_counter*z_counter) > 0) {
      invert_x = (args->invert_odd)&1;
      invert_y = (args->invert_odd)&2;
      invert_z = (args->invert_odd)&4;
    }else{
      invert_x = 0;
      invert_y = 0;
      invert_z = 0;
    }
    fprintf(stderr, "%d/%d/%d: invert_x: %c invert_y: %c invert_z: %c (x_counter=%d,y_counter=%d,z_counter=%d\n", x, y, z, (invert_x)?'y':'n', (invert_y)?'y':'n', (invert_z)?'y':'n', x_counter, y_counter, z_counter);
    /*Check to see if this is the middle block*/
    #ifdef DEBUG
    fprintf(stderr, "x=%d(s%d) y=%d(s%d) z=%d(s%d):\n", x, (args->N_neighbors[0]), y, (args->N_neighbors[1]), z, (args->N_neighbors[2]));
    #endif
    if((z == (args->N_neighbors[2]))
       && (y == (args->N_neighbors[1]))
       && (x == (args->N_neighbors[0]))) {
            #ifdef DEBUG
        fprintf(stderr, "center!\n");
            #endif
        if(center_d != NULL) {
    put_density_into_metadensity(meta_density, center_d, args->N_neighbors, x, y, z, invert_x, invert_y, invert_z, args->neighbor_offset);
        }
    }else{
            #ifdef DEBUG
        fprintf(stderr, "normal.\n");
            #endif
            /*We're not; copy in*/
        if(d != NULL) {
    put_density_into_metadensity(meta_density, d, args->N_neighbors, x, y, z, invert_x, invert_y, invert_z, args->neighbor_offset);
        }
    }
      /*We always put epsilon into the center, so do it on every loop*/
    put_density_into_metadensity(meta_epsilon, epsilon, args->N_neighbors, x, y, z, invert_x, invert_y, invert_z, args->neighbor_offset);
  }
      }
    }
  }
  /*if specified, we expand the calculation box here.*/
  if(args->expand_density_boxes != NULL) {
    fprintf(stderr, "Expanding meta_density box in %d, %d, %d.\n", args->expand_density_boxes[0], args->expand_density_boxes[1], args->expand_density_boxes[2]);
    struct density* t;
    t = expand_density_box_with_value(meta_density, args->expand_density_boxes, 0, &e);
    if(e != NULL) {
      fprintf(stderr, "ERROR: unable to expand meta_density out by %d,%d,%d grid sites: %s\n", args->expand_density_boxes[0], args->expand_density_boxes[1], args->expand_density_boxes[2], e->message);
      return 70;
    }
    free_density(meta_density);
    meta_density = t;
    fprintf(stderr, "Expanding meta_density box in %d, %d, %d.\n", args->expand_density_boxes[0], args->expand_density_boxes[1], args->expand_density_boxes[2]);
    t = expand_density_box(meta_epsilon, args->expand_density_boxes, &e);
    if(e != NULL) {
      fprintf(stderr, "ERROR: unable to expand meta_epsilon out by %d,%d,%d grid sites: %s\n", args->expand_density_boxes[0], args->expand_density_boxes[1], args->expand_density_boxes[2], e->message);
      return 70;
    }
    free_density(meta_epsilon);
    meta_epsilon = t;
  }
  /*write sweeps and slices for density*/
  write_slices(args, meta_density, "charge_density_slices", "Bohr", "e Bohr^2/Bohr^3");
  write_sweeps(args, meta_density, "charge_density_sweeps");
  write_slices(args, meta_epsilon, "epsilon_slices", "Bohr", "F m");
  write_sweeps(args, meta_epsilon, "epsilon_sweeps");

  /*A check to verify that this works properly.*/
  struct density* dcheck2 = new_density(d->N, d->dx, 0);
  get_density_from_metadensity(meta_density, dcheck2, args->expand_density_boxes, args->N_neighbors, 0, 0, 0, args->neighbor_offset);

  if(args->center_wf_list == NULL) {
  struct density* dcheck = new_density(d->N, d->dx, 0);
  get_density_from_metadensity(meta_density, dcheck, args->expand_density_boxes, args->N_neighbors, args->N_neighbors[0], args->N_neighbors[1], args->N_neighbors[2], args->neighbor_offset);
  fprintf(stderr, "CHECKING that the center density is zero:");
  for(int k=0; k<dcheck->N[2]; k++) {
    for(int j=0; j<dcheck->N[1]; j++) {
      for(int i=0; i<dcheck->N[0]; i++) {
  if(
     dcheck->storage[get_density_storage_index_(i, j, k, dcheck->N[0], dcheck->N[1])]
     != 0) {
    fprintf(stderr, "ERROR: dcheck[%d,%d,%d] (%g) != 0\n", 
      i, j, k,
      dcheck->storage[get_density_storage_index_(i, j, k, dcheck->N[0], dcheck->N[1])]);
  }
      }
    }
  }
  free_density(dcheck);
  }
  fprintf(stderr, "DONE\n");

  /*Density needs to be multiplied by qe and divided by grid spacing 
    (assuming all are same) in m*/
  if(((meta_density->dx)[0] != (meta_density->dx)[1])
     || ((meta_density->dx)[0] != (meta_density->dx)[2])) {
    fprintf(stderr, "ERROR: not all grid spacings are the same! (dx=%g dy=%g dz=%g)\n", (meta_density->dx)[0], (meta_density->dx)[1], (meta_density->dx)[2]);
    return 30;
  }
  printf("grid site distances are the same, and are %g (%gnm)\n", (meta_density->dx)[0], bohr_to_nm((meta_density->dx)[0]));
  /*This comes about from the following
   * the density is 1/dx[0]^3.  The differentiation causes rho to be *multiplied* by dx[0]^2.
   * So dx[0]^2/dx[0]^3 = 1/dx[0]
   */
  density_mul_by_const(meta_density, mag_q_e_in_C/(bohr_to_m((meta_density->dx)[0])));

  /**\note if we're calling get_potential_for_density directly, we're not using a dielectric term, so it's NULL*/
  struct density *phi;
  GError *err5 = NULL;
  if(args->use_inhomogeneous_dielectric) {
    phi = get_potential_for_density_inhomog_simple(meta_density, meta_epsilon, 1e-4, 0, args->save_inhomog_iterations, &err5);
    rhoprime = new_density(meta_density->N, meta_density->dx, 0);
  }else{
    phi = get_potential_for_density(meta_density, meta_epsilon, &rhoprime, NULL, NULL, &err5);
  }
  if(err5 != NULL) {
    fprintf(stderr, "ERROR getting potential: %s\n", err5->message);
    return 100;
  }
  if(phi == NULL) {
    fprintf(stderr, "Error allocating memory for get_potential_for_density!\n");
    return 5;
  }
  /*write sweeps and slices for phi*/
  write_slices(args, phi, "potential_slices", "Bohr", "V");
  write_sweeps(args, phi, "potential_sweeps");
  write_slices(args, rhoprime, "rhoprime_slices", "Bohr", "V/Bohr");
  write_sweeps(args, rhoprime, "rhoprime_sweeps");

  struct density* phi_center = new_density(d->N, d->dx, 0);
  if(phi_center == NULL) {
    fprintf(stderr, "ERROR: failed to get memory for central portion of the density!\n");
  }
  get_density_from_metadensity(phi, phi_center, args->expand_density_boxes, args->N_neighbors, args->N_neighbors[0], args->N_neighbors[1], args->N_neighbors[2], args->neighbor_offset);

  /*write sweeps and slices for phi*/
  write_slices(args, phi_center, "center_potential_slices", "Bohr", "V");
  write_sweeps(args, phi_center, "center_potential_sweeps");

  /**Convert the output potential to electron interaction energy in Hartree
   * per electron.
   * Multiplying by -1.0 gets us to eV; then we need to multiply by hartree/eV
  */
  density_mul_by_const(phi_center, -ev_in_hartree);

  struct dotcode_grid* g = dotcode_grid_from_density(phi_center, &e);
  if(e != NULL) {
    fprintf(stderr, "ERROR: unable to convert density to a dotcode grid for writing: %s\n", e->message);
    return 37;
  }
  
  e = dotcode_write_typed_file(DOTCODE_GRID, g, args->outfile);
  if(e != NULL) {
    char* pathstr = g_file_get_uri(args->outfile);
    fprintf(stderr, "ERROR: unable to write out the final potential (%s): %s\n", pathstr, e->message);
    g_free(pathstr);
    g_error_free(e);
    return 21;
  }
  
  if(d != meta_density) {
    free_density(d);
    free_density(meta_density);
  }else{
    free_density(d);
  }
  free_density(epsilon);
  free_density(meta_epsilon);
  /*free_density(phi);*/
  return 0;
}

void write_sweeps(struct args* args, struct density* d, char* subdir) {
  /*Output the sweeps if requested.*/
  if((!(args->sweep_x)) && (!(args->sweep_y)) && (!(args->sweep_z))) return;
  GFile* f = g_file_new_for_path(subdir);
  GError* e=NULL;
  g_file_make_directory_with_parents(f, NULL, &e);
  if((e!=NULL) && (e->code != G_IO_ERROR_EXISTS)){
    char* fn = g_file_get_uri(f);
    if(e == NULL) {
      fprintf(stderr, "ERROR: write_sweeps(%s): unable to make directory for writing sweep (%s)", subdir, fn);
    }else{
      fprintf(stderr, "ERROR: write_sweeps(%s): unable to make directory for writing sweep (%s): %s (%d/%d)", subdir, fn, e->message, e->code, G_IO_ERROR_EXISTS);
    }
    free(fn);
    return;
  }
  char* s;
  if(args->sweep_x) {
    e=NULL;
    s=NULL;
    write_density_sweep(X, d, &s, &e, f, args->pixel_strategy);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      free(s);
    }
  }
  if(args->sweep_y) {
    e=NULL;
    s=NULL;
    write_density_sweep(Y, d, &s, &e, f, args->pixel_strategy);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      free(s);
    }
  }
  if(args->sweep_z) {
    e=NULL;
    s=NULL;
    write_density_sweep(Z, d, &s, &e, f, args->pixel_strategy);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      free(s);
    }
  }
}

void write_slices(struct args* args, struct density* d, char* subdir, const char* coord_units, const char* density_units) {
  GError* e=NULL;
  /*Output the sweeps if requested.*/
  GFile* f = g_file_new_for_path(subdir);
  char* s;
  g_file_make_directory_with_parents(f, NULL, &e);
  if((e!=NULL) && (e->code != G_IO_ERROR_EXISTS)){
    char* fn = g_file_get_uri(f);
    if(e == NULL) {
      fprintf(stderr, "ERROR: write_slices(%s): unable to make directory for writing slices (%s)", subdir, fn);
    }else{
      fprintf(stderr, "ERROR: write_slices(%s): unable to make directory for writing slices (%s): %s (%d/%d)", subdir, fn, e->message, e->code, G_IO_ERROR_EXISTS);
    }
    free(fn);
    return;
  }
  if(args->save_ascii_values) {
    e=NULL;
    GString *ess = g_string_new("");
    g_string_printf(ess, "ascii_%s.pdata", subdir);
    fprintf(stderr, "Writing dump to %s: ", ess->str);
    GFile* av = g_file_get_child(f, ess->str);
    write_density_slice(X, NULL, d, &s, &e, av, args->sepchar, coord_units, density_units);
    if(e != NULL) {
      fprintf(stderr, "ERROR dumping ascii to file %s: %s\n", s, e->message);
    }
    fprintf(stderr, "%s\n", s);
    g_string_free(ess, TRUE);
    g_free(s);
  }
  if(args->povray) {
    e=NULL;
    GString *ess = g_string_new("");
    g_string_printf(ess, "ascii_%s.pov", subdir);
    fprintf(stderr, "Writing povray scene to %s: ", ess->str);
    GFile* av = g_file_get_child(f, ess->str);
    write_density_povray(d, av, args->pixel_strategy, &e);
    if(e != NULL) {
      fprintf(stderr, "ERROR dumping ascii to file %s: %s\n", s, e->message);
    }
    fprintf(stderr, "%s\n", s);
    g_string_free(ess, TRUE);
    g_free(s);
  }
  e=NULL;
  if((args->xslices == NULL) && (args->yslices == NULL) && (args->zslices == NULL)) {
    fprintf(stderr, "No slices to write. going back to main routine.\n");
    return;
  }
  GList *l;
  for(l=args->xslices; l!=NULL; l=l->next) {
    printf("Writing slice x=%d: ", *(int*)(l->data));
    s=NULL;
    e=NULL;
    write_density_slice(X, (int*)(l->data), d, &s, &e, f, '\t', coord_units, density_units);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      g_free(s);
    }
  }
  for(l=args->yslices; l!=NULL; l=l->next) {
    printf("Writing slice y=%d: ", *(int*)(l->data));
    s=NULL;
    e=NULL;
    write_density_slice(Y, (int*)(l->data), d, &s, &e, f, '\t', coord_units, density_units);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      g_free(s);
    }
  }
  for(l=args->zslices; l!=NULL; l=l->next) {
    printf("Writing slice z=%d: ", *(int*)(l->data));
    s=NULL;
    e=NULL;
    write_density_slice(Z, (int*)(l->data), d, &s, &e, f, '\t', coord_units, density_units);
    if(e!=NULL){
      printf("->%s. FAILED: %s.\n", s, e->message);
    }else{
      printf("->%s. done.\n", s);
    }
    if(s!=NULL) {
      g_free(s);
    }
  }
}