From 9d16bb46c887df26624897ef566e0029f7e47a85 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Wed, 20 Oct 2021 22:43:08 +0200 Subject: [PATCH 1/4] Include average neighbors and SIMD iterations per atom on stats Signed-off-by: Rafael Ravedutti --- src/stats.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/stats.c b/src/stats.c index b31b35b..71f7a00 100644 --- a/src/stats.c +++ b/src/stats.c @@ -14,11 +14,15 @@ void displayStatistics(Atom *atom, Parameter *param, Stats *stats, double *timer #ifdef COMPUTE_STATS double force_useful_volume = 1e-9 * ( (double)(atom->Nlocal * (param->ntimes + 1)) * (sizeof(MD_FLOAT) * 6 + sizeof(int)) + (double)(stats->total_force_neighs) * (sizeof(MD_FLOAT) * 3 + sizeof(int)) ); + 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); #endif printf("Statistics:\n"); printf("\tVector width: %d, Processor frequency: %.4f GHz\n", VECTOR_WIDTH, param->proc_freq); + printf("\tAverage neighbors per atom: %.4f\n", avg_neigh); + printf("\tAverage SIMD iterations per atom: %.4f\n", avg_simd); printf("\tTotal number of computed pair interactions: %lld\n", stats->total_force_neighs); printf("\tTotal number of SIMD iterations: %lld\n", stats->total_force_iters); printf("\tUseful read data volume for force computation: %.2fGB\n", force_useful_volume); From 2dac10469cc56fc0303cb79ef6c763c01c7bcce5 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 26 Oct 2021 00:40:39 +0200 Subject: [PATCH 2/4] Add EAM force field Signed-off-by: Rafael Ravedutti --- src/eam_utils.c | 254 +++++++++++++++++++++++++++++++++++++++ src/force_eam.c | 173 ++++++++++++++++++++++++++ src/includes/eam.h | 54 +++++++++ src/includes/parameter.h | 4 + src/main.c | 59 ++++++++- 5 files changed, 538 insertions(+), 6 deletions(-) create mode 100644 src/eam_utils.c create mode 100644 src/force_eam.c create mode 100644 src/includes/eam.h diff --git a/src/eam_utils.c b/src/eam_utils.c new file mode 100644 index 0000000..d52626d --- /dev/null +++ b/src/eam_utils.c @@ -0,0 +1,254 @@ +#include +#include + +#include +#include +#include +#include + +void initEam(Eam* eam, const char *input_file, int ntypes) { + eam->nmax = 0; + eam->fp = NULL; + eam->ntypes = ntypes; + eam->cutforcesq = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * sizeof(MD_FLOAT)); +} + +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 init_style(Eam* eam) { + // convert read-in file(s) to arrays and spline them + file2array(eam); + array2spline(eam); +} + +void read_file(Funcfl* file, const char* filename) { + FILE* fptr; + char line[MAXLINE]; + + fptr = fopen(filename, "r"); + if(fptr == NULL) { + printf("Can't open EAM Potential file: %s\n", filename); + exit(0); + } + + int tmp; + fgets(line, MAXLINE, fptr); + fgets(line, MAXLINE, fptr); + sscanf(line, "%d %lg", &tmp, &(file->mass)); + fgets(line, MAXLINE, fptr); + sscanf(line, "%d %lg %d %lg %lg", &file->nrho, &file->drho, &file->nr, &file->dr, &file->cut); + + //printf("Read: %lf %i %lf %i %lf %lf\n",file->mass,file->nrho,file->drho,file->nr,file->dr,file->cut); + file->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nhro + 1) * sizeof(MD_FLOAT)); + file->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT)); + file->zr = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT)); + grab(fptr, file->nrho, file->frho); + grab(fptr, file->nr, file->zr); + grab(fptr, file->nr, file->rhor); + for(int i = file->nrho; i > 0; i--) file->frho[i] = file->frho[i - 1]; + for(int i = file->nr; i > 0; i--) file->rhor[i] = file->rhor[i - 1]; + for(int i = file->nr; i > 0; i--) file->zr[i] = file->zr[i - 1]; + fclose(fptr); +} + +void file2array(Eam* eam) { + int i, j, k, m, n; + double sixth = 1.0 / 6.0; + + // determine max function params from all active funcfl files + // active means some element is pointing at it via map + int active; + double rmax, rhomax; + dr = drho = rmax = rhomax = 0.0; + active = 0; + Funcfl* file = eam->file; + dr = MAX(dr, file->dr); + drho = MAX(drho, file->drho); + rmax = MAX(rmax, (file->nr - 1) * file->dr); + rhomax = MAX(rhomax, (file->nrho - 1) * file->drho); + + // set nr,nrho from cutoff and spacings + // 0.5 is for round-off in divide + eam->nr = static_cast(rmax / dr + 0.5); + eam->nrho = static_cast(rhomax / drho + 0.5); + + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ + + // allocate frho arrays + // nfrho = # of funcfl files + 1 for zero array + eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nhro + 1) * sizeof(MD_FLOAT)); + + // interpolate each file's frho to a single grid and cutoff + double r, p, cof1, cof2, cof3, cof4; + n = 0; + + for(m = 1; m <= eam->nrho; m++) { + r = (m - 1) * drho; + p = r / file->drho + 1.0; + k = static_cast(p); + k = MIN(k, file->nrho - 2); + k = MAX(k, 2); + p -= k; + p = MIN(p, 2.0); + cof1 = -sixth * p * (p - 1.0) * (p - 2.0); + cof2 = 0.5 * (p * p - 1.0) * (p - 2.0); + cof3 = -0.5 * p * (p + 1.0) * (p - 2.0); + cof4 = sixth * p * (p * p - 1.0); + eam->frho[m] = cof1 * file->frho[k - 1] + cof2 * file->frho[k] + + cof3 * file->frho[k + 1] + cof4 * file->frho[k + 2]; + } + + + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ + + // allocate rhor arrays + // nrhor = # of funcfl files + eam->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT)); + + // interpolate each file's rhor to a single grid and cutoff + for(m = 1; m <= eam->nr; m++) { + r = (m - 1) * dr; + p = r / file->dr + 1.0; + k = static_cast(p); + k = MIN(k, file->nr - 2); + k = MAX(k, 2); + p -= k; + p = MIN(p, 2.0); + cof1 = -sixth * p * (p - 1.0) * (p - 2.0); + cof2 = 0.5 * (p * p - 1.0) * (p - 2.0); + cof3 = -0.5 * p * (p + 1.0) * (p - 2.0); + cof4 = sixth * p * (p * p - 1.0); + eam->rhor[m] = cof1 * file->rhor[k - 1] + cof2 * file->rhor[k] + + cof3 * file->rhor[k + 1] + cof4 * file->rhor[k + 2]; + //if(m==119)printf("BuildRho: %e %e %e %e %e %e\n",rhor[m],cof1,cof2,cof3,cof4,file->rhor[k]); + } + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for funcfl files, I,J mapping only depends on I + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of funcfl files + eam->z2r = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nr + 1) * sizeof(MD_FLOAT)); + + // create a z2r array for each file against other files, only for I >= J + // interpolate zri and zrj to a single grid and cutoff + double zri, zrj; + Funcfl* ifile = eam->file; + Funcfl* jfile = eam->file; + + for(m = 1; m <= eam->nr; m++) { + r = (m - 1) * dr; + p = r / ifile->dr + 1.0; + k = static_cast(p); + k = MIN(k, ifile->nr - 2); + k = MAX(k, 2); + p -= k; + p = MIN(p, 2.0); + cof1 = -sixth * p * (p - 1.0) * (p - 2.0); + cof2 = 0.5 * (p * p - 1.0) * (p - 2.0); + cof3 = -0.5 * p * (p + 1.0) * (p - 2.0); + cof4 = sixth * p * (p * p - 1.0); + zri = cof1 * ifile->zr[k - 1] + cof2 * ifile->zr[k] + + cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2]; + + p = r / jfile->dr + 1.0; + k = static_cast(p); + k = MIN(k, jfile->nr - 2); + k = MAX(k, 2); + p -= k; + p = MIN(p, 2.0); + cof1 = -sixth * p * (p - 1.0) * (p - 2.0); + cof2 = 0.5 * (p * p - 1.0) * (p - 2.0); + cof3 = -0.5 * p * (p + 1.0) * (p - 2.0); + cof4 = sixth * p * (p * p - 1.0); + zrj = cof1 * jfile->zr[k - 1] + cof2 * jfile->zr[k] + + cof3 * jfile->zr[k + 1] + cof4 * jfile->zr[k + 2]; + + eam->z2r[m] = 27.2 * 0.529 * zri * zrj; + } +} + +void array2spline(Eam* eam) { + rdr = 1.0 / eam->dr; + rdrho = 1.0 / eam->drho; + nrho_tot = (eam->nrho + 1) * 7 + 64; + nr_tot = (eam->nr + 1) * 7 + 64; + nrho_tot -= nrho_tot%64; + nr_tot -= nr_tot%64; + + int ntypes = eam->ntypes; + eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nrho_tot * sizeof(MD_FLOAT)); + eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT)); + eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT)); + interpolate(eam->nrho, eam->drho, eam->frho, eam->frho_spline); + interpolate(eam->nr, eam->dr, eam->rhor, eam->rhor_spline); + interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline); + + // replicate data for multiple types; + for(int tt = 0; tt < ntypes * ntypes; tt++) { + for(int k = 0; k < nrho_tot; k++) + eam->frho_spline[tt*nrho_tot + k] = eam->frho_spline[k]; + for(int k = 0; k < nr_tot; k++) + eam->rhor_spline[tt*nr_tot + k] = eam->rhor_spline[k]; + for(int k = 0; k < nr_tot; k++) + eam->z2r_spline[tt*nr_tot + k] = eam->z2r_spline[k]; + } +} + +void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) { + for(int m = 1; m <= n; m++) spline[m * 7 + 6] = f[m]; + + spline[1 * 7 + 5] = spline[2 * 7 + 6] - spline[1 * 7 + 6]; + spline[2 * 7 + 5] = 0.5 * (spline[3 * 7 + 6] - spline[1 * 7 + 6]); + spline[(n - 1) * 7 + 5] = 0.5 * (spline[n * 7 + 6] - spline[(n - 2) * 7 + 6]); + spline[n * 7 + 5] = spline[n * 7 + 6] - spline[(n - 1) * 7 + 6]; + + for(int m = 3; m <= n - 2; m++) + spline[m * 7 + 5] = ((spline[(m - 2) * 7 + 6] - spline[(m + 2) * 7 + 6]) + + 8.0 * (spline[(m + 1) * 7 + 6] - spline[(m - 1) * 7 + 6])) / 12.0; + + for(int m = 1; m <= n - 1; m++) { + spline[m * 7 + 4] = 3.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]) - + 2.0 * spline[m * 7 + 5] - spline[(m + 1) * 7 + 5]; + spline[m * 7 + 3] = spline[m * 7 + 5] + spline[(m + 1) * 7 + 5] - + 2.0 * (spline[(m + 1) * 7 + 6] - spline[m * 7 + 6]); + } + + spline[n * 7 + 4] = 0.0; + spline[n * 7 + 3] = 0.0; + + for(int m = 1; m <= n; m++) { + spline[m * 7 + 2] = spline[m * 7 + 5] / delta; + spline[m * 7 + 1] = 2.0 * spline[m * 7 + 4] / delta; + spline[m * 7 + 0] = 3.0 * spline[m * 7 + 3] / delta; + } +} + +void grab(FILE* fptr, MMD_int n, MMD_float* list) { + char* ptr; + char line[MAXLINE]; + int i = 0; + + while(i < n) { + fgets(line, MAXLINE, fptr); + ptr = strtok(line, " \t\n\r\f"); + list[i++] = atof(ptr); + while(ptr = strtok(NULL, " \t\n\r\f")) list[i++] = atof(ptr); + } +} diff --git a/src/force_eam.c b/src/force_eam.c new file mode 100644 index 0000000..5ab58da --- /dev/null +++ b/src/force_eam.c @@ -0,0 +1,173 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2021 RRZE, University Erlangen-Nuremberg + * + * This file is part of MD-Bench. + * + * MD-Bench is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A + * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License along + * with MD-Bench. If not, see . + * ======================================================================================= + */ +#include + +#include +#include +#include +#include +#include +#include +#include + +double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { + if(atom->nmax > eam->nmax) { + eam->nmax = atom->nmax; + if(eam->fp != NULL) { free(eam->fp); } + eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->nmax * sizeof(MD_FLOAT)); + } + + 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* rhor_spline = eam->rhor_spline; MD_FLOAT* fhro_spline = eam->fhro_spline; MD_FLOAT* z2r_spline = eam->z2r_spline; + double S = getTimeStamp(); + LIKWID_MARKER_START("force_eam_fp"); + + #pragma omp parallel for + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + MD_FLOAT xtmp = atom_x(i); + MD_FLOAT ytmp = atom_y(i); + MD_FLOAT ztmp = atom_z(i); + MD_FLOAT rhoi = 0; + const int type_i = atom->type[i]; + + #pragma ivdep + for(int k = 0; k < numneighs; k++) { + int j = neighs[k]; + MD_FLOAT delx = xtmp - atom_x(j); + MD_FLOAT dely = ytmp - atom_y(j); + MD_FLOAT delz = ztmp - atom_z(j); + MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; + const int type_j = atom->type[j]; + const int type_ij = type_i * atom->ntypes + type_j; + const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij]; + + if(rsq < cutforcesq) { + MD_FLOAT p = sqrt(rsq) * rdr + 1.0; + int m = static_cast(p); + m = m < nr - 1 ? m : nr - 1; + p -= m; + p = p < 1.0 ? p : 1.0; + + rhoi += ((rhor_spline[type_ij * nr_tot + m * 7 + 3] * p + + rhor_spline[type_ij * nr_tot + m * 7 + 4]) * p + + rhor_spline[type_ij * nr_tot + m * 7 + 5]) * p + + rhor_spline[type_ij * nr_tot + m * 7 + 6]; + } + } + + const int type_ii = type_i * type_i; + MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0; + int m = static_cast(p); + m = MAX(1, MIN(m, nrho - 1)); + p -= m; + p = MIN(p, 1.0); + fp[i] = (frho_spline[type_ii * nrho_tot + m * 7 + 0] * p + + frho_spline[type_ii * nrho_tot + m * 7 + 1]) * p + + frho_spline[type_ii * nrho_tot + m * 7 + 2]; + } + + LIKWID_MARKER_STOP("force_eam_fp"); + LIKWID_MARKER_START("force_eam"); + + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + MD_FLOAT xtmp = atom_x(i); + MD_FLOAT ytmp = atom_y(i); + MD_FLOAT ztmp = atom_z(i); + MD_FLOAT fix = 0; + MD_FLOAT fiy = 0; + MD_FLOAT fiz = 0; + const int type_i = atom->type[i]; + + #pragma ivdep + for(int k = 0; k < numneighs; k++) { + int j = neighs[k]; + MD_FLOAT delx = xtmp - atom_x(j); + MD_FLOAT dely = ytmp - atom_y(j); + MD_FLOAT delz = ztmp - atom_z(j); + MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; + const int type_j = atom->type[j]; + const int type_ij = type_i * atom->ntypes + type_j; + const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij]; + + if(rsq < cutforcesq) { + MMD_float r = sqrt(rsq); + MMD_float p = r * rdr + 1.0; + int m = static_cast(p); + m = m < nr - 1 ? m : nr - 1; + p -= m; + p = p < 1.0 ? p : 1.0; + + + // rhoip = derivative of (density at atom j due to atom i) + // rhojp = derivative of (density at atom i due to atom j) + // phi = pair potential energy + // phip = phi' + // z2 = phi * r + // z2p = (phi * r)' = (phi' r) + phi + // psip needs both fp[i] and fp[j] terms since r_ij appears in two + // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) + // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip + + MMD_float rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p + + rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p + + rhor_spline[type_ij * nr_tot + m * 7 + 2]; + + MMD_float z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p + + z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 2]; + + MMD_float z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p + + z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 6]; + + MMD_float recip = 1.0 / r; + MMD_float phi = z2 * recip; + MMD_float phip = z2p * recip - phi * recip; + MMD_float psip = fp[i] * rhoip + fp[j] * rhoip + phip; + MMD_float fpair = -psip * recip; + + fix += delx * fpair; + fiy += dely * fpair; + fiz += delz * fpair; + //fpair *= 0.5; + } + } + + fx[i] = fix; + fy[i] = fiy; + fz[i] = fiz; + addStat(stats->total_force_neighs, numneighs); + addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH); + } + + LIKWID_MARKER_STOP("force_eam"); + double E = getTimeStamp(); + return E-S; +} diff --git a/src/includes/eam.h b/src/includes/eam.h new file mode 100644 index 0000000..be23d41 --- /dev/null +++ b/src/includes/eam.h @@ -0,0 +1,54 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * This file is part of MD-Bench. + * + * MD-Bench is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A + * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License along + * with MD-Bench. If not, see . + * ======================================================================================= + */ +#include +#include + +#ifndef __EAM_H_ +#define __EAM_H_ +typedef struct { + int nrho, nr; + MD_FLOAT drho, dr, cut, mass; + MD_FLOAT *frho, *rhor, *zr; +} Funcfl; + +typedef struct { + MD_FLOAT* fp; + int nmax; + 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 init_eam(Eam* eam, int ntypes); +void coeff(Eam* eam, const char* arg); +void init_style(Eam* eam); +void read_file(Funcfl* file, const char* filename); +void file2array(Eam* eam); +void array2spline(Eam* eam); +void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline); +void grab(FILE* fptr, MMD_int n, MMD_float* list); +#endif diff --git a/src/includes/parameter.h b/src/includes/parameter.h index e95e1e8..e8eef2b 100644 --- a/src/includes/parameter.h +++ b/src/includes/parameter.h @@ -23,6 +23,9 @@ #ifndef __PARAMETER_H_ #define __PARAMETER_H_ +#define FF_LJ 0 +#define FF_EAM 1 + #if PRECISION == 1 #define MD_FLOAT float #else @@ -30,6 +33,7 @@ #endif typedef struct { + char* input_file; MD_FLOAT epsilon; MD_FLOAT sigma6; MD_FLOAT temp; diff --git a/src/main.c b/src/main.c index bff4422..00ae044 100644 --- a/src/main.c +++ b/src/main.c @@ -39,13 +39,17 @@ #include #include #include +#include #define HLINE "----------------------------------------------------------------------------\n" extern double computeForce(Parameter*, Atom*, Neighbor*, Stats*, int, int); +extern double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep); void init(Parameter *param) { + param->input_file = NULL; + param->force_field = FF_LJ; param->epsilon = 1.0; param->sigma6 = 1.0; param->rho = 0.8442; @@ -67,6 +71,7 @@ void init(Parameter *param) double setup( Parameter *param, + Eam *eam, Atom *atom, Neighbor *neighbor, Stats *stats) @@ -78,6 +83,7 @@ double setup( param->zprd = param->nz * param->lattice; S = getTimeStamp(); + if(param->force_field == FF_EAM) { initEam(eam, param->input_file, param->ntypes); } initAtom(atom); initNeighbor(neighbor, param); initPbc(); @@ -153,9 +159,24 @@ void printAtomState(Atom *atom) /* } */ } -int main (int argc, char** argv) +int str2ff(const char *string) +{ + if(strncmp(string, "lj", 2) == 0) return FF_LJ; + if(strncmp(string, "eam", 3) == 0) return FF_EAM; + return -1; +} + +const char* ff2str(int ff) +{ + if(ff == FF_LJ) { return "lj"; } + if(ff == FF_EAM) { return "eam"; } + return "invalid"; +} + +int main(int argc, char** argv) { double timer[NUMTIMER]; + Eam eam; Atom atom; Neighbor neighbor; Stats stats; @@ -172,6 +193,19 @@ int main (int argc, char** argv) for(int i = 0; i < argc; i++) { + if((strcmp(argv[i], "-f") == 0)) + { + if((param.force_field = str2ff(argv[++i])) < 0) { + fprintf(stderr, "Invalid force field!\n"); + exit(-1); + } + continue; + } + if((strcmp(argv[i], "-i") == 0)) + { + param.input_file = strdup(argv[++i]); + continue; + } if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0)) { param.ntimes = atoi(argv[++i]); @@ -192,7 +226,7 @@ int main (int argc, char** argv) param.nz = atoi(argv[++i]); continue; } - if((strcmp(argv[i], "-f") == 0)) + if((strcmp(argv[i], "--freq") == 0)) { param.proc_freq = atof(argv[++i]); continue; @@ -201,17 +235,24 @@ int main (int argc, char** argv) { printf("MD Bench: A minimalistic re-implementation of miniMD\n"); printf(HLINE); + printf("-f : force field (lj or eam), default lj\n"); + printf("-i : input file for EAM\n"); printf("-n / --nsteps : set number of timesteps for simulation\n"); printf("-nx/-ny/-nz : set linear dimension of systembox in x/y/z direction\n"); - printf("-f : processor frequency (GHz)\n"); + printf("--freq : processor frequency (GHz)\n"); printf(HLINE); exit(EXIT_SUCCESS); } } - setup(¶m, &atom, &neighbor, &stats); + setup(¶m, &eam, &atom, &neighbor, &stats); computeThermo(0, ¶m, &atom); - computeForce(¶m, &atom, &neighbor, &stats, 1, 0); + + if(param.force_field == FF_EAM) { + computeForceEam(&eam, &atom, &neighbor, &stats, 1, 0); + } else { + computeForce(¶m, &atom, &neighbor, &stats, 1, 0); + } timer[FORCE] = 0.0; timer[NEIGH] = 0.0; @@ -226,7 +267,12 @@ int main (int argc, char** argv) timer[NEIGH] += reneighbour(¶m, &atom, &neighbor); } - timer[FORCE] += computeForce(¶m, &atom, &neighbor, &stats, 0, n + 1); + if(param.force_field == FF_EAM) { + timer[FORCE] += computeForceEam(&eam, &atom, &neighbor, &stats, 0, n + 1); + } else { + timer[FORCE] += computeForce(¶m, &atom, &neighbor, &stats, 0, n + 1); + } + finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { @@ -238,6 +284,7 @@ int main (int argc, char** argv) computeThermo(-1, ¶m, &atom); printf(HLINE); + printf("Force field: %s\n", ff2str(param.force_field)); printf("Data layout for positions: %s\n", POS_DATA_LAYOUT); #if PRECISION == 1 printf("Using single precision floating point.\n"); From 40ddc9ad505e29a352d551fd0450485623be6748 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 26 Oct 2021 01:19:11 +0200 Subject: [PATCH 3/4] Fix errors introduced by last changes Signed-off-by: Rafael Ravedutti --- src/atom.c | 9 ------ src/eam_utils.c | 69 ++++++++++++++++++++++------------------ src/force_eam.c | 38 ++++++++++++---------- src/includes/atom.h | 2 -- src/includes/eam.h | 9 ++++-- src/includes/parameter.h | 1 + src/main-stub.c | 4 --- 7 files changed, 66 insertions(+), 66 deletions(-) diff --git a/src/atom.c b/src/atom.c index 4e28d7a..9964a14 100644 --- a/src/atom.c +++ b/src/atom.c @@ -41,14 +41,12 @@ void initAtom(Atom *atom) atom->Nlocal = 0; atom->Nghost = 0; atom->Nmax = 0; - #ifdef EXPLICIT_TYPES atom->type = NULL; atom->ntypes = 0; atom->epsilon = NULL; atom->sigma6 = NULL; atom->cutforcesq = NULL; atom->cutneighsq = NULL; - #endif } void createAtom(Atom *atom, Parameter *param) @@ -58,8 +56,6 @@ void createAtom(Atom *atom, Parameter *param) MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd; atom->Natoms = 4 * param->nx * param->ny * param->nz; atom->Nlocal = 0; - - #ifdef EXPLICIT_TYPES atom->ntypes = param->ntypes; atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); @@ -71,7 +67,6 @@ void createAtom(Atom *atom, Parameter *param) atom->cutneighsq[i] = param->cutneigh * param->cutneigh; atom->cutforcesq[i] = param->cutforce * param->cutforce; } - #endif MD_FLOAT alat = pow((4.0 / param->rho), (1.0 / 3.0)); int ilo = (int) (xlo / (0.5 * alat) - 1); @@ -142,9 +137,7 @@ void createAtom(Atom *atom, Parameter *param) atom->vx[atom->Nlocal] = vxtmp; atom->vy[atom->Nlocal] = vytmp; atom->vz[atom->Nlocal] = vztmp; - #ifdef EXPLICIT_TYPES atom->type[atom->Nlocal] = rand() % atom->ntypes; - #endif atom->Nlocal++; } } @@ -177,9 +170,7 @@ void growAtom(Atom *atom) atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); - #ifdef EXPLICIT_TYPES atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); - #endif } diff --git a/src/eam_utils.c b/src/eam_utils.c index d52626d..7d2ef77 100644 --- a/src/eam_utils.c +++ b/src/eam_utils.c @@ -1,16 +1,23 @@ #include #include +#include #include #include #include #include +#include + +#ifndef MAXLINE +#define MAXLINE 4096 +#endif void initEam(Eam* eam, const char *input_file, int 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); } void coeff(Eam* eam, const char* arg) { @@ -46,7 +53,7 @@ void read_file(Funcfl* file, const char* filename) { sscanf(line, "%d %lg %d %lg %lg", &file->nrho, &file->drho, &file->nr, &file->dr, &file->cut); //printf("Read: %lf %i %lf %i %lf %lf\n",file->mass,file->nrho,file->drho,file->nr,file->dr,file->cut); - file->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nhro + 1) * sizeof(MD_FLOAT)); + file->frho = (MD_FLOAT *) allocate(ALIGNMENT, (file->nrho + 1) * sizeof(MD_FLOAT)); file->rhor = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT)); file->zr = (MD_FLOAT *) allocate(ALIGNMENT, (file->nr + 1) * sizeof(MD_FLOAT)); grab(fptr, file->nrho, file->frho); @@ -66,18 +73,18 @@ void file2array(Eam* eam) { // active means some element is pointing at it via map int active; double rmax, rhomax; - dr = drho = rmax = rhomax = 0.0; + eam->dr = eam->drho = rmax = rhomax = 0.0; active = 0; Funcfl* file = eam->file; - dr = MAX(dr, file->dr); - drho = MAX(drho, file->drho); + eam->dr = MAX(eam->dr, file->dr); + eam->drho = MAX(eam->drho, file->drho); rmax = MAX(rmax, (file->nr - 1) * file->dr); rhomax = MAX(rhomax, (file->nrho - 1) * file->drho); // set nr,nrho from cutoff and spacings // 0.5 is for round-off in divide - eam->nr = static_cast(rmax / dr + 0.5); - eam->nrho = static_cast(rhomax / drho + 0.5); + eam->nr = (int)(rmax / eam->dr + 0.5); + eam->nrho = (int)(rhomax / eam->drho + 0.5); // ------------------------------------------------------------------ // setup frho arrays @@ -85,16 +92,16 @@ void file2array(Eam* eam) { // allocate frho arrays // nfrho = # of funcfl files + 1 for zero array - eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nhro + 1) * sizeof(MD_FLOAT)); + eam->frho = (MD_FLOAT *) allocate(ALIGNMENT, (eam->nrho + 1) * sizeof(MD_FLOAT)); // interpolate each file's frho to a single grid and cutoff double r, p, cof1, cof2, cof3, cof4; n = 0; for(m = 1; m <= eam->nrho; m++) { - r = (m - 1) * drho; + r = (m - 1) * eam->drho; p = r / file->drho + 1.0; - k = static_cast(p); + k = (int)(p); k = MIN(k, file->nrho - 2); k = MAX(k, 2); p -= k; @@ -118,9 +125,9 @@ void file2array(Eam* eam) { // interpolate each file's rhor to a single grid and cutoff for(m = 1; m <= eam->nr; m++) { - r = (m - 1) * dr; + r = (m - 1) * eam->dr; p = r / file->dr + 1.0; - k = static_cast(p); + k = (int)(p); k = MIN(k, file->nr - 2); k = MAX(k, 2); p -= k; @@ -153,9 +160,9 @@ void file2array(Eam* eam) { Funcfl* jfile = eam->file; for(m = 1; m <= eam->nr; m++) { - r = (m - 1) * dr; + r = (m - 1) * eam->dr; p = r / ifile->dr + 1.0; - k = static_cast(p); + k = (int)(p); k = MIN(k, ifile->nr - 2); k = MAX(k, 2); p -= k; @@ -168,7 +175,7 @@ void file2array(Eam* eam) { cof3 * ifile->zr[k + 1] + cof4 * ifile->zr[k + 2]; p = r / jfile->dr + 1.0; - k = static_cast(p); + k = (int)(p); k = MIN(k, jfile->nr - 2); k = MAX(k, 2); p -= k; @@ -185,33 +192,33 @@ void file2array(Eam* eam) { } void array2spline(Eam* eam) { - rdr = 1.0 / eam->dr; - rdrho = 1.0 / eam->drho; - nrho_tot = (eam->nrho + 1) * 7 + 64; - nr_tot = (eam->nr + 1) * 7 + 64; - nrho_tot -= nrho_tot%64; - nr_tot -= nr_tot%64; + eam->rdr = 1.0 / eam->dr; + eam->rdrho = 1.0 / eam->drho; + eam->nrho_tot = (eam->nrho + 1) * 7 + 64; + eam->nr_tot = (eam->nr + 1) * 7 + 64; + eam->nrho_tot -= eam->nrho_tot%64; + eam->nr_tot -= eam->nr_tot%64; int ntypes = eam->ntypes; - eam->frho_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nrho_tot * sizeof(MD_FLOAT)); - eam->rhor_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT)); - eam->z2r_spline = (MD_FLOAT *) allocate(ALIGNMENT, ntypes * ntypes * nr_tot * sizeof(MD_FLOAT)); + 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)); interpolate(eam->nrho, eam->drho, eam->frho, eam->frho_spline); interpolate(eam->nr, eam->dr, eam->rhor, eam->rhor_spline); interpolate(eam->nr, eam->dr, eam->z2r, eam->z2r_spline); // replicate data for multiple types; for(int tt = 0; tt < ntypes * ntypes; tt++) { - for(int k = 0; k < nrho_tot; k++) - eam->frho_spline[tt*nrho_tot + k] = eam->frho_spline[k]; - for(int k = 0; k < nr_tot; k++) - eam->rhor_spline[tt*nr_tot + k] = eam->rhor_spline[k]; - for(int k = 0; k < nr_tot; k++) - eam->z2r_spline[tt*nr_tot + k] = eam->z2r_spline[k]; + for(int k = 0; k < eam->nrho_tot; k++) + eam->frho_spline[tt*eam->nrho_tot + k] = eam->frho_spline[k]; + for(int k = 0; k < eam->nr_tot; k++) + eam->rhor_spline[tt*eam->nr_tot + k] = eam->rhor_spline[k]; + for(int k = 0; k < eam->nr_tot; k++) + eam->z2r_spline[tt*eam->nr_tot + k] = eam->z2r_spline[k]; } } -void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) { +void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline) { for(int m = 1; m <= n; m++) spline[m * 7 + 6] = f[m]; spline[1 * 7 + 5] = spline[2 * 7 + 6] - spline[1 * 7 + 6]; @@ -240,7 +247,7 @@ void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline) { } } -void grab(FILE* fptr, MMD_int n, MMD_float* list) { +void grab(FILE* fptr, int n, MD_FLOAT* list) { char* ptr; char line[MAXLINE]; int i = 0; diff --git a/src/force_eam.c b/src/force_eam.c index 5ab58da..1e5bb72 100644 --- a/src/force_eam.c +++ b/src/force_eam.c @@ -21,6 +21,7 @@ * ======================================================================================= */ #include +#include #include #include @@ -29,18 +30,21 @@ #include #include #include +#include double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { - if(atom->nmax > eam->nmax) { - eam->nmax = atom->nmax; + if(eam->nmax < atom->Nmax) { + eam->nmax = atom->Nmax; if(eam->fp != NULL) { free(eam->fp); } - eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->nmax * sizeof(MD_FLOAT)); + eam->fp = (MD_FLOAT *) allocate(ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT)); } 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* rhor_spline = eam->rhor_spline; MD_FLOAT* fhro_spline = eam->fhro_spline; MD_FLOAT* z2r_spline = eam->z2r_spline; + 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; double S = getTimeStamp(); LIKWID_MARKER_START("force_eam_fp"); @@ -67,7 +71,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i if(rsq < cutforcesq) { MD_FLOAT p = sqrt(rsq) * rdr + 1.0; - int m = static_cast(p); + int m = (int)(p); m = m < nr - 1 ? m : nr - 1; p -= m; p = p < 1.0 ? p : 1.0; @@ -81,7 +85,7 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i const int type_ii = type_i * type_i; MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0; - int m = static_cast(p); + int m = (int)(p); m = MAX(1, MIN(m, nrho - 1)); p -= m; p = MIN(p, 1.0); @@ -116,9 +120,9 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij]; if(rsq < cutforcesq) { - MMD_float r = sqrt(rsq); - MMD_float p = r * rdr + 1.0; - int m = static_cast(p); + MD_FLOAT r = sqrt(rsq); + MD_FLOAT p = r * rdr + 1.0; + int m = (int)(p); m = m < nr - 1 ? m : nr - 1; p -= m; p = p < 1.0 ? p : 1.0; @@ -134,24 +138,24 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip - MMD_float rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p + + MD_FLOAT rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p + rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p + rhor_spline[type_ij * nr_tot + m * 7 + 2]; - MMD_float z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p + + MD_FLOAT z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p + z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p + z2r_spline[type_ij * nr_tot + m * 7 + 2]; - MMD_float z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p + + MD_FLOAT z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p + z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p + z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p + z2r_spline[type_ij * nr_tot + m * 7 + 6]; - MMD_float recip = 1.0 / r; - MMD_float phi = z2 * recip; - MMD_float phip = z2p * recip - phi * recip; - MMD_float psip = fp[i] * rhoip + fp[j] * rhoip + phip; - MMD_float fpair = -psip * recip; + MD_FLOAT recip = 1.0 / r; + MD_FLOAT phi = z2 * recip; + MD_FLOAT phip = z2p * recip - phi * recip; + MD_FLOAT psip = fp[i] * rhoip + fp[j] * rhoip + phip; + MD_FLOAT fpair = -psip * recip; fix += delx * fpair; fiy += dely * fpair; diff --git a/src/includes/atom.h b/src/includes/atom.h index ed1a0f9..98d2117 100644 --- a/src/includes/atom.h +++ b/src/includes/atom.h @@ -30,14 +30,12 @@ typedef struct { MD_FLOAT *x, *y, *z; MD_FLOAT *vx, *vy, *vz; MD_FLOAT *fx, *fy, *fz; - #ifdef EXPLICIT_TYPES int *type; int ntypes; MD_FLOAT *epsilon; MD_FLOAT *sigma6; MD_FLOAT *cutforcesq; MD_FLOAT *cutneighsq; - #endif } Atom; extern void initAtom(Atom*); diff --git a/src/includes/eam.h b/src/includes/eam.h index be23d41..0ac250d 100644 --- a/src/includes/eam.h +++ b/src/includes/eam.h @@ -20,6 +20,8 @@ * with MD-Bench. If not, see . * ======================================================================================= */ +#include + #include #include @@ -34,6 +36,7 @@ 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; @@ -43,12 +46,12 @@ typedef struct { Funcfl* file; } Eam; -void init_eam(Eam* eam, int ntypes); +void initEam(Eam* eam, const char* input_file, int ntypes); void coeff(Eam* eam, const char* arg); void init_style(Eam* eam); void read_file(Funcfl* file, const char* filename); void file2array(Eam* eam); void array2spline(Eam* eam); -void interpolate(MMD_int n, MMD_float delta, MMD_float* f, MMD_float* spline); -void grab(FILE* fptr, MMD_int n, MMD_float* list); +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/includes/parameter.h b/src/includes/parameter.h index e8eef2b..ac6c737 100644 --- a/src/includes/parameter.h +++ b/src/includes/parameter.h @@ -33,6 +33,7 @@ #endif typedef struct { + int force_field; char* input_file; MD_FLOAT epsilon; MD_FLOAT sigma6; diff --git a/src/main-stub.c b/src/main-stub.c index 72521f2..07fbbe8 100644 --- a/src/main-stub.c +++ b/src/main-stub.c @@ -128,7 +128,6 @@ int main(int argc, const char *argv[]) { initAtom(atom); initStats(&stats); - #ifdef EXPLICIT_TYPES atom->ntypes = param.ntypes; atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); @@ -140,7 +139,6 @@ int main(int argc, const char *argv[]) { atom->cutneighsq[i] = param.cutneigh * param.cutneigh; atom->cutforcesq[i] = param.cutforce * param.cutforce; } - #endif DEBUG("Creating atoms...\n"); for(int i = 0; i < param.nx; ++i) { @@ -176,9 +174,7 @@ int main(int argc, const char *argv[]) { for(int jj = 0; jj < fac_y; ++jj) { for(int kk = 0; kk < fac_z; ++kk) { if(added_atoms < atoms_per_unit_cell) { - #ifdef EXPLICIT_TYPES atom->type[atom->Nlocal] = rand() % atom->ntypes; - #endif ADD_ATOM(ii * offset_x, jj * offset_y, kk * offset_z, vx, vy, vz); added_atoms++; } From 99d6a4bdd851a0475238d7a931a9242d48617b9c Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Tue, 26 Oct 2021 01:40:02 +0200 Subject: [PATCH 4/4] Fix Funcfl reference to avoid segfaults Signed-off-by: Rafael Ravedutti --- src/eam_utils.c | 35 +++++++++++++++++++++++++++++------ src/force_eam.c | 14 +++++++------- src/includes/eam.h | 2 +- 3 files changed, 37 insertions(+), 14 deletions(-) diff --git a/src/eam_utils.c b/src/eam_utils.c index 7d2ef77..8ffbad8 100644 --- a/src/eam_utils.c +++ b/src/eam_utils.c @@ -1,3 +1,25 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * This file is part of MD-Bench. + * + * MD-Bench is free software: you can redistribute it and/or modify it + * under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A + * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License along + * with MD-Bench. If not, see . + * ======================================================================================= + */ #include #include #include @@ -12,19 +34,20 @@ #define MAXLINE 4096 #endif -void initEam(Eam* eam, const char *input_file, int ntypes) { +void initEam(Eam* eam, const char* input_file, int 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); } void coeff(Eam* eam, const char* arg) { - read_file(eam->file, arg); + read_file(&eam->file, arg); int n = strlen(arg) + 1; int ntypes = eam->ntypes; - double cutmax = eam->file->cut; + double cutmax = eam->file.cut; for(int i=0; icutforcesq[i] = cutmax * cutmax; } @@ -75,7 +98,7 @@ void file2array(Eam* eam) { double rmax, rhomax; eam->dr = eam->drho = rmax = rhomax = 0.0; active = 0; - Funcfl* file = eam->file; + Funcfl* file = &eam->file; eam->dr = MAX(eam->dr, file->dr); eam->drho = MAX(eam->drho, file->drho); rmax = MAX(rmax, (file->nr - 1) * file->dr); @@ -156,8 +179,8 @@ void file2array(Eam* eam) { // create a z2r array for each file against other files, only for I >= J // interpolate zri and zrj to a single grid and cutoff double zri, zrj; - Funcfl* ifile = eam->file; - Funcfl* jfile = eam->file; + Funcfl* ifile = &eam->file; + Funcfl* jfile = &eam->file; for(m = 1; m <= eam->nr; m++) { r = (m - 1) * eam->dr; diff --git a/src/force_eam.c b/src/force_eam.c index 1e5bb72..ce32d91 100644 --- a/src/force_eam.c +++ b/src/force_eam.c @@ -139,17 +139,17 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip MD_FLOAT rhoip = (rhor_spline[type_ij * nr_tot + m * 7 + 0] * p + - rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p + - rhor_spline[type_ij * nr_tot + m * 7 + 2]; + rhor_spline[type_ij * nr_tot + m * 7 + 1]) * p + + rhor_spline[type_ij * nr_tot + m * 7 + 2]; MD_FLOAT z2p = (z2r_spline[type_ij * nr_tot + m * 7 + 0] * p + - z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p + - z2r_spline[type_ij * nr_tot + m * 7 + 2]; + z2r_spline[type_ij * nr_tot + m * 7 + 1]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 2]; MD_FLOAT z2 = ((z2r_spline[type_ij * nr_tot + m * 7 + 3] * p + - z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p + - z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p + - z2r_spline[type_ij * nr_tot + m * 7 + 6]; + z2r_spline[type_ij * nr_tot + m * 7 + 4]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 5]) * p + + z2r_spline[type_ij * nr_tot + m * 7 + 6]; MD_FLOAT recip = 1.0 / r; MD_FLOAT phi = z2 * recip; diff --git a/src/includes/eam.h b/src/includes/eam.h index 0ac250d..280f419 100644 --- a/src/includes/eam.h +++ b/src/includes/eam.h @@ -43,7 +43,7 @@ typedef struct { MD_FLOAT *frho, *rhor, *z2r; MD_FLOAT *rhor_spline, *frho_spline, *z2r_spline; MD_FLOAT *cutforcesq; - Funcfl* file; + Funcfl file; } Eam; void initEam(Eam* eam, const char* input_file, int ntypes);