From 53917a1287e51a73d79846c229df994c394cdf66 Mon Sep 17 00:00:00 2001
From: Han-Yi Chou <hchou10@illinois.edu>
Date: Wed, 17 Apr 2019 11:26:18 -0500
Subject: [PATCH] forgot to push some files
---
Kernel/gaussian_kernel.c | 91 +++++++++++++++++++++++++
Kernel/lenard_jones_kernel1.c | 101 +++++++++++++++++++++++++++
Kernel/onc_elec_kernel.c | 125 ++++++++++++++++++++++++++++++++++
Kernel/tabulated_pot_kernel.c | 100 +++++++++++++++++++++++++++
4 files changed, 417 insertions(+)
create mode 100644 Kernel/gaussian_kernel.c
create mode 100644 Kernel/lenard_jones_kernel1.c
create mode 100644 Kernel/onc_elec_kernel.c
create mode 100644 Kernel/tabulated_pot_kernel.c
diff --git a/Kernel/gaussian_kernel.c b/Kernel/gaussian_kernel.c
new file mode 100644
index 0000000..0cb1b4d
--- /dev/null
+++ b/Kernel/gaussian_kernel.c
@@ -0,0 +1,91 @@
+#include "gaussian_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct gaussian_kernel
+{
+ double sigma;
+ double amp;
+ double cutoff;
+};
+
+void* gaussian_create(char* name)
+{
+ void* src = malloc(sizeof(gaussian_kernel));
+ memset(src, 0, sizeof(gaussian_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, "Sigma") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((gaussian_kernel*)src)->sigma = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Amp") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((gaussian_kernel*)src)->amp = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Cutoff") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((gaussian_kernel*)src)->cutoff = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else
+ {
+ printf("Unrecognized pattern\n");
+ free(src);
+ return NULL;
+ }
+ }
+ fclose(infile);
+ return src;
+}
+
+double gaussian_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+ gaussian_kernel* kernel = (gaussian_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 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+ if(d > (kernel->cutoff)*(kernel->cutoff))
+ return 0;
+ else
+ return kernel->amp*exp(-d/((kernel->sigma)*(kernel->sigma)));
+}
+
+void gaussian_destroy(void* src)
+{
+ free(src);
+}
diff --git a/Kernel/lenard_jones_kernel1.c b/Kernel/lenard_jones_kernel1.c
new file mode 100644
index 0000000..d5eaa48
--- /dev/null
+++ b/Kernel/lenard_jones_kernel1.c
@@ -0,0 +1,101 @@
+#include "lenard_jones_kernel1.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct lenard_jones_kernel1
+{
+ double sigma;
+ double epsilon;
+ double cutoff;
+};
+
+void* lenard_jones_create1(char* name)
+{
+ void* src = malloc(sizeof(lenard_jones_kernel1));
+ 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, "Sigma") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((lenard_jones_kernel1*)src)->sigma = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Epsilon") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((lenard_jones_kernel1*)src)->epsilon = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Cutoff") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((lenard_jones_kernel1*)src)->cutoff = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else
+ {
+ printf("Unrecognized pattern\n");
+ free(src);
+ return NULL;
+ }
+ }
+ fclose(infile);
+ return src;
+}
+
+double lenard_jones_compute1(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+ lenard_jones_kernel1* kernel = (lenard_jones_kernel1*)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 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+ d = sqrt(d);
+ if(d > 0. && d <= kernel->cutoff*kernel->sigma)
+ {
+ d = kernel->sigma / d;
+ d = d*d*d*d*d*d;
+ //double e = -(kernel->c6)/d+(kernel->c12)/(d*d);
+ double e = kernel->epsilon*4.*d*(d-1.);
+ return (e>40.) ? 40. : e;
+ }
+ else if(d > kernel->cutoff*kernel->sigma)
+ return 0.;
+ else
+ {
+ return 40.;
+ }
+}
+
+void lenard_jones_destroy1(void* src)
+{
+ free(src);
+}
diff --git a/Kernel/onc_elec_kernel.c b/Kernel/onc_elec_kernel.c
new file mode 100644
index 0000000..3adadbc
--- /dev/null
+++ b/Kernel/onc_elec_kernel.c
@@ -0,0 +1,125 @@
+#include "onc_elec_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#define CoulombConst 332.15417
+
+struct onc_elec_kernel
+{
+ double kappa;
+ double Sz;
+ double z;
+ double q;
+};
+
+void* onc_elec_create(char* name)
+{
+ void* src = malloc(sizeof(onc_elec_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, "Kappa") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((onc_elec_kernel*)src)->kappa = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Sz") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((onc_elec_kernel*)src)->Sz = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Z") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((onc_elec_kernel*)src)->z = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else if(strcmp(tmp, "Charge") == 0)
+ {
+ tmp = strtok(NULL, " ");
+ if(tmp != NULL)
+ ((onc_elec_kernel*)src)->q = atof(tmp);
+ else
+ {
+ printf("Error formating\n");
+ free(src);
+ return NULL;
+ }
+ }
+ else
+ {
+ printf("Unrecognized pattern\n");
+ free(src);
+ return NULL;
+ }
+ }
+ fclose(infile);
+ return src;
+}
+
+static double computeDielectric(double d, onc_elec_kernel* pot)
+{
+ double z = pot->z;
+ double Sz = pot->Sz;
+ double tmp = exp(d/z);
+ return Sz*(1.-d*d/(z*z)*tmp/((tmp-1.)*(tmp-1.)));
+}
+/*
+static double computeDerivativeDielectric(double d, onc_elec_kernel* pot)
+{
+ double z = pot->z;
+ double Sz = pot->Sz;
+ double tmp = exp(d/z);
+ double tmp1 = 1./(tmp-1.);
+ return -Sz*d/(z*z)*tmp*tmp1*tmp1*(2.+d/z-2.*d/z*tmp*tmp1);
+}
+*/
+double onc_elec_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+ onc_elec_kernel* kernel = (onc_elec_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]);
+ double constant = CoulombConst;
+ double tmp = exp(-kernel->kappa*d);
+ if(d == 0.)
+ d = 1e-2;
+
+ double er = computeDielectric(d, kernel);
+ return constant*tmp/(er*d)*kernel->q;
+}
+
+void onc_elec_destroy(void* src)
+{
+ free(src);
+}
diff --git a/Kernel/tabulated_pot_kernel.c b/Kernel/tabulated_pot_kernel.c
new file mode 100644
index 0000000..3b28b99
--- /dev/null
+++ b/Kernel/tabulated_pot_kernel.c
@@ -0,0 +1,100 @@
+#include "tabulated_pot_kernel.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+struct tabulated_pot_kernel
+{
+ int num;
+ double* dist;
+ double* pot;
+};
+
+void* tabulated_pot_create(char* name)
+{
+ void* src = malloc(sizeof(tabulated_pot_kernel));
+ memset(src, 0, sizeof(tabulated_pot_kernel));
+ int capacity = 128;
+
+ ((tabulated_pot_kernel*)src)->dist = malloc(sizeof(double)*capacity);
+ ((tabulated_pot_kernel*)src)->pot = malloc(sizeof(double)*capacity);
+ 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(tmp == NULL)
+ return NULL;
+ ((tabulated_pot_kernel*)src)->dist[((tabulated_pot_kernel*)src)->num] = atof(tmp);
+ tmp = strtok(NULL, " ");
+ if(tmp == NULL)
+ return NULL;
+ ((tabulated_pot_kernel*)src)->pot[((tabulated_pot_kernel*)src)->num] = atof(tmp);
+ ((tabulated_pot_kernel*)src)->num++;
+ if(((tabulated_pot_kernel*)src)->num == capacity)
+ {
+ capacity *= 2;
+ ((tabulated_pot_kernel*)src)->dist = realloc(((tabulated_pot_kernel*)src)->dist, sizeof(double)*capacity);
+ ((tabulated_pot_kernel*)src)->pot = realloc(((tabulated_pot_kernel*)src)->pot, sizeof(double)*capacity);
+
+ }
+ }
+ ((tabulated_pot_kernel*)src)->dist = realloc(((tabulated_pot_kernel*)src)->dist, sizeof(double)*(((tabulated_pot_kernel*)src)->num));
+ ((tabulated_pot_kernel*)src)->pot = realloc(((tabulated_pot_kernel*)src)->pot, sizeof(double)*(((tabulated_pot_kernel*)src)->num));
+ fclose(infile);
+ return src;
+}
+
+static int binarySearch(double* dist, double r, int first, int last)
+{
+ int mid = (last+first)/2;
+ if(mid == 0)
+ return 0;
+ if(mid == first)
+ return mid;
+ else if(dist[mid] >= r && dist[mid-1] < r)
+ return mid;
+ else if(dist[mid] > r && dist[mid-1] > r)
+ return binarySearch(dist, r, first, mid);
+ else
+ return binarySearch(dist, r, mid, last);
+}
+
+static double interpolation(double r, tabulated_pot_kernel* tab)
+{
+ if(r <= tab->dist[0])
+ return tab->pot[0];
+ else if(r >= tab->dist[tab->num-1])
+ return tab->pot[tab->num-1];
+ int idx = binarySearch(tab->dist, r, 0, tab->num);
+ double y0 = tab->pot[(idx-1 < 0 ? 0 : idx-1)];
+ double y1 = tab->pot[idx];
+ double x0 = tab->dist[(idx-1 < 0 ? 0 : idx-1)];
+ double x1 = tab->dist[idx];
+ //printf("%lf %lf %lf\n", x0, r, x1);
+ return y0+(y1-y0)/(x1-x0)*(r-x0);
+}
+
+double tabulated_pot_compute(size_t i, size_t j, size_t k, double* origin, double** basis, void* src)
+{
+ tabulated_pot_kernel* kernel = (tabulated_pot_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 interpolation(d, kernel);
+}
+
+void tabulated_pot_destroy(void* src)
+{
+ free(((tabulated_pot_kernel*)src)->dist);
+ free(((tabulated_pot_kernel*)src)->pot);
+ free(src);
+}
--
GitLab