Axels requested revisions

src/USER-RANN/pair_rann.h

@@ -0,0 +1,174 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#ifdef PAIR_CLASS
+
+PairStyle(rann,PairRANN)
+
+#else
+
+#ifndef LMP_PAIR_RANN
+#define LMP_PAIR_RANN
+
+
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "pair.h"
+#include <map>
+
+
+namespace LAMMPS_NS {
+
+  namespace RANN {
+    //forward declarations
+    class Activation;
+    class Fingerprint;
+  }
+
+  class PairRANN : public Pair {
+   public:
+    //inherited functions
+    PairRANN(class LAMMPS *);
+    ~PairRANN();
+    void compute(int, int);
+    void settings(int, char **);
+    void coeff(int, char **);
+    void init_style();
+    double init_one(int, int);
+    void init_list(int , NeighList *);
+    void errorf(const char*,int,const char*);
+    int factorial(int);
+    //black magic for modular fingerprints and activations
+    class RANN::Activation ***activation;
+    class RANN::Fingerprint ***fingerprints;
+    typedef RANN::Fingerprint *(*FingerprintCreator)(PairRANN *);
+    typedef RANN::Activation *(*ActivationCreator)(PairRANN *);
+    typedef std::map<std::string,FingerprintCreator> FingerprintCreatorMap;
+    typedef std::map<std::string,ActivationCreator> ActivationCreatorMap;
+    FingerprintCreatorMap *fingerprint_map;
+    ActivationCreatorMap *activation_map;
+    RANN::Fingerprint * create_fingerprint(const char *);
+    RANN::Activation * create_activation(const char *);
+
+    //global variables
+    int nelements;                // # of elements (distinct from LAMMPS atom types since multiple atom types can be mapped to one element)
+    int nelementsp;        // nelements+1
+    char **elements;              // names of elements
+    char **elementsp;        // names of elements with "all" appended as the last "element"
+    double *mass;                 // mass of each element
+    double cutmax;        // max radial distance for neighbor lists
+    int *map;                     // mapping from atom types to elements
+    int *fingerprintcount;    // static variable used in initialization
+    int *fingerprintlength;       // # of input neurons defined by fingerprints of each element.
+    int *fingerprintperelement;   // # of fingerprints for each element
+    bool doscreen;//screening is calculated if any defined fingerprint uses it
+    bool allscreen;//all fingerprints use screening so screened neighbors can be completely ignored
+    bool dospin;
+    int res;//Resolution of function tables for cubic interpolation.
+    int memguess;
+    double *screening_min;
+    double *screening_max;
+    bool **weightdefined;
+    bool **biasdefined;
+    int nmax1;
+    int nmax2;
+    int fmax;
+    int fnmax;
+    //memory actively written to during each compute:
+    double *xn,*yn,*zn,*Sik,*dSikx,*dSiky,*dSikz,*dSijkx,*dSijky,*dSijkz,*sx,*sy,*sz,**dSijkxc,**dSijkyc,**dSijkzc,*dfeaturesx,*dfeaturesy,*dfeaturesz,*features;
+    double *layer,*sum,*sum1,**dlayerx,**dlayery,**dlayerz,**dlayersumx,**dlayersumy,**dlayersumz;
+    double **dsx,**dsy,**dsz,**dssumx,**dssumy,**dssumz;
+    int *tn,*jl;
+    bool *Bij;
+
+    struct Simulation{
+      int *id;
+    bool forces;
+    bool spins;
+    double **x;
+    double **f;
+    double **s;
+    double box[3][3];
+    double origin[3];
+    double **features;
+    double **dfx;
+    double **dfy;
+    double **dfz;
+    double **dsx;
+    double **dsy;
+    double **dsz;
+    int *ilist,*numneigh,**firstneigh,*type,inum,gnum;
+    };
+
+    struct NNarchitecture{
+    int layers;
+    int *dimensions;//vector of length layers with entries for neurons per layer
+    double **Weights;
+    double **Biases;
+    int *activations;//unused
+    int maxlayer;//longest layer (for memory allocation)
+    };
+
+    Simulation *sims;
+    NNarchitecture *net;//array of networks, 1 for each element.
+
+   private:
+    template <typename T> static RANN::Fingerprint *fingerprint_creator(PairRANN *);
+    template <typename T> static RANN::Activation *activation_creator(PairRANN *);
+    //new functions
+    void allocate(std::vector<std::string>);//called after reading element list, but before reading the rest of the potential
+    void read_file(char *);//read potential file
+    void strip_comments(char *,FILE *fp);
+    void read_atom_types(std::vector<std::string>,std::vector<std::string>);
+    void read_mass(std::vector<std::string>,std::vector<std::string>);
+    void read_fpe(std::vector<std::string>,std::vector<std::string>);//fingerprints per element. Count total fingerprints defined for each 1st element in element combinations
+    void read_fingerprints(std::vector<std::string>,std::vector<std::string>);
+    void read_fingerprint_constants(std::vector<std::string>,std::vector<std::string>);
+    void read_network_layers(std::vector<std::string>,std::vector<std::string>);//include input and output layer (hidden layers + 2)
+    void read_layer_size(std::vector<std::string>,std::vector<std::string>);
+    void read_weight(std::vector<std::string>,std::vector<std::string>,FILE*);//weights should be formatted as properly shaped matrices
+    void read_bias(std::vector<std::string>,std::vector<std::string>,FILE*);//biases should be formatted as properly shaped vectors
+    void read_activation_functions(std::vector<std::string>,std::vector<std::string>);
+    void read_screening(std::vector<std::string>,std::vector<std::string>);
+    bool check_potential();//after finishing reading potential file
+    void propagateforward(double *,double **,double **,int,int);//called by compute to get force and energy
+    void propagateforwardspin(double *,double **,double **,double**,int,int);//called by compute to get force and energy
+    void screen(int,int,int);
+    void cull_neighbor_list(int *,int,int);
+    void screen_neighbor_list(int *,int,int);
+  };
+
+}
+
+#endif
+#endif
+
+
+

src/USER-RANN/rann_activation.h

@@ -0,0 +1,73 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#ifndef LMP_RANN_ACTIVATION_H
+#define LMP_RANN_ACTIVATION_H
+
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Activation {
+    public:
+    Activation(class PairRANN *);
+    virtual ~Activation();
+    virtual double activation_function(double);
+    virtual double dactivation_function(double);
+    virtual double ddactivation_function(double);
+    bool empty;
+    const char *style;
+    };
+
+    Activation::Activation(PairRANN *_pair) {
+    empty = true;
+    style = "empty";
+    }
+
+    Activation::~Activation() {
+
+    }
+
+    //default is linear activation (no change).
+    double Activation::activation_function(double A) {
+    return A;
+    }
+
+    double Activation::dactivation_function(double A) {
+    return 1.0;
+    }
+
+    double Activation::ddactivation_function(double A) {
+    return 0.0;
+    }
+  }
+}
+
+
+#endif /* RANN_ACTIVATION_H_ */

src/USER-RANN/rann_activation_linear.h

@@ -0,0 +1,77 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef ACTIVATION_CLASS
+
+ActivationStyle(linear,Activation_linear)
+
+#else
+
+#ifndef LMP_RANN_ACTIVATION_LINEAR_H
+#define LMP_RANN_ACTIVATION_LINEAR_H
+
+#include "rann_activation.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Activation_linear : public Activation {
+    public:
+      Activation_linear(class PairRANN *);
+      double activation_function(double);
+      double dactivation_function(double);
+      double ddactivation_function(double);
+    };
+
+    Activation_linear::Activation_linear(PairRANN *_pair) : Activation(_pair) {
+      empty = false;
+      style = "linear";
+    }
+
+    double Activation_linear::activation_function(double A)
+    {
+      return A;
+    }
+
+    double Activation_linear::dactivation_function(double A)
+    {
+      return 1.0;
+    }
+
+    double Activation_linear::ddactivation_function(double) {
+      return 0.0;
+    }
+
+
+  }
+}
+
+#endif
+#endif /* ACTIVATION_LINEAR_H_ */

src/USER-RANN/rann_activation_sig_i.h

@@ -0,0 +1,77 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#ifdef ACTIVATION_CLASS
+
+ActivationStyle(sigI,Activation_sigI)
+
+#else
+
+#ifndef LMP_RANN_ACTIVATION_SIGI_H
+#define LMP_RANN_ACTIVATION_SIGI_H
+
+#include "rann_activation.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Activation_sigI : public Activation {
+    public:
+    Activation_sigI(class PairRANN *);
+    double activation_function(double);
+    double dactivation_function(double);
+    double ddactivation_function(double);
+    };
+
+    Activation_sigI::Activation_sigI(PairRANN *_pair) : Activation(_pair) {
+    empty = false;
+    style = "sigI";
+    }
+
+    double Activation_sigI::activation_function(double in) {
+    if (in>34)return in;
+    return 0.1*in + 0.9*log(exp(in) + 1);
+    }
+
+    double Activation_sigI::dactivation_function(double in) {
+    if (in>34)return 1;
+    return 0.1 + 0.9/(exp(in)+1)*exp(in);
+    }
+
+    double Activation_sigI::ddactivation_function(double in) {
+    if (in>34)return 0;
+    return 0.9*exp(in)/(exp(in)+1)/(exp(in)+1);;
+    }
+
+  }
+}
+
+#endif
+#endif /* ACTIVATION_SIGI_H_ */

src/USER-RANN/rann_fingerprint.cpp

@@ -0,0 +1,116 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#include "rann_fingerprint.h"
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint::Fingerprint(PairRANN *_pair)
+{
+  spin = false;
+  screen = false;
+  empty = true;
+  fullydefined = false;
+  n_body_type = 0;
+  style = "empty";
+  pair = _pair;
+}
+
+Fingerprint::~Fingerprint() {
+
+}
+
+bool Fingerprint::parse_values(std::string line,std::vector<std::string> line1) {
+  return false;
+}
+
+void Fingerprint::init(int *i,int id) {
+
+}
+
+void Fingerprint::allocate() {
+
+}
+
+void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+
+}
+
+void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+
+}
+
+void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *sx, double *sy, double *sz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+
+}
+
+void Fingerprint::compute_fingerprint(double *features,double *dfeaturesx,double *dfeaturesy,double * dfeaturesz,double *sx, double *sy, double *sz, double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+
+}
+
+void Fingerprint::write_values(FILE *fid) {
+
+}
+
+int Fingerprint::get_length() {
+  return 0;
+}
+
+//Smooth cutoff, goes from 1 to zero over the interval rc-dr to rc. Same as MEAM uses. Used by generateradialtable and generatexpcuttable.
+double Fingerprint::cutofffunction(double r,double rc, double dr) {
+  double out;
+  if (r < (rc -dr))out=1;
+  else if (r>rc)out=0;
+  else {
+	  out = 1-(rc-r)/dr;
+	  out *= out;
+	  out *= out;
+	  out = 1-out;
+	  out *= out;
+  }
+  return out;
+}
+
+void Fingerprint::generate_rinvssqrttable() {
+  int buf = 5;
+  int m;
+  double r1;
+  double cutmax = pair->cutmax;
+  int res = pair->res;
+  rinvsqrttable = new double[res+buf];
+  for (m=0;m<(res+buf);m++) {
+    r1 = cutmax*cutmax*(double)(m)/(double)(res);
+    rinvsqrttable[m] = 1/sqrt(r1);
+  }
+}
+
+
+
+

src/USER-RANN/rann_fingerprint.h

@@ -0,0 +1,72 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+
+----------------*/
+
+#ifndef LMP_RANN_FINGERPRINT_H
+#define LMP_RANN_FINGERPRINT_H
+
+#include "pair_rann.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Fingerprint {
+      public:
+      Fingerprint(PairRANN *);
+      virtual ~Fingerprint();
+      virtual bool parse_values(std::string,std::vector<std::string>);
+      virtual void write_values(FILE *);
+      virtual void init(int*,int);
+      virtual void allocate();
+      void init_screen(int);
+      virtual void compute_fingerprint(double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//no screen,no spin
+      virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//screen
+      virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//spin
+      virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//spin,screen
+      virtual int get_length();
+      virtual double cutofffunction(double,double, double);
+      virtual void generate_rinvssqrttable();
+      bool spin;
+      bool screen;
+      int n_body_type;//i-j vs. i-j-k vs. i-j-k-l, etc.
+      bool empty;
+      bool fullydefined;
+      int startingneuron;
+      int id;//based on ordering of fingerprints listed for i-j in potential file
+      const char *style;
+      int* atomtypes;
+      double *rinvsqrttable;
+      double rc;
+      PairRANN *pair;
+      };
+  }
+}
+
+
+#endif /* RANN_FINGERPRINT_H_ */

src/USER-RANN/rann_fingerprint_bond.cpp

@@ -0,0 +1,808 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#include "rann_fingerprint_bond.h"
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_bond::Fingerprint_bond(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 3;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha_k = new double [1];
+  alpha_k[0] = -1;
+  kmax = 0;
+  mlength = 0;
+  id = -1;
+  style = "bond";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  _pair->allscreen = false;
+}
+
+Fingerprint_bond::~Fingerprint_bond() {
+  delete [] alpha_k;
+  delete [] atomtypes;
+  delete [] expcuttable;
+  delete [] dfctable;
+  for (int i=0;i<(mlength*(mlength+1))>>1;i++) {
+    delete [] coeff[i];
+    delete [] coeffx[i];
+    delete [] coeffy[i];
+    delete [] coeffz[i];
+    delete [] Mf[i];
+  }
+  delete [] coeff;
+  delete [] coeffx;
+  delete [] coeffy;
+  delete [] coeffz;
+  delete [] Mf;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_bond::parse_values(std::string constant,std::vector<std::string> line1) {
+  int nwords,l;
+  nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alphak")==0) {
+    delete [] alpha_k;
+    alpha_k = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha_k[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("k")==0) {
+    kmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("m")==0) {
+    mlength = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for bond power");
+  if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_bond::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alphak:\n",style,id);
+  for (i=0;i<kmax;i++) {
+    fprintf(fid,"%f ",alpha_k[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:k:\n",style,id);
+  fprintf(fid,"%d\n",kmax);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:m:\n",style,id);
+  fprintf(fid,"%d\n",mlength);
+}
+
+void Fingerprint_bond::init(int *i,int _id) {
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  re = 0;
+  rc = 0;
+  mlength = 0;
+  kmax = 0;
+  alpha_k = new double [1];
+  alpha_k[0]=-1;
+  empty = false;
+  id = _id;
+}
+
+//number of neurons defined by this fingerprint
+int Fingerprint_bond::get_length() {
+  return mlength*kmax;
+}
+
+void Fingerprint_bond::allocate() {
+  generate_exp_cut_table();
+  generate_coefficients();
+  generate_rinvssqrttable();
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
+void Fingerprint_bond::generate_exp_cut_table() {
+  int m,n;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  expcuttable = new double [(res+buf)*(kmax)];
+  dfctable = new double [res+buf];
+  for (m=0;m<(res+buf);m++) {
+    r1 = cutmax*cutmax*(double)(m)/(double)(res);
+    for (n=0;n<(kmax);n++) {
+      expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+      dfctable[m]=0;
+    }
+    else{
+      dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop.
+void Fingerprint_bond::generate_coefficients() {      //calculates multinomial coefficient for each term
+  int p,mb,mc;
+  int n,p1,i1;
+  mb = mlength;
+  mc=(mb*(mb+1))>>1;
+  coeff = new int *[mc];
+  coeffx = new int *[mc];
+  coeffy = new int *[mc];
+  coeffz = new int *[mc];
+  for (p=0;p<mc;p++) {
+    coeff[p]=new int [mb];
+    coeffx[p]=new int [mb];
+    coeffy[p]=new int [mb];
+    coeffz[p]=new int [mb];
+  }
+  Mf = new int*[mc];
+  int *M = new int[mlength+1];
+  for (p=0;p<mlength+1;p++) {
+    M[p]=0;
+  }
+  for (p1=0;p1<mc;p1++) {
+    Mf[p1] = new int[mlength+1];
+    for (p=0;p<mlength+1;p++) {
+      Mf[p1][p]=0;
+    }
+  }
+  M[0] = 2;
+  Mf[0][0] = 2;
+  n = 1;
+  int m1 = 1;
+  bool go = true;
+  bool broke = false;
+  while (go) {
+    broke = false;
+    for (i1=0;i1<mlength-1;i1++) {
+      if (M[i1+1] == 0) {
+        M[i1+1]=M[i1+1]+1;
+        for (p1=0;p1<mlength+1;p1++) {
+          Mf[n][p1] = M[p1];
+        }
+        n = n+1;
+        broke = true;
+        break;
+      }
+    }
+    if (m1<mlength && !broke) {
+      M[m1]=M[m1]+1;
+      for (p1=m1+1;p1<mlength+1;p1++) {
+        M[p1]=0;
+      }
+      for (p1=0;p1<mlength+1;p1++) {
+        Mf[n][p1]=M[p1];
+      }
+      n=n+1;
+      broke = true;
+      m1 = m1+1;
+    }
+    if (!broke) {
+      go = false;
+    }
+  }
+  for (p=0;p<mb;p++) {
+    for (p1=0;p1<mc;p1++) {
+      if (p==0) {
+        coeffx[p1][p]=0;
+        coeffy[p1][p]=0;
+        coeffz[p1][p]=0;
+      }
+      else{
+        coeffx[p1][p]=coeffx[p1][p-1];
+        if (Mf[p1][p]==0) {
+          coeffx[p1][p]++;
+        }
+        coeffy[p1][p]=coeffy[p1][p-1];
+        if (Mf[p1][p]==1) {
+          coeffy[p1][p]++;
+        }
+        coeffz[p1][p]=coeffz[p1][p-1];
+        if (Mf[p1][p]==2) {
+          coeffz[p1][p]++;
+        }
+      }
+      coeff[p1][p] = pair->factorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]);
+    }
+  }
+}
+
+
+//Called by getproperties. Gets 3-body features and dfeatures
+void Fingerprint_bond::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i;
+  int *ilist,*numneigh;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  ilist = sim->ilist;
+  numneigh = sim->numneigh;
+  i = ilist[ii];
+  //select the more efficient algorithm for this particular potential and environment.
+  if (jnum*2>(mlength+1)*mlength*20) {
+    do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,sid,xn,yn,zn,tn,jnum,jl);
+  }
+  else{
+    do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,sid,xn,yn,zn,tn,jnum,jl);
+
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers
+void Fingerprint_bond::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  double cutmax = pair->cutmax;
+  int res = pair->res;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  int nelements=pair->nelements;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax+12];
+  int p = kmax;
+  int countmb=((mlength)*(mlength+1))>>1;
+  // calculate interpolation expr, rinvs and dfc, for each neighbor
+  for (jj = 0; jj < jnum; jj++) {
+    jtype = tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) {
+      expr[jj][0]=0;
+      continue;
+    }
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    expr[jj][0]=0;
+    continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+      expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+    }
+    double* q = &dfctable[m1-1];
+    double* ri = &rinvsqrttable[m1-1];
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
+
+    expr[jj][p]=delx*rinvs;
+    expr[jj][p+1]=dely*rinvs;
+    expr[jj][p+2]=delz*rinvs;
+    //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
+    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+      expr[jj][p] = 0.000001;
+    }
+    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+      expr[jj][p+1] = 0.000001;
+    }
+    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+      expr[jj][p+2] = 0.000001;
+    }
+    expr[jj][p+3] = -dfc*expr[jj][p];
+    expr[jj][p+4] = rinvs/expr[jj][p];
+    expr[jj][p+5] = rinvs*expr[jj][p];
+    expr[jj][p+6] = -dfc*expr[jj][p+1];
+    expr[jj][p+7] = rinvs/expr[jj][p+1];
+    expr[jj][p+8] = rinvs*expr[jj][p+1];
+    expr[jj][p+9] = -dfc*expr[jj][p+2];
+    expr[jj][p+10] = rinvs/expr[jj][p+2];
+    expr[jj][p+11] = rinvs*expr[jj][p+2];
+  }
+
+  int kb = kmax;
+  int mb = mlength;
+  count = startingneuron;
+  double Bb[mb];
+  double dBbx;
+  double dBby;
+  double dBbz;
+  double yprod;
+  for (mcount=0;mcount<countmb;mcount++) {
+    int *M = Mf[mcount];
+    int *_coeffx = coeffx[mcount];
+    int *_coeffy = coeffy[mcount];
+    int *_coeffz = coeffz[mcount];
+    int *_coeff = coeff[mcount];
+    a = mb+1;
+    for (a1=0;a1<mb;a1++) {
+      if (Mf[mcount][a1+1]==0) {
+        a = a1;
+        break;
+      }
+    }
+    for (n=0;n<kb;n++) {
+      for (a1=0;a1<mb;a1++) {
+        Bb[a1]=0;
+      }
+      ai = n;
+      double y1 = alpha_k[ai]/re;
+      //loop over jtype to get Bb
+      for (jj=0;jj<jnum;jj++) {
+        if (expr[jj][0]==0) {continue;}
+        jtype = tn[jj];
+        if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+          continue;
+        }
+        double yprod = expr[jj][ai];
+        double *y4 = &expr[jj][p];
+        for (a2=0;a2<a;a2++) {
+          yprod *= y4[M[a2+1]];
+        }
+        for (a2=a;a2<mb;a2++) {
+          Bb[a2]=Bb[a2]+yprod;
+          yprod *= y4[M[a2+1]];
+        }
+      }
+      if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
+        double Bg[mb];
+        for (a1=0;a1<mb;a1++) {
+          Bg[a1]=0;
+        }
+        ai = n;
+        double y1 = alpha_k[ai]/re;
+        //loop over ktype to get Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          double yprod = expr[jj][ai];
+          double *y4 = &expr[jj][p];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          for (a2=a;a2<mb;a2++) {
+            Bg[a2]=Bg[a2]+yprod;
+            yprod *= y4[M[a2+1]];
+          }
+        }
+        double B1;
+        //loop over ktype to get dBg*Bb
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bb[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //loop over jtype to get dBb*Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bg[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bb[a2]*Bg[a2]*_coeff[a2];
+        }
+      }
+      else{//Bb=Bg
+        double B1;
+        //loop over jtype to get 2*Bb*dBb
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = 2*Bb[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bb[a2]*Bb[a2]*_coeff[a2];
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
+void Fingerprint_bond::do3bodyfeatureset_doubleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,int ii, int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,ktype,kk,m,n;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int jtypes = atomtypes[1];
+  int ktypes = atomtypes[2];
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  double **x = sim->x;
+  int *type = sim->type;
+  int nelements = pair->nelements;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax];
+  double y[jnum][3];
+  double ri[jnum];
+  double dfc[jnum];
+  int kb = kmax;
+  int mb = mlength;
+  double c41[kmax];
+  double c51[kmax];
+  double c61[kmax];
+  double ct[kmax];
+  for (jj = 0; jj < jnum; jj++) {
+    jtype = tn[jj];
+    if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
+      expr[jj][0]=0;
+      continue;
+    }
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    	expr[jj][0]=0;
+    	continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere.
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+    	expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+    }
+    double* q = &dfctable[m1-1];
+    double* r2 = &rinvsqrttable[m1-1];
+    dfc[jj] = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
+    y[jj][0]=delx*ri[jj];
+    y[jj][1]=dely*ri[jj];
+    y[jj][2]=delz*ri[jj];
+  }
+  for (jj = 0; jj < jnum; jj++) {
+    if (expr[jj][0]==0)continue;
+    jtype = tn[jj];
+    if (jtypes != nelements && jtypes != jtype) {
+      continue;
+    }
+    for (n = 0;n<kmax;n++) {
+      ct[n] = 2*alpha_k[n]/re;
+      c41[n]=(-ct[n]+2*dfc[jj])*y[jj][0];
+      c51[n]=(-ct[n]+2*dfc[jj])*y[jj][1];
+      c61[n]= (-ct[n]+2*dfc[jj])*y[jj][2];
+    }
+    if (jtypes==ktypes) {
+      for (kk=jj+1;kk< jnum; kk++) {
+        if (expr[kk][0]==0)continue;
+        ktype = tn[kk];
+        if (ktypes != nelements && ktypes != ktype) {
+          continue;
+        }
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1  = 2*ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2  = 2*ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3  = 2*ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        double c10 = 2*ri[kk]*(y[jj][0]-dot*y[kk][0]);
+        double c11 = 2*ri[kk]*(y[jj][1]-dot*y[kk][1]);
+        double c12 = 2*ri[kk]*(y[jj][2]-dot*y[kk][2]);
+        //alternate formulation:
+//        double c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
+//        double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
+//        double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
+//        double c10 = 2*ri[kk]*y[jj][0]*(1-y[kk][0]*y[kk][0]);
+//        double c11 = 2*ri[kk]*y[jj][1]*(1-y[kk][1]*y[kk][1]);
+//        double c12 = 2*ri[kk]*y[jj][2]*(1-y[kk][2]*y[kk][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n];
+          double c5 = c51[n];
+          double c6 = c61[n];
+          double ct2 = -ct[n]+2*dfc[kk];
+          double c42 = ct2*y[kk][0];
+          double c52 = ct2*y[kk][1];
+          double c62 = ct2*y[kk][2];
+          //m=0
+          features[count]+=2*dot1;
+          dfeaturesx[jj*f+count]+=dot1*c4;
+          dfeaturesy[jj*f+count]+=dot1*c5;
+          dfeaturesz[jj*f+count]+=dot1*c6;
+          dfeaturesx[kk*f+count]+=dot1*c42;
+          dfeaturesy[kk*f+count]+=dot1*c52;
+          dfeaturesz[kk*f+count]+=dot1*c62;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          c42*=dot;
+          c52*=dot;
+          c62*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7;
+            dfeaturesy[jj*f+count]+=c8;
+            dfeaturesz[jj*f+count]+=c9;
+            dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
+            dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
+            dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
+            dot1*=dot;
+            features[count++]+=2*dot1;
+          }
+        }
+      }
+      kk=jj;
+      if (ktypes == nelements || ktypes == jtype) {
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1 = 2*ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2 = 2*ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3 = 2*ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        //alternate formulation:
+//        double c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
+//        double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
+//        double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n];
+          double c5 = c51[n];
+          double c6 = c61[n];
+          //m=0
+          features[count]+=dot1;
+          dfeaturesx[jj*f+count]+=dot1*c4;
+          dfeaturesy[jj*f+count]+=dot1*c5;
+          dfeaturesz[jj*f+count]+=dot1*c6;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7;
+            dfeaturesy[jj*f+count]+=c8;
+            dfeaturesz[jj*f+count]+=c9;
+            dot1*=dot;
+            features[count++]+=dot1;
+          }
+        }
+      }
+    }
+    else {
+      for (kk=0;kk<jnum; kk++) {
+        if (expr[kk][0]==0)continue;
+        ktype = tn[kk];
+        if (ktypes != nelements && ktypes != ktype) {
+          continue;
+        }
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1  = ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2  = ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3  = ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]);
+        double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]);
+        double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]);
+        //alternate formulation:
+//        double c1 = 2*ri[jj]*y[kk][0]*(1-y[jj][0]*y[jj][0]);
+//        double c2 = 2*ri[jj]*y[kk][1]*(1-y[jj][1]*y[jj][1]);
+//        double c3 = 2*ri[jj]*y[kk][2]*(1-y[jj][2]*y[jj][2]);
+//        double c10 = 2*ri[kk]*y[jj][0]*(1-y[kk][0]*y[kk][0]);
+//        double c11 = 2*ri[kk]*y[jj][1]*(1-y[kk][1]*y[kk][1]);
+//        double c12 = 2*ri[kk]*y[jj][2]*(1-y[kk][2]*y[kk][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n]/2;
+          double c5 = c51[n]/2;
+          double c6 = c61[n]/2;
+          double ct2 = -ct[n]/2+dfc[kk];
+          double c42 = ct2*y[kk][0];
+          double c52 = ct2*y[kk][1];
+          double c62 = ct2*y[kk][2];
+          //m=0
+          features[count]+=dot1;
+          dfeaturesx[jj*f+count]+=dot1*c4;
+          dfeaturesy[jj*f+count]+=dot1*c5;
+          dfeaturesz[jj*f+count]+=dot1*c6;
+          dfeaturesx[kk*f+count]+=dot1*c42;
+          dfeaturesy[kk*f+count]+=dot1*c52;
+          dfeaturesz[kk*f+count]+=dot1*c62;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          c42*=dot;
+          c52*=dot;
+          c62*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7;
+            dfeaturesy[jj*f+count]+=c8;
+            dfeaturesz[jj*f+count]+=c9;
+            dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
+            dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
+            dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
+            dot1*=dot;
+            features[count++]+=dot1;
+          }
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+

src/USER-RANN/rann_fingerprint_bond.h

@@ -0,0 +1,80 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(bond,Fingerprint_bond)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_BOND_H
+#define LMP_RANN_FINGERPRINT_BOND_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Fingerprint_bond : public Fingerprint {
+     public:
+      Fingerprint_bond(PairRANN *);
+      ~Fingerprint_bond();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,int,int,double *,double *,double *,int *,int,int *);
+      void do3bodyfeatureset_doubleneighborloop(double *,double *,double *,double *,int,int,double *,double *,double *,int *,int,int *);
+      void do3bodyfeatureset_singleneighborloop(double *,double *,double *,double *,int,int,double *,double *,double *,int *,int,int *);
+      void generate_exp_cut_table();
+      void generate_coefficients();
+      int get_length();
+
+      double *expcuttable;
+      double *dfctable;
+      double dr;
+      double *alpha_k;
+      double re;
+      int **coeff;
+      int **coeffx;
+      int **coeffy;
+      int **coeffz;
+      int kmax;
+      int mlength;
+      int **Mf;
+
+    };
+
+  }
+}
+
+
+#endif
+#endif /* FINGERPRINT_BOND_H_ */

src/USER-RANN/rann_fingerprint_bondscreened.cpp

@@ -0,0 +1,762 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#include "rann_fingerprint_bondscreened.h"
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_bondscreened::Fingerprint_bondscreened(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 3;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha_k = new double [1];
+  alpha_k[0] = -1;
+  kmax = 0;
+  mlength = 0;
+  id = -1;
+  style = "bondscreened";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  _pair->doscreen = true;
+  screen = true;
+}
+
+Fingerprint_bondscreened::~Fingerprint_bondscreened() {
+  delete [] alpha_k;
+  delete [] atomtypes;
+  delete [] expcuttable;
+  delete [] dfctable;
+  for (int i=0;i<(mlength*(mlength+1))>>1;i++) {
+    delete [] coeff[i];
+    delete [] coeffx[i];
+    delete [] coeffy[i];
+    delete [] coeffz[i];
+    delete [] Mf[i];
+  }
+  delete [] coeff;
+  delete [] coeffx;
+  delete [] coeffy;
+  delete [] coeffz;
+  delete [] Mf;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_bondscreened::parse_values(std::string constant,std::vector<std::string> line1) {
+  int nwords,l;
+  nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alphak")==0) {
+    delete [] alpha_k;
+    alpha_k = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha_k[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("k")==0) {
+    kmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("m")==0) {
+    mlength = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for bond power");
+  if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_bondscreened::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alphak:\n",style,id);
+  for (i=0;i<kmax;i++) {
+    fprintf(fid,"%f ",alpha_k[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:k:\n",style,id);
+  fprintf(fid,"%d\n",kmax);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:m:\n",style,id);
+  fprintf(fid,"%d\n",mlength);
+}
+
+void Fingerprint_bondscreened::init(int *i,int _id) {
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  re = 0;
+  rc = 0;
+  mlength = 0;
+  kmax = 0;
+  alpha_k = new double [1];
+  alpha_k[0]=-1;
+  empty = false;
+  id = _id;
+}
+
+//number of neurons defined by this fingerprint
+int Fingerprint_bondscreened::get_length() {
+  return mlength*kmax;
+}
+
+void Fingerprint_bondscreened::allocate() {
+  generate_exp_cut_table();
+  generate_coefficients();
+  generate_rinvssqrttable();
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
+void Fingerprint_bondscreened::generate_exp_cut_table() {
+    int m,n;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  expcuttable = new double [(res+buf)*(kmax)];
+  dfctable = new double [res+buf];
+  for (m=0;m<(res+buf);m++) {
+    r1 = cutmax*cutmax*(double)(m)/(double)(res);
+    for (n=0;n<(kmax);n++) {
+      expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+      dfctable[m]=0;
+    }
+    else{
+      dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+}
+
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop.
+void Fingerprint_bondscreened::generate_coefficients() {      //calculates multinomial coefficient for each term
+    int p,mb,mc;
+  int n,p1,i1;
+  mb = mlength;
+  mc=(mb*(mb+1))>>1;
+  coeff = new int *[mc];
+  coeffx = new int *[mc];
+  coeffy = new int *[mc];
+  coeffz = new int *[mc];
+  for (p=0;p<mc;p++) {
+    coeff[p]=new int [mb];
+    coeffx[p]=new int [mb];
+    coeffy[p]=new int [mb];
+    coeffz[p]=new int [mb];
+  }
+  Mf = new int*[mc];
+  int *M = new int[mlength+1];
+  for (p=0;p<mlength+1;p++) {
+    M[p]=0;
+  }
+  for (p1=0;p1<mc;p1++) {
+    Mf[p1] = new int[mlength+1];
+    for (p=0;p<mlength+1;p++) {
+      Mf[p1][p]=0;
+    }
+  }
+  M[0] = 2;
+  Mf[0][0] = 2;
+  n = 1;
+  int m1 = 1;
+  bool go = true;
+  bool broke = false;
+  while (go) {
+    broke = false;
+    for (i1=0;i1<mlength-1;i1++) {
+      if (M[i1+1] == 0) {
+        M[i1+1]=M[i1+1]+1;
+        for (p1=0;p1<mlength+1;p1++) {
+          Mf[n][p1] = M[p1];
+        }
+        n = n+1;
+        broke = true;
+        break;
+      }
+    }
+    if (m1<mlength && !broke) {
+      M[m1]=M[m1]+1;
+      for (p1=m1+1;p1<mlength+1;p1++) {
+        M[p1]=0;
+      }
+      for (p1=0;p1<mlength+1;p1++) {
+        Mf[n][p1]=M[p1];
+      }
+      n=n+1;
+      broke = true;
+      m1 = m1+1;
+    }
+    if (!broke) {
+      go = false;
+    }
+  }
+  for (p=0;p<mb;p++) {
+    for (p1=0;p1<mc;p1++) {
+      if (p==0) {
+        coeffx[p1][p]=0;
+        coeffy[p1][p]=0;
+        coeffz[p1][p]=0;
+      }
+      else{
+        coeffx[p1][p]=coeffx[p1][p-1];
+        if (Mf[p1][p]==0) {
+          coeffx[p1][p]++;
+        }
+        coeffy[p1][p]=coeffy[p1][p-1];
+        if (Mf[p1][p]==1) {
+          coeffy[p1][p]++;
+        }
+        coeffz[p1][p]=coeffz[p1][p-1];
+        if (Mf[p1][p]==2) {
+          coeffz[p1][p]++;
+        }
+      }
+      coeff[p1][p] = pair->factorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]);
+    }
+  }
+}
+
+
+//Called by getproperties. Gets 3-body features and dfeatures
+void Fingerprint_bondscreened::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i;
+  //select the more efficient algorithm for this particular potential and environment.
+  if (jnum*2>(mlength+1)*mlength*20) {
+    do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl);
+  }
+  else{
+    do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl);
+
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers
+void Fingerprint_bondscreened::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  double cutmax = pair->cutmax;
+  int res = pair->res;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  int nelements=pair->nelements;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax+12];
+  int p = kmax;
+  int countmb=((mlength)*(mlength+1))>>1;
+  // calculate interpolation expr, rinvs and dfc, for each neighbor
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+    jtype = tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) {
+      expr[jj][0]=0;
+      continue;
+    }
+    delx=xn[jj];
+    dely=yn[jj];
+    delz=zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    expr[jj][0]=0;
+    continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+      expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+      expr[jj][kk] *= Sik[jj];
+    }
+    double* q = &dfctable[m1-1];
+    double* ri = &rinvsqrttable[m1-1];
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
+
+    expr[jj][p]=delx*rinvs;
+    expr[jj][p+1]=dely*rinvs;
+    expr[jj][p+2]=delz*rinvs;
+    //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
+    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+      expr[jj][p] = 0.000001;
+    }
+    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+      expr[jj][p+1] = 0.000001;
+    }
+    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+      expr[jj][p+2] = 0.000001;
+    }
+    expr[jj][p+3] = -dfc*expr[jj][p]-dSikx[jj];
+    expr[jj][p+4] = rinvs/expr[jj][p];
+    expr[jj][p+5] = rinvs*expr[jj][p];
+    expr[jj][p+6] = -dfc*expr[jj][p+1]-dSiky[jj];
+    expr[jj][p+7] = rinvs/expr[jj][p+1];
+    expr[jj][p+8] = rinvs*expr[jj][p+1];
+    expr[jj][p+9] = -dfc*expr[jj][p+2]-dSikz[jj];
+    expr[jj][p+10] = rinvs/expr[jj][p+2];
+    expr[jj][p+11] = rinvs*expr[jj][p+2];
+  }
+
+  int kb = kmax;
+  int mb = mlength;
+  count = startingneuron;
+  double Bb[mb];
+  double dBbx;
+  double dBby;
+  double dBbz;
+  double dBbx1[mb];
+  double dBby1[mb];
+  double dBbz1[mb];
+  double yprod;
+  for (mcount=0;mcount<countmb;mcount++) {
+    int *M = Mf[mcount];
+    int *_coeffx = coeffx[mcount];
+    int *_coeffy =  coeffy[mcount];
+    int *_coeffz =  coeffz[mcount];
+    int *_coeff =  coeff[mcount];
+    a = mb+1;
+    for (a1=0;a1<mb;a1++) {
+      if (Mf[mcount][a1+1]==0) {
+        a = a1;
+        break;
+      }
+    }
+    for (n=0;n<kb;n++) {
+      for (a1=0;a1<mb;a1++) {
+        Bb[a1]=0;
+        dBbx1[a1] = 0;
+        dBby1[a1] = 0;
+        dBbz1[a1] = 0;
+      }
+      ai = n;
+      double y1 = alpha_k[ai]/re;
+      //loop over jtype to get Bb
+      for (jj=0;jj<jnum;jj++) {
+        if (Bij[jj]==false) {continue;}
+        if (expr[jj][0]==0) {continue;}
+        jtype = tn[jj];
+        if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+          continue;
+        }
+        double yprod = expr[jj][ai];
+        double *y4 = &expr[jj][p];
+        for (a2=0;a2<a;a2++) {
+          yprod *= y4[M[a2+1]];
+        }
+        for (a2=a;a2<mb;a2++) {
+          Bb[a2]=Bb[a2]+yprod;
+          yprod *= y4[M[a2+1]];
+        }
+      }
+      if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
+        double Bg[mb];
+        for (a1=0;a1<mb;a1++) {
+          Bg[a1]=0;
+        }
+        ai = n;
+        double y1 = alpha_k[ai]/re;
+        //loop over ktype to get Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          double yprod = expr[jj][ai];
+          double *y4 = &expr[jj][p];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          for (a2=a;a2<mb;a2++) {
+            Bg[a2]=Bg[a2]+yprod;
+            yprod *= y4[M[a2+1]];
+          }
+        }
+        double B1;
+        //loop over ktype to get dBg*Bb
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bb[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dBbx1[a2] -= dBbx;
+            dBby1[a2] -= dBby;
+            dBbz1[a2] -= dBbz;
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //loop over jtype to get dBb*Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bg[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dBbx1[a2] -= dBbx;
+            dBby1[a2] -= dBby;
+            dBbz1[a2] -= dBbz;
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bb[a2]*Bg[a2]*_coeff[a2];
+        }
+      }
+      else{//Bb=Bg
+        double B1;
+        //loop over jtype to get 2*Bb*dBb
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = 2*Bb[a2]*_coeff[a2]*yprod;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dBbx1[a2] -= dBbx;
+            dBby1[a2] -= dBby;
+            dBbz1[a2] -= dBbz;
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bb[a2]*Bb[a2]*_coeff[a2];
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
+void Fingerprint_bondscreened::do3bodyfeatureset_doubleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,ktype,kk,m,n;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int jtypes = atomtypes[1];
+  int ktypes = atomtypes[2];
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  int nelements = pair->nelements;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax];
+  double y[jnum][3];
+  double ri[jnum];
+  double dfc[jnum];
+  int kb = kmax;
+  int mb = mlength;
+  double Bijk[kb][mb];
+  double c41[kmax];
+  double c51[kmax];
+  double c61[kmax];
+  double ct[kmax];
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+    jtype = tn[jj];
+    if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
+    	expr[jj][0]=0;
+    	continue;
+    }
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    	expr[jj][0]=0;
+    	continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere.
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+    	expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+    	expr[jj][kk] *= Sik[jj];
+    }
+    double* q = &dfctable[m1-1];
+    double* r2 = &rinvsqrttable[m1-1];
+    dfc[jj] = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
+    y[jj][0]=delx*ri[jj];
+    y[jj][1]=dely*ri[jj];
+    y[jj][2]=delz*ri[jj];
+  }
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0)continue;
+    jtype =tn[jj];
+    if (jtypes != nelements && jtypes != jtype) {
+      continue;
+    }
+    for (n = 0;n<kmax;n++) {
+      ct[n] = alpha_k[n]/re;
+      c41[n]=(-ct[n]+dfc[jj])*y[jj][0]+dSikx[jj];
+      c51[n]=(-ct[n]+dfc[jj])*y[jj][1]+dSiky[jj];
+      c61[n]=(-ct[n]+dfc[jj])*y[jj][2]+dSikz[jj];
+    }
+    for (n=0;n<kb;n++) {for (m=0;m<mb;m++) {Bijk[n][m]=0;}}
+    for (kk=0;kk< jnum; kk++) {
+      if (Bij[kk]==false) {continue;}
+      if (expr[kk][0]==0)continue;
+      ktype = tn[kk];
+      if (ktypes != nelements && ktypes != ktype) {
+        continue;
+      }
+      count = startingneuron;
+      double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+      double c1  = ri[jj]*(y[kk][0]-dot*y[jj][0]);
+      double c2  = ri[jj]*(y[kk][1]-dot*y[jj][1]);
+      double c3  = ri[jj]*(y[kk][2]-dot*y[jj][2]);
+      double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]);
+      double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]);
+      double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]);
+      for (n=0;n<kb;n++) {
+        double dot1=expr[jj][n]*expr[kk][n];
+        double c4 = c41[n];
+        double c5 = c51[n];
+        double c6 = c61[n];
+        double ct2 = -ct[n]+dfc[kk];
+        double c42 = ct2*y[kk][0]+dSikx[kk];
+        double c52 = ct2*y[kk][1]+dSiky[kk];
+        double c62 = ct2*y[kk][2]+dSikz[kk];
+        //m=0
+        Bijk[n][0]+=dot1;
+        features[count]+=dot1;
+        dfeaturesx[jj*f+count]+=dot1*c4;
+        dfeaturesy[jj*f+count]+=dot1*c5;
+        dfeaturesz[jj*f+count]+=dot1*c6;
+        dfeaturesx[kk*f+count]+=dot1*c42;
+        dfeaturesy[kk*f+count]+=dot1*c52;
+        dfeaturesz[kk*f+count]+=dot1*c62;
+        c4*=dot;
+        c5*=dot;
+        c6*=dot;
+        c42*=dot;
+        c52*=dot;
+        c62*=dot;
+        count++;
+        for (m=1;m<mb;m++) {
+          double c7 = dot1*(m*c1+c4);
+          double c8 = dot1*(m*c2+c5);
+          double c9 = dot1*(m*c3+c6);
+          dfeaturesx[jj*f+count]+=c7;
+          dfeaturesy[jj*f+count]+=c8;
+          dfeaturesz[jj*f+count]+=c9;
+          dfeaturesx[kk*f+count]+=dot1*(m*c10+c42);
+          dfeaturesy[kk*f+count]+=dot1*(m*c11+c52);
+          dfeaturesz[kk*f+count]+=dot1*(m*c12+c62);
+          dot1*=dot;
+          features[count++]+=dot1;
+          Bijk[n][m] += dot1;
+        }
+      }
+    }
+    for (kk=0;kk<jnum;kk++) {
+      if (Bij[kk]==false) {continue;}
+      if (expr[kk][0]==0)continue;
+      ktype = tn[kk];
+      if (ktypes != nelements && ktypes != ktype) {
+        continue;
+      }
+      count = startingneuron;
+      for (n=0;n<kb;n++) {
+        for (m=0;m<mb;m++) {
+          dfeaturesx[kk*f+count]+=2*Bijk[n][m]*dSijkx[jj*jnum+kk];
+          dfeaturesy[kk*f+count]+=2*Bijk[n][m]*dSijky[jj*jnum+kk];
+          dfeaturesz[kk*f+count]+=2*Bijk[n][m]*dSijkz[jj*jnum+kk];
+          count++;
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+

src/USER-RANN/rann_fingerprint_bondscreened.h

@@ -0,0 +1,80 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(bondscreened,Fingerprint_bondscreened)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_BONDSCREENED_H
+#define LMP_RANN_FINGERPRINT_BONDSCREENED_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Fingerprint_bondscreened : public Fingerprint {
+     public:
+      Fingerprint_bondscreened(PairRANN *);
+      ~Fingerprint_bondscreened();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      void do3bodyfeatureset_doubleneighborloop(double *,double *,double *,double *,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      void do3bodyfeatureset_singleneighborloop(double *,double *,double *,double *,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      void generate_exp_cut_table();
+      void generate_coefficients();
+      int get_length();
+
+      double *expcuttable;
+      double *dfctable;
+      double dr;
+      double *alpha_k;
+      double re;
+      int **coeff;
+      int **coeffx;
+      int **coeffy;
+      int **coeffz;
+      int kmax;
+      int mlength;
+      int **Mf;
+
+    };
+  }
+
+}
+
+
+#endif
+#endif /* FINGERPRINT_BOND_H_ */

src/USER-RANN/rann_fingerprint_bondscreenedspin.cpp

@@ -0,0 +1,814 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#include "rann_fingerprint_bondscreenedspin.h"
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_bondscreenedspin::Fingerprint_bondscreenedspin(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 3;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha_k = new double [1];
+  alpha_k[0] = -1;
+  kmax = 0;
+  mlength = 0;
+  id = -1;
+  style = "bondscreenedspin";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  _pair->doscreen = true;
+  screen = true;
+  _pair->dospin = true;
+  spin = true;
+}
+
+Fingerprint_bondscreenedspin::~Fingerprint_bondscreenedspin() {
+  delete [] alpha_k;
+  delete [] atomtypes;
+  delete [] expcuttable;
+  delete [] dfctable;
+  for (int i=0;i<(mlength*(mlength+1))>>1;i++) {
+    delete [] coeff[i];
+    delete [] coeffx[i];
+    delete [] coeffy[i];
+    delete [] coeffz[i];
+    delete [] Mf[i];
+  }
+  delete [] coeff;
+  delete [] coeffx;
+  delete [] coeffy;
+  delete [] coeffz;
+  delete [] Mf;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_bondscreenedspin::parse_values(std::string constant,std::vector<std::string> line1) {
+  int nwords,l;
+  nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alphak")==0) {
+    delete [] alpha_k;
+    alpha_k = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha_k[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("k")==0) {
+    kmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("m")==0) {
+    mlength = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for bond power");
+  if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_bondscreenedspin::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alphak:\n",style,id);
+  for (i=0;i<kmax;i++) {
+    fprintf(fid,"%f ",alpha_k[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:k:\n",style,id);
+  fprintf(fid,"%d\n",kmax);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:m:\n",style,id);
+  fprintf(fid,"%d\n",mlength);
+}
+
+void Fingerprint_bondscreenedspin::init(int *i,int _id) {
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  re = 0;
+  rc = 0;
+  mlength = 0;
+  kmax = 0;
+  alpha_k = new double [1];
+  alpha_k[0]=-1;
+  empty = false;
+  id = _id;
+}
+
+//number of neurons defined by this fingerprint
+int Fingerprint_bondscreenedspin::get_length() {
+  return mlength*kmax;
+}
+
+void Fingerprint_bondscreenedspin::allocate() {
+  generate_exp_cut_table();
+  generate_coefficients();
+  generate_rinvssqrttable();
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
+void Fingerprint_bondscreenedspin::generate_exp_cut_table() {
+    int m,n;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  expcuttable = new double [(res+buf)*(kmax)];
+  dfctable = new double [res+buf];
+  for (m=0;m<(res+buf);m++) {
+    r1 = cutmax*cutmax*(double)(m)/(double)(res);
+    for (n=0;n<(kmax);n++) {
+      expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+      dfctable[m]=0;
+    }
+    else{
+      dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+}
+
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop.
+void Fingerprint_bondscreenedspin::generate_coefficients() {      //calculates multinomial coefficient for each term
+    int p,mb,mc;
+  int n,p1,i1;
+  mb = mlength;
+  mc=(mb*(mb+1))>>1;
+  coeff = new int *[mc];
+  coeffx = new int *[mc];
+  coeffy = new int *[mc];
+  coeffz = new int *[mc];
+  for (p=0;p<mc;p++) {
+    coeff[p]=new int [mb];
+    coeffx[p]=new int [mb];
+    coeffy[p]=new int [mb];
+    coeffz[p]=new int [mb];
+  }
+  Mf = new int*[mc];
+  int *M = new int[mlength+1];
+  for (p=0;p<mlength+1;p++) {
+    M[p]=0;
+  }
+  for (p1=0;p1<mc;p1++) {
+    Mf[p1] = new int[mlength+1];
+    for (p=0;p<mlength+1;p++) {
+      Mf[p1][p]=0;
+    }
+  }
+  M[0] = 2;
+  Mf[0][0] = 2;
+  n = 1;
+  int m1 = 1;
+  bool go = true;
+  bool broke = false;
+  while (go) {
+    broke = false;
+    for (i1=0;i1<mlength-1;i1++) {
+      if (M[i1+1] == 0) {
+        M[i1+1]=M[i1+1]+1;
+        for (p1=0;p1<mlength+1;p1++) {
+          Mf[n][p1] = M[p1];
+        }
+        n = n+1;
+        broke = true;
+        break;
+      }
+    }
+    if (m1<mlength && !broke) {
+      M[m1]=M[m1]+1;
+      for (p1=m1+1;p1<mlength+1;p1++) {
+        M[p1]=0;
+      }
+      for (p1=0;p1<mlength+1;p1++) {
+        Mf[n][p1]=M[p1];
+      }
+      n=n+1;
+      broke = true;
+      m1 = m1+1;
+    }
+    if (!broke) {
+      go = false;
+    }
+  }
+  for (p=0;p<mb;p++) {
+    for (p1=0;p1<mc;p1++) {
+      if (p==0) {
+        coeffx[p1][p]=0;
+        coeffy[p1][p]=0;
+        coeffz[p1][p]=0;
+      }
+      else{
+        coeffx[p1][p]=coeffx[p1][p-1];
+        if (Mf[p1][p]==0) {
+          coeffx[p1][p]++;
+        }
+        coeffy[p1][p]=coeffy[p1][p-1];
+        if (Mf[p1][p]==1) {
+          coeffy[p1][p]++;
+        }
+        coeffz[p1][p]=coeffz[p1][p-1];
+        if (Mf[p1][p]==2) {
+          coeffz[p1][p]++;
+        }
+      }
+      coeff[p1][p] = pair->factorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]);
+    }
+  }
+}
+
+
+//Called by getproperties. Gets 3-body features and dfeatures
+void Fingerprint_bondscreenedspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, double *dspinx, double *dspiny, double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i;
+  //select the more efficient algorithm for this particular potential and environment.
+  if (jnum*2>(mlength+1)*mlength*20) {
+    do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl);
+  }
+  else{
+    do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,sid,xn,yn,zn,tn,jnum,jl);
+
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers
+void Fingerprint_bondscreenedspin::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, double *dspinx, double *dspiny, double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  double cutmax = pair->cutmax;
+  int res = pair->res;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  int nelements=pair->nelements;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax+12];
+  int p = kmax;
+  int countmb=((mlength)*(mlength+1))>>1;
+  double *si = sim->s[i];
+  // calculate interpolation expr, rinvs and dfc, for each neighbor
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+    jtype = tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) {
+      expr[jj][0]=0;
+      continue;
+    }
+    delx=xn[jj];
+    dely=yn[jj];
+    delz=zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    expr[jj][0]=0;
+    continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+      expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+      expr[jj][kk] *= Sik[jj];
+    }
+    double* q = &dfctable[m1-1];
+    double* ri = &rinvsqrttable[m1-1];
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
+
+    expr[jj][p]=delx*rinvs;
+    expr[jj][p+1]=dely*rinvs;
+    expr[jj][p+2]=delz*rinvs;
+    //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
+    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+      expr[jj][p] = 0.000001;
+    }
+    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+      expr[jj][p+1] = 0.000001;
+    }
+    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+      expr[jj][p+2] = 0.000001;
+    }
+    expr[jj][p+3] = -dfc*expr[jj][p]-dSikx[jj];
+    expr[jj][p+4] = rinvs/expr[jj][p];
+    expr[jj][p+5] = rinvs*expr[jj][p];
+    expr[jj][p+6] = -dfc*expr[jj][p+1]-dSiky[jj];
+    expr[jj][p+7] = rinvs/expr[jj][p+1];
+    expr[jj][p+8] = rinvs*expr[jj][p+1];
+    expr[jj][p+9] = -dfc*expr[jj][p+2]-dSikz[jj];
+    expr[jj][p+10] = rinvs/expr[jj][p+2];
+    expr[jj][p+11] = rinvs*expr[jj][p+2];
+  }
+
+  int kb = kmax;
+  int mb = mlength;
+  count = startingneuron;
+  double Bb[mb];
+  double Bbs[mb];
+  double dBbx;
+  double dBby;
+  double dBbz;
+  double yprod;
+  for (mcount=0;mcount<countmb;mcount++) {
+    int *M = Mf[mcount];
+    int *_coeffx = coeffx[mcount];
+    int *_coeffy = coeffy[mcount];
+    int *_coeffz = coeffz[mcount];
+    int *_coeff = coeff[mcount];
+    a = mb+1;
+    for (a1=0;a1<mb;a1++) {
+      if (Mf[mcount][a1+1]==0) {
+        a = a1;
+        break;
+      }
+    }
+    for (n=0;n<kb;n++) {
+      for (a1=0;a1<mb;a1++) {
+        Bb[a1]=0;
+        Bbs[a1]=0;
+      }
+      ai = n;
+      double y1 = alpha_k[ai]/re;
+      //loop over jtype to get Bb
+      for (jj=0;jj<jnum;jj++) {
+        if (Bij[jj]==false) {continue;}
+        if (expr[jj][0]==0) {continue;}
+        jtype = tn[jj];
+        if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+          continue;
+        }
+        j = jl[jj];
+        double *sj = sim->s[j];
+        double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+        double yprod = expr[jj][ai];
+        double *y4 = &expr[jj][p];
+        for (a2=0;a2<a;a2++) {
+          yprod *= y4[M[a2+1]];
+        }
+        for (a2=a;a2<mb;a2++) {
+          Bb[a2]=Bb[a2]+yprod;
+          Bbs[a2]=Bbs[a2]+yprod*sp;
+          yprod *= y4[M[a2+1]];
+        }
+      }
+      if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
+        double Bg[mb];
+        double Bgs[mb];
+        for (a1=0;a1<mb;a1++) {
+          Bg[a1]=0;
+          Bgs[a1]=0;
+        }
+        ai = n;
+        double y1 = alpha_k[ai]/re;
+        //loop over ktype to get Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double yprod = expr[jj][ai];
+          double *y4 = &expr[jj][p];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          for (a2=a;a2<mb;a2++) {
+            Bg[a2]=Bg[a2]+yprod;
+            Bgs[a2]=Bgs[a2]+yprod*sp;
+            yprod *= y4[M[a2+1]];
+          }
+        }
+        double B1;
+        //loop over ktype to get dBg*Bb
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bbs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += yprod*Bbs[a2]*_coeff[a2]*si[0];
+            dspiny[ai] += yprod*Bbs[a2]*_coeff[a2]*si[1];
+            dspinz[ai] += yprod*Bbs[a2]*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //loop over jtype to get dBb*Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bgs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += Bgs[a2]*yprod*_coeff[a2]*si[0];
+            dspiny[ai] += Bgs[a2]*yprod*_coeff[a2]*si[1];
+            dspinz[ai] += Bgs[a2]*yprod*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bbs[a2]*Bgs[a2]*_coeff[a2];
+        }
+      }
+      else{//Bb=Bg
+        double B1;
+        //loop over jtype to get 2*Bb*dBb
+        for (jj=0;jj<jnum;jj++) {
+          if (Bij[jj]==false) {continue;}
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = 2*Bbs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[0];
+            dspiny[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[1];
+            dspinz[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+            for (kk=0;kk<jnum;kk++) {
+              if (Bij[kk]==false) {continue;}
+              dfeaturesx[n*mb+a2+count+kk*f]+=B1*dSijkx[jj*jnum+kk];
+              dfeaturesy[n*mb+a2+count+kk*f]+=B1*dSijky[jj*jnum+kk];
+              dfeaturesz[n*mb+a2+count+kk*f]+=B1*dSijkz[jj*jnum+kk];
+            }
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bbs[a2]*Bbs[a2]*_coeff[a2];
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
+void Fingerprint_bondscreenedspin::do3bodyfeatureset_doubleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz, double *dspinx, double *dspiny, double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij, int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,ktype,kk,m,n;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int jtypes = atomtypes[1];
+  int ktypes = atomtypes[2];
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  int nelements = pair->nelements;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax];
+  double y[jnum][3];
+  double ri[jnum];
+  double dfc[jnum];
+  int kb = kmax;
+  int mb = mlength;
+  double Bijk[kb][mb];
+  double c41[kmax];
+  double c51[kmax];
+  double c61[kmax];
+  double ct[kmax];
+  double *si = sim->s[i];
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+      jtype = tn[jj];
+      if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
+        expr[jj][0]=0;
+        continue;
+      }
+      delx = xn[jj];
+      dely = yn[jj];
+      delz = zn[jj];
+      rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq>rc*rc) {
+        expr[jj][0]=0;
+        continue;
+      }
+      double r1 = (rsq*((double)res)*cutinv2);
+      int m1 = (int)r1;
+      if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere.
+      r1 = r1-trunc(r1);
+      double *p0 = &expcuttable[(m1-1)*kmax];
+      double *p1 = &expcuttable[m1*kmax];
+      double *p2 = &expcuttable[(m1+1)*kmax];
+      double *p3 = &expcuttable[(m1+2)*kmax];
+      for (kk=0;kk<kmax;kk++) {
+        expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+        expr[jj][kk] *= Sik[jj];
+      }
+      double* q = &dfctable[m1-1];
+      double* r2 = &rinvsqrttable[m1-1];
+      dfc[jj] = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+      ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
+      y[jj][0]=delx*ri[jj];
+      y[jj][1]=dely*ri[jj];
+      y[jj][2]=delz*ri[jj];
+  }
+  for (jj = 0; jj < jnum; jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0)continue;
+    jtype =tn[jj];
+    if (jtypes != nelements && jtypes != jtype) {
+      continue;
+    }
+    j = jl[jj];
+    double *sj = sim->s[j];
+    double spj = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+    for (n = 0;n<kmax;n++) {
+      ct[n] = alpha_k[n]/re;
+      c41[n]=(-ct[n]+dfc[jj])*y[jj][0]+dSikx[jj];
+      c51[n]=(-ct[n]+dfc[jj])*y[jj][1]+dSiky[jj];
+      c61[n]=(-ct[n]+dfc[jj])*y[jj][2]+dSikz[jj];
+    }
+    for (n=0;n<kb;n++) {for (m=0;m<mb;m++) {Bijk[n][m]=0;}}
+    for (kk=0;kk< jnum; kk++) {
+      if (Bij[kk]==false) {continue;}
+      if (expr[kk][0]==0)continue;
+      ktype = tn[kk];
+      if (ktypes != nelements && ktypes != ktype) {
+        continue;
+      }
+      j = jl[kk];
+      double *sk = sim->s[j];
+      double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2];
+      count = startingneuron;
+      double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+      double c1  = ri[jj]*(y[kk][0]-dot*y[jj][0]);
+      double c2  = ri[jj]*(y[kk][1]-dot*y[jj][1]);
+      double c3  = ri[jj]*(y[kk][2]-dot*y[jj][2]);
+      double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]);
+      double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]);
+      double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]);
+      for (n=0;n<kb;n++) {
+        double dot1=expr[jj][n]*expr[kk][n];
+        double c4 = c41[n];
+        double c5 = c51[n];
+        double c6 = c61[n];
+        double ct2 = -ct[n]+dfc[kk];
+        double c42 = ct2*y[kk][0]+dSikx[kk];
+        double c52 = ct2*y[kk][1]+dSiky[kk];
+        double c62 = ct2*y[kk][2]+dSikz[kk];
+        //m=0
+        Bijk[n][0]+=dot1*spj*spk;
+        features[count]+=dot1*spj*spk;
+        dspinx[jj*f+count] += dot1*si[0]*spk;
+        dspiny[jj*f+count] += dot1*si[1]*spk;
+        dspinz[jj*f+count] += dot1*si[2]*spk;
+        dspinx[kk*f+count] += dot1*si[0]*spj;
+        dspiny[kk*f+count] += dot1*si[1]*spj;
+        dspinz[kk*f+count] += dot1*si[2]*spj;
+        dspinx[jnum*f+count] += dot1*(sj[0]*spk+sk[0]*spj);
+        dspiny[jnum*f+count] += dot1*(sj[1]*spk+sk[1]*spj);
+        dspinz[jnum*f+count] += dot1*(sj[2]*spk+sk[2]*spj);
+        dfeaturesx[jj*f+count]+=dot1*c4*spj*spk;
+        dfeaturesy[jj*f+count]+=dot1*c5*spj*spk;
+        dfeaturesz[jj*f+count]+=dot1*c6*spj*spk;
+        dfeaturesx[kk*f+count]+=dot1*c42*spj*spk;
+        dfeaturesy[kk*f+count]+=dot1*c52*spj*spk;
+        dfeaturesz[kk*f+count]+=dot1*c62*spj*spk;
+        c4*=dot;
+        c5*=dot;
+        c6*=dot;
+        c42*=dot;
+        c52*=dot;
+        c62*=dot;
+        count++;
+        for (m=1;m<mb;m++) {
+          double c7 = dot1*(m*c1+c4);
+          double c8 = dot1*(m*c2+c5);
+          double c9 = dot1*(m*c3+c6);
+          dfeaturesx[jj*f+count]+=c7*spj*spk;
+          dfeaturesy[jj*f+count]+=c8*spj*spk;
+          dfeaturesz[jj*f+count]+=c9*spj*spk;
+          dfeaturesx[kk*f+count]+=dot1*(m*c10+c42)*spj*spk;
+          dfeaturesy[kk*f+count]+=dot1*(m*c11+c52)*spj*spk;
+          dfeaturesz[kk*f+count]+=dot1*(m*c12+c62)*spj*spk;
+          dot1*=dot;
+          dspinx[jj*f+count] += dot1*si[0]*spk;
+          dspiny[jj*f+count] += dot1*si[1]*spk;
+          dspinz[jj*f+count] += dot1*si[2]*spk;
+          dspinx[kk*f+count] += dot1*si[0]*spj;
+          dspiny[kk*f+count] += dot1*si[1]*spj;
+          dspinz[kk*f+count] += dot1*si[2]*spj;
+          dspinx[jnum*f+count] += dot1*(sj[0]*spk+sk[0]*spj);
+          dspiny[jnum*f+count] += dot1*(sj[1]*spk+sk[1]*spj);
+          dspinz[jnum*f+count] += dot1*(sj[2]*spk+sk[2]*spj);
+          Bijk[n][m] += dot1*spj*spk;
+          features[count++]+=dot1*spj*spk;
+        }
+      }
+    }
+    for (kk=0;kk<jnum;kk++) {
+      if (Bij[kk]==false) {continue;}
+      if (expr[kk][0]==0)continue;
+      ktype = tn[kk];
+      if (ktypes != nelements && ktypes != ktype) {
+        continue;
+      }
+      count = startingneuron;
+      for (n=0;n<kb;n++) {
+        for (m=0;m<mb;m++) {
+          dfeaturesx[kk*f+count]+=2*Bijk[n][m]*dSijkx[jj*jnum+kk];
+          dfeaturesy[kk*f+count]+=2*Bijk[n][m]*dSijky[jj*jnum+kk];
+          dfeaturesz[kk*f+count]+=2*Bijk[n][m]*dSijkz[jj*jnum+kk];
+          count++;
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (Bij[jj]==false) {continue;}
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+

src/USER-RANN/rann_fingerprint_bondscreenedspin.h

@@ -0,0 +1,80 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(bondscreenedspin,Fingerprint_bondscreenedspin)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_BONDSCREENEDSPIN_H
+#define LMP_RANN_FINGERPRINT_BONDSCREENEDSPIN_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Fingerprint_bondscreenedspin : public Fingerprint {
+     public:
+      Fingerprint_bondscreenedspin(PairRANN *);
+      ~Fingerprint_bondscreenedspin();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      void do3bodyfeatureset_doubleneighborloop(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int *);
+      void do3bodyfeatureset_singleneighborloop(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      void generate_exp_cut_table();
+      void generate_coefficients();
+      int get_length();
+
+      double *expcuttable;
+      double *dfctable;
+      double dr;
+      double *alpha_k;
+      double re;
+      int **coeff;
+      int **coeffx;
+      int **coeffy;
+      int **coeffz;
+      int kmax;
+      int mlength;
+      int **Mf;
+
+    };
+
+  }
+}
+
+
+#endif
+#endif /* FINGERPRINT_BOND_H_ */

src/USER-RANN/rann_fingerprint_bondspin.cpp

@@ -0,0 +1,887 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#include "rann_fingerprint_bondspin.h"
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_bondspin::Fingerprint_bondspin(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 3;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha_k = new double [1];
+  alpha_k[0] = -1;
+  kmax = 0;
+  mlength = 0;
+  id = -1;
+  style = "bondspin";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  _pair->allscreen = false;
+  _pair->dospin = true;
+  spin = true;
+}
+
+Fingerprint_bondspin::~Fingerprint_bondspin() {
+  delete [] alpha_k;
+  delete [] atomtypes;
+  delete [] expcuttable;
+  delete [] dfctable;
+  for (int i=0;i<(mlength*(mlength+1))>>1;i++) {
+    delete [] coeff[i];
+    delete [] coeffx[i];
+    delete [] coeffy[i];
+    delete [] coeffz[i];
+    delete [] Mf[i];
+  }
+  delete [] coeff;
+  delete [] coeffx;
+  delete [] coeffy;
+  delete [] coeffz;
+  delete [] Mf;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_bondspin::parse_values(std::string constant,std::vector<std::string> line1) {
+  int nwords,l;
+  nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alphak")==0) {
+    delete [] alpha_k;
+    alpha_k = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha_k[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("k")==0) {
+    kmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("m")==0) {
+    mlength = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for bond power");
+  if (re!=0.0 && rc!=0.0 && alpha_k[0]!=-1 && dr!=0.0 && mlength!=0 && kmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_bondspin::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alphak:\n",style,id);
+  for (i=0;i<kmax;i++) {
+    fprintf(fid,"%f ",alpha_k[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:k:\n",style,id);
+  fprintf(fid,"%d\n",kmax);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:m:\n",style,id);
+  fprintf(fid,"%d\n",mlength);
+}
+
+void Fingerprint_bondspin::init(int *i,int _id) {
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  re = 0;
+  rc = 0;
+  mlength = 0;
+  kmax = 0;
+  alpha_k = new double [1];
+  alpha_k[0]=-1;
+  empty = false;
+  id = _id;
+}
+
+//number of neurons defined by this fingerprint
+int Fingerprint_bondspin::get_length() {
+  return mlength*kmax;
+}
+
+void Fingerprint_bondspin::allocate() {
+  generate_exp_cut_table();
+  generate_coefficients();
+  generate_rinvssqrttable();
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop and do3bodyfeatureset_doubleneighborloop.
+void Fingerprint_bondspin::generate_exp_cut_table() {
+  int m,n;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  expcuttable = new double [(res+buf)*(kmax)];
+  dfctable = new double [res+buf];
+  for (m=0;m<(res+buf);m++) {
+    r1 = cutmax*cutmax*(double)(m)/(double)(res);
+    for (n=0;n<(kmax);n++) {
+      expcuttable[n+m*(kmax)] = exp(-alpha_k[n]/re*sqrt(r1))*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+      dfctable[m]=0;
+    }
+    else{
+      dfctable[m]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+}
+
+//Generate table of complex functions for quick reference during compute. Used by do3bodyfeatureset_singleneighborloop.
+void Fingerprint_bondspin::generate_coefficients() {      //calculates multinomial coefficient for each term
+    int p,mb,mc;
+  int n,p1,i1;
+  mb = mlength;
+  mc=(mb*(mb+1))>>1;
+  coeff = new int *[mc];
+  coeffx = new int *[mc];
+  coeffy = new int *[mc];
+  coeffz = new int *[mc];
+  for (p=0;p<mc;p++) {
+    coeff[p]=new int [mb];
+    coeffx[p]=new int [mb];
+    coeffy[p]=new int [mb];
+    coeffz[p]=new int [mb];
+  }
+  Mf = new int*[mc];
+  int *M = new int[mlength+1];
+  for (p=0;p<mlength+1;p++) {
+    M[p]=0;
+  }
+  for (p1=0;p1<mc;p1++) {
+    Mf[p1] = new int[mlength+1];
+    for (p=0;p<mlength+1;p++) {
+      Mf[p1][p]=0;
+    }
+  }
+  M[0] = 2;
+  Mf[0][0] = 2;
+  n = 1;
+  int m1 = 1;
+  bool go = true;
+  bool broke = false;
+  while (go) {
+    broke = false;
+    for (i1=0;i1<mlength-1;i1++) {
+      if (M[i1+1] == 0) {
+        M[i1+1]=M[i1+1]+1;
+        for (p1=0;p1<mlength+1;p1++) {
+          Mf[n][p1] = M[p1];
+        }
+        n = n+1;
+        broke = true;
+        break;
+      }
+    }
+    if (m1<mlength && !broke) {
+      M[m1]=M[m1]+1;
+      for (p1=m1+1;p1<mlength+1;p1++) {
+        M[p1]=0;
+      }
+      for (p1=0;p1<mlength+1;p1++) {
+        Mf[n][p1]=M[p1];
+      }
+      n=n+1;
+      broke = true;
+      m1 = m1+1;
+    }
+    if (!broke) {
+      go = false;
+    }
+  }
+  for (p=0;p<mb;p++) {
+    for (p1=0;p1<mc;p1++) {
+      if (p==0) {
+        coeffx[p1][p]=0;
+        coeffy[p1][p]=0;
+        coeffz[p1][p]=0;
+      }
+      else{
+        coeffx[p1][p]=coeffx[p1][p-1];
+        if (Mf[p1][p]==0) {
+          coeffx[p1][p]++;
+        }
+        coeffy[p1][p]=coeffy[p1][p-1];
+        if (Mf[p1][p]==1) {
+          coeffy[p1][p]++;
+        }
+        coeffz[p1][p]=coeffz[p1][p-1];
+        if (Mf[p1][p]==2) {
+          coeffz[p1][p]++;
+        }
+      }
+      coeff[p1][p] = pair->factorial(p)/pair->factorial(coeffx[p1][p])/pair->factorial(coeffy[p1][p])/pair->factorial(coeffz[p1][p]);
+    }
+  }
+}
+
+
+//Called by getproperties. Gets 3-body features and dfeatures
+void Fingerprint_bondspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i;
+  int *ilist,*numneigh;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  ilist = sim->ilist;
+  numneigh = sim->numneigh;
+  i = ilist[ii];
+  //select the more efficient algorithm for this particular potential and environment.
+  if (jnum*2>(mlength+1)*mlength*20) {
+    do3bodyfeatureset_singleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,ii,sid,xn,yn,zn,tn,jnum,jl);
+  }
+  else{
+    do3bodyfeatureset_doubleneighborloop(features,dfeaturesx,dfeaturesy,dfeaturesz,dspinx,dspiny,dspinz,ii,sid,xn,yn,zn,tn,jnum,jl);
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for high neighbor numbers and small series of bond angle powers
+void Fingerprint_bondspin::do3bodyfeatureset_singleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,kk,m,n,mcount,a,a1,a2,ai;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int *type = sim->type;
+  double cutmax = pair->cutmax;
+  int res = pair->res;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  int nelements=pair->nelements;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax+12];
+  int p = kmax;
+  int countmb=((mlength)*(mlength+1))>>1;
+  double *si = sim->s[i];
+  // calculate interpolation expr, rinvs and dfc, for each neighbor
+  for (jj = 0; jj < jnum; jj++) {
+    jtype = tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype && atomtypes[2] != nelements && atomtypes[2] != jtype) {
+      expr[jj][0]=0;
+      continue;
+    }
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    expr[jj][0]=0;
+    continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+      expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+    }
+    double* q = &dfctable[m1-1];
+    double* ri = &rinvsqrttable[m1-1];
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double rinvs = ri[1] + 0.5 * r1*(ri[2] - ri[0] + r1*(2.0*ri[0] - 5.0*ri[1] + 4.0*ri[2] - ri[3] + r1*(3.0*(ri[1] - ri[2]) + ri[3] - ri[0])));
+
+    expr[jj][p]=delx*rinvs;
+    expr[jj][p+1]=dely*rinvs;
+    expr[jj][p+2]=delz*rinvs;
+    //Hack to avoid nan when x y or z component of radial vector is exactly 0. Shouldn't affect accuracy.
+    if (expr[jj][p]*expr[jj][p]<0.000000000001) {
+      expr[jj][p] = 0.000001;
+    }
+    if (expr[jj][p+1]*expr[jj][p+1]<0.000000000001) {
+      expr[jj][p+1] = 0.000001;
+    }
+    if (expr[jj][p+2]*expr[jj][p+2]<0.000000000001) {
+      expr[jj][p+2] = 0.000001;
+    }
+    expr[jj][p+3] = -dfc*expr[jj][p];
+    expr[jj][p+4] = rinvs/expr[jj][p];
+    expr[jj][p+5] = rinvs*expr[jj][p];
+    expr[jj][p+6] = -dfc*expr[jj][p+1];
+    expr[jj][p+7] = rinvs/expr[jj][p+1];
+    expr[jj][p+8] = rinvs*expr[jj][p+1];
+    expr[jj][p+9] = -dfc*expr[jj][p+2];
+    expr[jj][p+10] = rinvs/expr[jj][p+2];
+    expr[jj][p+11] = rinvs*expr[jj][p+2];
+  }
+
+  int kb = kmax;
+  int mb = mlength;
+  count = startingneuron;
+  double Bb[mb];
+  double Bbs[mb];
+  double dBbx;
+  double dBby;
+  double dBbz;
+  double yprod;
+  for (mcount=0;mcount<countmb;mcount++) {
+    int *M = Mf[mcount];
+    int *_coeffx = coeffx[mcount];
+    int *_coeffy = coeffy[mcount];
+    int *_coeffz = coeffz[mcount];
+    int *_coeff = coeff[mcount];
+    a = mb+1;
+    for (a1=0;a1<mb;a1++) {
+      if (Mf[mcount][a1+1]==0) {
+        a = a1;
+        break;
+      }
+    }
+    for (n=0;n<kb;n++) {
+      for (a1=0;a1<mb;a1++) {
+        Bb[a1]=0;
+        Bbs[a1]=0;
+      }
+      ai = n;
+      double y1 = alpha_k[ai]/re;
+      //loop over jtype to get Bb
+      for (jj=0;jj<jnum;jj++) {
+        if (expr[jj][0]==0) {continue;}
+        jtype = tn[jj];
+        if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+          continue;
+        }
+        j = jl[jj];
+        double *sj = sim->s[j];
+        double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+        double yprod = expr[jj][ai];
+        double *y4 = &expr[jj][p];
+        for (a2=0;a2<a;a2++) {
+          yprod *= y4[M[a2+1]];
+        }
+        for (a2=a;a2<mb;a2++) {
+          Bb[a2]=Bb[a2]+yprod;
+          Bbs[a2]=Bbs[a2]+yprod*sp;
+          yprod *= y4[M[a2+1]];
+        }
+      }
+      if (atomtypes[1]!=atomtypes[2]) {//Bb!=Bg
+        double Bg[mb];
+        double Bgs[mb];
+        for (a1=0;a1<mb;a1++) {
+          Bg[a1]=0;
+          Bgs[a1]=0;
+        }
+        ai = n;
+        double y1 = alpha_k[ai]/re;
+        //loop over ktype to get Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double yprod = expr[jj][ai];
+          double *y4 = &expr[jj][p];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          for (a2=a;a2<mb;a2++) {
+            Bg[a2]=Bg[a2]+yprod;
+            Bgs[a2]=Bgs[a2]+yprod*sp;
+            yprod *= y4[M[a2+1]];
+          }
+        }
+        double B1;
+        //loop over ktype to get dBg*Bb
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[2] != nelements && atomtypes[2] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bbs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += yprod*Bbs[a2]*_coeff[a2]*si[0];
+            dspiny[ai] += yprod*Bbs[a2]*_coeff[a2]*si[1];
+            dspinz[ai] += yprod*Bbs[a2]*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += yprod*Bbs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //loop over jtype to get dBb*Bg
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = Bgs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += Bgs[a2]*yprod*_coeff[a2]*si[0];
+            dspiny[ai] += Bgs[a2]*yprod*_coeff[a2]*si[1];
+            dspinz[ai] += Bgs[a2]*yprod*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += yprod*Bgs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bbs[a2]*Bgs[a2]*_coeff[a2];
+        }
+      }
+      else{//Bb=Bg
+        double B1;
+        //loop over jtype to get 2*Bb*dBb
+        for (jj=0;jj<jnum;jj++) {
+          if (expr[jj][0]==0) {continue;}
+          jtype = tn[jj];
+          if (atomtypes[1] != nelements && atomtypes[1] != jtype) {
+            continue;
+          }
+          j = jl[jj];
+          double *sj = sim->s[j];
+          double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+          double *y3 = &expr[jj][p+3];
+          double *y4 = &expr[jj][p];
+          ai = n;
+          yprod = expr[jj][ai];
+          for (a2=0;a2<a;a2++) {
+            yprod *= y4[M[a2+1]];
+          }
+          ai = n*(mb)+a+count+jj*f;
+          for (a2=a;a2<mb;a2++) {
+            B1 = 2*Bbs[a2]*_coeff[a2]*yprod*sp;
+            dBbx = -B1*(y1*y4[0]+y3[0]-_coeffx[a2]*y3[1]+a2*y3[2]);
+            dBby = -B1*(y1*y4[1]+y3[3]-_coeffy[a2]*y3[4]+a2*y3[5]);
+            dBbz = -B1*(y1*y4[2]+y3[6]-_coeffz[a2]*y3[7]+a2*y3[8]);
+            dfeaturesx[ai] += dBbx;
+            dfeaturesy[ai] += dBby;
+            dfeaturesz[ai] += dBbz;
+            dspinx[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[0];
+            dspiny[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[1];
+            dspinz[ai] += 2*Bbs[a2]*yprod*_coeff[a2]*si[2];
+            dspinx[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[0];
+            dspiny[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[1];
+            dspinz[ai-jj*f+jnum*f] += 2*yprod*Bbs[a2]*_coeff[a2]*sj[2];
+            yprod *= y4[M[a2+1]];
+            ai++;
+          }
+        }
+        //central atom derivative
+        for (a2=a;a2<mb;a2++) {
+          ai = n*(mb)+a2+count+jnum*f;
+          features[ai-jnum*f] += Bbs[a2]*Bbs[a2]*_coeff[a2];
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+
+//Called by do3bodyfeatureset. Algorithm for low neighbor numbers and large series of bond angle powers
+void Fingerprint_bondspin::do3bodyfeatureset_doubleneighborloop(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl) {
+  int i,j,jj,itype,jtype,ktype,kk,m,n;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  int jtypes = atomtypes[1];
+  int ktypes = atomtypes[2];
+  int count=0;
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  double **x = sim->x;
+  int *type = sim->type;
+  int nelements = pair->nelements;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  double cutinv2 = 1/cutmax/cutmax;
+  ilist = sim->ilist;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double expr[jnum][kmax];
+  double y[jnum][3];
+  double ri[jnum];
+  double dfc[jnum];
+  int kb = kmax;
+  int mb = mlength;
+  double c41[kmax];
+  double c51[kmax];
+  double c61[kmax];
+  double ct[kmax];
+  double *si = sim->s[i];
+  for (jj = 0; jj < jnum; jj++) {
+    jtype = tn[jj];
+    if (jtypes != nelements && jtypes != jtype && ktypes != nelements && ktypes != jtype) {
+    	expr[jj][0]=0;
+    	continue;
+    }
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq>rc*rc) {
+    	expr[jj][0]=0;
+    	continue;
+    }
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    if (!(m1>=1 && m1 <= res))pair->errorf(FLERR,"Neighbor list is invalid.");//usually results from nan somewhere.
+    r1 = r1-trunc(r1);
+    double *p0 = &expcuttable[(m1-1)*kmax];
+    double *p1 = &expcuttable[m1*kmax];
+    double *p2 = &expcuttable[(m1+1)*kmax];
+    double *p3 = &expcuttable[(m1+2)*kmax];
+    for (kk=0;kk<kmax;kk++) {
+    	expr[jj][kk] = p1[kk]+0.5*r1*(p2[kk]-p0[kk]+r1*(2.0*p0[kk]-5.0*p1[kk]+4.0*p2[kk]-p3[kk]+r1*(3.0*(p1[kk]-p2[kk])+p3[kk]-p0[kk])));
+    }
+    double* q = &dfctable[m1-1];
+    double* r2 = &rinvsqrttable[m1-1];
+    dfc[jj] = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    ri[jj] = r2[1] + 0.5 * r1*(r2[2] - r2[0] + r1*(2.0*r2[0] - 5.0*r2[1] + 4.0*r2[2] - r2[3] + r1*(3.0*(r2[1] - r2[2]) + r2[3] - r2[0])));
+    y[jj][0]=delx*ri[jj];
+    y[jj][1]=dely*ri[jj];
+    y[jj][2]=delz*ri[jj];
+  }
+  for (jj = 0; jj < jnum; jj++) {
+    if (expr[jj][0]==0)continue;
+    jtype = tn[jj];
+    if (jtypes != nelements && jtypes != jtype) {
+      continue;
+    }
+    j = jl[jj];
+    double *sj = sim->s[j];
+    double spj = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+    for (n = 0;n<kmax;n++) {
+      ct[n] = 2*alpha_k[n]/re;
+      c41[n]=(-ct[n]+2*dfc[jj])*y[jj][0];
+      c51[n]=(-ct[n]+2*dfc[jj])*y[jj][1];
+      c61[n]= (-ct[n]+2*dfc[jj])*y[jj][2];
+    }
+    if (jtypes==ktypes) {
+      for (kk=jj+1;kk< jnum; kk++) {
+        if (expr[kk][0]==0)continue;
+        ktype = tn[kk];
+        if (ktypes != nelements && ktypes != ktype) {
+          continue;
+        }
+        j = jl[kk];
+        double *sk = sim->s[j];
+        double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2];
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1  = 2*ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2  = 2*ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3  = 2*ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        double c10 = 2*ri[kk]*(y[jj][0]-dot*y[kk][0]);
+        double c11 = 2*ri[kk]*(y[jj][1]-dot*y[kk][1]);
+        double c12 = 2*ri[kk]*(y[jj][2]-dot*y[kk][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n];
+          double c5 = c51[n];
+          double c6 = c61[n];
+          double ct2 = -ct[n]+2*dfc[kk];
+          double c42 = ct2*y[kk][0];
+          double c52 = ct2*y[kk][1];
+          double c62 = ct2*y[kk][2];
+          //m=0
+          dspinx[jj*f+count] += 2*dot1*si[0]*spk;
+          dspiny[jj*f+count] += 2*dot1*si[1]*spk;
+          dspinz[jj*f+count] += 2*dot1*si[2]*spk;
+          dspinx[kk*f+count] += 2*dot1*si[0]*spj;
+          dspiny[kk*f+count] += 2*dot1*si[1]*spj;
+          dspinz[kk*f+count] += 2*dot1*si[2]*spj;
+          dspinx[jnum*f+count] += 2*dot1*(sj[0]*spk+sk[0]*spj);
+          dspiny[jnum*f+count] += 2*dot1*(sj[1]*spk+sk[1]*spj);
+          dspinz[jnum*f+count] += 2*dot1*(sj[2]*spk+sk[2]*spj);
+          features[count]+=2*dot1*spj*spk;
+          dfeaturesx[jj*f+count]+=dot1*c4*spj*spk;
+          dfeaturesy[jj*f+count]+=dot1*c5*spj*spk;
+          dfeaturesz[jj*f+count]+=dot1*c6*spj*spk;
+          dfeaturesx[kk*f+count]+=dot1*c42*spj*spk;
+          dfeaturesy[kk*f+count]+=dot1*c52*spj*spk;
+          dfeaturesz[kk*f+count]+=dot1*c62*spj*spk;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          c42*=dot;
+          c52*=dot;
+          c62*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7*spj*spk;
+            dfeaturesy[jj*f+count]+=c8*spj*spk;
+            dfeaturesz[jj*f+count]+=c9*spj*spk;
+            dfeaturesx[kk*f+count]+=dot1*(m*c10+c42)*spj*spk;
+            dfeaturesy[kk*f+count]+=dot1*(m*c11+c52)*spj*spk;
+            dfeaturesz[kk*f+count]+=dot1*(m*c12+c62)*spj*spk;
+            dot1*=dot;
+            dspinx[jj*f+count] += 2*dot1*si[0]*spk;
+            dspiny[jj*f+count] += 2*dot1*si[1]*spk;
+            dspinz[jj*f+count] += 2*dot1*si[2]*spk;
+            dspinx[kk*f+count] += 2*dot1*si[0]*spj;
+            dspiny[kk*f+count] += 2*dot1*si[1]*spj;
+            dspinz[kk*f+count] += 2*dot1*si[2]*spj;
+            dspinx[jnum*f+count] += 2*dot1*(sj[0]*spk+sk[0]*spj);
+            dspiny[jnum*f+count] += 2*dot1*(sj[1]*spk+sk[1]*spj);
+            dspinz[jnum*f+count] += 2*dot1*(sj[2]*spk+sk[2]*spj);
+            features[count++]+=2*dot1*spj*spk;
+          }
+        }
+      }
+      kk=jj;
+      if (ktypes == nelements || ktypes == jtype) {
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1 = 2*ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2 = 2*ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3 = 2*ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n];
+          double c5 = c51[n];
+          double c6 = c61[n];
+          //m=0
+          features[count]+=dot1*spj*spj;
+          dspinx[jj*f+count] += 2*dot1*si[0]*spj;
+          dspiny[jj*f+count] += 2*dot1*si[1]*spj;
+          dspinz[jj*f+count] += 2*dot1*si[2]*spj;
+          dspinx[jnum*f+count] += 2*dot1*sj[0]*spj;
+          dspiny[jnum*f+count] += 2*dot1*sj[1]*spj;
+          dspinz[jnum*f+count] += 2*dot1*sj[2]*spj;
+          dfeaturesx[jj*f+count]+=dot1*c4*spj*spj;
+          dfeaturesy[jj*f+count]+=dot1*c5*spj*spj;
+          dfeaturesz[jj*f+count]+=dot1*c6*spj*spj;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7*spj*spj;
+            dfeaturesy[jj*f+count]+=c8*spj*spj;
+            dfeaturesz[jj*f+count]+=c9*spj*spj;
+            dot1*=dot;
+            dspinx[jj*f+count] += 2*dot1*si[0]*spj;
+            dspiny[jj*f+count] += 2*dot1*si[1]*spj;
+            dspinz[jj*f+count] += 2*dot1*si[2]*spj;
+            dspinx[jnum*f+count] += 2*dot1*sj[0]*spj;
+            dspiny[jnum*f+count] += 2*dot1*sj[1]*spj;
+            dspinz[jnum*f+count] += 2*dot1*sj[2]*spj;
+            features[count++]+=dot1*spj*spj;
+          }
+        }
+      }
+    }
+    else {
+      for (kk=0;kk<jnum; kk++) {
+        if (expr[kk][0]==0)continue;
+        ktype = tn[kk];
+        if (ktypes != nelements && ktypes != ktype) {
+          continue;
+        }
+        j = jl[kk];
+        double *sk = sim->s[j];
+        double spk = si[0]*sk[0]+si[1]*sk[1]+si[2]*sk[2];
+        count = startingneuron;
+        double dot = (y[jj][0]*y[kk][0]+y[jj][1]*y[kk][1]+y[jj][2]*y[kk][2]);
+        double c1  = ri[jj]*(y[kk][0]-dot*y[jj][0]);
+        double c2  = ri[jj]*(y[kk][1]-dot*y[jj][1]);
+        double c3  = ri[jj]*(y[kk][2]-dot*y[jj][2]);
+        double c10 = ri[kk]*(y[jj][0]-dot*y[kk][0]);
+        double c11 = ri[kk]*(y[jj][1]-dot*y[kk][1]);
+        double c12 = ri[kk]*(y[jj][2]-dot*y[kk][2]);
+        for (n=0;n<kb;n++) {
+          double dot1=expr[jj][n]*expr[kk][n];
+          double c4 = c41[n]/2;
+          double c5 = c51[n]/2;
+          double c6 = c61[n]/2;
+          double ct2 = -ct[n]/2+dfc[kk];
+          double c42 = ct2*y[kk][0];
+          double c52 = ct2*y[kk][1];
+          double c62 = ct2*y[kk][2];
+          //m=0
+          features[count]+=dot1*spj*spk;
+          dspinx[jj*f+count] += dot1*si[0]*spk;
+          dspiny[jj*f+count] += dot1*si[1]*spk;
+          dspinz[jj*f+count] += dot1*si[2]*spk;
+          dspinx[kk*f+count] += dot1*si[0]*spj;
+          dspiny[kk*f+count] += dot1*si[1]*spj;
+          dspinz[kk*f+count] += dot1*si[2]*spj;
+          dspinx[jnum*f+count] += dot1*(sj[0]*spk+sk[0]*spj);
+          dspiny[jnum*f+count] += dot1*(sj[1]*spk+sk[1]*spj);
+          dspinz[jnum*f+count] += dot1*(sj[2]*spk+sk[2]*spj);
+          dfeaturesx[jj*f+count]+=dot1*c4*spj*spk;
+          dfeaturesy[jj*f+count]+=dot1*c5*spj*spk;
+          dfeaturesz[jj*f+count]+=dot1*c6*spj*spk;
+          dfeaturesx[kk*f+count]+=dot1*c42*spj*spk;
+          dfeaturesy[kk*f+count]+=dot1*c52*spj*spk;
+          dfeaturesz[kk*f+count]+=dot1*c62*spj*spk;
+          c4*=dot;
+          c5*=dot;
+          c6*=dot;
+          c42*=dot;
+          c52*=dot;
+          c62*=dot;
+          count++;
+          for (m=1;m<mb;m++) {
+            double c7 = dot1*(m*c1+c4);
+            double c8 = dot1*(m*c2+c5);
+            double c9 = dot1*(m*c3+c6);
+            dfeaturesx[jj*f+count]+=c7*spj*spk;
+            dfeaturesy[jj*f+count]+=c8*spj*spk;
+            dfeaturesz[jj*f+count]+=c9*spj*spk;
+            dfeaturesx[kk*f+count]+=dot1*(m*c10+c42)*spj*spk;
+            dfeaturesy[kk*f+count]+=dot1*(m*c11+c52)*spj*spk;
+            dfeaturesz[kk*f+count]+=dot1*(m*c12+c62)*spj*spk;
+            dot1*=dot;
+            dspinx[jj*f+count] += dot1*si[0]*spk;
+            dspiny[jj*f+count] += dot1*si[1]*spk;
+            dspinz[jj*f+count] += dot1*si[2]*spk;
+            dspinx[kk*f+count] += dot1*si[0]*spj;
+            dspiny[kk*f+count] += dot1*si[1]*spj;
+            dspinz[kk*f+count] += dot1*si[2]*spj;
+            dspinx[jnum*f+count] += dot1*(sj[0]*spk+sk[0]*spj);
+            dspiny[jnum*f+count] += dot1*(sj[1]*spk+sk[1]*spj);
+            dspinz[jnum*f+count] += dot1*(sj[2]*spk+sk[2]*spj);
+            features[count++]+=dot1*spj*spk;
+          }
+        }
+      }
+    }
+  }
+  for (jj=0;jj<jnum;jj++) {
+    if (expr[jj][0]==0) {continue;}
+    count = startingneuron;
+    for (n=0;n<kb;n++) {
+      for (m=0;m<mb;m++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+}
+

src/USER-RANN/rann_fingerprint_bondspin.h

@@ -0,0 +1,79 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(bondspin,Fingerprint_bondspin)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_BONDSPIN_H
+#define LMP_RANN_FINGERPRINT_BONDSPIN_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Fingerprint_bondspin : public Fingerprint {
+     public:
+      Fingerprint_bondspin(PairRANN *);
+      ~Fingerprint_bondspin();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);//spin
+      void do3bodyfeatureset_doubleneighborloop(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);
+      void do3bodyfeatureset_singleneighborloop(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);
+      void generate_exp_cut_table();
+      void generate_coefficients();
+      int get_length();
+
+      double *expcuttable;
+      double *dfctable;
+      double dr;
+      double *alpha_k;
+      double re;
+      int **coeff;
+      int **coeffx;
+      int **coeffy;
+      int **coeffz;
+      int kmax;
+      int mlength;
+      int **Mf;
+
+    };
+
+  }
+}
+
+
+#endif
+#endif /* FINGERPRINT_BOND_H_ */

src/USER-RANN/rann_fingerprint_radial.cpp

@@ -0,0 +1,239 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#include "rann_fingerprint_radial.h"
+
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_radial::Fingerprint_radial(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 2;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha = new double [1];
+  alpha[0] = -1;
+  nmax = 0;
+  omin = 0;
+  id = -1;
+  style = "radial";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  fullydefined = false;
+  _pair->allscreen = false;
+}
+
+Fingerprint_radial::~Fingerprint_radial()
+{
+  delete [] atomtypes;
+  delete [] radialtable;
+  delete [] alpha;
+  delete [] dfctable;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_radial::parse_values(std::string constant,std::vector<std::string> line1) {
+  int l;
+  int nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alpha")==0) {
+    delete [] alpha;
+    alpha = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("n")==0) {
+    nmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("o")==0) {
+    omin = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for radial power");
+  //code will run with default o=0 if o is never specified. All other values must be defined in potential file.
+  if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_radial::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alpha:\n",style,id);
+  for (i=0;i<(nmax-omin+1);i++) {
+    fprintf(fid,"%f ",alpha[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:o:\n",style,id);
+  fprintf(fid,"%d\n",omin);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:n:\n",style,id);
+  fprintf(fid,"%d\n",nmax);
+}
+
+//called after fingerprint is fully defined and tables can be computed.
+void Fingerprint_radial::allocate()
+{
+  int k,m;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  radialtable = new double [(res+buf)*get_length()];
+  dfctable = new double [res+buf];
+  for (k=0;k<(res+buf);k++) {
+    r1 = cutmax*cutmax*(double)(k)/(double)(res);
+    for (m=0;m<=(nmax-omin);m++) {
+      radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+        dfctable[k]=0;
+    }
+    else{
+    dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+  generate_rinvssqrttable();
+}
+
+//called after fingerprint is declared for i-j type, but before its parameters are read.
+void Fingerprint_radial::init(int *i,int _id)
+{
+  empty = false;
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  id = _id;
+}
+
+void Fingerprint_radial::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl)
+{
+  int nelements = pair->nelements;
+  int res = pair->res;
+  int i,j,jj,itype,jtype,l;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  //
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int count=0;
+//  double **x = sim->x;
+  int *type = sim->type;
+  ilist = sim->ilist;
+  double cutmax = pair->cutmax;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double cutinv2 = 1/cutmax/cutmax;
+  //loop over neighbors
+  for (jj = 0; jj < jnum; jj++) {
+    jtype =tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue;
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq > rc*rc)continue;
+    count = startingneuron;
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");}
+    if (radialtable[m1]==0) {continue;}
+    //cubic interpolation from tables
+    double *p1 = &radialtable[m1*(nmax-omin+1)];
+    double *p2 = &radialtable[(m1+1)*(nmax-omin+1)];
+    double *p3 = &radialtable[(m1+2)*(nmax-omin+1)];
+    double *p0 = &radialtable[(m1-1)*(nmax-omin+1)];
+    double *q = &dfctable[m1-1];
+    double *rinvs = &rinvsqrttable[m1-1];
+    r1 = r1-trunc(r1);
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0])));
+    for (l=0;l<=(nmax-omin);l++) {
+      double rt = p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l])));
+      features[count]+=rt;
+      rt *= (l+omin)/rsq+(-alpha[l]/re+dfc)*ri;
+      //update neighbor's features
+      dfeaturesx[jj*f+count]+=rt*delx;
+      dfeaturesy[jj*f+count]+=rt*dely;
+      dfeaturesz[jj*f+count]+=rt*delz;
+      //update atom's features
+      dfeaturesx[jnum*f+count]-=rt*delx;
+      dfeaturesy[jnum*f+count]-=rt*dely;
+      dfeaturesz[jnum*f+count]-=rt*delz;
+      count++;
+    }
+  }
+}
+
+int Fingerprint_radial::get_length()
+{
+  return nmax-omin+1;
+}

src/USER-RANN/rann_fingerprint_radial.h

@@ -0,0 +1,70 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(radial,Fingerprint_radial)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_RADIAL_H
+#define LMP_RANN_FINGERPRINT_RADIAL_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+
+    class Fingerprint_radial : public Fingerprint {
+     public:
+      Fingerprint_radial(PairRANN *);
+      ~Fingerprint_radial();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);
+      int get_length();
+
+      double *radialtable;
+      double *dfctable;
+      double dr;
+      double *alpha;
+      double re;
+      int nmax;//highest term
+      int omin;//lowest term
+
+    };
+
+  }
+}
+
+#endif
+#endif /* FINGERPRINT_RADIAL_H_ */

src/USER-RANN/rann_fingerprint_radialscreened.cpp

@@ -0,0 +1,253 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#include "rann_fingerprint_radialscreened.h"
+
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_radialscreened::Fingerprint_radialscreened(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 2;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha = new double [1];
+  alpha[0] = -1;
+  nmax = 0;
+  omin = 0;
+  id = -1;
+  style = "radialscreened";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  fullydefined = false;
+  _pair->doscreen = true;
+  screen = true;
+}
+
+Fingerprint_radialscreened::~Fingerprint_radialscreened()
+{
+  delete [] atomtypes;
+  delete [] radialtable;
+  delete [] alpha;
+  delete [] dfctable;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_radialscreened::parse_values(std::string constant,std::vector<std::string> line1) {
+  int l;
+  int nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alpha")==0) {
+    delete [] alpha;
+    alpha = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("n")==0) {
+    nmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("o")==0) {
+    omin = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for radial power");
+  //code will run with default o=0 if o is never specified. All other values must be defined in potential file.
+  if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_radialscreened::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alpha:\n",style,id);
+  for (i=0;i<(nmax-omin+1);i++) {
+    fprintf(fid,"%f ",alpha[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:o:\n",style,id);
+  fprintf(fid,"%d\n",omin);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:n:\n",style,id);
+  fprintf(fid,"%d\n",nmax);
+}
+
+//called after fingerprint is fully defined and tables can be computed.
+void Fingerprint_radialscreened::allocate()
+{
+  int k,m;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  radialtable = new double [(res+buf)*get_length()];
+  dfctable = new double [res+buf];
+  for (k=0;k<(res+buf);k++) {
+    r1 = cutmax*cutmax*(double)(k)/(double)(res);
+    for (m=0;m<=(nmax-omin);m++) {
+      radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+        dfctable[k]=0;
+    }
+    else{
+    dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+  generate_rinvssqrttable();
+}
+
+//called after fingerprint is declared for i-j type, but before its parameters are read.
+void Fingerprint_radialscreened::init(int *i,int _id)
+{
+  empty = false;
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  id = _id;
+}
+
+void Fingerprint_radialscreened::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl)
+{
+    int nelements = pair->nelements;
+    int res = pair->res;
+    int i,j,jj,itype,jtype,l,kk;
+    double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+    int *ilist,*jlist,*numneigh,**firstneigh;
+    //
+    PairRANN::Simulation *sim = &pair->sims[sid];
+    int count=0;
+    double **x = sim->x;
+    int *type = sim->type;
+    ilist = sim->ilist;
+    double cutmax = pair->cutmax;
+    i = ilist[ii];
+    itype = pair->map[type[i]];
+    int f = pair->net[itype].dimensions[0];
+    double cutinv2 = 1/cutmax/cutmax;
+    //loop over neighbors
+    for (jj = 0; jj < jnum; jj++) {
+      if (Bij[jj]==false) {continue;}
+      jtype = tn[jj];
+      if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue;
+      delx = xn[jj];
+      dely = yn[jj];
+      delz = zn[jj];
+      rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq > rc*rc)continue;
+      count = startingneuron;
+      double r1 = (rsq*((double)res)*cutinv2);
+      int m1 = (int)r1;
+      if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");}
+      if (radialtable[m1]==0) {continue;}
+      //cubic interpolation from tables
+      double *p1 = &radialtable[m1*(nmax-omin+1)];
+      double *p2 = &radialtable[(m1+1)*(nmax-omin+1)];
+      double *p3 = &radialtable[(m1+2)*(nmax-omin+1)];
+      double *p0 = &radialtable[(m1-1)*(nmax-omin+1)];
+      double *q = &dfctable[m1-1];
+      double *rinvs = &rinvsqrttable[m1-1];
+      r1 = r1-trunc(r1);
+      double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+      double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0])));
+      for (l=0;l<=(nmax-omin);l++) {
+        double rt = Sik[jj]*(p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l]))));
+        features[count]+=rt;
+        double rt1 = rt*((l+omin)/rsq+(-alpha[l]/re+dfc)*ri);
+        //update neighbor's features
+        dfeaturesx[jj*f+count]+=rt1*delx+rt*dSikx[jj];
+        dfeaturesy[jj*f+count]+=rt1*dely+rt*dSiky[jj];
+        dfeaturesz[jj*f+count]+=rt1*delz+rt*dSikz[jj];
+        for (kk=0;kk<jnum;kk++) {
+          if (Bij[kk]==false) {continue;}
+          dfeaturesx[kk*f+count]+=rt*dSijkx[jj*jnum+kk];
+          dfeaturesy[kk*f+count]+=rt*dSijky[jj*jnum+kk];
+          dfeaturesz[kk*f+count]+=rt*dSijkz[jj*jnum+kk];
+        }
+        count++;
+      }
+    }
+    for (jj=0;jj<jnum;jj++) {
+      if (Bij[jj]==false) {continue;}
+      count = startingneuron;
+      for (l=0;l<=(nmax-omin);l++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+
+int Fingerprint_radialscreened::get_length()
+{
+  return nmax-omin+1;
+}

src/USER-RANN/rann_fingerprint_radialscreened.h

@@ -0,0 +1,72 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(radialscreened,Fingerprint_radialscreened)
+
+#else
+#ifndef LMP_RANN_FINGERPRINT_RADIALSCREENED_H
+#define LMP_RANN_FINGERPRINT_RADIALSCREENED_H
+
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Fingerprint_radialscreened : public Fingerprint {
+     public:
+      Fingerprint_radialscreened(PairRANN *);
+      ~Fingerprint_radialscreened();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);
+      int get_length();
+
+      double *radialtable;
+      double *dfctable;
+      double dr;
+      double *alpha;
+      double re;
+      int nmax;//highest term
+      int omin;//lowest term
+
+    };
+}
+
+}
+
+#endif
+
+
+
+#endif /* FINGERPRINT_RADIALSCREENED_H_ */

src/USER-RANN/rann_fingerprint_radialscreenedspin.cpp

@@ -0,0 +1,264 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#include "rann_fingerprint_radialscreenedspin.h"
+
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_radialscreenedspin::Fingerprint_radialscreenedspin(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 2;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha = new double [1];
+  alpha[0] = -1;
+  nmax = 0;
+  omin = 0;
+  id = -1;
+  style = "radialscreenedspin";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  fullydefined = false;
+  _pair->doscreen = true;
+  screen = true;
+  _pair->dospin = true;
+  spin = true;
+}
+
+Fingerprint_radialscreenedspin::~Fingerprint_radialscreenedspin()
+{
+  delete [] atomtypes;
+  delete [] radialtable;
+  delete [] alpha;
+  delete [] dfctable;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_radialscreenedspin::parse_values(std::string constant,std::vector<std::string> line1) {
+  int l;
+  int nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alpha")==0) {
+    delete [] alpha;
+    alpha = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("n")==0) {
+    nmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("o")==0) {
+    omin = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for radial power");
+  //code will run with default o=0 if o is never specified. All other values must be defined in potential file.
+  if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_radialscreenedspin::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alpha:\n",style,id);
+  for (i=0;i<(nmax-omin+1);i++) {
+    fprintf(fid,"%f ",alpha[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:o:\n",style,id);
+  fprintf(fid,"%d\n",omin);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:n:\n",style,id);
+  fprintf(fid,"%d\n",nmax);
+}
+
+//called after fingerprint is fully defined and tables can be computed.
+void Fingerprint_radialscreenedspin::allocate()
+{
+  int k,m;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  radialtable = new double [(res+buf)*get_length()];
+  dfctable = new double [res+buf];
+  for (k=0;k<(res+buf);k++) {
+    r1 = cutmax*cutmax*(double)(k)/(double)(res);
+    for (m=0;m<=(nmax-omin);m++) {
+      radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+        dfctable[k]=0;
+    }
+    else{
+    dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+  generate_rinvssqrttable();
+}
+
+//called after fingerprint is declared for i-j type, but before its parameters are read.
+void Fingerprint_radialscreenedspin::init(int *i,int _id)
+{
+  empty = false;
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  id = _id;
+}
+
+void Fingerprint_radialscreenedspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,double *Sik, double *dSikx, double*dSiky, double *dSikz, double *dSijkx, double *dSijky, double *dSijkz, bool *Bij,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl)
+{
+    int nelements = pair->nelements;
+    int res = pair->res;
+    int i,j,jj,itype,jtype,l,kk;
+    double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+    int *ilist,*jlist,*numneigh,**firstneigh;
+    PairRANN::Simulation *sim = &pair->sims[sid];
+    int count=0;
+    double **x = sim->x;
+    int *type = sim->type;
+    ilist = sim->ilist;
+    double cutmax = pair->cutmax;
+    i = ilist[ii];
+    itype = pair->map[type[i]];
+    int f = pair->net[itype].dimensions[0];
+    double cutinv2 = 1/cutmax/cutmax;
+    double *si = sim->s[i];
+    //loop over neighbors
+    for (jj = 0; jj < jnum; jj++) {
+      if (Bij[jj]==false) {continue;}
+      jtype = tn[jj];
+      if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue;
+      delx = xn[jj];
+      dely = yn[jj];
+      delz = zn[jj];
+      rsq = delx*delx + dely*dely + delz*delz;
+      if (rsq > rc*rc)continue;
+      count = startingneuron;
+      double r1 = (rsq*((double)res)*cutinv2);
+      int m1 = (int)r1;
+      if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");}
+      if (radialtable[m1]==0) {continue;}
+      j=jl[jj];
+      double *sj = sim->s[j];
+      double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+      //cubic interpolation from tables
+      double *p1 = &radialtable[m1*(nmax-omin+1)];
+      double *p2 = &radialtable[(m1+1)*(nmax-omin+1)];
+      double *p3 = &radialtable[(m1+2)*(nmax-omin+1)];
+      double *p0 = &radialtable[(m1-1)*(nmax-omin+1)];
+      double *q = &dfctable[m1-1];
+      double *rinvs = &rinvsqrttable[m1-1];
+      r1 = r1-trunc(r1);
+      double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+      double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0])));
+      for (l=0;l<=(nmax-omin);l++) {
+        double rt = Sik[jj]*(p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l]))));
+        //update neighbor's features
+        dspinx[jj*f+count]+=rt*si[0];
+        dspiny[jj*f+count]+=rt*si[1];
+        dspinz[jj*f+count]+=rt*si[2];
+        dspinx[jnum*f+count]+=rt*sj[0];
+        dspiny[jnum*f+count]+=rt*sj[1];
+        dspinz[jnum*f+count]+=rt*sj[2];
+        rt *= sp;
+        features[count]+=rt;
+        double rt1 = rt*((l+omin)/rsq+(-alpha[l]/re+dfc)*ri);
+        dfeaturesx[jj*f+count]+=rt1*delx+rt*dSikx[jj];
+        dfeaturesy[jj*f+count]+=rt1*dely+rt*dSiky[jj];
+        dfeaturesz[jj*f+count]+=rt1*delz+rt*dSikz[jj];
+        for (kk=0;kk<jnum;kk++) {
+          if (Bij[kk]==false) {continue;}
+          dfeaturesx[kk*f+count]+=rt*dSijkx[jj*jnum+kk];
+          dfeaturesy[kk*f+count]+=rt*dSijky[jj*jnum+kk];
+          dfeaturesz[kk*f+count]+=rt*dSijkz[jj*jnum+kk];
+        }
+        count++;
+      }
+    }
+    for (jj=0;jj<jnum;jj++) {
+      if (Bij[jj]==false) {continue;}
+      count = startingneuron;
+      for (l=0;l<=(nmax-omin);l++) {
+        dfeaturesx[jnum*f+count]-=dfeaturesx[jj*f+count];
+        dfeaturesy[jnum*f+count]-=dfeaturesy[jj*f+count];
+        dfeaturesz[jnum*f+count]-=dfeaturesz[jj*f+count];
+        count++;
+      }
+    }
+  }
+
+int Fingerprint_radialscreenedspin::get_length()
+{
+  return nmax-omin+1;
+}

src/USER-RANN/rann_fingerprint_radialscreenedspin.h

@@ -0,0 +1,72 @@
+
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(radialscreenedspin,Fingerprint_radialscreenedspin)
+
+#else
+#ifndef LMP_RANN_FINGERPRINT_RADIALSCREENEDSPIN_H
+#define LMP_RANN_FINGERPRINT_RADIALSCREENEDSPIN_H
+
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Fingerprint_radialscreenedspin : public Fingerprint {
+     public:
+      Fingerprint_radialscreenedspin(PairRANN *);
+      ~Fingerprint_radialscreenedspin();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      virtual void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,double*,bool*,int,int,double*,double*,double*,int*,int,int*);//spin,screen
+      int get_length();
+
+      double *radialtable;
+      double *dfctable;
+      double dr;
+      double *alpha;
+      double re;
+      int nmax;//highest term
+      int omin;//lowest term
+
+    };
+  }
+
+}
+
+#endif
+
+
+
+#endif /* FINGERPRINT_RADIALSCREENED_H_ */

src/USER-RANN/rann_fingerprint_radialspin.cpp

@@ -0,0 +1,252 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+#include "rann_fingerprint_radialspin.h"
+
+
+
+using namespace LAMMPS_NS::RANN;
+
+Fingerprint_radialspin::Fingerprint_radialspin(PairRANN *_pair) : Fingerprint(_pair)
+{
+  n_body_type = 2;
+  dr = 0;
+  re = 0;
+  rc = 0;
+  alpha = new double [1];
+  alpha[0] = -1;
+  nmax = 0;
+  omin = 0;
+  id = -1;
+  style = "radialspin";
+  atomtypes = new int [n_body_type];
+  empty = true;
+  fullydefined = false;
+  _pair->allscreen = false;
+  _pair->dospin = true;
+  spin = true;
+}
+
+Fingerprint_radialspin::~Fingerprint_radialspin()
+{
+  delete [] atomtypes;
+  delete [] radialtable;
+  delete [] alpha;
+  delete [] dfctable;
+  delete [] rinvsqrttable;
+}
+
+bool Fingerprint_radialspin::parse_values(std::string constant,std::vector<std::string> line1) {
+  int l;
+  int nwords=line1.size();
+  if (constant.compare("re")==0) {
+    re = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("rc")==0) {
+    rc = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("alpha")==0) {
+    delete [] alpha;
+    alpha = new double [nwords];
+    for (l=0;l<nwords;l++) {
+      alpha[l]=strtod(line1[l].c_str(),NULL);
+    }
+  }
+  else if (constant.compare("dr")==0) {
+    dr = strtod(line1[0].c_str(),NULL);
+  }
+  else if (constant.compare("n")==0) {
+    nmax = strtol(line1[0].c_str(),NULL,10);
+  }
+  else if (constant.compare("o")==0) {
+    omin = strtol(line1[0].c_str(),NULL,10);
+  }
+  else pair->errorf(FLERR,"Undefined value for radial power");
+  //code will run with default o=0 if o is never specified. All other values must be defined in potential file.
+  if (re!=0 && rc!=0 && alpha!=0 && dr!=0 && nmax!=0)return true;
+  return false;
+}
+
+void Fingerprint_radialspin::write_values(FILE *fid) {
+  int i;
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:re:\n",style,id);
+  fprintf(fid,"%f\n",re);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:rc:\n",style,id);
+  fprintf(fid,"%f\n",rc);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:alpha:\n",style,id);
+  for (i=0;i<(nmax-omin+1);i++) {
+    fprintf(fid,"%f ",alpha[i]);
+  }
+  fprintf(fid,"\n");
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:dr:\n",style,id);
+  fprintf(fid,"%f\n",dr);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:o:\n",style,id);
+  fprintf(fid,"%d\n",omin);
+  fprintf(fid,"fingerprintconstants:");
+  fprintf(fid,"%s",pair->elementsp[atomtypes[0]]);
+  for (i=1;i<n_body_type;i++) {
+    fprintf(fid,"_%s",pair->elementsp[atomtypes[i]]);
+  }
+  fprintf(fid,":%s_%d:n:\n",style,id);
+  fprintf(fid,"%d\n",nmax);
+}
+
+//called after fingerprint is fully defined and tables can be computed.
+void Fingerprint_radialspin::allocate()
+{
+  int k,m;
+  double r1;
+  int buf = 5;
+  int res = pair->res;
+  double cutmax = pair->cutmax;
+  radialtable = new double [(res+buf)*get_length()];
+  dfctable = new double [res+buf];
+  for (k=0;k<(res+buf);k++) {
+    r1 = cutmax*cutmax*(double)(k)/(double)(res);
+    for (m=0;m<=(nmax-omin);m++) {
+      radialtable[k*(nmax-omin+1)+m]=pow(sqrt(r1)/re,m+omin)*exp(-alpha[m]*sqrt(r1)/re)*cutofffunction(sqrt(r1),rc,dr);
+    }
+    if (sqrt(r1)>=rc || sqrt(r1) <= (rc-dr)) {
+        dfctable[k]=0;
+    }
+    else{
+    dfctable[k]=-8*pow(1-(rc-sqrt(r1))/dr,3)/dr/(1-pow(1-(rc-sqrt(r1))/dr,4));
+    }
+  }
+  generate_rinvssqrttable();
+}
+
+//called after fingerprint is declared for i-j type, but before its parameters are read.
+void Fingerprint_radialspin::init(int *i,int _id)
+{
+  empty = false;
+  for (int j=0;j<n_body_type;j++) {atomtypes[j] = i[j];}
+  id = _id;
+}
+
+void Fingerprint_radialspin::compute_fingerprint(double * features,double * dfeaturesx,double *dfeaturesy,double *dfeaturesz,double * dspinx,double *dspiny,double *dspinz,int ii,int sid,double *xn,double *yn,double*zn,int *tn,int jnum,int *jl)
+{
+  int nelements = pair->nelements;
+  int res = pair->res;
+  int i,j,jj,itype,jtype,l;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  //
+  PairRANN::Simulation *sim = &pair->sims[sid];
+  int count=0;
+  int *type = sim->type;
+  ilist = sim->ilist;
+  double cutmax = pair->cutmax;
+  i = ilist[ii];
+  itype = pair->map[type[i]];
+  int f = pair->net[itype].dimensions[0];
+  double cutinv2 = 1/cutmax/cutmax;
+  double *si = sim->s[i];
+  firstneigh = sim->firstneigh;
+  jlist = firstneigh[i];
+  //loop over neighbors
+  for (jj = 0; jj < jnum; jj++) {
+    j = jl[jj];
+    jtype =tn[jj];
+    if (atomtypes[1] != nelements && atomtypes[1] != jtype)continue;
+    delx = xn[jj];
+    dely = yn[jj];
+    delz = zn[jj];
+    rsq = delx*delx + dely*dely + delz*delz;
+    if (rsq > rc*rc)continue;
+    count = startingneuron;
+    double r1 = (rsq*((double)res)*cutinv2);
+    int m1 = (int)r1;
+    if (m1>res || m1<1) {pair->errorf(FLERR,"invalid neighbor radius!");}
+    if (radialtable[m1]==0) {continue;}
+    double *sj = sim->s[j];
+    double sp = si[0]*sj[0]+si[1]*sj[1]+si[2]*sj[2];
+    //cubic interpolation from tables
+    double *p1 = &radialtable[m1*(nmax-omin+1)];
+    double *p2 = &radialtable[(m1+1)*(nmax-omin+1)];
+    double *p3 = &radialtable[(m1+2)*(nmax-omin+1)];
+    double *p0 = &radialtable[(m1-1)*(nmax-omin+1)];
+    double *q = &dfctable[m1-1];
+    double *rinvs = &rinvsqrttable[m1-1];
+    r1 = r1-trunc(r1);
+    double dfc = q[1] + 0.5 * r1*(q[2] - q[0] + r1*(2.0*q[0] - 5.0*q[1] + 4.0*q[2] - q[3] + r1*(3.0*(q[1] - q[2]) + q[3] - q[0])));
+    double ri = rinvs[1] + 0.5 * r1*(rinvs[2] - rinvs[0] + r1*(2.0*rinvs[0] - 5.0*rinvs[1] + 4.0*rinvs[2] - rinvs[3] + r1*(3.0*(rinvs[1] - rinvs[2]) + rinvs[3] - rinvs[0])));
+    for (l=0;l<=(nmax-omin);l++) {
+      double rt = p1[l]+0.5*r1*(p2[l]-p0[l]+r1*(2.0*p0[l]-5.0*p1[l]+4.0*p2[l]-p3[l]+r1*(3.0*(p1[l]-p2[l])+p3[l]-p0[l])));
+      dspinx[jj*f+count]+=rt*si[0];
+      dspiny[jj*f+count]+=rt*si[1];
+      dspinz[jj*f+count]+=rt*si[2];
+      dspinx[jnum*f+count]+=rt*sj[0];
+      dspiny[jnum*f+count]+=rt*sj[1];
+      dspinz[jnum*f+count]+=rt*sj[2];
+      rt *= sp;
+      features[count]+=rt;
+      rt *= (l+omin)/rsq+(-alpha[l]/re+dfc)*ri;
+      //update neighbor's features
+      dfeaturesx[jj*f+count]+=rt*delx;
+      dfeaturesy[jj*f+count]+=rt*dely;
+      dfeaturesz[jj*f+count]+=rt*delz;
+      //update atom's features
+      dfeaturesx[jnum*f+count]-=rt*delx;
+      dfeaturesy[jnum*f+count]-=rt*dely;
+      dfeaturesz[jnum*f+count]-=rt*delz;
+      count++;
+    }
+  }
+}
+
+int Fingerprint_radialspin::get_length()
+{
+  return nmax-omin+1;
+}

src/USER-RANN/rann_fingerprint_radialspin.h

@@ -0,0 +1,69 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+/*  ----------------------------------------------------------------------
+   Contributing authors: Christopher Barrett (MSU) barrett@me.msstate.edu
+                              Doyl Dickel (MSU) doyl@me.msstate.edu
+    ----------------------------------------------------------------------*/
+/*
+“The research described and the resulting data presented herein, unless
+otherwise noted, was funded under PE 0602784A, Project T53 "Military
+Engineering Applied Research", Task 002 under Contract No. W56HZV-17-C-0095,
+managed by the U.S. Army Combat Capabilities Development Command (CCDC) and
+the Engineer Research and Development Center (ERDC).  The work described in
+this document was conducted at CAVS, MSU.  Permission was granted by ERDC
+to publish this information. Any opinions, findings and conclusions or
+recommendations expressed in this material are those of the author(s) and
+do not necessarily reflect the views of the United States Army.​”
+
+DISTRIBUTION A. Approved for public release; distribution unlimited. OPSEC#4918
+ */
+
+#ifdef FINGERPRINT_CLASS
+
+FingerprintStyle(radialspin,Fingerprint_radialspin)
+
+#else
+
+#ifndef LMP_RANN_FINGERPRINT_RADIALSPIN_H
+#define LMP_RANN_FINGERPRINT_RADIALSPIN_H
+
+#include "rann_fingerprint.h"
+
+namespace LAMMPS_NS {
+  namespace RANN {
+    class Fingerprint_radialspin : public Fingerprint {
+     public:
+      Fingerprint_radialspin(PairRANN *);
+      ~Fingerprint_radialspin();
+      bool parse_values(std::string,std::vector<std::string>);
+      void write_values(FILE *);
+      void init(int*,int);
+      void allocate();
+      void compute_fingerprint(double*,double*,double*,double*,double*,double*,double*,int,int,double*,double*,double*,int*,int,int*);
+      int get_length();
+
+      double *radialtable;
+      double *dfctable;
+      double dr;
+      double *alpha;
+      double re;
+      int nmax;//highest term
+      int omin;//lowest term
+
+    };
+  }
+
+}
+
+#endif
+#endif /* FINGERPRINT_RADIAL_H_ */

src/USER-RANN/rann_list_activation.h

@@ -0,0 +1,2 @@
+#include "rann_activation_linear.h"
+#include "rann_activation_sig_i.h"

src/USER-RANN/rann_list_fingerprint.h

@@ -0,0 +1,8 @@
+#include "rann_fingerprint_bond.h"
+#include "rann_fingerprint_bondscreened.h"
+#include "rann_fingerprint_bondscreenedspin.h"
+#include "rann_fingerprint_bondspin.h"
+#include "rann_fingerprint_radial.h"
+#include "rann_fingerprint_radialscreened.h"
+#include "rann_fingerprint_radialscreenedspin.h"
+#include "rann_fingerprint_radialspin.h"