Fix a writer bug, the older one did not truncate the file. Add several...

Fix a writer bug, the older one did not truncate the file. Add several lenard-jones type interaction

GridDxFormat/grid_reader.c

@@ -232,7 +232,7 @@ int gridformat_writer(char* name, void* src)
     printf("Starting writing files %s\n", name);
     assert(src);
     grid* g = (grid*)src;
-    int fildes = open(name, O_RDWR | O_CREAT, 0664);
+    int fildes = open(name, O_RDWR | O_CREAT | O_TRUNC, 0664);
     if(fildes == -1)
     {
         perror("Error in opeing file for writing");
@@ -571,3 +571,15 @@ void* copy_grid_format(void* ptr)
     }
     return dest;   
 }
+
+void sum_grids(void* dest, int num, void** src)
+{
+    grid** in = (grid**)src;
+    grid* out = (grid*)dest;
+    for(size_t idx = 0; idx < out->size; idx++)
+    {
+        out->data[idx] = 0.;
+        for(size_t i = 0; i < (size_t)num; ++i)
+            out->data[idx] += in[i]->data[idx];
+    }
+}

Kernel/debye_huckel_kernel.c

@@ -64,7 +64,8 @@ double debye_huckel_compute(size_t i, size_t j, size_t k, double* origin, double
     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);   
+    double e = exp(-d/kernel->screen_length)/(d+kernel->shift);  
+    return (e > 100) ? 100: e;
 }
 
 void debye_huckel_destroy(void* src)

Kernel/lenard_jones_kernel.c

@@ -68,11 +68,11 @@ double lenard_jones_compute(size_t i, size_t j, size_t k, double* origin, double
     if(d > 0.)
     {
         double e = -(kernel->c6)/d+(kernel->c12)/(d*d);
-        return (e>20) ? 20 : e;
+        return (e>100) ? 100 : e;
     }
     else
     {
-        return 20;
+        return 100;
     }
 }
 

Kernel/lenard_jones_repulsion_kernel.c

+#include "lenard_jones_repulsion_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct lenard_jones_repulsion_kernel
+{
+    double c6; 
+    double c12;
+};
+
+void* lenard_jones_repulsion_create(char* name)
+{
+    void* src = malloc(sizeof(lenard_jones_repulsion_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_repulsion_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_repulsion_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_repulsion_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+    lenard_jones_repulsion_kernel* kernel = (lenard_jones_repulsion_kernel*)src;
+
+    double sigma2 = pow(kernel->c12/kernel->c6,0.33333333333333333333);
+    double epsilon = kernel->c6*kernel->c6/(4*kernel->c12);
+
+    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 d2 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+    double e;
+
+    if(d2<=pow(2.,0.33333333333333333)*sigma2)
+    {
+        d2=d2*d2*d2;
+        if(d2 > 0.)
+        {
+            e = (-(kernel->c6)/d2+(kernel->c12)/(d2*d2))+epsilon;
+            return (e>100) ? 100 : e;
+        }
+        else
+        {
+            return 100;
+        }
+    }
+    else
+        return 0;
+}
+
+void lenard_jones_repulsion_destroy(void* src)
+{
+    free(src);
+}

Kernel/lenard_jones_trunc_kernel.c

+#include "lenard_jones_trunc_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct lenard_jones_trunc_kernel
+{
+    double c6; 
+    double c12;
+    double scaling;
+};
+
+void* lenard_jones_trunc_create(char* name)
+{
+    void* src = malloc(sizeof(lenard_jones_trunc_kernel));
+    FILE* infile = fopen(name, "r");
+    if(infile == NULL)
+    {
+        perror("Error in opening file:");
+        free(src);
+        return NULL;
+    }
+    char buffer[512];
+    ((lenard_jones_trunc_kernel*)src)->c6 = 0.;
+    ((lenard_jones_trunc_kernel*)src)->c12 = 0.;
+    ((lenard_jones_trunc_kernel*)src)->scaling = 1.; 
+    while(fgets(buffer, 512, infile) != NULL)
+    {
+         char* tmp = strtok(buffer, " ");
+         if(strcmp(tmp, "C6") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((lenard_jones_trunc_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_trunc_kernel*)src)->c12 = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else if(strcmp(tmp, "Scaling") == 0)
+         {
+             tmp = strtok(NULL, " ");
+             if(tmp != NULL)
+                 ((lenard_jones_trunc_kernel*)src)->scaling = atof(tmp);
+             else
+             {
+                 printf("Error formating\n");
+                 free(src);
+                 return NULL;
+             }
+         }
+         else
+         {
+             printf("Unrecognized pattern\n");
+             free(src);
+             return NULL;
+         }
+    }
+    ((lenard_jones_trunc_kernel*)src)->c6 /= ((lenard_jones_trunc_kernel*)src)->scaling;
+    ((lenard_jones_trunc_kernel*)src)->c12 /= ((lenard_jones_trunc_kernel*)src)->scaling;
+    return src;
+}
+
+double lenard_jones_trunc_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+    lenard_jones_trunc_kernel* kernel = (lenard_jones_trunc_kernel*)src;
+
+    double sigma2 = pow(kernel->c12/kernel->c6,0.33333333333333333333);
+    double epsilon = kernel->c6*kernel->c6/(4*kernel->c12);
+
+    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 d2 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+    double e;
+    double rc2 = 2.5*2.5*sigma2;
+    //rc2 = rc2*rc2*rc2;
+
+    //d2 = d2*d2*d2;
+
+    if(d2<=pow(2.,0.33333333333333333)*sigma2)
+    {
+        d2=d2*d2*d2;
+        rc2=rc2*rc2*rc2;
+        if(d2 > 0.)
+        {
+            e = 2.*(-(kernel->c6)/d2+(kernel->c12)/(d2*d2))+epsilon-((-kernel->c6)/rc2+(kernel->c12)/(rc2*rc2));
+            return (e>200) ? 200 : e;
+        }
+        else
+        {
+            return 200;
+        }
+    }
+    else if(d2>rc2)
+        return 0;
+    else
+    {
+        d2=d2*d2*d2;
+        rc2=rc2*rc2*rc2;
+        return -(kernel->c6)/d2+(kernel->c12)/(d2*d2)-((-kernel->c6)/rc2+(kernel->c12)/(rc2*rc2));
+    }
+}
+
+void lenard_jones_trunc_destroy(void* src)
+{
+    free(src);
+}

inc/grid_reader.h

@@ -17,4 +17,5 @@ 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);
 void* resample_grid_data(void* src);
+void sum_grids(void* dest, int num, void** src);
 #endif

inc/lenard_jones_repulsion_kernel.h

+#ifdef KERNEL_TYPE
+ComputeType(lenard_jones_repulsion, lenard_jones_repulsion_compute);
+CreateType(lenard_jones_repulsion, lenard_jones_repulsion_create);
+DestroyType(lenard_jones_repulsion, lenard_jones_repulsion_destroy);
+#else
+#ifndef LENARD_JONES_REPULSION_H_
+#define LENARD_JONES_REPULSION_H_
+#include <stddef.h>
+typedef struct lenard_jones_repulsion_kernel lenard_jones_repulsion_kernel;
+void* lenard_jones_repulsion_create(char* name);
+double lenard_jones_repulsion_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src);
+void lenard_jones_repulsion_destroy(void*);
+#endif
+#endif

inc/lenard_jones_trunc_kernel.h

+#ifdef KERNEL_TYPE
+ComputeType(lenard_jones_trunc, lenard_jones_trunc_compute);
+CreateType(lenard_jones_trunc, lenard_jones_trunc_create);
+DestroyType(lenard_jones_trunc, lenard_jones_trunc_destroy);
+#else
+#ifndef LENARD_JONES_TRUNC_H_
+#define LENARD_JONES_TRUNC_H_
+#include <stddef.h>
+typedef struct lenard_jones_trunc_kernel lenard_jones_trunc_kernel;
+void* lenard_jones_trunc_create(char* name);
+double lenard_jones_trunc_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src);
+void lenard_jones_trunc_destroy(void*);
+#endif
+#endif

main.c

@@ -91,7 +91,7 @@ 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);
+    //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));