First implementation of the potential generation by using

convolution theorem with a given arbitrary potential kernel.

GridDxFormat/grid_reader.c

+#include <stdio.h>
+#include <errno.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <stddef.h>
+#include <string.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <unistd.h>
+#include <ctype.h>
+#include "grid_reader.h"
+typedef struct
+{
+    size_t size;
+    size_t nx, ny, nz;
+    double* data;
+    double*  origin;
+    double** basis;
+    
+} grid;
+
+void remove_new_line(char* str)
+{
+    while(*str != '\n')
+        ++str;
+    *str = '\0';
+}
+
+void remove_white_space(char* str)
+{
+    size_t len = strlen(str);
+    char* ptr = str;
+    //int count = 0;
+    while(*ptr != '\0')
+    {
+        if(*ptr == ' ' || *ptr == '\t')
+            ptr++;
+        else
+            break;
+    }
+    memmove(str, ptr, sizeof(char)*(len-(size_t)(ptr-str)+1));
+    len = strlen(str); //new length
+    ptr = str+len-1;
+    while(*ptr == ' ')
+        ptr--;
+        //truncated the trailing white space
+    *(ptr+1) = '\0';
+}
+
+void* gridformat_reader(const char* name)
+{
+    printf("Reading from file %s\n", name);
+    char buffer[512];
+    char delim[]=" \t";
+    int set_count = 0;
+    FILE* infile = fopen(name, "r");
+    grid* grids = (grid*)malloc(sizeof(grid));
+
+    grids->nx = 0;
+    grids->ny = 0;
+    grids->nz = 0;
+    grids->size = 0;
+    grids->origin = (double*)malloc(sizeof(double)*3);
+
+    grids->basis  = (double**)malloc(sizeof(double*)*3);
+    for(size_t i = 0; i < 3; ++i)
+        grids->basis[i] = (double*)malloc(sizeof(double)*3);
+
+    grids->data = NULL;
+
+     
+    if(infile == NULL)
+    {
+        perror("Error: ");
+        printf("Failure opening file %s\n", name);
+        return NULL;
+    }
+    size_t count = 0;
+    size_t count_data = 0;
+
+    while(fgets(buffer, 512, infile) != NULL)
+    {
+        char tmp[512];
+        char str[512];
+        strcpy(tmp, buffer);
+        remove_new_line(tmp);
+        remove_white_space(tmp); //strip the possible newline and puts '\0'
+        //skip comment
+        if(buffer[0] == '#')
+            continue;
+        strcpy(str, tmp);
+
+        if(strlen(str) == 0)
+            continue;
+        char* rest = strtok(tmp, " ");
+
+        if(rest == NULL)
+            continue; 
+        if(!strcmp(rest, "object") && set_count == 0)
+        {
+            strtok(str," ");
+            strtok(NULL, " ");
+            strtok(NULL, " ");
+            strtok(NULL, " ");
+            strtok(NULL, " ");
+
+            rest = strtok(NULL, " ");
+            grids->nx = (size_t)atoi(rest);
+
+            rest = strtok(NULL, " ");
+            grids->ny = (size_t)atoi(rest);
+
+            rest = strtok(NULL, " ");
+            grids->nz = (size_t)atoi(rest);
+
+            grids->size = grids->nx*grids->ny*grids->nz;
+            grids->data = (double*)malloc(sizeof(double)*(grids->size));
+            //grids->data = (double***)malloc(sizeof(double**)*(grids->nx));
+            if(grids->data == NULL)
+            {
+                perror("Error in malloc:");
+                //printf("Error message %s\n", strerror(errno));
+                return NULL;
+            }
+            /*
+            for(size_t i = 0; i < grids->nx; ++i)
+            {
+                
+                grids->data[i] = (double**)malloc(sizeof(double*)*(grids->ny));
+                if(grids->data[i] == NULL)
+                {
+                    perror("Error in mallc:");
+
+                    for(size_t j = 0; j < i; ++j)
+                        free(grids->data[j]);
+
+                    free(grids->data);
+
+                    return NULL;
+                }
+            }
+            for(size_t i = 0; i < grids->nx; ++i)
+            {
+                for(size_t j = 0; j < grids->ny; ++j)
+                {
+                    grids->data[i][j] = (double*)malloc(sizeof(double)*(grids->nz));
+
+                    if(grids->data[i][j] == NULL)
+                    {
+                        perror("Error in malloc:");
+
+                        for(size_t i1= 0; i1 < i; ++i1)
+                            for(size_t j1 = 0; j1 < grids->ny; ++j1)
+                                free(grids->data[i1][j1]);
+
+                        for(size_t j1 = 0; j1 < j; ++j1)
+                                free(grids->data[i][j1]);
+
+                        for(size_t i1 = 0; i1 < grids->nx; ++i1)
+                             free(grids->data[i1]);
+
+                        free(grids->data);
+                        return NULL;
+                    }
+                }
+            }*/
+            set_count = 1;
+        }
+        else if(!strcmp(rest, "origin"))
+        {
+            strtok(str," ");
+            grids->origin[0] = (double)atof(strtok(NULL, " "));
+            grids->origin[1] = (double)atof(strtok(NULL, " "));
+            grids->origin[2] = (double)atof(strtok(NULL, " "));
+        }
+        else if(!strcmp(rest, "delta")) 
+        {
+            strtok(str," ");
+            grids->basis[count][0] = (double)atof(strtok(NULL, " "));
+            grids->basis[count][1] = (double)atof(strtok(NULL, " "));
+            grids->basis[count][2] = (double)atof(strtok(NULL, " "));
+            ++count;
+        }
+        else if(isdigit(str[0]) || isdigit(str[1]))
+        {
+            rest = strtok(str, delim);
+            grids->data[count_data++] = (double)atof(rest);
+
+            rest = strtok(NULL, delim);
+            if(rest != NULL)
+                grids->data[count_data++] = (double)atof(rest);
+            else
+                continue;
+
+            rest = strtok(NULL, delim);
+            if(rest != NULL)
+                grids->data[count_data++] = (double)atof(rest);
+            else
+                continue;
+        }
+    }
+    if(count_data != grids->size)
+    {
+        printf("Invalid file, data and count is not consistent.\n");
+        grid_destroy(grids);
+        //free(grids);
+        return NULL;
+    }
+    printf("Sucessfully reading file %s\n", name);
+    
+    if(fclose(infile) != 0)
+        perror("Error in closing file:");
+
+    return grids;
+}
+
+void grid_destroy(void* src)
+{
+    free(((grid*)src)->data);
+    for(size_t i = 0; i < 3; ++i)
+        free(((grid*)src)->basis[i]);
+    free(((grid*)src)->origin);
+    free(((grid*)src)->basis);
+    free(src);
+}
+
+int gridformat_writer(char* name, void* src)
+{
+    ssize_t bytes_write;
+    printf("Starting writing files %s\n", name);
+    assert(src);
+    grid* g = (grid*)src;
+    int fildes = open(name, O_RDWR | O_CREAT, 0664);
+    if(fildes == -1)
+    {
+        perror("Error in opeing file for writing");
+        //exit(errno);
+        //return -1;
+        return errno;
+    }
+    bytes_write = write(fildes, "#OpenDx file generated from grids reader/writer program\n",56);
+    if(bytes_write < -1)
+    {
+        return errno;
+    }
+    else if(bytes_write < 56)
+    {
+        return 9999;
+    }
+
+    char buff[1024];
+    ssize_t len;
+    sprintf(buff, "object 1 class gridpositions counts %lu %lu %lu\n", g->nx, g->ny, g->nz);
+    len = strlen(buff);
+    bytes_write = write(fildes, buff, len);
+
+    if(bytes_write < -1)
+    {
+        return errno;
+    }
+    else if(bytes_write < len)
+    {
+        return 9999;
+    }
+
+    sprintf(buff, "origin %lf %lf %lf\n", g->origin[0], g->origin[1], g->origin[2]);
+    len = strlen(buff);
+    bytes_write = write(fildes, buff, len);
+
+    if(bytes_write < -1)
+    {
+        return errno;
+    }
+    else if(bytes_write < len)
+    {
+        return 9999;
+    }
+
+    for(size_t i = 0; i < 3; i++)
+    {
+        sprintf(buff, "delta %lf %lf %lf\n", g->basis[i][0], g->basis[i][1], g->basis[i][2]);
+        len = strlen(buff);
+        bytes_write = write(fildes, buff, len);
+        if(bytes_write < -1)
+        {
+            return errno;
+        }
+        else if(bytes_write < len)
+        {
+            return 9999;
+        }
+
+    }
+
+    sprintf(buff, "object 2 class gridconnections counts %lu %lu %lu\n", g->nx, g->ny, g->nz); 
+    len = strlen(buff);
+    bytes_write = write(fildes, buff, len);
+    if(bytes_write < -1)
+    {
+        return errno;
+    }
+    else if(bytes_write < len)
+    {
+        return 9999;
+    }
+
+    sprintf(buff, "object 3 class array type float rank 0 items %lu data follows\n", g->size);
+    len = strlen(buff);
+    bytes_write = write(fildes, buff, len);
+    if(bytes_write < -1)
+    {
+        return errno;
+    }
+    else if(bytes_write < len)
+    {
+        return 9999;
+    }
+
+    size_t lines = g->size/3;
+    for(size_t l = 0; l < lines; ++l)
+    {
+        sprintf(buff, "%.15lf %.15lf %.15lf\n", g->data[3*l], g->data[3*l+1], g->data[3*l+2]);
+        len = strlen(buff);
+        bytes_write = write(fildes, buff, len);
+        if(bytes_write < -1)
+        {
+            return errno;
+        }
+        else if(bytes_write < len)
+        {
+            return 9999;
+        }
+    }
+    if(g->size%3 == 1)
+    {
+        sprintf(buff, "%.15lf\n", g->data[g->size-1]);
+        len = strlen(buff);
+        bytes_write = write(fildes, buff, len);
+        if(bytes_write < -1)
+        {
+            return errno;
+        }
+        else if(bytes_write < len)
+        {
+            return 9999;
+        }
+    }
+    else if(g->size%3 == 2)
+    {
+        sprintf(buff, "%.15lf %.15lf\n", g->data[g->size-2], g->data[g->size-1]);
+        len = strlen(buff);
+        bytes_write = write(fildes, buff, len);
+        if(bytes_write < -1)
+        {
+            return errno;
+        }
+        else if(bytes_write < len)
+        {
+            return 9999;
+        }
+    }
+    if(close(fildes) == -1)
+    {
+        int err = errno;
+        perror("Error in closing files");
+        //exit(errno);
+        return err;
+    }
+    return 0;
+}
+
+void set_grid_data(void* src, compute_type compute, void* kernel)
+{
+    grid* g = (grid*)src; 
+    for(size_t i = 0; i < g->nx; ++i)
+    {
+        for(size_t j = 0; j < g->ny; ++j)
+        {
+            for(size_t k = 0; k < g->nz; ++k)
+            {
+                size_t idx = k + g->nz * j + g->nz * g->ny * i;
+                g->data[idx] = (*compute)(i, j, k, g->origin, g->basis, kernel);
+            }
+        }
+    } 
+}
+
+double* get_grid_origin(void* src)
+{
+    return ((grid*)src)->origin;
+}
+
+double** get_grid_basis(void* src)
+{
+    return ((grid*)src)->basis;
+}
+
+double* get_grid_data(void* src)
+{
+    return ((grid*)src)->data;
+}
+
+size_t get_nx(void* src)
+{
+    return ((grid*)src)->nx;
+}
+
+size_t get_ny(void* src)
+{
+    return ((grid*)src)->ny;
+}
+
+size_t get_nz(void* src)
+{
+    return ((grid*)src)->nz;
+}
+
+size_t get_size(void* src)
+{
+    return ((grid*)src)->size;
+}
+
+void* copy_grid_format(void* ptr)
+{
+    grid* src  = (grid*)ptr;
+    grid* dest = (grid*)malloc(sizeof(grid));
+    dest->nx = src->nx;
+    dest->ny = src->ny;
+    dest->nz = src->nz;
+    dest->size = src->size;
+    dest->data = (double*)malloc(sizeof(double)*dest->size);
+    dest->origin = (double*)malloc(sizeof(double)*3);
+
+    dest->basis  = (double**)malloc(sizeof(double*)*3);
+    for(size_t i = 0; i < 3; ++i)
+        dest->basis[i] = (double*)malloc(sizeof(double)*3);
+
+    for(size_t i = 0; i < 3; ++i)
+    {
+        for(size_t j = 0; j < 3; ++j)
+        {
+            dest->basis[i][j] = src->basis[i][j];
+        }
+        dest->origin[i] = src->origin[i];
+    }
+    return dest;   
+}

Kernel/debye_huckel_kernel.c

+#include "debye_huckel_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct debye_huckel_kernel
+{
+    double screen_length; 
+    double shift;
+};
+
+void* debye_huckel_create(char* name)
+{
+    void* src = malloc(sizeof(debye_huckel_kernel));
+    FILE* infile = fopen(name, "r");
+    if(infile == NULL)
+    {
+        perror("Error in opening file:");
+        free(src);
+        return NULL;
+    }
+    char buffer[512];
+    while(fgets(buffer, 512, infile) != NULL)
+    {
+         char* tmp = strtok(buffer, " ");
+         if(strcmp(tmp, "ScreenLength") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((debye_huckel_kernel*)src)->screen_length = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else if(strcmp(tmp, "Shift") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((debye_huckel_kernel*)src)->shift = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else
+         {
+             printf("Unrecognized pattern\n");
+             free(src);
+             return NULL;
+         }
+    }
+    return src;
+}
+
+double debye_huckel_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+    debye_huckel_kernel* kernel = (debye_huckel_kernel*)src;
+    double pos[3];
+    for(size_t dim = 0; dim < 3; ++dim)
+        pos[dim] = basis[dim][0]*i + basis[dim][1]*j + basis[dim][2]*k + origin[dim];
+    double d = sqrt(pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]);    
+    return exp(-d/kernel->screen_length)/(d+kernel->shift);   
+}
+
+void debye_huckel_destroy(void* src)
+{
+    free(src);
+}

LinkList/link_list.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stddef.h>
+#include "link_list.h"
+struct node
+{
+    void* key;
+    void* val;
+    struct node* next;
+    //my_link_list* ptr;
+};
+
+struct link_list
+{
+    //append one empty node
+    node* first;
+    //node last; 
+
+    default_constructor_type default_constructor_key, default_constructor_val;
+    copy_constructor_type copy_constructor_key, copy_constructor_val;
+    destructor_type destructor_key, destructor_val;
+    compare_type compare;
+};
+
+void* link_list_create(default_constructor_type def_key, default_constructor_type def_val,
+                          copy_constructor_type cpy_key,    copy_constructor_type cpy_val,
+                          destructor_type destruct_key,     destructor_type destruct_val, compare_type comp)
+{
+    void* src = malloc(sizeof(link_list));
+    link_list* ptr = (link_list*)src;
+    ptr->default_constructor_key = def_key;
+    ptr->default_constructor_val = def_val;
+    ptr->copy_constructor_key    = cpy_key;
+    ptr->copy_constructor_val    = cpy_val;
+    ptr->destructor_key          = destruct_key;
+    ptr->destructor_val          = destruct_val;
+    ptr->compare                 = comp;
+    ptr->first                   = (node*)malloc(sizeof(node));
+    ptr->first->next = NULL;
+    return src;
+}
+
+//on success return 0; otherwise return -1;
+int insert_list(void* list, int pos, void* key, void* val)
+{
+    link_list* ptr = (link_list*)list;
+    node* tmp = ptr->first;
+    int count = 0;
+    while(tmp->next != NULL && count != pos)
+    {
+        tmp = tmp->next;
+        ++count;
+    }
+    node* insert_item = (node*)malloc(sizeof(node));
+    
+    if(insert_item != NULL)
+    {
+        insert_item -> key = ptr->copy_constructor_key(key);
+        insert_item -> val = ptr->copy_constructor_val(val);
+        node* old = tmp->next;
+        tmp->next = insert_item;
+        insert_item->next = old;
+        return 0;
+    }
+    else
+    {
+        perror("Error malloc:");
+        return -1;
+    }
+}
+
+int append_list(void* list, void* key, void* val)
+{
+    link_list *ptr = (link_list*)list;
+    node* tmp = ptr->first;
+    while(tmp->next != NULL)
+    {   
+        tmp = tmp->next;
+    }
+    node* item = (node*)malloc(sizeof(node));
+    if(item == NULL)
+    {
+       perror("Error in malloc:");
+       return -1;
+    }
+    tmp -> next = item;
+    item -> val = ptr->copy_constructor_val(val);
+    item -> key = ptr->copy_constructor_key(key);
+    item->next = NULL;
+    return 0;
+    //return tmp; //return NULL if not found
+}
+
+void* find_list(void* list, void* key)
+{
+    link_list *ptr = (link_list*)list;
+    node* tmp = ptr->first->next;
+    while(tmp != NULL)
+    {
+        if(ptr->compare(tmp->key, key) == 0)
+            break;
+        else
+            tmp = tmp->next;
+    }
+    if(tmp != NULL)
+        return tmp->val; //return NULL if not found
+    else
+        return NULL;
+}
+
+void clear_list(void* list)
+{
+    link_list *ptr = (link_list*)list;
+    node* tmp = ptr->first->next;
+    while(tmp != NULL)
+    {
+        ptr->destructor_val(tmp->val);
+        ptr->destructor_key(tmp->key);
+        node* old = tmp;
+        tmp = tmp->next;
+        free(old);
+    }
+    free(ptr->first);
+    free(ptr);
+}
+
+void* string_default_constructor()
+{
+    return malloc(sizeof(char)*64);
+}
+
+void* string_copy_constructor(void* src)
+{
+    size_t len = strlen(src);
+    void* dest = malloc(sizeof(char)*(len+1));
+    strcpy(dest, src);
+    return dest;
+}
+
+void string_destructor(void* src)
+{
+    free(src);
+}
+
+int string_compare(void* a, void* b)
+{
+    return strcmp(a,b);
+}
+
+void string_print(void* src)
+{
+    link_list *ptr = (link_list*)src;
+    node* tmp = ptr->first->next;
+    while(tmp != NULL) 
+    {
+        printf("Key: %s, Val: %s\n", (char*)(tmp->key), (char*)(tmp->val));
+        tmp = tmp->next;
+    }
+}
+
+void string_functor_print(void* src)
+{
+    link_list *ptr = (link_list*)src;
+    node* tmp = ptr->first->next;
+    while(tmp != NULL)
+    {
+        printf("Key: %s, Val: %p\n", (char*)(tmp->key), (char*)(tmp->val));
+        tmp = tmp->next;
+    }
+}
+
+void* shallow_default_constructor(void* src)
+{
+    return src;
+}
+
+void* shallow_copy_constructor(void* src)
+{
+    return src;
+}
+
+void shallow_destructor(void* src)
+{
+}
+
+

fftw/fftw.c

+#include "fftw.h"
+fftw_plan fftwr2c_on_grid(fftw_complex *out, size_t nx, size_t ny, size_t nz, double* data)
+{
+     return fftw_plan_dft_r2c_3d(nx, ny, nz, data, out, FFTW_ESTIMATE);
+}
+
+fftw_plan fftwc2r_on_grid(double* data, size_t nx, size_t ny, size_t nz, fftw_complex *out)
+{
+     return fftw_plan_dft_c2r_3d(nx, ny, nz, out, data, FFTW_ESTIMATE);
+}
+

inc/debye_huckel_kernel.h

+#ifdef KERNEL_TYPE
+ComputeType(debye_huckel, debye_huckel_compute);
+CreateType(debye_huckel, debye_huckel_create);
+DestroyType(debye_huckel, debye_huckel_destroy);
+#else
+#ifndef DEBYE_HUCKEL_H_
+#define DEBYE_HUCKEL_H_
+#include <stddef.h>
+typedef struct debye_huckel_kernel debye_huckel_kernel;
+void* debye_huckel_create(char* name);
+double debye_huckel_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src);
+void debye_huckel_destroy(void*);
+#endif
+#endif

inc/fftw.h

+#ifndef FFTW_H_
+#define FFTW_H_
+#include <fftw3.h>
+fftw_plan fftwr2c_on_grid(fftw_complex *out, size_t nx, size_t ny, size_t nz, double* data);
+fftw_plan fftwc2r_on_grid(double* data, size_t nx, size_t ny, size_t nz, fftw_complex *out);
+#endif

inc/grid_reader.h

+#ifndef GRID_READER_H_
+#define GRID_READER_H_
+#include <stddef.h>
+void remove_new_line(char* str);
+void remove_white_space(char* str);
+void* gridformat_reader(const char* name);
+void grid_destroy(void* src);
+int gridformat_writer(char* name, void* g);
+double* get_grid_origin(void* src);
+double** get_grid_basis(void* src);
+double* get_grid_data(void* src);
+size_t get_nx(void* src);
+size_t get_ny(void* src);
+size_t get_nz(void* src);
+size_t get_size(void* src);
+void* copy_grid_format(void*);
+typedef double (*compute_type)(size_t,size_t,size_t,double*,double**, void*);
+void set_grid_data(void* src, compute_type compute, void* kernel);
+#endif

inc/kernel.h

+#ifndef KERNEL_H_
+#define KERNEL_H_
+typedef void* (*create_type)(char*);
+typedef void  (*destroy_type)(char*);
+#endif
+#include "debye_huckel_kernel.h"

inc/link_list.h

+#ifndef LINK_LIST_H_
+#define LINK_LIST_H_
+
+typedef struct node node;
+typedef struct link_list link_list;
+typedef void* (*default_constructor_type)();
+typedef void* (*copy_constructor_type)(void*);
+typedef void  (*destructor_type)(void*); 
+typedef int   (*compare_type)(void*, void*);
+typedef void  (*print_type)(void*);
+
+void* link_list_create(default_constructor_type def_key, default_constructor_type def_val,
+                         copy_constructor_type cpy_key,     copy_constructor_type cpy_val,
+                         destructor_type destruct_key,      destructor_type destruct_val, compare_type comp);
+
+
+int insert_list(void* list, int pos, void* key, void* val);
+int append_list(void* list, void* key, void* val);
+void* find_list(void* list, void* key);
+void clear_list(void* list);
+void* string_default_constructor();
+void* string_copy_constructor(void* src);
+void string_destructor(void* src);
+void string_print(void* src);
+void string_functor_print(void* src);
+int string_compare(void*, void*);
+void* shallow_default_constructor(void* src);
+void* shallow_copy_constructor(void* src);
+void shallow_destructor(void*);
+
+#endif

main.c

+#include "kernel.h"
+#include "link_list.h"
+#include "grid_reader.h"
+#include "fftw.h"
+#include <stdio.h>
+static void inner_product_grid(fftw_complex* dest, fftw_complex* src, size_t nx, size_t ny, size_t nz)
+{
+    double re, im;    
+    size_t size = nx * ny * (nz/2+1);
+    size_t V = nx*ny*nz;
+    for(size_t i = 0; i < size; ++i)
+    {
+        re = dest[i][0];
+        im = dest[i][1];
+        dest[i][0] = (re*src[i][0]-im*src[i][1])/V;
+        dest[i][1] = (re*src[i][1]+im*src[i][0])/V;
+    }
+}
+static void print_usage()
+{
+    printf("./exec input_density_dx_file potential_kernel_name kernel_init_file outpu_dx_file\n");
+}
+int main(int argc, char** argv)
+{
+
+    if(argc < 5)
+    {
+        printf("Invalid input arguments\n");
+        print_usage();
+        return -1;
+    }
+    link_list* compute_map = link_list_create(&shallow_default_constructor, &shallow_default_constructor,
+                                              &shallow_copy_constructor,    &shallow_copy_constructor,
+                                              &shallow_destructor,          &shallow_destructor, &string_compare);
+
+    link_list* create_map = link_list_create(&shallow_default_constructor, &shallow_default_constructor,
+                                              &shallow_copy_constructor,    &shallow_copy_constructor,
+                                              &shallow_destructor,          &shallow_destructor, &string_compare);
+
+    link_list* destroy_map = link_list_create(&shallow_default_constructor, &shallow_default_constructor,
+                                              &shallow_copy_constructor,    &shallow_copy_constructor,
+                                              &shallow_destructor,          &shallow_destructor, &string_compare);
+
+    #define CreateType(key, val) \
+    append_list(create_map, #key, &val)
+
+    #define ComputeType(key, val) \
+    append_list(compute_map,#key, &val)
+
+    #define DestroyType(key, val) \
+    append_list(destroy_map, #key, &val)
+
+    #define KERNEL_TYPE
+    #include "kernel.h"
+    #undef KERNEL_TYPE
+    #undef ComputeType 
+    #undef CreateType
+    #undef DestroyType
+
+    char* density_file  = argv[1];
+    char* kernel_name   = argv[2];
+    char* kernel_infile = argv[3];
+
+    void* density   = gridformat_reader(density_file);
+    void* potential = copy_grid_format(density);
+   
+    create_type creator = find_list(create_map, kernel_name);
+
+    void* kernel = (*creator)(kernel_infile);
+
+    compute_type functor = find_list(compute_map, kernel_name);
+    set_grid_data(potential, functor, kernel);
+    fftw_complex* density_out, *potential_out;
+    fftw_plan plan1, plan2;
+
+    density_out   = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (get_nx(density))*(get_ny(density))*(get_nz(density)/2+1));
+    potential_out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (get_nx(density))*(get_ny(density))*(get_nz(density)/2+1));
+
+    plan1 = fftwr2c_on_grid(density_out, get_nx(density), get_ny(density), get_nz(density), get_grid_data(density));
+
+    fftw_execute(plan1);
+    fftw_execute_dft_r2c(plan1, get_grid_data(potential), potential_out);
+
+    inner_product_grid(potential_out, density_out, (get_nx(density)), (get_ny(density)), (get_nz(density)));
+
+    plan2 = fftwc2r_on_grid(get_grid_data(potential), get_nx(density), get_ny(density), get_nz(density), potential_out);
+ 
+    fftw_execute(plan2);
+    gridformat_writer(argv[4], potential);
+   
+    //clean out all memory on the heap 
+    destroy_type destroy = find_list(destroy_map, kernel_name);
+    (*destroy)(kernel);
+
+    clear_list(compute_map);
+    clear_list(create_map);
+    clear_list(destroy_map);        
+
+    grid_destroy(density);
+    grid_destroy(potential);
+    fftw_free(density_out);
+    fftw_free(potential_out);
+
+    fftw_destroy_plan(plan1);
+    fftw_destroy_plan(plan2);
+    return 0;
+}