diff --git a/src/comm.cpp b/src/comm.cpp index 052de93793..3b4029e7d3 100644 --- a/src/comm.cpp +++ b/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) { diff --git a/src/comm.h b/src/comm.h index 90a215b4a6..860912a7d2 100644 --- a/src/comm.h +++ b/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 diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index eb53f89746..a9ce64399f 100644 --- a/src/comm_brick.cpp +++ b/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."); diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index 9d4c446057..2558559cf7 100644 --- a/src/comm_tiled.cpp +++ b/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."); diff --git a/src/info.cpp b/src/info.cpp index 25b9879408..74cbe6556f 100644 --- a/src/info.cpp +++ b/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) { diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 9964b23b60..0382624198 100644 --- a/src/neighbor.cpp +++ b/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;