diff --git a/doc/src/pair_atm.txt b/doc/src/pair_atm.txt index 2e961a1356..c2718973ad 100644 --- a/doc/src/pair_atm.txt +++ b/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 diff --git a/src/MANYBODY/pair_atm.cpp b/src/MANYBODY/pair_atm.cpp index 9afa36b3c0..e960398ab3 100644 --- a/src/MANYBODY/pair_atm.cpp +++ b/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); } /* ---------------------------------------------------------------------- diff --git a/src/MANYBODY/pair_atm.h b/src/MANYBODY/pair_atm.h index 60d22348c3..cc868df832 100644 --- a/src/MANYBODY/pair_atm.h +++ b/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();