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

src/MANYBODY/Install.csh

@@ -4,12 +4,14 @@ if ($1 == 1) then
 
   cp style_manybody.h ..
 
+  cp pair_airebo.cpp ..
   cp pair_eam.cpp ..
   cp pair_eam_alloy.cpp ..
   cp pair_eam_fs.cpp ..
   cp pair_sw.cpp ..
   cp pair_tersoff.cpp ..
 
+  cp pair_airebo.h ..
   cp pair_eam.h ..
   cp pair_eam_alloy.h ..
   cp pair_eam_fs.h ..
@@ -21,12 +23,14 @@ else if ($1 == 0) then
   rm ../style_manybody.h
   touch ../style_manybody.h
 
+  rm ../pair_airebo.cpp
   rm ../pair_eam.cpp
   rm ../pair_eam_alloy.cpp
   rm ../pair_eam_fs.cpp
   rm ../pair_sw.cpp
   rm ../pair_tersoff.cpp
 
+  rm ../pair_airebo.h
   rm ../pair_eam.h
   rm ../pair_eam_alloy.h
   rm ../pair_eam_fs.h

src/MANYBODY/pair_airebo.h

@@ -0,0 +1,117 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef PAIR_AIREBO_H
+#define PAIR_AIREBO_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairAIREBO : public Pair {
+ public:
+  PairAIREBO(class LAMMPS *);
+  virtual ~PairAIREBO();
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void write_restart(FILE *) {}
+  void read_restart(FILE *) {}
+  void write_restart_settings(FILE *) {}
+  void read_restart_settings(FILE *) {}
+  void single(int, int, int, int, double, double, double, int, One &) {}
+
+ private:
+  int me;
+  int ljflag,torflag;              // 0/1 if LJ,torsion terms included
+  int maxlocal;                    // size of numneigh, firstneigh arrays
+  int **pages;                     // neighbor list pages
+  int maxpage;                     // # of pages currently allocated
+  int pgsize;                      // size of neighbor page
+  int oneatom;                     // max # of neighbors for one atom
+  int npage;                       // current page in page list
+  int *map;
+  double cutlj;                    // user-specified LJ cutoff
+  double cutljrebosq;              // cut for when to compute
+                                   // REBO neighs of ghost atoms
+
+  double **cutljsq;                // LJ cutoffs for C,H types
+  double **lj1,**lj2,**lj3,**lj4;  // pre-computed LJ coeffs for C,H types
+  double cut3rebo;                 // maximum distance for 3rd REBO neigh
+
+  int *REBO_numneigh;              // # of pair neighbors for each atom
+  int **REBO_firstneigh;           // ptr to 1st neighbor of each atom
+  double *nC,*nH;                  // sum of weighing fns with REBO neighs
+
+  double smin,Nmin,Nmax,NCmin,NCmax,thmin,thmax;
+  double rcmin[2][2],rcmax[2][2],rcmaxsq[2][2],rcmaxp[2][2];
+  double Q[2][2],alpha[2][2],A[2][2],rho[2][2],BIJc[2][2][3],Beta[2][2][3];
+  double rcLJmin[2][2],rcLJmax[2][2],rcLJmaxsq[2][2],bLJmin[2][2],bLJmax[2][2];
+  double epsilon[2][2],sigma[2][2],epsilonT[2][2];
+
+  // spline coefficients
+
+  double gCdom[5],gC1[4][6],gC2[4][6],gHdom[4],gH[3][6];
+  double pCCdom[2][2],pCHdom[2][2],pCC[4][4][16],pCH[4][4][16];
+  double piCCdom[3][2],piCHdom[3][2],piHHdom[3][2];
+  double piCC[4][4][9][64],piCH[4][4][9][64],piHH[4][4][9][64];
+  double Tijdom[3][2],Tijc[4][4][9][64];
+
+  // spline knot values
+
+  double PCCf[5][5],PCCdfdx[5][5],PCCdfdy[5][5],PCHf[5][5];
+  double PCHdfdx[5][5],PCHdfdy[5][5];
+  double piCCf[5][5][10],piCCdfdx[5][5][10];
+  double piCCdfdy[5][5][10],piCCdfdz[5][5][10];
+  double piCHf[5][5][10],piCHdfdx[5][5][10];
+  double piCHdfdy[5][5][10],piCHdfdz[5][5][10];
+  double piHHf[5][5][10],piHHdfdx[5][5][10];
+  double piHHdfdy[5][5][10],piHHdfdz[5][5][10];
+  double Tf[5][5][10],Tdfdx[5][5][10],Tdfdy[5][5][10],Tdfdz[5][5][10];
+
+  void REBO_neigh();
+  void FREBO(int, double **);
+  void FLJ(int, double **);
+  void TORSION(int, double **);
+
+  double bondorder(int, int, double *, double, double, double **);
+  double bondorderLJ(int, int, double *, double, double,
+		     double *, double, double **);
+
+  double Sp(double, double, double, double &);
+  double Sp2(double, double, double, double &);
+
+  double gSpline(double, double, int, double *, double *);
+  double PijSpline(double, double, int, int, double *);
+  double piRCSpline(double, double, double, int, int, double *);
+  double TijSpline(double, double, double, double *);
+
+  double kronecker(int, int);
+
+  void add_pages(int);
+  void read_file(char *);
+
+  double Sp5th(double, double *, double *);
+  double Spbicubic(double, double, double *, double *);
+  double Sptricubic(double, double, double, double *, double *);
+  void spline_init();
+
+  void allocate();
+};
+
+}
+
+#endif
+

src/MANYBODY/style_manybody.h

@@ -12,6 +12,7 @@
 ------------------------------------------------------------------------- */
 
 #ifdef PairInclude
+#include "pair_airebo.h"
 #include "pair_eam.h"
 #include "pair_eam_alloy.h"
 #include "pair_eam_fs.h"
@@ -20,6 +21,7 @@
 #endif
 
 #ifdef PairClass
+PairStyle(airebo,PairAIREBO)
 PairStyle(eam,PairEAM)
 PairStyle(eam/alloy,PairEAMAlloy)
 PairStyle(eam/fs,PairEAMFS)

src/neighbor.cpp

@@ -44,6 +44,7 @@ using namespace LAMMPS_NS;
 #define SMALL 1.0e-6
 #define EXDELTA 1
 #define BIG 1.0e20
+#define CUT2BIN_RATIO 100
 
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
 #define MAX(a,b) ((a) > (b) ? (a) : (b))
@@ -239,6 +240,7 @@ void Neighbor::init()
   int i,j,m,n;
 
   ncalls = ndanger = 0;
+  dimension = domain->dimension;
   triclinic = domain->triclinic;
 
   // error check
@@ -520,19 +522,19 @@ void Neighbor::init()
       } else if (style == BIN) {
 	if (force->newton_pair == 0) {
 	  half_build = &Neighbor::granular_bin_no_newton;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_no_newton;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_no_newton;
 	} else if (triclinic) {
 	  half_build = &Neighbor::granular_bin_newton_tri;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_newton_tri;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_newton_tri;
 	} else {
 	  half_build = &Neighbor::granular_bin_newton;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_newton;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_newton;
@@ -547,19 +549,19 @@ void Neighbor::init()
       } else if (style == BIN) {
 	if (force->newton_pair == 0) {
 	  half_build = &Neighbor::respa_bin_no_newton;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_no_newton;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_no_newton;
 	} else if (triclinic) {
 	  half_build = &Neighbor::respa_bin_newton_tri;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_newton_tri;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_newton_tri;
 	} else {
 	  half_build = &Neighbor::respa_bin_newton;
-	  if (domain->dimension == 3)
+	  if (dimension == 3)
 	    half_stencil = &Neighbor::stencil_half_3d_newton;
 	  else
 	    half_stencil = &Neighbor::stencil_half_2d_newton;
@@ -582,7 +584,7 @@ void Neighbor::init()
 	    half_stencil = &Neighbor::stencil_none;
 	  } else {
 	    half_build = &Neighbor::half_bin_no_newton;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_no_newton;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_no_newton;
@@ -593,13 +595,13 @@ void Neighbor::init()
 	    half_stencil = &Neighbor::stencil_none;
 	  } else if (triclinic) {
 	    half_build = &Neighbor::half_bin_newton_tri;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_newton_tri;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_newton_tri;
 	  } else {
 	    half_build = &Neighbor::half_bin_newton;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_newton;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_newton;
@@ -612,7 +614,7 @@ void Neighbor::init()
 	    half_stencil = &Neighbor::stencil_none;
 	  } else {
 	    half_build = &Neighbor::half_bin_no_newton_multi;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_no_newton_multi;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_no_newton_multi;
@@ -623,13 +625,13 @@ void Neighbor::init()
 	    half_stencil = &Neighbor::stencil_none;
 	  } else if (triclinic) {
 	    half_build = &Neighbor::half_bin_newton_multi_tri;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_newton_multi_tri;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_newton_multi_tri;
 	  } else {
 	    half_build = &Neighbor::half_bin_newton_multi;
-	    if (domain->dimension == 3)
+	    if (dimension == 3)
 	      half_stencil = &Neighbor::stencil_half_3d_newton_multi;
 	    else
 	      half_stencil = &Neighbor::stencil_half_2d_newton_multi;
@@ -644,13 +646,13 @@ void Neighbor::init()
     if (style == NSQ) full_build = &Neighbor::full_nsq;
     else if (style == BIN) {
       full_build = &Neighbor::full_bin;
-      if (domain->dimension == 3)
+      if (dimension == 3)
 	full_stencil = &Neighbor::stencil_full_3d;
       else
 	full_stencil = &Neighbor::stencil_full_2d;
     } else {
       full_build = &Neighbor::full_bin_multi;
-      if (domain->dimension == 3)
+      if (dimension == 3)
 	full_stencil = &Neighbor::stencil_full_3d_multi;
       else
 	full_stencil = &Neighbor::stencil_full_2d_multi;
@@ -922,16 +924,19 @@ void Neighbor::build_full()
 
 /* ----------------------------------------------------------------------
    setup neighbor binning parameters
-   bin numbering is global: 0 = 0.0 to binsize, 1 = binsize to 2*binsize
+   bin numbering in each dimension is global:
-			    nbin-1 = bbox-binsize to binsize
+     0 = 0.0 to binsize, 1 = binsize to 2*binsize, etc
-			    nbin = bbox to bbox+binsize
+     nbin-1,nbin,etc = bbox-binsize to binsize, bbox to bbox+binsize, etc
-                            -1,-2,etc = -binsize to 0.0, -2*size to -size, etc
+     -1,-2,etc = -binsize to 0.0, -2*size to -size, etc
    code will work for any binsize
      since next(xyz) and stencil extend as far as necessary
      binsize = 1/2 of cutoff is roughly optimal
-   for orthogonal boxes, prd must be filled exactly by integer # of bins
+   for orthogonal boxes:
-     so procs on both sides of PBC see same bin boundary
+     a dim must be filled exactly by integer # of bins
-   for triclinic, tilted simulation box cannot contain integer # of bins
+     in periodic, procs on both sides of PBC must see same bin boundary
+     in non-periodic, coord2bin() still assumes this by use of nbin xyz
+   for triclinic boxes:
+     tilted simulation box cannot contain integer # of bins
      stencil & neigh list built differently to account for this
    mbinlo = lowest global bin any of my ghost atoms could fall into
    mbinhi = highest global bin any of my ghost atoms could fall into
@@ -979,36 +984,42 @@ void Neighbor::setup_bins()
   bbox[1] = bboxhi[1] - bboxlo[1];
   bbox[2] = bboxhi[2] - bboxlo[2];
 
+  // optimal bin size is roughly 1/2 the cutoff
   // for BIN style, binsize = 1/2 of max neighbor cutoff
   // for MULTI style, binsize = 1/2 of min neighbor cutoff
   // special case of all cutoffs = 0.0, binsize = box size
 
-  double binsize;
+  double binsize_optimal;
-  if (style == BIN) binsize = 0.5*cutneighmax;
+  if (style == BIN) binsize_optimal = 0.5*cutneighmax;
-  else binsize = 0.5*cutneighmin;
+  else binsize_optimal = 0.5*cutneighmin;
-  if (binsize == 0.0) binsize = bbox[0];
+  if (binsize_optimal == 0.0) binsize_optimal = bbox[0];
+  double binsizeinv = 1.0/binsize_optimal;
 
-  // test for too many global bins in any dimension due to huge domain
+  // test for too many global bins in any dimension due to huge global domain
 
-  double binsizeinv = 1.0/binsize;
   if (bbox[0]*binsizeinv > INT_MAX || bbox[1]*binsizeinv > INT_MAX ||
       bbox[2]*binsizeinv > INT_MAX)
     error->all("Domain too large for neighbor bins");
 
-  // divide box into bins
+  // create actual bins
-  // optimal size is roughly 1/2 the cutoff
+  // always have one bin even if cutoff > bbox
-  // use one bin even if cutoff >> bbox
+  // for 2d, nbinz = 1
 
   nbinx = static_cast<int> (bbox[0]*binsizeinv);
   nbiny = static_cast<int> (bbox[1]*binsizeinv);
-  if (domain->dimension == 3)
+  if (dimension == 3) nbinz = static_cast<int> (bbox[2]*binsizeinv);
-    nbinz = static_cast<int> (bbox[2]*binsizeinv);
   else nbinz = 1;
 
   if (nbinx == 0) nbinx = 1;
   if (nbiny == 0) nbiny = 1;
   if (nbinz == 0) nbinz = 1;
 
+  // compute actual bin size for nbins to fit into box exactly
+  // error if actual bin size << cutoff, since will create a zillion bins
+  // this happens when nbin = 1 and box size << cutoff
+  // typically due to non-periodic, flat system in a particular dim
+  // in that extreme case, should use NSQ not BIN neighbor style
+
   binsizex = bbox[0]/nbinx;
   binsizey = bbox[1]/nbiny;
   binsizez = bbox[2]/nbinz;
@@ -1016,7 +1027,12 @@ void Neighbor::setup_bins()
   bininvx = 1.0 / binsizex;
   bininvy = 1.0 / binsizey;
   bininvz = 1.0 / binsizez;
-
+  
+  if (binsize_optimal*bininvx > CUT2BIN_RATIO || 
+      binsize_optimal*bininvy > CUT2BIN_RATIO || 
+      binsize_optimal*bininvz > CUT2BIN_RATIO)
+    error->all("Cannot use neighbor bins - box size << cutoff");
+  
   // mbinlo/hi = lowest and highest global bins my ghost atoms could be in
   // coord = lowest and highest values of coords for my ghost atoms
   // static_cast(-1.5) = -1, so subract additional -1
@@ -1037,13 +1053,16 @@ void Neighbor::setup_bins()
   coord = bsubboxhi[1] + SMALL*bbox[1];
   mbinyhi = static_cast<int> ((coord-bboxlo[1])*bininvy);
 
-  coord = bsubboxlo[2] - SMALL*bbox[2];
+  if (dimension == 3) {
-  mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
+    coord = bsubboxlo[2] - SMALL*bbox[2];
-  if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
+    mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
-  coord = bsubboxhi[2] + SMALL*bbox[2];
+    if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
-  mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
+    coord = bsubboxhi[2] + SMALL*bbox[2];
+    mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
+  }
 
   // extend bins by 1 to insure stencil extent is included
+  // if 2d, only 1 bin in z
 
   mbinxlo = mbinxlo - 1;
   mbinxhi = mbinxhi + 1;
@@ -1053,15 +1072,12 @@ void Neighbor::setup_bins()
   mbinyhi = mbinyhi + 1;
   mbiny = mbinyhi - mbinylo + 1;
 
-  mbinzlo = mbinzlo - 1;
+  if (dimension == 3) {
-  mbinzhi = mbinzhi + 1;
+    mbinzlo = mbinzlo - 1;
+    mbinzhi = mbinzhi + 1;
+  } else mbinzlo = mbinzhi = 0;
   mbinz = mbinzhi - mbinzlo + 1;
 
-  // test for too many total local bins due to huge domain
-
-  if (1.0*mbinx*mbiny*mbinz > INT_MAX)
-    error->all("Domain too large for neighbor bins");
-
   // memory for bin ptrs
 
   mbins = mbinx*mbiny*mbinz;
@@ -1081,6 +1097,7 @@ void Neighbor::setup_bins()
   if (sy*binsizey < cutneighmax) sy++;
   int sz = static_cast<int> (cutneighmax*bininvz);
   if (sz*binsizez < cutneighmax) sz++;
+  if (dimension == 2) sz = 0;
 
   // allocate stencil memory and create stencil(s)
   // check half/full instead of half_every/full_every so stencils will be
@@ -1456,6 +1473,14 @@ void Neighbor::bin_atoms()
     bins[i] = binhead[ibin];
     binhead[ibin] = i;
   }
+
+  /*
+  for (i = nall-1; i >= 0; i--) {
+    ibin = coord2bin(x[i]);
+    bins[i] = binhead[ibin];
+    binhead[ibin] = i;
+  }
+  */
 }
 
 /* ----------------------------------------------------------------------

src/neighbor.h

@@ -120,6 +120,7 @@ class Neighbor : protected Pointers {
   double binsizex,binsizey,binsizez;  // bin sizes and inverse sizes
   double bininvx,bininvy,bininvz;
 
+  int dimension;                   // 2/3 for 2d/3d
   int triclinic;                   // 0 if domain is orthog, 1 if triclinic
   double *bboxlo,*bboxhi;          // copy of full domain bounding box