updated and slightly refactored tabulation for vashishta pair style

- tables are now dimensioned by nelements instead of ntypes
- tables are only created if used
- correctly identify max size of table
- add test for illegal cutoff for tabulation
- allocated memory for tables is accounted for
- add example input using 16-bit tables

examples/vashishta/in.indiumphosphide-tabulated

@@ -0,0 +1,76 @@
+# calculate the energy volume curve for InP zincblende
+
+# define volume range and filename
+
+variable	ndelta equal 100
+variable	volatom_min equal 20.0
+variable	volatom_max equal 29.0
+variable	evsvolfile string evsvol.dat
+
+# set up cell 
+
+units		metal
+ 
+boundary	p p p
+
+# setup loop variables for box volume
+
+variable	amin equal ${volatom_min}^(1/3)*2
+variable 	delta equal (${volatom_max}-${volatom_min})/${ndelta} 
+variable	scale equal (${delta}/v_volatom+1)^(1/3)
+
+# set up 8 atom InP zincblende unit cell
+
+lattice diamond ${amin}
+
+region		box prism &
+	 	0 1 &
+		0 1 &
+		0 1 &
+		0 0 0
+
+create_box	2 box
+
+create_atoms	1	box       &
+			basis 5 2 &
+			basis 6 2 &
+			basis 7 2 &
+			basis 8 2
+
+mass 		1 114.76
+mass 		2 30.98
+
+# choose potential
+
+pair_style	vashishta
+pair_coeff 	* * InP.vashishta In P
+pair_modify     table 16
+
+# setup neighbor style
+
+neighbor 	1.0 nsq
+neigh_modify 	once no every 1 delay 0 check yes
+
+# setup output
+
+thermo_style 	custom step temp pe press vol
+thermo_modify 	norm no
+variable 	volatom equal vol/atoms
+variable 	eatom equal pe/atoms
+print 		"# Volume [A^3/atom] Energy [eV/atom]" file ${evsvolfile}
+
+# loop over range of volumes
+
+label 		loop
+variable 	i loop ${ndelta}
+
+change_box 	all x scale ${scale} y scale ${scale} z scale ${scale} remap
+
+# calculate energy
+# no energy minimization needed for zincblende
+ 
+run 		0
+print 		"${volatom} ${eatom}" append ${evsvolfile}
+
+next 		i
+jump 		SELF loop

examples/vashishta/in.sio2-tabulated

@@ -0,0 +1,29 @@
+# test Vashishta potential for quartz
+
+units		metal
+boundary	p p p
+
+atom_style	atomic
+
+read_data	data.quartz
+
+replicate       4 4 4
+velocity	all create 2000.0 277387 mom yes
+displace_atoms	all move 0.05 0.9 0.4 units box
+
+pair_style 	vashishta
+pair_coeff	* *  SiO.1990.vashishta Si O
+pair_modify	table 16
+
+neighbor	0.3 bin
+neigh_modify	delay 10
+
+fix		1 all nve
+thermo		10
+timestep	0.001
+
+#dump		1 all cfg 10 *.cfg mass type xs ys zs vx vy vz fx fy fz
+#dump_modify	1 element Si O
+
+run		100
+

src/MANYBODY/pair_vashishta.cpp

@@ -97,11 +97,10 @@ void PairVashishta::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
 
       nTablebits = force->inumeric(FLERR,arg[iarg+1]);
-      if (nTablebits > sizeof(float)*CHAR_BIT) {
-        error->all(FLERR,"Too many total bits for bitmapped lookup table");
-      }
+      if ((nTablebits < 0) || (nTablebits > sizeof(int)*CHAR_BIT-1))
+        error->all(FLERR,"Illegal interpolation table size");
 
-      if(nTablebits == 0) {
+      if (nTablebits == 0) {
         useTable = false;
       } else {
         useTable = true;
@@ -112,11 +111,12 @@ void PairVashishta::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
       
       tabinner = force->numeric(FLERR,arg[iarg+1]);
+      if (tabinner <= 0.0)
+        error->all(FLERR,"Illegal inner cutoff for tabulation");
       iarg += 2;
     }
   }
-
-  createTable();
+  updateTables();
 }
 
 void PairVashishta::validateNeigh3Body() {
@@ -129,40 +129,34 @@ void PairVashishta::validateNeigh3Body() {
     }
 }
 
-void PairVashishta::createTable()
+void PairVashishta::updateTables()
 {
-  int ntypes = atom->ntypes+1;
-  tabinnersq = tabinner*tabinner;
+  memory->destroy(forceTable);
+  memory->destroy(potentialTable);
+  forceTable = NULL;
+  potentialTable = NULL;
 
-  int ntable = 1;
-  for (int i = 0; i < nTablebits; i++) ntable *= 2;
+  if (!useTable) return;
+
+  int ntable = 1<<nTablebits;
+  tabinnersq = tabinner*tabinner;
 
   deltaR2 = (cutmax*cutmax - tabinnersq) / (ntable-1);
   oneOverDeltaR2 = 1.0/deltaR2;
 
-  memory->destroy(forceTable);
-  memory->destroy(potentialTable);
-  memory->create(forceTable,ntypes,ntypes,ntable+1,"pair:vashishta:forceTable");
-  memory->create(potentialTable,ntypes,ntypes,ntable+1,"pair:vashishta:potentialTable");
+  memory->create(forceTable,nelements,nelements,ntable+1,"pair:vashishta:forceTable");
+  memory->create(potentialTable,nelements,nelements,ntable+1,"pair:vashishta:potentialTable");
 
-  for (int ii = 0; ii < ntypes; ii++) {
-    int i = map[ii];
-    for (int jj = 0; jj < ntypes; jj++) {
-      int j = map[jj];
-      if (i < 0 || j < 0 || ii == 0 || jj == 0) {
-        for(int tableIndex=0; tableIndex<=ntable; tableIndex++) {
-            forceTable[ii][jj][tableIndex] = 0;
-            potentialTable[ii][jj][tableIndex] = 0;
-        }
-      } else {
+  for (int i = 0; i < nelements; i++) {
+    for (int j = 0; j < nelements; j++) {
       int ijparam = elem2param[i][j][j];
-        for(int tableIndex=0; tableIndex<=ntable; tableIndex++) {
-            double rsq = tabinnersq + tableIndex*deltaR2;
+      for(int idx=0; idx <= ntable; idx++) {
+        double rsq = tabinnersq + idx*deltaR2;
         double fpair, eng;
-            twobody(&params[ijparam], rsq, fpair, 1, eng, false /* Don't ask for tabulated since we are generating it now*/ );
-            forceTable[ii][jj][tableIndex] = fpair;
-            potentialTable[ii][jj][tableIndex] = eng;
-        }
+        // call analytical version of two-body term function
+        twobody(&params[ijparam], rsq, fpair, 1, eng, false);
+        forceTable[i][j][idx] = fpair;
+        potentialTable[i][j][idx] = eng;
       }
     }
   }
@@ -642,7 +636,7 @@ void PairVashishta::setup_params()
     if (params[m].r0 > cutmax) cutmax = params[m].r0;
   }
 
-  createTable();
+  updateTables();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -660,12 +654,12 @@ void PairVashishta::twobody(Param *param, double rsq, double &fforce,
     int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2;
     double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; // double - int will only keep the 0.xxxx part
 
-    double force0 = forceTable[param->ielement+1][param->jelement+1][tableIndex];
-    double force1 = forceTable[param->ielement+1][param->jelement+1][tableIndex+1];
+    double force0 = forceTable[param->ielement][param->jelement][tableIndex];
+    double force1 = forceTable[param->ielement][param->jelement][tableIndex+1];
     fforce = (1.0 - fraction)*force0 + fraction*force1; // force is linearly interpolated between the two values
     if(evflag) {
-        double energy0 = potentialTable[param->ielement+1][param->jelement+1][tableIndex];
-        double energy1 = potentialTable[param->ielement+1][param->jelement+1][tableIndex+1];
+        double energy0 = potentialTable[param->ielement][param->jelement][tableIndex];
+        double energy1 = potentialTable[param->ielement][param->jelement][tableIndex+1];
         eng = (1.0 - fraction)*energy0 + fraction*energy1;
     }
   } else {
@@ -747,3 +741,15 @@ void PairVashishta::threebody(Param *paramij, Param *paramik, Param *paramijk,
 
   if (eflag) eng = facrad;
 }
+
+/* ----------------------------------------------------------------------
+ *  add memory usage for tabulation
+ * ---------------------------------------------------------------------- */
+
+double PairVashishta::memory_usage()
+{
+  double bytes = Pair::memory_usage();
+  if (useTable)
+    bytes += 2*nelements*nelements*sizeof(double)*(1<<nTablebits);
+  return bytes;
+}

src/MANYBODY/pair_vashishta.h

@@ -34,6 +34,7 @@ class PairVashishta : public Pair {
   void coeff(int, char **);
   virtual double init_one(int, int);
   virtual void init_style();
+  virtual double memory_usage();
 
  protected:
   struct Param {
@@ -71,7 +72,7 @@ class PairVashishta : public Pair {
   virtual void allocate();
   void read_file(char *);
   void setup_params();
-  void createTable();
+  void updateTables();
   void validateNeigh3Body();
   void twobody(Param *, double, double &, int, double &, bool tabulated = false);
   void threebody(Param *, Param *, Param *, double, double, double *, double *,

src/USER-OMP/pair_vashishta_omp.cpp

@@ -144,7 +144,7 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
         if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
       }
 
-      twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl, useTable);
+      twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl,useTable);
 
       fxtmp += delx*fpair;
       fytmp += dely*fpair;