From 70cc6aeb1926f3c948341b6ecc57e1cad39edd8d Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 26 Oct 2021 13:55:14 +0200 Subject: [PATCH] Add first working version of EAM Signed-off-by: Rafael Ravedutti --- src/eam_utils.c | 28 ++++++++++++---------------- src/force_eam.c | 8 ++++---- src/includes/eam.h | 10 ++++------ src/main.c | 2 +- src/stats.c | 2 +- 5 files changed, 22 insertions(+), 28 deletions(-) diff --git a/src/eam_utils.c b/src/eam_utils.c index 8ffbad8..fe0a51e 100644 --- a/src/eam_utils.c +++ b/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; icutforcesq[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)); diff --git a/src/force_eam.c b/src/force_eam.c index 54987d2..322647e 100644 --- a/src/force_eam.c +++ b/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); diff --git a/src/includes/eam.h b/src/includes/eam.h index 280f419..79e46c5 100644 --- a/src/includes/eam.h +++ b/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 diff --git a/src/main.c b/src/main.c index 43391d5..362dee9 100644 --- a/src/main.c +++ b/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(); diff --git a/src/stats.c b/src/stats.c index 71f7a00..defaeac 100644 --- a/src/stats.c +++ b/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);