Add lenard-jones potential kernel

Kernel/lenard_jones_kernel.c

+#include "lenard_jones_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct lenard_jones_kernel
+{
+    double c6; 
+    double c12;
+};
+
+void* lenard_jones_create(char* name)
+{
+    void* src = malloc(sizeof(lenard_jones_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, "C6") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((lenard_jones_kernel*)src)->c6 = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else if(strcmp(tmp, "C12") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((lenard_jones_kernel*)src)->c12 = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else
+         {
+             printf("Unrecognized pattern\n");
+             free(src);
+             return NULL;
+         }
+    }
+    return src;
+}
+
+double lenard_jones_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+    lenard_jones_kernel* kernel = (lenard_jones_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]);
+    d = d*d*d*d*d*d;
+    if(d > 0.)
+    {
+        double e = -(kernel->c6)/d+(kernel->c12)/(d*d);
+        return (e>1e4) ? 1e4 : e;
+    }
+    else
+    {
+        return 1e4;
+    }
+}
+
+void lenard_jones_destroy(void* src)
+{
+    free(src);
+}

Makefile

+CC=gcc
+INC=-I/usr/include -I./inc
+ifeq ($(dbg),1)
+	CFLAGS=-std=c99 -g -Og -Wall -Wextra
+else
+	CFLAGS=-std=c99 -O3 -Wall -Wextra
+endif
+CLIB=-lfftw3 -lm
+OBJ=main.o GridDxFormat/grid_reader.o Kernel/debye_huckel_kernel.o LinkList/link_list.o fftw/fftw.o Kernel/lenard_jones_kernel.o
+ifeq ($(dbg),1)
+	EXE=gen_pot-dbg
+else
+	EXE=gen_pot
+endif
+
+all: ${EXE}
+
+${EXE}:${OBJ}
+	${CC} ${CFLAGS} ${INC} ${OBJ} -o ${EXE} ${CLIB}
+
+main.o: main.c
+	${CC} ${CFLAGS} ${INC} -c main.c -o main.o
+
+GridDxFormat/grid_reader.o: GridDxFormat/grid_reader.c
+	${CC} ${CFLAGS} ${INC} -c GridDxFormat/grid_reader.c -o GridDxFormat/grid_reader.o
+
+Kernel/debye_huckel_kernel.o:Kernel/debye_huckel_kernel.c
+	${CC} ${CFLAGS} ${INC} -c Kernel/debye_huckel_kernel.c -o Kernel/debye_huckel_kernel.o
+
+Kernel/lenard_jones_kernel.o:Kernel/lenard_jones_kernel.c
+	${CC} ${CFLAGS} ${INC} -c Kernel/lenard_jones_kernel.c -o Kernel/lenard_jones_kernel.o
+
+LinkList/link_list.o:LinkList/link_list.c
+	${CC} ${CFLAGS} ${INC} -c LinkList/link_list.c -o LinkList/link_list.o
+
+fftw/fftw.o:fftw/fftw.c
+	${CC} ${CFLAGS} ${INC} -c fftw/fftw.c -o fftw/fftw.o
+
+clean:
+	rm -f ${OBJ} ${EXE}

inc/kernel.h

@@ -4,3 +4,4 @@ typedef void* (*create_type)(char*);
 typedef void  (*destroy_type)(char*);
 #endif
 #include "debye_huckel_kernel.h"
+#include "lenard_jones_kernel.h"

inc/lenard_jones_kernel.h

+#ifdef KERNEL_TYPE
+ComputeType(lenard_jones, lenard_jones_compute);
+CreateType(lenard_jones, lenard_jones_create);
+DestroyType(lenard_jones, lenard_jones_destroy);
+#else
+#ifndef LENARD_JONES_H_
+#define LENARD_JONES_H_
+#include <stddef.h>
+typedef struct lenard_jones_kernel lenard_jones_kernel;
+void* lenard_jones_create(char* name);
+double lenard_jones_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src);
+void lenard_jones_destroy(void*);
+#endif
+#endif

main.c

@@ -16,6 +16,20 @@ static void inner_product_grid(fftw_complex* dest, fftw_complex* src, size_t nx,
         dest[i][1] = (re*src[i][1]+im*src[i][0])/V;
     }
 }
+static void shift_zero_frequency(size_t nx, size_t ny, size_t nz, double* data)
+{
+    for(size_t i = 0; i < nx; ++i)
+    {
+        for(size_t j = 0; j < ny; ++j)
+        {
+            for(size_t k = 0; k < nz; ++k)
+            {
+                size_t idx = k + nz * (j + ny * i);
+                *(data+idx) = (*(data+idx))*2.*(0.5-(double)((i+j+k)%2));
+            }
+        }
+    }
+}
 static void print_usage()
 {
     printf("./exec input_density_dx_file potential_kernel_name kernel_init_file outpu_dx_file\n");
@@ -34,8 +48,8 @@ int main(int argc, char** argv)
                                               &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);
+                                             &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,
@@ -70,6 +84,10 @@ int main(int argc, char** argv)
 
     compute_type functor = find_list(compute_map, kernel_name);
     set_grid_data(potential, functor, kernel);
+    gridformat_writer("kernel.dx", potential);
+    //shift_zero_frequency(get_nx(density), get_ny(density), get_nz(density), get_grid_data(density));
+    //shift_zero_frequency(get_nx(potential), get_ny(potential), get_nz(potential), get_grid_data(potential));
+
     fftw_complex* density_out, *potential_out;
     fftw_plan plan1, plan2;