From d6e3280f884f61fb020dc3f94dd99c869f941f5c Mon Sep 17 00:00:00 2001
From: Han-Yi Chou <hchou10@illinois.edu>
Date: Mon, 8 Oct 2018 14:50:26 -0500
Subject: [PATCH] Fix a writer bug, the older one did not truncate the file.
Add several lenard-jones type interaction
---
GridDxFormat/grid_reader.c | 14 ++-
Kernel/debye_huckel_kernel.c | 3 +-
Kernel/lenard_jones_kernel.c | 4 +-
Kernel/lenard_jones_repulsion_kernel.c | 94 +++++++++++++++++++
Kernel/lenard_jones_trunc_kernel.c | 123 +++++++++++++++++++++++++
inc/grid_reader.h | 1 +
inc/lenard_jones_repulsion_kernel.h | 14 +++
inc/lenard_jones_trunc_kernel.h | 14 +++
main.c | 2 +-
9 files changed, 264 insertions(+), 5 deletions(-)
create mode 100644 Kernel/lenard_jones_repulsion_kernel.c
create mode 100644 Kernel/lenard_jones_trunc_kernel.c
create mode 100644 inc/lenard_jones_repulsion_kernel.h
create mode 100644 inc/lenard_jones_trunc_kernel.h
diff --git a/GridDxFormat/grid_reader.c b/GridDxFormat/grid_reader.c
index 990ab14..909a77a 100644
--- a/GridDxFormat/grid_reader.c
+++ b/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];
+ }
+}
diff --git a/Kernel/debye_huckel_kernel.c b/Kernel/debye_huckel_kernel.c
index 1e97d9e..b7f840f 100644
--- a/Kernel/debye_huckel_kernel.c
+++ b/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)
diff --git a/Kernel/lenard_jones_kernel.c b/Kernel/lenard_jones_kernel.c
index f789212..12a7e48 100644
--- a/Kernel/lenard_jones_kernel.c
+++ b/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;
}
}
diff --git a/Kernel/lenard_jones_repulsion_kernel.c b/Kernel/lenard_jones_repulsion_kernel.c
new file mode 100644
index 0000000..735c45e
--- /dev/null
+++ b/Kernel/lenard_jones_repulsion_kernel.c
@@ -0,0 +1,94 @@
+#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);
+}
diff --git a/Kernel/lenard_jones_trunc_kernel.c b/Kernel/lenard_jones_trunc_kernel.c
new file mode 100644
index 0000000..e5dad1c
--- /dev/null
+++ b/Kernel/lenard_jones_trunc_kernel.c
@@ -0,0 +1,123 @@
+#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);
+}
diff --git a/inc/grid_reader.h b/inc/grid_reader.h
index 9b97394..e5c04d8 100644
--- a/inc/grid_reader.h
+++ b/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
diff --git a/inc/lenard_jones_repulsion_kernel.h b/inc/lenard_jones_repulsion_kernel.h
new file mode 100644
index 0000000..b4071bc
--- /dev/null
+++ b/inc/lenard_jones_repulsion_kernel.h
@@ -0,0 +1,14 @@
+#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
diff --git a/inc/lenard_jones_trunc_kernel.h b/inc/lenard_jones_trunc_kernel.h
new file mode 100644
index 0000000..081e754
--- /dev/null
+++ b/inc/lenard_jones_trunc_kernel.h
@@ -0,0 +1,14 @@
+#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
diff --git a/main.c b/main.c
index 57d5934..915ab66 100644
--- a/main.c
+++ b/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));
--
GitLab