proposed SJP changes for 2nd cutoff_triple

doc/src/pair_atm.txt

@@ -10,16 +10,20 @@ pair_style atm command :h3
 
 [Syntax:]
 
-pair_style atm cutoff :pre
+SJP: add an arg
 
-cutoff = global cutoff for 3-body interactions (distance units) :ul
+pair_style atm cutoff cutoff_triple :pre
+
+cutoff = cutoff for each pair in 3-body interaction (distance units)
+cutoff_triple = additional cutoff applied to product of 3 pairwise distances (distance units) :ul
 
 [Examples:]
 
-pair_style atm 2.5
+SJP: is 1.5 a good value?
+pair_style atm 2.5 1.5
 pair_coeff * * * 0.072 :pre
 
-pair_style hybrid/overlay lj/cut 6.5 atm 2.5
+pair_style hybrid/overlay lj/cut 6.5 atm 2.5 1.5
 pair_coeff * * lj/cut 1.0 1.0
 pair_coeff 1 1 atm 1 0.064
 pair_coeff 1 1 atm 2 0.080
@@ -50,7 +54,10 @@ potential using the "pair_style hybrid/overlay"_pair_hybrid.html
 command as in the example above.
 
 The potential for a triplet of atom is calculated only if all 3
-distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff
+distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff.
+SJP: In addition, the product of the 3 distances r12*r23*r31 <
+cutoff_triple^3 is required, which limits the contribution of the
+potential to ???
 
 The following coefficients must be defined for each pair of atoms
 types via the "pair_coeff"_pair_coeff.html command as in the examples

src/MANYBODY/pair_atm.cpp

@@ -92,6 +92,10 @@ void PairATM::compute(int eflag, int vflag)
   int *type = atom->type;
 
   double cutoff_squared = cut_global*cut_global;
+  // SJP: new cutoff
+  double triple = cut_triple*cut_triple*cut_triple;
+  double cutoff_triple_sixth = triple*triple;
+
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -138,6 +142,10 @@ void PairATM::compute(int eflag, int vflag)
         rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
         if (rik2 > cutoff_squared) continue;
 
+        // SJP: add this line?
+
+        if (rij2*rjk2*rik2 > cutoff_triple_sixth) continue;
+
         nu_local = nu[type[i]][type[j]][type[k]];
         if (nu_local == 0.0) continue;
 
@@ -192,9 +200,14 @@ void PairATM::allocate()
 
 void PairATM::settings(int narg, char **arg)
 {
-  if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+  // SJP: change to 2 args, require triple <= global
+
+  if (narg != 2) error->all(FLERR,"Illegal pair_style command");
 
   cut_global = force->numeric(FLERR,arg[0]);
+  cut_triple = force->numeric(FLERR,arg[1]);
+
+  if (cut_triple > cut_global) error->all(FLERR,"Illegal pair_style command");
 }
 
 /* ----------------------------------------------------------------------
@@ -311,6 +324,8 @@ void PairATM::read_restart(FILE *fp)
 void PairATM::write_restart_settings(FILE *fp)
 {
   fwrite(&cut_global,sizeof(double),1,fp);
+  // SJP: 2nd arg
+  fwrite(&cut_triple,sizeof(double),1,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -319,9 +334,14 @@ void PairATM::write_restart_settings(FILE *fp)
 
 void PairATM::read_restart_settings(FILE *fp)
 {
+  // SJP: 2nd arg
   int me = comm->me;
-  if (me == 0) fread(&cut_global,sizeof(double),1,fp);
+  if (me == 0) {
+    fread(&cut_global,sizeof(double),1,fp);
+    fread(&cut_triple,sizeof(double),1,fp);
+  }
   MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&cut_triple,1,MPI_DOUBLE,0,world);
 }
 
 /* ----------------------------------------------------------------------

src/MANYBODY/pair_atm.h

@@ -39,7 +39,8 @@ class PairATM : public Pair {
   void read_restart_settings(FILE *);
 
  protected:
-  double cut_global, cutoff_squared;
+  // SJP: add 2nd cutoff
+  double cut_global,cut_triple;
   double ***nu;
 
   void allocate();