|
funkalicious 0.1
|
|
funkalicious 0.1
|
#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>
Include dependency graph for getChainPotential.c:Go to the source code of this file.
Functions | |
| 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) |
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
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
Definition at line 44 of file getChainPotential.c.
References bohr_to_m(), bohr_to_nm(), args::center_wf_list, density_mul_by_const(), DOTCODE_GRID, dotcode_grid1_get_from_gfile(), dotcode_grid1_get_property_grid(), dotcode_grid_from_density(), dotcode_write_typed_file(), dotcode_grid1::dx, density::dx, e, ev_in_hartree, expand_density_box(), expand_density_box_with_value(), args::expand_density_boxes, args::extra_neighbor, wavefunction::file, free_density(), get_density_for_list_of_wavefunctions(), get_density_from_metadensity(), get_density_storage_at_c_(), get_density_storage_index_(), get_potential_for_density(), get_potential_for_density_inhomog_simple(), args::invert_odd, load_wavefunctions(), mag_q_e_in_C, dotcode_grid1::N, density::N, args::N_neighbors, args::neighbor_offset, new_density(), args::outfile, process_args(), put_density_into_metadensity(), read_matdb_dotcode_gio(), args::save_inhomog_iterations, density::storage, args::use_inhomogeneous_dielectric, args::wf_list, write_slices(), and write_sweeps().
{
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;
}
Here is the call graph for this function:| void write_slices | ( | struct args * | args, |
| struct density * | d, | ||
| char * | subdir, | ||
| const char * | coord_units, | ||
| const char * | density_units | ||
| ) |
Definition at line 445 of file getChainPotential.c.
References e, args::pixel_strategy, args::povray, args::save_ascii_values, args::sepchar, write_density_povray(), write_density_slice(), X, args::xslices, Y, args::yslices, Z, and args::zslices.
Referenced by main().
{
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);
}
}
}
Here is the call graph for this function: Here is the caller graph for this function:Definition at line 387 of file getChainPotential.c.
References e, args::pixel_strategy, args::sweep_x, args::sweep_y, args::sweep_z, write_density_sweep(), X, Y, and Z.
Referenced by main().
{
/*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);
}
}
}
Here is the call graph for this function: Here is the caller graph for this function:
1.7.3 
