diff --git a/src/force_eam.c b/src/force_eam.c index c4d9fc7..b686abf 100644 --- a/src/force_eam.c +++ b/src/force_eam.c @@ -32,7 +32,7 @@ #include #include -double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { +double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats) { if(eam->nmax < atom->Nmax) { eam->nmax = atom->Nmax; if(eam->fp != NULL) { free(eam->fp); } diff --git a/src/force.c b/src/force_lj.c similarity index 92% rename from src/force.c rename to src/force_lj.c index 5ec92e7..73c0269 100644 --- a/src/force.c +++ b/src/force_lj.c @@ -26,13 +26,9 @@ #include #include #include +#include -double computeForce( - Parameter *param, - Atom *atom, - Neighbor *neighbor - ) -{ +double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) { int Nlocal = atom->Nlocal; int* neighs; MD_FLOAT* fx = atom->fx; @@ -97,10 +93,12 @@ double computeForce( 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"); double E = getTimeStamp(); - return E-S; } diff --git a/src/includes/parameter.h b/src/includes/parameter.h index b7c6017..bff46dd 100644 --- a/src/includes/parameter.h +++ b/src/includes/parameter.h @@ -23,9 +23,6 @@ #ifndef __PARAMETER_H_ #define __PARAMETER_H_ -#define FF_LJ 0 -#define FF_EAM 1 - #if PRECISION == 1 #define MD_FLOAT float #else diff --git a/src/force-tracing.c b/src/includes/tracing.h similarity index 61% rename from src/force-tracing.c rename to src/includes/tracing.h index b73305a..eab11ab 100644 --- a/src/force-tracing.c +++ b/src/includes/tracing.h @@ -20,13 +20,9 @@ * with MD-Bench. If not, see . * ======================================================================================= */ -#include - -#include #include #include #include -#include #if defined(MEM_TRACER) || defined(INDEX_TRACER) #include @@ -119,114 +115,4 @@ # define DIST_TRACE(l, e) #endif -double computeForceTracing(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { - MEM_TRACER_INIT; - INDEX_TRACER_INIT; - int Nlocal = atom->Nlocal; - int* neighs; - MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; - #ifndef EXPLICIT_TYPES - MD_FLOAT cutforcesq = param->cutforce * param->cutforce; - MD_FLOAT sigma6 = param->sigma6; - MD_FLOAT epsilon = param->epsilon; - #endif - - for(int i = 0; i < Nlocal; i++) { - fx[i] = 0.0; - fy[i] = 0.0; - fz[i] = 0.0; - } - - INDEX_TRACE_NATOMS(Nlocal, atom->Nghost, neighbor->maxneighs); - double S = getTimeStamp(); - LIKWID_MARKER_START("force"); - - for(int na = 0; na < (first_exec ? 1 : ATOMS_LOOP_RUNS); na++) { - #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 fix = 0; - MD_FLOAT fiy = 0; - MD_FLOAT fiz = 0; - - MEM_TRACE(atom_x(i), 'R'); - MEM_TRACE(atom_y(i), 'R'); - MEM_TRACE(atom_z(i), 'R'); - INDEX_TRACE_ATOM(i); - - #ifdef EXPLICIT_TYPES - const int type_i = atom->type[i]; - MEM_TRACE(atom->type(i), 'R'); - #endif - - #if defined(VARIANT) && VARIANT == stub && defined(NEIGHBORS_LOOP_RUNS) && NEIGHBORS_LOOP_RUNS > 1 - #define REPEAT_NEIGHBORS_LOOP - int nmax = first_exec ? 1 : NEIGHBORS_LOOP_RUNS; - for(int nn = 0; nn < (first_exec ? 1 : NEIGHBORS_LOOP_RUNS); nn++) { - #endif - - //DIST_TRACE_SORT(neighs, numneighs); - INDEX_TRACE(neighs, numneighs); - //DIST_TRACE(neighs, numneighs); - - 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; - - MEM_TRACE(neighs[k], 'R'); - MEM_TRACE(atom_x(j), 'R'); - MEM_TRACE(atom_y(j), 'R'); - MEM_TRACE(atom_z(j), 'R'); - - #ifdef EXPLICIT_TYPES - 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]; - const MD_FLOAT sigma6 = atom->sigma6[type_ij]; - const MD_FLOAT epsilon = atom->epsilon[type_ij]; - MEM_TRACE(atom->type(j), 'R'); - #endif - - if(rsq < cutforcesq) { - MD_FLOAT sr2 = 1.0 / rsq; - MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6; - MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon; - fix += delx * force; - fiy += dely * force; - fiz += delz * force; - } - } - - #ifdef REPEAT_NEIGHBORS_LOOP - } - #endif - - 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); - MEM_TRACE(fx[i], 'R'); - MEM_TRACE(fx[i], 'W'); - MEM_TRACE(fy[i], 'R'); - MEM_TRACE(fy[i], 'W'); - MEM_TRACE(fz[i], 'R'); - MEM_TRACE(fz[i], 'W'); - } - } - - LIKWID_MARKER_STOP("force"); - double E = getTimeStamp(); - - INDEX_TRACER_END; - MEM_TRACER_END; - return E-S; -} +extern void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timestep); diff --git a/src/includes/util.h b/src/includes/util.h index 5d6aeac..20d8806 100644 --- a/src/includes/util.h +++ b/src/includes/util.h @@ -33,5 +33,11 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#define FF_LJ 0 +#define FF_EAM 1 + extern double myrandom(int*); +extern void random_reset(int *seed, int ibase, double *coord); +extern int str2ff(const char *string); +extern const char* ff2str(int ff); #endif diff --git a/src/main-stub.c b/src/main-stub.c index 07fbbe8..4da90db 100644 --- a/src/main-stub.c +++ b/src/main-stub.c @@ -10,17 +10,21 @@ #include #include #include +#include #include #include +#include #define HLINE "----------------------------------------------------------------------------\n" #define LATTICE_DISTANCE 10.0 #define NEIGH_DISTANCE 1.0 -extern double computeForce(Parameter*, Atom*, Neighbor*, Stats*, int, int); +extern double computeForceLJ(Parameter*, Atom*, Neighbor*, Stats*); +extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*); void init(Parameter *param) { + param->input_file = NULL; param->epsilon = 1.0; param->sigma6 = 1.0; param->rho = 0.8442; @@ -39,13 +43,14 @@ void init(Parameter *param) { param->nstat = 100; param->temp = 1.44; param->every = 20; - param->proc_freq = 0.0; + param->proc_freq = 2.4; + param->eam_file = NULL; } // Show debug messages -//#define DEBUG(msg) printf(msg) +#define DEBUG(msg) printf(msg) // Do not show debug messages -#define DEBUG(msg) +//#define DEBUG(msg) #define ADD_ATOM(x, y, z, vx, vy, vz) atom_x(atom->Nlocal) = base_x + x * NEIGH_DISTANCE; \ atom_y(atom->Nlocal) = base_y + y * NEIGH_DISTANCE; \ @@ -56,6 +61,7 @@ void init(Parameter *param) { atom->Nlocal++ int main(int argc, const char *argv[]) { + Eam eam; Atom atom_data; Atom *atom = (Atom *)(&atom_data); Neighbor neighbor; @@ -71,6 +77,19 @@ int main(int argc, const 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], "-e") == 0)) + { + param.eam_file = strdup(argv[++i]); + continue; + } if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0)) { param.ntimes = atoi(argv[++i]); @@ -96,12 +115,12 @@ int main(int argc, const char *argv[]) { atoms_per_unit_cell = atoi(argv[++i]); continue; } - if((strcmp(argv[i], "-f") == 0)) + if((strcmp(argv[i], "--freq") == 0)) { param.proc_freq = atof(argv[++i]); continue; } - if((strcmp(argv[i], "-csv") == 0)) + if((strcmp(argv[i], "--csv") == 0)) { csv = 1; continue; @@ -110,16 +129,23 @@ int main(int argc, const 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("-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("-na : set number of atoms per unit cell\n"); - printf("-f : set CPU frequency (GHz) and display average cycles per atom and neighbors\n"); - printf("-csv: set output as CSV style\n"); + printf("--freq : set CPU frequency (GHz) and display average cycles per atom and neighbors\n"); + printf("--csv: set output as CSV style\n"); printf(HLINE); exit(EXIT_SUCCESS); } } + + if(param.force_field == FF_EAM) { + DEBUG("Initializing EAM parameters...\n"); + initEam(&eam, ¶m); + } + param.xprd = param.nx * LATTICE_DISTANCE; param.yprd = param.ny * LATTICE_DISTANCE; param.zprd = param.nz * LATTICE_DISTANCE; @@ -204,16 +230,29 @@ int main(int argc, const char *argv[]) { DEBUG("Initializing neighbor lists...\n"); initNeighbor(&neighbor, ¶m); DEBUG("Setting up neighbor lists...\n"); - setupNeighbor(); + setupNeighbor(¶m); DEBUG("Building neighbor lists...\n"); buildNeighbor(atom, &neighbor); DEBUG("Computing forces...\n"); - computeForce(¶m, atom, &neighbor, &stats, 1, 1); + if(param.force_field == FF_EAM) { + computeForceEam(&eam, ¶m, atom, &neighbor, &stats); + } else { + computeForceLJ(¶m, atom, &neighbor, &stats); + } double S, E; S = getTimeStamp(); for(int i = 0; i < param.ntimes; i++) { - computeForce(¶m, atom, &neighbor, &stats, 0, i + 1); + +#if defined(MEM_TRACER) || defined(INDEX_TRACER) + traceAddresses(¶m, atom, &neighbor, i + 1); +#endif + + if(param.force_field == FF_EAM) { + computeForceEam(&eam, ¶m, atom, &neighbor, &stats); + } else { + computeForceLJ(¶m, atom, &neighbor, &stats); + } } E = getTimeStamp(); double T_accum = E-S; diff --git a/src/main.c b/src/main.c index 95385de..ebb2dbd 100644 --- a/src/main.c +++ b/src/main.c @@ -41,12 +41,12 @@ #include #include #include +#include #define HLINE "----------------------------------------------------------------------------\n" -extern double computeForce(Parameter*, Atom*, Neighbor*); -extern double computeForceTracing(Parameter*, Atom*, Neighbor*, Stats*, int, int); -extern double computeForceEam(Eam* eam, Parameter*, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep); +extern double computeForceLJ(Parameter*, Atom*, Neighbor*, Stats*); +extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*); void init(Parameter *param) { @@ -166,20 +166,6 @@ void printAtomState(Atom *atom) /* } */ } -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]; @@ -266,17 +252,15 @@ int main(int argc, char** argv) setup(¶m, &eam, &atom, &neighbor, &stats); computeThermo(0, ¶m, &atom); - if(param.force_field == FF_EAM) { - computeForceEam(&eam, ¶m, &atom, &neighbor, &stats, 1, 0); - } else { -#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS) - computeForceTracing(¶m, &atom, &neighbor, &stats, 1, 0); -#else - computeForce(¶m, &atom, &neighbor); +#if defined(MEM_TRACER) || defined(INDEX_TRACER) + traceAddresses(¶m, &atom, &neighbor, n + 1); #endif + if(param.force_field == FF_EAM) { + timer[FORCE] = computeForceEam(&eam, ¶m, &atom, &neighbor, &stats); + } else { + timer[FORCE] = computeForceLJ(¶m, &atom, &neighbor, &stats); } - timer[FORCE] = 0.0; timer[NEIGH] = 0.0; timer[TOTAL] = getTimeStamp(); @@ -293,15 +277,16 @@ int main(int argc, char** argv) timer[NEIGH] += reneighbour(¶m, &atom, &neighbor); } - if(param.force_field == FF_EAM) { - timer[FORCE] += computeForceEam(&eam, ¶m, &atom, &neighbor, &stats, 0, n + 1); - } else { -#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS) - timer[FORCE] += computeForceTracing(¶m, &atom, &neighbor, &stats, 0, n + 1); -#else - timer[FORCE] += computeForce(¶m, &atom, &neighbor); +#if defined(MEM_TRACER) || defined(INDEX_TRACER) + traceAddresses(¶m, &atom, &neighbor, n + 1); #endif + + if(param.force_field == FF_EAM) { + timer[FORCE] += computeForceEam(&eam, ¶m, &atom, &neighbor, &stats); + } else { + timer[FORCE] += computeForceLJ(¶m, &atom, &neighbor, &stats); } + finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { diff --git a/src/thermo.c b/src/thermo.c index 70b897d..6f83663 100644 --- a/src/thermo.c +++ b/src/thermo.c @@ -25,6 +25,7 @@ #include #include +#include static int *steparr; static MD_FLOAT *tmparr; diff --git a/src/tracing.c b/src/tracing.c new file mode 100644 index 0000000..d91a97a --- /dev/null +++ b/src/tracing.c @@ -0,0 +1,73 @@ +/* + * ======================================================================================= + * + * 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 + +void traceAddresses(Parameter *param, Atom *atom, Neighbor *neighbor, int timestep) { + MEM_TRACER_INIT; + INDEX_TRACER_INIT; + int Nlocal = atom->Nlocal; + int* neighs; + MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; + + INDEX_TRACE_NATOMS(Nlocal, atom->Nghost, neighbor->maxneighs); + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + MEM_TRACE(atom_x(i), 'R'); + MEM_TRACE(atom_y(i), 'R'); + MEM_TRACE(atom_z(i), 'R'); + INDEX_TRACE_ATOM(i); + + #ifdef EXPLICIT_TYPES + MEM_TRACE(atom->type[i], 'R'); + #endif + + DIST_TRACE_SORT(neighs, numneighs); + INDEX_TRACE(neighs, numneighs); + DIST_TRACE(neighs, numneighs); + + for(int k = 0; k < numneighs; k++) { + MEM_TRACE(neighs[k], 'R'); + MEM_TRACE(atom_x(j), 'R'); + MEM_TRACE(atom_y(j), 'R'); + MEM_TRACE(atom_z(j), 'R'); + + #ifdef EXPLICIT_TYPES + MEM_TRACE(atom->type[j], 'R'); + #endif + } + + MEM_TRACE(fx[i], 'R'); + MEM_TRACE(fx[i], 'W'); + MEM_TRACE(fy[i], 'R'); + MEM_TRACE(fy[i], 'W'); + MEM_TRACE(fz[i], 'R'); + MEM_TRACE(fz[i], 'W'); + } + + INDEX_TRACER_END; + MEM_TRACER_END; +} diff --git a/src/util.c b/src/util.c index 1c271a3..1fa01e8 100644 --- a/src/util.c +++ b/src/util.c @@ -77,3 +77,17 @@ void random_reset(int *seed, int ibase, double *coord) for (i = 0; i < 5; i++) myrandom(seed); //save = 0; } + +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"; +}