moebiusband/MD-Bench
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add EAM without explicit types and update fp for PBC atoms

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2021-11-03 00:57:24 +01:00
parent 0f1e824507
commit ec556eb117
8 changed files with 120 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
src/atom.c
View File

@ -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; i<mbins; i++) { */
/* binpos[i] += binpos[i-1]; */
/* } */
/* double* new_x = (double*) malloc(Nmax * sizeof(double)); */
/* double* new_y = (double*) malloc(Nmax * sizeof(double)); */
/* double* new_z = (double*) malloc(Nmax * sizeof(double)); */
/* double* new_vx = (double*) malloc(Nmax * sizeof(double)); */
/* double* new_vy = (double*) malloc(Nmax * sizeof(double)); */
/* double* new_vz = (double*) malloc(Nmax * sizeof(double)); */
/* double* old_x = x; double* old_y = y; double* old_z = z; */
/* double* old_vx = vx; double* old_vy = vy; double* old_vz = vz; */
/* for(int mybin = 0; mybin<mbins; mybin++) { */
/* int start = mybin>0?binpos[mybin-1]:0; */
/* int count = binpos[mybin] - start; */
/* for(int k=0; k<count; k++) { */
/* int new_i = start + k; */
/* int old_i = bins[mybin * atoms_per_bin + k]; */
/* new_x[new_i] = old_x[old_i]; */
/* new_y[new_i] = old_y[old_i]; */
/* new_z[new_i] = old_z[old_i]; */
/* new_vx[new_i] = old_vx[old_i]; */
/* new_vy[new_i] = old_vy[old_i]; */
/* new_vz[new_i] = old_vz[old_i]; */
/* } */
/* } */
/* free(x); free(y); free(z); */
/* free(vx); free(vy); free(vz); */
/* x = new_x; y = new_y; z = new_z; */
/* vx = new_vx; vy = new_vy; vz = new_vz; */
/* } */

50
src/force_eam.c
View File

@ -32,7 +32,7 @@
#include <eam.h>
#include <util.h>
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;

1
src/includes/atom.h
View File

@ -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;

1
src/includes/neighbor.h
View File

@ -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

10
src/main.c
View File

@ -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(&param, &eam, &atom, &neighbor, &stats);
computeThermo(0, &param, &atom);
if(param.force_field == FF_EAM) {
computeForceEam(&eam, &atom, &neighbor, &stats, 1, 0);
computeForceEam(&eam, &param, &atom, &neighbor, &stats, 1, 0);
} else {
#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS)
computeForceTracing(&param, &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, &param, &atom, &neighbor, &stats, 0, n + 1);
} else {
#if defined(MEM_TRACER) || defined(INDEX_TRACER) || defined(COMPUTE_STATS)
timer[FORCE] += computeForceTracing(&param, &atom, &neighbor, &stats, 0, n + 1);

54
src/neighbor.c
View File

@ -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; i<mbins; i++) {
binpos[i] += binpos[i-1];
}
#ifdef AOS
double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
#else
double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_y = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_z = (double*) malloc(Nmax * sizeof(MD_FLOAT));
#endif
double* new_vx = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_vy = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* new_vz = (double*) malloc(Nmax * sizeof(MD_FLOAT));
double* old_x = atom->x; 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; mybin<mbins; mybin++) {
int start = mybin>0?binpos[mybin-1]:0;
int count = binpos[mybin] - start;
for(int k=0; k<count; k++) {
int new_i = start + k;
int old_i = bins[mybin * atoms_per_bin + k];
#ifdef AOS
new_x[new_i * 3 + 0] = old_x[old_i * 3 + 0];
new_x[new_i * 3 + 1] = old_x[old_i * 3 + 1];
new_x[new_i * 3 + 2] = old_x[old_i * 3 + 2];
#else
new_x[new_i] = old_x[old_i];
new_y[new_i] = old_y[old_i];
new_z[new_i] = old_z[old_i];
#endif
new_vx[new_i] = old_vx[old_i];
new_vy[new_i] = old_vy[old_i];
new_vz[new_i] = old_vz[old_i];
}
}
free(atom->x);
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;
}

24
src/pbc.c
View File

@ -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));