From 33388e99af42ed9079cab60116a6bf27f99ac9a5 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 3 Jun 2010 20:15:41 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4228
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/GRANULAR/pair_gran_hooke_history.cpp | 85 ++++++++++++------------
 src/GRANULAR/pair_gran_hooke_history.h   |  3 +
 src/pair.cpp                             |  3 +-
 3 files changed, 48 insertions(+), 43 deletions(-)

diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp
index 196934139d..ccdcf9883c 100644
--- a/src/GRANULAR/pair_gran_hooke_history.cpp
+++ b/src/GRANULAR/pair_gran_hooke_history.cpp
@@ -60,6 +60,11 @@ PairGranHookeHistory::~PairGranHookeHistory()
   if (allocated) {
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
+
+    delete [] onerad_dynamic;
+    delete [] onerad_frozen;
+    delete [] maxrad_dynamic;
+    delete [] maxrad_frozen;
   }
 }
 
@@ -280,6 +285,11 @@ void PairGranHookeHistory::allocate()
       setflag[i][j] = 0;
 
   cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+
+  onerad_dynamic = new double[n+1];
+  onerad_frozen = new double[n+1];
+  maxrad_dynamic = new double[n+1];
+  maxrad_frozen = new double[n+1];
 }
 
 /* ----------------------------------------------------------------------
@@ -379,45 +389,46 @@ void PairGranHookeHistory::init_style()
     fix_history->pair = this;
   }
 
-  // check for freeze Fix and set freeze_group_bit
+  // check for Fix freeze and set freeze_group_bit
 
   for (i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"freeze") == 0) break;
   if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit;
   else freeze_group_bit = 0;
 
-  // set cutoff by largest particles
-  // maxrad_dynamic = radius of largest dynamic particle, including inserted
-  // maxrad_frozen = radius of largest dynamic particle
-  // include frozen-dynamic interactions
-  // do not include frozen-frozen interactions
-  // include future inserted particles as dynamic
-  // cutforce was already set in pair::init(), but this sets it correctly
+  // check for Fix pour and set pour_type and pour_maxdiam
+
+  int pour_type = 0;
+  double pour_maxrad = 0.0;
+  for (i = 0; i < modify->nfix; i++)
+    if (strcmp(modify->fix[i]->style,"pour") == 0) break;
+  if (i < modify->nfix) {
+    pour_type = ((FixPour *) modify->fix[i])->ntype;
+    pour_maxrad = ((FixPour *) modify->fix[i])->radius_hi;
+  }
+
+  // set maxrad_dynamic and maxrad_frozen for each type
+  // include future Fix pour particles as dynamic
+
+  for (i = 1; i <= atom->ntypes; i++)
+    onerad_dynamic[i] = onerad_frozen[i] = 0.0;
+  if (pour_type) onerad_dynamic[pour_type] = pour_maxrad;
 
   double *radius = atom->radius;
   int *mask = atom->mask;
+  int *type = atom->type;
   int nlocal = atom->nlocal;
 
-  double maxrad_dynamic = 0.0;
-  for (i = 0; i < nlocal; i++)
-    if (!(mask[i] & freeze_group_bit))
-      maxrad_dynamic = MAX(maxrad_dynamic,radius[i]);
-  double mine = maxrad_dynamic;
-  MPI_Allreduce(&mine,&maxrad_dynamic,1,MPI_DOUBLE,MPI_MAX,world);
-
-  for (i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"pour") == 0)
-      maxrad_dynamic =
-	MAX(maxrad_dynamic,((FixPour *) modify->fix[i])->radius_hi);
-
-  double maxrad_frozen = 0.0;
   for (i = 0; i < nlocal; i++)
     if (mask[i] & freeze_group_bit)
-      maxrad_frozen = MAX(maxrad_frozen,radius[i]);
-  mine = maxrad_frozen;
-  MPI_Allreduce(&mine,&maxrad_frozen,1,MPI_DOUBLE,MPI_MAX,world);
-
-  cutforce = maxrad_dynamic + MAX(maxrad_dynamic,maxrad_frozen);
+      onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]);
+    else
+      onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
+  
+  MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,
+		MPI_DOUBLE,MPI_MAX,world);
+  MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,
+		MPI_DOUBLE,MPI_MAX,world);
 }
 
 /* ----------------------------------------------------------------------
@@ -439,23 +450,13 @@ double PairGranHookeHistory::init_one(int i, int j)
 {
   if (!allocated) allocate();
 
-  // return max diameter of any particle
-  // not used in granular neighbor calculation since particles have radius
-  // but will insure cutoff is big enough for any other neighbor lists built
-  // if no particles, return 1.0
+  // cutoff = sum of max I,J radii for
+  // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
 
-  double *radius = atom->radius;
-  int nlocal = atom->nlocal;
-
-  double maxrad = 0.0;
-  for (int m = 0; m < nlocal; m++) maxrad = MAX(maxrad,radius[m]);
-  double mine = 2.0*maxrad;
-
-  double maxdiam;
-  MPI_Allreduce(&mine,&maxdiam,1,MPI_DOUBLE,MPI_MAX,world);
-  if (maxdiam == 0.0) maxdiam = 1.0;
-
-  return maxdiam;
+  double cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
+  cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
+  cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
+  return cutoff;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h
index 7252721d12..f005aa4e7c 100644
--- a/src/GRANULAR/pair_gran_hooke_history.h
+++ b/src/GRANULAR/pair_gran_hooke_history.h
@@ -48,6 +48,9 @@ class PairGranHookeHistory : public Pair {
   int history;
   class FixShearHistory *fix_history;
 
+  double *onerad_dynamic,*onerad_frozen;
+  double *maxrad_dynamic,*maxrad_frozen;
+
   void allocate();
 };
 
diff --git a/src/pair.cpp b/src/pair.cpp
index 73daa5eeff..d72566b4e2 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -154,9 +154,10 @@ void Pair::init()
   // set cutsq for each I,J, used to neighbor
   // cutforce = max of all I,J cutoffs
 
-  double cut;
   cutforce = 0.0;
   etail = ptail = 0.0;
+  double cut;
+
   for (i = 1; i <= atom->ntypes; i++)
     for (j = i; j <= atom->ntypes; j++) {
       cut = init_one(i,j);