From a2303e5a15b27a40bf82b03bb1d5bf9f0c3c2b69 Mon Sep 17 00:00:00 2001
From: sergeylishchuk <sergey.lishchuk@gmail.com>
Date: Thu, 26 Jul 2018 11:54:04 +0100
Subject: [PATCH] more changes to ATM source and doc file

---
 doc/src/pair_atm.txt      | 14 +++++++-------
 src/MANYBODY/pair_atm.cpp | 40 +++++++++++++++++++++++++++++++--------
 src/MANYBODY/pair_atm.h   |  5 +++++
 3 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/doc/src/pair_atm.txt b/doc/src/pair_atm.txt
index e442c32814..6453b02bf8 100644
--- a/doc/src/pair_atm.txt
+++ b/doc/src/pair_atm.txt
@@ -41,17 +41,17 @@ gamma3 are as shown in this diagram:
 
 Note that for the interaction between a triplet of atoms I,J,K, there
 is no "central" atom.  The interaction is symmetric with respect to
-permutation of the three atoms.  In fact, it is computed 3 times by
-the code with each of I,J,K being atom 1.  Thus the nu value must be
-the same for all those permutations of the atom types of I,J,K as
-discussed below.
+permutation of the three atoms. Thus the nu value is
+the same for all those permutations of the atom types of I,J,K
+and needs to be specified only once, as discussed below.
 
 The {atm} potential is typically used in combination with a two-body
 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.
+The potential is calculated if atoms J and K are in the
+"neighbor list"_neighbor.html of atom I
+and the distances between atoms satisfy rIJ rJK rKI > cutoff^3.
 
 The following coefficients must be defined for each pair of atoms
 types via the "pair_coeff"_pair_coeff.html command as in the examples
@@ -88,7 +88,7 @@ pair_coeff 2 3 4 0.1 :pre
 Note that for a simulation with a single atom type, only a single
 entry is required, e.g.
 
-pair_coeff 1 1 1 0.25 :per
+pair_coeff 1 1 1 0.25 :pre
 
 For a simulation with two atom types, four pair_coeff commands will
 specify all possible nu values:
diff --git a/src/MANYBODY/pair_atm.cpp b/src/MANYBODY/pair_atm.cpp
index 2ba8912899..090a8352d5 100644
--- a/src/MANYBODY/pair_atm.cpp
+++ b/src/MANYBODY/pair_atm.cpp
@@ -91,8 +91,8 @@ void PairATM::compute(int eflag, int vflag)
   double **f = atom->f;
   int *type = atom->type;
 
-  double cutoffsq = cut_global*cut_global;
-
+  double cutoff_sixth = cut_global*cut_global;
+  cutoff_sixth = cutoff_sixth*cutoff_sixth*cutoff_sixth;
   inum = list->inum;
   ilist = list->ilist;
   numneigh = list->numneigh;
@@ -112,26 +112,32 @@ void PairATM::compute(int eflag, int vflag)
       j = jlist[jj];
       j &= NEIGHMASK;
       rij[0] = x[j][0] - xi;
+      if (rij[0] < 0.0) continue;
       rij[1] = x[j][1] - yi;
+      if (rij[0] == 0.0 and rij[1] < 0.0) continue;
       rij[2] = x[j][2] - zi;
+      if (rij[0] == 0.0 and rij[1] == 0.0 and rij[2] < 0.0) continue;
       rij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
 
       for (kk = jj+1; kk < jnum; kk++) {
         k = jlist[kk];
         k &= NEIGHMASK;
 
-        rik[0] = x[k][0] - xi;
-        rik[1] = x[k][1] - yi;
-        rik[2] = x[k][2] - zi;
-        rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
-
         rjk[0] = x[k][0] - x[j][0];
+        if (rjk[0] < 0.0) continue;
         rjk[1] = x[k][1] - x[j][1];
+        if (rjk[0] == 0.0 and rjk[1] < 0.0) continue;
         rjk[2] = x[k][2] - x[j][2];
+        if (rjk[0] == 0.0 and rjk[1] == 0.0 and rjk[2] < 0.0) continue;
         rjk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2];
 
+        rik[0] = x[k][0] - xi;
+        rik[1] = x[k][1] - yi;
+        rik[2] = x[k][2] - zi;
+        rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
+
         r6 = rij2*rik2*rjk2;
-        if (r6 > cutoffsq) continue;
+        if (r6 > cutoff_sixth) continue;
 
         nu_local = nu[type[i]][type[j]][type[k]];
         if (nu_local == 0.0) continue;
@@ -228,6 +234,9 @@ void PairATM::coeff(int narg, char **arg)
 
 void PairATM::init_style()
 {
+  if (force->newton_pair == 0)
+    error->all(FLERR,"Pair style ATM requires newton pair on");
+
   // need a full neighbor list
 
   int irequest = neighbor->request(this,instance_me);
@@ -258,6 +267,21 @@ double PairATM::init_one(int i, int j)
    proc 0 writes to restart file
 ------------------------------------------------------------------------- */
 
+// 
+// PLEASE HELP!
+// 
+// From the pair_hybrid documentation: 
+// 
+// "For the hybrid pair styles, the list of sub-styles and their respective
+// settings are written to binary restart files, so a pair_style command
+// does not need to specified in an input script that reads a restart file.
+// However, the coefficient information is not stored in the restart file.
+// Thus, pair_coeff commands need to be re-specified in the restart input
+// script."
+// 
+// As a result, the function "PairATM::write_restart" is not called from
+// "pair_hybrid.cpp". What is the best way to deal with this problem?
+// 
 void PairATM::write_restart(FILE *fp)
 {
   write_restart_settings(fp);
diff --git a/src/MANYBODY/pair_atm.h b/src/MANYBODY/pair_atm.h
index 45a6cda344..210351cfad 100644
--- a/src/MANYBODY/pair_atm.h
+++ b/src/MANYBODY/pair_atm.h
@@ -64,6 +64,11 @@ E: Incorrect args for pair coefficients
 
 Self-explanatory.  Check the input script or data file.
 
+E: Pair style ATM requires newton pair on
+
+See the newton command.  This is a restriction to use the ATM
+potential.
+
 E: All pair coeffs are not set
 
 All pair coefficients must be set in the data file or by the
-- 
GitLab