git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7184 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/ASPHERE/pair_line_lj.cpp

@@ -96,6 +96,7 @@ void PairLineLJ::compute(int eflag, int vflag)
   // grow discrete list if necessary and initialize
 
   if (nall > nmax) {
+    nmax = nall;
     memory->destroy(dnum);
     memory->destroy(dfirst);
     memory->create(dnum,nall,"pair:dnum");

src/ASPHERE/pair_line_lj.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairLineLJ : public Pair {
  public:
   PairLineLJ(class LAMMPS *);
-  ~PairLineLJ();
-  void compute(int, int);
+  virtual ~PairLineLJ();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);

src/ASPHERE/pair_resquared.cpp

@@ -133,6 +133,8 @@ void PairRESquared::compute(int eflag, int vflag)
       // compute if less than cutoff
 
       if (rsq < cutsq[itype][jtype]) {
+	fforce[0] = fforce[1] = fforce[2] = 0.0;
+
         switch (form[itype][jtype]) {
 
          case SPHERE_SPHERE:

src/ASPHERE/pair_tri_lj.cpp

@@ -100,6 +100,7 @@ void PairTriLJ::compute(int eflag, int vflag)
   // grow discrete list if necessary and initialize
 
   if (nall > nmax) {
+    nmax = nall;
     memory->destroy(dnum);
     memory->destroy(dfirst);
     memory->create(dnum,nall,"pair:dnum");
@@ -118,7 +119,7 @@ void PairTriLJ::compute(int eflag, int vflag)
     itype = type[i];
     jlist = firstneigh[i];
     jnum = numneigh[i];
-    
+
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
       j &= NEIGHMASK;

src/ASPHERE/pair_tri_lj.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairTriLJ : public Pair {
  public:
   PairTriLJ(class LAMMPS *);
-  ~PairTriLJ();
-  void compute(int, int);
+  virtual ~PairTriLJ();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);

src/CLASS2/angle_class2.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class AngleClass2 : public Angle {
  public:
   AngleClass2(class LAMMPS *);
-  ~AngleClass2();
-  void compute(int, int);
+  virtual ~AngleClass2();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *theta0,*k2,*k3,*k4;
   double *bb_k,*bb_r1,*bb_r2;
   double *ba_k1,*ba_k2,*ba_r1,*ba_r2;

src/CLASS2/bond_class2.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondClass2 : public Bond {
  public:
   BondClass2(class LAMMPS *);
-  ~BondClass2();
-  void compute(int, int);
+  virtual ~BondClass2();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *r0,*k2,*k3,*k4;
 
   void allocate();

src/CLASS2/improper_class2.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class ImproperClass2 : public Improper {
  public:
   ImproperClass2(class LAMMPS *);
-  ~ImproperClass2();
-  void compute(int, int);
+  virtual ~ImproperClass2();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k0,*chi0;
   double *aa_k1,*aa_k2,*aa_k3,*aa_theta0_1,*aa_theta0_2,*aa_theta0_3;
   int *setflag_i,*setflag_aa;

src/GPU/fix_gpu.cpp

@@ -43,7 +43,8 @@ extern double lmp_gpu_forces(double **f, double **tor, double *eatom,
 FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (lmp->cuda) error->all(FLERR,"Cannot use fix GPU with USER-CUDA mode enabled");
+  if (lmp->cuda) 
+    error->all(FLERR,"Cannot use fix GPU with USER-CUDA mode enabled");
 
   if (narg < 7) error->all(FLERR,"Illegal fix GPU command");
   if (strcmp(arg[1],"all") != 0) error->all(FLERR,"Illegal fix GPU command");
@@ -113,10 +114,8 @@ int FixGPU::setmask()
 
 void FixGPU::init()
 {
-  // Can only have 1 gpu fix that must be the first fix for a run
-  if ((void*)modify->fix[0] != (void*)this)
-    error->all(FLERR,"GPU is not the first fix for this run");
-  // Hybrid cannot be used with force/neigh option
+  // hybrid cannot be used with force/neigh option
+
   if (_gpu_mode == GPU_NEIGH)
     if (force->pair_match("hybrid",1) != NULL ||
 	force->pair_match("hybrid/overlay",1) != NULL)

src/GRANULAR/fix_pour.cpp

@@ -168,8 +168,10 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
   // assume grav = -magnitude at this point, enforce in init()
 
   int ifix;
-  for (ifix = 0; ifix < modify->nfix; ifix++)
+  for (ifix = 0; ifix < modify->nfix; ifix++) {
     if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
+    if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
+  }
   if (ifix == modify->nfix) 
     error->all(FLERR,"No fix gravity defined for fix pour");
   grav = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v;
@@ -267,8 +269,10 @@ void FixPour::init()
   // else insertion cannot work
 
   int ifix;
-  for (ifix = 0; ifix < modify->nfix; ifix++)
+  for (ifix = 0; ifix < modify->nfix; ifix++) {
     if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
+    if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
+  }
   if (ifix == modify->nfix) 
     error->all(FLERR,"No fix gravity defined for fix pour");
 

src/GRANULAR/fix_wall_gran.h

@@ -27,12 +27,12 @@ namespace LAMMPS_NS {
 class FixWallGran : public Fix {
  public:
   FixWallGran(class LAMMPS *, int, char **);
-  ~FixWallGran();
+  virtual ~FixWallGran();
   int setmask();
   void init();
   void setup(int);
-  void post_force(int);
-  void post_force_respa(int, int, int);
+  virtual void post_force(int);
+  virtual void post_force_respa(int, int, int);
 
   double memory_usage();
   void grow_arrays(int);
@@ -46,7 +46,7 @@ class FixWallGran : public Fix {
   int maxsize_restart();
   void reset_dt();
 
- private:
+ protected:
   int wallstyle,pairstyle,wiggle,wshear,axis;
   double kn,kt,gamman,gammat,xmu;
   double lo,hi,cylradius;

src/KSPACE/ewald.cpp

@@ -282,7 +282,7 @@ void Ewald::compute(int eflag, int vflag)
     for (k = 0; k < kcount; k++)
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + 
 			 sfacim_all[k]*sfacim_all[k]);
-    energy -= g_ewald*qsqsum/1.772453851 + 
+    energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e*scale;
   }

src/KSPACE/ewald.h

@@ -27,13 +27,13 @@ namespace LAMMPS_NS {
 class Ewald : public KSpace {
  public:
   Ewald(class LAMMPS *, int, char **);
-  ~Ewald();
+  virtual ~Ewald();
   void init();
   void setup();
-  void compute(int, int);
+  virtual void compute(int, int);
   double memory_usage();
 
- private:
+ protected:
   double precision;
   int kcount,kmax,kmax3d,kmax_created;
   double qqrd2e;
@@ -48,9 +48,9 @@ class Ewald : public KSpace {
   double *sfacrl,*sfacim,*sfacrl_all,*sfacim_all;
   double ***cs,***sn;
 
-  void eik_dot_r();
+  virtual void eik_dot_r();
   void coeffs();
-  void allocate();
+  virtual void allocate();
   void deallocate();
   void slabcorr(int);
 };

src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp

@@ -400,7 +400,8 @@ void PairLJCutCoulLongTIP4P::init_style()
     error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires newton pair on");
   if (!atom->q_flag)
     error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires atom attribute q");
-  if (strcmp(force->kspace_style,"pppm/tip4p") != 0)
+  if ( (strcmp(force->kspace_style,"pppm/tip4p") != 0) &&
+       (strcmp(force->kspace_style,"pppm/tip4p/proxy") != 0) )
     error->all(FLERR,"Pair style is incompatible with KSpace style");
   if (force->bond == NULL)
     error->all(FLERR,"Must use a bond style with TIP4P potential");

src/KSPACE/pppm.cpp

@@ -155,7 +155,8 @@ void PPPM::init()
 
   qdist = 0.0;
 
-  if (strcmp(force->kspace_style,"pppm/tip4p") == 0) {
+  if ( (strcmp(force->kspace_style,"pppm/tip4p") == 0) ||
+       (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
     if (force->pair == NULL)
       error->all(FLERR,"KSpace style is incompatible with Pair style");
     double *p_qdist = (double *) force->pair->extract("qdist",itmp);
@@ -184,6 +185,16 @@ void PPPM::init()
     alpha = qdist / (cos(0.5*theta) * blen);
   }
 
+  // if we have a /proxy pppm version check if the pair style is compatible
+
+  if ( (strcmp(force->kspace_style,"pppm/proxy") == 0) ||
+       (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
+    if (force->pair == NULL)
+      error->all(FLERR,"KSpace style is incompatible with Pair style");
+    if (strstr(force->pair_style,"pppm/") == NULL )
+      error->all(FLERR,"KSpace style is incompatible with Pair style");
+  }
+
   // compute qsum & qsqsum and warn if not charge-neutral
 
   qsum = qsqsum = 0.0;
@@ -709,7 +720,7 @@ void PPPM::compute(int eflag, int vflag)
     energy = energy_all;
    
     energy *= 0.5*volume;
-    energy -= g_ewald*qsqsum/1.772453851 +
+    energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e*scale;
   }
@@ -1515,8 +1526,7 @@ void PPPM::make_rho()
 
   // clear 3d density array
 
-  FFT_SCALAR *vec = &density_brick[nzlo_out][nylo_out][nxlo_out];
-  for (i = 0; i < ngrid; i++) vec[i] = ZEROF;
+  memset(&(density_brick[nzlo_out][nylo_out][nxlo_out]),0,ngrid*sizeof(FFT_SCALAR));
 
   // loop over my charges, add their contribution to nearby grid points
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@@ -1777,17 +1787,19 @@ void PPPM::compute_rho1d(const FFT_SCALAR &dx, const FFT_SCALAR &dy,
 			 const FFT_SCALAR &dz)
 {
   int k,l;
+  FFT_SCALAR r1,r2,r3;
 
   for (k = (1-order)/2; k <= order/2; k++) {
-    rho1d[0][k] = ZEROF;
-    rho1d[1][k] = ZEROF;
-    rho1d[2][k] = ZEROF;
+    r1 = r2 = r3 = ZEROF;
 
     for (l = order-1; l >= 0; l--) {
-      rho1d[0][k] = rho_coeff[l][k] + rho1d[0][k]*dx;
-      rho1d[1][k] = rho_coeff[l][k] + rho1d[1][k]*dy;
-      rho1d[2][k] = rho_coeff[l][k] + rho1d[2][k]*dz;
+      r1 = rho_coeff[l][k] + r1*dx;
+      r2 = rho_coeff[l][k] + r2*dy;
+      r3 = rho_coeff[l][k] + r3*dz;
     }
+    rho1d[0][k] = r1;
+    rho1d[1][k] = r2;
+    rho1d[2][k] = r3;
   }
 }
 

src/KSPACE/pppm_cg.cpp

@@ -178,7 +178,7 @@ void PPPMCG::compute(int eflag, int vflag)
     energy = energy_all;
    
     energy *= 0.5*volume;
-    energy -= g_ewald*qsqsum/1.772453851 +
+    energy -= g_ewald*qsqsum/MY_PIS +
       MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e*scale;
   }

src/MANYBODY/pair_adp.cpp

@@ -1022,8 +1022,7 @@ void PairADP::unpack_reverse_comm(int n, int *list, double *buf)
 
 double PairADP::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = Pair::memory_usage();
   bytes += 21 * nmax * sizeof(double);
   return bytes;
 }

src/MOLECULE/angle_charmm.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class AngleCharmm : public Angle {
  public:
   AngleCharmm(class LAMMPS *);
-  ~AngleCharmm();
-  void compute(int, int);
+  virtual ~AngleCharmm();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *k,*theta0,*k_ub,*r_ub;
 
   void allocate();

src/MOLECULE/angle_cosine.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class AngleCosine : public Angle {
  public:
   AngleCosine(class LAMMPS *);
-  ~AngleCosine();
-  void compute(int, int);
+  virtual ~AngleCosine();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *k;
 
   void allocate();

src/MOLECULE/angle_cosine_delta.h

@@ -28,7 +28,7 @@ namespace LAMMPS_NS {
 class AngleCosineDelta : public AngleCosineSquared {
  public:
   AngleCosineDelta(class LAMMPS *);
-  void compute(int, int);
+  virtual void compute(int, int);
   double single(int, int, int, int);
 };
 

src/MOLECULE/angle_cosine_periodic.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class AngleCosinePeriodic : public Angle {
  public:
   AngleCosinePeriodic(class LAMMPS *);
-  ~AngleCosinePeriodic();
-  void compute(int, int);
+  virtual ~AngleCosinePeriodic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *k; 
   int *multiplicity,*b;
 

src/MOLECULE/angle_harmonic.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class AngleHarmonic : public Angle {
  public:
   AngleHarmonic(class LAMMPS *);
-  ~AngleHarmonic();
-  void compute(int, int);
+  virtual ~AngleHarmonic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *k,*theta0;
 
   void allocate();

src/MOLECULE/angle_table.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class AngleTable : public Angle {
  public:
   AngleTable(class LAMMPS *);
-  ~AngleTable();
-  void compute(int, int);
+  virtual ~AngleTable();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double equilibrium_angle(int);
@@ -37,7 +37,7 @@ class AngleTable : public Angle {
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   int tabstyle,tablength;
   double *theta0;
 

src/MOLECULE/bond_fene.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class BondFENE : public Bond {
  public:
   BondFENE(class LAMMPS *);
-  ~BondFENE();
-  void compute(int, int);
+  virtual ~BondFENE();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
   double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondFENE : public Bond {
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double TWO_1_3;
   double *k,*r0,*epsilon,*sigma;
 

src/MOLECULE/bond_fene_expand.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class BondFENEExpand : public Bond {
  public:
   BondFENEExpand(class LAMMPS *);
-  ~BondFENEExpand();
-  void compute(int, int);
+  virtual ~BondFENEExpand();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
   double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondFENEExpand : public Bond {
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double TWO_1_3;
   double *k,*r0,*epsilon,*sigma,*shift;
 

src/MOLECULE/bond_harmonic.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondHarmonic : public Bond {
  public:
   BondHarmonic(class LAMMPS *);
-  ~BondHarmonic();
-  void compute(int, int);
+  virtual ~BondHarmonic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *k,*r0;
 
   void allocate();

src/MOLECULE/bond_morse.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondMorse : public Bond {
  public:
   BondMorse(class LAMMPS *);
-  ~BondMorse();
-  void compute(int, int);
+  virtual ~BondMorse();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *d0,*alpha,*r0;
 
   void allocate();

src/MOLECULE/bond_nonlinear.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondNonlinear : public Bond {
  public:
   BondNonlinear(class LAMMPS *);
-  ~BondNonlinear();
-  void compute(int, int);
+  virtual ~BondNonlinear();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *epsilon,*r0,*lamda;
 
   void allocate();

src/MOLECULE/bond_quartic.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class BondQuartic : public Bond {
  public:
   BondQuartic(class LAMMPS *);
-  ~BondQuartic();
-  void compute(int, int);
+  virtual ~BondQuartic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
   double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondQuartic : public Bond {
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double TWO_1_3;
   double *k,*b1,*b2,*rc,*u0;
 

src/MOLECULE/bond_table.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class BondTable : public Bond {
  public:
   BondTable(class LAMMPS *);
-  ~BondTable();
-  void compute(int, int);
+  virtual ~BondTable();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondTable : public Bond {
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   int tabstyle,tablength;
   double *r0;
 

src/MOLECULE/improper_cvff.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class ImproperCvff : public Improper {
  public:
   ImproperCvff(class LAMMPS *);
-  ~ImproperCvff();
-  void compute(int, int);
+  virtual ~ImproperCvff();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k;
   int *sign,*multiplicity;
 

src/MOLECULE/improper_harmonic.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class ImproperHarmonic : public Improper {
  public:
   ImproperHarmonic(class LAMMPS *);
-  ~ImproperHarmonic();
-  void compute(int, int);
+  virtual ~ImproperHarmonic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k,*chi;
 
   void allocate();

src/MOLECULE/improper_umbrella.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class ImproperUmbrella : public Improper {
  public:
   ImproperUmbrella(class LAMMPS *);
-  ~ImproperUmbrella();
-  void compute(int, int);
+  virtual ~ImproperUmbrella();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *kw, *w0, *C;
 
   void allocate();

src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp

@@ -162,7 +162,7 @@ void PairLJCutCoulLongTIP4POpt::eval()
   double r,r2inv,r6inv,forcecoul,forcelj,cforce;
   double factor_coul,factor_lj;
   double grij,expm2,prefactor,t,erfc,ddotf;
-  double xiM[3],xjM[3],v[6],xH1[3],xH2[3];
+  double v[6],xH1[3],xH2[3];
   double fdx,fdy,fdz,f1x,f1y,f1z,fOx,fOy,fOz,fHx,fHy,fHz;
   double *x1,*x2;
   int *ilist,*jlist,*numneigh,**firstneigh;

src/SHOCK/fix_wall_piston.cpp

@@ -24,8 +24,10 @@
 #include "random_mars.h"
 #include "force.h"
 #include "comm.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
@@ -203,7 +205,7 @@ void FixWallPiston::post_integrate()
   }
   else if (rampNL1flag) {
     paccelz = maxvz / tott;
-    angfreq = 6.283185 / (0.5 * tott);
+    angfreq = MY_2PI / (0.5 * tott);
 
     if (zloflag) {
       zlo = z0 + paccelz * (0.5*tt + 1.0/(angfreq*angfreq) - 1.0/(angfreq*angfreq)*cos(angfreq*t));
@@ -213,7 +215,7 @@ void FixWallPiston::post_integrate()
   }
   else if (rampNL2flag) {
     paccelz = maxvz / tott;
-    angfreq = 18.84956 / tott;
+    angfreq = 3.0*MY_2PI / tott;
 
     if (zloflag) {
       zlo = z0 + paccelz * (0.5*tt + 4.0/(3.0*angfreq*angfreq)*(1.0-cos(angfreq*t)) + 1.0/(6.0*angfreq*angfreq)*(1.0-cos(2.0*angfreq*t)));

src/USER-MISC/angle_cosine_shift.h

@@ -28,7 +28,7 @@ namespace LAMMPS_NS {
 class AngleCosineShift : public Angle {
  public:
   AngleCosineShift(class LAMMPS *);
-  ~AngleCosineShift();
+  virtual ~AngleCosineShift();
   virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
@@ -36,7 +36,7 @@ class AngleCosineShift : public Angle {
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   double *k;
   double *a;
   double *theta;

src/USER-MISC/angle_cosine_shift_exp.cpp

@@ -57,7 +57,7 @@ void AngleCosineShiftExp::compute(int eflag, int vflag)
   int i1,i2,i3,n,type;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double eangle,f1[3],f3[3],ff;
-  double rsq1,rsq2,r1,r2,c,s,cc,ss,a11,a12,a22;
+  double rsq1,rsq2,r1,r2,c,s,a11,a12,a22;
   double exp2,aa,uumin,cccpsss,cssmscc;            
 
   eangle = 0.0;

src/USER-MISC/angle_cosine_shift_exp.h

@@ -26,15 +26,15 @@ namespace LAMMPS_NS {
 class AngleCosineShiftExp : public Angle {
  public:
   AngleCosineShiftExp(class LAMMPS *);
-  ~AngleCosineShiftExp();
-  void compute(int, int);
+  virtual ~AngleCosineShiftExp();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_angle(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int);
 
- private:
+ protected:
   bool *doExpansion;
   double *umin,*a,*opt1;
   double *theta0;

src/USER-MISC/bond_harmonic_shift.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondHarmonicShift : public Bond {
  public:
   BondHarmonicShift(class LAMMPS *);
-  ~BondHarmonicShift();
-  void compute(int, int);
+  virtual ~BondHarmonicShift();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *k,*r0,*r1;
 
   void allocate();

src/USER-MISC/bond_harmonic_shift_cut.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class BondHarmonicShiftCut : public Bond {
  public:
   BondHarmonicShiftCut(class LAMMPS *);
-  ~BondHarmonicShiftCut();
-  void compute(int, int);
+  virtual ~BondHarmonicShiftCut();
+  virtual void compute(int, int);
   void coeff(int, char **);
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, double, int, int);
 
- private:
+ protected:
   double *k,*r0,*r1;
 
   void allocate();

src/angle.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "angle.h"
 #include "atom.h"
+#include "comm.h"
 #include "force.h"
 #include "math_const.h"
 #include "memory.h"
@@ -78,12 +79,12 @@ void Angle::ev_setup(int eflag, int vflag)
   if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
-    memory->create(eatom,maxeatom,"bond:eatom");
+    memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
   }
   if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
-    memory->create(vatom,maxvatom,6,"bond:vatom");
+    memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
   }
 
   // zero accumulators
@@ -216,7 +217,7 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
 
 double Angle::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = comm->nthreads*maxeatom * sizeof(double);
+  bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }

src/angle.h

@@ -20,6 +20,7 @@
 namespace LAMMPS_NS {
 
 class Angle : protected Pointers {
+  friend class ThrOMP;
  public:
   int allocated;
   int *setflag;

src/atom_vec_line.cpp

@@ -17,6 +17,7 @@
 #include "string.h"
 #include "atom_vec_line.h"
 #include "atom.h"
+#include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "force.h"
@@ -83,6 +84,8 @@ void AtomVecLine::grow(int n)
   if (n == 0) nmax += DELTA;
   else nmax = n;
   atom->nmax = nmax;
+  if (nmax < 0 || nmax > MAXSMALLINT)
+    error->one(FLERR,"Per-processor system is too big");
 
   tag = memory->grow(atom->tag,nmax,"atom:tag");
   type = memory->grow(atom->type,nmax,"atom:type");
@@ -90,12 +93,12 @@ void AtomVecLine::grow(int n)
   image = memory->grow(atom->image,nmax,"atom:image");
   x = memory->grow(atom->x,nmax,3,"atom:x");
   v = memory->grow(atom->v,nmax,3,"atom:v");
-  f = memory->grow(atom->f,nmax,3,"atom:f");
+  f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
 
   molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
   rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
   omega = memory->grow(atom->omega,nmax,3,"atom:omega");
-  torque = memory->grow(atom->torque,nmax,3,"atom:torque");
+  torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
   line = memory->grow(atom->line,nmax,"atom:line");
 
   if (atom->nextra_grow)
@@ -1115,12 +1118,12 @@ bigint AtomVecLine::memory_usage()
   if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
   if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
   if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
-  if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
+  if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
 
   if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
   if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
-  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
+  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3);
   if (atom->memcheck("line")) bytes += memory->usage(line,nmax);
 
   bytes += nmax_bonus*sizeof(Bonus);

src/atom_vec_tri.cpp

@@ -18,6 +18,7 @@
 #include "atom_vec_tri.h"
 #include "math_extra.h"
 #include "atom.h"
+#include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "force.h"
@@ -84,6 +85,8 @@ void AtomVecTri::grow(int n)
   if (n == 0) nmax += DELTA;
   else nmax = n;
   atom->nmax = nmax;
+  if (nmax < 0 || nmax > MAXSMALLINT)
+    error->one(FLERR,"Per-processor system is too big");
 
   tag = memory->grow(atom->tag,nmax,"atom:tag");
   type = memory->grow(atom->type,nmax,"atom:type");
@@ -91,12 +94,12 @@ void AtomVecTri::grow(int n)
   image = memory->grow(atom->image,nmax,"atom:image");
   x = memory->grow(atom->x,nmax,3,"atom:x");
   v = memory->grow(atom->v,nmax,3,"atom:v");
-  f = memory->grow(atom->f,nmax,3,"atom:f");
+  f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
 
   molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
   rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
   angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
-  torque = memory->grow(atom->torque,nmax,3,"atom:torque");
+  torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
   tri = memory->grow(atom->tri,nmax,"atom:tri");
 
   if (atom->nextra_grow)
@@ -1547,12 +1550,12 @@ bigint AtomVecTri::memory_usage()
   if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
   if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
   if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
-  if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
+  if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
 
   if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
   if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
   if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
-  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
+  if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3);
   if (atom->memcheck("tri")) bytes += memory->usage(tri,nmax);
 
   bytes += nmax_bonus*sizeof(Bonus);

src/bond.cpp

@@ -14,6 +14,7 @@
 #include "string.h"
 #include "bond.h"
 #include "atom.h"
+#include "comm.h"
 #include "force.h"
 #include "memory.h"
 #include "error.h"
@@ -80,12 +81,12 @@ void Bond::ev_setup(int eflag, int vflag)
   if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
-    memory->create(eatom,maxeatom,"bond:eatom");
+    memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
   }
   if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
-    memory->create(vatom,maxvatom,6,"bond:vatom");
+    memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
   }
 
   // zero accumulators
@@ -198,7 +199,7 @@ void Bond::ev_tally(int i, int j, int nlocal, int newton_bond,
 
 double Bond::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = comm->nthreads*maxeatom * sizeof(double);
+  bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }

src/bond.h

@@ -20,6 +20,7 @@
 namespace LAMMPS_NS {
 
 class Bond : protected Pointers {
+  friend class ThrOMP;
  public:
   int allocated;
   int *setflag;

src/dihedral.cpp

@@ -14,7 +14,7 @@
 #include "math.h"
 #include "dihedral.h"
 #include "atom.h"
-#include "atom.h"
+#include "comm.h"
 #include "force.h"
 #include "pair.h"
 #include "memory.h"
@@ -82,12 +82,12 @@ void Dihedral::ev_setup(int eflag, int vflag)
   if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
-    memory->create(eatom,maxeatom,"bond:eatom");
+    memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
   }
   if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
-    memory->create(vatom,maxvatom,6,"bond:vatom");
+    memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
   }
 
   // zero accumulators
@@ -242,7 +242,7 @@ void Dihedral::ev_tally(int i1, int i2, int i3, int i4,
 
 double Dihedral::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = comm->nthreads*maxeatom * sizeof(double);
+  bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }

src/domain.cpp

@@ -603,7 +603,7 @@ void Domain::minimum_image(double *delta)
    for triclinic, also add/subtract tilt factors in other dims as needed
 ------------------------------------------------------------------------- */
 
-void Domain::closest_image(double *xi, double *xj, double *xjimage)
+void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage)
 {
   double dx,dy,dz;
 

src/domain.h

@@ -93,16 +93,15 @@ class Domain : protected Pointers {
   virtual void set_local_box();
   virtual void reset_box();
   virtual void pbc();
-  virtual void remap(double *, int &);
-  virtual void remap(double *);
-  virtual void remap_near(double *, double *);
-  virtual void unmap(double *, int);
-  virtual void unmap(double *, int, double *);
-  virtual int minimum_image_check(double, double, double);
-  virtual void minimum_image(double &, double &, double &);
-  virtual void minimum_image(double *);
-  virtual void closest_image(double *, double *, double *);
-
+  void remap(double *, int &);
+  void remap(double *);
+  void remap_near(double *, double *);
+  void unmap(double *, int);
+  void unmap(double *, int, double *);
+  int minimum_image_check(double, double, double);
+  void minimum_image(double &, double &, double &);
+  void minimum_image(double *);
+  void closest_image(const double * const, const double * const, double * const);
   void set_lattice(int, char **);
   void add_region(int, char **);
   void delete_region(int, char **);

src/force.cpp

@@ -478,24 +478,61 @@ Improper *Force::new_improper(const char *style, const char *suffix, int &sflag)
    new kspace style 
 ------------------------------------------------------------------------- */
 
-void Force::create_kspace(int narg, char **arg)
+void Force::create_kspace(int narg, char **arg, const char *suffix)
 {
   delete [] kspace_style;
   if (kspace) delete kspace;
 
-  if (strcmp(arg[0],"none") == 0) kspace = NULL;
+  int sflag;
+  kspace = new_kspace(narg,arg,suffix,sflag);
+
+  if (sflag) {
+    char estyle[256];
+    sprintf(estyle,"%s/%s",arg[0],suffix);
+    int n = strlen(estyle) + 1;
+    kspace_style = new char[n];
+    strcpy(kspace_style,estyle);
+  } else {
+    int n = strlen(arg[0]) + 1;
+    kspace_style = new char[n];
+    strcpy(kspace_style,arg[0]);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   generate a kspace class
+------------------------------------------------------------------------- */
+
+KSpace *Force::new_kspace(int narg, char **arg, const char *suffix, int &sflag)
+{
+  if (suffix && lmp->suffix_enable) {
+    sflag = 1;
+    char estyle[256];
+    sprintf(estyle,"%s/%s",arg[0],suffix);
+
+    if (0) return NULL;
 
 #define KSPACE_CLASS
 #define KSpaceStyle(key,Class) \
-  else if (strcmp(arg[0],#key) == 0) kspace = new Class(lmp,narg-1,&arg[1]);
+  else if (strcmp(estyle,#key) == 0) return new Class(lmp,narg-1,&arg[1]);
+#include "style_kspace.h"
+#undef KSpaceStyle
+#undef KSPACE_CLASS
+
+  }
+
+  sflag = 0;
+
+  if (strcmp(arg[0],"none") == 0) return NULL;
+
+#define KSPACE_CLASS
+#define KSpaceStyle(key,Class) \
+  else if (strcmp(arg[0],#key) == 0) return  new Class(lmp,narg-1,&arg[1]);
 #include "style_kspace.h"
 #undef KSPACE_CLASS
 
   else error->all(FLERR,"Invalid kspace style");
-
-  int n = strlen(arg[0]) + 1;
-  kspace_style = new char[n];
-  strcpy(kspace_style,arg[0]);
+  return NULL;
 }
 
 /* ----------------------------------------------------------------------

src/force.h

@@ -86,7 +86,8 @@ class Force : protected Pointers {
   void create_improper(const char *, const char *suffix = NULL);
   class Improper *new_improper(const char *, const char *, int &);
 
-  void create_kspace(int, char **);
+  void create_kspace(int, char **, const char *suffix = NULL);
+  class KSpace *new_kspace(int, char **, const char *, int &);
 
   void set_special(int, char **);
   void bounds(char *, int, int &, int &);

src/improper.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "improper.h"
 #include "atom.h"
+#include "comm.h"
 #include "force.h"
 #include "memory.h"
 #include "error.h"
@@ -76,12 +77,12 @@ void Improper::ev_setup(int eflag, int vflag)
   if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
-    memory->create(eatom,maxeatom,"bond:eatom");
+    memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
   }
   if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
-    memory->create(vatom,maxvatom,6,"bond:vatom");
+    memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
   }
 
   // zero accumulators
@@ -236,7 +237,7 @@ void Improper::ev_tally(int i1, int i2, int i3, int i4,
 
 double Improper::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = comm->nthreads*maxeatom * sizeof(double);
+  bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }

src/improper.h

@@ -20,6 +20,7 @@
 namespace LAMMPS_NS {
 
 class Improper : protected Pointers {
+  friend class ThrOMP;
  public:
   int allocated;
   int *setflag;

src/input.cpp

@@ -1022,7 +1022,8 @@ void Input::improper_style()
 
 void Input::kspace_modify()
 {
-  if (force->kspace == NULL) error->all(FLERR,"KSpace style has not yet been set");
+  if (force->kspace == NULL) 
+    error->all(FLERR,"KSpace style has not yet been set");
   force->kspace->modify_params(narg,arg);
 }
 
@@ -1030,7 +1031,7 @@ void Input::kspace_modify()
 
 void Input::kspace_style()
 {
-  force->create_kspace(narg,arg);
+  force->create_kspace(narg,arg,lmp->suffix);
 }
 
 /* ---------------------------------------------------------------------- */

src/integrate.cpp

@@ -26,6 +26,7 @@ Integrate::Integrate(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 {
   elist_global = elist_atom = NULL;
   vlist_global = vlist_atom = NULL;
+  external_force_clear = 0;
 }
 
 /* ---------------------------------------------------------------------- */

src/integrate.h

@@ -33,6 +33,7 @@ class Integrate : protected Pointers {
  protected:
   int eflag,vflag;                  // flags for energy/virial computation
   int virial_style;                 // compute virial explicitly or implicitly
+  int external_force_clear;         // clear forces locally or externally
 
   int nelist_global,nelist_atom;    // # of PE,virial computes to check
   int nvlist_global,nvlist_atom;

src/kspace.h

@@ -19,6 +19,7 @@
 namespace LAMMPS_NS {
 
 class KSpace : protected Pointers {
+  friend class ThrOMP;
  public:
   double energy;
   double virial[6];

src/math_const.h

@@ -19,6 +19,8 @@ namespace LAMMPS_NS {
 namespace MathConst {
   static const double THIRD  = 1.0/3.0;
   static const double MY_PI  = 3.14159265358979323846; // pi
+  static const double MY_2PI = 6.28318530717958647692; // 2pi
+  static const double MY_3PI = 9.42477796076937971538; // 3pi
   static const double MY_PI2 = 1.57079632679489661923; // pi/2
   static const double MY_PI4 = 0.78539816339744830962; // pi/4
   static const double MY_PIS = 1.77245385090551602729; // sqrt(pi)

src/min.cpp

@@ -157,7 +157,8 @@ void Min::init()
   
   if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
     if (comm->me == 0) 
-      error->warning(FLERR,"Resetting reneighboring criteria during minimization");
+      error->warning(FLERR,
+		     "Resetting reneighboring criteria during minimization");
   }
 
   neighbor->every = 1;
@@ -227,9 +228,11 @@ void Min::setup()
   // remove these restriction eventually
 
   if (nextra_global && searchflag == 0)
-    error->all(FLERR,"Cannot use a damped dynamics min style with fix box/relax");
+    error->all(FLERR,
+	       "Cannot use a damped dynamics min style with fix box/relax");
   if (nextra_atom && searchflag == 0)
-    error->all(FLERR,"Cannot use a damped dynamics min style with per-atom DOF");
+    error->all(FLERR,
+	       "Cannot use a damped dynamics min style with per-atom DOF");
 
   // atoms may have migrated in comm->exchange()
 
@@ -510,6 +513,8 @@ void Min::force_clear()
 {
   int i;
 
+  if (external_force_clear) return;
+
   // clear global force array
   // nall includes ghosts only if either newton flag is set
 

src/min.h

@@ -49,6 +49,7 @@ class Min : protected Pointers {
  protected:
   int eflag,vflag;            // flags for energy/virial computation
   int virial_style;           // compute virial explicitly or implicitly
+  int external_force_clear;   // clear forces locally or externally
 
   double dmax;                // max dist to move any atom in one step
   int linestyle;              // 0 = backtrack, 1 = quadratic

src/neigh_full.cpp

@@ -108,7 +108,7 @@ void Neighbor::full_nsq(NeighList *list)
 
 /* ----------------------------------------------------------------------
    N^2 search for all neighbors
-   include neighbors of ghost atoms
+   include neighbors of ghost atoms, but no "special neighbors" for ghosts
    every neighbor pair appears in list of both atoms i and j
 ------------------------------------------------------------------------- */
 
@@ -187,12 +187,9 @@ void Neighbor::full_nsq_ghost(NeighList *list)
 	dely = ytmp - x[j][1];
 	delz = ztmp - x[j][2];
 	rsq = delx*delx + dely*dely + delz*delz;
-	if (rsq <= cutneighghostsq[itype][jtype]) {
-	  if (molecular) {
-	    which = find_special(special[i],nspecial[i],tag[j]);
-	    if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
-	  } else neighptr[n++] = j;
-	}
+
+	if (rsq <= cutneighghostsq[itype][jtype])
+	  neighptr[n++] = j;
       }
     }
 
@@ -201,7 +198,8 @@ void Neighbor::full_nsq_ghost(NeighList *list)
     numneigh[i] = n;
     npnt += n;
     if (n > oneatom || npnt >= pgsize)
-      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+      error->one(FLERR,
+		 "Neighbor list overflow, boost neigh_modify one or page");
   }
 
   list->inum = atom->nlocal;
@@ -295,7 +293,8 @@ void Neighbor::full_bin(NeighList *list)
     numneigh[i] = n;
     npnt += n;
     if (n > oneatom || npnt >= pgsize)
-      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+      error->one(FLERR,
+		 "Neighbor list overflow, boost neigh_modify one or page");
   }
 
   list->inum = inum;
@@ -304,7 +303,7 @@ void Neighbor::full_bin(NeighList *list)
 
 /* ----------------------------------------------------------------------
    binned neighbor list construction for all neighbors
-   include neighbors of ghost atoms
+   include neighbors of ghost atoms, but no "special neighbors" for ghosts
    every neighbor pair appears in list of both atoms i and j
 ------------------------------------------------------------------------- */
 
@@ -407,13 +406,9 @@ void Neighbor::full_bin_ghost(NeighList *list)
 	  dely = ytmp - x[j][1];
 	  delz = ztmp - x[j][2];
 	  rsq = delx*delx + dely*dely + delz*delz;
-	
-	  if (rsq <= cutneighghostsq[itype][jtype]) {
-	    if (molecular) {
-	      which = find_special(special[i],nspecial[i],tag[j]);
-	      if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
-	    } else neighptr[n++] = j;
-	  }
+
+	  if (rsq <= cutneighghostsq[itype][jtype])
+	    neighptr[n++] = j;
 	}
       }
     }
@@ -423,7 +418,8 @@ void Neighbor::full_bin_ghost(NeighList *list)
     numneigh[i] = n;
     npnt += n;
     if (n > oneatom || npnt >= pgsize)
-      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+      error->one(FLERR,
+		 "Neighbor list overflow, boost neigh_modify one or page");
   }
 
   list->inum = atom->nlocal;

src/neigh_list.cpp

@@ -15,17 +15,12 @@
 #include "neigh_list.h"
 #include "atom.h"
 #include "comm.h"
+#include "update.h"
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "memory.h"
 #include "error.h"
 
-
-
-#include "update.h"
-
-
-
 using namespace LAMMPS_NS;
 
 #define PGDELTA 1
@@ -255,6 +250,8 @@ void NeighList::print_attributes()
   printf("  %d = occasional\n",rq->occasional);
   printf("  %d = dnum\n",rq->dnum);
   printf("  %d = ghost\n",rq->ghost);
+  printf("  %d = cudable\n",rq->cudable);
+  printf("  %d = omp\n",rq->omp);
   printf("  %d = copy\n",rq->copy);
   printf("  %d = skip\n",rq->skip);
   printf("  %d = otherlist\n",rq->otherlist);

src/neigh_request.cpp

@@ -44,6 +44,8 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   // default is encode special bond flags
   // default is no auxiliary floating point values
   // default is no neighbors of ghosts
+  // default is no CUDA neighbor list build
+  // default is no multi-threaded neighbor list build
 
   occasional = 0;
   newton = 0;
@@ -51,6 +53,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
   dnum = 0;
   ghost = 0;
   cudable = 0;
+  omp = 0;
 
   // default is no copy or skip
 
@@ -101,6 +104,7 @@ int NeighRequest::identical(NeighRequest *other)
   if (dnum != other->dnum) same = 0;
   if (ghost != other->ghost) same = 0;
   if (cudable != other->cudable) same = 0;
+  if (omp != other->omp) same = 0;
 
   if (copy != other->copy) same = 0;
   if (same_skip(other) == 0) same = 0;
@@ -129,6 +133,7 @@ int NeighRequest::same_kind(NeighRequest *other)
   if (newton != other->newton) same = 0;
   if (ghost != other->ghost) same = 0;
   if (cudable != other->cudable) same = 0;
+  if (omp != other->omp) same = 0;
 
   return same;
 }
@@ -178,4 +183,5 @@ void NeighRequest::copy_request(NeighRequest *other)
   dnum = other->dnum;
   ghost = other->ghost;
   cudable = other->cudable;
+  omp = other->omp;
 }

src/neigh_request.h

@@ -76,6 +76,10 @@ class NeighRequest : protected Pointers {
 
   int cudable;
 
+  // 1 if using multi-threaded neighbor list build
+
+  int omp;
+
   // set by neighbor and pair_hybrid after all requests are made
   // these settings do not change kind value
 

src/neighbor.cpp

@@ -1698,9 +1698,8 @@ int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz)
    return 1 if should be excluded, 0 if included
 ------------------------------------------------------------------------- */
 
-int Neighbor::exclusion(int i, int j, int itype, int jtype, 
-			int *mask, int *molecule)
-{
+int Neighbor::exclusion(int i, int j, int itype, int jtype,
+			int *mask, int *molecule) const {
   int m;
 
   if (nex_type && ex_type[itype][jtype]) return 1;

src/neighbor.h

@@ -164,7 +164,7 @@ class Neighbor : protected Pointers {
   int coord2bin(double *);              // mapping atom coord to a bin
   int coord2bin(double *, int &, int &, int&); // ditto
 
-  int exclusion(int, int, int, int, int *, int *);  // test for pair exclusion
+  int exclusion(int, int, int, int, int *, int *) const;  // test for pair exclusion
   virtual void choose_build(int, class NeighRequest *);
   void choose_stencil(int, class NeighRequest *);
 

src/pair.cpp

@@ -281,12 +281,12 @@ void Pair::ev_setup(int eflag, int vflag)
   if (eflag_atom && atom->nmax > maxeatom) {
     maxeatom = atom->nmax;
     memory->destroy(eatom);
-    memory->create(eatom,maxeatom,"pair:eatom");
+    memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom");
   }
   if (vflag_atom && atom->nmax > maxvatom) {
     maxvatom = atom->nmax;
     memory->destroy(vatom);
-    memory->create(vatom,maxvatom,6,"pair:vatom");
+    memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom");
   }
 
   // zero accumulators
@@ -1151,7 +1151,7 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
 
 double Pair::memory_usage()
 {
-  double bytes = maxeatom * sizeof(double);
-  bytes += maxvatom*6 * sizeof(double);
+  double bytes = comm->nthreads*maxeatom * sizeof(double);
+  bytes += comm->nthreads*maxvatom*6 * sizeof(double);
   return bytes;
 }

src/pair.h

@@ -20,6 +20,7 @@ namespace LAMMPS_NS {
 
 class Pair : protected Pointers {
   friend class BondQuartic;
+  friend class BondQuarticOMP;
   friend class DihedralCharmm;
   friend class DihedralCharmmOMP;
   friend class FixGPU;

src/respa.cpp

@@ -599,6 +599,8 @@ void Respa::force_clear(int newtonflag)
 {
   int i;
 
+  if (external_force_clear) return;
+
   // clear global force array
   // nall includes ghosts only if newton flag is set
 

src/thermo.h

@@ -38,13 +38,13 @@ class Thermo : protected Pointers {
   int evaluate_keyword(char *, double *);
 
  private:
-  int me;
-
   char *line;
-  int nfield,nfield_initial;
   char **keyword;
   int *vtype;
 
+  int nfield,nfield_initial;
+  int me;
+
   char **format,**format_user;
   char *format_float_one_def,*format_float_multi_def;
   char *format_int_one_def,*format_int_multi_def;

src/verlet.cpp

@@ -308,6 +308,8 @@ void Verlet::force_clear()
 {
   int i;
 
+  if (external_force_clear) return;
+
   // clear force on all particles
   // if either newton flag is set, also include ghosts