Add first working version of EAM

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

src/eam_utils.c

@@ -34,28 +34,24 @@
 #define MAXLINE 4096
 #endif
 
-void initEam(Eam* eam, const char* input_file, int ntypes) {
+void initEam(Eam* eam, Parameter* param) {
+    int ntypes = param->ntypes;
     eam->nmax = 0;
     eam->fp = NULL;
-    eam->ntypes = ntypes;
-    eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT));
-    coeff(eam, input_file);
-    init_style(eam);
+    coeff(eam, param);
+    init_style(eam, param);
 }
 
-void coeff(Eam* eam, const char* arg) {
-    read_file(&eam->file, arg);
-    int n = strlen(arg) + 1;
-    int ntypes = eam->ntypes;
-    double cutmax = eam->file.cut;
-    for(int i=0; i<ntypes*ntypes; i++)
-        eam->cutforcesq[i] = cutmax * cutmax;
+void coeff(Eam* eam, Parameter* param) {
+    read_file(&eam->file, param->input_file);
+    param->cutforce = eam->file.cut;
+    param->cutneigh = param->cutforce + 0.3;
 }
 
-void init_style(Eam* eam) {
+void init_style(Eam* eam, Parameter* param) {
     // convert read-in file(s) to arrays and spline them
     file2array(eam);
-    array2spline(eam);
+    array2spline(eam, param);
 }
 
 void read_file(Funcfl* file, const char* filename) {
@@ -214,7 +210,7 @@ void file2array(Eam* eam) {
     }
 }
 
-void array2spline(Eam* eam) {
+void array2spline(Eam* eam, Parameter* param) {
     eam->rdr = 1.0 / eam->dr;
     eam->rdrho = 1.0 / eam->drho;
     eam->nrho_tot = (eam->nrho + 1) * 7 + 64;
@@ -222,7 +218,7 @@ void array2spline(Eam* eam) {
     eam->nrho_tot -= eam->nrho_tot%64;
     eam->nr_tot -= eam->nr_tot%64;
 
-    int ntypes = eam->ntypes;
+    int ntypes = param->ntypes;
     eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nrho_tot * sizeof(MD_FLOAT));
     eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));
     eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * eam->nr_tot * sizeof(MD_FLOAT));

src/force_eam.c

@@ -41,10 +41,10 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
 
     int Nlocal = atom->Nlocal;
     int* neighs;
-    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; MD_FLOAT* fp = eam->fp;
+    MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; int ntypes = atom->ntypes; MD_FLOAT* fp = eam->fp;
     MD_FLOAT* rhor_spline = eam->rhor_spline; MD_FLOAT* frho_spline = eam->frho_spline; MD_FLOAT* z2r_spline = eam->z2r_spline;
     int rdr = eam->rdr; int nr = eam->nr; int nr_tot = eam->nr_tot; int rdrho = eam->rdrho;
-    int nrho = eam->nrho; int nrho_tot = eam->nrho_tot; int ntypes = eam->ntypes;
+    int nrho = eam->nrho; int nrho_tot = eam->nrho_tot;
     double S = getTimeStamp();
     LIKWID_MARKER_START("force_eam_fp");
 
@@ -67,7 +67,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
             MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
             const int type_j = atom->type[j];
             const int type_ij = type_i * ntypes + type_j;
-            const MD_FLOAT cutforcesq = eam->cutforcesq[type_ij];
+            const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
 
             if(rsq < cutforcesq) {
                 MD_FLOAT p = sqrt(rsq) * rdr + 1.0;
@@ -115,7 +115,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i
             MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
             const int type_j = atom->type[j];
             const int type_ij = type_i * ntypes + type_j;
-            const MD_FLOAT cutforcesq = eam->cutforcesq[type_ij];
+            const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij];
 
             if(rsq < cutforcesq) {
                 MD_FLOAT r = sqrt(rsq);

src/includes/eam.h

@@ -36,22 +36,20 @@ typedef struct {
 typedef struct {
     MD_FLOAT* fp;
     int nmax;
-    int ntypes;
     int nrho, nr;
     int nrho_tot, nr_tot;
     MD_FLOAT dr, rdr, drho, rdrho;
     MD_FLOAT *frho, *rhor, *z2r;
     MD_FLOAT *rhor_spline, *frho_spline, *z2r_spline;
-    MD_FLOAT *cutforcesq;
     Funcfl file;
 } Eam;
 
-void initEam(Eam* eam, const char* input_file, int ntypes);
-void coeff(Eam* eam, const char* arg);
-void init_style(Eam* eam);
+void initEam(Eam* eam, Parameter* param);
+void coeff(Eam* eam, Parameter* param);
+void init_style(Eam* eam, Parameter *param);
 void read_file(Funcfl* file, const char* filename);
 void file2array(Eam* eam);
-void array2spline(Eam* eam);
+void array2spline(Eam* eam, Parameter* param);
 void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline);
 void grab(FILE* fptr, int n, MD_FLOAT* list);
 #endif

src/main.c

@@ -84,7 +84,7 @@ double setup(
     param->zprd = param->nz * param->lattice;
 
     S = getTimeStamp();
-    if(param->force_field == FF_EAM) { initEam(eam, param->input_file, param->ntypes); }
+    if(param->force_field == FF_EAM) { initEam(eam, param); }
     initAtom(atom);
     initNeighbor(neighbor, param);
     initPbc();

src/stats.c

@@ -17,7 +17,7 @@ void displayStatistics(Atom *atom, Parameter *param, Stats *stats, double *timer
     double avg_neigh = stats->total_force_neighs / (double)(atom->Nlocal * (param->ntimes + 1));
     double avg_simd = stats->total_force_iters / (double)(atom->Nlocal * (param->ntimes + 1));
 #ifdef EXPLICIT_TYPES
-    force_useful_volume += 1e-9 * (double)((atom.Nlocal * (param.ntimes + 1)) + stats.total_force_neighs) * sizeof(int);
+    force_useful_volume += 1e-9 * (double)((atom->Nlocal * (param->ntimes + 1)) + stats->total_force_neighs) * sizeof(int);
 #endif
     printf("Statistics:\n");
     printf("\tVector width: %d, Processor frequency: %.4f GHz\n", VECTOR_WIDTH, param->proc_freq);