lenard_jones_kernel1.c

#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);
}