From 4496e911258056274ecc9e1de25970dc80cef123 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Wed, 19 May 2021 23:51:02 +0200 Subject: [PATCH] Add version with explicit types for atoms Signed-off-by: Rafael Ravedutti --- src/atom.c | 29 +++++++++++++++++++++++++++++ src/force.c | 27 +++++++++++++++++++-------- src/includes/atom.h | 8 ++++++++ src/includes/parameter.h | 1 + src/main-stub.c | 23 ++++++++++++++++++++--- src/main.c | 7 ++++--- src/neighbor.c | 12 ++++++++++-- 7 files changed, 91 insertions(+), 16 deletions(-) diff --git a/src/atom.c b/src/atom.c index 687de55..4e28d7a 100644 --- a/src/atom.c +++ b/src/atom.c @@ -41,6 +41,14 @@ 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) @@ -50,6 +58,21 @@ 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)); + atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); + atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); + for(int i = 0; i < atom->ntypes * atom->ntypes; i++) { + atom->epsilon[i] = param->epsilon; + atom->sigma6[i] = param->sigma6; + 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); int ihi = (int) (xhi / (0.5 * alat) + 1); @@ -119,6 +142,9 @@ 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++; } } @@ -151,6 +177,9 @@ 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/force.c b/src/force.c index bd471c2..d1dda58 100644 --- a/src/force.c +++ b/src/force.c @@ -27,14 +27,16 @@ #include #include -double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) { +double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor, int ntimes) { + double S = getTimeStamp(); 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; - MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz; - double S, E; + #endif for(int i = 0; i < Nlocal; i++) { fx[i] = 0.0; @@ -42,7 +44,7 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) { fz[i] = 0.0; } -#pragma omp parallel for + #pragma omp parallel for for(int i = 0; i < Nlocal; i++) { neighs = &neighbor->neighbors[i * neighbor->maxneighs]; int numneighs = neighbor->numneigh[i]; @@ -52,16 +54,24 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) { MD_FLOAT fix = 0; MD_FLOAT fiy = 0; MD_FLOAT fiz = 0; + #ifdef EXPLICIT_TYPES + const int type_i = atom->type[i]; + #endif -// printf("%d: %d\n", i, numneighs); - - for(int n = 0; n < INTERNAL_LOOP_NTIMES; n++) { + for(int n = 0; n < ntimes; n++) { 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; + #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]; + #endif if(rsq < cutforcesq) { MD_FLOAT sr2 = 1.0 / rsq; @@ -79,5 +89,6 @@ double computeForce(Parameter *param, Atom *atom, Neighbor *neighbor) { fz[i] += fiz; } - return 0.0; + double E = getTimeStamp(); + return E-S; } diff --git a/src/includes/atom.h b/src/includes/atom.h index e090936..ed1a0f9 100644 --- a/src/includes/atom.h +++ b/src/includes/atom.h @@ -30,6 +30,14 @@ 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/parameter.h b/src/includes/parameter.h index 9d51b12..0ddec95 100644 --- a/src/includes/parameter.h +++ b/src/includes/parameter.h @@ -40,6 +40,7 @@ typedef struct { MD_FLOAT temp; MD_FLOAT rho; MD_FLOAT mass; + int ntypes; int ntimes; int nstat; int every; diff --git a/src/main-stub.c b/src/main-stub.c index ad2e580..a0973f1 100644 --- a/src/main-stub.c +++ b/src/main-stub.c @@ -16,12 +16,13 @@ #define LATTICE_DISTANCE 10.0 #define NEIGH_DISTANCE 1.0 -extern double computeForce( Parameter*, Atom*, Neighbor*); +extern double computeForce( Parameter*, Atom*, Neighbor*, int); void init(Parameter *param) { param->epsilon = 1.0; param->sigma6 = 1.0; param->rho = 0.8442; + param->ntypes = 4; param->ntimes = 200; param->nx = 4; param->ny = 4; @@ -49,6 +50,7 @@ void init(Parameter *param) { atom->vx[atom->Nlocal] = vy; \ atom->vy[atom->Nlocal] = vy; \ atom->vz[atom->Nlocal] = vz; \ + atom->type[atom->Nlocal] = rand() % atom->ntypes; \ atom->Nlocal++ int main(int argc, const char *argv[]) { @@ -123,6 +125,20 @@ int main(int argc, const char *argv[]) { DEBUG("Initializing atoms...\n"); initAtom(atom); + #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)); + atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); + atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT)); + for(int i = 0; i < atom->ntypes * atom->ntypes; i++) { + atom->epsilon[i] = param->epsilon; + atom->sigma6[i] = param->sigma6; + 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) { for(int j = 0; j < param.ny; ++j) { @@ -173,6 +189,7 @@ int main(int argc, const char *argv[]) { if(!csv) { printf("Number of timesteps: %d\n", param.ntimes); + printf("Number of times to compute the most internal loop: %d\n", INTERNAL_LOOP_NTIMES); printf("System size (unit cells): %dx%dx%d\n", param.nx, param.ny, param.nz); printf("Atoms per unit cell: %d\n", atoms_per_unit_cell); printf("Total number of atoms: %d\n", atom->Nlocal); @@ -188,13 +205,13 @@ int main(int argc, const char *argv[]) { DEBUG("Building neighbor lists...\n"); buildNeighbor(atom, &neighbor); DEBUG("Computing forces...\n"); - computeForce(¶m, atom, &neighbor); + computeForce(¶m, atom, &neighbor, 1); double S, E; S = getTimeStamp(); LIKWID_MARKER_START("force"); for(int i = 0; i < param.ntimes; i++) { - computeForce(¶m, atom, &neighbor); + computeForce(¶m, atom, &neighbor, INTERNAL_LOOP_NTIMES); } LIKWID_MARKER_STOP("force"); E = getTimeStamp(); diff --git a/src/main.c b/src/main.c index 6b09c8e..061da97 100644 --- a/src/main.c +++ b/src/main.c @@ -47,13 +47,14 @@ typedef enum { NUMTIMER } timertype; -extern double computeForce( Parameter*, Atom*, Neighbor*); +extern double computeForce( Parameter*, Atom*, Neighbor*, int); void init(Parameter *param) { param->epsilon = 1.0; param->sigma6 = 1.0; param->rho = 0.8442; + param->ntypes = 4; param->ntimes = 200; param->dt = 0.005; param->nx = 32; @@ -205,7 +206,7 @@ int main (int argc, char** argv) setup(¶m, &atom, &neighbor); computeThermo(0, ¶m, &atom); - computeForce(¶m, &atom, &neighbor); + computeForce(¶m, &atom, &neighbor, 1); timer[FORCE] = 0.0; timer[NEIGH] = 0.0; @@ -221,7 +222,7 @@ int main (int argc, char** argv) timer[NEIGH] += reneighbour(¶m, &atom, &neighbor); } - timer[FORCE] += computeForce(¶m, &atom, &neighbor); + timer[FORCE] += computeForce(¶m, &atom, &neighbor, 1); finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { diff --git a/src/neighbor.c b/src/neighbor.c index 72a6ee1..91c7ab8 100644 --- a/src/neighbor.c +++ b/src/neighbor.c @@ -196,7 +196,9 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) MD_FLOAT ytmp = atom_y(i); MD_FLOAT ztmp = atom_z(i); int ibin = coord2bin(xtmp, ytmp, ztmp); - + #ifdef EXPLICIT_TYPES + int type_i = atom->type[i]; + #endif for(int k = 0; k < nstencil; k++) { int jbin = ibin + stencil[k]; int* loc_bin = &bins[jbin * atoms_per_bin]; @@ -213,7 +215,13 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) MD_FLOAT delz = ztmp - atom_z(j); MD_FLOAT rsq = delx * delx + dely * dely + delz * delz; - if( rsq <= cutneighsq ) { + #ifdef EXPLICIT_TYPES + int type_j = atom->type[j]; + const MD_FLOAT cutoff = atom->cutneighsq[type_i * atom->ntypes + type_j]; + #else + const MD_FLOAT cutoff = cutneighsq; + #endif + if( rsq <= cutoff ) { neighptr[n++] = j; } }