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

src/ASPHERE/pair_gayberne.cpp

@@ -33,8 +33,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
-
 /* ---------------------------------------------------------------------- */
 
 PairGayBerne::PairGayBerne(LAMMPS *lmp) : Pair(lmp)

src/ASPHERE/pair_gayberne.h

@@ -39,6 +39,8 @@ class PairGayBerne : public Pair {
   void read_restart_settings(FILE *);
 
  protected:
+  enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
+
   double cut_global;
   double **cut;
 

src/ASPHERE/pair_resquared.cpp

@@ -33,8 +33,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
-
 /* ---------------------------------------------------------------------- */
 
 PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp), 

src/ASPHERE/pair_resquared.h

@@ -39,6 +39,8 @@ class PairRESquared : public Pair {
   void read_restart_settings(FILE *);
 
  protected:
+  enum{SPHERE_SPHERE,SPHERE_ELLIPSE,ELLIPSE_SPHERE,ELLIPSE_ELLIPSE};
+
   double cut_global;
   double **cut;
 

src/CLASS2/dihedral_class2.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class DihedralClass2 : public Dihedral {
  public:
   DihedralClass2(class LAMMPS *);
-  ~DihedralClass2();
-  void compute(int, int);
+  virtual ~DihedralClass2();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k1,*k2,*k3;
   double *phi1,*phi2,*phi3;
   double *mbt_f1,*mbt_f2,*mbt_f3,*mbt_r0;

src/CLASS2/pair_lj_class2.h

@@ -28,7 +28,7 @@ class PairLJClass2 : public Pair {
  public:
   PairLJClass2(class LAMMPS *);
   virtual ~PairLJClass2();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);

src/CLASS2/pair_lj_class2_coul_cut.h

@@ -28,7 +28,7 @@ class PairLJClass2CoulCut : public Pair {
  public:
   PairLJClass2CoulCut(class LAMMPS *);
   virtual ~PairLJClass2CoulCut();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();

src/CLASS2/pair_lj_class2_coul_long.h

@@ -28,7 +28,7 @@ class PairLJClass2CoulLong : public Pair {
  public:
   PairLJClass2CoulLong(class LAMMPS *);
   virtual ~PairLJClass2CoulLong();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();

src/COLLOID/pair_colloid.cpp

@@ -30,8 +30,6 @@
 
 using namespace LAMMPS_NS;
 
-enum{SMALL_SMALL,SMALL_LARGE,LARGE_LARGE};
-
 /* ---------------------------------------------------------------------- */
 
 PairColloid::PairColloid(LAMMPS *lmp) : Pair(lmp) {}

src/COLLOID/pair_colloid.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairColloid : public Pair {
  public:
   PairColloid(class LAMMPS *);
-  ~PairColloid();
-  void compute(int, int);
+  virtual ~PairColloid();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -38,7 +38,9 @@ class PairColloid : public Pair {
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
 
- private:
+  enum {SMALL_SMALL,SMALL_LARGE,LARGE_LARGE};
+
+ protected:
   double cut_global;
   double **cut;
   double **a12,**d1,**d2,**diameter,**a1,**a2,**offset;

src/COLLOID/pair_lubricate.h

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

src/COLLOID/pair_yukawa_colloid.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairYukawaColloid : public PairYukawa {
  public:
   PairYukawaColloid(class LAMMPS *);
-  ~PairYukawaColloid() {}
-  void compute(int, int);
+  virtual ~PairYukawaColloid() {}
+  virtual void compute(int, int);
   void init_style();
   double init_one(int, int);
   double single(int, int, int, int, double, double, double, double &);

src/DIPOLE/pair_dipole_cut.h

@@ -26,9 +26,9 @@ namespace LAMMPS_NS {
 
 class PairDipoleCut : public Pair {
  public:
-    PairDipoleCut(class LAMMPS *);
-  ~PairDipoleCut();
-  void compute(int, int);
+  PairDipoleCut(class LAMMPS *);
+  virtual ~PairDipoleCut();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -38,7 +38,7 @@ class PairDipoleCut : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
 
- private:
+ protected:
   double cut_lj_global,cut_coul_global;
   double **cut_lj,**cut_ljsq;
   double **cut_coul,**cut_coulsq;

src/GRANULAR/fix_pour.h

@@ -26,8 +26,11 @@ namespace LAMMPS_NS {
 
 class FixPour : public Fix {
   friend class PairGranHertzHistory;
+  friend class PairGranHertzHistoryOMP;
   friend class PairGranHooke;
+  friend class PairGranHookeOMP;
   friend class PairGranHookeHistory;
+  friend class PairGranHookeHistoryOMP;
   friend class PairGranHookeCuda;
 
  public:

src/GRANULAR/fix_wall_gran.cpp

@@ -221,8 +221,12 @@ void FixWallGran::init()
     pairstyle = HOOKE;
   else if (force->pair_match("gran/hooke/history",1))
     pairstyle = HOOKE_HISTORY;
+  else if (force->pair_match("gran/hooke/history/omp",1))
+    pairstyle = HOOKE_HISTORY;
   else if (force->pair_match("gran/hertz/history",1))
     pairstyle = HERTZ_HISTORY;
+  else if (force->pair_match("gran/hertz/history/omp",1))
+    pairstyle = HERTZ_HISTORY;
   else error->all(FLERR,"Fix wall/gran is incompatible with Pair style");
 }
 

src/GRANULAR/pair_gran_hertz_history.h

@@ -27,7 +27,7 @@ namespace LAMMPS_NS {
 class PairGranHertzHistory : public PairGranHookeHistory {
  public:
   PairGranHertzHistory(class LAMMPS *);
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
 };
 

src/GRANULAR/pair_gran_hooke.h

@@ -27,7 +27,7 @@ namespace LAMMPS_NS {
 class PairGranHooke : public PairGranHookeHistory {
  public:
   PairGranHooke(class LAMMPS *);
-  void compute(int, int);
+  virtual void compute(int, int);
 };
 
 }

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -46,6 +46,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
   no_virial_fdotr_compute = 1;
   history = 1;
   fix_history = NULL;
+  suffix = NULL;
 
   laststep = -1;
 }
@@ -55,6 +56,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
 PairGranHookeHistory::~PairGranHookeHistory()
 {
   if (fix_history) modify->delete_fix("SHEAR_HISTORY");
+  if (suffix) delete[] suffix;
 
   if (allocated) {
     memory->destroy(setflag);
@@ -393,7 +395,7 @@ void PairGranHookeHistory::init_style()
     fixarg[0] = (char *) "SHEAR_HISTORY";
     fixarg[1] = (char *) "all";
     fixarg[2] = (char *) "SHEAR_HISTORY";
-    modify->add_fix(3,fixarg);
+    modify->add_fix(3,fixarg,suffix);
     delete [] fixarg;
     fix_history = (FixShearHistory *) modify->fix[modify->nfix-1];
     fix_history->pair = this;

src/GRANULAR/pair_gran_hooke_history.h

@@ -49,7 +49,8 @@ class PairGranHookeHistory : public Pair {
 
   bigint laststep;
   class FixShearHistory *fix_history;
-  int shearupdate;
+
+  char *suffix;
 
   double *onerad_dynamic,*onerad_frozen;
   double *maxrad_dynamic,*maxrad_frozen;

src/KSPACE/pair_born_coul_long.h

@@ -28,7 +28,7 @@ class PairBornCoulLong : public Pair {
  public:
   PairBornCoulLong(class LAMMPS *);
   virtual ~PairBornCoulLong();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();

src/KSPACE/pair_lj_charmm_coul_long.h

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

src/MANYBODY/fix_qeq_comb.h

@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
 class FixQEQComb : public Fix {
  public:
   FixQEQComb(class LAMMPS *, int, char **);
-  ~FixQEQComb();
+  virtual ~FixQEQComb();
   int setmask();
-  void init();
+  virtual void init();
   void setup(int);
-  void post_force(int);
+  virtual void post_force(int);
   void post_force_respa(int,int,int);
   double memory_usage();
 
- private:
+ protected:
   int me,firstflag;
   double precision;
   int nlevels_respa;

src/MANYBODY/pair_adp.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairADP : public Pair {
  public:
   PairADP(class LAMMPS *);
-  ~PairADP();
-  void compute(int, int);
+  virtual ~PairADP();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -40,7 +40,7 @@ class PairADP : public Pair {
   void unpack_reverse_comm(int, int *, double *);
   double memory_usage();
   
- private:
+ protected:
   int nmax;                   // allocated size of per-atom arrays
   double cutforcesq,cutmax;
 

src/MANYBODY/pair_eam.h

@@ -46,7 +46,7 @@ class PairEAM : public Pair {
 
   PairEAM(class LAMMPS *);
   virtual ~PairEAM();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   virtual void coeff(int, char **);
   void init_style();

src/MANYBODY/pair_eim.h

@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
 class PairEIM : public Pair {
  public:
   PairEIM(class LAMMPS *);
-  ~PairEIM();
-  void compute(int, int);
+  virtual ~PairEIM();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -41,7 +41,7 @@ class PairEIM : public Pair {
   void unpack_reverse_comm(int, int *, double *);
   double memory_usage();
 
- private:
+ protected:
   double **cutforcesq,cutmax;
   int nmax;
   double *rho,*fp;

src/MANYBODY/pair_sw.h

@@ -27,14 +27,14 @@ namespace LAMMPS_NS {
 class PairSW : public Pair {
  public:
   PairSW(class LAMMPS *);
-  ~PairSW();
-  void compute(int, int);
+  virtual ~PairSW();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
   void init_style();
 
- private:
+ protected:
   struct Param {
     double epsilon,sigma;
     double littlea,lambda,gamma,costheta;

src/MOLECULE/angle_hybrid.cpp

@@ -207,8 +207,10 @@ void AngleHybrid::settings(int narg, char **arg)
   // one exception is 1st arg of style "table", which is non-numeric
   // need a better way to skip these exceptions
 
+  int dummy;
   nstyles = 0;
   i = 0;
+
   while (i < narg) {
     for (m = 0; m < nstyles; m++)
       if (strcmp(arg[i],keywords[m]) == 0) 
@@ -217,7 +219,7 @@ void AngleHybrid::settings(int narg, char **arg)
       error->all(FLERR,"Angle style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[i],"none") == 0) 
       error->all(FLERR,"Angle style hybrid cannot have none as an argument");
-    styles[nstyles] = force->new_angle(arg[i]);
+    styles[nstyles] = force->new_angle(arg[i],lmp->suffix,dummy);
     keywords[nstyles] = new char[strlen(arg[i])+1];
     strcpy(keywords[nstyles],arg[i]);
     istyle = i;
@@ -321,14 +323,14 @@ void AngleHybrid::read_restart(FILE *fp)
 
   allocate();
   
-  int n;
+  int n,dummy;
   for (int m = 0; m < nstyles; m++) {
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
-    styles[m] = force->new_angle(keywords[m]);
+    styles[m] = force->new_angle(keywords[m],lmp->suffix,dummy);
   }
 }
 

src/MOLECULE/dihedral_charmm.h

@@ -28,14 +28,14 @@ namespace LAMMPS_NS {
 class DihedralCharmm : public Dihedral {
  public:
   DihedralCharmm(class LAMMPS *);
-  ~DihedralCharmm();
-  void compute(int, int);
+  virtual ~DihedralCharmm();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k,*weight,*cos_shift,*sin_shift;
   int *multiplicity,*shift;
   double **lj14_1,**lj14_2,**lj14_3,**lj14_4;

src/MOLECULE/dihedral_harmonic.h

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

src/MOLECULE/dihedral_helix.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class DihedralHelix : public Dihedral {
  public:
   DihedralHelix(class LAMMPS *);
-  ~DihedralHelix();
-  void compute(int, int);
+  virtual ~DihedralHelix();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *aphi,*bphi,*cphi;
 
   void allocate();

src/MOLECULE/dihedral_hybrid.cpp

@@ -165,16 +165,19 @@ void DihedralHybrid::settings(int narg, char **arg)
   styles = new Dihedral*[nstyles];
   keywords = new char*[nstyles];
 
+  int dummy;
+
   for (int m = 0; m < nstyles; m++) {
     for (int i = 0; i < m; i++)
       if (strcmp(arg[m],arg[i]) == 0) 
 	error->all(FLERR,"Dihedral style hybrid cannot use "
 		   "same dihedral style twice");
     if (strcmp(arg[m],"hybrid") == 0) 
-      error->all(FLERR,"Dihedral style hybrid cannot have hybrid as an argument");
+      error->all(FLERR,
+		 "Dihedral style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[m],"none") == 0) 
       error->all(FLERR,"Dihedral style hybrid cannot have none as an argument");
-    styles[m] = force->new_dihedral(arg[m]);
+    styles[m] = force->new_dihedral(arg[m],lmp->suffix,dummy);
     keywords[m] = new char[strlen(arg[m])+1];
     strcpy(keywords[m],arg[m]);
   }
@@ -269,14 +272,14 @@ void DihedralHybrid::read_restart(FILE *fp)
 
   allocate();
 
-  int n;
+  int n,dummy;
   for (int m = 0; m < nstyles; m++) {
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
-    styles[m] = force->new_dihedral(keywords[m]);
+    styles[m] = force->new_dihedral(keywords[m],lmp->suffix,dummy);
   }
 }
 

src/MOLECULE/dihedral_multi_harmonic.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class DihedralMultiHarmonic : public Dihedral {
  public:
   DihedralMultiHarmonic(class LAMMPS *);
-  ~DihedralMultiHarmonic();
-  void compute(int, int);
+  virtual ~DihedralMultiHarmonic();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *a1,*a2,*a3,*a4,*a5;
 
   void allocate();

src/MOLECULE/dihedral_opls.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class DihedralOPLS : public Dihedral {
  public:
   DihedralOPLS(class LAMMPS *);
-  ~DihedralOPLS();
-  void compute(int, int);
+  virtual ~DihedralOPLS();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   double *k1,*k2,*k3,*k4;
 
   void allocate();

src/MOLECULE/improper_hybrid.cpp

@@ -165,15 +165,19 @@ void ImproperHybrid::settings(int narg, char **arg)
   styles = new Improper*[nstyles];
   keywords = new char*[nstyles];
 
+  int dummy;
+
   for (int m = 0; m < nstyles; m++) {
     for (int i = 0; i < m; i++)
       if (strcmp(arg[m],arg[i]) == 0) 
-	error->all(FLERR,"Improper style hybrid cannot use same improper style twice");
+	error->all(FLERR,
+		   "Improper style hybrid cannot use same improper style twice");
     if (strcmp(arg[m],"hybrid") == 0) 
-      error->all(FLERR,"Improper style hybrid cannot have hybrid as an argument");
+      error->all(FLERR,
+		 "Improper style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[m],"none") == 0) 
       error->all(FLERR,"Improper style hybrid cannot have none as an argument");
-    styles[m] = force->new_improper(arg[m]);
+    styles[m] = force->new_improper(arg[m],lmp->suffix,dummy);
     keywords[m] = new char[strlen(arg[m])+1];
     strcpy(keywords[m],arg[m]);
   }
@@ -256,14 +260,14 @@ void ImproperHybrid::read_restart(FILE *fp)
 
   allocate();
   
-  int n;
+  int n,dummy;
   for (int m = 0; m < nstyles; m++) {
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
-    styles[m] = force->new_improper(keywords[m]);
+    styles[m] = force->new_improper(keywords[m],lmp->suffix,dummy);
   }
 }
 

src/MOLECULE/pair_hbond_dreiding_morse.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairHbondDreidingMorse : public PairHbondDreidingLJ {
  public:
   PairHbondDreidingMorse(class LAMMPS *);
-  ~PairHbondDreidingMorse() {}
-  void compute(int, int);
+  virtual ~PairHbondDreidingMorse() {};
+  virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
   double single(int, int, int, int, double, double, double, double &);

src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.h

@@ -27,7 +27,7 @@ namespace LAMMPS_NS {
 class PairLJCharmmCoulCharmmImplicit : public PairLJCharmmCoulCharmm {
  public:
   PairLJCharmmCoulCharmmImplicit(class LAMMPS *);
-  void compute(int, int);
+  virtual void compute(int, int);
   double single(int, int, int, int, double, double, double, double &);
 };
 

src/PERI/fix_peri_neigh.cpp

@@ -260,8 +260,8 @@ void FixPeriNeigh::setup(int vflag)
   double **x0 = atom->x0;
   double half_lc = 0.5*(domain->lattice->xlattice);
   double vfrac_scale;
-  PairPeriLPS *pairlps = dynamic_cast<PairPeriLPS*>(anypair);
-  PairPeriPMB *pairpmb = dynamic_cast<PairPeriPMB*>(anypair);
+  PairPeriLPS *pairlps = static_cast<PairPeriLPS*>(anypair);
+  PairPeriPMB *pairpmb = static_cast<PairPeriPMB*>(anypair);
   for (i = 0; i < nlocal; i++) {
     double xtmp0 = x0[i][0];
     double ytmp0 = x0[i][1];

src/PERI/fix_peri_neigh.h

@@ -26,7 +26,9 @@ namespace LAMMPS_NS {
 
 class FixPeriNeigh : public Fix {
   friend class PairPeriPMB;
+  friend class PairPeriPMBOMP;
   friend class PairPeriLPS;
+  friend class PairPeriLPSOMP;
   friend class ComputeDamageAtom;
 
  public:

src/PERI/pair_peri_lps.cpp

@@ -174,7 +174,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
 
         double kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) /
 	  (3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]);
-        rk = (kshort * vfrac[j]) * (dr / sqrt(cutsq[itype][jtype]));
+        rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]);
 
         if (r > 0.0) fpair = -(rk/r);
         else fpair = 0.0;
@@ -214,9 +214,9 @@ void PairPeriLPS::compute(int eflag, int vflag)
     comm->forward_comm_fix(modify->fix[ifix_peri]);
 
   // Volume-dependent part of the energy
-  for (i = 0; i < nlocal; i++) {   
-    itype = type[i];
-    if (eflag) {
+  if (eflag) {
+    for (i = 0; i < nlocal; i++) {   
+      itype = type[i];
       if (eflag_global)
 	eng_vdwl += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]);
       if (eflag_atom)
@@ -266,7 +266,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
       delz0 = ztmp0 - x0[j][2];						
       if (periodic) domain->minimum_image(delx0,dely0,delz0);   
       jtype = type[j];
-      delta = sqrt(cutsq[itype][jtype]);
+      delta = cut[itype][jtype];
       r = sqrt(rsq);
       dr = r - r0[i][jj];
 
@@ -316,8 +316,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
       if (first)
          s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch);
       else
-         s0_new[i] = MAX(s0_new[i],
-                         s00[itype][jtype] - (alpha[itype][jtype] * stretch));
+         s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - (alpha[itype][jtype] * stretch));
 
       first = false;
     }
@@ -372,11 +371,11 @@ void PairPeriLPS::coeff(int narg, char **arg)
   force->bounds(arg[0],atom->ntypes,ilo,ihi);
   force->bounds(arg[1],atom->ntypes,jlo,jhi);
 
-  double bulkmodulus_one = atof(arg[2]);	
-  double shearmodulus_one = atof(arg[3]);		
-  double cut_one = atof(arg[4]);			
-  double s00_one = atof(arg[5]);			
-  double alpha_one = atof(arg[6]);			
+  double bulkmodulus_one = atof(arg[2]);
+  double shearmodulus_one = atof(arg[3]);
+  double cut_one = atof(arg[4]);
+  double s00_one = atof(arg[5]);
+  double alpha_one = atof(arg[6]);
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
@@ -402,10 +401,11 @@ double PairPeriLPS::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
 
-  bulkmodulus[j][i] = bulkmodulus[i][j];		
-  shearmodulus[j][i] = shearmodulus[i][j];		
-  s00[j][i] = s00[i][j];		
-  alpha[j][i] = alpha[i][j];		
+  bulkmodulus[j][i] = bulkmodulus[i][j];
+  shearmodulus[j][i] = shearmodulus[i][j];
+  s00[j][i] = s00[i][j];
+  alpha[j][i] = alpha[i][j];
+  cut[j][i] = cut[i][j];
 
   return cut[i][j];
 }
@@ -689,7 +689,7 @@ void PairPeriLPS::compute_dilatation()
       if (fabs(dr) < 2.2204e-016) dr = 0.0;
 
       jtype = type[j];
-      delta = sqrt(cutsq[itype][jtype]);
+      delta = cut[itype][jtype];
 
       // scale vfrac[j] if particle j near the horizon
 

src/PERI/pair_peri_lps.h

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

src/PERI/pair_peri_pmb.cpp

@@ -161,7 +161,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
         dr = r - d_ij;
 
         rk = (15.0 * kspring[itype][jtype] * vfrac[j]) * 
-	  (dr / sqrt(cutsq[itype][jtype]));
+	  (dr / cut[itype][jtype]);
         if (r > 0.0) fpair = -(rk/r);
         else fpair = 0.0;
 
@@ -182,7 +182,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
 
   // grow bond forces array if necessary
 
-  if (nlocal > nmax) {
+  if (atom->nmax > nmax) {
     memory->destroy(s0_new);
     nmax = atom->nmax;
     memory->create(s0_new,nmax,"pair:s0_new");
@@ -223,7 +223,7 @@ void PairPeriPMB::compute(int eflag, int vflag)
       if (periodic) domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-      delta = sqrt(cutsq[itype][jtype]);
+      delta = cut[itype][jtype];
       r = sqrt(rsq);
       dr = r - r0[i][jj];
 
@@ -345,6 +345,7 @@ double PairPeriPMB::init_one(int i, int j)
   kspring[j][i] = kspring[i][j];
   alpha[j][i] = alpha[i][j];
   s00[j][i] = s00[i][j];
+  cut[j][i] = cut[i][j];
 
   return cut[i][j];
 }

src/PERI/pair_peri_pmb.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairPeriPMB : public Pair {
  public:
   PairPeriPMB(class LAMMPS *);
-  ~PairPeriPMB();
-  void compute(int, int);
+  virtual ~PairPeriPMB();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -38,7 +38,7 @@ class PairPeriPMB : public Pair {
   void write_restart_settings(FILE *) {}
   void read_restart_settings(FILE *) {}
   double single(int, int, int, int, double, double, double, double &);
-  double memory_usage();
+  virtual double memory_usage();
 
  protected:
   int ifix_peri;

src/USER-AWPMD/atom_vec_wavepacket.cpp

@@ -84,7 +84,7 @@ void AtomVecWavepacket::grow(int n)
   spin = memory->grow(atom->spin,nmax,"atom:spin");
   eradius = memory->grow(atom->eradius,nmax,"atom:eradius");
   ervel = memory->grow(atom->ervel,nmax,"atom:ervel");
-  erforce = memory->grow(atom->erforce,nmax,"atom:erforce");
+  erforce = memory->grow(atom->erforce,nmax*comm->nthreads,"atom:erforce");
 
   cs = memory->grow(atom->cs,2*nmax,"atom:cs");
   csforce = memory->grow(atom->csforce,2*nmax,"atom:csforce");
@@ -997,7 +997,7 @@ bigint AtomVecWavepacket::memory_usage()
   if (atom->memcheck("spin")) bytes += memory->usage(spin,nmax);
   if (atom->memcheck("eradius")) bytes += memory->usage(eradius,nmax);
   if (atom->memcheck("ervel")) bytes += memory->usage(ervel,nmax);
-  if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax);
+  if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax*comm->nthreads);
 
   if (atom->memcheck("ervelforce")) bytes += memory->usage(ervelforce,nmax);
   if (atom->memcheck("cs")) bytes += memory->usage(cs,2*nmax);

src/USER-EFF/atom_vec_electron.cpp

@@ -78,7 +78,7 @@ void AtomVecElectron::grow(int n)
   spin = memory->grow(atom->spin,nmax,"atom:spin");
   eradius = memory->grow(atom->eradius,nmax,"atom:eradius");
   ervel = memory->grow(atom->ervel,nmax,"atom:ervel");
-  erforce = memory->grow(atom->erforce,nmax,"atom:erforce");
+  erforce = memory->grow(atom->erforce,nmax*comm->nthreads,"atom:erforce");
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
@@ -845,7 +845,7 @@ bigint AtomVecElectron::memory_usage()
   if (atom->memcheck("spin")) bytes += memory->usage(spin,nmax);
   if (atom->memcheck("eradius")) bytes += memory->usage(eradius,nmax);
   if (atom->memcheck("ervel")) bytes += memory->usage(ervel,nmax);
-  if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax);
+  if (atom->memcheck("erforce")) bytes += memory->usage(erforce,nmax*comm->nthreads);
   
   return bytes;
 }

src/USER-MISC/dihedral_cosine_shift_exp.h

@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
 class DihedralCosineShiftExp : public Dihedral {
  public:
   DihedralCosineShiftExp(class LAMMPS *);
-  ~DihedralCosineShiftExp();
-  void compute(int, int);
+  virtual ~DihedralCosineShiftExp();
+  virtual void compute(int, int);
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
 
- private:
+ protected:
   bool *doExpansion;
   double *umin,*a,*opt1;
   double *sint;

src/USER-MISC/pair_cdeam.h

@@ -35,7 +35,7 @@ public:
   virtual ~PairCDEAM();
 
   /// Calculates the energies and forces for all atoms in the system.
-  void compute(int, int);
+  virtual void compute(int, int);
   
   /// Parses the pair_coeff command parameters for this pair style.
   void coeff(int, char **);
@@ -62,7 +62,7 @@ public:
       v = (v + hcoeff[i]) * x;
     }
     return v + hcoeff[0];
-  }
+  };
   
   // Calculates the derivative of the h(x) polynomial.
   inline double evalHprime(double x) const {
@@ -71,7 +71,7 @@ public:
       v = (v + (double)i * hcoeff[i]) * x;
     }
     return v + hcoeff[1];
-  }
+  };
   
   // We have two versions of the CD-EAM potential. The user can choose which one he wants to use:
   //
@@ -124,7 +124,7 @@ public:
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;
-  }
+  };
   
   // Converts a density value to an index value to be used in a spline table lookup.
   inline EAMTableIndex rhoToTableIndex(double rho) const {
@@ -135,43 +135,43 @@ public:
     index.p -= index.m;
     index.p = index.p <= 1.0 ? index.p : 1.0;
     return index;
-  }
+  };
   
   // Computes the derivative of rho(r)
   inline double RhoPrimeOfR(const EAMTableIndex& index, int itype, int jtype) const {
     const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m];
     return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
-  }
+  };
   
   // Computes rho(r)
   inline double RhoOfR(const EAMTableIndex& index, int itype, int jtype) const {
     const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m];
     return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
-  }
+  };
   
   // Computes the derivative of F(rho)
   inline double FPrimeOfRho(const EAMTableIndex& index, int itype) const {
     const double* coeff = frho_spline[type2frho[itype]][index.m];
     return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
-  }
+  };
   
   // Computes F(rho)
   inline double FofRho(const EAMTableIndex& index, int itype) const {
     const double* coeff = frho_spline[type2frho[itype]][index.m];
     return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
-  }
+  };
   
   // Computes the derivative of z2(r)
   inline double Z2PrimeOfR(const EAMTableIndex& index, int itype, int jtype) const {
     const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
     return (coeff[0]*index.p + coeff[1])*index.p + coeff[2];
-  }
-  
+  };
+
   // Computes z2(r)
   inline double Z2OfR(const EAMTableIndex& index, int itype, int jtype) const {
     const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
     return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
-  }
+  };
   
   // Computes pair potential V_ij(r).
   inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR) const {
@@ -180,7 +180,7 @@ public:
     const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m];
     const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6];
     return z2 * oneOverR;
-  }
+  };
   
   // Computes pair potential V_ij(r) and its derivative.
   inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR, double& phid) const {
@@ -194,7 +194,7 @@ public:
     const double phi = z2 * oneOverR;
     phid = z2p * oneOverR - phi * oneOverR;
     return phi;
-  }
+  };
   
   // Parameters
   
@@ -224,7 +224,7 @@ public:
      PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {}
    };
  
-};
+}
 
 #endif
 #endif

src/USER-MISC/pair_dipole_sf.h

@@ -26,9 +26,9 @@ namespace LAMMPS_NS {
 
 class PairDipoleSF : public Pair {
  public:
-    PairDipoleSF(class LAMMPS *);
-  ~PairDipoleSF();
-  void compute(int, int);
+  PairDipoleSF(class LAMMPS *);
+  virtual ~PairDipoleSF();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
@@ -38,7 +38,7 @@ class PairDipoleSF : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
 
- private:
+ protected:
   double cut_lj_global,cut_coul_global;
   double **cut_lj,**cut_ljsq;
   double **cut_coul,**cut_coulsq;

src/USER-MISC/pair_edip.cpp

@@ -554,15 +554,17 @@ void PairEDIP::allocateGrids(void)
 
 void PairEDIP::allocatePreLoops(void)
 {
-  memory->create(preInvR_ij,leadDimInteractionList,"edip:preInvR_ij");
-  memory->create(preExp3B_ij,leadDimInteractionList,"edip:preExp3B_ij");
-  memory->create(preExp3BDerived_ij,leadDimInteractionList,
+  int nthreads = comm->nthreads;
+  
+  memory->create(preInvR_ij,nthreads*leadDimInteractionList,"edip:preInvR_ij");
+  memory->create(preExp3B_ij,nthreads*leadDimInteractionList,"edip:preExp3B_ij");
+  memory->create(preExp3BDerived_ij,nthreads*leadDimInteractionList,
 		 "edip:preExp3BDerived_ij");
-  memory->create(preExp2B_ij,leadDimInteractionList,"edip:preExp2B_ij");
-  memory->create(preExp2BDerived_ij,leadDimInteractionList,
+  memory->create(preExp2B_ij,nthreads*leadDimInteractionList,"edip:preExp2B_ij");
+  memory->create(preExp2BDerived_ij,nthreads*leadDimInteractionList,
 		 "edip:preExp2BDerived_ij");
-  memory->create(prePow2B_ij,leadDimInteractionList,"edip:prePow2B_ij");
-  memory->create(preForceCoord,5*leadDimInteractionList,"edip:preForceCoord");
+  memory->create(prePow2B_ij,nthreads*leadDimInteractionList,"edip:prePow2B_ij");
+  memory->create(preForceCoord,5*nthreads*leadDimInteractionList,"edip:preForceCoord");
 }
 
 /* ----------------------------------------------------------------------

src/USER-MISC/pair_edip.h

@@ -27,14 +27,14 @@ namespace LAMMPS_NS {
 class PairEDIP : public Pair {
  public:
   PairEDIP(class LAMMPS *);
-  ~PairEDIP();
-  void compute(int, int);
+  virtual ~PairEDIP();
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
   void init_style();
 
- private:
+ protected:
   struct Param {
     double A, B;
     double cutoffA, cutoffC, cutsq;

src/USER-MISC/pair_lj_sf.cpp

@@ -92,12 +92,8 @@ void PairLJShiftedForce::compute(int eflag, int vflag)
     
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
-      
-      if (j < nall) factor_lj = 1.0;
-      else {
-	factor_lj = special_lj[j/nall];
-	j %= nall;
-      }
+      factor_lj = special_lj[sbmask(j)];
+      j &= NEIGHMASK;
       
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];

src/USER-SPH/atom_vec_meso.cpp

@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "atom_vec_meso.h"
 #include "atom.h"
+#include "comm.h"
 #include "domain.h"
 #include "modify.h"
 #include "fix.h"
@@ -68,12 +69,12 @@ void AtomVecMeso::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");
 
 	rho = memory->grow(atom->rho, nmax, "atom:rho");
-	drho = memory->grow(atom->drho, nmax, "atom:drho");
+	drho = memory->grow(atom->drho, nmax*comm->nthreads, "atom:drho");
 	e = memory->grow(atom->e, nmax, "atom:e");
-	de = memory->grow(atom->de, nmax, "atom:de");
+	de = memory->grow(atom->de, nmax*comm->nthreads, "atom:de");
 	vest = memory->grow(atom->vest, nmax, 3, "atom:vest");
 	cv = memory->grow(atom->cv, nmax, "atom:cv");
 
@@ -846,19 +847,19 @@ bigint AtomVecMeso::memory_usage() {
 	if (atom->memcheck("image"))
 		bytes += memory->usage(image, nmax);
 	if (atom->memcheck("x"))
-		bytes += memory->usage(x, nmax);
+		bytes += memory->usage(x, nmax, 3);
 	if (atom->memcheck("v"))
-		bytes += memory->usage(v, nmax);
+		bytes += memory->usage(v, nmax, 3);
 	if (atom->memcheck("f"))
-		bytes += memory->usage(f, nmax);
+		bytes += memory->usage(f, nmax*comm->nthreads, 3);
 	if (atom->memcheck("rho"))
 		bytes += memory->usage(rho, nmax);
 	if (atom->memcheck("drho"))
-		bytes += memory->usage(drho, nmax);
+		bytes += memory->usage(drho, nmax*comm->nthreads);
 	if (atom->memcheck("e"))
 		bytes += memory->usage(e, nmax);
 	if (atom->memcheck("de"))
-		bytes += memory->usage(de, nmax);
+		bytes += memory->usage(de, nmax*comm->nthreads);
 	if (atom->memcheck("cv"))
 		bytes += memory->usage(cv, nmax);
 	if (atom->memcheck("vest"))

src/atom_vec_ellipsoid.cpp

@@ -90,7 +90,7 @@ void AtomVecEllipsoid::grow(int n)
 
   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");
   ellipsoid = memory->grow(atom->ellipsoid,nmax,"atom:ellipsoid");
 
   if (atom->nextra_grow)
@@ -1250,7 +1250,7 @@ bigint AtomVecEllipsoid::memory_usage()
 
   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("ellipsoid")) bytes += memory->usage(ellipsoid,nmax);
   
   bytes += nmax_bonus*sizeof(Bonus);

src/atom_vec_sphere.cpp

@@ -102,7 +102,7 @@ void AtomVecSphere::grow(int n)
   radius = memory->grow(atom->radius,nmax,"atom:radius");
   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");
 
   if (atom->nextra_grow)
     for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@@ -1053,7 +1053,7 @@ bigint AtomVecSphere::memory_usage()
   if (atom->memcheck("radius")) bytes += memory->usage(radius,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);
 
   return bytes;
 }

src/bond_hybrid.cpp

@@ -206,8 +206,10 @@ void BondHybrid::settings(int narg, char **arg)
   // one exception is 1st arg of style "table", which is non-numeric
   // need a better way to skip these exceptions
 
+  int dummy;
   nstyles = 0;
   i = 0;
+
   while (i < narg) {
     for (m = 0; m < nstyles; m++)
       if (strcmp(arg[i],keywords[m]) == 0) 
@@ -216,7 +218,7 @@ void BondHybrid::settings(int narg, char **arg)
       error->all(FLERR,"Bond style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[i],"none") == 0) 
       error->all(FLERR,"Bond style hybrid cannot have none as an argument");
-    styles[nstyles] = force->new_bond(arg[i]);
+    styles[nstyles] = force->new_bond(arg[i],lmp->suffix,dummy);
     keywords[nstyles] = new char[strlen(arg[i])+1];
     strcpy(keywords[nstyles],arg[i]);
     istyle = i;
@@ -319,14 +321,14 @@ void BondHybrid::read_restart(FILE *fp)
 
   allocate();
 
-  int n;
+  int n,dummy;
   for (int m = 0; m < nstyles; m++) {
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
-    styles[m] = force->new_bond(keywords[m]);
+    styles[m] = force->new_bond(keywords[m],lmp->suffix,dummy);
   }
 }
 

src/compute_pair.cpp

@@ -37,6 +37,7 @@ ComputePair::ComputePair(LAMMPS *lmp, int narg, char **arg) :
   timeflag = 1;
 
   int n = strlen(arg[3]) + 1;
+  if (lmp->suffix) n += strlen(lmp->suffix) + 1;
   pstyle = new char[n];
   strcpy(pstyle,arg[3]);
 
@@ -46,8 +47,17 @@ ComputePair::ComputePair(LAMMPS *lmp, int narg, char **arg) :
     if (strcmp(arg[4],"ecoul") == 0) evalue = ECOUL;
   } else evalue = EPAIR;
 
+  // check if pair style with and without suffix exists
+
   pair = force->pair_match(pstyle,1);
-  if (!pair) error->all(FLERR,"Unrecognized pair style in compute pair command");
+  if (!pair && lmp->suffix) {
+    strcat(pstyle,"/");
+    strcat(pstyle,lmp->suffix);
+    pair = force->pair_match(pstyle,1);
+  }
+
+  if (!pair) 
+    error->all(FLERR,"Unrecognized pair style in compute pair command");
   npair = pair->nextra;
 
   if (npair) {
@@ -75,7 +85,8 @@ void ComputePair::init()
   // recheck for pair style in case it has been deleted
 
   pair = force->pair_match(pstyle,1);
-  if (!pair) error->all(FLERR,"Unrecognized pair style in compute pair command");
+  if (!pair) 
+    error->all(FLERR,"Unrecognized pair style in compute pair command");
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_box_relax.cpp

@@ -901,7 +901,7 @@ void FixBoxRelax::compute_press_target()
 
 /* ----------------------------------------------------------------------
    compute PV and strain energy for access to the user
-/* ---------------------------------------------------------------------- */
+   ---------------------------------------------------------------------- */
 
 double FixBoxRelax::compute_scalar()
 {

src/fix_gravity.h

@@ -32,11 +32,11 @@ class FixGravity : public Fix {
   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 compute_scalar();
 
- private:
+ protected:
   int style;
   double magnitude,dt;
   double phi,theta,phigrad,thetagrad;

src/fix_nve_sphere.h

@@ -27,13 +27,13 @@ namespace LAMMPS_NS {
 class FixNVESphere : public FixNVE {
  public:
   FixNVESphere(class LAMMPS *, int, char **);
-  ~FixNVESphere() {}
+  virtual ~FixNVESphere() {}
   int setmask();
   void init();
-  void initial_integrate(int);
-  void final_integrate();
+  virtual void initial_integrate(int);
+  virtual void final_integrate();
 
- private:
+ protected:
   int extra;
 };
 

src/fix_shear_history.h

@@ -35,7 +35,7 @@ class FixShearHistory : public Fix {
   int setmask();
   void init();
   void setup_pre_exchange();
-  void pre_exchange();
+  virtual void pre_exchange();
 
   double memory_usage();
   void grow_arrays(int);
@@ -48,7 +48,7 @@ class FixShearHistory : public Fix {
   int size_restart(int);
   int maxsize_restart();
 
- private:
+ protected:
   int *npartner;                // # of touching partners of each atom
   int **partner;                // tags for the partners
   double ***shearpartner;       // 3 shear values with the partner

src/force.cpp

@@ -112,7 +112,7 @@ void Force::init()
    create a pair style, called from input script or restart file
 ------------------------------------------------------------------------- */
 
-void Force::create_pair(const char *style, char *suffix)
+void Force::create_pair(const char *style, const char *suffix)
 {
   delete [] pair_style;
   if (pair) delete pair;
@@ -137,7 +137,7 @@ void Force::create_pair(const char *style, char *suffix)
    generate a pair class, first with suffix appended
 ------------------------------------------------------------------------- */
 
-Pair *Force::new_pair(const char *style, char *suffix, int &sflag)
+Pair *Force::new_pair(const char *style, const char *suffix, int &sflag)
 {
   if (suffix && lmp->suffix_enable) {
     sflag = 1;
@@ -185,20 +185,7 @@ Pair *Force::pair_match(const char *word, int exact)
   if (exact && strcmp(pair_style,word) == 0) return pair;
   else if (!exact && strstr(pair_style,word)) return pair;
 
-  else if (strcmp(pair_style,"hybrid") == 0) {
-    PairHybrid *hybrid = (PairHybrid *) pair;
-    count = 0;
-    for (int i = 0; i < hybrid->nstyles; i++) {
-      if (exact && strcmp(hybrid->keywords[i],word) == 0)
-	return hybrid->styles[i];
-      if (!exact && strstr(hybrid->keywords[i],word)) {
-	iwhich = i;
-	count++;
-      }
-    }
-    if (!exact && count == 1) return hybrid->styles[iwhich];
-
-  } else if (strcmp(pair_style,"hybrid/overlay") == 0) {
+  else if (strstr(pair_style,"hybrid/overlay")) {
     PairHybridOverlay *hybrid = (PairHybridOverlay *) pair;
     count = 0;
     for (int i = 0; i < hybrid->nstyles; i++) {
@@ -210,6 +197,19 @@ Pair *Force::pair_match(const char *word, int exact)
       }
     }
     if (!exact && count == 1) return hybrid->styles[iwhich];
+
+  } else if (strstr(pair_style,"hybrid")) {
+    PairHybrid *hybrid = (PairHybrid *) pair;
+    count = 0;
+    for (int i = 0; i < hybrid->nstyles; i++) {
+      if (exact && strcmp(hybrid->keywords[i],word) == 0)
+	return hybrid->styles[i];
+      if (!exact && strstr(hybrid->keywords[i],word)) {
+	iwhich = i;
+	count++;
+      }
+    }
+    if (!exact && count == 1) return hybrid->styles[iwhich];
   }
 
   return NULL;
@@ -219,23 +219,50 @@ Pair *Force::pair_match(const char *word, int exact)
    create a bond style, called from input script or restart file
 ------------------------------------------------------------------------- */
 
-void Force::create_bond(const char *style)
+void Force::create_bond(const char *style, const char *suffix)
 {
   delete [] bond_style;
   if (bond) delete bond;
 
-  bond = new_bond(style);
-  int n = strlen(style) + 1;
-  bond_style = new char[n];
-  strcpy(bond_style,style);
+  int sflag;
+  bond = new_bond(style,suffix,sflag);
+
+  if (sflag) {
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+    int n = strlen(estyle) + 1;
+    bond_style = new char[n];
+    strcpy(bond_style,estyle);
+  } else {
+    int n = strlen(style) + 1;
+    bond_style = new char[n];
+    strcpy(bond_style,style);
+  }
 }
 
 /* ----------------------------------------------------------------------
-   generate a bond class
+   generate a bond class, fist with suffix appended
 ------------------------------------------------------------------------- */
 
-Bond *Force::new_bond(const char *style)
+Bond *Force::new_bond(const char *style, const char *suffix, int &sflag)
 {
+  if (suffix && lmp->suffix_enable) {
+    sflag = 1;
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+
+    if (0) return NULL;
+
+#define BOND_CLASS
+#define BondStyle(key,Class) \
+    else if (strcmp(estyle,#key) == 0) return new Class(lmp);
+#include "style_bond.h"
+#undef BondStyle
+#undef BOND_CLASS
+  }
+
+  sflag = 0;
+
   if (strcmp(style,"none") == 0) return NULL;
 
 #define BOND_CLASS
@@ -267,23 +294,51 @@ Bond *Force::bond_match(const char *style)
    create an angle style, called from input script or restart file
 ------------------------------------------------------------------------- */
 
-void Force::create_angle(const char *style)
+void Force::create_angle(const char *style, const char *suffix)
 {
   delete [] angle_style;
   if (angle) delete angle;
 
-  angle = new_angle(style);
-  int n = strlen(style) + 1;
-  angle_style = new char[n];
-  strcpy(angle_style,style);
+  int sflag;
+  angle = new_angle(style,suffix,sflag);
+
+  if (sflag) {
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+    int n = strlen(estyle) + 1;
+    angle_style = new char[n];
+    strcpy(angle_style,estyle);
+  } else {
+    int n = strlen(style) + 1;
+    angle_style = new char[n];
+    strcpy(angle_style,style);
+  }
 }
 
 /* ----------------------------------------------------------------------
    generate an angle class
 ------------------------------------------------------------------------- */
 
-Angle *Force::new_angle(const char *style)
+Angle *Force::new_angle(const char *style, const char *suffix, int &sflag)
 {
+  if (suffix && lmp->suffix_enable) {
+    sflag = 1;
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+
+    if (0) return NULL;
+
+#define ANGLE_CLASS
+#define AngleStyle(key,Class) \
+    else if (strcmp(estyle,#key) == 0) return new Class(lmp);
+#include "style_angle.h"
+#undef AngleStyle
+#undef ANGLE_CLASS
+
+  }
+
+  sflag = 0;
+
   if (strcmp(style,"none") == 0) return NULL;
 
 #define ANGLE_CLASS
@@ -300,29 +355,58 @@ Angle *Force::new_angle(const char *style)
    create a dihedral style, called from input script or restart file
 ------------------------------------------------------------------------- */
 
-void Force::create_dihedral(const char *style)
+void Force::create_dihedral(const char *style, const char *suffix)
 {
   delete [] dihedral_style;
   if (dihedral) delete dihedral;
 
-  dihedral = new_dihedral(style);
-  int n = strlen(style) + 1;
-  dihedral_style = new char[n];
-  strcpy(dihedral_style,style);
+  int sflag;
+  dihedral = new_dihedral(style,suffix,sflag);
+
+  if (sflag) {
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+    int n = strlen(estyle) + 1;
+    dihedral_style = new char[n];
+    strcpy(dihedral_style,estyle);
+  } else {
+    int n = strlen(style) + 1;
+    dihedral_style = new char[n];
+    strcpy(dihedral_style,style);
+  }
 }
 
 /* ----------------------------------------------------------------------
    generate a dihedral class
 ------------------------------------------------------------------------- */
 
-Dihedral *Force::new_dihedral(const char *style)
+Dihedral *Force::new_dihedral(const char *style, const char *suffix, int &sflag)
 {
+  if (suffix && lmp->suffix_enable) {
+    sflag = 1;
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+
+    if (0) return NULL;
+
+#define DIHEDRAL_CLASS
+#define DihedralStyle(key,Class) \
+    else if (strcmp(estyle,#key) == 0) return new Class(lmp);
+#include "style_dihedral.h"
+#undef DihedralStyle
+#undef DIHEDRAL_CLASS
+
+  }
+
+  sflag = 0;
+
   if (strcmp(style,"none") == 0) return NULL;
 
 #define DIHEDRAL_CLASS
 #define DihedralStyle(key,Class) \
   else if (strcmp(style,#key) == 0) return new Class(lmp);
 #include "style_dihedral.h"
+#undef DihedralStyle
 #undef DIHEDRAL_CLASS
 
   else error->all(FLERR,"Invalid dihedral style");
@@ -333,23 +417,51 @@ Dihedral *Force::new_dihedral(const char *style)
    create an improper style, called from input script or restart file
 ------------------------------------------------------------------------- */
 
-void Force::create_improper(const char *style)
+void Force::create_improper(const char *style, const char *suffix)
 {
   delete [] improper_style;
   if (improper) delete improper;
 
-  improper = new_improper(style);
-  int n = strlen(style) + 1;
-  improper_style = new char[n];
-  strcpy(improper_style,style);
+  int sflag;
+  improper = new_improper(style,suffix,sflag);
+
+  if (sflag) {
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+    int n = strlen(estyle) + 1;
+    improper_style = new char[n];
+    strcpy(improper_style,estyle);
+  } else {
+    int n = strlen(style) + 1;
+    improper_style = new char[n];
+    strcpy(improper_style,style);
+  }
 }
 
 /* ----------------------------------------------------------------------
    generate a improper class
 ------------------------------------------------------------------------- */
 
-Improper *Force::new_improper(const char *style)
+Improper *Force::new_improper(const char *style, const char *suffix, int &sflag)
 {
+  if (suffix && lmp->suffix_enable) {
+    sflag = 1;
+    char estyle[256];
+    sprintf(estyle,"%s/%s",style,suffix);
+
+    if (0) return NULL;
+
+#define IMPROPER_CLASS
+#define ImproperStyle(key,Class) \
+    else if (strcmp(estyle,#key) == 0) return new Class(lmp);
+#include "style_improper.h"
+#undef ImproperStyle
+#undef IMPROPER_CLASS
+
+  }
+
+  sflag = 0;
+
   if (strcmp(style,"none") == 0) return NULL;
 
 #define IMPROPER_CLASS

src/force.h

@@ -69,22 +69,22 @@ class Force : protected Pointers {
   ~Force();
   void init();
 
-  void create_pair(const char *, char *suffix = NULL);
-  class Pair *new_pair(const char *, char *, int &);
+  void create_pair(const char *, const char *suffix = NULL);
+  class Pair *new_pair(const char *, const char *, int &);
   class Pair *pair_match(const char *, int);
 
-  void create_bond(const char *);
-  class Bond *new_bond(const char *);
+  void create_bond(const char *, const char *suffix = NULL);
+  class Bond *new_bond(const char *, const char *, int &);
   class Bond *bond_match(const char *); 
 
-  void create_angle(const char *);
-  class Angle *new_angle(const char *);
+  void create_angle(const char *, const char *suffix = NULL);
+  class Angle *new_angle(const char *, const char *, int &);
 
-  void create_dihedral(const char *);
-  class Dihedral *new_dihedral(const char *);
+  void create_dihedral(const char *, const char *suffix = NULL);
+  class Dihedral *new_dihedral(const char *, const char *, int &);
 
-  void create_improper(const char *);
-  class Improper *new_improper(const char *);
+  void create_improper(const char *, const char *suffix = NULL);
+  class Improper *new_improper(const char *, const char *, int &);
 
   void create_kspace(int, char **);
 

src/input.cpp

@@ -834,7 +834,7 @@ void Input::angle_style()
   if (narg < 1) error->all(FLERR,"Illegal angle_style command");
   if (atom->avec->angles_allow == 0) 
     error->all(FLERR,"Angle_style command when no angles allowed");
-  force->create_angle(arg[0]);
+  force->create_angle(arg[0],lmp->suffix);
   if (force->angle) force->angle->settings(narg-1,&arg[1]);
 }
 
@@ -875,7 +875,7 @@ void Input::bond_style()
   if (narg < 1) error->all(FLERR,"Illegal bond_style command");
   if (atom->avec->bonds_allow == 0) 
     error->all(FLERR,"Bond_style command when no bonds allowed");
-  force->create_bond(arg[0]);
+  force->create_bond(arg[0],lmp->suffix);
   if (force->bond) force->bond->settings(narg-1,&arg[1]);
 }
 
@@ -937,7 +937,7 @@ void Input::dihedral_style()
   if (narg < 1) error->all(FLERR,"Illegal dihedral_style command");
   if (atom->avec->dihedrals_allow == 0) 
     error->all(FLERR,"Dihedral_style command when no dihedrals allowed");
-  force->create_dihedral(arg[0]);
+  force->create_dihedral(arg[0],lmp->suffix);
   if (force->dihedral) force->dihedral->settings(narg-1,&arg[1]);
 }
 
@@ -1014,7 +1014,7 @@ void Input::improper_style()
   if (narg < 1) error->all(FLERR,"Illegal improper_style command");
   if (atom->avec->impropers_allow == 0) 
     error->all(FLERR,"Improper_style command when no impropers allowed");
-  force->create_improper(arg[0]);
+  force->create_improper(arg[0],lmp->suffix);
   if (force->improper) force->improper->settings(narg-1,&arg[1]);
 }
 

src/pair.h

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

src/pair_born.h

@@ -27,8 +27,9 @@ namespace LAMMPS_NS {
 class PairBorn : public Pair {
  public:
   PairBorn(class LAMMPS *);
-  ~PairBorn();
-  void compute(int, int);
+  virtual ~PairBorn();
+
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -39,7 +40,7 @@ class PairBorn : public Pair {
   double single(int, int, int, int, double, double, double, double &);
   void *extract(char *, int &);
 
- private:
+ protected:
   double cut_global;
   double **cut;
   double **a,**rho,**sigma,**c, **d;

src/pair_buck.h

@@ -28,7 +28,7 @@ class PairBuck : public Pair {
  public:
   PairBuck(class LAMMPS *);
   virtual ~PairBuck();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);

src/pair_coul_debye.h

@@ -27,13 +27,13 @@ namespace LAMMPS_NS {
 class PairCoulDebye : public PairCoulCut {
  public:
   PairCoulDebye(class LAMMPS *);
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
 
- private:
+ protected:
   double kappa;
 };
 

src/pair_dpd_tstat.h

@@ -34,7 +34,7 @@ class PairDPDTstat : public PairDPD {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
 
- private:
+ protected:
   double t_start,t_stop;
 };
 

src/pair_gauss.h

@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
 class PairGauss : public Pair {
  public:
   PairGauss(class LAMMPS *);
-  ~PairGauss();
-  void compute(int,int);
+  virtual ~PairGauss();
+  virtual void compute(int,int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);  
@@ -39,7 +39,7 @@ class PairGauss : public Pair {
   double single(int, int, int, int, double, double, double, double &);
   void *extract(char *, int &);
  
- private: 
+ protected: 
   double cut_global;
   double **cut;
   double **a,**b;

src/pair_hybrid.cpp

@@ -227,9 +227,9 @@ void PairHybrid::settings(int narg, char **arg)
   // need a better way to skip these exceptions
 
   int dummy;
-
   nstyles = 0;
   i = 0;
+
   while (i < narg) {
     for (m = 0; m < nstyles; m++)
       if (strcmp(arg[i],keywords[m]) == 0) 
@@ -238,7 +238,7 @@ void PairHybrid::settings(int narg, char **arg)
       error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[i],"none") == 0) 
       error->all(FLERR,"Pair style hybrid cannot have none as an argument");
-    styles[nstyles] = force->new_pair(arg[i],NULL,dummy);
+    styles[nstyles] = force->new_pair(arg[i],lmp->suffix,dummy);
     keywords[nstyles] = new char[strlen(arg[i])+1];
     strcpy(keywords[nstyles],arg[i]);
     istyle = i;
@@ -582,7 +582,7 @@ void PairHybrid::read_restart(FILE *fp)
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
-    styles[m] = force->new_pair(keywords[m],NULL,dummy);
+    styles[m] = force->new_pair(keywords[m],lmp->suffix,dummy);
     styles[m]->read_restart_settings(fp);
   }
 }

src/pair_lj96_cut.h

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

src/pair_lj_cubic.h

@@ -27,8 +27,9 @@ namespace LAMMPS_NS {
 class PairLJCubic : public Pair {
  public:
   PairLJCubic(class LAMMPS *);
-  ~PairLJCubic();
-  void compute(int, int);
+  virtual ~PairLJCubic();
+
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -38,7 +39,7 @@ class PairLJCubic : public Pair {
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
 
- private:
+ protected:
   double **cut,**cut_inner,**cut_inner_sq;
   double **epsilon,**sigma;
   double **lj1,**lj2,**lj3,**lj4;

src/pair_lj_cut_coul_debye.h

@@ -28,7 +28,7 @@ class PairLJCutCoulDebye : public PairLJCutCoulCut {
  public:
   PairLJCutCoulDebye(class LAMMPS *);
   virtual ~PairLJCutCoulDebye() {}
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);

src/pair_lj_expand.h

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

src/pair_lj_smooth.h

@@ -28,7 +28,7 @@ class PairLJSmooth : public Pair {
  public:
   PairLJSmooth(class LAMMPS *);
   virtual ~PairLJSmooth();
-  void compute(int, int);
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);

src/pair_morse.h

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

src/pair_soft.h

@@ -29,8 +29,9 @@ class PairSoft : public Pair {
 
  public:
   PairSoft(class LAMMPS *);
-  ~PairSoft();
-  void compute(int, int);
+  virtual ~PairSoft();
+
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -41,7 +42,7 @@ class PairSoft : public Pair {
   double single(int, int, int, int, double, double, double, double &);
   void *extract(char *, int &);
 
- private:
+ protected:
   double PI;
   double cut_global;
   double **prefactor;

src/pair_table.h

@@ -27,8 +27,9 @@ namespace LAMMPS_NS {
 class PairTable : public Pair {
  public:
   PairTable(class LAMMPS *);
-  ~PairTable();
-  void compute(int, int);
+  virtual ~PairTable();
+
+  virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   double init_one(int, int);
@@ -39,8 +40,9 @@ class PairTable : public Pair {
   double single(int, int, int, int, double, double, double, double &);
   void *extract(char *, int &);
 
- private:
+ protected:
   int tabstyle,tablength;
+  enum {LOOKUP=0, LINEAR=1, SPLINE=2, BITMAP=3};
   struct Table {
     int ninput,rflag,fpflag,match,ntablebits;
     int nshiftbits,nmask;

src/set.cpp

@@ -30,7 +30,10 @@
 #include "math_extra.h"
 #include "error.h"
 
+#include "math_const.h"
+
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
 enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
@@ -42,13 +45,6 @@ enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
 
 /* ---------------------------------------------------------------------- */
 
-Set::Set(LAMMPS *lmp) : Pointers(lmp)
-{
-  PI = 4.0*atan(1.0);
-}
-
-/* ---------------------------------------------------------------------- */
-
 void Set::command(int narg, char **arg)
 {
   if (domain->box_exist == 0)
@@ -415,11 +411,11 @@ void Set::set(int keyword)
 
     else if (keyword == DENSITY) {
       if (atom->radius_flag && atom->radius[i] > 0.0)
-	atom->rmass[i] = 4.0*PI/3.0 * 
+	atom->rmass[i] = 4.0*MY_PI*THIRD * 
 	  atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
       else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
 	double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
-	atom->rmass[i] = 4.0*PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
+	atom->rmass[i] = 4.0*MY_PI*THIRD * shape[0]*shape[1]*shape[2] * dvalue;
       } else atom->rmass[i] = dvalue;
 
     // reset any or all of 3 image flags
@@ -450,7 +446,7 @@ void Set::set(int keyword)
       if (atom->ellipsoid[i] < 0)
 	error->one(FLERR,"Cannot set quaternion for atom that is not an ellipsoid");
       double *quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
-      double theta2 = 0.5 * PI * wvalue/180.0;
+      double theta2 = MY_PI2 * wvalue/180.0;
       double sintheta2 = sin(theta2);
       quat[0] = cos(theta2);
       quat[1] = xvalue * sintheta2;
@@ -554,8 +550,8 @@ void Set::setrandom(int keyword)
 	  s = random->uniform();
 	  t1 = sqrt(1.0-s);
 	  t2 = sqrt(s);
-	  theta1 = 2.0*PI*random->uniform();
-	  theta2 = 2.0*PI*random->uniform();
+	  theta1 = 2.0*MY_PI*random->uniform();
+	  theta2 = 2.0*MY_PI*random->uniform();
 	  quat[0] = cos(theta2)*t2;
 	  quat[1] = sin(theta1)*t1;
 	  quat[2] = cos(theta1)*t1;
@@ -572,7 +568,7 @@ void Set::setrandom(int keyword)
 		       "that is not an ellipsoid");
 	  quat = bonus[ellipsoid[i]].quat;
 	  random->reset(seed,x[i]);
-	  theta2 = PI*random->uniform();
+	  theta2 = MY_PI*random->uniform();
 	  quat[0] = cos(theta2);
 	  quat[1] = 0.0;
 	  quat[2] = 0.0;

src/set.h

@@ -26,7 +26,7 @@ namespace LAMMPS_NS {
 
 class Set : protected Pointers {
  public:
-  Set(class LAMMPS *);
+  Set(class LAMMPS *lmp) : Pointers(lmp) {};
   void command(int, char **);
 
  private:

src/variable.cpp

@@ -32,7 +32,10 @@
 #include "memory.h"
 #include "error.h"
 
+#include "math_const.h"
+
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define VARDELTA 4
 #define MAXLEVEL 4
@@ -88,8 +91,6 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
   precedence[MULTIPLY] = precedence[DIVIDE] = 6;
   precedence[CARAT] = 7;
   precedence[UNARY] = precedence[NOT] = 8;
-
-  PI = 4.0*atan(1.0);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1767,7 +1768,7 @@ double Variable::collapse_tree(Tree *tree)
     tree->type = VALUE;
     if (arg3 == 0.0) error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
-    double omega = 2.0*PI/arg3;
+    double omega = 2.0*MY_PI/arg3;
     tree->value = arg1 + arg2*sin(omega*delta*update->dt);
     return tree->value;
   }
@@ -1781,7 +1782,7 @@ double Variable::collapse_tree(Tree *tree)
     tree->type = VALUE;
     if (arg3 == 0.0) error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
-    double omega = 2.0*PI/arg3;
+    double omega = 2.0*MY_PI/arg3;
     tree->value = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
     return tree->value;
   }
@@ -1997,7 +1998,7 @@ double Variable::eval_tree(Tree *tree, int i)
     arg3 = eval_tree(tree->right,i);
     if (arg3 == 0.0) error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
-    double omega = 2.0*PI/arg3;
+    double omega = 2.0*MY_PI/arg3;
     arg = arg1 + arg2*sin(omega*delta*update->dt);
     return arg;
   }
@@ -2008,7 +2009,7 @@ double Variable::eval_tree(Tree *tree, int i)
     arg3 = eval_tree(tree->right,i);
     if (arg3 == 0.0) error->one(FLERR,"Invalid math function in variable formula");
     double delta = update->ntimestep - update->beginstep;
-    double omega = 2.0*PI/arg3;
+    double omega = 2.0*MY_PI/arg3;
     arg = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
     return arg;
   }
@@ -2385,7 +2386,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       if (value3 == 0.0)
 	error->all(FLERR,"Invalid math function in variable formula");
       double delta = update->ntimestep - update->beginstep;
-      double omega = 2.0*PI/value3;
+      double omega = 2.0*MY_PI/value3;
       double value = value1 + value2*sin(omega*delta*update->dt);
       argstack[nargstack++] = value;
     }
@@ -2399,7 +2400,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
       if (value3 == 0.0)
 	error->all(FLERR,"Invalid math function in variable formula");
       double delta = update->ntimestep - update->beginstep;
-      double omega = 2.0*PI/value3;
+      double omega = 2.0*MY_PI/value3;
       double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
       argstack[nargstack++] = value;
     }
@@ -3058,7 +3059,7 @@ int Variable::is_constant(char *word)
 
 double Variable::constant(char *word)
 {
-  if (strcmp(word,"PI") == 0) return PI;
+  if (strcmp(word,"PI") == 0) return MY_PI;
   return 0.0;
 }
 

src/variable.h

@@ -35,7 +35,6 @@ class Variable : protected Pointers {
   double evaluate_boolean(char *);
 
  private:
-  int me;
   int nvar;                // # of defined variables
   int maxvar;              // max # of variables arrays can hold
   char **names;            // name of each variable
@@ -44,13 +43,13 @@ class Variable : protected Pointers {
   int *which;              // next available value for each variable
   int *pad;                // 1 = pad loop/uloop variables with 0s, 0 = no pad
   char ***data;            // str value of each variable's values
-  double PI;
 
   class RanMars *randomequal;   // random number generator for equal-style vars
   class RanMars *randomatom;    // random number generator for atom-style vars
 
   int precedence[16];      // precedence level of math operators
                            // set length to include OR in enum
+  int me;
 
   struct Tree {            // parse tree for atom-style variables
     double value;