refactor communication cutoff estimator. move it to Comm class. and adjust heuristics.

src/comm.cpp

@@ -20,6 +20,7 @@
 #include "atom_vec.h"
 #include "force.h"
 #include "pair.h"
+#include "bond.h"
 #include "modify.h"
 #include "fix.h"
 #include "compute.h"
@@ -585,6 +586,41 @@ void Comm::set_proc_grid(int outflag)
   }
 }
 
+/* ----------------------------------------------------------------------
+   determine suitable communication cutoff.
+   it is the maximum of the user specified value and estimates based on
+   the maximum neighbor list cutoff and largest bond equilibrium length.
+   we use the 1.5x the bond equilibrium distance as cutoff, if only a
+   bond style exists or only bond and angle styles exists.  If dihedrals
+   or impropers are present we multiply by 2.0. This plus the
+   "neighbor list skin" will become the default communication cutoff, if
+   no pair style is defined and thus avoids all kinds of unexpected behavior
+   for such systems.  If a pair style exists, the result is the maximum of
+   the bond based cutoff and the largest pair cutoff and the user
+   specified communication cutoff.
+------------------------------------------------------------------------- */
+
+double Comm::get_comm_cutoff()
+{
+  double maxcommcutoff = 0.0;
+  if (force->bond) {
+    int n = atom->nbondtypes;
+    for (int i = 1; i <= n; ++i)
+      maxcommcutoff = MAX(maxcommcutoff,force->bond->equilibrium_distance(i));
+
+    if (force->dihedral || force->improper) {
+      maxcommcutoff *= 2.0;
+    } else {
+      maxcommcutoff *=1.5;
+    }
+    maxcommcutoff += neighbor->skin;
+  }
+  maxcommcutoff = MAX(maxcommcutoff,neighbor->cutneighmax);
+  maxcommcutoff = MAX(maxcommcutoff,cutghostuser);
+
+  return maxcommcutoff;
+}
+
 /* ----------------------------------------------------------------------
    determine which proc owns atom with coord x[3] based on current decomp
    x will be in box (orthogonal) or lamda coords (triclinic)
@@ -963,11 +999,6 @@ rendezvous_all2all(int n, char *inbuf, int insize, int inorder, int *procs,
     return 0;    // all nout_rvous are 0, no 2nd irregular
   }
 
-
-
-
-
-
   // create procs and outbuf for All2all if necesary
 
   if (!outorder) {

src/comm.h

@@ -70,6 +70,8 @@ class Comm : protected Pointers {
   void set_processors(int, char **);      // set 3d processor grid attributes
   virtual void set_proc_grid(int outflag = 1); // setup 3d grid of procs
 
+  double get_comm_cutoff();     // determine communication cutoff
+
   virtual void setup() = 0;                      // setup 3d comm pattern
   virtual void forward_comm(int dummy = 0) = 0;  // forward comm of atom coords
   virtual void reverse_comm() = 0;               // reverse comm of forces

src/comm_brick.cpp

@@ -175,7 +175,7 @@ void CommBrick::setup()
   int ntypes = atom->ntypes;
   double *prd,*sublo,*subhi;
 
-  double cut = MAX(neighbor->cutneighmax,cutghostuser);
+  double cut = get_comm_cutoff();
   if ((cut == 0.0) && (me == 0))
     error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
                    "will be generated. Atoms may get lost.");

src/comm_tiled.cpp

@@ -157,7 +157,7 @@ void CommTiled::setup()
   // set cutoff for comm forward and comm reverse
   // check that cutoff < any periodic box length
 
-  double cut = MAX(neighbor->cutneighmax,cutghostuser);
+  double cut = get_comm_cutoff();
   if ((cut == 0.0) && (me == 0))
     error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
                    "will be generated. Atoms may get lost.");

src/info.cpp

@@ -366,7 +366,7 @@ void Info::command(int narg, char **arg)
     if (comm->mode == 0) {
       fprintf(out,"Communication mode = single\n");
       fprintf(out,"Communication cutoff = %g\n",
-              MAX(comm->cutghostuser,neighbor->cutneighmax));
+              comm->get_comm_cutoff());
     }
 
     if (comm->mode == 1) {

src/neighbor.cpp

@@ -35,7 +35,6 @@
 #include "comm.h"
 #include "force.h"
 #include "pair.h"
-#include "bond.h"
 #include "domain.h"
 #include "group.h"
 #include "modify.h"
@@ -277,31 +276,6 @@ void Neighbor::init()
   // cutneigh = force cutoff + skin if cutforce > 0, else cutneigh = 0
   // cutneighghost = pair cutghost if it requests it, else same as cutneigh
 
-  // also consider bonded interactions for estimating the the neighborlist
-  // and communication cutoff. we use the 1.5x the bond equilibrium distance
-  // as cutoff, if only a bond style exists. if also an angle style exists we
-  // multiply by 2.5, for dihedral or improper we multiply by 3.5. (1,2, or 3
-  // bonds plus half a bond length total stretch).
-  // this plus "skin" will become the default communication cutoff, if no
-  // pair style is defined. otherwise the maximum of the largest pairwise
-  // cutoff of this is used.
-
-  double maxbondcutoff = 0.0;
-  if (force->bond) {
-    n = atom->nbondtypes;
-    for (i = 1; i <= n; ++i) {
-       double bondcutoff = force->bond->equilibrium_distance(i);
-       maxbondcutoff = MAX(bondcutoff,maxbondcutoff);
-    }
-    if (force->dihedral || force->improper) {
-      maxbondcutoff *= 3.5;
-    } else if (force->angle) {
-      maxbondcutoff *=2.5;
-    } else {
-      maxbondcutoff *=1.5;
-    }
-  }
-
   triggersq = 0.25*skin*skin;
   boxcheck = 0;
   if (domain->box_change && (domain->xperiodic || domain->yperiodic ||
@@ -319,7 +293,7 @@ void Neighbor::init()
 
   double cutoff,delta,cut;
   cutneighmin = BIG;
-  cutneighmax = maxbondcutoff;
+  cutneighmax = 0.0;
 
   for (i = 1; i <= n; i++) {
     cuttype[i] = cuttypesq[i] = 0.0;