diff --git a/src/ASPHERE/compute_erotate_asphere.cpp b/src/ASPHERE/compute_erotate_asphere.cpp index 612619f810..3a0322e17a 100644 --- a/src/ASPHERE/compute_erotate_asphere.cpp +++ b/src/ASPHERE/compute_erotate_asphere.cpp @@ -15,14 +15,13 @@ #include "compute_erotate_asphere.h" #include "math_extra.h" #include "atom.h" +#include "update.h" #include "force.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 - /* ---------------------------------------------------------------------- */ ComputeERotateAsphere::ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) : @@ -63,7 +62,7 @@ void ComputeERotateAsphere::init() double ComputeERotateAsphere::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **quat = atom->quat; double **angmom = atom->angmom; diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index ca40b40876..b12b35e704 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -19,6 +19,7 @@ #include "compute_temp_asphere.h" #include "math_extra.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "modify.h" @@ -29,9 +30,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) : @@ -150,11 +148,10 @@ void ComputeTempAsphere::dof_compute() double ComputeTempAsphere::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; if (tempbias) { - if (!(tbias->invoked & INVOKED_SCALAR)) - double tmp = tbias->compute_scalar(); + if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } @@ -213,10 +210,10 @@ void ComputeTempAsphere::compute_vector() { int i; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; if (tempbias) { - if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index b943712b01..de23f8d291 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -42,6 +42,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : { if (narg < 6) error->all("Illegal fix pour command"); + time_depend = 1; + if (!atom->radius_flag || !atom->rmass_flag) error->all("Fix pour requires atom attributes radius, rmass"); diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 1724a6b0ca..b371821ed4 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -45,6 +45,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : { if (narg < 4) error->all("Illegal fix wall/gran command"); + time_depend = 1; + if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag) error->all("Fix wall/gran requires atom attributes radius, omega, torque"); @@ -199,10 +201,8 @@ void FixWallGran::init() // same initialization as in pair_gran_history::init_style() xkkt = xkk * 2.0/7.0; - dt = update->dt; double gammas = 0.5*gamman; - gamman_dl = gamman/dt; - gammas_dl = gammas/dt; + dt = update->dt; // set pairstyle from granular pair style @@ -324,9 +324,11 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; - double fn,fs,ft,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3,rinv; + double fn,fs,ft,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3,rinv,rsqinv; r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -334,16 +336,12 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*dx + vr2*dy + vr3*dz; - vn1 = dx*vnnr / rsq; - vn2 = dy*vnnr / rsq; - vn3 = dz*vnnr / rsq; + vn1 = dx*vnnr * rsqinv; + vn2 = dy*vnnr * rsqinv; + vn3 = dz*vnnr * rsqinv; // tangential component @@ -353,20 +351,15 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, // relative rotational velocity - wr1 = radius*omega[0]; - wr2 = radius*omega[1]; - wr3 = radius*omega[2]; + wr1 = radius*omega[0] * rinv; + wr2 = radius*omega[1] * rinv; + wr3 = radius*omega[2] * rinv; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // normal forces = Hookian contact + normal velocity damping xmeff = mass; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radius-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radius-r)*rinv - damp; // relative velocities @@ -379,31 +372,26 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz, // force normalization fn = xmu * fabs(ccel*r); - fs = xmeff*gammas_dl*vrel; + fs = xmeff*gammas*vrel; if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; - // shear friction forces + // tangential force due to tangential velocity damping fs1 = -ft*vtr1; fs2 = -ft*vtr2; fs3 = -ft*vtr3; - // force components + // forces & torques - ccelx = dx*ccel + fs1; - ccely = dy*ccel + fs2; - ccelz = dz*ccel + fs3; + fx = dx*ccel + fs1; + fy = dy*ccel + fs2; + fz = dz*ccel + fs3; - // forces + f[0] += fx; + f[1] += fy; + f[2] += fz; - f[0] += ccelx; - f[1] += ccely; - f[2] += ccelz; - - // torques - - rinv = 1/r; tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -421,10 +409,12 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; - double fn,fs,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3; - double shrmag,rsht,rinv; + double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; + double shrmag,rsht,rinv,rsqinv; r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -432,16 +422,12 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*dx + vr2*dy + vr3*dz; - vn1 = dx*vnnr / rsq; - vn2 = dy*vnnr / rsq; - vn3 = dz*vnnr / rsq; + vn1 = dx*vnnr * rsqinv; + vn2 = dy*vnnr * rsqinv; + vn3 = dz*vnnr * rsqinv; // tangential component @@ -451,20 +437,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, // relative rotational velocity - wr1 = radius*omega[0]; - wr2 = radius*omega[1]; - wr3 = radius*omega[2]; + wr1 = radius*omega[0] * rinv; + wr2 = radius*omega[1] * rinv; + wr3 = radius*omega[2] * rinv; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // normal forces = Hookian contact + normal velocity damping xmeff = mass; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radius-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radius-r)*rinv - damp; // relative velocities @@ -476,60 +457,54 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz, // shear history effects - shear[0] += vtr1; - shear[1] += vtr2; - shear[2] += vtr3; + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - // rotate shear displacements correctly + // rotate shear displacements rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; - rsht = rsht/rsq; + rsht = rsht*rsqinv; shear[0] -= rsht*dx; shear[1] -= rsht*dy; shear[2] -= rsht*dz; - // tangential forces + // tangential forces = shear + tangential velocity damping - fs1 = - (xkkt*shear[0] + xmeff*gammas_dl*vtr1); - fs2 = - (xkkt*shear[1] + xmeff*gammas_dl*vtr2); - fs3 = - (xkkt*shear[2] + xmeff*gammas_dl*vtr3); + fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1); + fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2); + fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3); - // force normalization + // rescale frictional displacements and forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fn = xmu * fabs(ccel*r); - // shrmag is magnitude of shearwall - // rescale frictional displacements and forces if needed - if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) - - xmeff*gammas_dl*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) - - xmeff*gammas_dl*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) - - xmeff*gammas_dl*vtr3/xkkt; + 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; } else fs1 = fs2 = fs3 = 0.0; } - ccelx = dx*ccel + fs1; - ccely = dy*ccel + fs2; - ccelz = dz*ccel + fs3; + // forces & torques - // forces + fx = dx*ccel + fs1; + fy = dy*ccel + fs2; + fz = dz*ccel + fs3; - f[0] += ccelx; - f[1] += ccely; - f[2] += ccelz; + f[0] += fx; + f[1] += fy; + f[2] += fz; - // torques - - rinv = 1/r; tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -547,10 +522,12 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, { double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel; - double fn,fs,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3; - double shrmag,rsht,rhertz,rinv; + double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; + double shrmag,rsht,rhertz,rinv,rsqinv; r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -558,10 +535,6 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*dx + vr2*dy + vr3*dz; @@ -577,20 +550,15 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, // relative rotational velocity - wr1 = radius*omega[0]; - wr2 = radius*omega[1]; - wr3 = radius*omega[2]; + wr1 = radius*omega[0] * rinv; + wr2 = radius*omega[1] * rinv; + wr3 = radius*omega[2] * rinv; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // normal forces = Hertzian contact + normal velocity damping xmeff = mass; - damp = xmeff*gamman_dl*vnnr/rsq; - ccel = xkk*(radius-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radius-r)*rinv - damp; rhertz = sqrt(radius - r); ccel = rhertz * ccel; @@ -604,60 +572,54 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz, // shear history effects - shear[0] += vtr1; - shear[1] += vtr2; - shear[2] += vtr3; + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - // rotate shear displacements correctly + // rotate shear displacements rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; - rsht = rsht/rsq; + rsht = rsht*rsqinv; shear[0] -= rsht*dx; shear[1] -= rsht*dy; shear[2] -= rsht*dz; - // tangential forces + // tangential forces = shear + tangential velocity damping - fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas_dl*vtr1); - fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas_dl*vtr2); - fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas_dl*vtr3); + fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1); + fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2); + fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3); - // force normalization + // rescale frictional displacements and forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fn = xmu * fabs(ccel*r); - // shrmag is magnitude of shearwall - // rescale frictional displacements and forces if needed - if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) - - xmeff*gammas_dl*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) - - xmeff*gammas_dl*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) - - xmeff*gammas_dl*vtr3/xkkt; + 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; } else fs1 = fs2 = fs3 = 0.0; } - ccelx = dx*ccel + fs1; - ccely = dy*ccel + fs2; - ccelz = dz*ccel + fs3; + // forces & torques - // forces + fx = dx*ccel + fs1; + fy = dy*ccel + fs2; + fz = dz*ccel + fs3; - f[0] += ccelx; - f[1] += ccely; - f[2] += ccelz; + f[0] += fx; + f[1] += fy; + f[2] += fz; - // torques - - rinv = 1/r; tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -778,7 +740,4 @@ int FixWallGran::size_restart(int nlocal) void FixWallGran::reset_dt() { dt = update->dt; - double gammas = 0.5*gamman; - gamman_dl = gamman/dt; - gammas_dl = gammas/dt; } diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h index 8c8fb35403..11108def0c 100644 --- a/src/GRANULAR/fix_wall_gran.h +++ b/src/GRANULAR/fix_wall_gran.h @@ -40,9 +40,9 @@ class FixWallGran : public Fix { private: int wallstyle,pairstyle,wiggle,wshear,axis; - double xkk,xkkt,gamman,xmu; + double xkk,xkkt,gamman,gammas,xmu; double lo,hi,cylradius; - double dt,gamman_dl,gammas_dl; + double dt; double amplitude,period,omega,time_origin,vshear; int *touch; diff --git a/src/GRANULAR/pair_gran_hertzian.cpp b/src/GRANULAR/pair_gran_hertzian.cpp index fc52c19e64..94597fc512 100644 --- a/src/GRANULAR/pair_gran_hertzian.cpp +++ b/src/GRANULAR/pair_gran_hertzian.cpp @@ -42,7 +42,7 @@ void PairGranHertzian::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -108,6 +108,8 @@ void PairGranHertzian::compute(int eflag, int vflag) } else { r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -115,16 +117,12 @@ void PairGranHertzian::compute(int eflag, int vflag) vr2 = v[i][1] - v[j][1]; vr3 = v[i][2] - v[j][2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; // tangential component @@ -134,24 +132,19 @@ void PairGranHertzian::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; + 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; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // 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_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radsum-r)*rinv - damp; rhertz = sqrt(radsum - r); - ccel = rhertz * ccel; + ccel *= rhertz; // relative velocities @@ -162,31 +155,29 @@ void PairGranHertzian::compute(int eflag, int vflag) vrel = sqrt(vrel); // shear history effects - // shrmag = magnitude of shear touch[jj] = 1; shear = &allshear[3*jj]; - shear[0] += vtr1; - shear[1] += vtr2; - shear[2] += vtr3; + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - // rotate shear displacements correctly + // rotate shear displacements rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht /= rsq; + rsht *= rsqinv; shear[0] -= rsht*delx; shear[1] -= rsht*dely; shear[2] -= rsht*delz; - // tangential forces + // tangential forces = shear + tangential velocity damping - fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas_dl*vtr1); - fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas_dl*vtr2); - fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas_dl*vtr3); + fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1); + fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2); + fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3); - // force normalization // rescale frictional displacements and forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); @@ -194,12 +185,12 @@ void PairGranHertzian::compute(int eflag, int vflag) if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) - - xmeff*gammas_dl*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) - - xmeff*gammas_dl*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) - - xmeff*gammas_dl*vtr3/xkkt; + 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 *= fn/fs; fs2 *= fn/fs; fs3 *= fn/fs; @@ -219,7 +210,6 @@ void PairGranHertzian::compute(int eflag, int vflag) f[i][1] += fy; f[i][2] += fz; - rinv = 1/r; tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); diff --git a/src/GRANULAR/pair_gran_history.cpp b/src/GRANULAR/pair_gran_history.cpp index f0650f3238..d014469e35 100644 --- a/src/GRANULAR/pair_gran_history.cpp +++ b/src/GRANULAR/pair_gran_history.cpp @@ -69,7 +69,7 @@ void PairGranHistory::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -135,6 +135,8 @@ void PairGranHistory::compute(int eflag, int vflag) } else { r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -142,16 +144,12 @@ void PairGranHistory::compute(int eflag, int vflag) vr2 = v[i][1] - v[j][1]; vr3 = v[i][2] - v[j][2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; // tangential component @@ -161,22 +159,17 @@ void PairGranHistory::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; + 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; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // 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_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radsum-r)*rinv - damp; // relative velocities @@ -187,31 +180,29 @@ void PairGranHistory::compute(int eflag, int vflag) vrel = sqrt(vrel); // shear history effects - // shrmag = magnitude of shear touch[jj] = 1; shear = &allshear[3*jj]; - shear[0] += vtr1; - shear[1] += vtr2; - shear[2] += vtr3; + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - // rotate shear displacements correctly + // rotate shear displacements rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht /= rsq; + rsht *= rsqinv; shear[0] -= rsht*delx; shear[1] -= rsht*dely; shear[2] -= rsht*delz; - // tangential forces + // tangential forces = shear + tangential velocity damping - fs1 = - (xkkt*shear[0] + xmeff*gammas_dl*vtr1); - fs2 = - (xkkt*shear[1] + xmeff*gammas_dl*vtr2); - fs3 = - (xkkt*shear[2] + xmeff*gammas_dl*vtr3); + fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1); + fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2); + fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3); - // force normalization // rescale frictional displacements and forces if needed fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); @@ -219,12 +210,12 @@ void PairGranHistory::compute(int eflag, int vflag) if (fs > fn) { if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) - - xmeff*gammas_dl*vtr1/xkkt; - shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) - - xmeff*gammas_dl*vtr2/xkkt; - shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) - - xmeff*gammas_dl*vtr3/xkkt; + 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 *= fn/fs; fs2 *= fn/fs; fs3 *= fn/fs; @@ -244,7 +235,6 @@ void PairGranHistory::compute(int eflag, int vflag) f[i][1] += fy; f[i][2] += fz; - rinv = 1/r; tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); @@ -297,15 +287,6 @@ void PairGranHistory::settings(int narg, char **arg) gamman = atof(arg[1]); xmu = atof(arg[2]); dampflag = atoi(arg[3]); - - // granular styles do not use pair_coeff, so set setflag for everything now - - if (!allocated) allocate(); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - setflag[i][j] = 1; } /* ---------------------------------------------------------------------- @@ -314,7 +295,22 @@ void PairGranHistory::settings(int narg, char **arg) void PairGranHistory::coeff(int narg, char **arg) { - error->all("Granular pair styles do not use pair_coeff settings"); + if (narg > 2) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -344,11 +340,9 @@ void PairGranHistory::init_style() } xkkt = xkk * 2.0/7.0; - dt = update->dt; double gammas = 0.5*gamman; if (dampflag == 0) gammas = 0.0; - gamman_dl = gamman/dt; - gammas_dl = gammas/dt; + dt = update->dt; // if shear history is stored: // check if newton flag is valid @@ -512,8 +506,4 @@ void *PairGranHistory::extract(char *str) void PairGranHistory::reset_dt() { dt = update->dt; - double gammas = 0.5*gamman; - if (dampflag == 0) gammas = 0.0; - gamman_dl = gamman/dt; - gammas_dl = gammas/dt; } diff --git a/src/GRANULAR/pair_gran_history.h b/src/GRANULAR/pair_gran_history.h index 9509d6636f..cccbc161ce 100644 --- a/src/GRANULAR/pair_gran_history.h +++ b/src/GRANULAR/pair_gran_history.h @@ -38,8 +38,8 @@ class PairGranHistory : public Pair { protected: double xkk,xkkt,xmu; int dampflag; - double gamman; - double dt,gamman_dl,gammas_dl; + double gamman,gammas; + double dt; int freeze_group_bit; int history; class FixShearHistory *fix_history; diff --git a/src/GRANULAR/pair_gran_no_history.cpp b/src/GRANULAR/pair_gran_no_history.cpp index 54fb62a594..dbeedcf045 100644 --- a/src/GRANULAR/pair_gran_no_history.cpp +++ b/src/GRANULAR/pair_gran_no_history.cpp @@ -41,7 +41,7 @@ void PairGranNoHistory::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -91,6 +91,8 @@ void PairGranNoHistory::compute(int eflag, int vflag) if (rsq < radsum*radsum) { r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; // relative translational velocity @@ -98,16 +100,12 @@ void PairGranNoHistory::compute(int eflag, int vflag) vr2 = v[i][1] - v[j][1]; vr3 = v[i][2] - v[j][2]; - vr1 *= dt; - vr2 *= dt; - vr3 *= dt; - // normal component vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr / rsq; - vn2 = dely*vnnr / rsq; - vn3 = delz*vnnr / rsq; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; // tangential component @@ -117,22 +115,17 @@ void PairGranNoHistory::compute(int eflag, int vflag) // relative rotational velocity - wr1 = radi*omega[i][0] + radj*omega[j][0]; - wr2 = radi*omega[i][1] + radj*omega[j][1]; - wr3 = radi*omega[i][2] + radj*omega[j][2]; + 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; - wr1 *= dt/r; - wr2 *= dt/r; - wr3 *= dt/r; - - // normal damping term - // this definition of DAMP includes the extra 1/r term + // 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_dl*vnnr/rsq; - ccel = xkk*(radsum-r)/r - damp; + damp = xmeff*gamman*vnnr*rsqinv; + ccel = xkk*(radsum-r)*rinv - damp; // relative velocities @@ -145,11 +138,11 @@ void PairGranNoHistory::compute(int eflag, int vflag) // force normalization fn = xmu * fabs(ccel*r); - fs = xmeff*gammas_dl*vrel; + fs = xmeff*gammas*vrel; if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; - // shear friction forces + // tangential force due to tangential velocity damping fs1 = -ft*vtr1; fs2 = -ft*vtr2; @@ -164,7 +157,6 @@ void PairGranNoHistory::compute(int eflag, int vflag) f[i][1] += fy; f[i][2] += fz; - rinv = 1/r; tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); diff --git a/src/MOLECULE/Install.csh b/src/MOLECULE/Install.csh index 5a5ec8ebee..7442989420 100644 --- a/src/MOLECULE/Install.csh +++ b/src/MOLECULE/Install.csh @@ -34,6 +34,8 @@ if ($1 == 1) then cp dihedral_multi_harmonic.cpp .. cp dihedral_opls.cpp .. cp dump_bond.cpp .. + cp fix_bond_break.cpp .. + cp fix_bond_create.cpp .. cp fix_bond_swap.cpp .. cp improper.cpp .. cp improper_cvff.cpp .. @@ -69,6 +71,8 @@ if ($1 == 1) then cp dihedral_multi_harmonic.h .. cp dihedral_opls.h .. cp dump_bond.h .. + cp fix_bond_break.h .. + cp fix_bond_create.h .. cp fix_bond_swap.h .. # cp improper.h .. cp improper_cvff.h .. @@ -109,6 +113,8 @@ else if ($1 == 0) then rm ../dihedral_multi_harmonic.cpp rm ../dihedral_opls.cpp rm ../dump_bond.cpp + rm ../fix_bond_break.cpp + rm ../fix_bond_create.cpp rm ../fix_bond_swap.cpp rm ../improper.cpp rm ../improper_cvff.cpp @@ -144,6 +150,8 @@ else if ($1 == 0) then rm ../dihedral_multi_harmonic.h rm ../dihedral_opls.h rm ../dump_bond.h + rm ../fix_bond_break.h + rm ../fix_bond_create.h rm ../fix_bond_swap.h # rm ../improper.h rm ../improper_cvff.h diff --git a/src/MOLECULE/angle_hybrid.cpp b/src/MOLECULE/angle_hybrid.cpp index 23cd6ad468..04bb495a62 100644 --- a/src/MOLECULE/angle_hybrid.cpp +++ b/src/MOLECULE/angle_hybrid.cpp @@ -103,7 +103,6 @@ void AngleHybrid::compute(int eflag, int vflag) else evflag = 0; for (m = 0; m < nstyles; m++) { - if (styles[m] == NULL) continue; neighbor->nanglelist = nanglelist[m]; neighbor->anglelist = anglelist[m]; @@ -168,6 +167,8 @@ void AngleHybrid::settings(int narg, char **arg) error->all("Angle style hybrid cannot use same angle style twice"); if (strcmp(arg[m],"hybrid") == 0) error->all("Angle style hybrid cannot have hybrid as an argument"); + if (strcmp(arg[m],"none") == 0) + error->all("Angle style hybrid cannot have none as an argument"); styles[m] = force->new_angle(arg[m]); keywords[m] = new char[strlen(arg[m])+1]; strcpy(keywords[m],arg[m]); @@ -185,12 +186,18 @@ void AngleHybrid::coeff(int which, int narg, char **arg) int ilo,ihi; force->bounds(arg[0],atom->nangletypes,ilo,ihi); - // 2nd arg = angle style name (harmonic, etc) + // 2nd arg = angle sub-style name + // allow for "none" as valid sub-style name int m; for (m = 0; m < nstyles; m++) if (strcmp(arg[1],keywords[m]) == 0) break; - if (m == nstyles) error->all("Angle coeff for hybrid has invalid style"); + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[1],"none") == 0) none = 1; + else error->all("Angle coeff for hybrid has invalid style"); + } // move 1st arg to 2nd arg // just copy ptrs, since arg[] points into original input line @@ -199,15 +206,19 @@ void AngleHybrid::coeff(int which, int narg, char **arg) // invoke sub-style coeff() starting with 1st arg - if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]); + if (!none) styles[m]->coeff(which,narg-1,&arg[1]); // set setflag and which type maps to which sub-style - // if sub-style is NULL for "none", still set setflag + // if sub-style is none: set hybrid setflag, wipe out map for (int i = ilo; i <= ihi; i++) { - map[i] = m; - if (styles[m] == NULL) setflag[i] = 1; - else setflag[i] = styles[m]->setflag[i]; + if (none) { + setflag[i] = 1; + map[i] = -1; + } else { + setflag[i] = styles[m]->setflag[i]; + map[i] = m; + } } } @@ -217,6 +228,7 @@ void AngleHybrid::coeff(int which, int narg, char **arg) double AngleHybrid::equilibrium_angle(int i) { + if (map[i] < 0) error->one("Invoked angle equil angle on angle style none"); return styles[map[i]]->equilibrium_angle(i); } @@ -265,8 +277,8 @@ void AngleHybrid::read_restart(FILE *fp) double AngleHybrid::single(int type, int i1, int i2, int i3) { - if (styles[map[type]]) return styles[map[type]]->single(type,i1,i2,i3); - return 0.0; + if (map[type] < 0) error->one("Invoked angle single on angle style none"); + return styles[map[type]]->single(type,i1,i2,i3); } /* ---------------------------------------------------------------------- diff --git a/src/MOLECULE/bond_hybrid.cpp b/src/MOLECULE/bond_hybrid.cpp index 06642f42c2..7076f2e9d7 100644 --- a/src/MOLECULE/bond_hybrid.cpp +++ b/src/MOLECULE/bond_hybrid.cpp @@ -102,7 +102,6 @@ void BondHybrid::compute(int eflag, int vflag) else evflag = 0; for (m = 0; m < nstyles; m++) { - if (styles[m] == NULL) continue; neighbor->nbondlist = nbondlist[m]; neighbor->bondlist = bondlist[m]; @@ -167,6 +166,8 @@ void BondHybrid::settings(int narg, char **arg) error->all("Bond style hybrid cannot use same bond style twice"); if (strcmp(arg[m],"hybrid") == 0) error->all("Bond style hybrid cannot have hybrid as an argument"); + if (strcmp(arg[m],"none") == 0) + error->all("Bond style hybrid cannot have none as an argument"); styles[m] = force->new_bond(arg[m]); keywords[m] = new char[strlen(arg[m])+1]; strcpy(keywords[m],arg[m]); @@ -184,12 +185,18 @@ void BondHybrid::coeff(int narg, char **arg) int ilo,ihi; force->bounds(arg[0],atom->nbondtypes,ilo,ihi); - // 2nd arg = bond style name (harmonic, fene, etc) + // 2nd arg = bond sub-style name + // allow for "none" as valid sub-style name int m; for (m = 0; m < nstyles; m++) if (strcmp(arg[1],keywords[m]) == 0) break; - if (m == nstyles) error->all("Bond coeff for hybrid has invalid style"); + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[1],"none") == 0) none = 1; + else error->all("Bond coeff for hybrid has invalid style"); + } // move 1st arg to 2nd arg // just copy ptrs, since arg[] points into original input line @@ -198,14 +205,15 @@ void BondHybrid::coeff(int narg, char **arg) // invoke sub-style coeff() starting with 1st arg - if (styles[m]) styles[m]->coeff(narg-1,&arg[1]); + if (!none) styles[m]->coeff(narg-1,&arg[1]); // set setflag and which type maps to which sub-style - // if sub-style is NULL for "none", still set setflag + // if sub-style is none: set hybrid setflag, wipe out map for (int i = ilo; i <= ihi; i++) { - map[i] = m; setflag[i] = 1; + if (none) map[i] = -1; + else map[i] = m; } } @@ -223,6 +231,7 @@ void BondHybrid::init_style() double BondHybrid::equilibrium_distance(int i) { + if (map[i] < 0) error->one("Invoked bond equil distance on bond style none"); return styles[map[i]]->equilibrium_distance(i); } @@ -271,8 +280,8 @@ void BondHybrid::read_restart(FILE *fp) double BondHybrid::single(int type, double rsq, int i, int j) { - if (styles[map[type]]) return styles[map[type]]->single(type,rsq,i,j); - return 0.0; + if (map[type] < 0) error->one("Invoked bond single on bond style none"); + return styles[map[type]]->single(type,rsq,i,j); } /* ---------------------------------------------------------------------- diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index 75130c8383..d540d54a01 100755 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -245,7 +245,7 @@ void BondQuartic::init_style() if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) - error->all("Must use special bonds = 1,1,1 with bond style quartic"); + error->all("Bond style quartic requires special_bonds = 1,1,1"); } /* ---------------------------------------------------------------------- diff --git a/src/MOLECULE/dihedral_hybrid.cpp b/src/MOLECULE/dihedral_hybrid.cpp index a956180ac1..60fda615bd 100644 --- a/src/MOLECULE/dihedral_hybrid.cpp +++ b/src/MOLECULE/dihedral_hybrid.cpp @@ -105,7 +105,6 @@ void DihedralHybrid::compute(int eflag, int vflag) else evflag = 0; for (m = 0; m < nstyles; m++) { - if (styles[m] == NULL) continue; neighbor->ndihedrallist = ndihedrallist[m]; neighbor->dihedrallist = dihedrallist[m]; @@ -170,6 +169,8 @@ void DihedralHybrid::settings(int narg, char **arg) error->all("Dihedral style hybrid cannot use same dihedral style twice"); if (strcmp(arg[m],"hybrid") == 0) error->all("Dihedral style hybrid cannot have hybrid as an argument"); + if (strcmp(arg[m],"none") == 0) + error->all("Dihedral style hybrid cannot have none as an argument"); styles[m] = force->new_dihedral(arg[m]); keywords[m] = new char[strlen(arg[m])+1]; strcpy(keywords[m],arg[m]); @@ -187,12 +188,18 @@ void DihedralHybrid::coeff(int which, int narg, char **arg) int ilo,ihi; force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi); - // 2nd arg = dihedral style name (harmonic, etc) + // 2nd arg = dihedral sub-style name + // allow for "none" as valid sub-style name int m; for (m = 0; m < nstyles; m++) if (strcmp(arg[1],keywords[m]) == 0) break; - if (m == nstyles) error->all("Dihedral coeff for hybrid has invalid style"); + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[1],"none") == 0) none = 1; + else error->all("Dihedral coeff for hybrid has invalid style"); + } // move 1st arg to 2nd arg // just copy ptrs, since arg[] points into original input line @@ -201,15 +208,19 @@ void DihedralHybrid::coeff(int which, int narg, char **arg) // invoke sub-style coeff() starting with 1st arg - if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]); + if (!none) styles[m]->coeff(which,narg-1,&arg[1]); // set setflag and which type maps to which sub-style - // if sub-style is NULL for "none", still set setflag + // if sub-style is none: set hybrid setflag, wipe out map for (int i = ilo; i <= ihi; i++) { - map[i] = m; - if (styles[m] == NULL) setflag[i] = 1; - else setflag[i] = styles[m]->setflag[i]; + if (none) { + setflag[i] = 1; + map[i] = -1; + } else { + setflag[i] = styles[m]->setflag[i]; + map[i] = m; + } } } diff --git a/src/MOLECULE/fix_bond_swap.cpp b/src/MOLECULE/fix_bond_swap.cpp index 36bd71de53..b42ebe32b3 100644 --- a/src/MOLECULE/fix_bond_swap.cpp +++ b/src/MOLECULE/fix_bond_swap.cpp @@ -77,6 +77,8 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) : nmax = 0; alist = NULL; + + naccept = foursome = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/MOLECULE/improper_hybrid.cpp b/src/MOLECULE/improper_hybrid.cpp index 46fb9adf63..4a836d7df2 100644 --- a/src/MOLECULE/improper_hybrid.cpp +++ b/src/MOLECULE/improper_hybrid.cpp @@ -105,7 +105,6 @@ void ImproperHybrid::compute(int eflag, int vflag) else evflag = 0; for (m = 0; m < nstyles; m++) { - if (styles[m] == NULL) continue; neighbor->nimproperlist = nimproperlist[m]; neighbor->improperlist = improperlist[m]; @@ -170,6 +169,8 @@ void ImproperHybrid::settings(int narg, char **arg) error->all("Improper style hybrid cannot use same improper style twice"); if (strcmp(arg[m],"hybrid") == 0) error->all("Improper style hybrid cannot have hybrid as an argument"); + if (strcmp(arg[m],"none") == 0) + error->all("Improper style hybrid cannot have none as an argument"); styles[m] = force->new_improper(arg[m]); keywords[m] = new char[strlen(arg[m])+1]; strcpy(keywords[m],arg[m]); @@ -187,12 +188,18 @@ void ImproperHybrid::coeff(int which, int narg, char **arg) int ilo,ihi; force->bounds(arg[0],atom->nimpropertypes,ilo,ihi); - // 2nd arg = improper style name (harmonic, etc) + // 2nd arg = improper sub-style name + // allow for "none" as valid sub-style name int m; for (m = 0; m < nstyles; m++) if (strcmp(arg[1],keywords[m]) == 0) break; - if (m == nstyles) error->all("Improper coeff for hybrid has invalid style"); + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[1],"none") == 0) none = 1; + else error->all("Improper coeff for hybrid has invalid style"); + } // move 1st arg to 2nd arg // just copy ptrs, since arg[] points into original input line @@ -201,15 +208,19 @@ void ImproperHybrid::coeff(int which, int narg, char **arg) // invoke sub-style coeff() starting with 1st arg - if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]); + if (!none) styles[m]->coeff(which,narg-1,&arg[1]); // set setflag and which type maps to which sub-style - // if sub-style is NULL for "none", still set setflag + // if sub-style is none: set hybrid setflag, wipe out map for (int i = ilo; i <= ihi; i++) { - map[i] = m; - if (styles[m] == NULL) setflag[i] = 1; - else setflag[i] = styles[m]->setflag[i]; + if (none) { + setflag[i] = 1; + map[i] = -1; + } else { + setflag[i] = styles[m]->setflag[i]; + map[i] = m; + } } } diff --git a/src/MOLECULE/style_molecule.h b/src/MOLECULE/style_molecule.h index 4ee90a09f3..68e98d6764 100644 --- a/src/MOLECULE/style_molecule.h +++ b/src/MOLECULE/style_molecule.h @@ -90,10 +90,14 @@ DumpStyle(bond,DumpBond) #endif #ifdef FixInclude +#include "fix_bond_break.h" +#include "fix_bond_create.h" #include "fix_bond_swap.h" #endif #ifdef FixClass +FixStyle(bond/break,FixBondBreak) +FixStyle(bond/create,FixBondCreate) FixStyle(bond/swap,FixBondSwap) #endif diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp index c573c97529..8df2ec603a 100644 --- a/src/PERI/compute_damage_atom.cpp +++ b/src/PERI/compute_damage_atom.cpp @@ -18,6 +18,7 @@ #include "string.h" #include "compute_damage_atom.h" #include "atom.h" +#include "update.h" #include "modify.h" #include "comm.h" #include "force.h" @@ -72,6 +73,8 @@ void ComputeDamageAtom::init() void ComputeDamageAtom::compute_peratom() { + invoked_peratom = update->ntimestep; + // grow damage array if necessary if (atom->nlocal > nmax) { diff --git a/src/USER-ACKLAND/compute_ackland_atom.cpp b/src/USER-ACKLAND/compute_ackland_atom.cpp index 80db2b40b9..10a9af1868 100644 --- a/src/USER-ACKLAND/compute_ackland_atom.cpp +++ b/src/USER-ACKLAND/compute_ackland_atom.cpp @@ -19,8 +19,8 @@ #include "string.h" #include "compute_ackland_atom.h" #include "atom.h" -#include "modify.h" #include "update.h" +#include "modify.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" @@ -34,7 +34,6 @@ using namespace LAMMPS_NS; enum{UNKNOWN,BCC,FCC,HCP,ICO}; -#define INVOKED_PERATOM 4 /* ---------------------------------------------------------------------- */ @@ -102,7 +101,7 @@ void ComputeAcklandAtom::compute_peratom() int *ilist,*jlist,*numneigh,**firstneigh; int chi[8]; - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; // grow structure array if necessary diff --git a/src/atom.cpp b/src/atom.cpp index 7cc2115625..9e1b39e044 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -51,6 +51,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0; nbonds = nangles = ndihedrals = nimpropers = 0; bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; + extra_bond_per_atom = 0; firstgroupname = NULL; diff --git a/src/atom.h b/src/atom.h index 3690678841..f4cf022b01 100644 --- a/src/atom.h +++ b/src/atom.h @@ -34,6 +34,7 @@ class Atom : protected Pointers { int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes; int nbonds,nangles,ndihedrals,nimpropers; int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom; + int extra_bond_per_atom; int firstgroup; // store atoms in this group first, -1 if unset int nfirst; // # of atoms in first group on this proc @@ -51,8 +52,9 @@ class Atom : protected Pointers { double *radius,*density,*rmass,*vfrac,*s0; double **x0; - int maxspecial; - int **nspecial,**special; + int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs + int **special; // IDs of 1-2,1-3,1-4 neighs of each atom + int maxspecial; // special[nlocal][maxspecial] int *num_bond; int **bond_type; diff --git a/src/compute.cpp b/src/compute.cpp index 6472ca3ccd..4cd76b95d3 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -66,16 +66,17 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) id_pre = NULL; timeflag = 0; - invoked = 0; comm_forward = comm_reverse = 0; + invoked_scalar = invoked_vector = invoked_peratom = -1; + // set modify defaults extra_dof = domain->dimension; dynamic = 0; // setup list of timesteps - + ntime = maxtime = 0; tlist = NULL; } @@ -175,3 +176,12 @@ int Compute::matchstep(int ntimestep) } return 0; } + +/* ---------------------------------------------------------------------- + clean out list of timesteps to call the compute on +------------------------------------------------------------------------- */ + +void Compute::clearstep() +{ + ntime = 0; +} diff --git a/src/compute.h b/src/compute.h index 1c1cb5490b..233e61f9ea 100644 --- a/src/compute.h +++ b/src/compute.h @@ -54,7 +54,10 @@ class Compute : protected Pointers { int maxtime; // max # of entries time list can hold int *tlist; // time list of steps the Compute is called on - int invoked; // set when Compute is invoked, to avoid re-invoking + int invoked_flag; // 1 if invoked or accessed this step, 0 if not + int invoked_scalar; // last timestep on which compute_scalar() was invoked + int invoked_vector; // ditto for compute_vector() + int invoked_peratom; // ditto for compute_peratom() double dof; // degrees-of-freedom for temperature @@ -85,6 +88,7 @@ class Compute : protected Pointers { void addstep(int); int matchstep(int); + void clearstep(); virtual double memory_usage() {return 0.0;} diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp index cf4b2d63be..dc51ba7ad1 100644 --- a/src/compute_centro_atom.cpp +++ b/src/compute_centro_atom.cpp @@ -14,8 +14,8 @@ #include "string.h" #include "compute_centro_atom.h" #include "atom.h" -#include "modify.h" #include "update.h" +#include "modify.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" @@ -27,8 +27,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) : @@ -59,6 +57,9 @@ ComputeCentroAtom::~ComputeCentroAtom() void ComputeCentroAtom::init() { + if (force->pair == NULL) + error->all("Compute centro/atom requires a pair style be defined"); + int count = 0; for (int i = 0; i < modify->ncompute; i++) if (strcmp(modify->compute[i]->style,"centro/atom") == 0) count++; @@ -91,7 +92,7 @@ void ComputeCentroAtom::compute_peratom() int *ilist,*jlist,*numneigh,**firstneigh; double pairs[66]; - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; // grow centro array if necessary diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index f9e398f453..f121f04042 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -15,6 +15,7 @@ #include "stdlib.h" #include "compute_coord_atom.h" #include "atom.h" +#include "update.h" #include "modify.h" #include "neighbor.h" #include "neigh_list.h" @@ -27,8 +28,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : @@ -90,7 +89,7 @@ void ComputeCoordAtom::compute_peratom() double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *ilist,*jlist,*numneigh,**firstneigh; - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; // grow coordination array if necessary diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp index 447b065342..62b41cbfa8 100644 --- a/src/compute_displace_atom.cpp +++ b/src/compute_displace_atom.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_displace_atom.h" #include "atom.h" +#include "update.h" #include "domain.h" #include "modify.h" #include "fix.h" @@ -23,8 +24,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) : @@ -75,7 +74,7 @@ void ComputeDisplaceAtom::init() void ComputeDisplaceAtom::compute_peratom() { - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; // grow local displacement array if necessary diff --git a/src/compute_erotate_sphere.cpp b/src/compute_erotate_sphere.cpp index 98b0e1591f..91fb7a3e9f 100644 --- a/src/compute_erotate_sphere.cpp +++ b/src/compute_erotate_sphere.cpp @@ -14,6 +14,7 @@ #include "mpi.h" #include "compute_erotate_sphere.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "group.h" @@ -21,8 +22,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 - #define INERTIA 0.4 // moment of inertia for sphere /* ---------------------------------------------------------------------- */ @@ -75,7 +74,7 @@ void ComputeERotateSphere::init() double ComputeERotateSphere::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **omega = atom->omega; double *radius = atom->radius; diff --git a/src/compute_group_group.cpp b/src/compute_group_group.cpp index f5d42f09ce..79a6f63b90 100644 --- a/src/compute_group_group.cpp +++ b/src/compute_group_group.cpp @@ -18,6 +18,7 @@ #include "mpi.h" #include "compute_group_group.h" #include "atom.h" +#include "update.h" #include "force.h" #include "pair.h" #include "neighbor.h" @@ -28,9 +29,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) : @@ -62,7 +60,7 @@ ComputeGroupGroup::~ComputeGroupGroup() void ComputeGroupGroup::init() { if (force->pair == NULL) - error->all("Compute group/group requires pair style be defined"); + error->all("Compute group/group requires a pair style be defined"); pair = force->pair; cutsq = force->pair->cutsq; @@ -85,8 +83,7 @@ void ComputeGroupGroup::init_list(int id, NeighList *ptr) double ComputeGroupGroup::compute_scalar() { - invoked |= INVOKED_SCALAR; - invoked |= INVOKED_VECTOR; + invoked_scalar = invoked_vector = update->ntimestep; interact(); return scalar; @@ -96,8 +93,7 @@ double ComputeGroupGroup::compute_scalar() void ComputeGroupGroup::compute_vector() { - invoked |= INVOKED_SCALAR; - invoked |= INVOKED_VECTOR; + invoked_scalar = invoked_vector = update->ntimestep; interact(); } diff --git a/src/compute_ke.cpp b/src/compute_ke.cpp index c9e35d3fcb..7b77783c68 100644 --- a/src/compute_ke.cpp +++ b/src/compute_ke.cpp @@ -14,6 +14,7 @@ #include "mpi.h" #include "compute_ke.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "group.h" @@ -21,8 +22,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 - /* ---------------------------------------------------------------------- */ ComputeKE::ComputeKE(LAMMPS *lmp, int narg, char **arg) : @@ -45,7 +44,7 @@ void ComputeKE::init() double ComputeKE::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **v = atom->v; double *rmass = atom->rmass; diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp index e2417a4e6d..0596e68787 100644 --- a/src/compute_ke_atom.cpp +++ b/src/compute_ke_atom.cpp @@ -14,6 +14,7 @@ #include "string.h" #include "compute_ke_atom.h" #include "atom.h" +#include "update.h" #include "modify.h" #include "comm.h" #include "force.h" @@ -22,8 +23,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputeKEAtom::ComputeKEAtom(LAMMPS *lmp, int narg, char **arg) : @@ -60,7 +59,7 @@ void ComputeKEAtom::init() void ComputeKEAtom::compute_peratom() { - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; // grow ke array if necessary diff --git a/src/compute_pe.cpp b/src/compute_pe.cpp index 4788785dc5..694b06dd33 100644 --- a/src/compute_pe.cpp +++ b/src/compute_pe.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_pe.h" #include "atom.h" +#include "update.h" #include "force.h" #include "pair.h" #include "bond.h" @@ -28,8 +29,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 - /* ---------------------------------------------------------------------- */ ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) : @@ -73,7 +72,9 @@ ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) : double ComputePE::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; + if (update->eflag_global != invoked_scalar) + error->all("Energy was not tallied on needed timestep"); double one = 0.0; if (pairflag && force->pair) diff --git a/src/compute_pe_atom.cpp b/src/compute_pe_atom.cpp index 20fd75318c..31ac624480 100755 --- a/src/compute_pe_atom.cpp +++ b/src/compute_pe_atom.cpp @@ -14,6 +14,7 @@ #include "string.h" #include "compute_pe_atom.h" #include "atom.h" +#include "update.h" #include "comm.h" #include "force.h" #include "pair.h" @@ -26,8 +27,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) : @@ -76,7 +75,9 @@ void ComputePEAtom::compute_peratom() { int i; - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; + if (update->eflag_atom != invoked_peratom) + error->all("Per-atom energy was not tallied on needed timestep"); // grow local energy array if necessary diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp index 495e2b561a..9e1e490133 100644 --- a/src/compute_pressure.cpp +++ b/src/compute_pressure.cpp @@ -16,6 +16,7 @@ #include "stdlib.h" #include "compute_pressure.h" #include "atom.h" +#include "update.h" #include "domain.h" #include "modify.h" #include "fix.h" @@ -30,9 +31,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) : @@ -159,13 +157,16 @@ void ComputePressure::init() double ComputePressure::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; + if (update->vflag_global != invoked_scalar) + error->all("Virial was not tallied on needed timestep"); // invoke temperature it it hasn't been already double t; if (keflag) { - if (temperature->invoked & INVOKED_SCALAR) t = temperature->scalar; + if (temperature->invoked_scalar == update->ntimestep) + t = temperature->scalar; else t = temperature->compute_scalar(); } @@ -197,13 +198,15 @@ double ComputePressure::compute_scalar() void ComputePressure::compute_vector() { - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; + if (update->vflag_global != invoked_vector) + error->all("Virial was not tallied on needed timestep"); // invoke temperature it it hasn't been already double *ke_tensor; if (keflag) { - if (!(temperature->invoked & INVOKED_VECTOR)) + if (temperature->invoked_vector != update->ntimestep) temperature->compute_vector(); ke_tensor = temperature->vector; } diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index f34b2cc4fe..f2b725177d 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -27,10 +27,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 -#define INVOKED_PERATOM 4 - enum{SUM,MINN,MAXX}; enum{X,V,F,COMPUTE,FIX,VARIABLE}; @@ -224,7 +220,7 @@ void ComputeReduce::init() double ComputeReduce::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double one = compute_one(0); if (mode == SUM) @@ -240,7 +236,7 @@ double ComputeReduce::compute_scalar() void ComputeReduce::compute_vector() { - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; for (int m = 0; m < nvalues; m++) onevec[m] = compute_one(m); if (mode == SUM) @@ -288,7 +284,7 @@ double ComputeReduce::compute_one(int m) // invoke compute if not previously invoked } else if (which[m] == COMPUTE) { - if (!(modify->compute[n]->invoked & INVOKED_PERATOM)) + if (modify->compute[n]->invoked_peratom != update->ntimestep) modify->compute[n]->compute_peratom(); if (j == 0) { diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index 3f1703704b..9019bdfec5 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_stress_atom.h" #include "atom.h" +#include "update.h" #include "comm.h" #include "force.h" #include "pair.h" @@ -29,8 +30,6 @@ using namespace LAMMPS_NS; -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) : @@ -85,7 +84,9 @@ void ComputeStressAtom::compute_peratom() { int i,j; - invoked |= INVOKED_PERATOM; + invoked_peratom = update->ntimestep; + if (update->vflag_atom != invoked_peratom) + error->all("Per-atom virial was not tallied on needed timestep"); // grow local stress array if necessary diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index c43b7e868a..01c8327354 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_temp.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "modify.h" @@ -24,9 +25,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) : @@ -75,7 +73,7 @@ void ComputeTemp::dof_compute() double ComputeTemp::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; @@ -109,7 +107,7 @@ void ComputeTemp::compute_vector() { int i; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index 7ded6c967b..fbbce90680 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -16,6 +16,7 @@ #include "string.h" #include "compute_temp_com.h" #include "atom.h" +#include "update.h" #include "force.h" #include "group.h" #include "modify.h" @@ -29,9 +30,6 @@ using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) : @@ -86,7 +84,7 @@ double ComputeTempCOM::compute_scalar() { double vthermal[3]; - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); @@ -125,7 +123,7 @@ void ComputeTempCOM::compute_vector() int i; double vthermal[3]; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index ad496ee13c..8fab755cce 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -20,6 +20,7 @@ #include "compute_temp_deform.h" #include "domain.h" #include "atom.h" +#include "update.h" #include "force.h" #include "modify.h" #include "fix.h" @@ -31,8 +32,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp /* ---------------------------------------------------------------------- */ @@ -103,7 +102,7 @@ double ComputeTempDeform::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; @@ -152,7 +151,7 @@ void ComputeTempDeform::compute_vector() { double lamda[3],vstream[3],vthermal[3]; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 079b969b3c..490a5356d1 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -15,6 +15,7 @@ #include "stdlib.h" #include "compute_temp_partial.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "modify.h" @@ -25,9 +26,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : @@ -95,7 +93,7 @@ int ComputeTempPartial::dof_remove(int i) double ComputeTempPartial::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; @@ -130,7 +128,7 @@ void ComputeTempPartial::compute_vector() { int i; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index fe79202310..395bcb81e8 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -16,6 +16,7 @@ #include "string.h" #include "compute_temp_ramp.h" #include "atom.h" +#include "update.h" #include "force.h" #include "group.h" #include "modify.h" @@ -30,9 +31,6 @@ using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : @@ -151,7 +149,7 @@ double ComputeTempRamp::compute_scalar() { double fraction,vramp,vthermal[3]; - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; @@ -193,7 +191,7 @@ void ComputeTempRamp::compute_vector() int i; double fraction,vramp,vthermal[3]; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index 69475e5bbc..f2e06be3d6 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_temp_region.h" #include "atom.h" +#include "update.h" #include "force.h" #include "modify.h" #include "domain.h" @@ -24,9 +25,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) : @@ -77,7 +75,7 @@ int ComputeTempRegion::dof_remove(int i) double ComputeTempRegion::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; @@ -122,7 +120,7 @@ void ComputeTempRegion::compute_vector() { int i; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index cbd5f1bd08..d338e78fb6 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -15,6 +15,7 @@ #include "string.h" #include "compute_temp_sphere.h" #include "atom.h" +#include "update.h" #include "force.h" #include "domain.h" #include "modify.h" @@ -24,9 +25,6 @@ using namespace LAMMPS_NS; -#define INVOKED_SCALAR 1 -#define INVOKED_VECTOR 2 - #define INERTIA 0.4 // moment of inertia for sphere /* ---------------------------------------------------------------------- */ @@ -136,11 +134,10 @@ void ComputeTempSphere::dof_compute() double ComputeTempSphere::compute_scalar() { - invoked |= INVOKED_SCALAR; + invoked_scalar = update->ntimestep; if (tempbias) { - if (!(tbias->invoked & INVOKED_SCALAR)) - double tmp = tbias->compute_scalar(); + if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } @@ -186,10 +183,10 @@ void ComputeTempSphere::compute_vector() { int i; - invoked |= INVOKED_VECTOR; + invoked_vector = update->ntimestep; if (tempbias) { - if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector(); + if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index a567a6cfd2..8c2668f37a 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -201,7 +201,10 @@ void DeleteAtoms::delete_overlap(int narg, char **arg) lmp->init(); // error check on cutoff + // if no pair style, neighbor list will be empty + if (force->pair == NULL) + error->all("Delete_atoms requires a pair style be defined"); if (cut > neighbor->cutneighmax) error->all("Delete_atoms cutoff > neighbor cutoff"); diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 02f644e5a4..43d1f9f035 100644 --- a/src/delete_bonds.cpp +++ b/src/delete_bonds.cpp @@ -21,7 +21,6 @@ #include "neighbor.h" #include "comm.h" #include "force.h" -#include "pair.h" #include "group.h" #include "special.h" #include "error.h" diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 775a52a0a9..6d04a3933f 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -32,7 +32,7 @@ using namespace LAMMPS_NS; // customize by adding keyword to 1st enum -enum{TAG,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ, +enum{ID,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ, VX,VY,VZ,FX,FY,FZ, Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ, @@ -40,8 +40,6 @@ enum{TAG,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ, enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE}; -#define INVOKED_PERATOM 4 // same as in computes - /* ---------------------------------------------------------------------- */ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) : @@ -271,10 +269,14 @@ int DumpCustom::count() // invoke Computes for per-atom dump quantities // only if not already invoked - if (ncompute) - for (i = 0; i < ncompute; i++) - if (!(compute[i]->invoked & INVOKED_PERATOM)) + if (ncompute) { + int ntimestep = update->ntimestep; + for (i = 0; i < ncompute; i++) { + if (compute[i]->invoked_peratom != ntimestep) compute[i]->compute_peratom(); + compute[i]->invoked_flag = 1; + } + } // evaluate atom-style Variables for per-atom dump quantities @@ -322,7 +324,7 @@ int DumpCustom::count() // customize by adding to if statement - if (thresh_array[ithresh] == TAG) { + if (thresh_array[ithresh] == ID) { int *tag = atom->tag; for (i = 0; i < nlocal; i++) dchoose[i] = tag[i]; ptr = dchoose; @@ -658,8 +660,8 @@ void DumpCustom::parse_fields(int narg, char **arg) for (int iarg = 5; iarg < narg; iarg++) { i = iarg-5; - if (strcmp(arg[iarg],"tag") == 0) { - pack_choice[i] = &DumpCustom::pack_tag; + if (strcmp(arg[iarg],"id") == 0) { + pack_choice[i] = &DumpCustom::pack_id; vtype[i] = INT; } else if (strcmp(arg[iarg],"mol") == 0) { if (!atom->molecule_flag) @@ -1051,7 +1053,7 @@ int DumpCustom::modify_param(int narg, char **arg) // set keyword type of threshhold // customize by adding to if statement - if (strcmp(arg[1],"tag") == 0) thresh_array[nthresh] = TAG; + if (strcmp(arg[1],"id") == 0) thresh_array[nthresh] = ID; else if (strcmp(arg[1],"mol") == 0) thresh_array[nthresh] = MOL; else if (strcmp(arg[1],"type") == 0) thresh_array[nthresh] = TYPE; else if (strcmp(arg[1],"mass") == 0) thresh_array[nthresh] = MASS; @@ -1314,7 +1316,7 @@ void DumpCustom::pack_variable(int n) /* ---------------------------------------------------------------------- */ -void DumpCustom::pack_tag(int n) +void DumpCustom::pack_id(int n) { int *tag = atom->tag; int nlocal = atom->nlocal; diff --git a/src/dump_custom.h b/src/dump_custom.h index 3d3243a88a..da0fbfcc1a 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -88,7 +88,7 @@ class DumpCustom : public Dump { typedef void (DumpCustom::*FnPtrPack)(int); FnPtrPack *pack_choice; // ptrs to pack functions - void pack_tag(int); + void pack_id(int); void pack_molecule(int); void pack_type(int); void pack_mass(int); diff --git a/src/fix.cpp b/src/fix.cpp index 2fe254298c..f58077b509 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -53,6 +53,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) virial_flag = 0; no_change_box = 0; time_integrate = 0; + time_depend = 0; restart_pbc = 0; scalar_flag = vector_flag = peratom_flag = 0; @@ -65,17 +66,21 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // mask settings - same as in modify.cpp INITIAL_INTEGRATE = 1; - PRE_EXCHANGE = 2; - PRE_NEIGHBOR = 4; - POST_FORCE = 8; - FINAL_INTEGRATE = 16; - END_OF_STEP = 32; - THERMO_ENERGY = 64; - INITIAL_INTEGRATE_RESPA = 128; - POST_FORCE_RESPA = 256; - FINAL_INTEGRATE_RESPA = 512; - MIN_POST_FORCE = 1024; - MIN_ENERGY = 2048; + POST_INTEGRATE = 2; + PRE_EXCHANGE = 4; + PRE_NEIGHBOR = 8; + PRE_FORCE = 16; + POST_FORCE = 32; + FINAL_INTEGRATE = 64; + END_OF_STEP = 128; + THERMO_ENERGY = 256; + INITIAL_INTEGRATE_RESPA = 512; + POST_INTEGRATE_RESPA = 1024; + PRE_FORCE_RESPA = 2048; + POST_FORCE_RESPA = 4096; + FINAL_INTEGRATE_RESPA = 8192; + MIN_POST_FORCE = 16384; + MIN_ENERGY = 32768; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix.h b/src/fix.h index f7b62b6a90..b0a675a1fc 100644 --- a/src/fix.h +++ b/src/fix.h @@ -34,13 +34,14 @@ class Fix : protected Pointers { int virial_flag; // 1 if Fix contributes to virial, 0 if not int no_change_box; // 1 if cannot swap ortho <-> triclinic int time_integrate; // 1 if fix performs time integration, 0 if no + int time_depend; // 1 if fix is timestep dependent, 0 if not int restart_pbc; // 1 if fix moves atoms (except integrate) // so that write_restart must remap to PBC int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int size_vector; // N = size of global vector - int scalar_vector_freq; // frequency compute s/v data is available at + int scalar_vector_freq; // frequency s/v data is available at int extscalar; // 0/1 if scalar is intensive/extensive int extvector; // 0/1/-1 if vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component @@ -57,9 +58,11 @@ class Fix : protected Pointers { double virial[6]; // accumlated virial double **vatom; // accumulated per-atom virial - int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR; // mask settings - int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY; - int INITIAL_INTEGRATE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA; + int INITIAL_INTEGRATE,POST_INTEGRATE; // mask settings + int PRE_EXCHANGE,PRE_NEIGHBOR; + int PRE_FORCE,POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY; + int INITIAL_INTEGRATE_RESPA,POST_INTEGRATE_RESPA; + int PRE_FORCE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA; int MIN_POST_FORCE,MIN_ENERGY; Fix(class LAMMPS *, int, char **); @@ -73,8 +76,10 @@ class Fix : protected Pointers { virtual void setup(int) {} virtual void min_setup(int) {} virtual void initial_integrate(int) {} + virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} + virtual void pre_force(int) {} virtual void post_force(int) {} virtual void final_integrate() {} virtual void end_of_step() {} @@ -91,6 +96,8 @@ class Fix : protected Pointers { virtual int maxsize_restart() {return 0;} virtual void initial_integrate_respa(int, int, int) {} + virtual void post_integrate_respa(int, int) {} + virtual void pre_force_respa(int, int, int) {} virtual void post_force_respa(int, int, int) {} virtual void final_integrate_respa(int) {} diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp index a5a1389f8e..596cf38acb 100644 --- a/src/fix_ave_atom.cpp +++ b/src/fix_ave_atom.cpp @@ -27,8 +27,6 @@ using namespace LAMMPS_NS; enum{X,V,F,COMPUTE,FIX,VARIABLE}; -#define INVOKED_PERATOM 4 // same as in computes - /* ---------------------------------------------------------------------- */ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : @@ -36,6 +34,8 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : { if (narg < 7) error->all("Illegal fix ave/atom command"); + time_depend = 1; + nevery = atoi(arg[3]); nrepeat = atoi(arg[4]); peratom_freq = atoi(arg[5]); @@ -167,6 +167,7 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : atom->add_callback(0); // zero the array since dump may access it on timestep 0 + // zero the array since a variable may access it before first run int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) @@ -185,9 +186,9 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : nvalid -= (nrepeat-1)*nevery; if (nvalid < update->ntimestep) nvalid += peratom_freq; - // add nvalid to ALL computes that store invocation times + // add nvalid to all computes that store invocation times // since don't know a priori which are invoked by this fix - // once in end_of_step() can just set timestep for ones actually invoked + // once in end_of_step() can set timestep for ones actually invoked modify->addstep_compute_all(nvalid); } @@ -268,7 +269,8 @@ void FixAveAtom::end_of_step() // skip if not step which requires doing something - if (update->ntimestep != nvalid) return; + int ntimestep = update->ntimestep; + if (ntimestep != nvalid) return; // zero if first step @@ -307,7 +309,8 @@ void FixAveAtom::end_of_step() } else if (which[m] == COMPUTE) { Compute *compute = modify->compute[n]; - if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom(); + if (compute->invoked_peratom != ntimestep) compute->compute_peratom(); + compute->invoked_flag = 1; if (j == 0) { double *compute_scalar = compute->scalar_atom; @@ -351,7 +354,7 @@ void FixAveAtom::end_of_step() } irepeat = 0; - nvalid = update->ntimestep+peratom_freq - (nrepeat-1)*nevery; + nvalid = ntimestep+peratom_freq - (nrepeat-1)*nevery; modify->addstep_compute(nvalid); // average the final result for the Nfreq timestep diff --git a/src/fix_ave_force.cpp b/src/fix_ave_force.cpp index 242e81326c..d1176c657e 100644 --- a/src/fix_ave_force.cpp +++ b/src/fix_ave_force.cpp @@ -41,6 +41,9 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) : else yvalue = atof(arg[4]); if (strcmp(arg[5],"NULL") == 0) zflag = 0; else zvalue = atof(arg[5]); + + foriginal_all[0] = foriginal_all[1] = + foriginal_all[2] = foriginal_all[3] = 0.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index eef3178ee6..753520361c 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -38,7 +38,6 @@ enum{BOX,LATTICE,REDUCED}; enum{ONE,RUNNING,WINDOW}; #define BIG 1000000000 -#define INVOKED_PERATOM 4 // same as in computes /* ---------------------------------------------------------------------- */ @@ -50,6 +49,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : MPI_Comm_rank(world,&me); no_change_box = 1; + time_depend = 1; nevery = atoi(arg[3]); nrepeat = atoi(arg[4]); @@ -150,7 +150,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : // optional args normflag = ALL; - scaleflag = BOX; + scaleflag = LATTICE; fp = NULL; ave = ONE; @@ -315,9 +315,9 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : nvalid -= (nrepeat-1)*nevery; if (nvalid < update->ntimestep) nvalid += nfreq; - // add nvalid to ALL computes that store invocation times + // add nvalid to all computes that store invocation times // since don't know a priori which are invoked by this fix - // once in end_of_step() can just set timestep for ones actually invoked + // once in end_of_step() can set timestep for ones actually invoked modify->addstep_compute_all(nvalid); } @@ -417,7 +417,8 @@ void FixAveSpatial::end_of_step() // skip if not step which requires doing something - if (update->ntimestep != nvalid) return; + int ntimestep = update->ntimestep; + if (ntimestep != nvalid) return; // if computing the first sample, setup layers // compute current origin = boundary for some layer @@ -621,7 +622,8 @@ void FixAveSpatial::end_of_step() } else if (which[m] == COMPUTE) { Compute *compute = modify->compute[n]; - if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom(); + if (compute->invoked_peratom != ntimestep) compute->compute_peratom(); + compute->invoked_flag = 1; double *scalar = compute->scalar_atom; double **vector = compute->vector_atom; int jm1 = j - 1; @@ -699,7 +701,7 @@ void FixAveSpatial::end_of_step() } irepeat = 0; - nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery; + nvalid = ntimestep+nfreq - (nrepeat-1)*nevery; modify->addstep_compute(nvalid); // time average across samples @@ -782,7 +784,7 @@ void FixAveSpatial::end_of_step() // output result to file if (fp && me == 0) { - fprintf(fp,"%d %d\n",update->ntimestep,nlayers); + fprintf(fp,"%d %d\n",ntimestep,nlayers); for (m = 0; m < nlayers; m++) { fprintf(fp," %d %g %g",m+1,coord[m],count_total[m]/norm); for (i = 0; i < nvalues; i++) fprintf(fp," %g",values_total[m][i]/norm); @@ -803,6 +805,7 @@ double FixAveSpatial::compute_vector(int n) int ivalue = n % nvalues; int ilayer = n / nvalues; if (ilayer >= nlayers) return 0.0; + if (values_total == NULL) return 0.0; return values_total[ilayer][ivalue]/norm; } diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 5712b1b205..e5eaa35e2e 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -32,9 +32,6 @@ using namespace LAMMPS_NS; enum{COMPUTE,FIX,VARIABLE}; enum{ONE,RUNNING,WINDOW}; -#define INVOKED_SCALAR 1 // same as in computes -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : @@ -42,6 +39,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : { if (narg < 7) error->all("Illegal fix ave/time command"); + time_depend = 1; + MPI_Comm_rank(world,&me); nevery = atoi(arg[3]); @@ -248,9 +247,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : nvalid -= (nrepeat-1)*nevery; if (nvalid < update->ntimestep) nvalid += nfreq; - // add nvalid to ALL computes that store invocation times + // add nvalid to all computes that store invocation times // since don't know a priori which are invoked by this fix - // once in end_of_step() can just set timestep for ones actually invoked + // once in end_of_step() can set timestep for ones actually invoked modify->addstep_compute_all(nvalid); } @@ -330,7 +329,8 @@ void FixAveTime::end_of_step() // skip if not step which requires doing something - if (update->ntimestep != nvalid) return; + int ntimestep = update->ntimestep; + if (ntimestep != nvalid) return; // zero if first step @@ -351,11 +351,13 @@ void FixAveTime::end_of_step() Compute *compute = modify->compute[m]; if (argindex[i] == 0) { - if (compute->invoked & INVOKED_SCALAR) vector[i] += compute->scalar; + if (compute->invoked_scalar == ntimestep) vector[i] += compute->scalar; else vector[i] += compute->compute_scalar(); + compute->invoked_flag = 1; } else { - if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector(); + if (compute->invoked_vector != ntimestep) compute->compute_vector(); vector[i] += compute->vector[argindex[i]-1]; + compute->invoked_flag = 1; } // access fix fields, guaranteed to be ready @@ -383,7 +385,7 @@ void FixAveTime::end_of_step() } irepeat = 0; - nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery; + nvalid = ntimestep+nfreq - (nrepeat-1)*nevery; modify->addstep_compute(nvalid); // average the final result for the Nfreq timestep @@ -395,7 +397,7 @@ void FixAveTime::end_of_step() // if ave = RUNNING, combine with all previous Nfreq timestep values // if ave = WINDOW, combine with nwindow most recent Nfreq timestep values - if (update->ntimestep >= startstep) { + if (ntimestep >= startstep) { if (ave == ONE) { for (i = 0; i < nvalues; i++) vector_total[i] = vector[i]; norm = 1; @@ -424,7 +426,7 @@ void FixAveTime::end_of_step() // output result to file if (fp && me == 0) { - fprintf(fp,"%d",update->ntimestep); + fprintf(fp,"%d",ntimestep); for (i = 0; i < nvalues; i++) fprintf(fp," %g",vector_total[i]/norm); fprintf(fp,"\n"); fflush(fp); diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 7aaec977b4..c62328b9f3 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -46,6 +46,7 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) box_change = 1; no_change_box = 1; + time_depend = 1; nevery = atoi(arg[3]); if (nevery <= 0) error->all("Illegal fix deform command"); diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index 64f656e9d4..8ff2c49092 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -37,6 +37,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : if (narg < 7) error->all("Illegal fix deposit command"); restart_global = 1; + time_depend = 1; // required args diff --git a/src/fix_dt_reset.cpp b/src/fix_dt_reset.cpp index 358e770699..0fe6f0f101 100644 --- a/src/fix_dt_reset.cpp +++ b/src/fix_dt_reset.cpp @@ -43,6 +43,7 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) : { if (narg < 7) error->all("Illegal fix dt/reset command"); + time_depend = 1; scalar_flag = 1; vector_flag = 1; size_vector = 1; diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 576b1d14b3..9dda46ea30 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -32,6 +32,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : { if (narg < 5) error->all("Illegal fix gravity command"); + time_depend = 1; + magnitude = atof(arg[3]); if (strcmp(arg[4],"chute") == 0) { diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index ae9672388f..fd45bd60d2 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -38,6 +38,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : { if (narg < 4) error->all("Illegal fix indent command"); + time_depend = 1; scalar_flag = 1; vector_flag = 1; size_vector = 3; diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 44f1ce2b14..a02a6737ab 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -40,6 +40,8 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : { if (narg < 7) error->all("Illegal fix langevin command"); + time_depend = 1; + t_start = atof(arg[3]); t_stop = atof(arg[4]); t_period = atof(arg[5]); diff --git a/src/fix_nve_limit.cpp b/src/fix_nve_limit.cpp index 0fd132f26b..3cbd46eb30 100644 --- a/src/fix_nve_limit.cpp +++ b/src/fix_nve_limit.cpp @@ -37,6 +37,8 @@ FixNVELimit::FixNVELimit(LAMMPS *lmp, int narg, char **arg) : extscalar = 1; xlimit = atof(arg[3]); + + ncount = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 6a806c1f6c..ed7e616d75 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -62,8 +62,6 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) : if (t_period <= 0.0) error->all("Fix nvt period must be > 0.0"); t_freq = 1.0 / t_period; - eta = eta_dot = 0.0; - // create a new compute temp style // id = fix-ID + temp, compute group = fix group @@ -86,6 +84,10 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(3,newarg); delete [] newarg; tflag = 1; + + // Nose/Hoover temp init + + eta = eta_dot = 0.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_recenter.cpp b/src/fix_recenter.cpp index e2f1406f86..e5f243597d 100644 --- a/src/fix_recenter.cpp +++ b/src/fix_recenter.cpp @@ -28,6 +28,8 @@ using namespace LAMMPS_NS; +enum{BOX,LATTICE,FRACTION}; + /* ---------------------------------------------------------------------- */ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) : @@ -52,7 +54,7 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) : // optional args group2bit = groupbit; - scaleflag = 1; + scaleflag = LATTICE; int iarg = 6; while (iarg < narg) { @@ -62,9 +64,9 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) : group2bit = group->bitmask[igroup2]; iarg += 2; } else if (strcmp(arg[iarg],"units") == 0) { - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = 2; + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE; + else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = FRACTION; else error->all("Illegal fix recenter command"); iarg += 2; } else error->all("Illegal fix recenter command"); @@ -72,11 +74,11 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) : // scale xcom,ycom,zcom - if (scaleflag == 1 && domain->lattice == NULL) + if (scaleflag == LATTICE && domain->lattice == NULL) error->all("Use of fix recenter with undefined lattice"); double xscale,yscale,zscale; - if (scaleflag == 1) { + if (scaleflag == LATTICE) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; @@ -140,7 +142,7 @@ void FixRecenter::initial_integrate(int vflag) double xtarget,ytarget,ztarget; double *bboxlo,*bboxhi; - if (scaleflag == 2) { + if (scaleflag == FRACTION) { if (domain->triclinic == 0) { bboxlo = domain->boxlo; bboxhi = domain->boxhi; @@ -151,15 +153,18 @@ void FixRecenter::initial_integrate(int vflag) } if (xinitflag) xtarget = xinit; - else if (scaleflag == 2) xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]); + else if (scaleflag == FRACTION) + xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]); else xtarget = xcom; if (yinitflag) ytarget = yinit; - else if (scaleflag == 2) ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]); + else if (scaleflag == FRACTION) + ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]); else ytarget = ycom; if (zinitflag) ztarget = zinit; - else if (scaleflag == 2) ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]); + else if (scaleflag == FRACTION) + ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]); else ztarget = zcom; // current COM diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp index ecc2c78609..a55c828daf 100644 --- a/src/fix_shake.cpp +++ b/src/fix_shake.cpp @@ -64,7 +64,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) : grow_arrays(atom->nmax); atom->add_callback(0); - // set comm size needed by this Fix + // set comm size needed by this fix comm_forward = 3; diff --git a/src/fix_spring_self.cpp b/src/fix_spring_self.cpp index 65c2d41b11..2a196d958a 100644 --- a/src/fix_spring_self.cpp +++ b/src/fix_spring_self.cpp @@ -72,6 +72,8 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) : xoriginal[i][2] = x[i][2] + zbox*zprd; } else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; } + + espring = 0.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp index c857bae82c..4b580ac928 100644 --- a/src/fix_wall_lj126.cpp +++ b/src/fix_wall_lj126.cpp @@ -35,6 +35,7 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) : { if (narg < 8) error->all("Illegal fix wall/lj126 command"); + time_depend = 1; scalar_flag = 1; vector_flag = 1; size_vector = 3; diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp index 833785f958..27161da6fc 100644 --- a/src/fix_wall_lj93.cpp +++ b/src/fix_wall_lj93.cpp @@ -31,6 +31,7 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) : { if (narg < 8) error->all("Illegal fix wall/lj93 command"); + time_depend = 1; scalar_flag = 1; vector_flag = 1; size_vector = 3; diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index a1aa470f9b..ef6fc33770 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -16,7 +16,6 @@ #include "atom.h" #include "modify.h" #include "domain.h" -#include "comm.h" #include "error.h" using namespace LAMMPS_NS; @@ -52,29 +51,14 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : int FixWallReflect::setmask() { int mask = 0; - mask |= INITIAL_INTEGRATE; + mask |= POST_INTEGRATE; + mask |= POST_INTEGRATE_RESPA; return mask; } /* ---------------------------------------------------------------------- */ -void FixWallReflect::init() -{ - // warn if any integrate fix comes after this one - - int after = 0; - int flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(id,modify->fix[i]->id) == 0) after = 1; - else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1; - } - if (flag && comm->me == 0) - error->warning("Fix wall/reflect should come after all other integration fixes"); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallReflect::initial_integrate(int vflag) +void FixWallReflect::post_integrate() { double xlo = domain->boxlo[0]; double xhi = domain->boxhi[0]; diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h index 53c588ad72..294f0672b2 100644 --- a/src/fix_wall_reflect.h +++ b/src/fix_wall_reflect.h @@ -22,8 +22,7 @@ class FixWallReflect : public Fix { public: FixWallReflect(class LAMMPS *, int, char **); int setmask(); - void init(); - void initial_integrate(int); + void post_integrate(); private: int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag; diff --git a/src/fix_wiggle.cpp b/src/fix_wiggle.cpp index 6904a564b0..93037ca8aa 100644 --- a/src/fix_wiggle.cpp +++ b/src/fix_wiggle.cpp @@ -30,6 +30,8 @@ FixWiggle::FixWiggle(LAMMPS *lmp, int narg, char **arg) : { if (narg != 6) error->all("Illegal fix wiggle command"); + time_depend = 1; + // parse command-line args if (strcmp(arg[3],"x") == 0) axis = 0; diff --git a/src/force.cpp b/src/force.cpp index 400840b143..f5cb1aee2c 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -56,6 +56,7 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp) special_lj[1] = special_lj[2] = special_lj[3] = 0.0; special_coul[1] = special_coul[2] = special_coul[3] = 0.0; special_dihedral = 0; + special_extra = 0; dielectric = 1.0; @@ -353,54 +354,72 @@ void Force::set_special(int narg, char **arg) { if (narg == 0) error->all("Illegal special_bonds command"); - if (strcmp(arg[0],"charmm") == 0) { - if (narg != 1) error->all("Illegal special_bonds command"); - special_lj[1] = 0.0; - special_lj[2] = 0.0; - special_lj[3] = 0.0; - special_coul[1] = 0.0; - special_coul[2] = 0.0; - special_coul[3] = 0.0; - special_dihedral = 0; - return; + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"amber") == 0) { + if (iarg+1 > narg) error->all("Illegal special_bonds command"); + special_lj[1] = 0.0; + special_lj[2] = 0.0; + special_lj[3] = 0.5; + special_coul[1] = 0.0; + special_coul[2] = 0.0; + special_coul[3] = 5.0/6.0; + iarg += 1; + } else if (strcmp(arg[iarg],"charmm") == 0) { + if (iarg+1 > narg) error->all("Illegal special_bonds command"); + special_lj[1] = 0.0; + special_lj[2] = 0.0; + special_lj[3] = 0.0; + special_coul[1] = 0.0; + special_coul[2] = 0.0; + special_coul[3] = 0.0; + iarg += 1; + } else if (strcmp(arg[iarg],"fene") == 0) { + if (iarg+1 > narg) error->all("Illegal special_bonds command"); + special_lj[1] = 0.0; + special_lj[2] = 1.0; + special_lj[3] = 1.0; + special_coul[1] = 0.0; + special_coul[2] = 1.0; + special_coul[3] = 1.0; + iarg += 1; + } else if (strcmp(arg[iarg],"lj/coul") == 0) { + if (iarg+4 > narg) error->all("Illegal special_bonds command"); + special_lj[1] = special_coul[1] = atof(arg[iarg+1]); + special_lj[2] = special_coul[2] = atof(arg[iarg+2]); + special_lj[3] = special_coul[3] = atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"lj") == 0) { + if (iarg+4 > narg) error->all("Illegal special_bonds command"); + special_lj[1] = atof(arg[iarg+1]); + special_lj[2] = atof(arg[iarg+2]); + special_lj[3] = atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"coul") == 0) { + if (iarg+4 > narg) error->all("Illegal special_bonds command"); + special_coul[1] = atof(arg[iarg+1]); + special_coul[2] = atof(arg[iarg+2]); + special_coul[3] = atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"dihedral") == 0) { + if (iarg+2 > narg) error->all("Illegal special_bonds command"); + if (strcmp(arg[iarg+1],"no") == 0) special_dihedral = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) special_dihedral = 1; + else error->all("Illegal special_bonds command"); + iarg += 2; + } else if (strcmp(arg[iarg],"extra") == 0) { + if (iarg+2 > narg) error->all("Illegal special_bonds command"); + special_extra = atoi(arg[iarg+1]); + iarg += 2; + } else error->all("Illegal special_bonds command"); } - if (strcmp(arg[0],"amber") == 0) { - if (narg != 1) error->all("Illegal special_bonds command"); - special_lj[1] = 0.0; - special_lj[2] = 0.0; - special_lj[3] = 0.5; - special_coul[1] = 0.0; - special_coul[2] = 0.0; - special_coul[3] = 5.0/6.0; - special_dihedral = 0; - return; - } + for (int i = 1; i <= 3; i++) + if (special_lj[i] < 0.0 || special_lj[i] > 1.0 || + special_coul[i] < 0.0 || special_coul[i] > 1.0) + error->all("Illegal special_bonds command"); - int iarg; - if (strcmp(arg[0],"dihedral") == 0) { - special_dihedral = 1; - iarg = 1; - } else if (strcmp(arg[0],"explicit") == 0) { - special_dihedral = 0; - iarg = 1; - } else { - special_dihedral = 0; - iarg = 0; - } - - if (narg-iarg == 3) { - special_lj[1] = special_coul[1] = atof(arg[iarg+0]); - special_lj[2] = special_coul[2] = atof(arg[iarg+1]); - special_lj[3] = special_coul[3] = atof(arg[iarg+2]); - } else if (narg-iarg == 6) { - special_lj[1] = atof(arg[iarg+0]); - special_lj[2] = atof(arg[iarg+1]); - special_lj[3] = atof(arg[iarg+2]); - special_coul[1] = atof(arg[iarg+3]); - special_coul[2] = atof(arg[iarg+4]); - special_coul[3] = atof(arg[iarg+5]); - } else error->all("Illegal special_bonds command"); + if (special_extra < 0) error->all("Illegal special_bonds command"); } /* ---------------------------------------------------------------------- diff --git a/src/force.h b/src/force.h index 5cc62080d0..96a8ab8b4b 100644 --- a/src/force.h +++ b/src/force.h @@ -54,6 +54,7 @@ class Force : protected Pointers { double special_coul[4]; // 1-2, 1-3, 1-4 prefactors for Coulombics int special_dihedral; // 0 if defined dihedrals are ignored // 1 if only weight 1,4 atoms if in a dihedral + int special_extra; // extra space for added bonds Force(class LAMMPS *); ~Force(); diff --git a/src/input.cpp b/src/input.cpp index 4872e86d68..ba582b14fe 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -1080,13 +1080,15 @@ void Input::shape() void Input::special_bonds() { - // store 1-3,1-4 and dihedral flag values before change + // store 1-3,1-4 and dihedral/extra flag values before change + // change in 1-2 coeffs will not change the special list double lj2 = force->special_lj[2]; double lj3 = force->special_lj[3]; double coul2 = force->special_coul[2]; double coul3 = force->special_coul[3]; int dihedral = force->special_dihedral; + int extra = force->special_extra; force->set_special(narg,arg); @@ -1095,7 +1097,7 @@ void Input::special_bonds() if (domain->box_exist && atom->molecular) { if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] || coul2 != force->special_coul[2] || coul3 != force->special_coul[3] || - dihedral != force->special_dihedral) { + dihedral != force->special_dihedral || extra != force->special_extra) { Special special(lmp); special.build(); } diff --git a/src/integrate.cpp b/src/integrate.cpp index 1225a56a8b..372a5f8227 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -13,6 +13,7 @@ #include "stdlib.h" #include "integrate.h" +#include "update.h" #include "modify.h" #include "compute.h" @@ -104,6 +105,8 @@ void Integrate::ev_set(int ntimestep) if (elist_atom[i]->matchstep(ntimestep)) break; if (i < nelist_atom) eflag_atom = 2; + if (eflag_global) update->eflag_global = update->ntimestep; + if (eflag_atom) update->eflag_atom = update->ntimestep; eflag = eflag_global + eflag_atom; int vflag_global = 0; @@ -116,5 +119,7 @@ void Integrate::ev_set(int ntimestep) if (vlist_atom[i]->matchstep(ntimestep)) break; if (i < nvlist_atom) vflag_atom = 4; + if (vflag_global) update->vflag_global = update->ntimestep; + if (vflag_atom) update->vflag_atom = update->ntimestep; vflag = vflag_global + vflag_atom; } diff --git a/src/min.cpp b/src/min.cpp index bf141fd8f7..d134d75e5f 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -14,6 +14,7 @@ #include "stdlib.h" #include "string.h" #include "min.h" +#include "update.h" #include "modify.h" #include "compute.h" #include "error.h" @@ -116,6 +117,8 @@ void Min::ev_set(int ntimestep) if (elist_atom[i]->matchstep(ntimestep)) break; if (i < nelist_atom) eflag_atom = 2; + if (eflag_global) update->eflag_global = update->ntimestep; + if (eflag_atom) update->eflag_atom = update->ntimestep; eflag = eflag_global + eflag_atom; int vflag_global = 0; @@ -128,5 +131,7 @@ void Min::ev_set(int ntimestep) if (vlist_atom[i]->matchstep(ntimestep)) break; if (i < nvlist_atom) vflag_atom = 4; + if (vflag_global) update->vflag_global = update->ntimestep; + if (vflag_atom) update->vflag_atom = update->ntimestep; vflag = vflag_global + vflag_atom; } diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 887b725e90..a1dde826e9 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -170,8 +170,8 @@ void MinCG::run() // set output->next values to this timestep // call eng_force to insure vflag is set when forces computed // output->write does final output for thermo, dump, restart files - // add ntimestep to ALL computes that store invocation times - // since just hardwired call to thermo/dumps and they may not be ready + // add ntimestep to all computes that store invocation times + // since just hardwired call to thermo/dumps and computes may not be ready if (niter < update->nsteps) { niter++; diff --git a/src/modify.cpp b/src/modify.cpp index bece86e77c..f6ff32c45b 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -36,37 +36,45 @@ using namespace LAMMPS_NS; // mask settings - same as in fix.cpp -#define INITIAL_INTEGRATE 1 -#define PRE_EXCHANGE 2 -#define PRE_NEIGHBOR 4 -#define POST_FORCE 8 -#define FINAL_INTEGRATE 16 -#define END_OF_STEP 32 -#define THERMO_ENERGY 64 -#define INITIAL_INTEGRATE_RESPA 128 -#define POST_FORCE_RESPA 256 -#define FINAL_INTEGRATE_RESPA 512 -#define MIN_POST_FORCE 1024 -#define MIN_ENERGY 2048 +#define INITIAL_INTEGRATE 1 +#define POST_INTEGRATE 2 +#define PRE_EXCHANGE 4 +#define PRE_NEIGHBOR 8 +#define PRE_FORCE 16 +#define POST_FORCE 32 +#define FINAL_INTEGRATE 64 +#define END_OF_STEP 128 +#define THERMO_ENERGY 256 +#define INITIAL_INTEGRATE_RESPA 512 +#define POST_INTEGRATE_RESPA 1024 +#define PRE_FORCE_RESPA 2048 +#define POST_FORCE_RESPA 4096 +#define FINAL_INTEGRATE_RESPA 8192 +#define MIN_POST_FORCE 16384 +#define MIN_ENERGY 32768 /* ---------------------------------------------------------------------- */ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) { nfix = maxfix = 0; - n_initial_integrate = 0; + n_initial_integrate = n_post_integrate = 0; n_pre_exchange = n_pre_neighbor = 0; - n_post_force = n_final_integrate = n_end_of_step = n_thermo_energy = 0; - n_initial_integrate_respa = n_post_force_respa = n_final_integrate_respa = 0; + n_pre_force = n_post_force = 0; + n_final_integrate = n_end_of_step = n_thermo_energy = 0; + n_initial_integrate_respa = n_post_integrate_respa = 0; + n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; n_min_post_force = n_min_energy = 0; fix = NULL; fmask = NULL; - list_initial_integrate = NULL; + list_initial_integrate = list_post_integrate = NULL; list_pre_exchange = list_pre_neighbor = NULL; - list_post_force = list_final_integrate = list_end_of_step = NULL; + list_pre_force = list_post_force = NULL; + list_final_integrate = list_end_of_step = NULL; list_thermo_energy = NULL; - list_initial_integrate_respa = list_post_force_respa = NULL; + list_initial_integrate_respa = list_post_integrate_respa = NULL; + list_pre_force_respa = list_post_force_respa = NULL; list_final_integrate_respa = NULL; list_min_post_force = list_min_energy = NULL; @@ -101,13 +109,17 @@ Modify::~Modify() memory->sfree(compute); delete [] list_initial_integrate; + delete [] list_post_integrate; delete [] list_pre_exchange; delete [] list_pre_neighbor; + delete [] list_pre_force; delete [] list_post_force; delete [] list_final_integrate; delete [] list_end_of_step; delete [] list_thermo_energy; delete [] list_initial_integrate_respa; + delete [] list_post_integrate_respa; + delete [] list_pre_force_respa; delete [] list_post_force_respa; delete [] list_final_integrate_respa; delete [] list_min_post_force; @@ -134,8 +146,10 @@ void Modify::init() // create lists of fixes to call at each stage of run list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate); + list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate); list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange); list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor); + list_init(PRE_FORCE,n_pre_force,list_pre_force); list_init(POST_FORCE,n_post_force,list_post_force); list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate); list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step); @@ -143,8 +157,12 @@ void Modify::init() list_init(INITIAL_INTEGRATE_RESPA, n_initial_integrate_respa,list_initial_integrate_respa); + list_init(POST_INTEGRATE_RESPA, + n_post_integrate_respa,list_post_integrate_respa); list_init(POST_FORCE_RESPA, n_post_force_respa,list_post_force_respa); + list_init(PRE_FORCE_RESPA, + n_pre_force_respa,list_pre_force_respa); list_init(FINAL_INTEGRATE_RESPA, n_final_integrate_respa,list_final_integrate_respa); @@ -161,9 +179,21 @@ void Modify::init() for (i = 0; i < nfix; i++) fix[i]->init(); // init each compute - // add initial step to ALL computes that store invocation times + // reset their invoked flags + // else may not be invoked on initial timestep of this run + // if were invoked on last step of previous run + // would be bad if state of system changed between runs + // add initial timestep to all computes that store invocation times + // since any of them may be invoked by initial thermo + // do not clear stored invocation times within compute, + // b/c some may be holdovers from previous run, like for ave fixes - for (i = 0; i < ncompute; i++) compute[i]->init(); + for (i = 0; i < ncompute; i++) { + compute[i]->init(); + compute[i]->invoked_scalar = -1; + compute[i]->invoked_vector = -1; + compute[i]->invoked_peratom = -1; + } modify->addstep_compute_all(update->ntimestep); // set global flag if any fix has its restart_pbc flag set @@ -222,6 +252,16 @@ void Modify::initial_integrate(int vflag) fix[list_initial_integrate[i]]->initial_integrate(vflag); } +/* ---------------------------------------------------------------------- + post_integrate call only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::post_integrate() +{ + for (int i = 0; i < n_post_integrate; i++) + fix[list_post_integrate[i]]->post_integrate(); +} + /* ---------------------------------------------------------------------- pre_exchange call only for relevant fixes ------------------------------------------------------------------------- */ @@ -243,7 +283,17 @@ void Modify::pre_neighbor() } /* ---------------------------------------------------------------------- - force adjustment call only for relevant fixes + pre_force call only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::pre_force(int vflag) +{ + for (int i = 0; i < n_pre_force; i++) + fix[list_pre_force[i]]->pre_force(vflag); +} + +/* ---------------------------------------------------------------------- + post_force call only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_force(int vflag) @@ -300,7 +350,27 @@ void Modify::initial_integrate_respa(int vflag, int ilevel, int flag) } /* ---------------------------------------------------------------------- - rRESPA force adjustment call only for relevant fixes + rRESPA post_integrate call only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::post_integrate_respa(int ilevel, int flag) +{ + for (int i = 0; i < n_post_integrate_respa; i++) + fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,flag); +} + +/* ---------------------------------------------------------------------- + rRESPA pre_force call only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::pre_force_respa(int vflag, int ilevel, int iloop) +{ + for (int i = 0; i < n_pre_force_respa; i++) + fix[list_pre_force_respa[i]]->pre_force_respa(vflag,ilevel,iloop); +} + +/* ---------------------------------------------------------------------- + rRESPA post_force call only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_force_respa(int vflag, int ilevel, int iloop) @@ -614,39 +684,40 @@ int Modify::find_compute(char *id) } /* ---------------------------------------------------------------------- - clear invoked flag of all computes - always called before computes are invoked + loop only over computes that store invocation times + clear their invoked flag for this step + called before computes are invoked ------------------------------------------------------------------------- */ void Modify::clearstep_compute() { - for (int icompute = 0; icompute < ncompute; icompute++) - compute[icompute]->invoked = 0; + for (int icompute = 0; icompute < n_timeflag; icompute++) + compute[list_timeflag[icompute]]->invoked_flag = 0; } /* ---------------------------------------------------------------------- loop only over computes that store invocation times - if it was invoked, schedule the next invocation - always called after computes are invoked + if its invoked flag set on this timestep, schedule the next invocation + called after computes are invoked ------------------------------------------------------------------------- */ -void Modify::addstep_compute(int ntimestep) +void Modify::addstep_compute(int newstep) { for (int icompute = 0; icompute < n_timeflag; icompute++) - if (compute[list_timeflag[icompute]]->invoked) - compute[list_timeflag[icompute]]->addstep(ntimestep); + if (compute[list_timeflag[icompute]]->invoked_flag) + compute[list_timeflag[icompute]]->addstep(newstep); } /* ---------------------------------------------------------------------- - loop only all computes + loop over all computes schedule the next invocation for those that store invocation times + called when not sure what computes will be needed on newstep ------------------------------------------------------------------------- */ -void Modify::addstep_compute_all(int ntimestep) +void Modify::addstep_compute_all(int newstep) { - for (int icompute = 0; icompute < ncompute; icompute++) - if (compute[icompute]->timeflag) - compute[icompute]->addstep(ntimestep); + for (int icompute = 0; icompute < n_timeflag; icompute++) + compute[list_timeflag[icompute]]->addstep(newstep); } /* ---------------------------------------------------------------------- diff --git a/src/modify.h b/src/modify.h index 5ed6a12d6f..e969192c0b 100644 --- a/src/modify.h +++ b/src/modify.h @@ -22,9 +22,11 @@ namespace LAMMPS_NS { class Modify : protected Pointers { public: int nfix,maxfix; - int n_initial_integrate,n_pre_decide,n_pre_exchange,n_pre_neighbor; - int n_post_force,n_final_integrate,n_end_of_step,n_thermo_energy; - int n_initial_integrate_respa,n_post_force_respa,n_final_integrate_respa; + int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor; + int n_pre_force,n_post_force; + int n_final_integrate,n_end_of_step,n_thermo_energy; + int n_initial_integrate_respa,n_post_integrate_respa; + int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa; int n_min_post_force,n_min_energy; int nfix_restart_peratom; @@ -41,15 +43,19 @@ class Modify : protected Pointers { void init(); void setup(int); void initial_integrate(int); + void post_integrate(); void pre_decide(); void pre_exchange(); void pre_neighbor(); + void pre_force(int); void post_force(int); void final_integrate(); void end_of_step(); double thermo_energy(); void initial_integrate_respa(int,int,int); + void post_integrate_respa(int,int); + void pre_force_respa(int,int,int); void post_force_respa(int,int,int); void final_integrate_respa(int); @@ -78,12 +84,13 @@ class Modify : protected Pointers { double memory_usage(); private: - // lists of fixes to apply at different times - int *list_initial_integrate,*list_pre_decide; + // lists of fixes to apply at different stages of timestep + int *list_initial_integrate,*list_post_integrate; int *list_pre_exchange,*list_pre_neighbor; - int *list_post_force,*list_final_integrate,*list_end_of_step; - int *list_thermo_energy; - int *list_initial_integrate_respa,*list_post_force_respa; + int *list_pre_force,*list_post_force; + int *list_final_integrate,*list_end_of_step,*list_thermo_energy; + int *list_initial_integrate_respa,*list_post_integrate_respa; + int *list_pre_force_respa,*list_post_force_respa; int *list_final_integrate_respa; int *list_min_post_force,*list_min_energy; diff --git a/src/pair.h b/src/pair.h index 3eae3f4041..cffcc62cb6 100644 --- a/src/pair.h +++ b/src/pair.h @@ -46,7 +46,7 @@ class Pair : protected Pointers { 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 - class NeighList *listgranhistory; // granular history list used by some + class NeighList *listgranhistory; // granular history list used by some pairs class NeighList *listinner; // rRESPA lists used by some pairs class NeighList *listmiddle; class NeighList *listouter; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index c3115e0b7f..33aa5cf5db 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -193,11 +193,14 @@ void PairHybrid::settings(int narg, char **arg) // count sub-styles by skipping numeric args // one exception is 1st arg of style "table", which is non-numeric word + // one exception is 1st two args of style "lj/coul", which are non-numeric + // need a better way to skip these exceptions nstyles = 0; i = 0; while (i < narg) { if (strcmp(arg[i],"table") == 0) i++; + if (strcmp(arg[i],"lj/coul") == 0) i += 2; i++; while (i < narg && !isalpha(arg[i][0])) i++; nstyles++; @@ -304,9 +307,9 @@ void PairHybrid::coeff(int narg, char **arg) } // set setflag and which type pairs map to which sub-style - // if sub-style is none: set hybrid subflag, wipe out map + // if sub-style is none: set hybrid setflag, wipe out map // else: set hybrid setflag & map only if substyle setflag is set - // previous mappings are wiped out + // previous mappings are wiped out int count = 0; for (int i = ilo; i <= ihi; i++) { diff --git a/src/pair_hybrid_overlay.cpp b/src/pair_hybrid_overlay.cpp index 42f1f8e4b1..44675c94bb 100644 --- a/src/pair_hybrid_overlay.cpp +++ b/src/pair_hybrid_overlay.cpp @@ -67,8 +67,8 @@ void PairHybridOverlay::coeff(int narg, char **arg) // set setflag and which type pairs map to which sub-style // if sub-style is none: set hybrid subflag, wipe out map // else: set hybrid setflag & map only if substyle setflag is set - // multiple mappings are allowed, - // but don't map again if already mapped to this sub-style + // multiple mappings are allowed, + // but don't map again if already mapped to this sub-style int count = 0; for (int i = ilo; i <= ihi; i++) { diff --git a/src/read_data.cpp b/src/read_data.cpp index bfeb3abb80..a2ae8cef8c 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -38,6 +38,8 @@ using namespace LAMMPS_NS; #define CHUNK 1024 #define DELTA 4 +#define NSECTIONS 22 // change when add to header::section_keywords + #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) @@ -82,14 +84,17 @@ void ReadData::command(int narg, char **arg) if (screen) fprintf(screen,"Scanning data file ...\n"); open(arg[0]); header(0); - scan(&atom->bond_per_atom,&atom->angle_per_atom, - &atom->dihedral_per_atom,&atom->improper_per_atom); + scan(atom->bond_per_atom,atom->angle_per_atom, + atom->dihedral_per_atom,atom->improper_per_atom); fclose(fp); + atom->bond_per_atom += atom->extra_bond_per_atom; } + MPI_Bcast(&atom->bond_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->angle_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->dihedral_per_atom,1,MPI_INT,0,world); MPI_Bcast(&atom->improper_per_atom,1,MPI_INT,0,world); + } else atom->bond_per_atom = atom->angle_per_atom = atom->dihedral_per_atom = atom->improper_per_atom = 0; @@ -122,9 +127,9 @@ void ReadData::command(int narg, char **arg) domain->set_local_box(); // read rest of file in free format + // if add a section keyword, add to header::section_keywords and NSECTIONS int atomflag = 0; - parse_keyword(1,1); while (strlen(keyword)) { if (strcmp(keyword,"Atoms") == 0) { @@ -242,7 +247,7 @@ void ReadData::command(int narg, char **arg) } else { char str[128]; - sprintf(str,"Unknown section in data file: %s",keyword); + sprintf(str,"Unknown identifier in data file: %s",keyword); error->all(str); } @@ -281,6 +286,15 @@ void ReadData::header(int flag) int n; char *ptr; + char *section_keywords[NSECTIONS] = + {"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers", + "Masses","Shapes","Dipoles", + "Pair Coeffs","Bond Coeffs","Angle Coeffs", + "Dihedral Coeffs","Improper Coeffs", + "BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs", + "EndBondTorsion Coeffs","AngleTorsion Coeffs", + "AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"}; + // skip 1st line of file if (me == 0) { @@ -322,6 +336,7 @@ void ReadData::header(int flag) else if (strstr(line,"angles")) sscanf(line,"%d",&atom->nangles); else if (strstr(line,"dihedrals")) sscanf(line,"%d",&atom->ndihedrals); else if (strstr(line,"impropers")) sscanf(line,"%d",&atom->nimpropers); + else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes); else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes); else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes); @@ -329,6 +344,10 @@ void ReadData::header(int flag) sscanf(line,"%d",&atom->ndihedraltypes); else if (strstr(line,"improper types")) sscanf(line,"%d",&atom->nimpropertypes); + + else if (strstr(line,"extra bond per atom")) + sscanf(line,"%d",&atom->extra_bond_per_atom); + else if (strstr(line,"xlo xhi")) sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]); else if (strstr(line,"ylo yhi")) @@ -341,6 +360,17 @@ void ReadData::header(int flag) } else break; } + // check that exiting string is a valid section keyword + + parse_keyword(1,flag); + for (n = 0; n < NSECTIONS; n++) + if (strcmp(keyword,section_keywords[n]) == 0) break; + if (n == NSECTIONS) { + char str[128]; + sprintf(str,"Unknown identifier in data file: %s",keyword); + error->all(str); + } + // error check on consistency of header values if ((atom->nbonds || atom->nbondtypes) && @@ -914,15 +944,14 @@ void ReadData::impropercoeffs(int which) proc 0 scans the data file for topology maximums ------------------------------------------------------------------------- */ -void ReadData::scan(int *pbond_per, int *pangle_per, - int *pdihedral_per, int *pimproper_per) +void ReadData::scan(int &bond_per_atom, int &angle_per_atom, + int &dihedral_per_atom, int &improper_per_atom) { int i,tmp1,tmp2,atom1,atom2,atom3,atom4; - int bond_per,angle_per,dihedral_per,improper_per; char *eof; int natoms = static_cast (atom->natoms); - bond_per = angle_per = dihedral_per = improper_per = 0; + bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; // allocate topology counting vector // initially, array length = 1 to natoms @@ -931,8 +960,6 @@ void ReadData::scan(int *pbond_per, int *pangle_per, int cmax = natoms + 1; int *count = (int *) memory->smalloc(cmax*sizeof(int),"read_data:count"); - parse_keyword(1,0); - while (strlen(keyword)) { if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes); @@ -1050,9 +1077,9 @@ void ReadData::scan(int *pbond_per, int *pangle_per, count[atom1]++; count[atom2]++; } - for (i = 1; i < cmax; i++) bond_per = MAX(bond_per,count[i]); - if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per); - if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per); + for (i = 1; i < cmax; i++) bond_per_atom = MAX(bond_per_atom,count[i]); + if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per_atom); + if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per_atom); } else if (strcmp(keyword,"Angles") == 0) { @@ -1077,9 +1104,9 @@ void ReadData::scan(int *pbond_per, int *pangle_per, count[atom2]++; count[atom3]++; } - for (i = 1; i < cmax; i++) angle_per = MAX(angle_per,count[i]); - if (screen) fprintf(screen," %d = max angles/atom\n",angle_per); - if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per); + for (i = 1; i < cmax; i++) angle_per_atom = MAX(angle_per_atom,count[i]); + if (screen) fprintf(screen," %d = max angles/atom\n",angle_per_atom); + if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per_atom); } else if (strcmp(keyword,"Dihedrals") == 0) { @@ -1108,11 +1135,12 @@ void ReadData::scan(int *pbond_per, int *pangle_per, count[atom3]++; count[atom4]++; } - for (i = 1; i < cmax; i++) dihedral_per = MAX(dihedral_per,count[i]); + for (i = 1; i < cmax; i++) + dihedral_per_atom = MAX(dihedral_per_atom,count[i]); if (screen) - fprintf(screen," %d = max dihedrals/atom\n",dihedral_per); + fprintf(screen," %d = max dihedrals/atom\n",dihedral_per_atom); if (logfile) - fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per); + fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per_atom); } else if (strcmp(keyword,"Impropers") == 0) { for (i = 1; i < cmax; i++) count[i] = 0; @@ -1140,11 +1168,12 @@ void ReadData::scan(int *pbond_per, int *pangle_per, count[atom3]++; count[atom4]++; } - for (i = 1; i < cmax; i++) improper_per = MAX(improper_per,count[i]); + for (i = 1; i < cmax; i++) + improper_per_atom = MAX(improper_per_atom,count[i]); if (screen) - fprintf(screen," %d = max impropers/atom\n",improper_per); + fprintf(screen," %d = max impropers/atom\n",improper_per_atom); if (logfile) - fprintf(logfile," %d = max impropers/atom\n",improper_per); + fprintf(logfile," %d = max impropers/atom\n",improper_per_atom); } else { char str[128]; @@ -1161,17 +1190,11 @@ void ReadData::scan(int *pbond_per, int *pangle_per, // error check that topology was specified in file - if ((atom->nbonds && !bond_per) || (atom->nangles && !angle_per) || - (atom->ndihedrals && !dihedral_per) || - (atom->nimpropers && !improper_per)) + if ((atom->nbonds && !bond_per_atom) || + (atom->nangles && !angle_per_atom) || + (atom->ndihedrals && !dihedral_per_atom) || + (atom->nimpropers && !improper_per_atom)) error->one("Needed topology not in data file"); - - // return values - - *pbond_per = bond_per; - *pangle_per = angle_per; - *pdihedral_per = dihedral_per; - *pimproper_per = improper_per; } /* ---------------------------------------------------------------------- diff --git a/src/read_data.h b/src/read_data.h index c18b23002f..bff0c4d58d 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -33,7 +33,7 @@ class ReadData : protected Pointers { char **arg; void open(char *); - void scan(int *, int *, int *, int*); + void scan(int &, int &, int &, int &); int reallocate(int **, int, int); void header(int); void parse_keyword(int, int); diff --git a/src/respa.cpp b/src/respa.cpp index d8df6e9243..1e1350fdc6 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -415,6 +415,8 @@ void Respa::recurse(int ilevel) for (int iloop = 0; iloop < loop[ilevel]; iloop++) { modify->initial_integrate_respa(vflag,ilevel,0); + if (modify->n_post_integrate_respa) + modify->post_integrate_respa(ilevel,iloop); if (ilevel) recurse(ilevel-1); @@ -427,6 +429,9 @@ void Respa::recurse(int ilevel) if (ilevel == nlevels-1) { modify->initial_integrate_respa(vflag,ilevel,1); + if (modify->n_post_integrate_respa) + modify->post_integrate_respa(ilevel,iloop); + int nflag = neighbor->decide(); if (nflag) { if (modify->n_pre_exchange) modify->pre_exchange(); @@ -446,6 +451,7 @@ void Respa::recurse(int ilevel) neighbor->build(); timer->stamp(TIME_NEIGHBOR); } + } else if (ilevel == 0) { timer->stamp(); comm->communicate(); @@ -453,6 +459,8 @@ void Respa::recurse(int ilevel) } force_clear(newton[ilevel]); + if (modify->n_pre_force_respa) + modify->pre_force_respa(vflag,ilevel,iloop); timer->stamp(); if (level_bond == ilevel && force->bond) { diff --git a/src/special.cpp b/src/special.cpp index 7b66f60ad0..4ae54d98f6 100644 --- a/src/special.cpp +++ b/src/special.cpp @@ -619,11 +619,13 @@ void Special::combine() } // compute global maxspecial, must be at least 1 + // add in extra factor from special_bonds command // allocate correct special array with same nmax, new maxspecial // previously allocated one must be destroyed - // must force AtomVec class to update its ptr to special + // must make AtomVec class update its ptr to special MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world); + atom->maxspecial += force->special_extra; atom->maxspecial = MAX(atom->maxspecial,1); if (me == 0) { @@ -694,7 +696,7 @@ void Special::combine() /* ---------------------------------------------------------------------- trim list of 1-4 neighbors by checking defined dihedrals - delete a 1-4 neigh if those atoms are not end points of a defined dihedral + delete a 1-4 neigh if they are not end atoms of a defined dihedral ------------------------------------------------------------------------- */ void Special::dihedral_trim() @@ -827,15 +829,6 @@ void Special::dihedral_trim() } else for (i = 0; i < nlocal; i++) nspecial[i][2] = nspecial[i][1]; - // reset atom->maxspecial - - int maxspecial = 0; - for (i = 0; i < nlocal; i++) - maxspecial = MAX(maxspecial,nspecial[i][2]); - - MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world); - atom->maxspecial = MAX(atom->maxspecial,1); - // stats on new 1-4 neighbor counts onefourcount = 0.0; diff --git a/src/style_molecule.h b/src/style_molecule.h index 4ee90a09f3..68e98d6764 100644 --- a/src/style_molecule.h +++ b/src/style_molecule.h @@ -90,10 +90,14 @@ DumpStyle(bond,DumpBond) #endif #ifdef FixInclude +#include "fix_bond_break.h" +#include "fix_bond_create.h" #include "fix_bond_swap.h" #endif #ifdef FixClass +FixStyle(bond/break,FixBondBreak) +FixStyle(bond/create,FixBondCreate) FixStyle(bond/swap,FixBondSwap) #endif diff --git a/src/thermo.cpp b/src/thermo.cpp index d22b2adb84..59ad878a1f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -60,9 +60,6 @@ enum{INT,FLOAT}; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) -#define INVOKED_SCALAR 1 // same as in computes -#define INVOKED_VECTOR 2 - /* ---------------------------------------------------------------------- */ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) @@ -73,8 +70,6 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) style = new char[n]; strcpy(style,arg[0]); - thermoflag = 1; - // set thermo_modify defaults normuserflag = 0; @@ -261,6 +256,7 @@ void Thermo::compute(int flag) int i; firststep = flag; + int ntimestep = update->ntimestep; // check for lost atoms // turn off normflag if natoms = 0 to avoid divide by 0 @@ -274,12 +270,13 @@ void Thermo::compute(int flag) for (i = 0; i < ncompute; i++) { if (compute_which[i] == 0) { - if (!(computes[i]->invoked & INVOKED_SCALAR)) - double tmp = computes[i]->compute_scalar(); + if (computes[i]->invoked_scalar != ntimestep) + computes[i]->compute_scalar(); } else { - if (!(computes[i]->invoked & INVOKED_VECTOR)) + if (computes[i]->invoked_vector != ntimestep) computes[i]->compute_vector(); } + computes[i]->invoked_flag = 1; } // if lineflag = MULTILINE, prepend step/cpu header line @@ -289,7 +286,7 @@ void Thermo::compute(int flag) double cpu; if (flag) cpu = timer->elapsed(TIME_LOOP); else cpu = 0.0; - loc = sprintf(&line[loc],format_multi,update->ntimestep,cpu); + loc = sprintf(&line[loc],format_multi,ntimestep,cpu); } // add each thermo value to line with its specific format @@ -301,7 +298,7 @@ void Thermo::compute(int flag) } // kludge for RedStorm timing issue - // if (update->ntimestep == 100) return; + // if (ntimestep == 100) return; // print line to screen and logfile @@ -819,89 +816,203 @@ void Thermo::create_compute(char *id, char *cstyle, char *extra) } /* ---------------------------------------------------------------------- - compute a single thermodyanmic value - word is any supported keyword in custom list + compute a single thermodyanmic value, word is any keyword in custom list called when a variable is evaluated by Variable class return value as double in answer - return 0 if OK, 1 if str is invalid - customize a new keyword by adding to if statement and error tests + return 0 if str is recoginzed keyword, 1 if unrecognized + customize a new keyword by adding to if statement ------------------------------------------------------------------------- */ int Thermo::evaluate_keyword(char *word, double *answer) { - // check if Compute pointers exist for keywords that need them - // error if thermo style does not use the Compute - - if (strcmp(word,"temp") == 0 || strcmp(word,"press") == 0 || - strcmp(word,"ke") == 0 || strcmp(word,"etotal") == 0 || - strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 || - strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 || - strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0) - if (!temperature) - error->all("Variable uses compute via thermo keyword that thermo does not"); - - if (strcmp(word,"press") == 0 || - strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 || - strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 || - strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0) - if (!pressure) - error->all("Variable uses compute via thermo keyword that thermo does not"); - - if (strcmp(word,"pe") == 0 || strcmp(word,"etotal") == 0 || - strcmp(word,"enthalpy") == 0) - if (!pe) - error->all("Variable uses compute via thermo keyword that thermo does not"); - - // set compute_pe invocation flag for keywords that use energy - // but don't call compute_pe explicitly - - if (strcmp(word,"evdwl") == 0 || strcmp(word,"ecoul") == 0 || - strcmp(word,"epair") == 0 || strcmp(word,"ebond") == 0 || - strcmp(word,"eangle") == 0 || strcmp(word,"edihed") == 0 || - strcmp(word,"eimp") == 0 || strcmp(word,"emol") == 0 || - strcmp(word,"elong") == 0 || strcmp(word,"etail") == 0) { - if (!pe) - error->all("Variable uses compute via thermo keyword that thermo does not"); - pe->invoked |= INVOKED_SCALAR; - } - - // toggle thermoflag off/on - // so individual compute routines know they are not being called from - // Thermo::compute() which pre-calls Computes - - thermoflag = 0; - - // invoke the lo-level thermo routine to compute the variable value + // invoke a lo-level thermo routine to compute the variable value + // if keyword requires a compute, is error if thermo doesn't use the compute + // if in middle of run and needed compute is not current, invoke it + // if inbetween runs and needed compute is not current, error + // set invoked flag for pe and pressure for keywords that use them + // this insures tallying on future needed steps via clearstep/addstep + // for keywords that use pe indirectly (evdwl, ebond, etc): + // check if energy was tallied on this timestep and set pe->invoked_flag if (strcmp(word,"step") == 0) { compute_step(); dvalue = ivalue; + } else if (strcmp(word,"atoms") == 0) { compute_atoms(); dvalue = ivalue; - } - else if (strcmp(word,"cpu") == 0) compute_cpu(); - else if (strcmp(word,"temp") == 0) compute_temp(); + } else if (strcmp(word,"cpu") == 0) { + if (update->whichflag < 0) firststep = 0; + compute_cpu(); - else if (strcmp(word,"press") == 0) compute_press(); - else if (strcmp(word,"pe") == 0) compute_pe(); - else if (strcmp(word,"ke") == 0) compute_ke(); - else if (strcmp(word,"etotal") == 0) compute_etotal(); - else if (strcmp(word,"enthalpy") == 0) compute_enthalpy(); - - else if (strcmp(word,"evdwl") == 0) compute_evdwl(); - else if (strcmp(word,"ecoul") == 0) compute_ecoul(); - else if (strcmp(word,"epair") == 0) compute_epair(); - else if (strcmp(word,"ebond") == 0) compute_ebond(); - else if (strcmp(word,"eangle") == 0) compute_eangle(); - else if (strcmp(word,"edihed") == 0) compute_edihed(); - else if (strcmp(word,"eimp") == 0) compute_eimp(); - else if (strcmp(word,"emol") == 0) compute_emol(); - else if (strcmp(word,"elong") == 0) compute_elong(); - else if (strcmp(word,"etail") == 0) compute_etail(); - - else if (strcmp(word,"vol") == 0) compute_vol(); + } else if (strcmp(word,"temp") == 0) { + if (!temperature) + error->all("Thermo keyword in variable requires thermo to use/init temp"); + if (temperature->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else temperature->compute_scalar(); + } + compute_temp(); + + } else if (strcmp(word,"press") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_scalar(); + } + pressure->invoked_flag = 1; + compute_press(); + + } else if (strcmp(word,"pe") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (pe->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pe->compute_scalar(); + } + pe->invoked_flag = 1; + compute_pe(); + + } else if (strcmp(word,"ke") == 0) { + if (!temperature) + error->all("Thermo keyword in variable requires thermo to use/init temp"); + if (temperature->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else temperature->compute_scalar(); + } + compute_ke(); + + } else if (strcmp(word,"etotal") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (pe->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pe->compute_scalar(); + } + if (!temperature) + error->all("Thermo keyword in variable requires thermo to use/init temp"); + if (temperature->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else temperature->compute_scalar(); + } + pe->invoked_flag = 1; + compute_etotal(); + + } else if (strcmp(word,"enthalpy") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (pe->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pe->compute_scalar(); + } + if (!temperature) + error->all("Thermo keyword in variable requires thermo to use/init temp"); + if (temperature->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else temperature->compute_scalar(); + } + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_scalar(); + } + pe->invoked_flag = 1; + pressure->invoked_flag = 1; + compute_enthalpy(); + + } else if (strcmp(word,"evdwl") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_evdwl(); + + } else if (strcmp(word,"ecoul") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_ecoul(); + + } else if (strcmp(word,"epair") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_epair(); + + } else if (strcmp(word,"ebond") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_ebond(); + + } else if (strcmp(word,"eangle") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_eangle(); + + } else if (strcmp(word,"edihed") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_edihed(); + + } else if (strcmp(word,"eimp") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_eimp(); + + } else if (strcmp(word,"emol") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_emol(); + + } else if (strcmp(word,"elong") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_elong(); + + } else if (strcmp(word,"etail") == 0) { + if (!pe) + error->all("Thermo keyword in variable requires thermo to use/init pe"); + if (update->eflag_global != update->ntimestep) + error->all("Energy was not tallied on needed timestep"); + pe->invoked_flag = 1; + compute_etail(); + + } else if (strcmp(word,"vol") == 0) compute_vol(); else if (strcmp(word,"lx") == 0) compute_lx(); else if (strcmp(word,"ly") == 0) compute_ly(); else if (strcmp(word,"lz") == 0) compute_lz(); @@ -917,16 +1028,73 @@ int Thermo::evaluate_keyword(char *word, double *answer) else if (strcmp(word,"xz") == 0) compute_xz(); else if (strcmp(word,"yz") == 0) compute_yz(); - else if (strcmp(word,"pxx") == 0) compute_pxx(); - else if (strcmp(word,"pyy") == 0) compute_pyy(); - else if (strcmp(word,"pzz") == 0) compute_pzz(); - else if (strcmp(word,"pxy") == 0) compute_pxy(); - else if (strcmp(word,"pxz") == 0) compute_pxz(); - else if (strcmp(word,"pyz") == 0) compute_pyz(); + else if (strcmp(word,"pxx") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pxx(); - else return 1; + } else if (strcmp(word,"pyy") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pyy(); - thermoflag = 1; + } else if (strcmp(word,"pzz") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pzz(); + + } else if (strcmp(word,"pxy") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pxy(); + + } else if (strcmp(word,"pxz") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pxz(); + + } else if (strcmp(word,"pyz") == 0) { + if (!pressure) + error->all("Thermo keyword in variable requires thermo to use/init press"); + if (pressure->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable thermo keyword is not current"); + else pressure->compute_vector(); + } + pressure->invoked_flag = 1; + compute_pyz(); + + } else return 1; *answer = dvalue; return 0; @@ -934,8 +1102,6 @@ int Thermo::evaluate_keyword(char *word, double *answer) /* ---------------------------------------------------------------------- extraction of Compute, Fix, Variable results - ignore thermoflag: - 3 methods only called by Thermo::compute(), not by variable evaluation compute/fix are normalized by atoms if returning extensive value variable value is not normalized (formula should normalize if desired) ------------------------------------------------------------------------- */ @@ -986,8 +1152,8 @@ void Thermo::compute_variable() /* ---------------------------------------------------------------------- one method for every keyword thermo can output - thermoflag = 1 if called by Thermo::compute() - thermoflag = 0 if called by evaluate_keyword() via Variable + called by compute() or evaluate_keyword() + compute will have already been called set ivalue/dvalue if value is integer/double customize a new keyword by adding a method ------------------------------------------------------------------------- */ @@ -1001,8 +1167,7 @@ void Thermo::compute_step() void Thermo::compute_atoms() { - if (thermoflag) ivalue = static_cast (natoms); - else ivalue = static_cast (atom->natoms); + ivalue = static_cast (natoms); } /* ---------------------------------------------------------------------- */ @@ -1017,27 +1182,21 @@ void Thermo::compute_cpu() void Thermo::compute_temp() { - if (thermoflag || temperature->invoked & INVOKED_SCALAR) - dvalue = temperature->scalar; - else dvalue = temperature->compute_scalar(); + dvalue = temperature->scalar; } /* ---------------------------------------------------------------------- */ void Thermo::compute_press() { - if (thermoflag || pressure->invoked & INVOKED_SCALAR) - dvalue = pressure->scalar; - else dvalue = pressure->compute_scalar(); + dvalue = pressure->scalar; } /* ---------------------------------------------------------------------- */ void Thermo::compute_pe() { - if (thermoflag || pe->invoked & INVOKED_SCALAR) - dvalue = pe->scalar; - else dvalue = pe->compute_scalar(); + dvalue = pe->scalar; if (normflag) dvalue /= natoms; } @@ -1045,10 +1204,7 @@ void Thermo::compute_pe() void Thermo::compute_ke() { - if (thermoflag || temperature->invoked & INVOKED_SCALAR) - dvalue = temperature->scalar; - else dvalue = temperature->compute_scalar(); - + dvalue = temperature->scalar; dvalue *= 0.5 * temperature->dof * force->boltz; if (normflag) dvalue /= natoms; } @@ -1058,12 +1214,7 @@ void Thermo::compute_ke() void Thermo::compute_etotal() { compute_pe(); - - double ke; - if (thermoflag || temperature->invoked & INVOKED_SCALAR) - ke = temperature->scalar; - else ke = temperature->compute_scalar(); - + double ke = temperature->scalar; ke *= 0.5 * temperature->dof * force->boltz; if (normflag) ke /= natoms; dvalue += ke; @@ -1307,8 +1458,6 @@ void Thermo::compute_yz() void Thermo::compute_pxx() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[0]; } @@ -1316,8 +1465,6 @@ void Thermo::compute_pxx() void Thermo::compute_pyy() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[1]; } @@ -1325,8 +1472,6 @@ void Thermo::compute_pyy() void Thermo::compute_pzz() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[2]; } @@ -1334,8 +1479,6 @@ void Thermo::compute_pzz() void Thermo::compute_pxy() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[3]; } @@ -1343,8 +1486,6 @@ void Thermo::compute_pxy() void Thermo::compute_pxz() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[4]; } @@ -1352,7 +1493,5 @@ void Thermo::compute_pxz() void Thermo::compute_pyz() { - if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR)) - pressure->compute_vector(); dvalue = pressure->vector[5]; } diff --git a/src/thermo.h b/src/thermo.h index cc14576421..3b7bc60036 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -60,7 +60,6 @@ class Thermo : protected Pointers { int ivalue; // integer value to print double dvalue,natoms; // dvalue = double value to print int ifield; // which field in thermo output is being computed - int thermoflag; // 1 when called by compute(), 0 from variable eval int *field2index; // which compute,fix,variable calcs this field int *argindex; // index into compute,fix scalar,vector diff --git a/src/update.cpp b/src/update.cpp index 5abc8acb61..c28a596ed4 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -16,6 +16,9 @@ #include "update.h" #include "neighbor.h" #include "force.h" +#include "modify.h" +#include "fix.h" +#include "compute.h" #include "output.h" #include "memory.h" #include "error.h" @@ -42,6 +45,8 @@ Update::Update(LAMMPS *lmp) : Pointers(lmp) firststep = laststep = 0; beginstep = endstep = 0; + eflag_global = vflag_global = -1; + unit_style = NULL; set_units("lj"); @@ -207,6 +212,10 @@ void Update::create_minimize(int narg, char **arg) /* ---------------------------------------------------------------------- reset timestep from input script do not allow dump files or a restart to be defined + do not allow any timestep-dependent fixes to be defined + reset eflag/vflag global so nothing will think eng/virial are current + reset invoked flags of computes, so nothing will think they are current + clear timestep list of computes that store one ------------------------------------------------------------------------- */ void Update::reset_timestep(int narg, char **arg) @@ -219,6 +228,21 @@ void Update::reset_timestep(int narg, char **arg) if (output->restart && output->last_restart >= 0) error->all("Cannot reset timestep with restart file already written"); + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->time_depend) + error->all("Cannot reset timestep with a time-dependent fix defined"); + + eflag_global = vflag_global = -1; + + for (int i = 0; i < modify->ncompute; i++) { + modify->compute[i]->invoked_scalar = -1; + modify->compute[i]->invoked_vector = -1; + modify->compute[i]->invoked_peratom = -1; + } + + for (int i = 0; i < modify->ncompute; i++) + if (modify->compute[i]->timeflag) modify->compute[i]->clearstep(); + ntimestep = atoi(arg[0]); if (ntimestep < 0) error->all("Timestep must be >= 0"); } diff --git a/src/update.h b/src/update.h index 1e5847cbe8..a417df4aab 100644 --- a/src/update.h +++ b/src/update.h @@ -30,6 +30,9 @@ class Update : protected Pointers { int first_update; // 0 before initial update, 1 after int max_eval; // max force evaluations for minimizer + int eflag_global,eflag_atom; // timestep global/peratom eng is tallied + int vflag_global,vflag_atom; // ditto for virial + char *unit_style; class Integrate *integrate; diff --git a/src/variable.cpp b/src/variable.cpp index 94f140e4e0..4a12a6f5b9 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -43,10 +43,6 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY, SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN, CEIL,FLOOR,ROUND,VALUE,ATOMARRAY,TYPEARRAY}; -#define INVOKED_SCALAR 1 // same as in computes -#define INVOKED_VECTOR 2 -#define INVOKED_PERATOM 4 - /* ---------------------------------------------------------------------- */ Variable::Variable(LAMMPS *lmp) : Pointers(lmp) @@ -643,15 +639,15 @@ double Variable::evaluate(char *str, Tree **tree) char *id = new char[n]; strcpy(id,&word[2]); - if (update->first_update == 0) - error->all("Compute in variable formula before initial run"); - int icompute = modify->find_compute(id); if (icompute < 0) error->all("Invalid compute ID in variable formula"); Compute *compute = modify->compute[icompute]; - + compute->invoked_flag = 1; delete [] id; + if (domain->box_exist == 0) + error->all("Variable evalulation before simulation box is defined"); + // parse zero or one or two trailing brackets // point i beyond last bracket // nbracket = # of bracket pairs @@ -674,9 +670,12 @@ double Variable::evaluate(char *str, Tree **tree) if (nbracket == 0 && compute->scalar_flag) { - if (!(compute->invoked & INVOKED_SCALAR)) - value1 = compute->compute_scalar(); - else value1 = compute->scalar; + if (compute->invoked_scalar != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_scalar(); + } + value1 = compute->scalar; if (tree) { Tree *newtree = new Tree(); newtree->type = VALUE; @@ -691,7 +690,11 @@ double Variable::evaluate(char *str, Tree **tree) if (index1 > compute->size_vector) error->all("Compute vector in variable formula is too small"); - if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector(); + if (compute->invoked_vector != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_vector(); + } value1 = compute->vector[index1-1]; if (tree) { Tree *newtree = new Tree(); @@ -708,8 +711,11 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all("Per-atom compute in equal-style variable formula"); - if (!(compute->invoked & INVOKED_PERATOM)) - compute->compute_peratom(); + if (compute->invoked_peratom != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_peratom(); + } Tree *newtree = new Tree(); newtree->type = ATOMARRAY; newtree->array = compute->scalar_atom; @@ -722,8 +728,11 @@ double Variable::evaluate(char *str, Tree **tree) } else if (nbracket == 1 && index1 > 0 && compute->peratom_flag && compute->size_peratom == 0) { - if (!(compute->invoked & INVOKED_PERATOM)) - compute->compute_peratom(); + if (compute->invoked_peratom != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_peratom(); + } peratom2global(1,NULL,compute->scalar_atom,1,index1, tree,treestack,ntreestack,argstack,nargstack); @@ -736,8 +745,11 @@ double Variable::evaluate(char *str, Tree **tree) error->all("Per-atom compute in equal-style variable formula"); if (index2 > compute->size_peratom) error->all("Compute vector in variable formula is too small"); - if (!(compute->invoked & INVOKED_PERATOM)) - compute->compute_peratom(); + if (compute->invoked_peratom != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_peratom(); + } Tree *newtree = new Tree(); newtree->type = ATOMARRAY; newtree->array = &compute->vector_atom[0][index2-1]; @@ -752,8 +764,11 @@ double Variable::evaluate(char *str, Tree **tree) if (index2 > compute->size_peratom) error->all("Compute vector in variable formula is too small"); - if (!(compute->invoked & INVOKED_PERATOM)) - compute->compute_peratom(); + if (compute->invoked_peratom != update->ntimestep) { + if (update->whichflag < 0) + error->all("Compute used in variable is not current"); + else compute->compute_peratom(); + } peratom2global(1,NULL,&compute->vector_atom[0][index2-1], compute->size_peratom,index1, tree,treestack,ntreestack,argstack,nargstack); @@ -769,15 +784,14 @@ double Variable::evaluate(char *str, Tree **tree) char *id = new char[n]; strcpy(id,&word[2]); - if (update->first_update == 0) - error->all("Fix in variable formula before initial run"); - int ifix = modify->find_fix(id); if (ifix < 0) error->all("Invalid fix ID in variable formula"); Fix *fix = modify->fix[ifix]; - delete [] id; + if (domain->box_exist == 0) + error->all("Variable evaluation before simulation box is defined"); + // parse zero or one or two trailing brackets // point i beyond last bracket // nbracket = # of bracket pairs @@ -800,7 +814,8 @@ double Variable::evaluate(char *str, Tree **tree) if (nbracket == 0 && fix->scalar_flag) { - if (update->ntimestep % fix->scalar_vector_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->scalar_vector_freq) error->all("Fix in variable not computed at compatible time"); value1 = fix->compute_scalar(); if (tree) { @@ -817,7 +832,8 @@ double Variable::evaluate(char *str, Tree **tree) if (index1 > fix->size_vector) error->all("Fix vector in variable formula is too small"); - if (update->ntimestep % fix->scalar_vector_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->scalar_vector_freq) error->all("Fix in variable not computed at compatible time"); value1 = fix->compute_vector(index1-1); if (tree) { @@ -835,7 +851,8 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all("Per-atom fix in equal-style variable formula"); - if (update->ntimestep % fix->peratom_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->peratom_freq) error->all("Fix in variable not computed at compatible time"); Tree *newtree = new Tree(); newtree->type = ATOMARRAY; @@ -849,7 +866,8 @@ double Variable::evaluate(char *str, Tree **tree) } else if (nbracket == 1 && index1 > 0 && fix->peratom_flag && fix->size_peratom == 0) { - if (update->ntimestep % fix->peratom_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->peratom_freq) error->all("Fix in variable not computed at compatible time"); peratom2global(1,NULL,fix->scalar_atom,1,index1, tree,treestack,ntreestack,argstack,nargstack); @@ -863,7 +881,8 @@ double Variable::evaluate(char *str, Tree **tree) error->all("Per-atom fix in equal-style variable formula"); if (index2 > fix->size_peratom) error->all("Fix vector in variable formula is too small"); - if (update->ntimestep % fix->peratom_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->peratom_freq) error->all("Fix in variable not computed at compatible time"); Tree *newtree = new Tree(); newtree->type = ATOMARRAY; @@ -879,7 +898,8 @@ double Variable::evaluate(char *str, Tree **tree) if (index2 > fix->size_peratom) error->all("Fix vector in variable formula is too small"); - if (update->ntimestep % fix->peratom_freq) + if (update->whichflag >= 0 && + update->ntimestep % fix->peratom_freq) error->all("Fix in variable not computed at compatible time"); peratom2global(1,NULL,&fix->vector_atom[0][index2-1], fix->size_peratom,index1, @@ -984,6 +1004,9 @@ double Variable::evaluate(char *str, Tree **tree) i = int_between_brackets(str,i,id,1); i++; + if (domain->box_exist == 0) + error->all("Variable evaluation before simulation box is defined"); + // ID between brackets exists: atom value // empty brackets: atom vector @@ -998,8 +1021,9 @@ double Variable::evaluate(char *str, Tree **tree) // ---------------- } else { - if (update->first_update == 0) - error->all("Thermo keyword in variable formula before initial run"); + if (domain->box_exist == 0) + error->all("Variable evaluation before simulation box is defined"); + int flag = output->thermo->evaluate_keyword(word,&value1); if (flag) error->all("Invalid thermo keyword in variable formula"); if (tree) { diff --git a/src/verlet.cpp b/src/verlet.cpp index d271a83cdf..c2f995ea0a 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -137,6 +137,7 @@ void Verlet::iterate(int n) // initial time integration modify->initial_integrate(vflag); + if (modify->n_post_integrate) modify->post_integrate(); // regular communication vs neighbor list rebuild @@ -168,6 +169,7 @@ void Verlet::iterate(int n) // force computations force_clear(); + if (modify->n_pre_force) modify->pre_force(vflag); timer->stamp();