diff --git a/src/atom.c b/src/atom.c index 9964a14..1ed54e9 100644 --- a/src/atom.c +++ b/src/atom.c @@ -172,49 +172,3 @@ void growAtom(Atom *atom) atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT)); atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int)); } - - -/* void sortAtom() */ -/* { */ -/* binatoms(neighbor); */ -/* int* binpos = neighbor->bincount; */ -/* int* bins = neighbor->bins; */ - -/* int mbins = neighbor->mbins; */ -/* int atoms_per_bin = neighbor->atoms_per_bin; */ - -/* for(int i=1; i0?binpos[mybin-1]:0; */ -/* int count = binpos[mybin] - start; */ - -/* for(int k=0; kfrho[m] = cof1 * file->frho[k - 1] + cof2 * file->frho[k] + - cof3 * file->frho[k + 1] + cof4 * file->frho[k + 2]; + cof3 * file->frho[k + 1] + cof4 * file->frho[k + 2]; } @@ -161,7 +161,7 @@ void file2array(Eam* eam) { 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]; + 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]); } diff --git a/src/force_eam.c b/src/force_eam.c index 322647e..a284e2e 100644 --- a/src/force_eam.c +++ b/src/force_eam.c @@ -32,7 +32,7 @@ #include #include -double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { +double computeForceEam(Eam* eam, Parameter* param, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep) { if(eam->nmax < atom->Nmax) { eam->nmax = atom->Nmax; if(eam->fp != NULL) { free(eam->fp); } @@ -46,8 +46,8 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i 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"); + LIKWID_MARKER_START("force_eam_fp"); #pragma omp parallel for for(int i = 0; i < Nlocal; i++) { neighs = &neighbor->neighbors[i * neighbor->maxneighs]; @@ -56,8 +56,9 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i MD_FLOAT ytmp = atom_y(i); MD_FLOAT ztmp = atom_z(i); MD_FLOAT rhoi = 0; +#ifdef EXPLICIT_TYPES const int type_i = atom->type[i]; - +#endif #pragma ivdep for(int k = 0; k < numneighs; k++) { int j = neighs[k]; @@ -65,36 +66,58 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i 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 * ntypes + type_j; const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij]; - +#else + const MD_FLOAT cutforcesq = param->cutforce * param->cutforce; +#endif if(rsq < cutforcesq) { MD_FLOAT p = sqrt(rsq) * rdr + 1.0; int m = (int)(p); m = m < nr - 1 ? m : nr - 1; p -= m; p = p < 1.0 ? p : 1.0; +#ifdef EXPLICIT_TYPES 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]; +#else + rhoi += ((rhor_spline[m * 7 + 3] * p + + rhor_spline[m * 7 + 4]) * p + + rhor_spline[m * 7 + 5]) * p + + rhor_spline[m * 7 + 6]; +#endif } } +#ifdef EXPLICIT_TYPES const int type_ii = type_i * type_i; +#endif MD_FLOAT p = 1.0 * rhoi * rdrho + 1.0; int m = (int)(p); m = MAX(1, MIN(m, nrho - 1)); p -= m; p = MIN(p, 1.0); +#ifdef EXPLICIT_TYPES 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]; +#else + fp[i] = (frho_spline[m * 7 + 0] * p + frho_spline[m * 7 + 1]) * p + frho_spline[m * 7 + 2]; +#endif } - LIKWID_MARKER_STOP("force_eam_fp"); - LIKWID_MARKER_START("force_eam"); + LIKWID_MARKER_STOP("force_eam_fp"); + + // We still need to update fp for PBC atoms + for(int i = 0; i < atom->Nghost; i++) { + fp[Nlocal + i] = fp[atom->border_map[i]]; + } + + LIKWID_MARKER_START("force_eam"); for(int i = 0; i < Nlocal; i++) { neighs = &neighbor->neighbors[i * neighbor->maxneighs]; int numneighs = neighbor->numneigh[i]; @@ -104,7 +127,9 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i MD_FLOAT fix = 0; MD_FLOAT fiy = 0; MD_FLOAT fiz = 0; +#ifdef EXPLICIT_TYPES const int type_i = atom->type[i]; +#endif #pragma ivdep for(int k = 0; k < numneighs; k++) { @@ -113,9 +138,13 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i 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 * ntypes + type_j; const MD_FLOAT cutforcesq = atom->cutforcesq[type_ij]; +#else + const MD_FLOAT cutforcesq = param->cutforce * param->cutforce; +#endif if(rsq < cutforcesq) { MD_FLOAT r = sqrt(rsq); @@ -136,6 +165,7 @@ 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 +#ifdef EXPLICIT_TYPES 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]; @@ -148,6 +178,14 @@ double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, i 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]; +#else + MD_FLOAT rhoip = (rhor_spline[m * 7 + 0] * p + rhor_spline[m * 7 + 1]) * p + rhor_spline[m * 7 + 2]; + MD_FLOAT z2p = (z2r_spline[m * 7 + 0] * p + z2r_spline[m * 7 + 1]) * p + z2r_spline[m * 7 + 2]; + MD_FLOAT z2 = ((z2r_spline[m * 7 + 3] * p + + z2r_spline[m * 7 + 4]) * p + + z2r_spline[m * 7 + 5]) * p + + z2r_spline[m * 7 + 6]; +#endif MD_FLOAT recip = 1.0 / r; MD_FLOAT phi = z2 * recip; diff --git a/src/includes/atom.h b/src/includes/atom.h index 98d2117..89c01ea 100644 --- a/src/includes/atom.h +++ b/src/includes/atom.h @@ -30,6 +30,7 @@ typedef struct { MD_FLOAT *x, *y, *z; MD_FLOAT *vx, *vy, *vz; MD_FLOAT *fx, *fy, *fz; + int *border_map; int *type; int ntypes; MD_FLOAT *epsilon; diff --git a/src/includes/neighbor.h b/src/includes/neighbor.h index c9bad95..8ea52f0 100644 --- a/src/includes/neighbor.h +++ b/src/includes/neighbor.h @@ -37,4 +37,5 @@ extern void initNeighbor(Neighbor*, Parameter*); extern void setupNeighbor(); extern void binatoms(Atom*); extern void buildNeighbor(Atom*, Neighbor*); +extern void sortAtom(Atom*); #endif diff --git a/src/main.c b/src/main.c index 6dbb648..5c20a7f 100644 --- a/src/main.c +++ b/src/main.c @@ -46,7 +46,7 @@ extern double computeForce(Parameter*, Atom*, Neighbor*); extern double computeForceTracing(Parameter*, Atom*, Neighbor*, Stats*, int, int); -extern double computeForceEam(Eam* eam, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep); +extern double computeForceEam(Eam* eam, Parameter*, Atom *atom, Neighbor *neighbor, Stats *stats, int first_exec, int timestep); void init(Parameter *param) { @@ -89,7 +89,7 @@ double setup( S = getTimeStamp(); initAtom(atom); initNeighbor(neighbor, param); - initPbc(); + initPbc(atom); initStats(stats); setupNeighbor(); createAtom(atom, param); @@ -115,7 +115,7 @@ double reneighbour( updateAtomsPbc(atom, param); setupPbc(atom, param); updatePbc(atom, param); - /* sortAtom(); */ + //sortAtom(atom); buildNeighbor(atom, neighbor); LIKWID_MARKER_STOP("reneighbour"); E = getTimeStamp(); @@ -257,7 +257,7 @@ int main(int argc, char** argv) setup(¶m, &eam, &atom, &neighbor, &stats); computeThermo(0, ¶m, &atom); if(param.force_field == FF_EAM) { - computeForceEam(&eam, &atom, &neighbor, &stats, 1, 0); + 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); @@ -284,7 +284,7 @@ int main(int argc, char** argv) } if(param.force_field == FF_EAM) { - timer[FORCE] += computeForceEam(&eam, &atom, &neighbor, &stats, 0, n + 1); + 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); diff --git a/src/neighbor.c b/src/neighbor.c index 91c7ab8..4fa3416 100644 --- a/src/neighbor.c +++ b/src/neighbor.c @@ -340,3 +340,57 @@ void binatoms(Atom *atom) } } } + +void sortAtom(Atom* atom) { + binatoms(atom); + int Nmax = atom->Nmax; + int* binpos = bincount; + + for(int i=1; ix; double* old_y = atom->y; double* old_z = atom->z; + double* old_vx = atom->vx; double* old_vy = atom->vy; double* old_vz = atom->vz; + + for(int mybin = 0; mybin0?binpos[mybin-1]:0; + int count = binpos[mybin] - start; + for(int k=0; kx); + atom->x = new_x; +#ifndef AOS + free(atom->y); + free(atom->z); + atom->y = new_y; atom->z = new_z; +#endif + free(atom->vx); free(atom->vy); free(atom->vz); + atom->vx = new_vx; atom->vy = new_vy; atom->vz = new_vz; +} diff --git a/src/pbc.c b/src/pbc.c index f2bbb19..3610701 100644 --- a/src/pbc.c +++ b/src/pbc.c @@ -30,16 +30,15 @@ #define DELTA 20000 static int NmaxGhost; -static int *BorderMap; static int *PBCx, *PBCy, *PBCz; -static void growPbc(); +static void growPbc(Atom*); /* exported subroutines */ -void initPbc() +void initPbc(Atom* atom) { NmaxGhost = 0; - BorderMap = NULL; + atom->border_map = NULL; PBCx = NULL; PBCy = NULL; PBCz = NULL; } @@ -47,15 +46,16 @@ void initPbc() /* uses mapping created in setupPbc */ void updatePbc(Atom *atom, Parameter *param) { + int *border_map = atom->border_map; int nlocal = atom->Nlocal; MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; MD_FLOAT zprd = param->zprd; for(int i = 0; i < atom->Nghost; i++) { - atom_x(nlocal + i) = atom_x(BorderMap[i]) + PBCx[i] * xprd; - atom_y(nlocal + i) = atom_y(BorderMap[i]) + PBCy[i] * yprd; - atom_z(nlocal + i) = atom_z(BorderMap[i]) + PBCz[i] * zprd; + atom_x(nlocal + i) = atom_x(border_map[i]) + PBCx[i] * xprd; + atom_y(nlocal + i) = atom_y(border_map[i]) + PBCy[i] * yprd; + atom_z(nlocal + i) = atom_z(border_map[i]) + PBCz[i] * zprd; } } @@ -95,7 +95,7 @@ void updateAtomsPbc(Atom *atom, Parameter *param) * that are then enforced in updatePbc */ #define ADDGHOST(dx,dy,dz) \ Nghost++; \ - BorderMap[Nghost] = i; \ + border_map[Nghost] = i; \ PBCx[Nghost] = dx; \ PBCy[Nghost] = dy; \ PBCz[Nghost] = dz; \ @@ -103,6 +103,7 @@ void updateAtomsPbc(Atom *atom, Parameter *param) void setupPbc(Atom *atom, Parameter *param) { + int *border_map = atom->border_map; MD_FLOAT xprd = param->xprd; MD_FLOAT yprd = param->yprd; MD_FLOAT zprd = param->zprd; @@ -115,7 +116,8 @@ void setupPbc(Atom *atom, Parameter *param) growAtom(atom); } if (Nghost + 7 >= NmaxGhost) { - growPbc(); + growPbc(atom); + border_map = atom->border_map; } MD_FLOAT x = atom_x(i); @@ -158,12 +160,12 @@ void setupPbc(Atom *atom, Parameter *param) } /* internal subroutines */ -void growPbc() +void growPbc(Atom* atom) { int nold = NmaxGhost; NmaxGhost += DELTA; - BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + atom->border_map = (int*) reallocate(atom->border_map, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));