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

Add version with explicit types for atoms

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2021-05-19 23:51:02 +02:00
parent f7f7ae2002
commit 4496e91125
7 changed files with 91 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
src/atom.c
View File

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

27
src/force.c
View File

@ -27,14 +27,16 @@
#include <parameter.h>
#include <atom.h>
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;
}

8
src/includes/atom.h
View File

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

1
src/includes/parameter.h
View File

@ -40,6 +40,7 @@ typedef struct {
MD_FLOAT temp;
MD_FLOAT rho;
MD_FLOAT mass;
int ntypes;
int ntimes;
int nstat;
int every;

23
src/main-stub.c
View File

@ -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(&param, atom, &neighbor);
computeForce(&param, atom, &neighbor, 1);
double S, E;
S = getTimeStamp();
LIKWID_MARKER_START("force");
for(int i = 0; i < param.ntimes; i++) {
computeForce(&param, atom, &neighbor);
computeForce(&param, atom, &neighbor, INTERNAL_LOOP_NTIMES);
}
LIKWID_MARKER_STOP("force");
E = getTimeStamp();

7
src/main.c
View File

@ -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(&param, &atom, &neighbor);
computeThermo(0, &param, &atom);
computeForce(&param, &atom, &neighbor);
computeForce(&param, &atom, &neighbor, 1);
timer[FORCE] = 0.0;
timer[NEIGH] = 0.0;
@ -221,7 +222,7 @@ int main (int argc, char** argv)
timer[NEIGH] += reneighbour(&param, &atom, &neighbor);
}
timer[FORCE] += computeForce(&param, &atom, &neighbor);
timer[FORCE] += computeForce(&param, &atom, &neighbor, 1);
finalIntegrate(&param, &atom);
if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) {

12
src/neighbor.c
View File

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