diff --git a/lammps/atom.c b/lammps/atom.c
index d923b0b..46e4110 100644
--- a/lammps/atom.c
+++ b/lammps/atom.c
@@ -55,7 +55,7 @@ void initAtom(Atom *atom) {
atom->sigma6 = NULL;
atom->cutforcesq = NULL;
atom->cutneighsq = NULL;
- atom->radius = 1.0;
+ atom->radius = NULL;
atom->av = NULL;
atom->r = NULL;
}
@@ -144,9 +144,9 @@ void createAtom(Atom *atom, Parameter *param) {
atom_x(atom->Nlocal) = xtmp;
atom_y(atom->Nlocal) = ytmp;
atom_z(atom->Nlocal) = ztmp;
- atom->vx[atom->Nlocal] = vxtmp;
- atom->vy[atom->Nlocal] = vytmp;
- atom->vz[atom->Nlocal] = vztmp;
+ atom_vx(atom->Nlocal) = vxtmp;
+ atom_vy(atom->Nlocal) = vytmp;
+ atom_vz(atom->Nlocal) = vztmp;
atom->type[atom->Nlocal] = rand() % atom->ntypes;
atom->Nlocal++;
}
@@ -220,9 +220,9 @@ int readAtom_pdb(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
- atom->vx[atom_id] = 0.0;
- atom->vy[atom_id] = 0.0;
- atom->vz[atom_id] = 0.0;
+ atom_vx(atom_id) = 0.0;
+ atom_vy(atom_id) = 0.0;
+ atom_vz(atom_id) = 0.0;
occupancy = atof(strtok(NULL, " "));
charge = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
@@ -299,9 +299,9 @@ int readAtom_gro(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
- atom->vx[atom_id] = atof(strtok(NULL, " "));
- atom->vy[atom_id] = atof(strtok(NULL, " "));
- atom->vz[atom_id] = atof(strtok(NULL, " "));
+ atom_vx(atom_id) = atof(strtok(NULL, " "));
+ atom_vy(atom_id) = atof(strtok(NULL, " "));
+ atom_vz(atom_id) = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
atom->Natoms++;
atom->Nlocal++;
@@ -396,9 +396,9 @@ int readAtom_dmp(Atom* atom, Parameter* param) {
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
- atom->vx[atom_id] = atof(strtok(NULL, " "));
- atom->vy[atom_id] = atof(strtok(NULL, " "));
- atom->vz[atom_id] = atof(strtok(NULL, " "));
+ atom_vx(atom_id) = atof(strtok(NULL, " "));
+ atom_vy(atom_id) = atof(strtok(NULL, " "));
+ atom_vz(atom_id) = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id], atom->ntypes);
read_atoms++;
}
@@ -441,20 +441,22 @@ void growAtom(Atom *atom) {
#ifdef AOS
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
+ atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
#else
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->y = (MD_FLOAT*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
atom->z = (MD_FLOAT*) reallocate(atom->z, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
+ atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
+ 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->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));
#endif
- atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
- 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->radius = (int *) reallocate(atom->radius, ALIGNMENT, atom->Nmax * sizeof(int), nold * sizeof(int));
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
index 8162481..6ba1c70 100644
--- a/lammps/force_dem.c
+++ b/lammps/force_dem.c
@@ -20,6 +20,7 @@
* with MD-Bench. If not, see .
* =======================================================================================
*/
+#include
#include
#include
@@ -34,6 +35,8 @@
double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
int Nlocal = atom->Nlocal;
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
@@ -81,54 +84,54 @@ double computeForceDEMFullNeigh_plain_c(Parameter *param, Atom *atom, Neighbor *
#endif
if(rsq < cutforcesq) {
- MD_FLOAT sq = sqrt(rsq);
+ MD_FLOAT r = sqrt(rsq);
// penetration depth
- MD_FLOAT p = irad + jrad - sq;
- if(p < 0) { continue; }
+ 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;
- // 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;
- // 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 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);
+ MD_FLOAT vr = sqrt(delvx * delvx + delvy * delvy + delvz * delvz);
- // 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 / r;
+ MD_FLOAT ny = dely / r;
+ MD_FLOAT nz = delz / r;
- // normal distance
- MD_FLOAT nx = delx / sq;
- MD_FLOAT ny = dely / sq;
- MD_FLOAT nz = delz / sq;
+ // normal contact velocity
+ MD_FLOAT nvx = delvx / vr;
+ MD_FLOAT nvy = delvy / vr;
+ MD_FLOAT nvz = delvz / vr;
- // 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;
+ // forces
+ fix += k_s * p * nx - k_dn * nvx;
+ fiy += k_s * p * ny - k_dn * nvy;
+ fiz += k_s * p * nz - k_dn * nvz;
+ // tangential force
+ //fix += MIN(kdt * vtsq, kf * fnx) * tx;
+ //fiy += MIN(kdt * vtsq, kf * fny) * ty;
+ //fiz += MIN(kdt * vtsq, kf * fnz) * tz;
+ // 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 {
diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h
index 50b15a1..4c473d8 100644
--- a/lammps/includes/atom.h
+++ b/lammps/includes/atom.h
@@ -38,7 +38,7 @@ typedef struct {
MD_FLOAT *cutforcesq;
MD_FLOAT *cutneighsq;
// DEM
- MD_FLOAT radius;
+ MD_FLOAT *radius;
MD_FLOAT *av;
MD_FLOAT *r;
} Atom;
@@ -56,6 +56,9 @@ extern void growAtom(Atom*);
#define atom_x(i) atom->x[(i) * 3 + 0]
#define atom_y(i) atom->x[(i) * 3 + 1]
#define atom_z(i) atom->x[(i) * 3 + 2]
+#define atom_vx(i) atom->vx[(i) * 3 + 0]
+#define atom_vy(i) atom->vx[(i) * 3 + 1]
+#define atom_vz(i) atom->vx[(i) * 3 + 2]
#define atom_fx(i) atom->fx[(i) * 3 + 0]
#define atom_fy(i) atom->fx[(i) * 3 + 1]
#define atom_fz(i) atom->fx[(i) * 3 + 2]
@@ -64,6 +67,9 @@ extern void growAtom(Atom*);
#define atom_x(i) atom->x[i]
#define atom_y(i) atom->y[i]
#define atom_z(i) atom->z[i]
+#define atom_vx(i) atom->vx[i]
+#define atom_vy(i) atom->vy[i]
+#define atom_vz(i) atom->vz[i]
#define atom_fx(i) atom->fx[i]
#define atom_fy(i) atom->fy[i]
#define atom_fz(i) atom->fz[i]
diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h
index 5c8eba7..497f4ce 100644
--- a/lammps/includes/parameter.h
+++ b/lammps/includes/parameter.h
@@ -58,6 +58,10 @@ typedef struct {
MD_FLOAT xprd, yprd, zprd;
double proc_freq;
char* eam_file;
+ // DEM
+ MD_FLOAT k_s;
+ MD_FLOAT k_dn;
+ MD_FLOAT gx, gy, gz;
} Parameter;
void initParameter(Parameter*);
diff --git a/lammps/main-stub.c b/lammps/main-stub.c
index 2e8f1ff..eb6aeee 100644
--- a/lammps/main-stub.c
+++ b/lammps/main-stub.c
@@ -227,9 +227,9 @@ int main(int argc, const char *argv[]) {
atom_x(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
atom_y(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
atom_z(atom->Nlocal) = (MD_FLOAT)(i) * 0.00001;
- atom->vx[atom->Nlocal] = 0.0;
- atom->vy[atom->Nlocal] = 0.0;
- atom->vz[atom->Nlocal] = 0.0;
+ atom_vx(atom->Nlocal) = 0.0;
+ atom_vy(atom->Nlocal) = 0.0;
+ atom_vz(atom->Nlocal) = 0.0;
atom->Nlocal++;
}
diff --git a/lammps/main.c b/lammps/main.c
index f770fa1..3551ce1 100644
--- a/lammps/main.c
+++ b/lammps/main.c
@@ -102,26 +102,21 @@ double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) {
}
void initialIntegrate(Parameter *param, Atom *atom) {
- MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
-
for(int i = 0; i < atom->Nlocal; i++) {
- vx[i] += param->dtforce * atom_fx(i);
- vy[i] += param->dtforce * atom_fy(i);
- vz[i] += param->dtforce * atom_fz(i);
- atom_x(i) = atom_x(i) + param->dt * vx[i];
- atom_y(i) = atom_y(i) + param->dt * vy[i];
- atom_z(i) = atom_z(i) + param->dt * vz[i];
+ atom_vx(i) += param->dtforce * atom_fx(i);
+ atom_vy(i) += param->dtforce * atom_fy(i);
+ atom_vz(i) += param->dtforce * atom_fz(i);
+ atom_x(i) = atom_x(i) + param->dt * atom_vx(i);
+ atom_y(i) = atom_y(i) + param->dt * atom_vy(i);
+ atom_z(i) = atom_z(i) + param->dt * atom_vz(i);
}
}
void finalIntegrate(Parameter *param, Atom *atom) {
- MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
- MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
-
for(int i = 0; i < atom->Nlocal; i++) {
- vx[i] += param->dtforce * atom_fx(i);
- vy[i] += param->dtforce * atom_fy(i);
- vz[i] += param->dtforce * atom_fz(i);
+ atom_vx(i) += param->dtforce * atom_fx(i);
+ atom_vy(i) += param->dtforce * atom_fy(i);
+ atom_vz(i) += param->dtforce * atom_fz(i);
}
}
diff --git a/lammps/neighbor.c b/lammps/neighbor.c
index 8900a04..826a00c 100644
--- a/lammps/neighbor.c
+++ b/lammps/neighbor.c
@@ -349,14 +349,15 @@ void sortAtom(Atom* atom) {
#ifdef AOS
MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
+ MD_FLOAT* new_vx = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
#else
MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_y = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_z = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
-#endif
MD_FLOAT* new_vx = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_vy = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
MD_FLOAT* new_vz = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+#endif
MD_FLOAT* old_x = atom->x; MD_FLOAT* old_y = atom->y; MD_FLOAT* old_z = atom->z;
MD_FLOAT* old_vx = atom->vx; MD_FLOAT* old_vy = atom->vy; MD_FLOAT* old_vz = atom->vz;
@@ -370,24 +371,32 @@ void sortAtom(Atom* atom) {
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];
+ new_vx[new_i * 3 + 0] = old_vx[old_i * 3 + 0];
+ new_vx[new_i * 3 + 1] = old_vx[old_i * 3 + 1];
+ new_vx[new_i * 3 + 2] = old_vx[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];
+#endif
}
}
free(atom->x);
+ free(atom->vx);
atom->x = new_x;
+ atom->vx = new_vx;
#ifndef AOS
free(atom->y);
free(atom->z);
- atom->y = new_y; atom->z = new_z;
+ free(atom->vy);
+ free(atom->vz);
+ atom->y = new_y;
+ atom->z = new_z;
+ atom->vy = new_vy;
+ atom->vz = new_vz;
#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/lammps/parameter.c b/lammps/parameter.c
index 90fc3cf..e16f538 100644
--- a/lammps/parameter.c
+++ b/lammps/parameter.c
@@ -54,6 +54,12 @@ void initParameter(Parameter *param) {
param->v_out_every = 5;
param->half_neigh = 0;
param->proc_freq = 2.4;
+ // DEM
+ param->k_s = 1.0;
+ param->k_dn = 1.0;
+ param->gx = 0.0;
+ param->gy = 0.0;
+ param->gz = 0.0;
}
void readParameter(Parameter *param, const char *filename) {
diff --git a/lammps/thermo.c b/lammps/thermo.c
index 73e94ae..0fbbf41 100644
--- a/lammps/thermo.c
+++ b/lammps/thermo.c
@@ -70,12 +70,8 @@ void setupThermo(Parameter *param, int natoms)
void computeThermo(int iflag, Parameter *param, Atom *atom)
{
MD_FLOAT t = 0.0, p;
- MD_FLOAT* vx = atom->vx;
- MD_FLOAT* vy = atom->vy;
- MD_FLOAT* vz = atom->vz;
-
for(int i = 0; i < atom->Nlocal; i++) {
- t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
+ t += (atom_vx(i) * atom_vx(i) + atom_vy(i) * atom_vy(i) + atom_vz(i) * atom_vz(i)) * param->mass;
}
t = t * t_scale;
@@ -100,12 +96,11 @@ void adjustThermo(Parameter *param, Atom *atom)
{
/* zero center-of-mass motion */
MD_FLOAT vxtot = 0.0; MD_FLOAT vytot = 0.0; MD_FLOAT vztot = 0.0;
- MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
for(int i = 0; i < atom->Nlocal; i++) {
- vxtot += vx[i];
- vytot += vy[i];
- vztot += vz[i];
+ vxtot += atom_vx(i);
+ vytot += atom_vy(i);
+ vztot += atom_vz(i);
}
vxtot = vxtot / atom->Natoms;
@@ -113,24 +108,24 @@ void adjustThermo(Parameter *param, Atom *atom)
vztot = vztot / atom->Natoms;
for(int i = 0; i < atom->Nlocal; i++) {
- vx[i] -= vxtot;
- vy[i] -= vytot;
- vz[i] -= vztot;
+ atom_vx(i) -= vxtot;
+ atom_vy(i) -= vytot;
+ atom_vz(i) -= vztot;
}
t_act = 0;
MD_FLOAT t = 0.0;
for(int i = 0; i < atom->Nlocal; i++) {
- t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
+ t += (atom_vx(i) * atom_vx(i) + atom_vy(i) * atom_vy(i) + atom_vz(i) * atom_vz(i)) * param->mass;
}
t *= t_scale;
MD_FLOAT factor = sqrt(param->temp / t);
for(int i = 0; i < atom->Nlocal; i++) {
- vx[i] *= factor;
- vy[i] *= factor;
- vz[i] *= factor;
+ atom_vx(i) *= factor;
+ atom_vy(i) *= factor;
+ atom_vz(i) *= factor;
}
}