diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index a418e9352d..1afecdf94d 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -276,3 +276,149 @@ void PairGranHertzHistory::settings(int narg, char **arg) kn /= force->nktv2p; kt /= force->nktv2p; } + +/* ---------------------------------------------------------------------- */ + +double PairGranHertzHistory::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double radi,radj,radsum; + double r,rinv,rsqinv,delx,dely,delz; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; + double meff,damp,ccel,polyhertz; + double vtr1,vtr2,vtr3,vrel,shrmag,rsht; + double fs1,fs2,fs3,fs,fn; + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + + if (rsq >= radsum*radsum) { + fforce = 0.0; + svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + return 0.0; + } + + r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; + + // relative translational velocity + + double **v = atom->v; + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + double **x = atom->x; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + + vnnr = vr1*delx + vr2*dely + vr3*delz; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + double **omega = atom->omega; + wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; + wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; + wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; + + // normal force = Hertzian contact + normal velocity damping + + double *rmass = atom->rmass; + double *mass = atom->mass; + int *mask = atom->mask; + + if (rmass) { + 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]; + } else { + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; + polyhertz = sqrt((radsum-r)*radi*radj / radsum); + ccel *= polyhertz; + + // relative velocities + + vtr1 = vt1 - (delz*wr2-dely*wr3); + vtr2 = vt2 - (delx*wr3-delz*wr1); + vtr3 = vt3 - (dely*wr1-delx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // shear history effects + // neighprev = index of found neigh on previous call + // search entire jnum list of neighbors of I for neighbor J + // start from neighprev, since will typically be next neighbor + // reset neighprev to 0 as necessary + + int *jlist = list->firstneigh[i]; + int jnum = list->numneigh[i]; + int *touch = list->listgranhistory->firstneigh[i]; + double *allshear = list->listgranhistory->firstdouble[i]; + + for (int jj = 0; jj < jnum; jj++) { + neighprev++; + if (neighprev >= jnum) neighprev = 0; + if (touch[neighprev] == j) break; + } + + double *shear = &allshear[3*neighprev]; + shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + + shear[2]*shear[2]); + + // rotate shear displacements + + rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; + rsht *= rsqinv; + + // tangential forces = shear + tangential velocity damping + + 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 + + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + fn = xmu * fabs(ccel*r); + + if (fs > fn) { + if (shrmag != 0.0) { + fs1 *= fn/fs; + fs2 *= fn/fs; + fs3 *= fn/fs; + fs *= fn/fs; + } else fs1 = fs2 = fs3 = fs = 0.0; + } + + // set all forces and return no energy + + fforce = ccel; + svector[0] = fs1; + svector[1] = fs2; + svector[2] = fs3; + svector[3] = fs; + return 0.0; +} diff --git a/src/GRANULAR/pair_gran_hertz_history.h b/src/GRANULAR/pair_gran_hertz_history.h index 5c5877a161..1b86354fb1 100644 --- a/src/GRANULAR/pair_gran_hertz_history.h +++ b/src/GRANULAR/pair_gran_hertz_history.h @@ -29,6 +29,7 @@ class PairGranHertzHistory : public PairGranHookeHistory { PairGranHertzHistory(class LAMMPS *); virtual void compute(int, int); void settings(int, char **); + double single(int, int, int, int, double, double, double, double &); }; } diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 95651e7606..4f2c00e270 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -191,3 +191,111 @@ void PairGranHooke::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } + +/* ---------------------------------------------------------------------- */ + +double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double radi,radj,radsum,r,rinv,rsqinv; + double delx,dely,delz; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; + double meff,damp,ccel; + double fn,fs,ft,fs1,fs2,fs3; + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + + // zero out forces if caller requests non-touching pair outside cutoff + + if (rsq >= radsum*radsum) { + fforce = 0.0; + svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + return 0.0; + } + + r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; + + // relative translational velocity + + double **v = atom->v; + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + double **x = atom->x; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + + vnnr = vr1*delx + vr2*dely + vr3*delz; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + double **omega = atom->omega; + wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; + wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; + wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; + + // normal forces = Hookian contact + normal velocity damping + + double *rmass = atom->rmass; + double *mass = atom->mass; + int *mask = atom->mask; + + if (rmass) { + 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]; + } else { + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; + + // relative velocities + + vtr1 = vt1 - (delz*wr2-dely*wr3); + vtr2 = vt2 - (delx*wr3-delz*wr1); + vtr3 = vt3 - (dely*wr1-delx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // force normalization + + fn = xmu * fabs(ccel*r); + fs = meff*gammat*vrel; + if (vrel != 0.0) ft = MIN(fn,fs) / vrel; + else ft = 0.0; + + // set all forces and return no energy + + fforce = ccel; + svector[0] = -ft*vtr1; + svector[1] = -ft*vtr2; + svector[2] = -ft*vtr3; + svector[3] = sqrt(svector[0]*svector[0] + + svector[1]*svector[1] + + svector[2]*svector[2]); + return 0.0; +} diff --git a/src/GRANULAR/pair_gran_hooke.h b/src/GRANULAR/pair_gran_hooke.h index d5cff796f2..20c2c3c1bd 100644 --- a/src/GRANULAR/pair_gran_hooke.h +++ b/src/GRANULAR/pair_gran_hooke.h @@ -28,6 +28,7 @@ class PairGranHooke : public PairGranHookeHistory { public: PairGranHooke(class LAMMPS *); virtual void compute(int, int); + double single(int, int, int, int, double, double, double, double &); }; } diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 847ed84a58..7e8cbc22fc 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -42,19 +42,24 @@ using namespace LAMMPS_NS; PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) { - single_enable = 0; + single_enable = 1; no_virial_fdotr_compute = 1; history = 1; fix_history = NULL; suffix = NULL; + single_extra = 4; + svector = new double[4]; + laststep = -1; + neighprev = 0; } /* ---------------------------------------------------------------------- */ PairGranHookeHistory::~PairGranHookeHistory() { + delete [] svector; if (fix_history) modify->delete_fix("SHEAR_HISTORY"); if (suffix) delete[] suffix; @@ -545,3 +550,147 @@ void PairGranHookeHistory::reset_dt() { dt = update->dt; } + +/* ---------------------------------------------------------------------- */ + +double PairGranHookeHistory::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double radi,radj,radsum; + double r,rinv,rsqinv,delx,dely,delz; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; + double meff,damp,ccel,polyhertz; + double vtr1,vtr2,vtr3,vrel,shrmag,rsht; + double fs1,fs2,fs3,fs,fn; + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + + if (rsq >= radsum*radsum) { + fforce = 0.0; + svector[0] = svector[1] = svector[2] = svector[3] = 0.0; + return 0.0; + } + + r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; + + // relative translational velocity + + double **v = atom->v; + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + double **x = atom->x; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + + vnnr = vr1*delx + vr2*dely + vr3*delz; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + double **omega = atom->omega; + wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; + wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; + wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; + + // normal force = Hertzian contact + normal velocity damping + + double *rmass = atom->rmass; + double *mass = atom->mass; + int *mask = atom->mask; + + if (rmass) { + 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]; + } else { + meff = mass[itype]*mass[jtype] / (mass[itype]+mass[jtype]); + if (mask[i] & freeze_group_bit) meff = mass[jtype]; + if (mask[j] & freeze_group_bit) meff = mass[itype]; + } + + damp = meff*gamman*vnnr*rsqinv; + ccel = kn*(radsum-r)*rinv - damp; + + // relative velocities + + vtr1 = vt1 - (delz*wr2-dely*wr3); + vtr2 = vt2 - (delx*wr3-delz*wr1); + vtr3 = vt3 - (dely*wr1-delx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // shear history effects + // neighprev = index of found neigh on previous call + // search entire jnum list of neighbors of I for neighbor J + // start from neighprev, since will typically be next neighbor + // reset neighprev to 0 as necessary + + int *jlist = list->firstneigh[i]; + int jnum = list->numneigh[i]; + int *touch = list->listgranhistory->firstneigh[i]; + double *allshear = list->listgranhistory->firstdouble[i]; + + for (int jj = 0; jj < jnum; jj++) { + neighprev++; + if (neighprev >= jnum) neighprev = 0; + if (touch[neighprev] == j) break; + } + + double *shear = &allshear[3*neighprev]; + shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + + shear[2]*shear[2]); + + // rotate shear displacements + + rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; + rsht *= rsqinv; + + // tangential forces = shear + tangential velocity damping + + 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 + + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + fn = xmu * fabs(ccel*r); + + if (fs > fn) { + if (shrmag != 0.0) { + fs1 *= fn/fs; + fs2 *= fn/fs; + fs3 *= fn/fs; + fs *= fn/fs; + } else fs1 = fs2 = fs3 = fs = 0.0; + } + + // set all forces and return no energy + + fforce = ccel; + svector[0] = fs1; + svector[1] = fs2; + svector[2] = fs3; + svector[3] = fs; + return 0.0; +} diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index e645c12bfe..2e88fc2d46 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -38,6 +38,7 @@ class PairGranHookeHistory : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + double single(int, int, int, int, double, double, double, double &); void reset_dt(); protected: @@ -55,6 +56,8 @@ class PairGranHookeHistory : public Pair { double *onerad_dynamic,*onerad_frozen; double *maxrad_dynamic,*maxrad_frozen; + int neighprev; + void allocate(); }; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 8c35f01702..af1c415e79 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -13,6 +13,7 @@ #include "math.h" #include "string.h" +#include "stdlib.h" #include "compute_pair_local.h" #include "atom.h" #include "update.h" @@ -29,6 +30,8 @@ using namespace LAMMPS_NS; #define DELTA 10000 +enum{DIST,ENG,FORCE,FX,FY,FZ,PN}; + /* ---------------------------------------------------------------------- */ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : @@ -41,18 +44,32 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : if (nvalues == 1) size_local_cols = 0; else size_local_cols = nvalues; - dflag = eflag = fflag = -1; - nvalues = 0; + pstyle = new int[nvalues]; + pindex = new int[nvalues]; - int i; + nvalues = 0; for (int iarg = 3; iarg < narg; iarg++) { - i = iarg-3; - if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++; - else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++; - else if (strcmp(arg[iarg],"force") == 0) fflag = nvalues++; - else error->all(FLERR,"Invalid keyword in compute pair/local command"); + if (strcmp(arg[iarg],"dist") == 0) pstyle[nvalues++] = DIST; + else if (strcmp(arg[iarg],"eng") == 0) pstyle[nvalues++] = ENG; + else if (strcmp(arg[iarg],"force") == 0) pstyle[nvalues++] = FORCE; + else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX; + else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY; + else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ; + else if (arg[iarg][0] == 'p') { + int n = atoi(&arg[iarg][1]); + if (n <= 0) error->all(FLERR, + "Invalid keyword in compute pair/local command"); + pstyle[nvalues] = PN; + pindex[nvalues++] = n-1; + } else error->all(FLERR,"Invalid keyword in compute pair/local command"); } + // set singleflag if need to call pair->single() + + singleflag = 0; + for (int i = 0; i < nvalues; i++) + if (pstyle[i] != DIST) singleflag = 1; + nmax = 0; vector = NULL; array = NULL; @@ -64,17 +81,24 @@ ComputePairLocal::~ComputePairLocal() { memory->destroy(vector); memory->destroy(array); + delete [] pstyle; + delete [] pindex; } /* ---------------------------------------------------------------------- */ void ComputePairLocal::init() { - if (force->pair == NULL) + if (singleflag && force->pair == NULL) error->all(FLERR,"No pair style is defined for compute pair/local"); - if (force->pair->single_enable == 0) + if (singleflag && force->pair->single_enable == 0) error->all(FLERR,"Pair style does not support compute pair/local"); + for (int i = 0; i < nvalues; i++) + if (pstyle[i] == PN && pindex[i] >= force->pair->single_extra) + error->all(FLERR,"Pair style does not have single field" + " requested by compute pair/local"); + // need an occasional half neighbor list int irequest = neighbor->request((void *) this); @@ -101,7 +125,7 @@ void ComputePairLocal::compute_local() ncount = compute_pairs(0); if (ncount > nmax) reallocate(ncount); size_local_rows = ncount; - ncount = compute_pairs(1); + compute_pairs(1); } /* ---------------------------------------------------------------------- @@ -117,7 +141,7 @@ int ComputePairLocal::compute_pairs(int flag) double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,eng,fpair,factor_coul,factor_lj; int *ilist,*jlist,*numneigh,**firstneigh; - double *dbuf,*ebuf,*fbuf; + double *ptr; double **x = atom->x; int *type = atom->type; @@ -138,23 +162,13 @@ int ComputePairLocal::compute_pairs(int flag) // loop over neighbors of my atoms // skip if I or J are not in group - - if (flag) { - if (nvalues == 1) { - if (dflag >= 0) dbuf = vector; - if (eflag >= 0) ebuf = vector; - if (fflag >= 0) fbuf = vector; - } else { - if (dflag >= 0) dbuf = &array[0][dflag]; - if (eflag >= 0) ebuf = &array[0][eflag]; - if (fflag >= 0) fbuf = &array[0][fflag]; - } - } + // for flag = 0, just count pair interactions within force cutoff + // for flag = 1, calculate requested output fields Pair *pair = force->pair; double **cutsq = force->pair->cutsq; - m = n = 0; + m = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (!(mask[i] & groupbit)) continue; @@ -183,13 +197,37 @@ int ComputePairLocal::compute_pairs(int flag) if (rsq >= cutsq[itype][jtype]) continue; if (flag) { - if (dflag >= 0) dbuf[n] = sqrt(rsq); - if (eflag >= 0 || fflag >= 0) { + if (singleflag) eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); - if (eflag >= 0) ebuf[n] = eng; - if (fflag >= 0) fbuf[n] = sqrt(rsq)*fpair; + + if (nvalues == 1) ptr = &vector[m]; + else ptr = array[m]; + + for (n = 0; n < nvalues; n++) { + switch (pstyle[n]) { + case DIST: + ptr[n] = sqrt(rsq); + break; + case ENG: + ptr[n] = eng; + break; + case FORCE: + ptr[n] = sqrt(rsq)*fpair; + break; + case FX: + ptr[n] = delx*fpair; + break; + case FY: + ptr[n] = dely*fpair; + break; + case FZ: + ptr[n] = delz*fpair; + break; + case PN: + ptr[n] = pair->svector[pindex[n]]; + break; + } } - n += nvalues; } m++; diff --git a/src/compute_pair_local.h b/src/compute_pair_local.h index e6706c193f..aaa8cdf4c9 100644 --- a/src/compute_pair_local.h +++ b/src/compute_pair_local.h @@ -37,6 +37,10 @@ class ComputePairLocal : public Compute { int nvalues,dflag,eflag,fflag; int ncount; + int *pstyle; // style of each requested output + int *pindex; // for pI, index of the output (0 to M-1) + int singleflag; + int nmax; double *vector; double **array; diff --git a/src/pair.cpp b/src/pair.cpp index 88e938aea2..3f60f94926 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -58,6 +58,8 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) nextra = 0; pvector = NULL; + single_extra = 0; + svector = NULL; // pair_modify settings diff --git a/src/pair.h b/src/pair.h index 6a88af5e07..9c0c74a1a4 100644 --- a/src/pair.h +++ b/src/pair.h @@ -53,6 +53,9 @@ class Pair : protected Pointers { int nextra; // # of extra quantities pair style calculates double *pvector; // vector of extra pair quantities + int single_extra; // number of extra single values calculated + double *svector; // vector of extra single quantities + class NeighList *list; // standard neighbor list used by most pairs class NeighList *listhalf; // half list used by some pairs class NeighList *listfull; // full list used by some pairs diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 206782be27..71a1366001 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -45,6 +45,8 @@ PairHybrid::~PairHybrid() delete [] keywords; } + delete [] svector; + if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -265,19 +267,32 @@ void PairHybrid::settings(int narg, char **arg) if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse); } - // single_enable = 0 if any sub-style = 0 + // single_enable = 1 if any sub-style is set // respa_enable = 1 if any sub-style is set // no_virial_fdotr_compute = 1 if any sub-style is set // ghostneigh = 1 if any sub-style is set + single_enable = 0; for (m = 0; m < nstyles; m++) - if (styles[m]->single_enable == 0) single_enable = 0; + if (styles[m]->single_enable) single_enable = 1; for (m = 0; m < nstyles; m++) if (styles[m]->respa_enable) respa_enable = 1; for (m = 0; m < nstyles; m++) if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1; for (m = 0; m < nstyles; m++) if (styles[m]->ghostneigh) ghostneigh = 1; + + // single_extra = min of all sub-style single_extra + // allocate svector + + single_extra = styles[0]->single_extra; + for (m = 1; m < nstyles; m++) + single_extra = MIN(single_extra,styles[m]->single_extra); + + if (single_extra) { + delete [] svector; + svector = new double[single_extra]; + } } /* ---------------------------------------------------------------------- @@ -613,10 +628,17 @@ double PairHybrid::single(int i, int j, int itype, int jtype, for (int m = 0; m < nmap[itype][jtype]; m++) { if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) { if (styles[map[itype][jtype][m]]->single_enable == 0) - error->all(FLERR,"Pair hybrid sub-style does not support single call"); + error->one(FLERR,"Pair hybrid sub-style does not support single call"); + esum += styles[map[itype][jtype][m]]-> single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone); fforce += fone; + + // copy substyle extra values into hybrid's svector + + if (single_extra && styles[map[itype][jtype][m]]->single_extra) + for (m = 0; m < single_extra; m++) + svector[m] = styles[map[itype][jtype][m]]->svector[m]; } } @@ -666,7 +688,8 @@ void *PairHybrid::extract(char *str, int &dim) double *p_newvalue = (double *) ptr; double newvalue = *p_newvalue; if (cutptr && newvalue != cutvalue) - error->all(FLERR,"Coulomb cutoffs of pair hybrid sub-styles do not match"); + error->all(FLERR, + "Coulomb cutoffs of pair hybrid sub-styles do not match"); cutptr = ptr; cutvalue = newvalue; } else if (ptr) return ptr; diff --git a/src/random_mars.cpp b/src/random_mars.cpp index 2430eabd56..0dd5990362 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -27,7 +27,7 @@ RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp) double s,t; if (seed <= 0 || seed > 900000000) - error->all(FLERR,"Invalid seed for Marsaglia random # generator"); + error->one(FLERR,"Invalid seed for Marsaglia random # generator"); save = 0; u = new double[97+1]; diff --git a/src/random_park.cpp b/src/random_park.cpp index 6867a61a56..3e53aee20a 100644 --- a/src/random_park.cpp +++ b/src/random_park.cpp @@ -40,7 +40,7 @@ using namespace LAMMPS_NS; RanPark::RanPark(LAMMPS *lmp, int seed_init) : Pointers(lmp) { if (seed_init <= 0) - error->all(FLERR,"Invalid seed for Park random # generator"); + error->one(FLERR,"Invalid seed for Park random # generator"); seed = seed_init; save = 0; }