diff --git a/lammps/atom.c b/lammps/atom.c
index d6fac51..d923b0b 100644
--- a/lammps/atom.c
+++ b/lammps/atom.c
@@ -41,8 +41,7 @@
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#endif
-void initAtom(Atom *atom)
-{
+void initAtom(Atom *atom) {
atom->x = NULL; atom->y = NULL; atom->z = NULL;
atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
@@ -56,10 +55,12 @@ void initAtom(Atom *atom)
atom->sigma6 = NULL;
atom->cutforcesq = NULL;
atom->cutneighsq = NULL;
+ atom->radius = 1.0;
+ atom->av = NULL;
+ atom->r = NULL;
}
-void createAtom(Atom *atom, Parameter *param)
-{
+void createAtom(Atom *atom, Parameter *param) {
MD_FLOAT xlo = 0.0; MD_FLOAT xhi = param->xprd;
MD_FLOAT ylo = 0.0; MD_FLOAT yhi = param->yprd;
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd;
@@ -453,4 +454,8 @@ void growAtom(Atom *atom) {
atom->vy = (MD_FLOAT*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->vz = (MD_FLOAT*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->type = (int *) reallocate(atom->type, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
+ #ifdef DEM
+ atom->av = (MD_FLOAT*) reallocate(atom->av, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
+ atom->r = (MD_FLOAT*) reallocate(atom->r, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 4, nold * sizeof(MD_FLOAT) * 4);
+ #endif
}
diff --git a/lammps/force_dem.c b/lammps/force_dem.c
new file mode 100644
index 0000000..8162481
--- /dev/null
+++ b/lammps/force_dem.c
@@ -0,0 +1,232 @@
+/*
+ * =======================================================================================
+ *
+ * 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
+#include
+#include
+
+// TODO: Joint common files for gromacs and lammps variants
+#include "../gromacs/includes/simd.h"
+
+double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
+ int Nlocal = atom->Nlocal;
+ int* neighs;
+#ifndef EXPLICIT_TYPES
+ MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
+#endif
+
+ for(int i = 0; i < Nlocal; i++) {
+ atom_fx(i) = 0.0;
+ atom_fy(i) = 0.0;
+ atom_fz(i) = 0.0;
+ }
+
+ double S = getTimeStamp();
+ LIKWID_MARKER_START("force");
+
+#pragma omp parallel for
+ for(int i = 0; i < Nlocal; i++) {
+ neighs = &neighbor->neighbors[i * neighbor->maxneighs];
+ int numneighs = neighbor->numneigh[i];
+ MD_FLOAT irad = atom->radius[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;
+
+#ifdef EXPLICIT_TYPES
+ const int type_i = atom->type[i];
+#endif
+
+ for(int k = 0; k < numneighs; k++) {
+ int j = neighs[k];
+ MD_FLOAT jrad = atom->radius[j];
+ MD_FLOAT xj = atom_x(j);
+ MD_FLOAT yj = atom_y(j);
+ MD_FLOAT zj = atom_z(j);
+ MD_FLOAT delx = xtmp - xj;
+ MD_FLOAT dely = ytmp - yj;
+ MD_FLOAT delz = ztmp - zj;
+ 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];
+#endif
+
+ if(rsq < cutforcesq) {
+ MD_FLOAT sq = sqrt(rsq);
+ // penetration depth
+ MD_FLOAT p = irad + jrad - sq;
+ if(p < 0) { continue; }
+
+ // contact position
+ MD_FLOAT cterm = jrad / (irad + jrad);
+ MD_FLOAT cx = xj + cterm * delx;
+ MD_FLOAT cy = yj + cterm * dely;
+ MD_FLOAT cz = zj + cterm * delz;
+
+ // delta contact and particle position
+ MD_FLOAT delcx = cx - xtmp;
+ MD_FLOAT delcy = cy - ytmp;
+ MD_FLOAT delcz = cz - ztmp;
+
+ // contact velocity
+ MD_FLOAT cvx = (atom_vx(i) + atom_avx(i) * delcx) - (atom_vx(j) + atom_avx(j) * (cx - xj));
+ MD_FLOAT cvy = (atom_vy(i) + atom_avy(i) * delcy) - (atom_vy(j) + atom_avy(j) * (cy - yj));
+ MD_FLOAT cvz = (atom_vz(i) + atom_avz(i) * delcz) - (atom_vz(j) + atom_avz(j) * (cz - zj));
+ MD_FLOAT vsq = sqrt(cvx * cvx + cvy * cvy + cvz * cvz);
+
+ // normal distance
+ MD_FLOAT nx = delx / sq;
+ MD_FLOAT ny = dely / sq;
+ MD_FLOAT nz = delz / sq;
+
+ // normal contact velocity
+ MD_FLOAT vnx = cvx / vsq;
+ MD_FLOAT vny = cvy / vsq;
+ MD_FLOAT vnz = cvz / vsq;
+
+ // forces
+ MD_FLOAT fnx = ks * p * nx - kdn * vn;
+ MD_FLOAT fny = ks * p * ny - kdn * vn;
+ MD_FLOAT fnz = ks * p * nz - kdn * vn;
+ MD_FLOAT ftx = MIN(kdt * vtsq, kf * fnx) * tx;
+ MD_FLOAT fty = MIN(kdt * vtsq, kf * fny) * ty;
+ MD_FLOAT ftz = MIN(kdt * vtsq, kf * fnz) * tz;
+ fix += fnx + ftx;
+ fiy += fny + fty;
+ fiz += fnz + ftz;
+
+ // torque
+ //MD_FLOAT taux = delcx * ftx;
+ //MD_FLOAT tauy = delcy * fty;
+ //MD_FLOAT tauz = delcz * ftz;
+
+#ifdef USE_REFERENCE_VERSION
+ addStat(stats->atoms_within_cutoff, 1);
+ } else {
+ addStat(stats->atoms_outside_cutoff, 1);
+#endif
+ }
+ }
+
+ atom_fx(i) += fix;
+ atom_fy(i) += fiy;
+ atom_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;
+}
+
+double computeForceDEMHalfNeigh(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
+ int Nlocal = atom->Nlocal;
+ int* neighs;
+#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++) {
+ atom_fx(i) = 0.0;
+ atom_fy(i) = 0.0;
+ atom_fz(i) = 0.0;
+ }
+
+ double S = getTimeStamp();
+ LIKWID_MARKER_START("forceLJ-halfneigh");
+
+ 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;
+
+#ifdef EXPLICIT_TYPES
+ const int type_i = atom->type[i];
+#endif
+
+ // Pragma required to vectorize the inner loop
+#ifdef ENABLE_OMP_SIMD
+ #pragma omp simd reduction(+: fix,fiy,fiz)
+#endif
+ 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;
+ 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;
+
+ // We do not need to update forces for ghost atoms
+ if(j < Nlocal) {
+ atom_fx(j) -= delx * force;
+ atom_fy(j) -= dely * force;
+ atom_fz(j) -= delz * force;
+ }
+ }
+ }
+
+ atom_fx(i) += fix;
+ atom_fy(i) += fiy;
+ atom_fz(i) += fiz;
+
+ addStat(stats->total_force_neighs, numneighs);
+ addStat(stats->total_force_iters, (numneighs + VECTOR_WIDTH - 1) / VECTOR_WIDTH);
+ }
+
+ LIKWID_MARKER_STOP("forceLJ-halfneigh");
+ double E = getTimeStamp();
+ return E-S;
+}
diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h
index cf161c9..50b15a1 100644
--- a/lammps/includes/atom.h
+++ b/lammps/includes/atom.h
@@ -37,6 +37,10 @@ typedef struct {
MD_FLOAT *sigma6;
MD_FLOAT *cutforcesq;
MD_FLOAT *cutneighsq;
+ // DEM
+ MD_FLOAT radius;
+ MD_FLOAT *av;
+ MD_FLOAT *r;
} Atom;
extern void initAtom(Atom*);