From 73ed317254622957db472359d6f63acf9987019e Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Mon, 22 Jun 2020 11:38:02 -0600
Subject: [PATCH 01/18] Added compute_mliap.cpp/h to MLIAP package

---
 src/MLIAP/compute_mliap.cpp | 356 ++++++++++++++++++++++++++++++++++++
 src/MLIAP/compute_mliap.h   |  87 +++++++++
 2 files changed, 443 insertions(+)
 create mode 100644 src/MLIAP/compute_mliap.cpp
 create mode 100644 src/MLIAP/compute_mliap.h

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
new file mode 100644
index 00000000000..80cfbca0588
--- /dev/null
+++ b/src/MLIAP/compute_mliap.cpp
@@ -0,0 +1,356 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <cstring>
+#include <cstdlib>
+#include "mliap_model_linear.h"
+#include "mliap_model_quadratic.h"
+#include "mliap_descriptor_snap.h"
+#include "compute_mliap.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "pair.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{SCALAR,VECTOR,ARRAY};
+
+ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg), cutsq(NULL), list(NULL), mliap(NULL),
+  mliap_peratom(NULL), mliapall(NULL)
+{
+
+  array_flag = 1;
+  extarray = 0;
+
+  int ntypes = atom->ntypes;
+
+  if (narg < 4)
+    error->all(FLERR,"Illegal compute mliap command");
+
+  // process keywords
+
+  int iarg = 0;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"model") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
+      if (strcmp(arg[iarg+1],"linear") == 0) {
+        if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
+        model = new MLIAPModelLinear(lmp,arg[iarg+2]);
+        iarg += 3;
+      } else if (strcmp(arg[iarg+1],"quadratic") == 0) {
+        if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
+        model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
+        iarg += 3;
+      } else error->all(FLERR,"Illegal compute mliap command");
+    } else if (strcmp(arg[iarg],"descriptor") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
+      if (strcmp(arg[iarg+1],"sna") == 0) {
+        if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
+        descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]);
+        iarg += 3;
+      } else error->all(FLERR,"Illegal compute mliap command");
+    }
+  }
+
+  nparams = model->nparams;
+  nperdim = nparams;
+  ndims_force = 3;
+  ndims_virial = 6;
+  yoffset = nperdim;
+  zoffset = 2*nperdim;
+  natoms = atom->natoms;
+  size_array_rows = 1+ndims_force*natoms+ndims_virial;
+  size_array_cols = nperdim*atom->ntypes+1;
+  lastcol = size_array_cols-1;
+
+  ndims_peratom = ndims_force;
+  size_peratom = ndims_peratom*nperdim*atom->ntypes;
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeMLIAP::~ComputeMLIAP()
+{
+  memory->destroy(mliap);
+  memory->destroy(mliapall);
+  memory->destroy(mliap_peratom);
+  memory->destroy(cutsq);
+
+  memory->destroy(map);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeMLIAP::init()
+{
+  if (force->pair == NULL)
+    error->all(FLERR,"Compute mliap requires a pair style be defined");
+
+  if (descriptor->get_cutmax() > force->pair->cutforce)
+    error->all(FLERR,"Compute mliap cutoff is longer than pairwise cutoff");
+
+  // need an occasional full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->occasional = 1;
+
+  int count = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (strcmp(modify->compute[i]->style,"mliap") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute mliap");
+
+  // allocate memory for global array
+
+  memory->create(mliap,size_array_rows,size_array_cols,
+                 "mliap:mliap");
+  memory->create(mliapall,size_array_rows,size_array_cols,
+                 "mliap:mliapall");
+  array = mliapall;
+
+  // find compute for reference energy
+
+  char *id_pe = (char *) "thermo_pe";
+  int ipe = modify->find_compute(id_pe);
+  if (ipe == -1)
+    error->all(FLERR,"compute thermo_pe does not exist.");
+  c_pe = modify->compute[ipe];
+
+  // add compute for reference virial tensor
+
+  char *id_virial = (char *) "mliap_press";
+  char **newarg = new char*[5];
+  newarg[0] = id_virial;
+  newarg[1] = (char *) "all";
+  newarg[2] = (char *) "pressure";
+  newarg[3] = (char *) "NULL";
+  newarg[4] = (char *) "virial";
+  modify->add_compute(5,newarg);
+  delete [] newarg;
+
+  int ivirial = modify->find_compute(id_virial);
+  if (ivirial == -1)
+    error->all(FLERR,"compute mliap_press does not exist.");
+  c_virial = modify->compute[ivirial];
+
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeMLIAP::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeMLIAP::compute_array()
+{
+  int ntotal = atom->nlocal + atom->nghost;
+
+  invoked_array = update->ntimestep;
+
+  // grow mliap_peratom array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->destroy(mliap_peratom);
+    nmax = atom->nmax;
+    memory->create(mliap_peratom,nmax,size_peratom,
+                   "mliap:mliap_peratom");
+  }
+
+  if (gamma_max < list->inum) {
+    memory->grow(descriptors,list->inum,ndescriptors,"PairMLIAP:descriptors");
+    memory->grow(gamma,nparams,list->inum,ndescriptors,"PairMLIAP:gamma");
+    gamma_max = list->inum;
+  }
+
+  // clear global array
+
+  for (int irow = 0; irow < size_array_rows; irow++)
+    for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
+      mliap[irow][icoeff] = 0.0;
+
+  // clear local peratom array
+
+  for (int i = 0; i < ntotal; i++)
+    for (int icoeff = 0; icoeff < size_peratom; icoeff++) {
+      mliap_peratom[i][icoeff] = 0.0;
+    }
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+
+  // compute descriptors, if needed
+
+  if (model->nonlinearflag)
+    descriptor->forward(map, list, descriptors);
+
+  // ***********THIS IS NOT RIGHT**********************
+  // This whole idea is flawed. The gamma matrix is too big to
+  // store. Instead, we should generate the A matrix,
+  // just as ComputeSNAP does, and then pass it to 
+  // the model, which can evaluate gradients of E, F, sigma,
+  // w.r.t. model parameters.
+
+  // calculate descriptor contributions to parameter gradients
+  // and gamma = double gradient w.r.t. parameters and descriptors
+
+  // i.e. gamma = d2E/d\sigma.dB_i
+  // sigma is a parameter and B_i is a descriptor of atom i
+  // for SNAP, this is a sparse nparams*natoms*ndescriptors matrix,
+  // but in general it could be fully dense. 
+ 
+  // *******Not implemented yet*****************
+  // This should populate the energy row and gamma
+  // For the linear model energy row will look just like the Bi accumulation
+  // in ComputeSNAP i.e. accumulating the intput descriptors vector, 
+  // while gamma will be just 1's and 0's
+  // For the quadratic model, the energy row will be similar,
+  // while gamma will be 1's, 0's and Bi's  
+
+  //   model->param_gradient(list, descriptors, mliap[0], gamma);
+
+  // calculate descriptor gradient contributions to parameter gradients
+
+  // *******Not implemented yet*****************
+  // This will just take gamma and multiply it with
+  // descriptor gradient contributions i.e. dblist
+  // this will resemble snadi accumualation in ComputeSNAP
+  // descriptor->param_backward(list, gamma, snadi);
+
+  // accumulate descriptor gradient contributions to global array
+
+  for (int itype = 0; itype < atom->ntypes; itype++) {
+    const int typeoffset_local = ndims_peratom*nperdim*itype;
+    const int typeoffset_global = nperdim*itype;
+    for (int icoeff = 0; icoeff < nperdim; icoeff++) {
+      int irow = 1;
+      for (int i = 0; i < ntotal; i++) {
+        double *snadi = mliap_peratom[i]+typeoffset_local;
+        int iglobal = atom->tag[i];
+        int irow = 3*(iglobal-1)+1;
+        mliap[irow][icoeff+typeoffset_global] += snadi[icoeff];
+        mliap[irow+1][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
+        mliap[irow+2][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
+      }
+    }
+  }
+
+ // accumulate forces to global array
+
+  for (int i = 0; i < atom->nlocal; i++) {
+    int iglobal = atom->tag[i];
+    int irow = 3*(iglobal-1)+1;
+    mliap[irow][lastcol] = atom->f[i][0];
+    mliap[irow+1][lastcol] = atom->f[i][1];
+    mliap[irow+2][lastcol] = atom->f[i][2];
+  }
+
+  // accumulate bispectrum virial contributions to global array
+
+  dbdotr_compute();
+
+  // sum up over all processes
+
+  MPI_Allreduce(&mliap[0][0],&mliapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
+
+  // assign energy to last column
+
+  int irow = 0;
+  double reference_energy = c_pe->compute_scalar();
+  mliapall[irow++][lastcol] = reference_energy;
+
+  // assign virial stress to last column
+  // switch to Voigt notation
+
+  c_virial->compute_vector();
+  irow += 3*natoms;
+  mliapall[irow++][lastcol] = c_virial->vector[0];
+  mliapall[irow++][lastcol] = c_virial->vector[1];
+  mliapall[irow++][lastcol] = c_virial->vector[2];
+  mliapall[irow++][lastcol] = c_virial->vector[5];
+  mliapall[irow++][lastcol] = c_virial->vector[4];
+  mliapall[irow++][lastcol] = c_virial->vector[3];
+
+}
+
+/* ----------------------------------------------------------------------
+   compute global virial contributions via summing r_i.dB^j/dr_i over
+   own & ghost atoms
+------------------------------------------------------------------------- */
+
+void ComputeMLIAP::dbdotr_compute()
+{
+  double **x = atom->x;
+  int irow0 = 1+ndims_force*natoms;
+
+  // sum over bispectrum contributions to forces
+  // on all particles including ghosts
+
+  int nall = atom->nlocal + atom->nghost;
+  for (int i = 0; i < nall; i++)
+    for (int itype = 0; itype < atom->ntypes; itype++) {
+      const int typeoffset_local = ndims_peratom*nperdim*itype;
+      const int typeoffset_global = nperdim*itype;
+      double *snadi = mliap_peratom[i]+typeoffset_local;
+      for (int icoeff = 0; icoeff < nperdim; icoeff++) {
+        double dbdx = snadi[icoeff];
+        double dbdy = snadi[icoeff+yoffset];
+        double dbdz = snadi[icoeff+zoffset];
+        int irow = irow0;
+        mliap[irow++][icoeff+typeoffset_global] += dbdx*x[i][0];
+        mliap[irow++][icoeff+typeoffset_global] += dbdy*x[i][1];
+        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][2];
+        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][1];
+        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][0];
+        mliap[irow++][icoeff+typeoffset_global] += dbdy*x[i][0];
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage
+------------------------------------------------------------------------- */
+
+double ComputeMLIAP::memory_usage()
+{
+
+  double bytes = size_array_rows*size_array_cols *
+    sizeof(double);                                     // mliap
+  bytes += size_array_rows*size_array_cols *
+    sizeof(double);                                     // mliapall
+  bytes += nmax*size_peratom * sizeof(double);          // mliap_peratom
+  int n = atom->ntypes+1;
+  bytes += n*sizeof(int);        // map
+
+  return bytes;
+}
diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h
new file mode 100644
index 00000000000..bc84880b5e2
--- /dev/null
+++ b/src/MLIAP/compute_mliap.h
@@ -0,0 +1,87 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(mliap,ComputeMLIAP)
+
+#else
+
+#ifndef LMP_COMPUTE_MLIAP_H
+#define LMP_COMPUTE_MLIAP_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeMLIAP : public Compute {
+ public:
+  ComputeMLIAP(class LAMMPS *, int, char **);
+  ~ComputeMLIAP();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int natoms, nmax, size_peratom, lastcol;
+  int nperdim, yoffset, zoffset;
+  int ndims_peratom, ndims_force, ndims_virial;
+  double **cutsq;
+  class NeighList *list;
+  double **mliap, **mliapall;
+  double **mliap_peratom;
+  int *map;  // map types to [0,nelements)
+  int nelements;
+
+  double*** gamma;             // gammas for all atoms in list
+  double** descriptors;        // descriptors for all atoms in list
+  int ndescriptors;            // number of descriptors 
+  int gamma_max;               // number of atoms allocated for beta, descriptors
+  int nparams;                 // number of model paramters per element
+
+  class MLIAPModel* model;
+  class MLIAPDescriptor* descriptor;
+
+  Compute *c_pe;
+  Compute *c_virial;
+
+  void dbdotr_compute();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute snap requires a pair style be defined
+
+Self-explanatory.
+
+E: Compute snap cutoff is longer than pairwise cutoff
+
+UNDOCUMENTED
+
+W: More than one compute snad/atom
+
+Self-explanatory.
+
+*/

From 9c7005a91cdda15d2142cec1d328a856c7fc4f43 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Tue, 30 Jun 2020 19:14:01 -0600
Subject: [PATCH 02/18] It compiles, but not yet working

---
 src/MLIAP/compute_mliap.cpp         |  60 ++++++++----
 src/MLIAP/compute_mliap.h           |   7 +-
 src/MLIAP/mliap_descriptor.h        |   7 +-
 src/MLIAP/mliap_descriptor_snap.cpp | 145 ++++++++++++++++++++++------
 src/MLIAP/mliap_descriptor_snap.h   |   4 +-
 src/MLIAP/mliap_model.h             |   2 +
 src/MLIAP/mliap_model_linear.cpp    |  64 +++++++++++-
 src/MLIAP/mliap_model_linear.h      |   2 +
 src/MLIAP/mliap_model_quadratic.cpp | 107 ++++++++++++++++++++
 src/MLIAP/mliap_model_quadratic.h   |   2 +
 src/MLIAP/pair_mliap.cpp            |   5 +-
 11 files changed, 348 insertions(+), 57 deletions(-)

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index 80cfbca0588..e86dce732ef 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
 enum{SCALAR,VECTOR,ARRAY};
 
 ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), cutsq(NULL), list(NULL), mliap(NULL),
+  Compute(lmp, narg, arg), list(NULL), mliap(NULL),
   mliap_peratom(NULL), mliapall(NULL)
 {
 
@@ -46,6 +46,11 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 4)
     error->all(FLERR,"Illegal compute mliap command");
 
+  // set flags for required keywords
+
+  int modelflag = 0;
+  int descriptorflag = 0;
+
   // process keywords
 
   int iarg = 0;
@@ -62,6 +67,7 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
         model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
         iarg += 3;
       } else error->all(FLERR,"Illegal compute mliap command");
+      modelflag = 1;
     } else if (strcmp(arg[iarg],"descriptor") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
       if (strcmp(arg[iarg+1],"sna") == 0) {
@@ -69,10 +75,17 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
         descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]);
         iarg += 3;
       } else error->all(FLERR,"Illegal compute mliap command");
-    }
+      descriptorflag = 1;
+    } else
+      error->all(FLERR,"Illegal compute mliap command");
   }
 
+  if (modelflag == 0 || descriptorflag == 0)
+    error->all(FLERR,"Illegal compute_style command");
+
   nparams = model->nparams;
+  nelements = model->nelements;
+  gamma_nnz = model->get_gamma_nnz();
   nperdim = nparams;
   ndims_force = 3;
   ndims_virial = 6;
@@ -96,9 +109,17 @@ ComputeMLIAP::~ComputeMLIAP()
   memory->destroy(mliap);
   memory->destroy(mliapall);
   memory->destroy(mliap_peratom);
-  memory->destroy(cutsq);
 
   memory->destroy(map);
+
+  memory->destroy(descriptors);
+  memory->destroy(gamma_row_index);
+  memory->destroy(gamma_col_index);
+  memory->destroy(gamma);
+  memory->destroy(egradient);
+
+  delete model;
+  delete descriptor;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -108,7 +129,7 @@ void ComputeMLIAP::init()
   if (force->pair == NULL)
     error->all(FLERR,"Compute mliap requires a pair style be defined");
 
-  if (descriptor->get_cutmax() > force->pair->cutforce)
+  if (descriptor->cutmax > force->pair->cutforce)
     error->all(FLERR,"Compute mliap cutoff is longer than pairwise cutoff");
 
   // need an occasional full neighbor list
@@ -187,8 +208,11 @@ void ComputeMLIAP::compute_array()
   }
 
   if (gamma_max < list->inum) {
-    memory->grow(descriptors,list->inum,ndescriptors,"PairMLIAP:descriptors");
-    memory->grow(gamma,nparams,list->inum,ndescriptors,"PairMLIAP:gamma");
+    memory->grow(descriptors,list->inum,ndescriptors,"ComputeMLIAP:descriptors");
+    memory->grow(gamma_row_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_row_index");
+    memory->grow(gamma_col_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_col_index");
+    memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
+    memory->grow(egradient,nelements*nparams,"ComputeMLIAP:egradient");
     gamma_max = list->inum;
   }
 
@@ -214,19 +238,12 @@ void ComputeMLIAP::compute_array()
   if (model->nonlinearflag)
     descriptor->forward(map, list, descriptors);
 
-  // ***********THIS IS NOT RIGHT**********************
-  // This whole idea is flawed. The gamma matrix is too big to
-  // store. Instead, we should generate the A matrix,
-  // just as ComputeSNAP does, and then pass it to 
-  // the model, which can evaluate gradients of E, F, sigma,
-  // w.r.t. model parameters.
-
   // calculate descriptor contributions to parameter gradients
   // and gamma = double gradient w.r.t. parameters and descriptors
 
-  // i.e. gamma = d2E/d\sigma.dB_i
-  // sigma is a parameter and B_i is a descriptor of atom i
-  // for SNAP, this is a sparse nparams*natoms*ndescriptors matrix,
+  // i.e. gamma = d2E/dsigma_l.dB_k
+  // sigma_l is a parameter and B_k is a descriptor of atom i
+  // for SNAP, this is a sparse natoms*nparams*ndescriptors matrix,
   // but in general it could be fully dense. 
  
   // *******Not implemented yet*****************
@@ -237,16 +254,23 @@ void ComputeMLIAP::compute_array()
   // For the quadratic model, the energy row will be similar,
   // while gamma will be 1's, 0's and Bi's  
 
-  //   model->param_gradient(list, descriptors, mliap[0], gamma);
+  model->param_gradient(map, list, descriptors, gamma_row_index, 
+                        gamma_col_index, gamma, egradient);
+
 
   // calculate descriptor gradient contributions to parameter gradients
 
   // *******Not implemented yet*****************
   // This will just take gamma and multiply it with
   // descriptor gradient contributions i.e. dblist
-  // this will resemble snadi accumualation in ComputeSNAP
+  // this will resemble snadi accumulation in ComputeSNAP
+
   // descriptor->param_backward(list, gamma, snadi);
 
+  descriptor->param_backward(map, list, gamma_nnz, gamma_row_index, 
+                             gamma_col_index, gamma, mliap_peratom,
+                             yoffset, zoffset);
+
   // accumulate descriptor gradient contributions to global array
 
   for (int itype = 0; itype < atom->ntypes; itype++) {
diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h
index bc84880b5e2..14f4a4bdaed 100644
--- a/src/MLIAP/compute_mliap.h
+++ b/src/MLIAP/compute_mliap.h
@@ -37,18 +37,21 @@ class ComputeMLIAP : public Compute {
   int natoms, nmax, size_peratom, lastcol;
   int nperdim, yoffset, zoffset;
   int ndims_peratom, ndims_force, ndims_virial;
-  double **cutsq;
   class NeighList *list;
   double **mliap, **mliapall;
   double **mliap_peratom;
   int *map;  // map types to [0,nelements)
   int nelements;
 
-  double*** gamma;             // gammas for all atoms in list
   double** descriptors;        // descriptors for all atoms in list
   int ndescriptors;            // number of descriptors 
   int gamma_max;               // number of atoms allocated for beta, descriptors
   int nparams;                 // number of model paramters per element
+  int gamma_nnz;               // number of non-zero entries in gamma
+  double** gamma;              // gamma element
+  int** gamma_row_index;       // row (parameter) index 
+  int** gamma_col_index;       // column (descriptor) index 
+  double* egradient;           // energy gradient w.r.t. parameters
 
   class MLIAPModel* model;
   class MLIAPDescriptor* descriptor;
diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h
index cd42cb3be6b..f7f0db250db 100644
--- a/src/MLIAP/mliap_descriptor.h
+++ b/src/MLIAP/mliap_descriptor.h
@@ -24,15 +24,16 @@ class MLIAPDescriptor : protected Pointers  {
   ~MLIAPDescriptor();
   virtual void forward(int*, class NeighList*, double**)=0;
   virtual void backward(class PairMLIAP*, class NeighList*, double**, int)=0;
+  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 
+                              double**, int, int)=0;
   virtual void init()=0;
-  virtual double get_cutoff(int, int)=0;
-  virtual double get_cutmax()=0;
   virtual double memory_usage()=0;
 
   int ndescriptors;              // number of descriptors
   int nelements;                 // # of unique elements
   char **elements;               // names of unique elements
-
+  double **cutsq;                // nelem x nelem rcutsq values 
+  double cutmax;                 // maximum cutoff needed
 protected:
 
 };
diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp
index eedeeb44835..de3d7d3c9b6 100644
--- a/src/MLIAP/mliap_descriptor_snap.cpp
+++ b/src/MLIAP/mliap_descriptor_snap.cpp
@@ -58,6 +58,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()
     delete[] elements;
     memory->destroy(radelem);
     memory->destroy(wjelem);
+    memory->destroy(cutsq);
   }
 
   delete snaptr;
@@ -111,16 +112,13 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
       int jtype = type[j];
       int jelem = map[jtype];
 
-      //      printf("i = %d j = %d itype = %d jtype = %d cutsq[i][j] = %g rsq = %g\n",i,j,itype,jtype,cutsq[itype][jtype],rsq);
-
-      double rcutsqtmp = get_cutoff(ielem, jelem);
-      if (rsq < rcutsqtmp*rcutsqtmp) {
+      if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
         snaptr->wj[ninside] = wjelem[jelem];
-        snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
+        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
         snaptr->element[ninside] = jelem; // element index for chem snap
         ninside++;
       }
@@ -196,13 +194,13 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
       int jtype = type[j];
       int jelem = pairmliap->map[jtype];
 
-      if (rsq < pairmliap->cutsq[itype][jtype]&&rsq>1e-20) {
+      if (rsq < cutsq[itype][jtype]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
         snaptr->wj[ninside] = wjelem[jelem];
-        snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
+        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
         snaptr->element[ninside] = jelem; // element index for chem snap
         ninside++;
       }
@@ -242,7 +240,7 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
       f[j][1] -= fij[1];
       f[j][2] -= fij[2];
 
-      // add in gloabl and per-atom virial contributions
+      // add in global and per-atom virial contributions
       // this is optional and has no effect on force calculation
       
       if (vflag)
@@ -256,6 +254,105 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
 
 }
 
+/* ----------------------------------------------------------------------
+   compute forces for each atom
+   ---------------------------------------------------------------------- */
+
+void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, 
+                                         int gamma_nnz, int **gamma_row_index, 
+                                         int **gamma_col_index, double **gamma, double **snadi,
+                                         int yoffset, int zoffset)
+{
+  int i,j,jnum,ninside;
+  double delx,dely,delz,evdwl,rsq;
+  double fij[3];
+  int *jlist,*numneigh,**firstneigh;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  for (int ii = 0; ii < list->inum; ii++) {
+    i = list->ilist[ii];
+
+    const double xtmp = x[i][0];
+    const double ytmp = x[i][1];
+    const double ztmp = x[i][2];
+    const int itype = type[i];
+    int ielem = 0;
+    if (chemflag)
+      ielem = map[itype];
+    const double radi = radelem[ielem];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    // insure rij, inside, wj, and rcutij are of size jnum
+
+    snaptr->grow_rij(jnum);
+
+    // rij[][3] = displacements between atom I and those neighbors
+    // inside = indices of neighbors of I within cutoff
+    // wj = weights for neighbors of I within cutoff
+    // rcutij = cutoffs for neighbors of I within cutoff
+    // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
+
+    ninside = 0;
+    for (int jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      delx = x[j][0] - xtmp;
+      dely = x[j][1] - ytmp;
+      delz = x[j][2] - ztmp;
+      rsq = delx*delx + dely*dely + delz*delz;
+      int jtype = type[j];
+      int jelem = map[jtype];
+
+      if (rsq < cutsq[itype][jtype]) {
+        snaptr->rij[ninside][0] = delx;
+        snaptr->rij[ninside][1] = dely;
+        snaptr->rij[ninside][2] = delz;
+        snaptr->inside[ninside] = j;
+        snaptr->wj[ninside] = wjelem[jelem];
+        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
+        snaptr->element[ninside] = jelem; // element index for chem snap
+        ninside++;
+      }
+    }
+
+    snaptr->compute_ui(ninside, ielem);
+    snaptr->compute_zi();
+    snaptr->compute_bi(ielem);
+
+    for (int jj = 0; jj < ninside; jj++) {
+      const int j = snaptr->inside[jj];
+      snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                             snaptr->rcutij[jj], jj, snaptr->element[jj]);
+      snaptr->compute_dbidrj();
+
+      // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj
+
+        for (int inz = 0; inz < gamma_nnz; inz++) {
+          const int l = gamma_row_index[ii][inz];
+          const int k = gamma_col_index[ii][inz];
+          snadi[i][l]         += gamma[ii][inz]*snaptr->dblist[k][0];
+          snadi[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1];
+          snadi[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2];
+          snadi[j][l]         -= gamma[ii][inz]*snaptr->dblist[k][0];
+          snadi[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1];
+          snadi[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2];
+        }
+
+    }
+  }
+
+}
+
 /* ----------------------------------------------------------------------
    set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */
@@ -270,7 +367,6 @@ void MLIAPDescriptorSNAP::init()
   snaptr->init();
 
   ndescriptors = snaptr->ncoeff;
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -417,31 +513,22 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
       !elementsflag || !radelemflag || !wjelemflag)
     error->all(FLERR,"Incorrect SNAP parameter file");
 
-}
+  // construct cutsq
 
-/* ----------------------------------------------------------------------
-   provide cutoff distance for two elements
-------------------------------------------------------------------------- */
-
-double MLIAPDescriptorSNAP::get_cutoff(int ielem, int jelem)
-{
-  return (radelem[ielem] + radelem[jelem])*rcutfac;
-}
-
-/* ----------------------------------------------------------------------
-   calculate maximum cutoff distance
-------------------------------------------------------------------------- */
-
-double MLIAPDescriptorSNAP::get_cutmax()
-{
   double cut;
-  double cutmax = 0.0;
-  for(int ielem = 0; ielem <= nelements; ielem++) {
+  cutmax = 0.0;
+  memory->create(cutsq,nelements,nelements,"mliap/descriptor/snap:cutsq");
+  for (int ielem = 0; ielem < nelements; ielem++) {
     cut = 2.0*radelem[ielem]*rcutfac;
     if (cut > cutmax) cutmax = cut;
-    return cutmax;
+    cutsq[ielem][ielem] = cut*cut;
+    printf("ielem = %d ielem = %d cutsq[i][i] = %g\n",ielem, ielem, cutsq[ielem][ielem]);
+    for(int jelem = ielem+1; jelem < nelements; jelem++) {
+      cut = (radelem[ielem]+radelem[jelem])*rcutfac;
+      cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut*cut;
+      printf("ielem = %d jelem = %d cutsq[i][j] = %g\n",ielem, jelem, cutsq[ielem][jelem]);
+    }
   }
-  return cutmax;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h
index 15691fabfef..4843fbbcc80 100644
--- a/src/MLIAP/mliap_descriptor_snap.h
+++ b/src/MLIAP/mliap_descriptor_snap.h
@@ -24,9 +24,9 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor  {
   ~MLIAPDescriptorSNAP();
   virtual void forward(int*, class NeighList*, double**);
   virtual void backward(class PairMLIAP*, class NeighList*, double**, int);
+  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 
+                              double**, int, int);
   virtual void init();
-  virtual double get_cutoff(int, int);
-  virtual double get_cutmax();
   virtual double memory_usage();
 
   double rcutfac;                // declared public to workaround gcc 4.9
diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h
index aeb16cb299a..06b23707a98 100644
--- a/src/MLIAP/mliap_model.h
+++ b/src/MLIAP/mliap_model.h
@@ -23,6 +23,8 @@ class MLIAPModel : protected Pointers {
   MLIAPModel(LAMMPS*, char*);
   ~MLIAPModel();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int)=0;
+  virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*)=0;
+  virtual int get_gamma_nnz()=0;
   virtual void init();
   virtual double memory_usage();
   int nelements;                 // # of unique elements
diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp
index fe3665938fd..63e8220d141 100644
--- a/src/MLIAP/mliap_model_linear.cpp
+++ b/src/MLIAP/mliap_model_linear.cpp
@@ -40,7 +40,8 @@ MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
 MLIAPModelLinear::~MLIAPModelLinear(){}
 
 /* ----------------------------------------------------------------------
-   Calculate model gradients w.r.t descriptors for each atom dE(B_i)/dB_i
+   Calculate model gradients w.r.t descriptors 
+   for each atom beta_i = dE(B_i)/dB_i
    ---------------------------------------------------------------------- */
 
 void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double **descriptors, double **beta, int eflag)
@@ -77,3 +78,64 @@ void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double **
   }
 }
 
+/* ----------------------------------------------------------------------
+   Calculate model double gradients w.r.t descriptors and parameters
+   for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, 
+   where sigma_l is a parameter, B_k a descriptor, 
+   and atom subscript i is omitted
+
+   gamma is in CSR format:
+      nnz = number of non-zero values
+      gamma_row_index[inz] = l indices, 0 <= l < nparams 
+      gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors
+      gamma[i][inz] = non-zero values, 0 <= inz < nnz
+
+   egradient is derivative of energy w.r.t. parameters
+   ---------------------------------------------------------------------- */
+
+void MLIAPModelLinear::param_gradient(int *map, NeighList* list, 
+                                         double **descriptors, 
+                                         int **gamma_row_index, int **gamma_col_index, 
+                                         double **gamma, double *egradient)
+{
+  int i;
+  int *type = atom->type;
+
+  // zero out energy gradients
+
+  for (int l = 0; l < nelements*nparams; l++)
+    egradient[l] = 0.0;
+    
+  for (int ii = 0; ii < list->inum; ii++) {
+
+    i = list->ilist[ii];
+    const int itype = type[i];
+    const int ielem = map[itype];
+    const int elemoffset = nparams*ielem;
+
+    int l = elemoffset+1;
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
+      gamma[ii][icoeff] = 1.0;
+      gamma_row_index[ii][icoeff] = l++;
+      gamma_col_index[ii][icoeff] = icoeff;
+    }
+
+    // gradient of energy of atom I w.r.t. parameters
+    
+    l = elemoffset;
+    egradient[l++] += 1.0;
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++)
+      egradient[l++] += descriptors[ii][icoeff];
+    
+  }
+
+}
+/* ----------------------------------------------------------------------
+   count the number of non-zero entries in gamma matrix
+   ---------------------------------------------------------------------- */
+
+int MLIAPModelLinear::get_gamma_nnz()
+{
+  int inz = ndescriptors;
+  return inz;
+}
diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h
index 53c7c36cefc..b7646aab458 100644
--- a/src/MLIAP/mliap_model_linear.h
+++ b/src/MLIAP/mliap_model_linear.h
@@ -23,6 +23,8 @@ class MLIAPModelLinear : public MLIAPModel {
   MLIAPModelLinear(LAMMPS*, char*);
   ~MLIAPModelLinear();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
+  virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);
+  virtual int get_gamma_nnz();
 
 protected:
 };
diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp
index 94d89724167..554d1a06a76 100644
--- a/src/MLIAP/mliap_model_quadratic.cpp
+++ b/src/MLIAP/mliap_model_quadratic.cpp
@@ -101,3 +101,110 @@ void MLIAPModelQuadratic::gradient(PairMLIAP* pairmliap, NeighList* list, double
   }
 }
 
+
+/* ----------------------------------------------------------------------
+   Calculate model double gradients w.r.t descriptors and parameters
+   for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, 
+   where sigma_l is a parameter, B_k a descriptor, 
+   and atom subscript i is omitted
+
+   gamma is in CSR format:
+      nnz = number of non-zero values
+      gamma_row_index[inz] = l indices, 0 <= l < nparams 
+      gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors
+      gamma[i][inz] = non-zero values, 0 <= inz < nnz
+
+   egradient is derivative of energy w.r.t. parameters
+   ---------------------------------------------------------------------- */
+
+void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list, 
+                                         double **descriptors, 
+                                         int **gamma_row_index, int **gamma_col_index, 
+                                         double **gamma, double *egradient)
+{
+  int i;
+  int *type = atom->type;
+
+  // zero out energy gradients
+
+  for (int l = 0; l < nelements*nparams; l++)
+    egradient[l] = 0.0;
+    
+  for (int ii = 0; ii < list->inum; ii++) {
+
+    i = list->ilist[ii];
+    const int itype = type[i];
+    const int ielem = map[itype];
+    const int elemoffset = nparams*ielem;
+
+    // linear contributions
+
+    int l = elemoffset+1;
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
+      gamma[ii][icoeff] = 1.0;
+      gamma_row_index[ii][icoeff] = l++;
+      gamma_col_index[ii][icoeff] = icoeff;
+    }
+
+    // quadratic contributions
+
+    int inz = ndescriptors;
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
+      double bveci = descriptors[ii][icoeff];
+      gamma[ii][inz] = bveci;
+      gamma_row_index[ii][inz] = l++;
+      gamma_col_index[ii][inz] = icoeff;
+      inz++;
+      for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
+        double bvecj = descriptors[ii][jcoeff];
+        gamma[ii][inz] = bvecj; // derivative w.r.t. B[icoeff]
+        gamma_row_index[ii][inz] = l;
+        gamma_col_index[ii][inz] = icoeff;
+        inz++;
+        gamma[ii][inz] = bveci; // derivative w.r.t. B[jcoeff]
+        gamma_row_index[ii][inz] = l;
+        gamma_col_index[ii][inz] = jcoeff;
+        inz++;
+        l++;
+      }
+    }
+
+    // gradient of energy of atom I w.r.t. parameters
+    
+    l = elemoffset;
+    egradient[l++] += 1.0;
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++)
+      egradient[l++] += descriptors[ii][icoeff];
+    
+    // quadratic contributions
+    
+    for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
+      double bveci = descriptors[ii][icoeff];
+      egradient[l++] += 0.5*bveci*bveci;
+      for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
+        double bvecj = descriptors[ii][jcoeff];
+        egradient[l++] += bveci*bvecj;
+      }
+    }
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+   count the number of non-zero entries in gamma matrix
+   ---------------------------------------------------------------------- */
+
+int MLIAPModelQuadratic::get_gamma_nnz()
+{
+  int inz = ndescriptors;
+  for (int icoeff = 0; icoeff < ndescriptors; icoeff++) {
+    inz++;
+    for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) {
+        inz++;
+        inz++;
+    }
+  }
+
+  return inz;
+}
+
diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h
index 6ca0697919c..4cfc669bbb0 100644
--- a/src/MLIAP/mliap_model_quadratic.h
+++ b/src/MLIAP/mliap_model_quadratic.h
@@ -23,6 +23,8 @@ class MLIAPModelQuadratic : public MLIAPModel {
   MLIAPModelQuadratic(LAMMPS*, char*);
   ~MLIAPModelQuadratic();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
+  virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);
+  virtual int get_gamma_nnz();
 
 protected:
 };
diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp
index b97dd51d67d..129fcdd6b78 100644
--- a/src/MLIAP/pair_mliap.cpp
+++ b/src/MLIAP/pair_mliap.cpp
@@ -310,7 +310,9 @@ void PairMLIAP::init_style()
 double PairMLIAP::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
-  return descriptor->get_cutoff(map[i],map[j]);
+  printf("itype = %d jtype = %d map[i] = %d map[j] = %d cutsq = %g\n",
+         i,j,map[i],map[j],descriptor->cutsq[map[i]][map[j]]);
+  return sqrt(descriptor->cutsq[map[i]][map[j]]);
 }
 
 /* ----------------------------------------------------------------------
@@ -323,7 +325,6 @@ double PairMLIAP::memory_usage()
 
   int n = atom->ntypes+1;
   bytes += n*n*sizeof(int);      // setflag
-  bytes += n*n*sizeof(double);   // cutsq
   bytes += beta_max*ndescriptors*sizeof(double); // descriptors
   bytes += beta_max*ndescriptors*sizeof(double); // beta
 

From aee3d9f2cdde443f4d1871c408486fad81625eae Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Wed, 1 Jul 2020 15:03:45 -0600
Subject: [PATCH 03/18] Generated working example, not quite correct

---
 examples/mliap/in.mliap.snap.compute | 95 ++++++++++++++++++++++++++++
 src/MLIAP/compute_mliap.cpp          | 66 +++++++++++++------
 src/MLIAP/mliap_descriptor_snap.cpp  | 67 +++++++++++++-------
 src/MLIAP/mliap_model.cpp            | 10 +++
 src/MLIAP/mliap_model.h              |  1 +
 src/MLIAP/mliap_model_linear.cpp     | 13 +++-
 src/MLIAP/mliap_model_linear.h       |  1 +
 src/MLIAP/mliap_model_quadratic.cpp  |  9 +++
 src/MLIAP/mliap_model_quadratic.h    |  1 +
 9 files changed, 217 insertions(+), 46 deletions(-)
 create mode 100644 examples/mliap/in.mliap.snap.compute

diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute
new file mode 100644
index 00000000000..813e568387c
--- /dev/null
+++ b/examples/mliap/in.mliap.snap.compute
@@ -0,0 +1,95 @@
+# Demonstrate bispectrum computes
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	ny equal ${nrep}
+variable 	nz equal ${nrep}
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+create_box	2 box
+create_atoms	2 box
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_coeff 	* * ${zblz} ${zblz}
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 0
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string &
+"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	vb all snav/atom ${snap_options}
+compute 	db all snad/atom ${snap_options}
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+group 		snapgroup2 type 2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_30 equal c_db[2][30]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(B_{000}^i, all i of type 2) 
+#                4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom &
+		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 &
+		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] 
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index e86dce732ef..69c00fc94c5 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -35,14 +35,13 @@ enum{SCALAR,VECTOR,ARRAY};
 
 ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg), list(NULL), mliap(NULL),
-  mliap_peratom(NULL), mliapall(NULL)
+  mliap_peratom(NULL), mliapall(NULL), map(NULL), 
+  descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL),
+  gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL)
 {
-
   array_flag = 1;
   extarray = 0;
 
-  int ntypes = atom->ntypes;
-
   if (narg < 4)
     error->all(FLERR,"Illegal compute mliap command");
 
@@ -53,19 +52,23 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
 
   // process keywords
 
-  int iarg = 0;
+  int iarg = 3;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"model") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");
       if (strcmp(arg[iarg+1],"linear") == 0) {
-        if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
-        model = new MLIAPModelLinear(lmp,arg[iarg+2]);
-        iarg += 3;
+	if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command");
+	int ntmp1 = atoi(arg[iarg+2]);
+	int ntmp2 = atoi(arg[iarg+3]);
+        model = new MLIAPModelLinear(lmp,ntmp1,ntmp2);
+        iarg += 4;
       } else if (strcmp(arg[iarg+1],"quadratic") == 0) {
-        if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command");
-        model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);
-        iarg += 3;
+	if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command");
+	int ntmp1 = atoi(arg[iarg+2]);
+	int ntmp2 = atoi(arg[iarg+3]);
+        model = new MLIAPModelQuadratic(lmp,ntmp1,ntmp2);
+        iarg += 4;
       } else error->all(FLERR,"Illegal compute mliap command");
       modelflag = 1;
     } else if (strcmp(arg[iarg],"descriptor") == 0) {
@@ -83,6 +86,7 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   if (modelflag == 0 || descriptorflag == 0)
     error->all(FLERR,"Illegal compute_style command");
 
+  ndescriptors = descriptor->ndescriptors;
   nparams = model->nparams;
   nelements = model->nelements;
   gamma_nnz = model->get_gamma_nnz();
@@ -100,6 +104,8 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   size_peratom = ndims_peratom*nperdim*atom->ntypes;
 
   nmax = 0;
+  gamma_max = 0;
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -147,14 +153,32 @@ void ComputeMLIAP::init()
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one compute mliap");
 
+  // initialize model and descriptor
+
+  model->init();
+  descriptor->init();
+
+  // consistency checks
+
+  if (descriptor->ndescriptors != model->ndescriptors)
+    error->all(FLERR,"Incompatible model and descriptor definitions");
+  if (descriptor->nelements != model->nelements)
+    error->all(FLERR,"Incompatible model and descriptor definitions");
+  if (nelements != atom->ntypes)
+    error->all(FLERR,"nelements must equal ntypes");
+
   // allocate memory for global array
 
+  printf("size_array_rows = %d size_array_cols = %d\n",size_array_rows,size_array_cols);
   memory->create(mliap,size_array_rows,size_array_cols,
                  "mliap:mliap");
   memory->create(mliapall,size_array_rows,size_array_cols,
                  "mliap:mliapall");
   array = mliapall;
 
+  printf("nelements = %d nparams = %d\n",nelements,nparams);
+  memory->create(egradient,nelements*nparams,"ComputeMLIAP:egradient");
+
   // find compute for reference energy
 
   char *id_pe = (char *) "thermo_pe";
@@ -203,19 +227,11 @@ void ComputeMLIAP::compute_array()
   if (atom->nmax > nmax) {
     memory->destroy(mliap_peratom);
     nmax = atom->nmax;
+    printf("nmax = %d size_peratom = %d\n",nmax,size_peratom);
     memory->create(mliap_peratom,nmax,size_peratom,
                    "mliap:mliap_peratom");
   }
 
-  if (gamma_max < list->inum) {
-    memory->grow(descriptors,list->inum,ndescriptors,"ComputeMLIAP:descriptors");
-    memory->grow(gamma_row_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_row_index");
-    memory->grow(gamma_col_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_col_index");
-    memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
-    memory->grow(egradient,nelements*nparams,"ComputeMLIAP:egradient");
-    gamma_max = list->inum;
-  }
-
   // clear global array
 
   for (int irow = 0; irow < size_array_rows; irow++)
@@ -233,6 +249,15 @@ void ComputeMLIAP::compute_array()
 
   neighbor->build_one(list);
 
+  if (gamma_max < list->inum) {
+    memory->grow(descriptors,list->inum,ndescriptors,"ComputeMLIAP:descriptors");
+    memory->grow(gamma_row_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_row_index");
+    memory->grow(gamma_col_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_col_index");
+    memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
+    gamma_max = list->inum;
+  }
+  printf("gamma_max %d %d %d %d %d %d %p\n",gamma_max, ndescriptors, gamma_nnz, nelements, nparams, list->inum, list);
+
   // compute descriptors, if needed
 
   if (model->nonlinearflag)
@@ -312,6 +337,7 @@ void ComputeMLIAP::compute_array()
   int irow = 0;
   double reference_energy = c_pe->compute_scalar();
   mliapall[irow++][lastcol] = reference_energy;
+  printf("Reference energy = %g %g %g %d %d\n",reference_energy,mliapall[irow-1][lastcol],array[irow-1][lastcol],irow-1,lastcol);
 
   // assign virial stress to last column
   // switch to Voigt notation
diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp
index de3d7d3c9b6..ffb4ab65ff5 100644
--- a/src/MLIAP/mliap_descriptor_snap.cpp
+++ b/src/MLIAP/mliap_descriptor_snap.cpp
@@ -45,6 +45,13 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename):
   wjelem = NULL;
   snaptr = NULL;
   read_paramfile(paramfilename);
+
+  snaptr = new SNA(lmp, rfac0, twojmax,
+                   rmin0, switchflag, bzeroflag,
+                   chemflag, bnormflag, wselfallflag, nelements);
+
+  ndescriptors = snaptr->ncoeff;
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -85,8 +92,9 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
     const double ytmp = x[i][1];
     const double ztmp = x[i][2];
     const int itype = type[i];
-    const int ielem = map[itype];
-    const double radi = radelem[ielem];
+    // element map not yet implemented
+    // const int ielem = map[itype];
+    const int ielem = itype-1;
 
     jlist = list->firstneigh[i];
     jnum = list->numneigh[i];
@@ -110,16 +118,18 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
       delz = x[j][2] - ztmp;
       rsq = delx*delx + dely*dely + delz*delz;
       int jtype = type[j];
-      int jelem = map[jtype];
+      // element map not yet implemented
+      // const int jelem = map[jtype];
+      const int jelem = jtype-1;
 
       if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
-        snaptr->wj[ninside] = wjelem[jelem];
-        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
-        snaptr->element[ninside] = jelem; // element index for chem snap
+	snaptr->wj[ninside] = wjelem[jelem];
+	snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
+	snaptr->element[ninside] = jelem; // element index for chem snap
         ninside++;
       }
     }
@@ -129,6 +139,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
     else
       snaptr->compute_ui(ninside, 0);
     snaptr->compute_zi();
+
     if (chemflag)
       snaptr->compute_bi(ielem);
     else
@@ -168,7 +179,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
     const double ztmp = x[i][2];
     const int itype = type[i];
     const int ielem = pairmliap->map[itype];
-    const double radi = radelem[ielem];
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -194,7 +204,7 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
       int jtype = type[j];
       int jelem = pairmliap->map[jtype];
 
-      if (rsq < cutsq[itype][jtype]) {
+      if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
@@ -286,8 +296,9 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
     const int itype = type[i];
     int ielem = 0;
     if (chemflag)
-      ielem = map[itype];
-    const double radi = radelem[ielem];
+    // element map not yet implemented
+    //   ielem = map[itype];
+      const int ielem = itype-1;
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -311,28 +322,43 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
       delz = x[j][2] - ztmp;
       rsq = delx*delx + dely*dely + delz*delz;
       int jtype = type[j];
-      int jelem = map[jtype];
-
-      if (rsq < cutsq[itype][jtype]) {
+      // element map not yet implemented
+      // int jelem = map[jtype];
+      const int jelem = jtype-1;
+      if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
-        snaptr->wj[ninside] = wjelem[jelem];
-        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
+	snaptr->wj[ninside] = wjelem[jelem];
+	snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
         snaptr->element[ninside] = jelem; // element index for chem snap
         ninside++;
       }
     }
 
-    snaptr->compute_ui(ninside, ielem);
+    if(chemflag)
+      snaptr->compute_ui(ninside, ielem);
+    else
+      snaptr->compute_ui(ninside, 0);
+
+
     snaptr->compute_zi();
-    snaptr->compute_bi(ielem);
+    if(chemflag)
+      snaptr->compute_bi(ielem);
+    else
+      snaptr->compute_bi(0);
 
     for (int jj = 0; jj < ninside; jj++) {
       const int j = snaptr->inside[jj];
+
+    if(chemflag)
       snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
                              snaptr->rcutij[jj], jj, snaptr->element[jj]);
+    else
+      snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                             snaptr->rcutij[jj], jj, 0);
+
       snaptr->compute_dbidrj();
 
       // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj
@@ -359,14 +385,7 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
 
 void MLIAPDescriptorSNAP::init()
 {
-
-  snaptr = new SNA(lmp, rfac0, twojmax,
-                   rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
-
   snaptr->init();
-
-  ndescriptors = snaptr->ncoeff;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/MLIAP/mliap_model.cpp b/src/MLIAP/mliap_model.cpp
index 94aafc93a83..5ab96e19d4b 100644
--- a/src/MLIAP/mliap_model.cpp
+++ b/src/MLIAP/mliap_model.cpp
@@ -40,6 +40,16 @@ MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp)
 
 /* ---------------------------------------------------------------------- */
 
+MLIAPModel::MLIAPModel(LAMMPS* lmp, int nelements_in, int nparams_in) : Pointers(lmp)
+{
+  nelements = nelements_in;
+  nparams = nparams_in;
+  coeffelem = NULL;
+  nonlinearflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
 MLIAPModel::~MLIAPModel()
 {
   memory->destroy(coeffelem);
diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h
index 06b23707a98..57f6f762149 100644
--- a/src/MLIAP/mliap_model.h
+++ b/src/MLIAP/mliap_model.h
@@ -21,6 +21,7 @@ namespace LAMMPS_NS {
 class MLIAPModel : protected Pointers {
 public:
   MLIAPModel(LAMMPS*, char*);
+  MLIAPModel(LAMMPS*, int, int);
   ~MLIAPModel();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int)=0;
   virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*)=0;
diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp
index 63e8220d141..a6d098ed265 100644
--- a/src/MLIAP/mliap_model_linear.cpp
+++ b/src/MLIAP/mliap_model_linear.cpp
@@ -31,7 +31,14 @@ using namespace LAMMPS_NS;
 MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : 
   MLIAPModel(lmp, coefffilename)
 {
-  nonlinearflag = 0;
+  ndescriptors = nparams - 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) : 
+  MLIAPModel(lmp, nelements_in, nparams_in)
+{
   ndescriptors = nparams - 1;
 }
 
@@ -110,7 +117,9 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list,
 
     i = list->ilist[ii];
     const int itype = type[i];
-    const int ielem = map[itype];
+    // element map not yet implemented
+    // const int ielem = map[itype];
+    const int ielem = itype-1;
     const int elemoffset = nparams*ielem;
 
     int l = elemoffset+1;
diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h
index b7646aab458..db4cb265146 100644
--- a/src/MLIAP/mliap_model_linear.h
+++ b/src/MLIAP/mliap_model_linear.h
@@ -21,6 +21,7 @@ namespace LAMMPS_NS {
 class MLIAPModelLinear : public MLIAPModel {
 public:
   MLIAPModelLinear(LAMMPS*, char*);
+  MLIAPModelLinear(LAMMPS*, int, int);
   ~MLIAPModelLinear();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
   virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);
diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp
index 554d1a06a76..062e66959a1 100644
--- a/src/MLIAP/mliap_model_quadratic.cpp
+++ b/src/MLIAP/mliap_model_quadratic.cpp
@@ -37,6 +37,15 @@ MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) :
 
 /* ---------------------------------------------------------------------- */
 
+MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) : 
+  MLIAPModel(lmp, nelements_in, nparams_in)
+{
+  nonlinearflag = 1;
+  ndescriptors = sqrt(2*nparams)-1;
+}
+
+/* ---------------------------------------------------------------------- */
+
 MLIAPModelQuadratic::~MLIAPModelQuadratic(){}
 
 /* ----------------------------------------------------------------------
diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h
index 4cfc669bbb0..b597fc76646 100644
--- a/src/MLIAP/mliap_model_quadratic.h
+++ b/src/MLIAP/mliap_model_quadratic.h
@@ -21,6 +21,7 @@ namespace LAMMPS_NS {
 class MLIAPModelQuadratic : public MLIAPModel {
 public:
   MLIAPModelQuadratic(LAMMPS*, char*);
+  MLIAPModelQuadratic(LAMMPS*, int, int);
   ~MLIAPModelQuadratic();
   virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int);
   virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);

From c57564252c2d4137146bde176636631124619daa Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Wed, 1 Jul 2020 15:07:37 -0600
Subject: [PATCH 04/18] Added descriptor file for compute mliap

---
 examples/mliap/compute.mliap.descriptor | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)
 create mode 100644 examples/mliap/compute.mliap.descriptor

diff --git a/examples/mliap/compute.mliap.descriptor b/examples/mliap/compute.mliap.descriptor
new file mode 100644
index 00000000000..7fa4f1bcb09
--- /dev/null
+++ b/examples/mliap/compute.mliap.descriptor
@@ -0,0 +1,19 @@
+# LAMMPS SNAP parameters for MLIAP compute example
+
+# required
+rcutfac 1.0
+twojmax 2
+
+# elements
+
+nelems 2
+elems A B 
+radelems 2.3 2.0
+welems 1.0 0.96
+
+# optional
+
+rfac0 0.99363
+rmin0 0
+bzeroflag 0
+

From 7350dd61d57a1fb3a2be6a1f43d5ac39dba5bdd8 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Wed, 1 Jul 2020 17:06:55 -0600
Subject: [PATCH 05/18] Tweaked input to match snap/in.snap.compute

---
 examples/mliap/compute.mliap.descriptor | 2 +-
 examples/mliap/in.mliap.snap.compute    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/mliap/compute.mliap.descriptor b/examples/mliap/compute.mliap.descriptor
index 7fa4f1bcb09..c3fc432886d 100644
--- a/examples/mliap/compute.mliap.descriptor
+++ b/examples/mliap/compute.mliap.descriptor
@@ -16,4 +16,4 @@ welems 1.0 0.96
 rfac0 0.99363
 rmin0 0
 bzeroflag 0
-
+switchflag 0
diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute
index 813e568387c..bcdb08bc11a 100644
--- a/examples/mliap/in.mliap.snap.compute
+++ b/examples/mliap/in.mliap.snap.compute
@@ -84,7 +84,7 @@ thermo 		100
 
 thermo_style	custom &
 		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 &
-		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] 
+		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[13][12] c_snap[7][12] 
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]

From 2a4e51fa386a4d3a7382a461976dd4b53ee4c88a Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Wed, 1 Jul 2020 17:11:02 -0600
Subject: [PATCH 06/18] Fixed a few problems, still not quite matching compute
 snap

---
 src/MLIAP/compute_mliap.cpp | 27 +++++++++------------------
 1 file changed, 9 insertions(+), 18 deletions(-)

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index 69c00fc94c5..46c062ecaaf 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -260,8 +260,7 @@ void ComputeMLIAP::compute_array()
 
   // compute descriptors, if needed
 
-  if (model->nonlinearflag)
-    descriptor->forward(map, list, descriptors);
+  descriptor->forward(map, list, descriptors);
 
   // calculate descriptor contributions to parameter gradients
   // and gamma = double gradient w.r.t. parameters and descriptors
@@ -271,27 +270,12 @@ void ComputeMLIAP::compute_array()
   // for SNAP, this is a sparse natoms*nparams*ndescriptors matrix,
   // but in general it could be fully dense. 
  
-  // *******Not implemented yet*****************
-  // This should populate the energy row and gamma
-  // For the linear model energy row will look just like the Bi accumulation
-  // in ComputeSNAP i.e. accumulating the intput descriptors vector, 
-  // while gamma will be just 1's and 0's
-  // For the quadratic model, the energy row will be similar,
-  // while gamma will be 1's, 0's and Bi's  
-
   model->param_gradient(map, list, descriptors, gamma_row_index, 
                         gamma_col_index, gamma, egradient);
 
 
   // calculate descriptor gradient contributions to parameter gradients
 
-  // *******Not implemented yet*****************
-  // This will just take gamma and multiply it with
-  // descriptor gradient contributions i.e. dblist
-  // this will resemble snadi accumulation in ComputeSNAP
-
-  // descriptor->param_backward(list, gamma, snadi);
-
   descriptor->param_backward(map, list, gamma_nnz, gamma_row_index, 
                              gamma_col_index, gamma, mliap_peratom,
                              yoffset, zoffset);
@@ -328,6 +312,14 @@ void ComputeMLIAP::compute_array()
 
   dbdotr_compute();
 
+  // copy descriptor gradient contributions to global array
+
+  for (int itype = 0; itype < atom->ntypes; itype++) {
+    const int typeoffset_global = nperdim*itype;
+    for (int icoeff = 0; icoeff < nperdim; icoeff++)
+      mliap[0][icoeff+typeoffset_global] = egradient[icoeff+typeoffset_global];
+  }
+
   // sum up over all processes
 
   MPI_Allreduce(&mliap[0][0],&mliapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
@@ -337,7 +329,6 @@ void ComputeMLIAP::compute_array()
   int irow = 0;
   double reference_energy = c_pe->compute_scalar();
   mliapall[irow++][lastcol] = reference_energy;
-  printf("Reference energy = %g %g %g %d %d\n",reference_energy,mliapall[irow-1][lastcol],array[irow-1][lastcol],irow-1,lastcol);
 
   // assign virial stress to last column
   // switch to Voigt notation

From bc365117676611bb44b5dd582788c19248d549f8 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 14:32:47 -0600
Subject: [PATCH 07/18] Correctly reproduced examples/in.snap.compute, not yet
 for quadratic case

---
 src/MLIAP/compute_mliap.cpp         | 12 ++++------
 src/MLIAP/mliap_descriptor_snap.cpp | 35 +++++++++++++----------------
 src/thermo.cpp                      |  1 +
 3 files changed, 21 insertions(+), 27 deletions(-)

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index 46c062ecaaf..b32078bdc13 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -93,8 +93,8 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   nperdim = nparams;
   ndims_force = 3;
   ndims_virial = 6;
-  yoffset = nperdim;
-  zoffset = 2*nperdim;
+  yoffset = nperdim*atom->ntypes;
+  zoffset = 2*yoffset;
   natoms = atom->natoms;
   size_array_rows = 1+ndims_force*natoms+ndims_virial;
   size_array_cols = nperdim*atom->ntypes+1;
@@ -169,14 +169,12 @@ void ComputeMLIAP::init()
 
   // allocate memory for global array
 
-  printf("size_array_rows = %d size_array_cols = %d\n",size_array_rows,size_array_cols);
   memory->create(mliap,size_array_rows,size_array_cols,
                  "mliap:mliap");
   memory->create(mliapall,size_array_rows,size_array_cols,
                  "mliap:mliapall");
   array = mliapall;
 
-  printf("nelements = %d nparams = %d\n",nelements,nparams);
   memory->create(egradient,nelements*nparams,"ComputeMLIAP:egradient");
 
   // find compute for reference energy
@@ -227,7 +225,6 @@ void ComputeMLIAP::compute_array()
   if (atom->nmax > nmax) {
     memory->destroy(mliap_peratom);
     nmax = atom->nmax;
-    printf("nmax = %d size_peratom = %d\n",nmax,size_peratom);
     memory->create(mliap_peratom,nmax,size_peratom,
                    "mliap:mliap_peratom");
   }
@@ -256,7 +253,6 @@ void ComputeMLIAP::compute_array()
     memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
     gamma_max = list->inum;
   }
-  printf("gamma_max %d %d %d %d %d %d %p\n",gamma_max, ndescriptors, gamma_nnz, nelements, nparams, list->inum, list);
 
   // compute descriptors, if needed
 
@@ -283,7 +279,7 @@ void ComputeMLIAP::compute_array()
   // accumulate descriptor gradient contributions to global array
 
   for (int itype = 0; itype < atom->ntypes; itype++) {
-    const int typeoffset_local = ndims_peratom*nperdim*itype;
+    const int typeoffset_local = nperdim*itype;
     const int typeoffset_global = nperdim*itype;
     for (int icoeff = 0; icoeff < nperdim; icoeff++) {
       int irow = 1;
@@ -360,7 +356,7 @@ void ComputeMLIAP::dbdotr_compute()
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++)
     for (int itype = 0; itype < atom->ntypes; itype++) {
-      const int typeoffset_local = ndims_peratom*nperdim*itype;
+      const int typeoffset_local = nperdim*itype;
       const int typeoffset_global = nperdim*itype;
       double *snadi = mliap_peratom[i]+typeoffset_local;
       for (int icoeff = 0; icoeff < nperdim; icoeff++) {
diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp
index ffb4ab65ff5..bdbeaff43eb 100644
--- a/src/MLIAP/mliap_descriptor_snap.cpp
+++ b/src/MLIAP/mliap_descriptor_snap.cpp
@@ -228,9 +228,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
     // add to Fi, subtract from Fj
 
     snaptr->compute_yi(beta[ii]);
-    //for (int q=0; q<snaptr->idxu_max*2; q++){
-    //  fprintf(screen, "%i %f\n",q, snaptr->ylist_r[q]);
-    //}
 
     for (int jj = 0; jj < ninside; jj++) {
       int j = snaptr->inside[jj];
@@ -294,11 +291,11 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
     const double ytmp = x[i][1];
     const double ztmp = x[i][2];
     const int itype = type[i];
-    int ielem = 0;
-    if (chemflag)
+    int ielem = itype-1;
+    // element map not yet implemented
+    // if (chemflag)
     // element map not yet implemented
     //   ielem = map[itype];
-      const int ielem = itype-1;
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -322,14 +319,17 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
       delz = x[j][2] - ztmp;
       rsq = delx*delx + dely*dely + delz*delz;
       int jtype = type[j];
+      int jelem = jtype-1;
       // element map not yet implemented
-      // int jelem = map[jtype];
-      const int jelem = jtype-1;
+      // if (chemflag)
+      // jelem = map[jtype];
+
       if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
+      // element map not yet implemented
 	snaptr->wj[ninside] = wjelem[jelem];
 	snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
         snaptr->element[ninside] = jelem; // element index for chem snap
@@ -337,14 +337,13 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
       }
     }
 
-    if(chemflag)
+    if (chemflag)
       snaptr->compute_ui(ninside, ielem);
     else
       snaptr->compute_ui(ninside, 0);
 
-
     snaptr->compute_zi();
-    if(chemflag)
+    if (chemflag)
       snaptr->compute_bi(ielem);
     else
       snaptr->compute_bi(0);
@@ -352,12 +351,12 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
     for (int jj = 0; jj < ninside; jj++) {
       const int j = snaptr->inside[jj];
 
-    if(chemflag)
-      snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                             snaptr->rcutij[jj], jj, snaptr->element[jj]);
-    else
-      snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                             snaptr->rcutij[jj], jj, 0);
+      if(chemflag)
+        snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+      else
+        snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                               snaptr->rcutij[jj],jj, 0);
 
       snaptr->compute_dbidrj();
 
@@ -541,11 +540,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
     cut = 2.0*radelem[ielem]*rcutfac;
     if (cut > cutmax) cutmax = cut;
     cutsq[ielem][ielem] = cut*cut;
-    printf("ielem = %d ielem = %d cutsq[i][i] = %g\n",ielem, ielem, cutsq[ielem][ielem]);
     for(int jelem = ielem+1; jelem < nelements; jelem++) {
       cut = (radelem[ielem]+radelem[jelem])*rcutfac;
       cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut*cut;
-      printf("ielem = %d jelem = %d cutsq[i][j] = %g\n",ielem, jelem, cutsq[ielem][jelem]);
     }
   }
 }
diff --git a/src/thermo.cpp b/src/thermo.cpp
index ea476dda04c..339b641bb2c 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -315,6 +315,7 @@ void Thermo::header()
 
   std::string hdr;
   for (int i = 0; i < nfield; i++) hdr +=  keyword[i] + std::string(" ");
+  hdr += "\n";
   if (me == 0) utils::logmesg(lmp,hdr);
 }
 

From bf8043fdb8c4c91266085924c5dc98c182187c60 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 14:34:11 -0600
Subject: [PATCH 08/18] Correctly reproduced examples/in.snap.compute, not yet
 for quadratic case

---
 examples/mliap/in.mliap.quadratic.compute | 95 +++++++++++++++++++++++
 examples/mliap/in.mliap.snap.compute      | 10 +--
 2 files changed, 100 insertions(+), 5 deletions(-)
 create mode 100644 examples/mliap/in.mliap.quadratic.compute

diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute
new file mode 100644
index 00000000000..0f9ba702909
--- /dev/null
+++ b/examples/mliap/in.mliap.quadratic.compute
@@ -0,0 +1,95 @@
+# Demonstrate MLIAP quadratic compute
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	ny equal ${nrep}
+variable 	nz equal ${nrep}
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+create_box	2 box
+create_atoms	2 box
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_coeff 	* * ${zblz} ${zblz}
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 1
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string &
+"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	vb all snav/atom ${snap_options}
+compute 	db all snad/atom ${snap_options}
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+group 		snapgroup2 type 2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_100 equal c_db[2][100]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(0.5*(B_{222}^i)^2, all i of type 2) 
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j) 
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom &
+		pe            pxy            c_bsum2[20]   c_vbsum[220]    v_db_2_100 &
+		c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][37] 
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute
index bcdb08bc11a..4b88a3e52d9 100644
--- a/examples/mliap/in.mliap.snap.compute
+++ b/examples/mliap/in.mliap.snap.compute
@@ -65,7 +65,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_30 equal c_db[2][30]
+variable	db_2_25 equal c_db[2][25]
 
 # set up compute snap generating global array
 
@@ -77,14 +77,14 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(B_{000}^i, all i of type 2) 
-#                4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
 thermo_style	custom &
-		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 &
-		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[13][12] c_snap[7][12] 
+		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 &
+		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]

From cae9788d4296e1ad931f3a4b131ce9f2052e0d97 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 14:43:48 -0600
Subject: [PATCH 09/18] Changed the compute examples to expose effect of
 different list orderings

---
 examples/snap/in.snap.compute           | 10 +++++-----
 examples/snap/in.snap.compute.quadratic | 10 +++++-----
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/examples/snap/in.snap.compute b/examples/snap/in.snap.compute
index e698a134a60..b0c7314882e 100644
--- a/examples/snap/in.snap.compute
+++ b/examples/snap/in.snap.compute
@@ -70,7 +70,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_30 equal c_db[2][30]
+variable	db_2_25 equal c_db[2][25]
 
 # set up compute snap generating global array
 
@@ -82,14 +82,14 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(B_{000}^i, all i of type 2) 
-#                4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
 thermo_style	custom &
-		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 &
-		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] 
+		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 &
+		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10] 
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
diff --git a/examples/snap/in.snap.compute.quadratic b/examples/snap/in.snap.compute.quadratic
index 00e46bd3a82..c31ee4286dd 100644
--- a/examples/snap/in.snap.compute.quadratic
+++ b/examples/snap/in.snap.compute.quadratic
@@ -70,7 +70,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_120 equal c_db[2][120]
+variable	db_2_100 equal c_db[2][100]
 
 # set up compute snap generating global array
 
@@ -80,16 +80,16 @@ fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
 thermo 		100
 
 # test output:   1: total potential energy
-#                2: xy component of stress tensor
+#                2: xz component of stress tensor
 #                3: Sum(0.5*(B_{222}^i)^2, all i of type 2) 
 #                4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j) 
-#                5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
 thermo_style	custom &
-		pe            pxy            c_bsum2[20]   c_vbsum[240]   v_db_2_120 & 
-		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
+		pe            pxy            c_bsum2[20]   c_vbsum[220]   v_db_2_100 & 
+		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]

From db43aadf09b67e8643fad51084b9879eeb0387f6 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 15:05:14 -0600
Subject: [PATCH 10/18] Got compute working for quadratic model

---
 examples/mliap/in.mliap.quadratic.compute | 2 +-
 src/MLIAP/mliap_model_quadratic.cpp       | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute
index 0f9ba702909..6b346392fa5 100644
--- a/examples/mliap/in.mliap.quadratic.compute
+++ b/examples/mliap/in.mliap.quadratic.compute
@@ -84,7 +84,7 @@ thermo 		100
 
 thermo_style	custom &
 		pe            pxy            c_bsum2[20]   c_vbsum[220]    v_db_2_100 &
-		c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][37] 
+		c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp
index 062e66959a1..ae97b30c476 100644
--- a/src/MLIAP/mliap_model_quadratic.cpp
+++ b/src/MLIAP/mliap_model_quadratic.cpp
@@ -143,7 +143,9 @@ void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list,
 
     i = list->ilist[ii];
     const int itype = type[i];
-    const int ielem = map[itype];
+    // element map not yet implemented
+    // const int ielem = map[itype];
+    const int ielem = itype-1;
     const int elemoffset = nparams*ielem;
 
     // linear contributions

From 07afe1c66df75e43c5c1ef71e4fb6fb0196c6f60 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 16:14:15 -0600
Subject: [PATCH 11/18] Added basic element map

---
 src/MLIAP/compute_mliap.cpp         |  7 +++++++
 src/MLIAP/mliap_descriptor_snap.cpp | 20 ++++----------------
 src/MLIAP/mliap_model_linear.cpp    |  4 +---
 src/MLIAP/mliap_model_quadratic.cpp |  4 +---
 4 files changed, 13 insertions(+), 22 deletions(-)

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index b32078bdc13..74e0a34ee5b 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -106,6 +106,13 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   nmax = 0;
   gamma_max = 0;
 
+  // create a minimal map, placeholder for more general map
+
+  memory->create(map,atom->ntypes+1,"compute_mliap:map");
+
+  for (int i = 1; i <= atom->ntypes; i++)
+    map[i] = i-1;
+
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp
index bdbeaff43eb..87b6c171668 100644
--- a/src/MLIAP/mliap_descriptor_snap.cpp
+++ b/src/MLIAP/mliap_descriptor_snap.cpp
@@ -92,9 +92,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
     const double ytmp = x[i][1];
     const double ztmp = x[i][2];
     const int itype = type[i];
-    // element map not yet implemented
-    // const int ielem = map[itype];
-    const int ielem = itype-1;
+    const int ielem = map[itype];
 
     jlist = list->firstneigh[i];
     jnum = list->numneigh[i];
@@ -118,9 +116,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
       delz = x[j][2] - ztmp;
       rsq = delx*delx + dely*dely + delz*delz;
       int jtype = type[j];
-      // element map not yet implemented
-      // const int jelem = map[jtype];
-      const int jelem = jtype-1;
+      const int jelem = map[jtype];
 
       if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
@@ -291,11 +287,7 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
     const double ytmp = x[i][1];
     const double ztmp = x[i][2];
     const int itype = type[i];
-    int ielem = itype-1;
-    // element map not yet implemented
-    // if (chemflag)
-    // element map not yet implemented
-    //   ielem = map[itype];
+    const int ielem = map[itype];
 
     jlist = firstneigh[i];
     jnum = numneigh[i];
@@ -319,17 +311,13 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
       delz = x[j][2] - ztmp;
       rsq = delx*delx + dely*dely + delz*delz;
       int jtype = type[j];
-      int jelem = jtype-1;
-      // element map not yet implemented
-      // if (chemflag)
-      // jelem = map[jtype];
+      const int jelem = map[jtype];
 
       if (rsq < cutsq[ielem][jelem]) {
         snaptr->rij[ninside][0] = delx;
         snaptr->rij[ninside][1] = dely;
         snaptr->rij[ninside][2] = delz;
         snaptr->inside[ninside] = j;
-      // element map not yet implemented
 	snaptr->wj[ninside] = wjelem[jelem];
 	snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
         snaptr->element[ninside] = jelem; // element index for chem snap
diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp
index a6d098ed265..f597d10b332 100644
--- a/src/MLIAP/mliap_model_linear.cpp
+++ b/src/MLIAP/mliap_model_linear.cpp
@@ -117,9 +117,7 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list,
 
     i = list->ilist[ii];
     const int itype = type[i];
-    // element map not yet implemented
-    // const int ielem = map[itype];
-    const int ielem = itype-1;
+    const int ielem = map[itype];
     const int elemoffset = nparams*ielem;
 
     int l = elemoffset+1;
diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp
index ae97b30c476..062e66959a1 100644
--- a/src/MLIAP/mliap_model_quadratic.cpp
+++ b/src/MLIAP/mliap_model_quadratic.cpp
@@ -143,9 +143,7 @@ void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list,
 
     i = list->ilist[ii];
     const int itype = type[i];
-    // element map not yet implemented
-    // const int ielem = map[itype];
-    const int ielem = itype-1;
+    const int ielem = map[itype];
     const int elemoffset = nparams*ielem;
 
     // linear contributions

From 28bbb6afbc9189c02ef2ab3f3ff11162e8935b8a Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Thu, 2 Jul 2020 17:09:10 -0600
Subject: [PATCH 12/18] Doc page for compute mliap

---
 doc/src/compute_mliap.rst | 154 ++++++++++++++++++++++++++++++++++++++
 doc/src/pair_mliap.rst    |   8 +-
 2 files changed, 158 insertions(+), 4 deletions(-)
 create mode 100644 doc/src/compute_mliap.rst

diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst
new file mode 100644
index 00000000000..c04a647993e
--- /dev/null
+++ b/doc/src/compute_mliap.rst
@@ -0,0 +1,154 @@
+.. index:: compute mliap
+
+compute mliap command
+====================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute mliap
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute mliap model linear nelems nparams descriptor sna InP.mliap.descriptor
+
+Description
+"""""""""""
+
+Compute style *mliap* provides a general interface to the gradient
+of machine-learning interatomic potentials w.r.t. model parameters. 
+It is used primarily for calculating the gradient of energy, force, and
+stress components w.r.t. model parameters, which is useful when training
+:doc:`MLIAP pair_style <pair_mliap>` to match target data.
+It provides separate 
+definitions of the interatomic potential functional form (*model*)
+and the geometric quantities that characterize the atomic positions
+(*descriptor*). By defining *model* and *descriptor* separately, 
+it is possible to use many different models with a given descriptor,
+or many different descriptors with a given model. Currently, the 
+compute supports just two models, *linear* and *quadratic*,
+and one descriptor, *sna*, the SNAP descriptor used by 
+:doc:`pair_style snap <pair_snap>`, including the linear, quadratic,
+and chem variants. Work is currently underway to extend
+the interface to handle neural network energy models,
+and it is also straightforward to add new descriptor styles.
+
+The compute *mliap* command must be followed by two keywords
+*model* and *descriptor* in either order. 
+
+The *model* keyword is followed by a model style, currently limited to
+either *linear* or *quadratic*. In both cases,
+this is followed by two arguments. *nelems* is the number of elements.
+It must be equal to the number of LAMMPS atom types. *nparams*
+is the number of parameters per element for this model i.e.
+the number of paramter gradients for each element. Note these definitions
+are identical to those of *nelems* and *nparams* in the 
+:doc:`pair_style mliap <pair_mliap>` model file.
+ 
+The *descriptor* keyword is followed by a descriptor style, and additional arguments.
+Currently the only descriptor style is *sna*, indicating the bispectrum component 
+descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of 
+:doc:`pair_style snap <pair_snap>`.
+The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped 
+out by the LAMMPS command parser. A single additional argument specifies the descriptor filename 
+containing the parameters and setting used by the SNAP descriptor. 
+The descriptor filename usually ends in the *.mliap.descriptor* extension.
+The format of this file is identical to descriptor file in the 
+:doc:`pair_style mliap <pair_mliap>`, and is described in detail
+there. 
+
+.. note::
+
+The number of LAMMPS atom types (and the value of *nelems* in the model)
+must match the value of *nelems* in the descriptor file. 
+
+Compute *mliap* calculates a global array containing gradient information.
+The number of columns in the array is the :math:`nelems \times nparams + 1`.
+The first row of the array contain the derivative of potential energy w.r.t to
+each parameter and each element. The last six rows
+of the array contain the corresponding derivatives of the
+virial stress tensor, listed in Voigt notation: *pxx*, *pyy*, *pzz*,
+*pyz*, *pxz*, *pxy*. In between 3\*\ *N* rows containing the derivatives
+of the force components. 
+
+The element in the last column of each row contains
+the potential energy, force, or stress, according to the row.
+These quantities correspond to the user-specified reference potential
+that must be subtracted from the target data when fitting SNAP.
+The potential energy calculation uses the built in compute *thermo_pe*.
+The stress calculation uses a compute called *snap_press* that is
+automatically created behind the scenes, according to the following
+command:
+
+.. code-block:: LAMMPS
+
+   compute snap_press all pressure NULL virial
+
+See section below on output for a detailed explanation of the data
+layout in the global array.
+
+Atoms not in the group do not contribute to this compute. 
+Neighbor atoms not in the group do not contribute to this compute.
+
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (i.e. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently.
+
+.. note::
+
+   If you have a bonded system, then the settings of
+   :doc:`special_bonds <special_bonds>` command can remove pairwise
+   interactions between atoms in the same bond, angle, or dihedral.  This
+   is the default setting for the :doc:`special_bonds <special_bonds>`
+   command, and means those pairwise interactions do not appear in the
+   neighbor list.  Because this fix uses the neighbor list, it also means
+   those pairs will not be included in the calculation.  One way to get
+   around this, is to write a dump file, and use the :doc:`rerun <rerun>`
+   command to compute the bispectrum components for snapshots in the dump
+   file.  The rerun script can use a :doc:`special_bonds <special_bonds>`
+   command that includes all pairs in the neighbor list.
+
+----------
+
+**Output info:**
+
+Compute *mliap* evaluates a global array.
+The columns are arranged into
+*nelems* blocks, listed in order of element *I*\ . Each block
+contains one column for each of the *nparams* model parameters. 
+A final column contains the corresponding energy, force component
+on an atom, or virial stress component. The rows of the array appear
+in the following order:
+
+* 1 row: Derivatives of potential energy w.r.t. each parameter of each element.
+* 3\*\ *N* rows: Derivatives of force components. x, y, and z components of 
+force on atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
+* 6 rows: Derivatives of virial stress tensor  w.r.t. each parameter of each element.
+The ordering of the rows follows Voigt notation: *pxx*, *pyy*, *pzz*,
+*pyz*, *pxz*, *pxy*.
+
+These values can be accessed by any command that uses per-atom values
+from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
+page for an overview of LAMMPS output options.
+
+Restrictions
+""""""""""""
+
+This compute is part of the MLIAP package.  It is only enabled if
+LAMMPS was built with that package.  In addition, building LAMMPS with the MLIAP package
+requires building LAMMPS with the SNAP package.
+See the :doc:`Build package <Build_package>` doc page for more info.
+doc page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style snap <pair_mliap>`
+
+**Default:** none
diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst
index b569a269106..c65c4c66d54 100644
--- a/doc/src/pair_mliap.rst
+++ b/doc/src/pair_mliap.rst
@@ -47,9 +47,9 @@ where N is the number of LAMMPS atom types.
 The *model* keyword is followed by a model style, currently limited to
 either *linear* or *quadratic*. In both cases,
 this is followed by a single argument specifying the model filename containing the 
-linear or quadratic coefficients for a set of elements. 
+parameters for a set of elements. 
 The model filename usually ends in the *.mliap.model* extension.
-It may contain coefficients for many elements. The only requirement is that it
+It may contain parameters for many elements. The only requirement is that it
 contain at least those element names appearing in the
 *pair_coeff* command.
 
@@ -58,10 +58,10 @@ but follows a strict format after that. The first non-blank non-comment
 line must contain two integers:
 
 * nelems  = Number of elements
-* ncoeff = Number of coefficients
+* nparams = Number of parameters
 
 This is followed by one block for each of the *nelem* elements.
-Each block consists of *ncoeff* coefficients, one per line.
+Each block consists of *nparams* parameters, one per line.
 Note that this format is similar, but not identical to that used
 for the :doc:`pair_style snap <pair_snap>` coefficient file.
 Specifically, the line containing the element weight and radius is omitted, 

From 9d0aee7426a1b9efded94f1b49714833d0c31ec7 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Fri, 3 Jul 2020 13:43:51 -0600
Subject: [PATCH 13/18] Added doc page for compute mliap and updated examples

---
 doc/src/compute_mliap.rst                     |  80 ++++---
 doc/src/compute_sna_atom.rst                  |   4 +-
 doc/src/pair_mliap.rst                        |  21 +-
 .../log.03Jul20.mliap.quadratic.compute.g++.1 | 199 +++++++++++++++++
 .../log.03Jul20.mliap.quadratic.compute.g++.4 | 200 ++++++++++++++++++
 .../log.03Jul20.mliap.snap.compute.g++.1      | 199 +++++++++++++++++
 .../log.03Jul20.mliap.snap.compute.g++.4      | 200 ++++++++++++++++++
 examples/snap/in.snap.compute.quadratic       |   4 +-
 ...e.g++.1 => log.03Jul20.snap.compute.g++.1} |  26 ++-
 ...e.g++.4 => log.03Jul20.snap.compute.g++.4} |  30 ++-
 ... log.03Jul20.snap.compute.quadratic.g++.1} |  26 ++-
 ... log.03Jul20.snap.compute.quadratic.g++.4} |  30 ++-
 12 files changed, 921 insertions(+), 98 deletions(-)
 create mode 100644 examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1
 create mode 100644 examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4
 create mode 100644 examples/mliap/log.03Jul20.mliap.snap.compute.g++.1
 create mode 100644 examples/mliap/log.03Jul20.mliap.snap.compute.g++.4
 rename examples/snap/{log.27Nov18.snap.compute.g++.1 => log.03Jul20.snap.compute.g++.1} (87%)
 rename examples/snap/{log.27Nov18.snap.compute.g++.4 => log.03Jul20.snap.compute.g++.4} (86%)
 rename examples/snap/{log.27Nov18.snap.compute.quadratic.g++.1 => log.03Jul20.snap.compute.quadratic.g++.1} (87%)
 rename examples/snap/{log.27Nov18.snap.compute.quadratic.g++.4 => log.03Jul20.snap.compute.quadratic.g++.4} (86%)

diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst
index c04a647993e..91253b29c2e 100644
--- a/doc/src/compute_mliap.rst
+++ b/doc/src/compute_mliap.rst
@@ -1,21 +1,36 @@
 .. index:: compute mliap
 
 compute mliap command
-====================
+=====================
 
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   compute mliap
+   compute ID group-ID mliap ... keyword values ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* mliap = style name of this compute command
+* two keyword/value pairs must be appended
+* keyword = *model* or *descriptor*
+
+  .. parsed-literal::
+
+       *model* values = style Nelems Nparams
+         style = *linear* or *quadratic*
+         Nelems = number of elements
+         Nparams = number of parameters per element 
+       *descriptor* values = style filename
+         style = *sna*
+         filename = name of file containing descriptor definitions
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   compute mliap model linear nelems nparams descriptor sna InP.mliap.descriptor
+   compute mliap model linear 2 31 descriptor sna Ta06A.mliap.descriptor
 
 Description
 """""""""""
@@ -24,7 +39,7 @@ Compute style *mliap* provides a general interface to the gradient
 of machine-learning interatomic potentials w.r.t. model parameters. 
 It is used primarily for calculating the gradient of energy, force, and
 stress components w.r.t. model parameters, which is useful when training
-:doc:`MLIAP pair_style <pair_mliap>` to match target data.
+:doc:`mliap pair_style <pair_mliap>` models to match target data.
 It provides separate 
 definitions of the interatomic potential functional form (*model*)
 and the geometric quantities that characterize the atomic positions
@@ -46,7 +61,7 @@ either *linear* or *quadratic*. In both cases,
 this is followed by two arguments. *nelems* is the number of elements.
 It must be equal to the number of LAMMPS atom types. *nparams*
 is the number of parameters per element for this model i.e.
-the number of paramter gradients for each element. Note these definitions
+the number of parameter gradients for each element. Note these definitions
 are identical to those of *nelems* and *nparams* in the 
 :doc:`pair_style mliap <pair_mliap>` model file.
  
@@ -58,43 +73,44 @@ The \'p\' in SNAP is dropped, because keywords that match pair_styles are silent
 out by the LAMMPS command parser. A single additional argument specifies the descriptor filename 
 containing the parameters and setting used by the SNAP descriptor. 
 The descriptor filename usually ends in the *.mliap.descriptor* extension.
-The format of this file is identical to descriptor file in the 
+The format of this file is identical to the descriptor file in the 
 :doc:`pair_style mliap <pair_mliap>`, and is described in detail
 there. 
 
 .. note::
 
-The number of LAMMPS atom types (and the value of *nelems* in the model)
-must match the value of *nelems* in the descriptor file. 
+   The number of LAMMPS atom types (and the value of *nelems* in the model)
+   must match the value of *nelems* in the descriptor file. 
 
 Compute *mliap* calculates a global array containing gradient information.
-The number of columns in the array is the :math:`nelems \times nparams + 1`.
-The first row of the array contain the derivative of potential energy w.r.t to
+The number of columns in the array is :math:`nelems \times nparams + 1`.
+The first row of the array contain the derivative of potential energy w.r.t. to
 each parameter and each element. The last six rows
 of the array contain the corresponding derivatives of the
 virial stress tensor, listed in Voigt notation: *pxx*, *pyy*, *pzz*,
-*pyz*, *pxz*, *pxy*. In between 3\*\ *N* rows containing the derivatives
-of the force components. 
+*pyz*, *pxz*, *pxy*. In between the energy and stress rows are
+the 3\*\ *N* rows containing the derivatives of the force components. 
+See section below on output for a detailed description of how 
+rows and columns are ordered. 
 
 The element in the last column of each row contains
 the potential energy, force, or stress, according to the row.
 These quantities correspond to the user-specified reference potential
-that must be subtracted from the target data when fitting SNAP.
+that must be subtracted from the target data when training a model.
 The potential energy calculation uses the built in compute *thermo_pe*.
-The stress calculation uses a compute called *snap_press* that is
+The stress calculation uses a compute called *mliap_press* that is
 automatically created behind the scenes, according to the following
 command:
 
 .. code-block:: LAMMPS
 
-   compute snap_press all pressure NULL virial
+   compute mliap_press all pressure NULL virial
 
 See section below on output for a detailed explanation of the data
 layout in the global array.
 
 Atoms not in the group do not contribute to this compute. 
 Neighbor atoms not in the group do not contribute to this compute.
-
 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (i.e. each time a snapshot of atoms
 is dumped).  Thus it can be inefficient to compute/dump this quantity
@@ -102,17 +118,19 @@ too frequently.
 
 .. note::
 
-   If you have a bonded system, then the settings of
-   :doc:`special_bonds <special_bonds>` command can remove pairwise
+   If the user-specified reference potentials includes bonded and
+   non-bonded pairwise interactions, then the settings of
+   :doc:`special_bonds <special_bonds>` command can remove pairwise 
    interactions between atoms in the same bond, angle, or dihedral.  This
    is the default setting for the :doc:`special_bonds <special_bonds>`
    command, and means those pairwise interactions do not appear in the
    neighbor list.  Because this fix uses the neighbor list, it also means
-   those pairs will not be included in the calculation.  One way to get
-   around this, is to write a dump file, and use the :doc:`rerun <rerun>`
-   command to compute the bispectrum components for snapshots in the dump
-   file.  The rerun script can use a :doc:`special_bonds <special_bonds>`
-   command that includes all pairs in the neighbor list.
+   those pairs will not be included in the calculation. The :doc:`rerun <rerun>`
+   command is not an option here, since the reference potential is required
+   for the last column of the global array. A work-around is to prevent
+   pairwise interactions from being removed by explicitly adding a 
+   *tiny* positive value for every pairwise interaction that would otherwise be
+   set to zero in the :doc:`special_bonds <special_bonds>` command.
 
 ----------
 
@@ -127,15 +145,14 @@ on an atom, or virial stress component. The rows of the array appear
 in the following order:
 
 * 1 row: Derivatives of potential energy w.r.t. each parameter of each element.
-* 3\*\ *N* rows: Derivatives of force components. x, y, and z components of 
-force on atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
-* 6 rows: Derivatives of virial stress tensor  w.r.t. each parameter of each element.
-The ordering of the rows follows Voigt notation: *pxx*, *pyy*, *pzz*,
-*pyz*, *pxz*, *pxy*.
+* 3\*\ *N* rows: Derivatives of force components. x, y, and z components of force on atom *i* appearing in consecutive rows. The atoms are sorted based on atom ID.
+* 6 rows: Derivatives of virial stress tensor  w.r.t. each parameter of each element. The ordering of the rows follows Voigt notation: *pxx*, *pyy*, *pzz*, *pyz*, *pxz*, *pxy*.
 
-These values can be accessed by any command that uses per-atom values
+These values can be accessed by any command that uses a global array
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
-page for an overview of LAMMPS output options.
+page for an overview of LAMMPS output options. To see how this command
+can be used within a Python workflow to train machine-learning interatomic
+potentials, see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
 
 Restrictions
 """"""""""""
@@ -144,11 +161,10 @@ This compute is part of the MLIAP package.  It is only enabled if
 LAMMPS was built with that package.  In addition, building LAMMPS with the MLIAP package
 requires building LAMMPS with the SNAP package.
 See the :doc:`Build package <Build_package>` doc page for more info.
-doc page for more info.
 
 Related commands
 """"""""""""""""
 
-:doc:`pair_style snap <pair_mliap>`
+:doc:`pair_style mliap <pair_mliap>`
 
 **Default:** none
diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst
index 4e42c4523b8..cf9339a30d0 100644
--- a/doc/src/compute_sna_atom.rst
+++ b/doc/src/compute_sna_atom.rst
@@ -391,7 +391,9 @@ of :math:`K N_{elem}^3` columns.
 
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
-page for an overview of LAMMPS output options.
+page for an overview of LAMMPS output options. To see how this command
+can be used within a Python workflow to train SNAP potentials, 
+see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
 
 Restrictions
 """"""""""""
diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst
index c65c4c66d54..d4fd897b954 100644
--- a/doc/src/pair_mliap.rst
+++ b/doc/src/pair_mliap.rst
@@ -8,7 +8,19 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   pair_style mliap
+   pair_style mliap ... keyword values ...
+
+* two keyword/value pairs must be appended
+* keyword = *model* or *descriptor*
+
+  .. parsed-literal::
+
+       *model* values = style filename
+         style = *linear* or *quadratic*
+         filename = name of file containing model definitions
+       *descriptor* values = style filename
+         style = *sna*
+         filename = name of file containing descriptor definitions
 
 Examples
 """"""""
@@ -23,7 +35,7 @@ Description
 """""""""""
 
 Pair style *mliap* provides a general interface to families of 
-machine-learning interatomic potentials. It provides separate 
+machine-learning interatomic potentials. It allows separate 
 definitions of the interatomic potential functional form (*model*)
 and the geometric quantities that characterize the atomic positions
 (*descriptor*). By defining *model* and *descriptor* separately, 
@@ -34,6 +46,9 @@ and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap <pa
 and chem variants. Work is currently underway to extend
 the interface to handle neural network energy models,
 and it is also straightforward to add new descriptor styles.
+In order to train a model, it is useful to know the gradient or derivative
+of energy, force, and stress w.r.t. model parameters. This information
+can be accessed using the related :doc:`compute mliap <compute_mliap>` command.
 
 The pair_style *mliap* command must be followed by two keywords
 *model* and *descriptor* in either order. A single
@@ -131,6 +146,6 @@ See the :doc:`Build package <Build_package>` doc page for more info.
 Related commands
 """"""""""""""""
 
-:doc:`pair_style snap  <pair_snap>`,
+:doc:`pair_style snap  <pair_snap>`, :doc:`compute mliap <compute_mliap>`
 
 **Default:** none
diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1
new file mode 100644
index 00000000000..d5226cec890
--- /dev/null
+++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1
@@ -0,0 +1,199 @@
+LAMMPS (15 Jun 2020)
+# Demonstrate MLIAP quadratic compute
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	nx equal 1
+variable 	ny equal ${nrep}
+variable 	ny equal 1
+variable 	nz equal ${nrep}
+variable 	nz equal 1
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+lattice         bcc 2
+Lattice spacing in x,y,z = 2 2 2
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+region		box block 0 1 0 ${ny} 0 ${nz}
+region		box block 0 1 0 1 0 ${nz}
+region		box block 0 1 0 1 0 1
+create_box	2 box
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
+  1 by 1 by 1 MPI processor grid
+create_atoms	2 box
+Created 2 atoms
+  create_atoms CPU = 0.000 seconds
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_style 	zbl 4 ${zblcutouter}
+pair_style 	zbl 4 4.8
+pair_coeff 	* * ${zblz} ${zblz}
+pair_coeff 	* * 73 ${zblz}
+pair_coeff 	* * 73 73
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 1
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+compute 	vb all snav/atom ${snap_options}
+compute 	vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+compute 	db all snad/atom ${snap_options}
+compute 	db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+0 atoms in group snapgroup1
+group 		snapgroup2 type 2
+2 atoms in group snapgroup2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_100 equal c_db[2][100]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
+SNAP keyword rcutfac 1.0 
+SNAP keyword twojmax 2 
+SNAP keyword nelems 2 
+SNAP keyword elems A 
+SNAP keyword radelems 2.3 
+SNAP keyword welems 1.0 
+SNAP keyword rfac0 0.99363 
+SNAP keyword rmin0 0 
+SNAP keyword bzeroflag 0 
+SNAP keyword switchflag 0 
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[220]    v_db_2_100 		c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
+run             0
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.8
+  ghost atom cutoff = 6.8
+  binsize = 3.4, bins = 1 1 1
+  5 neighbor lists, perpetual/occasional/extra = 1 4 0
+  (1) pair zbl, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) compute sna/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (3) compute snav/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (4) compute snad/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (5) compute mliap, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 22.45 | 22.45 | 22.45 Mbytes
+PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 
+   322.86952    1505558.1 4.2492771e+08     -4952473     28484035    322.86952    1505558.1 4.2492771e+08     -4952473     28484035 
+Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
+
+100.0% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 1e-06      |            |       |100.00
+
+Nlocal:    2 ave 2 max 2 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    853 ave 853 max 853 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    330 ave 330 max 330 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:  660 ave 660 max 660 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 660
+Ave neighs/atom = 330
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:00
diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4
new file mode 100644
index 00000000000..dfa6483b536
--- /dev/null
+++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4
@@ -0,0 +1,200 @@
+LAMMPS (15 Jun 2020)
+# Demonstrate MLIAP quadratic compute
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	nx equal 1
+variable 	ny equal ${nrep}
+variable 	ny equal 1
+variable 	nz equal ${nrep}
+variable 	nz equal 1
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+lattice         bcc 2
+Lattice spacing in x,y,z = 2 2 2
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+region		box block 0 1 0 ${ny} 0 ${nz}
+region		box block 0 1 0 1 0 ${nz}
+region		box block 0 1 0 1 0 1
+create_box	2 box
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
+  1 by 2 by 2 MPI processor grid
+create_atoms	2 box
+Created 2 atoms
+  create_atoms CPU = 0.001 seconds
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_style 	zbl 4 ${zblcutouter}
+pair_style 	zbl 4 4.8
+pair_coeff 	* * ${zblz} ${zblz}
+pair_coeff 	* * 73 ${zblz}
+pair_coeff 	* * 73 73
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 1
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+compute 	vb all snav/atom ${snap_options}
+compute 	vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+compute 	db all snad/atom ${snap_options}
+compute 	db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 1 bzeroflag 0 switchflag 0
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+0 atoms in group snapgroup1
+group 		snapgroup2 type 2
+2 atoms in group snapgroup2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_100 equal c_db[2][100]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21
+SNAP keyword rcutfac 1.0 
+SNAP keyword twojmax 2 
+SNAP keyword nelems 2 
+SNAP keyword elems A 
+SNAP keyword radelems 2.3 
+SNAP keyword welems 1.0 
+SNAP keyword rfac0 0.99363 
+SNAP keyword rmin0 0 
+SNAP keyword bzeroflag 0 
+SNAP keyword switchflag 0 
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[220]    v_db_2_100 		c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42]
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
+run             0
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.8
+  ghost atom cutoff = 6.8
+  binsize = 3.4, bins = 1 1 1
+  5 neighbor lists, perpetual/occasional/extra = 1 4 0
+  (1) pair zbl, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) compute sna/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (3) compute snav/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (4) compute snad/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (5) compute mliap, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
+Per MPI rank memory allocation (min/avg/max) = 22.18 | 22.18 | 22.18 Mbytes
+PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 
+   322.86952    1505558.1 4.2492771e+08     -4952473     28484035    322.86952    1505558.1 4.2492771e+08     -4952473     28484035 
+Loop time of 2e-06 on 4 procs for 0 steps with 2 atoms
+
+100.0% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 2e-06      |            |       |100.00
+
+Nlocal:    0.5 ave 1 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Nghost:    734.5 ave 735 max 734 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Neighs:    82.5 ave 177 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 1 1
+FullNghs:  165 ave 330 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+
+Total # of neighbors = 660
+Ave neighs/atom = 330
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:00
diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1
new file mode 100644
index 00000000000..083d451afab
--- /dev/null
+++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1
@@ -0,0 +1,199 @@
+LAMMPS (15 Jun 2020)
+# Demonstrate bispectrum computes
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	nx equal 1
+variable 	ny equal ${nrep}
+variable 	ny equal 1
+variable 	nz equal ${nrep}
+variable 	nz equal 1
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+lattice         bcc 2
+Lattice spacing in x,y,z = 2 2 2
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+region		box block 0 1 0 ${ny} 0 ${nz}
+region		box block 0 1 0 1 0 ${nz}
+region		box block 0 1 0 1 0 1
+create_box	2 box
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
+  1 by 1 by 1 MPI processor grid
+create_atoms	2 box
+Created 2 atoms
+  create_atoms CPU = 0.001 seconds
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_style 	zbl 4 ${zblcutouter}
+pair_style 	zbl 4 4.8
+pair_coeff 	* * ${zblz} ${zblz}
+pair_coeff 	* * 73 ${zblz}
+pair_coeff 	* * 73 73
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 0
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+compute 	vb all snav/atom ${snap_options}
+compute 	vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+compute 	db all snad/atom ${snap_options}
+compute 	db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+0 atoms in group snapgroup1
+group 		snapgroup2 type 2
+2 atoms in group snapgroup2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_25 equal c_db[2][25]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
+SNAP keyword rcutfac 1.0 
+SNAP keyword twojmax 2 
+SNAP keyword nelems 2 
+SNAP keyword elems A 
+SNAP keyword radelems 2.3 
+SNAP keyword welems 1.0 
+SNAP keyword rfac0 0.99363 
+SNAP keyword rmin0 0 
+SNAP keyword bzeroflag 0 
+SNAP keyword switchflag 0 
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(B_{000}^i, all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
+run             0
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.8
+  ghost atom cutoff = 6.8
+  binsize = 3.4, bins = 1 1 1
+  5 neighbor lists, perpetual/occasional/extra = 1 4 0
+  (1) pair zbl, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) compute sna/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (3) compute snav/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (4) compute snad/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (5) compute mliap, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 10.73 | 10.73 | 10.73 Mbytes
+PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 
+   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
+Loop time of 0 on 1 procs for 0 steps with 2 atoms
+
+0.0% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 0          |            |       |  0.00
+
+Nlocal:    2 ave 2 max 2 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    853 ave 853 max 853 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    330 ave 330 max 330 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:  660 ave 660 max 660 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 660
+Ave neighs/atom = 330
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:00
diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4
new file mode 100644
index 00000000000..ebedaee049d
--- /dev/null
+++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4
@@ -0,0 +1,200 @@
+LAMMPS (15 Jun 2020)
+# Demonstrate bispectrum computes
+
+# initialize simulation
+
+variable 	nsteps index 0
+variable 	nrep equal 1
+variable 	a equal 2.0
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable 	nx equal ${nrep}
+variable 	nx equal 1
+variable 	ny equal ${nrep}
+variable 	ny equal 1
+variable 	nz equal ${nrep}
+variable 	nz equal 1
+
+boundary	p p p
+
+atom_modify	map hash
+lattice         bcc $a
+lattice         bcc 2
+Lattice spacing in x,y,z = 2 2 2
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+region		box block 0 1 0 ${ny} 0 ${nz}
+region		box block 0 1 0 1 0 ${nz}
+region		box block 0 1 0 1 0 1
+create_box	2 box
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
+  1 by 2 by 2 MPI processor grid
+create_atoms	2 box
+Created 2 atoms
+  create_atoms CPU = 0.001 seconds
+
+mass 		* 180.88
+
+displace_atoms 	all random 0.1 0.1 0.1 123456
+
+# set up reference potential
+
+variable 	zblcutinner equal 4
+variable 	zblcutouter equal 4.8
+variable 	zblz equal 73
+pair_style 	zbl ${zblcutinner} ${zblcutouter}
+pair_style 	zbl 4 ${zblcutouter}
+pair_style 	zbl 4 4.8
+pair_coeff 	* * ${zblz} ${zblz}
+pair_coeff 	* * 73 ${zblz}
+pair_coeff 	* * 73 73
+
+# choose SNA parameters
+
+variable 	twojmax equal 2
+variable 	rcutfac equal 1.0
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	radelem1 equal 2.3
+variable 	radelem2 equal 2.0
+variable	wj1 equal 1.0
+variable	wj2 equal 0.96
+variable	quadratic equal 0
+variable	bzero equal 0
+variable	switch equal 0
+variable 	snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
+1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch}
+1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+
+# set up per-atom computes
+
+compute 	b all sna/atom ${snap_options}
+compute 	b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+compute 	vb all snav/atom ${snap_options}
+compute 	vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+compute 	db all snad/atom ${snap_options}
+compute 	db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
+
+# perform sums over atoms
+
+group 		snapgroup1 type 1
+0 atoms in group snapgroup1
+group 		snapgroup2 type 2
+2 atoms in group snapgroup2
+compute         bsum1 snapgroup1 reduce sum c_b[*]
+compute         bsum2 snapgroup2 reduce sum c_b[*]
+# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
+# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
+compute		vbsum all reduce sum c_vb[*]
+# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
+variable	db_2_25 equal c_db[2][25]
+
+# set up compute snap generating global array
+
+compute   	snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6
+SNAP keyword rcutfac 1.0 
+SNAP keyword twojmax 2 
+SNAP keyword nelems 2 
+SNAP keyword elems A 
+SNAP keyword radelems 2.3 
+SNAP keyword welems 1.0 
+SNAP keyword rfac0 0.99363 
+SNAP keyword rmin0 0 
+SNAP keyword bzeroflag 0 
+SNAP keyword switchflag 0 
+fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
+
+thermo 		100
+
+# test output:   1: total potential energy
+#                2: xy component of stress tensor
+#                3: Sum(B_{000}^i, all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#
+#                followed by 5 counterparts from compute snap
+
+thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
+thermo_modify 	norm no
+
+# dump 		mydump_db all custom 1000 dump_db id c_db[*]
+# dump_modify 	mydump_db sort id
+
+# Run MD
+
+run             ${nsteps}
+run             0
+Neighbor list info ...
+  update every 1 steps, delay 10 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.8
+  ghost atom cutoff = 6.8
+  binsize = 3.4, bins = 1 1 1
+  5 neighbor lists, perpetual/occasional/extra = 1 4 0
+  (1) pair zbl, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+  (2) compute sna/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (3) compute snav/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (4) compute snad/atom, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+  (5) compute mliap, occasional
+      attributes: full, newton on
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
+Per MPI rank memory allocation (min/avg/max) = 10.69 | 10.69 | 10.7 Mbytes
+PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 
+   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
+Loop time of 3.25e-06 on 4 procs for 0 steps with 2 atoms
+
+115.4% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0          | 0          | 0          |   0.0 |  0.00
+Output  | 0          | 0          | 0          |   0.0 |  0.00
+Modify  | 0          | 0          | 0          |   0.0 |  0.00
+Other   |            | 3.25e-06   |            |       |100.00
+
+Nlocal:    0.5 ave 1 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Nghost:    734.5 ave 735 max 734 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+Neighs:    82.5 ave 177 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 1 1
+FullNghs:  165 ave 330 max 0 min
+Histogram: 2 0 0 0 0 0 0 0 0 2
+
+Total # of neighbors = 660
+Ave neighs/atom = 330
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:00
diff --git a/examples/snap/in.snap.compute.quadratic b/examples/snap/in.snap.compute.quadratic
index c31ee4286dd..e03d4af3bfe 100644
--- a/examples/snap/in.snap.compute.quadratic
+++ b/examples/snap/in.snap.compute.quadratic
@@ -80,9 +80,9 @@ fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector
 thermo 		100
 
 # test output:   1: total potential energy
-#                2: xz component of stress tensor
+#                2: xy component of stress tensor
 #                3: Sum(0.5*(B_{222}^i)^2, all i of type 2) 
-#                4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j) 
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j) 
 #                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
diff --git a/examples/snap/log.27Nov18.snap.compute.g++.1 b/examples/snap/log.03Jul20.snap.compute.g++.1
similarity index 87%
rename from examples/snap/log.27Nov18.snap.compute.g++.1
rename to examples/snap/log.03Jul20.snap.compute.g++.1
index b189c552da7..8b9fe80a07b 100644
--- a/examples/snap/log.27Nov18.snap.compute.g++.1
+++ b/examples/snap/log.03Jul20.snap.compute.g++.1
@@ -1,6 +1,4 @@
-LAMMPS (20 Nov 2019)
-OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
-  using 1 OpenMP thread(s) per MPI task
+LAMMPS (15 Jun 2020)
 # Demonstrate bispectrum computes
 
 # initialize simulation
@@ -30,11 +28,11 @@ region		box block 0 1 0 ${ny} 0 ${nz}
 region		box block 0 1 0 1 0 ${nz}
 region		box block 0 1 0 1 0 1
 create_box	2 box
-Created orthogonal box = (0 0 0) to (2 2 2)
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
   1 by 1 by 1 MPI processor grid
 create_atoms	2 box
 Created 2 atoms
-  create_atoms CPU = 0.000478029 secs
+  create_atoms CPU = 0.000 seconds
 
 mass 		* 180.88
 
@@ -105,7 +103,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_30 equal c_db[2][30]
+variable	db_2_25 equal c_db[2][25]
 
 # set up compute snap generating global array
 
@@ -118,12 +116,12 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(B_{000}^i, all i of type 2)
-#                4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
-thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
+thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
@@ -166,11 +164,11 @@ Neighbor list info ...
       stencil: full/bin/3d
       bin: standard
 Per MPI rank memory allocation (min/avg/max) = 10.06 | 10.06 | 10.06 Mbytes
-PotEng Pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] 
-   322.86952    1505558.1    364182.88    381.32218   -855.04473    322.86952    1505558.1    364182.88    381.32218   -855.04473 
-Loop time of 9.53674e-07 on 1 procs for 0 steps with 2 atoms
+PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10] 
+   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
+Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
 
-104.9% CPU use with 1 MPI tasks x 1 OpenMP threads
+200.0% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
@@ -180,7 +178,7 @@ Neigh   | 0          | 0          | 0          |   0.0 |  0.00
 Comm    | 0          | 0          | 0          |   0.0 |  0.00
 Output  | 0          | 0          | 0          |   0.0 |  0.00
 Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 9.537e-07  |            |       |100.00
+Other   |            | 1e-06      |            |       |100.00
 
 Nlocal:    2 ave 2 max 2 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
diff --git a/examples/snap/log.27Nov18.snap.compute.g++.4 b/examples/snap/log.03Jul20.snap.compute.g++.4
similarity index 86%
rename from examples/snap/log.27Nov18.snap.compute.g++.4
rename to examples/snap/log.03Jul20.snap.compute.g++.4
index f979123b2a7..ff68bf75d5e 100644
--- a/examples/snap/log.27Nov18.snap.compute.g++.4
+++ b/examples/snap/log.03Jul20.snap.compute.g++.4
@@ -1,6 +1,4 @@
-LAMMPS (20 Nov 2019)
-OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
-  using 1 OpenMP thread(s) per MPI task
+LAMMPS (15 Jun 2020)
 # Demonstrate bispectrum computes
 
 # initialize simulation
@@ -30,11 +28,11 @@ region		box block 0 1 0 ${ny} 0 ${nz}
 region		box block 0 1 0 1 0 ${nz}
 region		box block 0 1 0 1 0 1
 create_box	2 box
-Created orthogonal box = (0 0 0) to (2 2 2)
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
   1 by 2 by 2 MPI processor grid
 create_atoms	2 box
 Created 2 atoms
-  create_atoms CPU = 0.000610113 secs
+  create_atoms CPU = 0.001 seconds
 
 mass 		* 180.88
 
@@ -105,7 +103,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_30 equal c_db[2][30]
+variable	db_2_25 equal c_db[2][25]
 
 # set up compute snap generating global array
 
@@ -118,12 +116,12 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(B_{000}^i, all i of type 2)
-#                4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
-thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[60]    v_db_2_30 		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10]
+thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10]
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
@@ -165,13 +163,13 @@ Neighbor list info ...
       pair build: full/bin/atomonly
       stencil: full/bin/3d
       bin: standard
-WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:936)
-Per MPI rank memory allocation (min/avg/max) = 8.211 | 8.254 | 8.295 Mbytes
-PotEng Pxy c_bsum2[1] c_vbsum[60] v_db_2_30 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] 
-   322.86952    1505558.1    364182.88    381.32218   -855.04473    322.86952    1505558.1    364182.88    381.32218   -855.04473 
-Loop time of 2.38419e-06 on 4 procs for 0 steps with 2 atoms
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
+Per MPI rank memory allocation (min/avg/max) = 9.945 | 9.988 | 10.03 Mbytes
+PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10] 
+   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
+Loop time of 2.75e-06 on 4 procs for 0 steps with 2 atoms
 
-104.9% CPU use with 4 MPI tasks x 1 OpenMP threads
+90.9% CPU use with 4 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
@@ -181,7 +179,7 @@ Neigh   | 0          | 0          | 0          |   0.0 |  0.00
 Comm    | 0          | 0          | 0          |   0.0 |  0.00
 Output  | 0          | 0          | 0          |   0.0 |  0.00
 Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 2.384e-06  |            |       |100.00
+Other   |            | 2.75e-06   |            |       |100.00
 
 Nlocal:    0.5 ave 1 max 0 min
 Histogram: 2 0 0 0 0 0 0 0 0 2
diff --git a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.1 b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.1
similarity index 87%
rename from examples/snap/log.27Nov18.snap.compute.quadratic.g++.1
rename to examples/snap/log.03Jul20.snap.compute.quadratic.g++.1
index c2054023ede..a52810063cf 100644
--- a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.1
+++ b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.1
@@ -1,6 +1,4 @@
-LAMMPS (20 Nov 2019)
-OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
-  using 1 OpenMP thread(s) per MPI task
+LAMMPS (15 Jun 2020)
 # Demonstrate bispectrum computes
 
 # initialize simulation
@@ -30,11 +28,11 @@ region		box block 0 1 0 ${ny} 0 ${nz}
 region		box block 0 1 0 1 0 ${nz}
 region		box block 0 1 0 1 0 1
 create_box	2 box
-Created orthogonal box = (0 0 0) to (2 2 2)
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
   1 by 1 by 1 MPI processor grid
 create_atoms	2 box
 Created 2 atoms
-  create_atoms CPU = 0.000473976 secs
+  create_atoms CPU = 0.001 seconds
 
 mass 		* 180.88
 
@@ -105,7 +103,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_120 equal c_db[2][120]
+variable	db_2_100 equal c_db[2][100]
 
 # set up compute snap generating global array
 
@@ -118,12 +116,12 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
-#                4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
-thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[240]   v_db_2_120 		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
+thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[220]   v_db_2_100 		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
@@ -166,11 +164,11 @@ Neighbor list info ...
       stencil: full/bin/3d
       bin: standard
 Per MPI rank memory allocation (min/avg/max) = 21.78 | 21.78 | 21.78 Mbytes
-PotEng Pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40] 
-   322.86952    1505558.1 4.2492771e+08    7860489.6    -17625699    322.86952    1505558.1 4.2492771e+08    7860489.6    -17625699 
-Loop time of 2.14577e-06 on 1 procs for 0 steps with 2 atoms
+PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40] 
+   322.86952    1505558.1 4.2492771e+08     -4952473     28484035    322.86952    1505558.1 4.2492771e+08     -4952473     28484035 
+Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms
 
-93.2% CPU use with 1 MPI tasks x 1 OpenMP threads
+200.0% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
@@ -180,7 +178,7 @@ Neigh   | 0          | 0          | 0          |   0.0 |  0.00
 Comm    | 0          | 0          | 0          |   0.0 |  0.00
 Output  | 0          | 0          | 0          |   0.0 |  0.00
 Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 2.146e-06  |            |       |100.00
+Other   |            | 1e-06      |            |       |100.00
 
 Nlocal:    2 ave 2 max 2 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
diff --git a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.4 b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.4
similarity index 86%
rename from examples/snap/log.27Nov18.snap.compute.quadratic.g++.4
rename to examples/snap/log.03Jul20.snap.compute.quadratic.g++.4
index 819397f7455..ba674cd4fa8 100644
--- a/examples/snap/log.27Nov18.snap.compute.quadratic.g++.4
+++ b/examples/snap/log.03Jul20.snap.compute.quadratic.g++.4
@@ -1,6 +1,4 @@
-LAMMPS (20 Nov 2019)
-OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:93)
-  using 1 OpenMP thread(s) per MPI task
+LAMMPS (15 Jun 2020)
 # Demonstrate bispectrum computes
 
 # initialize simulation
@@ -30,11 +28,11 @@ region		box block 0 1 0 ${ny} 0 ${nz}
 region		box block 0 1 0 1 0 ${nz}
 region		box block 0 1 0 1 0 1
 create_box	2 box
-Created orthogonal box = (0 0 0) to (2 2 2)
+Created orthogonal box = (0.0 0.0 0.0) to (2.0 2.0 2.0)
   1 by 2 by 2 MPI processor grid
 create_atoms	2 box
 Created 2 atoms
-  create_atoms CPU = 0.000118971 secs
+  create_atoms CPU = 0.001 seconds
 
 mass 		* 180.88
 
@@ -105,7 +103,7 @@ compute         bsum2 snapgroup2 reduce sum c_b[*]
 # fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
 compute		vbsum all reduce sum c_vb[*]
 # fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
-variable	db_2_120 equal c_db[2][120]
+variable	db_2_100 equal c_db[2][100]
 
 # set up compute snap generating global array
 
@@ -118,12 +116,12 @@ thermo 		100
 # test output:   1: total potential energy
 #                2: xy component of stress tensor
 #                3: Sum(0.5*(B_{222}^i)^2, all i of type 2)
-#                4: xy component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
-#                5: z component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
+#                4: xz component of Sum(Sum(r_j*(0.5*(dB_{222}^i)^2/dR[j]), all i of type 2), all j)
+#                5: y component of -Sum(d(0.5*(B_{222}^i)^2/dR[2]), all i of type 2)
 #
 #                followed by 5 counterparts from compute snap
 
-thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[240]   v_db_2_120 		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40]
+thermo_style	custom 		pe            pxy            c_bsum2[20]   c_vbsum[220]   v_db_2_100 		c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40]
 thermo_modify 	norm no
 
 # dump 		mydump_db all custom 1000 dump_db id c_db[*]
@@ -165,13 +163,13 @@ Neighbor list info ...
       pair build: full/bin/atomonly
       stencil: full/bin/3d
       bin: standard
-WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:936)
-Per MPI rank memory allocation (min/avg/max) = 19.7 | 19.74 | 19.78 Mbytes
-PotEng Pxy c_bsum2[20] c_vbsum[240] v_db_2_120 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[13][40] c_snap[7][40] 
-   322.86952    1505558.1 4.2492771e+08    7860489.6    -17625699    322.86952    1505558.1 4.2492771e+08    7860489.6    -17625699 
-Loop time of 2.80142e-06 on 4 procs for 0 steps with 2 atoms
+WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964)
+Per MPI rank memory allocation (min/avg/max) = 21.43 | 21.47 | 21.52 Mbytes
+PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][41] c_snap[13][41] c_snap[1][40] c_snap[12][40] c_snap[6][40] 
+   322.86952    1505558.1 4.2492771e+08     -4952473     28484035    322.86952    1505558.1 4.2492771e+08     -4952473     28484035 
+Loop time of 2.25e-06 on 4 procs for 0 steps with 2 atoms
 
-107.1% CPU use with 4 MPI tasks x 1 OpenMP threads
+111.1% CPU use with 4 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
@@ -181,7 +179,7 @@ Neigh   | 0          | 0          | 0          |   0.0 |  0.00
 Comm    | 0          | 0          | 0          |   0.0 |  0.00
 Output  | 0          | 0          | 0          |   0.0 |  0.00
 Modify  | 0          | 0          | 0          |   0.0 |  0.00
-Other   |            | 2.801e-06  |            |       |100.00
+Other   |            | 2.25e-06   |            |       |100.00
 
 Nlocal:    0.5 ave 1 max 0 min
 Histogram: 2 0 0 0 0 0 0 0 0 2

From feec9673d892f1143fa115b0ec4620a85858718f Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Fri, 3 Jul 2020 14:03:00 -0600
Subject: [PATCH 14/18] Fixed html and spelling  warnings

---
 doc/src/Commands_compute.rst                | 1 +
 doc/src/compute.rst                         | 3 ++-
 doc/utils/sphinx-config/false_positives.txt | 1 +
 3 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst
index a54c4b4e03b..4d2ceabb4d2 100644
--- a/doc/src/Commands_compute.rst
+++ b/doc/src/Commands_compute.rst
@@ -80,6 +80,7 @@ KOKKOS, o = USER-OMP, t = OPT.
    * :doc:`ke/eff <compute_ke_eff>`
    * :doc:`ke/rigid <compute_ke_rigid>`
    * :doc:`mesont <compute_mesont>`
+   * :doc:`mliap <compute_mliap>`
    * :doc:`momentum <compute_momentum>`
    * :doc:`msd <compute_msd>`
    * :doc:`msd/chunk <compute_msd_chunk>`
diff --git a/doc/src/compute.rst b/doc/src/compute.rst
index 754cc3e1aa4..b74d2f70ba7 100644
--- a/doc/src/compute.rst
+++ b/doc/src/compute.rst
@@ -225,6 +225,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`ke/atom/eff <compute_ke_atom_eff>` - per-atom translational and radial kinetic energy in the electron force field model
 * :doc:`ke/eff <compute_ke_eff>` - kinetic energy of a group of nuclei and electrons in the electron force field model
 * :doc:`ke/rigid <compute_ke_rigid>` - translational kinetic energy of rigid bodies
+* :doc:`mliap <compute_mliap>` - gradients of energy and forces w.r.t. model parameters and related quantities for training machine learning interatomic potentials
 * :doc:`momentum <compute_momentum>` - translational momentum
 * :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms
 * :doc:`msd/chunk <compute_msd_chunk>` - mean-squared displacement for each chunk
@@ -273,7 +274,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`smd/ulsph/strain/rate <compute_smd_ulsph_strain_rate>` -
 * :doc:`smd/ulsph/stress <compute_smd_ulsph_stress>` - per-particle Cauchy stress tensor and von Mises equivalent stress in Smooth Mach Dynamics
 * :doc:`smd/vol <compute_smd_vol>` - per-particle volumes and their sum in Smooth Mach Dynamics
-* :doc:`snap <compute_sna_atom>` - bispectrum components and related quantities for a group of atoms
+* :doc:`snap <compute_sna_atom>` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials
 * :doc:`sna/atom <compute_sna_atom>` - bispectrum components for each atom
 * :doc:`snad/atom <compute_sna_atom>` - derivative of bispectrum components for each atom
 * :doc:`snav/atom <compute_sna_atom>` - virial contribution from bispectrum components for each atom
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 2bc3c2a0f68..9f188bd2f8d 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -2081,6 +2081,7 @@ Novint
 np
 Npair
 Npairs
+nparams
 nparticle
 npernode
 nph

From 7fe2df423f6214ae5143804b1713e52317128ee6 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Fri, 3 Jul 2020 14:22:43 -0600
Subject: [PATCH 15/18] Added a package README

---
 src/MLIAP/README | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 src/MLIAP/README

diff --git a/src/MLIAP/README b/src/MLIAP/README
new file mode 100644
index 00000000000..77e38d3c622
--- /dev/null
+++ b/src/MLIAP/README
@@ -0,0 +1,30 @@
+This package provides a general interface to families of 
+machine-learning interatomic potentials (MLIAPs). This interface consists
+of a mliap pair style and a mliap compute. 
+
+The mliap pair style is used when running simulations with energies and 
+forces calculated by an MLIAP. The interface allows separate 
+definitions of the interatomic potential functional form (*model*)
+and the geometric quantities that characterize the atomic positions
+(*descriptor*). By defining *model* and *descriptor* separately, 
+it is possible to use many different models with a given descriptor,
+or many different descriptors with a given model. Currently, the pair_style 
+supports just two models, *linear* and *quadratic*,
+and one descriptor, *sna*, the SNAP descriptor, including the 
+linear, quadratic, and chem variants. Work is currently underway to extend
+the interface to handle neural network energy models,
+and it is also straightforward to add new descriptor styles.
+
+The mliap compute style provides gradients of the energy, force,
+and stress tensor w.r.t. model parameters.
+These are useful when training MLIAPs to match target data.
+Any *model or *descriptor* that has been implemented for the
+*mliap* pair style can also be accessed by the *mliap* compute.
+In addition to the energy, force, and stress gradients, w.r.t.
+each *model* parameter, the compute also calculates the energy,
+force, and stress contributions from a user-specified 
+reference potential. 
+
+To see how this command
+can be used within a Python workflow to train machine-learning interatomic
+potentials, see the examples in FitSNAP https://github.com/FitSNAP/FitSNAP>.

From 37b2778d4b377d64cbc885919f4c19851536dfdc Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sat, 4 Jul 2020 13:34:35 -0400
Subject: [PATCH 16/18] Remove debug output

---
 src/MLIAP/pair_mliap.cpp | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp
index 14d4cce7538..a42e7dc231b 100644
--- a/src/MLIAP/pair_mliap.cpp
+++ b/src/MLIAP/pair_mliap.cpp
@@ -310,8 +310,6 @@ void PairMLIAP::init_style()
 double PairMLIAP::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
-  printf("itype = %d jtype = %d map[i] = %d map[j] = %d cutsq = %g\n",
-         i,j,map[i],map[j],descriptor->cutsq[map[i]][map[j]]);
   return sqrt(descriptor->cutsq[map[i]][map[j]]);
 }
 

From 759797733ddc0d1361eb670ef781050009761779 Mon Sep 17 00:00:00 2001
From: Aidan Thompson <athomps@sandia.gov>
Date: Sun, 5 Jul 2020 18:58:32 -0600
Subject: [PATCH 17/18] Generalized the variable and function names

---
 src/MLIAP/compute_mliap.cpp         | 120 +++++++++++-------------
 src/MLIAP/compute_mliap.h           |   8 +-
 src/MLIAP/mliap_descriptor.h        |   8 +-
 src/MLIAP/mliap_descriptor_snap.cpp | 138 ++++++++++++++++++++++++----
 src/MLIAP/mliap_descriptor_snap.h   |   8 +-
 src/MLIAP/pair_mliap.cpp            |   4 +-
 src/SNAP/compute_snap.cpp           |  19 +---
 7 files changed, 197 insertions(+), 108 deletions(-)

diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp
index 74e0a34ee5b..655e70c6f9d 100644
--- a/src/MLIAP/compute_mliap.cpp
+++ b/src/MLIAP/compute_mliap.cpp
@@ -35,7 +35,7 @@ enum{SCALAR,VECTOR,ARRAY};
 
 ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg), list(NULL), mliap(NULL),
-  mliap_peratom(NULL), mliapall(NULL), map(NULL), 
+  gradforce(NULL), mliapall(NULL), map(NULL), 
   descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL),
   gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL)
 {
@@ -90,18 +90,16 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
   nparams = model->nparams;
   nelements = model->nelements;
   gamma_nnz = model->get_gamma_nnz();
-  nperdim = nparams;
   ndims_force = 3;
   ndims_virial = 6;
-  yoffset = nperdim*atom->ntypes;
+  yoffset = nparams*nelements;
   zoffset = 2*yoffset;
   natoms = atom->natoms;
   size_array_rows = 1+ndims_force*natoms+ndims_virial;
-  size_array_cols = nperdim*atom->ntypes+1;
+  size_array_cols = nparams*nelements+1;
   lastcol = size_array_cols-1;
 
-  ndims_peratom = ndims_force;
-  size_peratom = ndims_peratom*nperdim*atom->ntypes;
+  size_gradforce = ndims_force*nparams*nelements;
 
   nmax = 0;
   gamma_max = 0;
@@ -121,7 +119,7 @@ ComputeMLIAP::~ComputeMLIAP()
 {
   memory->destroy(mliap);
   memory->destroy(mliapall);
-  memory->destroy(mliap_peratom);
+  memory->destroy(gradforce);
 
   memory->destroy(map);
 
@@ -186,7 +184,7 @@ void ComputeMLIAP::init()
 
   // find compute for reference energy
 
-  char *id_pe = (char *) "thermo_pe";
+  std::string id_pe = std::string("thermo_pe");
   int ipe = modify->find_compute(id_pe);
   if (ipe == -1)
     error->all(FLERR,"compute thermo_pe does not exist.");
@@ -194,15 +192,9 @@ void ComputeMLIAP::init()
 
   // add compute for reference virial tensor
 
-  char *id_virial = (char *) "mliap_press";
-  char **newarg = new char*[5];
-  newarg[0] = id_virial;
-  newarg[1] = (char *) "all";
-  newarg[2] = (char *) "pressure";
-  newarg[3] = (char *) "NULL";
-  newarg[4] = (char *) "virial";
-  modify->add_compute(5,newarg);
-  delete [] newarg;
+  std::string id_virial = std::string("mliap_press");
+  std::string pcmd = id_virial + " all pressure NULL virial";
+  modify->add_compute(pcmd);
 
   int ivirial = modify->find_compute(id_virial);
   if (ivirial == -1)
@@ -227,28 +219,28 @@ void ComputeMLIAP::compute_array()
 
   invoked_array = update->ntimestep;
 
-  // grow mliap_peratom array if necessary
+  // grow gradforce array if necessary
 
   if (atom->nmax > nmax) {
-    memory->destroy(mliap_peratom);
+    memory->destroy(gradforce);
     nmax = atom->nmax;
-    memory->create(mliap_peratom,nmax,size_peratom,
-                   "mliap:mliap_peratom");
+    memory->create(gradforce,nmax,size_gradforce,
+                   "mliap:gradforce");
   }
 
-  // clear global array
-
-  for (int irow = 0; irow < size_array_rows; irow++)
-    for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
-      mliap[irow][icoeff] = 0.0;
-
-  // clear local peratom array
+  // clear gradforce array
 
   for (int i = 0; i < ntotal; i++)
-    for (int icoeff = 0; icoeff < size_peratom; icoeff++) {
-      mliap_peratom[i][icoeff] = 0.0;
+    for (int j = 0; j < size_gradforce; j++) {
+      gradforce[i][j] = 0.0;
     }
 
+  // clear global array
+
+  for (int irow = 0; irow < size_array_rows; irow++)
+    for (int jcol = 0; jcol < size_array_cols; jcol++)
+      mliap[irow][jcol] = 0.0;
+
   // invoke full neighbor list (will copy or build if necessary)
 
   neighbor->build_one(list);
@@ -261,9 +253,9 @@ void ComputeMLIAP::compute_array()
     gamma_max = list->inum;
   }
 
-  // compute descriptors, if needed
+  // compute descriptors
 
-  descriptor->forward(map, list, descriptors);
+  descriptor->compute_descriptors(map, list, descriptors);
 
   // calculate descriptor contributions to parameter gradients
   // and gamma = double gradient w.r.t. parameters and descriptors
@@ -279,29 +271,28 @@ void ComputeMLIAP::compute_array()
 
   // calculate descriptor gradient contributions to parameter gradients
 
-  descriptor->param_backward(map, list, gamma_nnz, gamma_row_index, 
-                             gamma_col_index, gamma, mliap_peratom,
+  descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index, 
+                             gamma_col_index, gamma, gradforce,
                              yoffset, zoffset);
 
   // accumulate descriptor gradient contributions to global array
 
-  for (int itype = 0; itype < atom->ntypes; itype++) {
-    const int typeoffset_local = nperdim*itype;
-    const int typeoffset_global = nperdim*itype;
-    for (int icoeff = 0; icoeff < nperdim; icoeff++) {
+  for (int ielem = 0; ielem < nelements; ielem++) {
+    const int elemoffset = nparams*ielem;
+    for (int jparam = 0; jparam < nparams; jparam++) {
       int irow = 1;
       for (int i = 0; i < ntotal; i++) {
-        double *snadi = mliap_peratom[i]+typeoffset_local;
+        double *gradforcei = gradforce[i]+elemoffset;
         int iglobal = atom->tag[i];
         int irow = 3*(iglobal-1)+1;
-        mliap[irow][icoeff+typeoffset_global] += snadi[icoeff];
-        mliap[irow+1][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
-        mliap[irow+2][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
+        mliap[irow][jparam+elemoffset] += gradforcei[jparam];
+        mliap[irow+1][jparam+elemoffset] += gradforcei[jparam+yoffset];
+        mliap[irow+2][jparam+elemoffset] += gradforcei[jparam+zoffset];
       }
     }
   }
 
- // accumulate forces to global array
+ // copy forces to global array
 
   for (int i = 0; i < atom->nlocal; i++) {
     int iglobal = atom->tag[i];
@@ -315,25 +306,25 @@ void ComputeMLIAP::compute_array()
 
   dbdotr_compute();
 
-  // copy descriptor gradient contributions to global array
+  // copy energy gradient contributions to global array
 
-  for (int itype = 0; itype < atom->ntypes; itype++) {
-    const int typeoffset_global = nperdim*itype;
-    for (int icoeff = 0; icoeff < nperdim; icoeff++)
-      mliap[0][icoeff+typeoffset_global] = egradient[icoeff+typeoffset_global];
+  for (int ielem = 0; ielem < nelements; ielem++) {
+    const int elemoffset = nparams*ielem;
+    for (int jparam = 0; jparam < nparams; jparam++)
+      mliap[0][jparam+elemoffset] = egradient[jparam+elemoffset];
   }
 
   // sum up over all processes
 
   MPI_Allreduce(&mliap[0][0],&mliapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
 
-  // assign energy to last column
+  // copy energy to last column
 
   int irow = 0;
   double reference_energy = c_pe->compute_scalar();
   mliapall[irow++][lastcol] = reference_energy;
 
-  // assign virial stress to last column
+  // copy virial stress to last column
   // switch to Voigt notation
 
   c_virial->compute_vector();
@@ -362,21 +353,20 @@ void ComputeMLIAP::dbdotr_compute()
 
   int nall = atom->nlocal + atom->nghost;
   for (int i = 0; i < nall; i++)
-    for (int itype = 0; itype < atom->ntypes; itype++) {
-      const int typeoffset_local = nperdim*itype;
-      const int typeoffset_global = nperdim*itype;
-      double *snadi = mliap_peratom[i]+typeoffset_local;
-      for (int icoeff = 0; icoeff < nperdim; icoeff++) {
-        double dbdx = snadi[icoeff];
-        double dbdy = snadi[icoeff+yoffset];
-        double dbdz = snadi[icoeff+zoffset];
+    for (int ielem = 0; ielem < nelements; ielem++) {
+      const int elemoffset = nparams*ielem;
+      double *gradforcei = gradforce[i]+elemoffset;
+      for (int jparam = 0; jparam < nparams; jparam++) {
+        double dbdx = gradforcei[jparam];
+        double dbdy = gradforcei[jparam+yoffset];
+        double dbdz = gradforcei[jparam+zoffset];
         int irow = irow0;
-        mliap[irow++][icoeff+typeoffset_global] += dbdx*x[i][0];
-        mliap[irow++][icoeff+typeoffset_global] += dbdy*x[i][1];
-        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][2];
-        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][1];
-        mliap[irow++][icoeff+typeoffset_global] += dbdz*x[i][0];
-        mliap[irow++][icoeff+typeoffset_global] += dbdy*x[i][0];
+        mliap[irow++][jparam+elemoffset] += dbdx*x[i][0];
+        mliap[irow++][jparam+elemoffset] += dbdy*x[i][1];
+        mliap[irow++][jparam+elemoffset] += dbdz*x[i][2];
+        mliap[irow++][jparam+elemoffset] += dbdz*x[i][1];
+        mliap[irow++][jparam+elemoffset] += dbdz*x[i][0];
+        mliap[irow++][jparam+elemoffset] += dbdy*x[i][0];
       }
     }
 }
@@ -392,7 +382,7 @@ double ComputeMLIAP::memory_usage()
     sizeof(double);                                     // mliap
   bytes += size_array_rows*size_array_cols *
     sizeof(double);                                     // mliapall
-  bytes += nmax*size_peratom * sizeof(double);          // mliap_peratom
+  bytes += nmax*size_gradforce * sizeof(double);          // gradforce
   int n = atom->ntypes+1;
   bytes += n*sizeof(int);        // map
 
diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h
index 14f4a4bdaed..74b5fc1f93c 100644
--- a/src/MLIAP/compute_mliap.h
+++ b/src/MLIAP/compute_mliap.h
@@ -34,12 +34,12 @@ class ComputeMLIAP : public Compute {
   double memory_usage();
 
  private:
-  int natoms, nmax, size_peratom, lastcol;
-  int nperdim, yoffset, zoffset;
-  int ndims_peratom, ndims_force, ndims_virial;
+  int natoms, nmax, size_gradforce, lastcol;
+  int yoffset, zoffset;
+  int ndims_force, ndims_virial;
   class NeighList *list;
   double **mliap, **mliapall;
-  double **mliap_peratom;
+  double **gradforce;
   int *map;  // map types to [0,nelements)
   int nelements;
 
diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h
index f7f0db250db..072a5b4729c 100644
--- a/src/MLIAP/mliap_descriptor.h
+++ b/src/MLIAP/mliap_descriptor.h
@@ -22,9 +22,11 @@ class MLIAPDescriptor : protected Pointers  {
 public:
   MLIAPDescriptor(LAMMPS*);
   ~MLIAPDescriptor();
-  virtual void forward(int*, class NeighList*, double**)=0;
-  virtual void backward(class PairMLIAP*, class NeighList*, double**, int)=0;
-  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 
+  virtual void compute_descriptors(int*, class NeighList*, double**)=0;
+  virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int)=0;
+  virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, 
+                              double**, int, int)=0;
+  virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, 
                               double**, int, int)=0;
   virtual void init()=0;
   virtual double memory_usage()=0;
diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp
index b1542f25d45..c2879e942a7 100644
--- a/src/MLIAP/mliap_descriptor_snap.cpp
+++ b/src/MLIAP/mliap_descriptor_snap.cpp
@@ -76,7 +76,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()
    compute descriptors for each atom
    ---------------------------------------------------------------------- */
 
-void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptors)
+void MLIAPDescriptorSNAP::compute_descriptors(int* map, NeighList* list, double **descriptors)
 {
   int i,j,jnum,ninside;
   double delx,dely,delz,rsq;
@@ -151,7 +151,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
    compute forces for each atom
    ---------------------------------------------------------------------- */
 
-void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double **beta, int vflag)
+void MLIAPDescriptorSNAP::compute_forces(PairMLIAP* pairmliap, NeighList* list, double **beta, int vflag)
 {
   int i,j,jnum,ninside;
   double delx,dely,delz,rsq;
@@ -256,12 +256,12 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
 }
 
 /* ----------------------------------------------------------------------
-   compute forces for each atom
+   compute force gradient for each atom
    ---------------------------------------------------------------------- */
 
-void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, 
+void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, 
                                          int gamma_nnz, int **gamma_row_index, 
-                                         int **gamma_col_index, double **gamma, double **snadi,
+                                         int **gamma_col_index, double **gamma, double **gradforce,
                                          int yoffset, int zoffset)
 {
   int i,j,jnum,ninside;
@@ -345,20 +345,126 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
                                snaptr->rcutij[jj],jj, 0);
 
       snaptr->compute_dbidrj();
-
+      
       // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj
+      
+      for (int inz = 0; inz < gamma_nnz; inz++) {
+        const int l = gamma_row_index[ii][inz];
+        const int k = gamma_col_index[ii][inz];
+        gradforce[i][l]         += gamma[ii][inz]*snaptr->dblist[k][0];
+        gradforce[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1];
+        gradforce[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2];
+        gradforce[j][l]         -= gamma[ii][inz]*snaptr->dblist[k][0];
+        gradforce[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1];
+        gradforce[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2];
+      }
+      
+    }
+  }
 
-        for (int inz = 0; inz < gamma_nnz; inz++) {
-          const int l = gamma_row_index[ii][inz];
-          const int k = gamma_col_index[ii][inz];
-          snadi[i][l]         += gamma[ii][inz]*snaptr->dblist[k][0];
-          snadi[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1];
-          snadi[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2];
-          snadi[j][l]         -= gamma[ii][inz]*snaptr->dblist[k][0];
-          snadi[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1];
-          snadi[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2];
-        }
+}
+
+/* ----------------------------------------------------------------------
+   compute descriptor gradients for each neighbor atom
+   ---------------------------------------------------------------------- */
+
+void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list, 
+                                         int gamma_nnz, int **gamma_row_index, 
+                                         int **gamma_col_index, double **gamma, double **graddesc,
+                                         int yoffset, int zoffset)
+{
+  int i,j,jnum,ninside;
+  double delx,dely,delz,evdwl,rsq;
+  double fij[3];
+  int *jlist,*numneigh,**firstneigh;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  for (int ii = 0; ii < list->inum; ii++) {
+    i = list->ilist[ii];
+
+    const double xtmp = x[i][0];
+    const double ytmp = x[i][1];
+    const double ztmp = x[i][2];
+    const int itype = type[i];
+    const int ielem = map[itype];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    // insure rij, inside, wj, and rcutij are of size jnum
+
+    snaptr->grow_rij(jnum);
+
+    // rij[][3] = displacements between atom I and those neighbors
+    // inside = indices of neighbors of I within cutoff
+    // wj = weights for neighbors of I within cutoff
+    // rcutij = cutoffs for neighbors of I within cutoff
+    // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
+
+    ninside = 0;
+    for (int jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      delx = x[j][0] - xtmp;
+      dely = x[j][1] - ytmp;
+      delz = x[j][2] - ztmp;
+      rsq = delx*delx + dely*dely + delz*delz;
+      int jtype = type[j];
+      const int jelem = map[jtype];
+
+      if (rsq < cutsq[ielem][jelem]) {
+        snaptr->rij[ninside][0] = delx;
+        snaptr->rij[ninside][1] = dely;
+        snaptr->rij[ninside][2] = delz;
+        snaptr->inside[ninside] = j;
+	snaptr->wj[ninside] = wjelem[jelem];
+	snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
+        snaptr->element[ninside] = jelem; // element index for chem snap
+        ninside++;
+      }
+    }
+
+    if (chemflag)
+      snaptr->compute_ui(ninside, ielem);
+    else
+      snaptr->compute_ui(ninside, 0);
+
+    snaptr->compute_zi();
+    if (chemflag)
+      snaptr->compute_bi(ielem);
+    else
+      snaptr->compute_bi(0);
+
+    for (int jj = 0; jj < ninside; jj++) {
+      const int j = snaptr->inside[jj];
+
+      if(chemflag)
+        snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+      else
+        snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
+                               snaptr->rcutij[jj],jj, 0);
+
+      snaptr->compute_dbidrj();
+
+      // Accumulate dB_k^i/dRi, dB_k^i/dRj
 
+      for (int k = 0; k < ndescriptors; k++) {
+        graddesc[i][k] = snaptr->dblist[k][0];
+        graddesc[i][k] = snaptr->dblist[k][1];
+        graddesc[i][k] = snaptr->dblist[k][2];
+        graddesc[j][k] = -snaptr->dblist[k][0];
+        graddesc[j][k] = -snaptr->dblist[k][1];
+        graddesc[j][k] = -snaptr->dblist[k][2];
+      } 
     }
   }
 
diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h
index 4843fbbcc80..a929e011493 100644
--- a/src/MLIAP/mliap_descriptor_snap.h
+++ b/src/MLIAP/mliap_descriptor_snap.h
@@ -22,9 +22,11 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor  {
 public:
   MLIAPDescriptorSNAP(LAMMPS*, char*);
   ~MLIAPDescriptorSNAP();
-  virtual void forward(int*, class NeighList*, double**);
-  virtual void backward(class PairMLIAP*, class NeighList*, double**, int);
-  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 
+  virtual void compute_descriptors(int*, class NeighList*, double**);
+  virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int);
+  virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, 
+                              double**, int, int);
+  virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, 
                               double**, int, int);
   virtual void init();
   virtual double memory_usage();
diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp
index a42e7dc231b..2dd70fc6c23 100644
--- a/src/MLIAP/pair_mliap.cpp
+++ b/src/MLIAP/pair_mliap.cpp
@@ -86,7 +86,7 @@ void PairMLIAP::compute(int eflag, int vflag)
   // compute descriptors, if needed
 
   if (model->nonlinearflag || eflag)
-    descriptor->forward(map, list, descriptors);
+    descriptor->compute_descriptors(map, list, descriptors);
 
   // compute E_i and beta_i = dE_i/dB_i for all i in list
 
@@ -94,7 +94,7 @@ void PairMLIAP::compute(int eflag, int vflag)
 
   // calculate force contributions beta_i*dB_i/dR_j
 
-  descriptor->backward(this, list, beta, vflag);
+  descriptor->compute_forces(this, list, beta, vflag);
 
   // calculate stress
 
diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp
index c9e10065710..613fc1a862d 100644
--- a/src/SNAP/compute_snap.cpp
+++ b/src/SNAP/compute_snap.cpp
@@ -210,14 +210,9 @@ void ComputeSnap::init()
                  "snap:snapall");
   array = snapall;
 
-  // INCOMPLETE: modify->find_compute()
-  // was called 223960 times by snappy Ta example
-  // that is over 600 times per config?
-  // how is this possible???
-
   // find compute for reference energy
 
-  char *id_pe = (char *) "thermo_pe";
+  std::string id_pe = std::string("thermo_pe");
   int ipe = modify->find_compute(id_pe);
   if (ipe == -1)
     error->all(FLERR,"compute thermo_pe does not exist.");
@@ -225,15 +220,9 @@ void ComputeSnap::init()
 
   // add compute for reference virial tensor
 
-  char *id_virial = (char *) "snap_press";
-  char **newarg = new char*[5];
-  newarg[0] = id_virial;
-  newarg[1] = (char *) "all";
-  newarg[2] = (char *) "pressure";
-  newarg[3] = (char *) "NULL";
-  newarg[4] = (char *) "virial";
-  modify->add_compute(5,newarg);
-  delete [] newarg;
+  std::string id_virial = std::string("snap_press");
+  std::string pcmd = id_virial + " all pressure NULL virial";
+  modify->add_compute(pcmd);
 
   int ivirial = modify->find_compute(id_virial);
   if (ivirial == -1)

From fe12ea273457ca3ca5daeaa992589bc60826942d Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 5 Jul 2020 22:46:48 -0400
Subject: [PATCH 18/18] simplify compute creation

---
 src/SNAP/compute_snap.cpp | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp
index 613fc1a862d..7eb216a63d4 100644
--- a/src/SNAP/compute_snap.cpp
+++ b/src/SNAP/compute_snap.cpp
@@ -212,19 +212,16 @@ void ComputeSnap::init()
 
   // find compute for reference energy
 
-  std::string id_pe = std::string("thermo_pe");
-  int ipe = modify->find_compute(id_pe);
+  int ipe = modify->find_compute("thermo_pe");
   if (ipe == -1)
     error->all(FLERR,"compute thermo_pe does not exist.");
   c_pe = modify->compute[ipe];
 
   // add compute for reference virial tensor
 
-  std::string id_virial = std::string("snap_press");
-  std::string pcmd = id_virial + " all pressure NULL virial";
-  modify->add_compute(pcmd);
+  modify->add_compute("snap_press all pressure NULL virial");
 
-  int ivirial = modify->find_compute(id_virial);
+  int ivirial = modify->find_compute("snap_press");
   if (ivirial == -1)
     error->all(FLERR,"compute snap_press does not exist.");
   c_virial = modify->compute[ivirial];