Skip to content
Snippets Groups Projects
Commit 53917a12 authored by hchou10's avatar hchou10
Browse files

forgot to push some files

parent 141591ad
No related branches found
No related tags found
No related merge requests found
Pipeline #168412 canceled
#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);
}
#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);
}
#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);
}
#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);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment