Apply first changes to DEM kernel

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

lammps/force_dem.c

@@ -37,9 +37,7 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
     int* neighs;
     MD_FLOAT k_s = param->k_s;
     MD_FLOAT k_dn = param->k_dn;
-#ifndef EXPLICIT_TYPES
     MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
-#endif
 
     for(int i = 0; i < Nlocal; i++) {
         atom_fx(i) = 0.0;
@@ -50,7 +48,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
     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];
@@ -62,10 +59,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
         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];
@@ -77,32 +70,11 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
             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 r = sqrt(rsq);
-                // penetration depth
                 MD_FLOAT p = irad + jrad - r;
-                if(p >= 0) {
-                    // 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
+                if(p > 0) {
-                    //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 delvx = atom_vx(i) - atom_vx(j);
                     MD_FLOAT delvy = atom_vy(i) - atom_vy(j);
                     MD_FLOAT delvz = atom_vz(i) - atom_vz(j);
@@ -119,9 +91,28 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
                     MD_FLOAT nvz = delvz / vr;
 
                     // forces
-                    fix += k_s * p * nx - k_dn * nvx;
+                    atom_fx(i) += k_s * p * nx - k_dn * nvx;
-                    fiy += k_s * p * ny - k_dn * nvy;
+                    atom_fy(i) += k_s * p * ny - k_dn * nvy;
-                    fiz += k_s * p * nz - k_dn * nvz;
+                    atom_fz(i) += k_s * p * nz - k_dn * nvz;
+                    atom_fx(j) += -k_s * p * nx - k_dn * nvx;
+                    atom_fy(j) += -k_s * p * ny - k_dn * nvy;
+                    atom_fz(j) += -k_s * p * nz - k_dn * nvz;
+
+                    // 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));
 
                     // tangential force
                     //fix += MIN(kdt * vtsq, kf * fnx) * tx;
@@ -140,10 +131,6 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
             }
         }
 
-        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);
     }
@@ -152,84 +139,3 @@ double computeForceDemFullNeigh(Parameter *param, Atom *atom, Neighbor *neighbor
     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;
-}

lammps/main.c

@@ -50,7 +50,6 @@ extern double computeForceLJFullNeigh_simd(Parameter*, Atom*, Neighbor*, Stats*)
 extern double computeForceLJHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*);
 extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*);
 extern double computeForceDemFullNeigh(Parameter*, Atom*, Neighbor*, Stats*);
-extern double computeForceDemHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*);
 
 #ifdef USE_SIMD_KERNEL
 #   define KERNEL_NAME              "SIMD"
@@ -135,7 +134,8 @@ double computeForce(Eam *eam, Parameter *param, Atom *atom, Neighbor *neighbor,
         return computeForceEam(eam, param, atom, neighbor, stats);
     } else if(param->force_field == FF_DEM) {
         if(param->half_neigh) {
-            return computeForceDemHalfNeigh(param, atom, neighbor, stats);
+            fprintf(stderr, "Error: DEM cannot use half neighbor-lists!\n");
+            return 0.0;
         } else {
             return computeForceDemFullNeigh(param, atom, neighbor, stats);
         }

lammps/parameter.c

@@ -130,6 +130,9 @@ void readParameter(Parameter *param, const char *filename) {
         }
     }
 
+    // Update dtforce
+    param->dtforce = 0.5 * param->dt;
+
     // Update sigma6 parameter
     MD_FLOAT s2 = param->sigma * param->sigma;
     param->sigma6 = s2 * s2 * s2;