From cd5a34021d03bffefeb3ad256e9eaf6de41f143b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 21 May 2010 15:23:59 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4172
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/neigh_full.cpp       | 18 +++++++++++++++---
 src/neigh_half_bin.cpp   | 20 ++++++++++++++++----
 src/neigh_half_multi.cpp | 20 ++++++++++++++++----
 src/neigh_half_nsq.cpp   | 17 ++++++++++++++---
 src/neigh_respa.cpp      | 37 ++++++++++++++++++++++++++++++-------
 src/neighbor.cpp         | 37 -------------------------------------
 src/neighbor.h           | 36 ++++++++++++++++++++++++++++++++++--
 7 files changed, 125 insertions(+), 60 deletions(-)

diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp
index 31c753579c..de5d576dfc 100644
--- a/src/neigh_full.cpp
+++ b/src/neigh_full.cpp
@@ -30,6 +30,10 @@ void Neighbor::full_nsq(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -81,7 +85,7 @@ void Neighbor::full_nsq(NeighList *list)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -116,6 +120,10 @@ void Neighbor::full_bin(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -170,7 +178,7 @@ void Neighbor::full_bin(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -208,6 +216,10 @@ void Neighbor::full_multi(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -268,7 +280,7 @@ void Neighbor::full_multi(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp
index 570b87307e..3432b4a0d1 100644
--- a/src/neigh_half_bin.cpp
+++ b/src/neigh_half_bin.cpp
@@ -37,6 +37,10 @@ void Neighbor::half_bin_no_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -93,7 +97,7 @@ void Neighbor::half_bin_no_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -130,6 +134,10 @@ void Neighbor::half_bin_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -188,7 +196,7 @@ void Neighbor::half_bin_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -209,7 +217,7 @@ void Neighbor::half_bin_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -246,6 +254,10 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -309,7 +321,7 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp
index bbfff83655..b61a53114d 100644
--- a/src/neigh_half_multi.cpp
+++ b/src/neigh_half_multi.cpp
@@ -39,6 +39,10 @@ void Neighbor::half_multi_no_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -101,7 +105,7 @@ void Neighbor::half_multi_no_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -140,6 +144,10 @@ void Neighbor::half_multi_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -199,7 +207,7 @@ void Neighbor::half_multi_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -227,7 +235,7 @@ void Neighbor::half_multi_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -266,6 +274,10 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -337,7 +349,7 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp
index 1f3bf89ae8..0447c3ffd4 100644
--- a/src/neigh_half_nsq.cpp
+++ b/src/neigh_half_nsq.cpp
@@ -31,6 +31,12 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
+  // loop over each atom, storing neighbors
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -81,7 +87,7 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -111,8 +117,13 @@ void Neighbor::half_nsq_newton(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr;
 
-  double **x = atom->x;
+  // loop over each atom, storing neighbors
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
   int *tag = atom->tag;
+
+  double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
   int *molecule = atom->molecule;
@@ -180,7 +191,7 @@ void Neighbor::half_nsq_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp
index daed5d6793..a1e0194a33 100644
--- a/src/neigh_respa.cpp
+++ b/src/neigh_respa.cpp
@@ -32,6 +32,12 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr,*neighptr_inner,*neighptr_middle;
 
+  // loop over each atom, storing neighbors
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -122,7 +128,7 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -184,8 +190,13 @@ void Neighbor::respa_nsq_newton(NeighList *list)
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *neighptr,*neighptr_inner,*neighptr_middle;
 
-  double **x = atom->x;
+  // loop over each atom, storing neighbors
+
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
   int *tag = atom->tag;
+
+  double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
   int *molecule = atom->molecule;
@@ -292,7 +303,7 @@ void Neighbor::respa_nsq_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -361,6 +372,10 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -456,7 +471,7 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -525,6 +540,10 @@ void Neighbor::respa_bin_newton(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -623,7 +642,7 @@ void Neighbor::respa_bin_newton(NeighList *list)
       rsq = delx*delx + dely*dely + delz*delz;
 
       if (rsq <= cutneighsq[itype][jtype]) {
-	if (molecular) which = find_special(i,j);
+	if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	else which = 0;
 	if (which == 0) neighptr[n++] = j;
 	else if (which > 0) neighptr[n++] = which*nall + j;
@@ -655,7 +674,7 @@ void Neighbor::respa_bin_newton(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
@@ -724,6 +743,10 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
 
   // loop over each atom, storing neighbors
 
+  int **special = atom->special;
+  int **nspecial = atom->nspecial;
+  int *tag = atom->tag;
+
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
@@ -827,7 +850,7 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
 	rsq = delx*delx + dely*dely + delz*delz;
 
 	if (rsq <= cutneighsq[itype][jtype]) {
-	  if (molecular) which = find_special(i,j);
+	  if (molecular) which = find_special(special[i],nspecial[i],tag[j]);
 	  else which = 0;
 	  if (which == 0) neighptr[n++] = j;
 	  else if (which > 0) neighptr[n++] = which*nall + j;
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 1aa2743687..4cd449659f 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1443,43 +1443,6 @@ void Neighbor::modify_params(int narg, char **arg)
   }
 }
 
-/* ----------------------------------------------------------------------
-   determine if atom j is in special list of atom i
-   if it is not, return 0
-   if it is and special flag is 0 (both coeffs are 0.0), return -1
-   if it is and special flag is 1 (both coeffs are 1.0), return 0
-   if it is and special flag is 2 (otherwise), return 1,2,3
-     for which neighbor it is (and which coeff it maps to)
-------------------------------------------------------------------------- */
-
-int Neighbor::find_special(int i, int j)
-{
-  int *list = atom->special[i];
-  int n1 = atom->nspecial[i][0];
-  int n2 = atom->nspecial[i][1];
-  int n3 = atom->nspecial[i][2];
-  int tag = atom->tag[j];
-
-  for (int i = 0; i < n3; i++) {
-    if (list[i] == tag) {
-      if (i < n1) {
-	if (special_flag[1] == 0) return -1;
-	else if (special_flag[1] == 1) return 0;
-	else return 1;
-      } else if (i < n2) {
-	if (special_flag[2] == 0) return -1;
-	else if (special_flag[2] == 1) return 0;
-	else return 2;
-      } else {
-	if (special_flag[3] == 0) return -1;
-	else if (special_flag[3] == 1) return 0;
-	else return 3;
-      }
-    }
-  }
-  return 0;
-}
-
 /* ----------------------------------------------------------------------
    bin owned and ghost atoms
 ------------------------------------------------------------------------- */
diff --git a/src/neighbor.h b/src/neighbor.h
index 5f5951fc37..6a21950af6 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -149,11 +149,44 @@ class Neighbor : protected Pointers {
   void bin_atoms();                     // bin all atoms
   double bin_distance(int, int, int);   // distance between binx
   int coord2bin(double *);              // mapping atom coord to a bin
-  int find_special(int, int);           // look for special neighbor
+
   int exclusion(int, int, int, int, int *, int *);  // test for pair exclusion
   void choose_build(int, class NeighRequest *);
   void choose_stencil(int, class NeighRequest *);
 
+  // find_special: determine if atom j is in special list of atom i
+  // if it is not, return 0
+  // if it is and special flag is 0 (both coeffs are 0.0), return -1
+  // if it is and special flag is 1 (both coeffs are 1.0), return 0
+  // if it is and special flag is 2 (otherwise), return 1,2,3
+  //   for which neighbor it is (and which coeff it maps to)
+
+  inline int find_special(const int *list, const int *nspecial, 
+			  const int tag) const {
+    const int n1 = nspecial[0];
+    const int n2 = nspecial[1];
+    const int n3 = nspecial[2];
+
+    for (int i = 0; i < n3; i++) {
+      if (list[i] == tag) {
+	if (i < n1) {
+	  if (special_flag[1] == 0) return -1;
+	  else if (special_flag[1] == 1) return 0;
+	  else return 1;
+	} else if (i < n2) {
+	  if (special_flag[2] == 0) return -1;
+	  else if (special_flag[2] == 1) return 0;
+	  else return 2;
+	} else {
+	  if (special_flag[3] == 0) return -1;
+	  else if (special_flag[3] == 1) return 0;
+	  else return 3;
+	}
+      }
+    }
+    return 0;
+  };
+
   // pairwise build functions
 
   typedef void (Neighbor::*PairPtr)(class NeighList *);
@@ -236,7 +269,6 @@ class Neighbor : protected Pointers {
   BondPtr improper_build;             // ptr to improper list functions
   void improper_all();                // improper list with all impropers
   void improper_partial();            // exclude certain impropers
-
 };
 
 }
-- 
GitLab