diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 685d222539..2aadb57e10 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -274,8 +274,8 @@ void FixPour::init() // check if a shear history fix exists fix_history = NULL; - if (force->pair_match("gran/hooke/history") || - force->pair_match("gran/hertz/history")) + if (force->pair_match("gran/hooke/history",1) || + force->pair_match("gran/hertz/history",1)) for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"SHEAR_HISTORY") == 0) fix_history = (FixShearHistory *) modify->fix[i]; diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index b371821ed4..5c8ffae419 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER}; -enum{NO_HISTORY,HISTORY,HERTZIAN}; +enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY}; #define BIG 1.0e20 @@ -190,25 +190,17 @@ int FixWallGran::setmask() void FixWallGran::init() { - // set local values from Pair values - - if (force->pair == NULL) - error->all("Fix wall/gran is incompatible with Pair style"); - double *p_xkk = (double *) force->pair->extract("xkk"); - if (!p_xkk) error->all("Fix wall/gran is incompatible with Pair style"); - xkk = *p_xkk; - - // same initialization as in pair_gran_history::init_style() - - xkkt = xkk * 2.0/7.0; - double gammas = 0.5*gamman; dt = update->dt; // set pairstyle from granular pair style - if (force->pair_match("gran/no_history")) pairstyle = NO_HISTORY; - else if (force->pair_match("gran/history")) pairstyle = HISTORY; - else if (force->pair_match("gran/hertzian")) pairstyle = HERTZIAN; + if (force->pair_match("gran/hooke",1)) + pairstyle = HOOKE; + else if (force->pair_match("gran/hooke/history",1)) + pairstyle = HOOKE_HISTORY; + else if (force->pair_match("gran/hertz/history",1)) + pairstyle = HERTZ_HISTORY; + else error->all("Fix wall/gran is incompatible with Pair style"); } /* ---------------------------------------------------------------------- */ @@ -295,21 +287,21 @@ void FixWallGran::post_force(int vflag) rsq = dx*dx + dy*dy + dz*dz; if (rsq > radius[i]*radius[i]) { - if (pairstyle != NO_HISTORY) { + if (pairstyle != HOOKE) { shear[i][0] = 0.0; shear[i][1] = 0.0; shear[i][2] = 0.0; } } else { - if (pairstyle == NO_HISTORY) - no_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],rmass[i]); - else if (pairstyle == HISTORY) - history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],rmass[i],shear[i]); - else if (pairstyle == HERTZIAN) - hertzian(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],rmass[i],shear[i]); + if (pairstyle == HOOKE) + hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], + radius[i],rmass[i]); + else if (pairstyle == HOOKE_HISTORY) + hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], + radius[i],rmass[i],shear[i]); + else if (pairstyle == HERTZ_HISTORY) + hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], + radius[i],rmass[i],shear[i]); } } } @@ -317,13 +309,13 @@ void FixWallGran::post_force(int vflag) /* ---------------------------------------------------------------------- */ -void FixWallGran::no_history(double rsq, double dx, double dy, double dz, - double *vwall, double *v, - double *f, double *omega, double *torque, - double radius, double mass) +void FixWallGran::hooke(double rsq, double dx, double dy, double dz, + double *vwall, double *v, + double *f, double *omega, double *torque, + double radius, double mass) { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; + double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel; double fn,fs,ft,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3,rinv,rsqinv; r = sqrt(rsq); @@ -357,9 +349,9 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, // normal forces = Hookian contact + normal velocity damping - xmeff = mass; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radius-r)*rinv - damp; + meff = mass; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radius-r)*rinv - damp; // relative velocities @@ -372,7 +364,7 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, // force normalization fn = xmu * fabs(ccel*r); - fs = xmeff*gammas*vrel; + fs = meff*gammat*vrel; if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; @@ -402,13 +394,13 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, /* ---------------------------------------------------------------------- */ -void FixWallGran::history(double rsq, double dx, double dy, double dz, - double *vwall, double *v, - double *f, double *omega, double *torque, - double radius, double mass, double *shear) +void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, + double *vwall, double *v, + double *f, double *omega, double *torque, + double radius, double mass, double *shear) { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; + double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel; double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; double shrmag,rsht,rinv,rsqinv; @@ -443,9 +435,9 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, // normal forces = Hookian contact + normal velocity damping - xmeff = mass; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radius-r)*rinv - damp; + meff = mass; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radius-r)*rinv - damp; // relative velocities @@ -472,9 +464,9 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, // tangential forces = shear + tangential velocity damping - fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1); - fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2); - fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3); + fs1 = - (kt*shear[0] + meff*gammat*vtr1); + fs2 = - (kt*shear[1] + meff*gammat*vtr2); + fs3 = - (kt*shear[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -483,15 +475,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) - - xmeff*gammas*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) - - xmeff*gammas*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) - - xmeff*gammas*vtr3/xkkt; - fs1 = fs1 * fn / fs ; - fs2 = fs2 * fn / fs; - fs3 = fs3 * fn / fs; + shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + meff*gammat*vtr1/kt; + shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + meff*gammat*vtr2/kt; + shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + meff*gammat*vtr3/kt; + fs1 *= fn/fs ; + fs2 *= fn/fs; + fs3 *= fn/fs; } else fs1 = fs2 = fs3 = 0.0; } @@ -515,15 +507,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, /* ---------------------------------------------------------------------- */ -void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, - double *vwall, double *v, - double *f, double *omega, double *torque, - double radius, double mass, double *shear) +void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, + double *vwall, double *v, + double *f, double *omega, double *torque, + double radius, double mass, double *shear) { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; + double wr1,wr2,wr3,meff,damp,ccel,vtr1,vtr2,vtr3,vrel; double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; - double shrmag,rsht,rhertz,rinv,rsqinv; + double shrmag,rsht,polyhertz,rinv,rsqinv; r = sqrt(rsq); rinv = 1.0/r; @@ -556,11 +548,11 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, // normal forces = Hertzian contact + normal velocity damping - xmeff = mass; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radius-r)*rinv - damp; - rhertz = sqrt(radius - r); - ccel = rhertz * ccel; + meff = mass; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radius-r)*rinv - damp; + polyhertz = sqrt((radius-r)*radius); + ccel *= polyhertz; // relative velocities @@ -587,9 +579,9 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, // tangential forces = shear + tangential velocity damping - fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1); - fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2); - fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3); + fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1); + fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2); + fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -598,15 +590,15 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) - - xmeff*gammas*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) - - xmeff*gammas*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) - - xmeff*gammas*vtr3/xkkt; - fs1 = fs1 * fn / fs ; - fs2 = fs2 * fn / fs; - fs3 = fs3 * fn / fs; + shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + meff*gammat*vtr1/kt; + shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + meff*gammat*vtr2/kt; + shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + meff*gammat*vtr3/kt; + fs1 *= fn/fs ; + fs2 *= fn/fs; + fs3 *= fn/fs; } else fs1 = fs2 = fs3 = 0.0; } diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 11108def0c..f2db1263a5 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -40,7 +40,7 @@ class FixWallGran : public Fix { private: int wallstyle,pairstyle,wiggle,wshear,axis; - double xkk,xkkt,gamman,gammas,xmu; + double kn,kt,gamman,gammat,xmu; double lo,hi,cylradius; double dt; double amplitude,period,omega,time_origin,vshear; @@ -48,14 +48,14 @@ class FixWallGran : public Fix { int *touch; double **shear; - void no_history(double, double, double, double, double *, - double *, double *, double *, double *, double, double); - void history(double, double, double, double, double *, - double *, double *, double *, double *, double, double, - double *); - void hertzian(double, double, double, double, double *, - double *, double *, double *, double *, double, double, - double *); + void hooke(double, double, double, double, double *, + double *, double *, double *, double *, double, double); + void hooke_history(double, double, double, double, double *, + double *, double *, double *, double *, double, double, + double *); + void hertz_history(double, double, double, double, double *, + double *, double *, double *, double *, double, double, + double *); }; } diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 9592e8d5b1..1563a09b35 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -16,12 +16,14 @@ ------------------------------------------------------------------------- */ #include "math.h" +#include "stdlib.h" #include "stdio.h" #include "string.h" #include "pair_gran_hertz_history.h" #include "atom.h" #include "force.h" #include "neigh_list.h" +#include "error.h" using namespace LAMMPS_NS; @@ -47,9 +49,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag) double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; - double xmeff,damp,ccel,tor1,tor2,tor3; + double meff,damp,ccel,tor1,tor2,tor3; double fn,fs,fs1,fs2,fs3; - double shrmag,rsht,rhertz; + double shrmag,rsht,polyhertz; int *ilist,*jlist,*numneigh,**firstneigh; int *touch,**firsttouch; double *shear,*allshear,**firstshear; @@ -139,13 +141,13 @@ void PairGranHertzHistory::compute(int eflag, int vflag) // normal force = Hertzian contact + normal velocity damping - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radsum-r)*rinv - damp; - rhertz = sqrt(radsum - r); - ccel *= rhertz; + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; + polyhertz = sqrt((radsum-r)*radi*radj / radsum); + ccel *= polyhertz; // relative velocities @@ -175,9 +177,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag) // tangential forces = shear + tangential velocity damping - fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1); - fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2); - fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3); + fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1); + fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2); + fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -186,20 +188,16 @@ void PairGranHertzHistory::compute(int eflag, int vflag) if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) - - xmeff*gammas*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) - - xmeff*gammas*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) - - xmeff*gammas*vtr3/xkkt; + shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + meff*gammat*vtr1/kt; + shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + meff*gammat*vtr2/kt; + shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + meff*gammat*vtr3/kt; fs1 *= fn/fs; fs2 *= fn/fs; fs3 *= fn/fs; - } else { - fs1 = 0.0; - fs2 = 0.0; - fs3 = 0.0; - } + } else fs1 = fs2 = fs3 = 0.0; } // forces & torques @@ -233,3 +231,33 @@ void PairGranHertzHistory::compute(int eflag, int vflag) } } } + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairGranHertzHistory::settings(int narg, char **arg) +{ + if (narg != 6) error->all("Illegal pair_style command"); + + kn = atof(arg[0]); + if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; + else kt = atof(arg[1]); + + gamman = atof(arg[2]); + if (strcmp(arg[3],"NULL") == 0) gammat = 0.5 * gamman; + else gammat = atof(arg[3]); + + xmu = atof(arg[4]); + dampflag = atoi(arg[5]); + if (dampflag == 0) gammat = 0.0; + + if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || + xmu < 0.0 || xmu > 1.0 || dampflag < 0 || dampflag > 1) + error->all("Illegal pair_style command"); + + // convert Kn and Kt from pressure units to force/distance^2 + + kn /= force->nktv2p; + kt /= force->nktv2p; +} diff --git a/src/GRANULAR/pair_gran_hertz_history.h b/src/GRANULAR/pair_gran_hertz_history.h index 08f388f5fd..a5799a82f5 100644 --- a/src/GRANULAR/pair_gran_hertz_history.h +++ b/src/GRANULAR/pair_gran_hertz_history.h @@ -22,6 +22,7 @@ class PairGranHertzHistory : public PairGranHookeHistory { public: PairGranHertzHistory(class LAMMPS *); void compute(int, int); + void settings(int, char **); }; } diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 4ff0c1d8be..3309583841 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -45,7 +45,7 @@ void PairGranHooke::compute(int eflag, int vflag) double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; - double xmeff,damp,ccel,tor1,tor2,tor3; + double meff,damp,ccel,tor1,tor2,tor3; double fn,fs,ft,fs1,fs2,fs3; int *ilist,*jlist,*numneigh,**firstneigh; @@ -121,11 +121,11 @@ void PairGranHooke::compute(int eflag, int vflag) // normal forces = Hookian contact + normal velocity damping - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radsum-r)*rinv - damp; + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; // relative velocities @@ -138,7 +138,7 @@ void PairGranHooke::compute(int eflag, int vflag) // force normalization fn = xmu * fabs(ccel*r); - fs = xmeff*gammas*vrel; + fs = meff*gammat*vrel; if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 0a38177874..69691e1d92 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -73,7 +73,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; - double xmeff,damp,ccel,tor1,tor2,tor3; + double meff,damp,ccel,tor1,tor2,tor3; double fn,fs,fs1,fs2,fs3; double shrmag,rsht; int *ilist,*jlist,*numneigh,**firstneigh; @@ -165,11 +165,11 @@ void PairGranHookeHistory::compute(int eflag, int vflag) // normal forces = Hookian contact + normal velocity damping - xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); - if (mask[i] & freeze_group_bit) xmeff = rmass[j]; - if (mask[j] & freeze_group_bit) xmeff = rmass[i]; - damp = xmeff*gamman*vnnr*rsqinv; - ccel = xkk*(radsum-r)*rinv - damp; + meff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]); + if (mask[i] & freeze_group_bit) meff = rmass[j]; + if (mask[j] & freeze_group_bit) meff = rmass[i]; + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; // relative velocities @@ -199,9 +199,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag) // tangential forces = shear + tangential velocity damping - fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1); - fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2); - fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3); + fs1 = - (kt*shear[0] + meff*gammat*vtr1); + fs2 = - (kt*shear[1] + meff*gammat*vtr2); + fs3 = - (kt*shear[2] + meff*gammat*vtr3); // rescale frictional displacements and forces if needed @@ -210,20 +210,16 @@ void PairGranHookeHistory::compute(int eflag, int vflag) if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) - - xmeff*gammas*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) - - xmeff*gammas*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) - - xmeff*gammas*vtr3/xkkt; + shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - + meff*gammat*vtr1/kt; + shear[1] = (fn/fs) * (shear[1] + meff*gammat*vtr2/kt) - + meff*gammat*vtr2/kt; + shear[2] = (fn/fs) * (shear[2] + meff*gammat*vtr3/kt) - + meff*gammat*vtr3/kt; fs1 *= fn/fs; fs2 *= fn/fs; fs3 *= fn/fs; - } else { - fs1 = 0.0; - fs2 = 0.0; - fs3 = 0.0; - } + } else fs1 = fs2 = fs3 = 0.0; } // forces & torques @@ -281,12 +277,23 @@ void PairGranHookeHistory::allocate() void PairGranHookeHistory::settings(int narg, char **arg) { - if (narg != 4) error->all("Illegal pair_style command"); + if (narg != 6) error->all("Illegal pair_style command"); - xkk = atof(arg[0]); - gamman = atof(arg[1]); - xmu = atof(arg[2]); - dampflag = atoi(arg[3]); + kn = atof(arg[0]); + if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0; + else kt = atof(arg[1]); + + gamman = atof(arg[2]); + if (strcmp(arg[3],"NULL") == 0) gammat = 0.5 * gamman; + else gammat = atof(arg[3]); + + xmu = atof(arg[4]); + dampflag = atoi(arg[5]); + if (dampflag == 0) gammat = 0.0; + + if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || + xmu < 0.0 || xmu > 1.0 || dampflag < 0 || dampflag > 1) + error->all("Illegal pair_style command"); } /* ---------------------------------------------------------------------- @@ -339,9 +346,6 @@ void PairGranHookeHistory::init_style() neighbor->requests[irequest]->dnum = 3; } - xkkt = xkk * 2.0/7.0; - double gammas = 0.5*gamman; - if (dampflag == 0) gammas = 0.0; dt = update->dt; // if shear history is stored: @@ -466,8 +470,10 @@ void PairGranHookeHistory::read_restart(FILE *fp) void PairGranHookeHistory::write_restart_settings(FILE *fp) { - fwrite(&xkk,sizeof(double),1,fp); + fwrite(&kn,sizeof(double),1,fp); + fwrite(&kt,sizeof(double),1,fp); fwrite(&gamman,sizeof(double),1,fp); + fwrite(&gammat,sizeof(double),1,fp); fwrite(&xmu,sizeof(double),1,fp); fwrite(&dampflag,sizeof(int),1,fp); } @@ -479,30 +485,23 @@ void PairGranHookeHistory::write_restart_settings(FILE *fp) void PairGranHookeHistory::read_restart_settings(FILE *fp) { if (comm->me == 0) { - fread(&xkk,sizeof(double),1,fp); + fread(&kn,sizeof(double),1,fp); + fread(&kt,sizeof(double),1,fp); fread(&gamman,sizeof(double),1,fp); + fread(&gammat,sizeof(double),1,fp); fread(&xmu,sizeof(double),1,fp); fread(&dampflag,sizeof(int),1,fp); } - MPI_Bcast(&xkk,1,MPI_DOUBLE,0,world); + MPI_Bcast(&kn,1,MPI_DOUBLE,0,world); + MPI_Bcast(&kt,1,MPI_DOUBLE,0,world); MPI_Bcast(&gamman,1,MPI_DOUBLE,0,world); + MPI_Bcast(&gammat,1,MPI_DOUBLE,0,world); MPI_Bcast(&xmu,1,MPI_DOUBLE,0,world); MPI_Bcast(&dampflag,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- */ -void *PairGranHookeHistory::extract(char *str) -{ - if (strcmp(str,"xkk") == 0) return (void *) &xkk; - else if (strcmp(str,"gamman") == 0) return (void *) &gamman; - else if (strcmp(str,"xmu") == 0) return (void *) &xmu; - else if (strcmp(str,"dampflag") == 0) return (void *) &dampflag; - return NULL; -} - -/* ---------------------------------------------------------------------- */ - void PairGranHookeHistory::reset_dt() { dt = update->dt; diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index f340575bda..e88e168b90 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -23,7 +23,7 @@ class PairGranHookeHistory : public Pair { PairGranHookeHistory(class LAMMPS *); ~PairGranHookeHistory(); virtual void compute(int, int); - void settings(int, char **); + virtual void settings(int, char **); void coeff(int, char **); void init_style(); void init_list(int, class NeighList *); @@ -32,13 +32,11 @@ class PairGranHookeHistory : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); - void *extract(char *); void reset_dt(); protected: - double xkk,xkkt,xmu; + double kn,kt,gamman,gammat,xmu; int dampflag; - double gamman,gammas; double dt; int freeze_group_bit; int history; diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index 2067c57998..61bed32869 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -143,7 +143,7 @@ void FixPeriNeigh::setup(int vflag) // scan neighbor list to set maxpartner - Pair *anypair = force->pair_match("peri"); + Pair *anypair = force->pair_match("peri",0); double **cutsq = anypair->cutsq; for (ii = 0; ii < inum; ii++) { diff --git a/src/force.cpp b/src/force.cpp index f5cb1aee2c..a1cb6402d9 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -152,23 +152,30 @@ Pair *Force::new_pair(const char *style) } /* ---------------------------------------------------------------------- - return ptr to current pair class or hybrid sub-class if it contains word + return ptr to current pair class or hybrid sub-class + if exact, then style name must be exact match to word + if not exact, style name must contain word else return NULL ------------------------------------------------------------------------- */ -Pair *Force::pair_match(const char *word) +Pair *Force::pair_match(const char *word, int exact) { - if (strstr(pair_style,word)) return pair; + if (exact && strcmp(pair_style,word) == 0) return pair; + else if (!exact && strstr(pair_style,word)) return pair; else if (strcmp(pair_style,"hybrid") == 0) { PairHybrid *hybrid = (PairHybrid *) pair; for (int i = 0; i < hybrid->nstyles; i++) { - if (strstr(hybrid->keywords[i],word)) + if (exact && strcmp(hybrid->keywords[i],word) == 0) + return hybrid->styles[i]; + else if (!exact && strstr(hybrid->keywords[i],word)) return hybrid->styles[i]; } } else if (strcmp(pair_style,"hybrid/overlay") == 0) { PairHybridOverlay *hybrid = (PairHybridOverlay *) pair; for (int i = 0; i < hybrid->nstyles; i++) { - if (strstr(hybrid->keywords[i],word)) + if (exact && strcmp(hybrid->keywords[i],word) == 0) + return hybrid->styles[i]; + else if (!exact && strstr(hybrid->keywords[i],word)) return hybrid->styles[i]; } } diff --git a/src/force.h b/src/force.h index 96a8ab8b4b..84948fd271 100644 --- a/src/force.h +++ b/src/force.h @@ -62,7 +62,7 @@ class Force : protected Pointers { void create_pair(const char *); class Pair *new_pair(const char *); - class Pair *pair_match(const char *); + class Pair *pair_match(const char *, int); void create_bond(const char *); class Bond *new_bond(const char *); diff --git a/src/pair.cpp b/src/pair.cpp index 9d4a2ea94d..1b34d55f9e 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -792,18 +792,18 @@ void Pair::write_file(int narg, char **arg) // if pair style = soft, set prefactor to final value - Pair *spair = force->pair_match("soft"); + Pair *spair = force->pair_match("soft",1); if (spair) ((PairSoft *) spair)->prefactor[itype][jtype] = ((PairSoft *) spair)->prestop[itype][jtype]; - // if pair style = EAM, swap in dummy fp vector + // if pair style = any of EAM, swap in dummy fp vector double eamfp[2]; eamfp[0] = eamfp[1] = 0.0; double *eamfp_hold; - Pair *epair = force->pair_match("eam"); + Pair *epair = force->pair_match("eam",0); if (epair) epair->swap_eam(eamfp,&eamfp_hold); // if atom style defines charge, swap in dummy q vec