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

src/delete_atoms.cpp

@@ -24,8 +24,7 @@
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "comm.h"
+#include "memory.h"
-#include "force.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -55,15 +54,14 @@ void DeleteAtoms::command(int narg, char **arg)
   // allocate and initialize deletion list
   
   int nlocal = atom->nlocal;
-  int *dlist = new int[nlocal];
+  dlist = (int *) memory->smalloc(nlocal*sizeof(int),"delete_atoms:dlist");
-
   for (int i = 0; i < nlocal; i++) dlist[i] = 0;
 
   // delete the atoms
 
-  if (strcmp(arg[0],"group") == 0) delete_group(narg,arg,dlist);
+  if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
-  else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg,dlist);
+  else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
-  else if (strcmp(arg[0],"overlap") == 0) delete_overlap(narg,arg,dlist);
+  else if (strcmp(arg[0],"overlap") == 0) delete_overlap(narg,arg);
   else error->all("Illegal delete_atoms command");
 
   // delete local atoms in list
@@ -80,7 +78,7 @@ void DeleteAtoms::command(int narg, char **arg)
     } else i++;
   }
   atom->nlocal = nlocal;
-  delete [] dlist;
+  memory->sfree(dlist);
 
   // if non-molecular system, reset atom tags to be contiguous
   // set all atom IDs to 0, call tag_extend()
@@ -122,7 +120,7 @@ void DeleteAtoms::command(int narg, char **arg)
    group will still exist
 ------------------------------------------------------------------------- */
 
-void DeleteAtoms::delete_group(int narg, char **arg, int *dlist)
+void DeleteAtoms::delete_group(int narg, char **arg)
 {
   if (narg != 2) error->all("Illegal delete_atoms command");
 
@@ -141,7 +139,7 @@ void DeleteAtoms::delete_group(int narg, char **arg, int *dlist)
    delete all atoms in region
 ------------------------------------------------------------------------- */
 
-void DeleteAtoms::delete_region(int narg, char **arg, int *dlist)
+void DeleteAtoms::delete_region(int narg, char **arg)
 {
   if (narg != 2) error->all("Illegal delete_atoms command");
   
@@ -162,34 +160,33 @@ void DeleteAtoms::delete_region(int narg, char **arg, int *dlist)
    no guarantee that minimium number of atoms will be deleted
 ------------------------------------------------------------------------- */
 
-void DeleteAtoms::delete_overlap(int narg, char **arg, int *dlist)
+void DeleteAtoms::delete_overlap(int narg, char **arg)
 {
-  if (narg < 2) error->all("Illegal delete_atoms command");
+  if (narg < 4) error->all("Illegal delete_atoms command");
     
-  // read args including optional type info
+  // read args
 
   double cut = atof(arg[1]);
   double cutsq = cut*cut;
 
-  int typeflag,type1,type2;
+  int igroup1 = group->find(arg[2]);
-  if (narg == 2) typeflag = 0;
+  int igroup2 = group->find(arg[3]);
-  else if (narg == 3) {
+  if (igroup1 < 0 || igroup2 < 0)
-    typeflag = 1;
+    error->all("Could not find delete_atoms group ID");
-    type1 = atoi(arg[2]);
+
-  } else if (narg == 4) {
+  int group1bit = group->bitmask[igroup1];
-    typeflag = 2;
+  int group2bit = group->bitmask[igroup2];
-    type1 = atoi(arg[2]);
-    type2 = atoi(arg[3]);
-  } else error->all("Illegal delete_atoms command");
 
   if (comm->me == 0 && screen)
     fprintf(screen,"System init for delete_atoms ...\n");
 
-  // request a half neighbor list for use by this command
+  // request a full neighbor list for use by this command
 
   int irequest = neighbor->request((void *) this);
   neighbor->requests[irequest]->pair = 0;
   neighbor->requests[irequest]->command = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
   neighbor->requests[irequest]->occasional = 1;
 
   // init entire system since comm->borders and neighbor->build is done
@@ -218,30 +215,13 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *dlist)
   NeighList *list = neighbor->lists[irequest];
   neighbor->build_one(irequest);
 
-  // create an atom map if one doesn't exist already
+  // double loop over owned atoms and their full neighbor list
-
-  int mapflag = 0;
-  if (atom->map_style == 0) {
-    mapflag = 1;
-    atom->map_style = 1;
-    atom->map_init();
-    atom->map_set();
-  }
-
-  // double loop over owned atoms and their neighbors
   // at end of loop, there are no more overlaps
-  // criteria for i,j to undergo a deletion event:
+  // only ever delete owned atom I, never J even if owned
-  //   weighting factor != 0.0 for this pair
-  //     could be 0 and still be in neigh list for long-range Coulombics
-  //   local copy of j (map[tag[j]]) has not already been deleted
-  //   distance between i,j is less than cutoff
-  //   i,j are of valid types
-  // if all criteria met, delete i and skip to next i in outer loop
-  //   unless j is ghost and newton_pair off and tag[j] < tag[i]
-  //   then rely on other proc to delete
 
   int *tag = atom->tag;
   int *type = atom->type;
+  int *mask = atom->mask;
   double **x = atom->x;
   int nlocal = atom->nlocal;
   int nall = atom->nlocal + atom->nghost;
@@ -270,17 +250,16 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *dlist)
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
 
+      // if weighting factors are 0, skip this pair
+      // could be 0 and still be in neigh list for long-range Coulombics
+      // want consistency with non-charged pairs which wouldn't be in list
+
       if (j >= nall) {
 	if (special_coul[j/nall] == 0.0 && special_lj[j/nall] == 0.0) continue;
 	j %= nall;
       }
 
-      if (j < nlocal) {
+      // only consider deletion if I,J distance < cutoff
-	if (dlist[j]) continue;
-      } else {
-	m = atom->map(tag[j]);
-	if (m < nlocal && dlist[m]) continue;
-      }
 
       delx = xtmp - x[j][0];
       dely = ytmp - x[j][1];
@@ -288,24 +267,30 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *dlist)
       rsq = delx*delx + dely*dely + delz*delz;
       if (rsq >= cutsq) continue;
 
-      if (typeflag) {
+      // only consider deletion if I,J are in groups 1,2 respectively
-	jtype = type[j];
+      // true whether J is owned or ghost atom
-	if (typeflag == 1 && itype != type1 && jtype != type1) continue;
-	if (typeflag == 2 && !(itype == type1 && jtype == type2) && 
-	    !(itype == type2 && jtype == type1)) continue;
-      }
 
-      if (j >= nlocal && newton_pair == 0 && tag[j] < tag[i]) continue;
+      if (!(mask[i] & group1bit)) continue;
+      if (!(mask[j] & group2bit)) continue;
+
+      // J is owned atom:
+      //   delete atom I if atom J has not already been deleted
+      // J is ghost atom:
+      //   delete atom I if J,I is not a candidate deletion pair
+      //     due to being in groups 1,2 respectively
+      //   if they are candidate pair, then either:
+      //      another proc owns J and could delete J
+      //      J is a ghost of another of my owned atoms, and I could delete J
+      //   test on tags of I,J insures that only I or J is deleted
+
+      if (j < nlocal) {
+	if (dlist[j]) continue;
+      } else if ((mask[i] & group2bit) && (mask[j] & group1bit)) {
+	if (tag[i] > tag[j]) continue;
+      }
 
       dlist[i] = 1;
       break;
     }
   }
-
-  // delete temporary atom map
-
-  if (mapflag) {
-    atom->map_delete();
-    atom->map_style = 0;
-  }
 }

src/delete_atoms.h

@@ -24,9 +24,11 @@ class DeleteAtoms : protected Pointers {
   void command(int, char **);
 
  private:
-  void delete_group(int, char **, int *);
+  int *dlist;
-  void delete_region(int, char **, int *);
+
-  void delete_overlap(int, char **, int *);
+  void delete_group(int, char **);
+  void delete_region(int, char **);
+  void delete_overlap(int, char **);
 };
 
 }

src/finish.cpp

@@ -294,6 +294,7 @@ void Finish::end(int flag)
   }
 
   // find a non-skip neighbor list containing half the pairwise interactions
+  // count neighbors in that list for stats purposes
 
   for (m = 0; m < neighbor->old_nrequest; m++)
     if ((neighbor->old_requests[m]->half || neighbor->old_requests[m]->gran ||
@@ -302,9 +303,13 @@ void Finish::end(int flag)
 	neighbor->old_requests[m]->skip == 0) break;
 
   int nneigh = 0;
-  if (m < neighbor->old_nrequest)
+  if (m < neighbor->old_nrequest) {
-    for (i = 0; i < atom->nlocal; i++)
+    int inum = neighbor->lists[m]->inum;
-      nneigh += neighbor->lists[m]->numneigh[i];
+    int *ilist = neighbor->lists[m]->ilist;
+    int *numneigh = neighbor->lists[m]->numneigh;
+    for (int ii = 0; ii < inum; ii++)
+      nneigh += numneigh[ilist[ii]];
+  }
 
   tmp = nneigh;
   stats(1,&tmp,&ave,&max,&min,10,histo);

src/pair_hybrid.cpp

@@ -342,10 +342,11 @@ void PairHybrid::init_style()
 
   for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();
 
-  // create skip lists for each neigh request
+  // create skip lists for each pair neigh request
   // any kind of list can have its skip flag set at this stage
 
   for (i = 0; i < neighbor->nrequest; i++) {
+    if (!neighbor->requests[i]->pair) continue;
 
     // find associated sub-style
 
@@ -458,6 +459,7 @@ void PairHybrid::modify_requests()
   int i,j;
   NeighRequest *irq,*jrq;
 
+  // loop over pair requests only
   // if list is skip list and not copy, look for non-skip list of same kind
   // if one exists, point at that one
   // else make new non-skip request of same kind and point at that one
@@ -471,6 +473,8 @@ void PairHybrid::modify_requests()
   //   which invokes this routine
 
   for (i = 0; i < neighbor->nrequest; i++) {
+    if (!neighbor->requests[i]->pair) continue;
+
     irq = neighbor->requests[i];
     if (irq->skip == 0 || irq->copy) continue;
     if (irq->half_from_full) {
@@ -479,6 +483,7 @@ void PairHybrid::modify_requests()
     }
 
     for (j = 0; j < neighbor->nrequest; j++) {
+      if (!neighbor->requests[j]->pair) continue;
       jrq = neighbor->requests[j];
       if (irq->same_kind(jrq) && jrq->skip == 0) break;
     }

src/pair_hybrid_overlay.cpp

@@ -101,15 +101,18 @@ void PairHybridOverlay::modify_requests()
   int i,j;
   NeighRequest *irq,*jrq;
 
-  // loop over pair requests
+  // loop over pair requests only
   // if a previous list is same kind with same skip attributes
   // then make this one a copy list of that one
   // works whether both lists are no-skip or yes-skip
   // will not point a list at a copy list, but at copy list's parent
 
   for (i = 0; i < neighbor->nrequest; i++) {
+    if (!neighbor->requests[i]->pair) continue;
+
     irq = neighbor->requests[i];
     for (j = 0; j < i; j++) {
+      if (!neighbor->requests[j]->pair) continue;
       jrq = neighbor->requests[j];
       if (irq->same_kind(jrq) && irq->same_skip(jrq)) {
 	irq->copy = 1;