provide keyword/value option to compute ackland/atom for selecting legacy or current variant of implementation

doc/src/compute_ackland_atom.txt

@@ -10,19 +10,29 @@ compute ackland/atom command :h3
 
 [Syntax:]
 
-compute ID group-ID ackland/atom :pre
+compute ID group-ID ackland/atom keyword/value :pre
 
-ID, group-ID are documented in "compute"_compute.html command
-ackland/atom = style name of this compute command :ul
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+ackland/atom = style name of this compute command :l
+
+zero or more keyword/value pairs may be appended :l
+keyword = {legacy} :l
+  {legacy} yes/no = use ({yes}) or do not use ({no}) legacy ackland algorithm implementation :pre 
+:ule
 
 [Examples:]
 
-compute 1 all ackland/atom :pre
+compute 1 all ackland/atom
+compute 1 all ackland/atom legacy yes :pre
 
 [Description:]
 
 Defines a computation that calculates the local lattice structure
 according to the formulation given in "(Ackland)"_#Ackland.
+Historically, LAMMPS had two, slightly different implementations of
+the algorithm from the paper. With the {legacy} keyword, it is
+possible to switch between the pre-2015 ({legacy yes}) and post-2015
+implemention ({legacy no}). The post-2015 variant is the default. 
 
 In contrast to the "centro-symmetry
 parameter"_compute_centro_atom.html this method is stable against
@@ -66,7 +76,8 @@ integers defined above.
 
 "compute centro/atom"_compute_centro_atom.html
 
-[Default:] none
+[Default:] 
+The keyword {legacy} defaults to {no}.
 
 :line
 

src/USER-MISC/compute_ackland_atom.cpp

@@ -14,6 +14,7 @@
 /* ----------------------------------------------------------------------
    Contributing author:  G. Ziegenhain, gerolf@ziegenhain.com
                          Copyright (C) 2007
+   Updated algorithm by: Brian Barnes, brian.c.barnes11.civ@mail.mil
 ------------------------------------------------------------------------- */
 
 #include <cmath>
@@ -40,7 +41,8 @@ enum{UNKNOWN,BCC,FCC,HCP,ICO};
 ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute ackland/atom command");
+  if ((narg < 3) || (narg > 5))
+    error->all(FLERR,"Illegal compute ackland/atom command");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -48,10 +50,26 @@ ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
   nmax = 0;
   structure = NULL;
   maxneigh = 0;
+  legacy = 0;
   distsq = NULL;
   nearest = NULL;
   nearest_n0 = NULL;
   nearest_n1 = NULL;
+
+  int iarg = 3;
+  while (narg > iarg) {
+    if (strcmp("legacy",arg[iarg]) == 0) {
+      ++iarg;
+      if (iarg >= narg)
+        error->all(FLERR,"Invalid compute ackland/atom command");
+      if (strcmp("yes",arg[iarg]) == 0)
+        legacy = 1;
+      else if (strcmp("no",arg[iarg]) == 0)
+        legacy = 0;
+      else error->all(FLERR,"Invalid compute ackland/atom command");
+    }
+    ++iarg;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -231,6 +249,53 @@ void ComputeAcklandAtom::compute_peratom()
           else chi[7]++;
         }
       }
+      if (legacy) {
+
+        // This is the original implementation by Gerolf Ziegenhain
+        // Deviations from the different lattice structures
+
+        double delta_bcc = 0.35*chi[4]/(double)(chi[5]+chi[6]-chi[4]);
+        double delta_cp = fabs(1.-(double)chi[6]/24.);
+        double delta_fcc = 0.61*(fabs((double)(chi[0]+chi[1]-6.))+
+                                 (double)chi[2])/6.0;
+        double delta_hcp = (fabs((double)chi[0]-3.)+
+                            fabs((double)chi[0]+(double)chi[1]+
+                                 (double)chi[2]+(double)chi[3]-9.0))/12.0;
+
+        // Identification of the local structure according to the reference
+
+        if (chi[0] == 7)       { delta_bcc = 0.; }
+        else if (chi[0] == 6)  { delta_fcc = 0.; }
+        else if (chi[0] <= 3)  { delta_hcp = 0.; }
+
+        if (chi[7] > 0.)
+          structure[i] = UNKNOWN;
+        else
+          if (chi[4] < 3.)
+            {
+              if (n1 > 13 || n1 < 11)
+                structure[i] = UNKNOWN;
+              else
+                structure[i] = ICO;
+            } else
+            if (delta_bcc <= delta_cp)
+              {
+                if (n1 < 11)
+                  structure[i] = UNKNOWN;
+                else
+                  structure[i] = BCC;
+              } else
+              if (n1 > 12 || n1 < 11)
+                structure[i] = UNKNOWN;
+              else
+                if (delta_fcc < delta_hcp)
+                  structure[i] = FCC;
+                else
+                  structure[i] = HCP;
+
+      } else {
+
+        // This is the updated implementation by Brian Barnes
         
         if (chi[7] > 0 || n0 < 11) structure[i] = UNKNOWN;
         else if (chi[0] == 7) structure[i] = BCC;
@@ -252,8 +317,11 @@ void ComputeAcklandAtom::compute_peratom()
           double delta_fcc = 0.61*(fabs((double)(chi[0]+chi[1]-6))
                                    +(double)chi[2])/6.0;
 
-        double delta_hcp = (fabs((double)chi[0]-3.)+fabs((double)chi[0]
-                            +(double)chi[1]+(double)chi[2]+(double)chi[3]
+          double delta_hcp = (fabs((double)chi[0]-3.)
+                              +fabs((double)chi[0]
+                                    +(double)chi[1]
+                                    +(double)chi[2]
+                                    +(double)chi[3]
                                     -9.0))/12.0;
 
           // Identification of the local structure according to the reference
@@ -276,6 +344,7 @@ void ComputeAcklandAtom::compute_peratom()
             }
           }
         }
+      }
     } else structure[i] = 0.0;
   }
 }

src/USER-MISC/compute_ackland_atom.h

@@ -34,7 +34,7 @@ class ComputeAcklandAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax,maxneigh;
+  int nmax,maxneigh,legacy;
   double *distsq;
   int *nearest, *nearest_n0, *nearest_n1;
   double *structure;