// clang-format off
/* -*- c++ -*- ----------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   https://www.lammps.org/ 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 "pair_rann.h"

#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_special.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "tokenizer.h"
#include "update.h"

#include <cmath>
#include <cstring>

#include "rann_activation_linear.h"
#include "rann_activation_sig_i.h"
#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"

#define MAXLINE 1024

using namespace LAMMPS_NS;

static const char cite_ml_rann_package[] =
  "ML-RANN package:\n\n"
  "@Article{Nitol2021,\n"
  " author = {Nitol, Mashroor S and Dickel, Doyl E and Barrett, Christopher D},\n"
  " title = {Artificial neural network potential for pure zinc},\n"
  " journal = {Computational Materials Science},\n"
  " year =    2021,\n"
  " volume =  188,\n"
  " pages =   {110207}\n"
  "}\n\n";


PairRANN::PairRANN(LAMMPS *lmp) : Pair(lmp)
{
  if (lmp->citeme) lmp->citeme->add(cite_ml_rann_package);

  //initialize ints and bools
  single_enable = 0;
  restartinfo = 0;
  one_coeff = 1;
  manybody_flag = 1;
  allocated = 0;
  nelements = -1;
  nelementsp = -1;
  comm_forward = 0;
  comm_reverse = 0;
  res = 10000;
  cutmax = 0;
  dospin = false;
  memguess = 0;
  nmax1 = 0;
  nmax2 = 0;
  fmax = 0;
  fnmax = 0;
  //at least one of the following two flags will change during fingerprint definition:
  doscreen = false;
  allscreen = true;

  //null init for arrays with sizes not yet determined.
  elements = nullptr;
  mass = nullptr;
  elementsp = nullptr;
  map  = nullptr;
  fingerprintcount = nullptr;
  fingerprintlength = nullptr;
  fingerprintperelement = nullptr;
  screening_min = nullptr;
  screening_max = nullptr;
  weightdefined = nullptr;
  biasdefined = nullptr;
  xn = nullptr;
  yn = nullptr;
  zn = nullptr;
  Sik = nullptr;
  dSikx = nullptr;
  dSiky = nullptr;
  dSikz = nullptr;
  dSijkx = nullptr;
  dSijky = nullptr;
  dSijkz = nullptr;
  sx = nullptr;
  sy = nullptr;
  sz = nullptr;
  dSijkxc = nullptr;
  dSijkyc = nullptr;
  dSijkzc = nullptr;
  dfeaturesx = nullptr;
  dfeaturesy = nullptr;
  dfeaturesz = nullptr;
  features = nullptr;
  layer = nullptr;
  sum = nullptr;
  sum1 = nullptr;
  dlayerx = nullptr;
  dlayery = nullptr;
  dlayerz = nullptr;
  dlayersumx = nullptr;
  dlayersumy = nullptr;
  dlayersumz = nullptr;
  dsx = nullptr;
  dsy = nullptr;
  dsz = nullptr;
  dssumx = nullptr;
  dssumy = nullptr;
  dssumz = nullptr;
  tn = nullptr;
  jl = nullptr;
  Bij = nullptr;
  sims = nullptr;
  net = nullptr;
  activation = nullptr;
  fingerprints = nullptr;
}

PairRANN::~PairRANN()
{
  deallocate();
}

void PairRANN::deallocate()
{
  //clear memory
  delete[] mass;
  for (int i=0;i<nelements;i++) {delete [] elements[i];}
  delete[] elements;
  for (int i=0;i<nelementsp;i++) {delete [] elementsp[i];}
  delete[] elementsp;
  for (int i=0;i<=nelements;i++) {
    if (net[i].layers>0) {
      for (int j=0;j<net[i].layers-1;j++) {
        delete[] net[i].Weights[j];
        delete[] net[i].Biases[j];
        delete activation[i][j];
      }
      delete[] activation[i];
      delete[] net[i].Weights;
      delete[] net[i].Biases;
      delete[] weightdefined[i];
      delete[] biasdefined[i];
    }
    delete[] net[i].dimensions;
  }
  delete[] net;
  delete[] map;
  for (int i=0;i<nelementsp;i++) {
    if (fingerprintlength[i]>0) {
      for (int j=0;j<fingerprintperelement[i];j++) {
        delete fingerprints[i][j];
      }
      delete[] fingerprints[i];
    }
  }
  delete[] fingerprints;
  delete[] activation;
  delete[] weightdefined;
  delete[] biasdefined;
  delete[] fingerprintcount;
  delete[] fingerprintperelement;
  delete[] fingerprintlength;
  delete[] screening_min;
  delete[] screening_max;
  memory->destroy(xn);
  memory->destroy(yn);
  memory->destroy(zn);
  memory->destroy(tn);
  memory->destroy(jl);
  memory->destroy(features);
  memory->destroy(dfeaturesx);
  memory->destroy(dfeaturesy);
  memory->destroy(dfeaturesz);
  memory->destroy(layer);
  memory->destroy(sum);
  memory->destroy(sum1);
  memory->destroy(dlayerx);
  memory->destroy(dlayery);
  memory->destroy(dlayerz);
  memory->destroy(dlayersumx);
  memory->destroy(dlayersumy);
  memory->destroy(dlayersumz);
  memory->destroy(Sik);
  memory->destroy(Bij);
  memory->destroy(dSikx);
  memory->destroy(dSiky);
  memory->destroy(dSikz);
  memory->destroy(dSijkx);
  memory->destroy(dSijky);
  memory->destroy(dSijkz);
  memory->destroy(dSijkxc);
  memory->destroy(dSijkyc);
  memory->destroy(dSijkzc);
  memory->destroy(sx);
  memory->destroy(sy);
  memory->destroy(sz);
  memory->destroy(dsx);
  memory->destroy(dsy);
  memory->destroy(dsz);
  memory->destroy(dssumx);
  memory->destroy(dssumy);
  memory->destroy(dssumz);
  memory->destroy(setflag);
  memory->destroy(cutsq);
}

void PairRANN::allocate(const std::vector<std::string> &elementwords)
{
  int i,n;
  n = atom->ntypes;
  memory->create(setflag,n+1,n+1,"pair:setflag");
  memory->create(cutsq,n+1,n+1,"pair:cutsq");
  cutmax = 0;
  nmax1 = 100;
  nmax2 = 20;
  fmax = 0;
  fnmax = 0;
  nelementsp=nelements+1;
  //initialize arrays
  elements = new char *[nelements];
  elementsp = new char *[nelementsp];//elements + 'all'
  mass = new double[nelements];
  net = new NNarchitecture[nelementsp];
  weightdefined = new bool*[nelementsp];
  biasdefined = new bool *[nelementsp];
  activation = new RANN::Activation**[nelementsp];
  fingerprints = new RANN::Fingerprint**[nelementsp];
  fingerprintlength = new int[nelementsp];
  fingerprintperelement = new int[nelementsp];
  fingerprintcount = new int[nelementsp];
  screening_min = new double[nelements*nelements*nelements];
  screening_max = new double[nelements*nelements*nelements];
  for (i=0;i<nelements;i++) {
    for (int j =0;j<nelements;j++) {
      for (int k=0;k<nelements;k++) {
        screening_min[i*nelements*nelements+j*nelements+k] = 0.8;//default values. Custom values may be read from potential file later.
        screening_max[i*nelements*nelements+j*nelements+k] = 2.8;//default values. Custom values may be read from potential file later.
      }
    }
  }
  for (i=0;i<=nelements;i++) {
    fingerprintlength[i]=0;
    fingerprintperelement[i] = -1;
    fingerprintcount[i] = 0;
    if (i<nelements) {
      mass[i]=-1.0;
      elements[i]= utils::strdup(elementwords[i]);
    }
    elementsp[i] = utils::strdup(elementwords[i]);
    net[i].layers = 0;
    net[i].dimensions = new int[1];
    net[i].dimensions[0]=0;
  }
}

void PairRANN::settings(int narg, char ** /*arg*/)
{
  //read pair_style command in input file
  if (narg > 0) error->one(FLERR,"Illegal pair_style command");
}

void PairRANN::coeff(int narg, char **arg)
{
  int i,j;
  deallocate();//clear allocation from any previous coeff
  map = new int[atom->ntypes+1];
  if (narg != 3 + atom->ntypes) error->one(FLERR,"Incorrect args for pair coefficients");
  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->one(FLERR,"Incorrect args for pair coefficients");
  nelements = -1;
  read_file(arg[2]);
  // read args that map atom types to elements in potential file
  // map[i] = which element the Ith atom type is, -1 if NULL
  for (i = 3; i < narg; i++) {
    if (strcmp(arg[i],"NULL") == 0) {
      map[i-2] = -1;
      continue;
    }
    for (j = 0; j < nelements; j++) {
      if (strcmp(arg[i],elements[j]) == 0) break;
    }
    if (j < nelements) map[i-2] = j;
    else error->one(FLERR,"No matching element in NN potential file");
  }
  // clear setflag since coeff() called once with I,J = * *
  int n = atom->ntypes;
  for (i = 1; i <= n; i++) {
    for (j = i; j <= n; j++) {
      setflag[i][j] = 0;
    }
  }
  // set setflag i,j for type pairs where both are mapped to elements
  // set mass of atom type if i = j
  int count = 0;
  for (i = 1; i <= n; i++) {
    for (j = i; j <= n; j++) {
      if (map[i] >= 0 && map[j] >= 0) {
        setflag[i][j] = 1;
        if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
        count++;
      }
    }
  }
  if (count == 0) error->one(FLERR,"Incorrect args for pair coefficients");
  for (i=0;i<nelementsp;i++) {
    for (j=0;j<fingerprintperelement[i];j++) {
      fingerprints[i][j]->allocate();
    }
  }
  allocated=1;
}

void PairRANN::read_file(char *filename)
{
  FILE *fp;
  int eof = 0;
  std::string line,line1;
  const int longline = 4096;
  int linenum=0;
  char linetemp[longline];
  std::string strtemp;
  char *ptr;
  std::vector<std::string> linev,line1v;
  fp = utils::open_potential(filename,lmp,nullptr);
  if (fp == nullptr) {error->one(FLERR,"Cannot open RANN potential file");}
  ptr=fgets(linetemp,longline,fp);
  linenum++;
  strtemp=utils::trim_comment(linetemp);
  while (strtemp.empty()) {
          ptr=fgets(linetemp,longline,fp);
          strtemp=utils::trim_comment(linetemp);
          linenum++;
  }
  line=strtemp;
  while (eof == 0) {
    ptr=fgets(linetemp,longline,fp);
    linenum++;
    if (ptr == nullptr) {
      fclose(fp);
      if (check_potential()) {
        error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
      }
      else {
        eof = 1;
        break;
      }
    }
    strtemp=utils::trim_comment(linetemp);
    while (strtemp.empty()) {
        ptr=fgets(linetemp,longline,fp);
        strtemp=utils::trim_comment(linetemp);
        linenum++;
    }
    line1=linetemp;
    Tokenizer values = Tokenizer(line,": ,\t_\n");
    Tokenizer values1 = Tokenizer(line1,": ,\t_\n");
    linev = values.as_vector();
    line1v = values1.as_vector();
    if (linev[0]=="atomtypes") read_atom_types(line1v,filename,linenum);
    else if (linev[0]=="mass") read_mass(linev,line1v,filename,linenum);
    else if (linev[0]=="fingerprintsperelement") read_fpe(linev,line1v,filename,linenum);
    else if (linev[0]=="fingerprints") read_fingerprints(linev,line1v,filename,linenum);
    else if (linev[0]=="fingerprintconstants") read_fingerprint_constants(linev,line1v,filename,linenum);
    else if (linev[0]=="networklayers") read_network_layers(linev,line1v,filename,linenum);
    else if (linev[0]=="layersize") read_layer_size(linev,line1v,filename,linenum);
    else if (linev[0]=="weight") read_weight(linev,line1v,fp,filename,&linenum);
    else if (linev[0]=="bias") read_bias(linev,line1v,fp,filename,&linenum);
    else if (linev[0]=="activationfunctions") read_activation_functions(linev,line1v,filename,linenum);
    else if (linev[0]=="screening") read_screening(linev,line1v,filename,linenum);
    else if (linev[0]=="calibrationparameters") continue;//information on how the network was trained
    else error->one(FLERR,"Could not understand file syntax: unknown keyword");
    ptr=fgets(linetemp,longline,fp);
    linenum++;
    strtemp=utils::trim_comment(linetemp);
    while (strtemp.empty()) {
        ptr=fgets(linetemp,longline,fp);
        strtemp=utils::trim_comment(linetemp);
        linenum++;
    }
    if (ptr == nullptr) {
      if (check_potential()) {
        error->one(FLERR,"Invalid syntax in potential file, values are inconsistent or missing");
      }
      else {
        eof = 1;
        break;
      }
    }
    line=linetemp;
  }
}

void PairRANN::read_atom_types(std::vector<std::string> line,char *filename,int linenum) {
  int nwords = line.size();
  if (nwords < 1) error->one(filename,linenum,"Incorrect syntax for atom types");
  nelements = nwords;
  line.push_back("all");
  allocate(line);
}

void PairRANN::read_mass(const std::vector<std::string> &line1, const std::vector<std::string> &line2, const char *filename,int linenum) {
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before mass in potential file.");
  for (int i=0;i<nelements;i++) {
    if (line1[1].compare(elements[i])==0) {
      mass[i]=utils::numeric(filename,linenum,line2[0].c_str(),1,lmp);
      return;
    }
  }
  error->one(filename,linenum-1,"mass element not found in atom types.");
}

void PairRANN::read_fpe(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i;
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints per element in potential file.");
  for (i=0;i<nelementsp;i++) {
    if (line[1].compare(elementsp[i])==0) {
      fingerprintperelement[i] = utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
      fingerprints[i] = new RANN::Fingerprint *[fingerprintperelement[i]];
      for (int j=0;j<fingerprintperelement[i];j++) {
        fingerprints[i][j]=new RANN::Fingerprint(this);
      }
      return;
    }
  }
  error->one(filename,linenum-1,"fingerprint-per-element element not found in atom types");
}

void PairRANN::read_fingerprints(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int nwords1,nwords,i,j,k,i1,*atomtypes;
  bool found;
  nwords1 = line1.size();
  nwords = line.size();
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
  atomtypes = new int[nwords-1];
  for (i=1;i<nwords;i++) {
    found = false;
    for (j=0;j<nelementsp;j++) {
      if (line[i].compare(elementsp[j])==0) {
        atomtypes[i-1]=j;
        found = true;
        break;
      }
    }
    if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
  }
  i = atomtypes[0];
  k = 0;
  if (fingerprintperelement[i]==-1) {error->one(filename,linenum-1,"fingerprint per element must be defined before fingerprints");}
  while (k<nwords1) {
    i1 = fingerprintcount[i];
    if (i1>=fingerprintperelement[i]) {error->one(filename,linenum,"more fingerprints defined than fingerprint per element");}
    delete fingerprints[i][i1];
    fingerprints[i][i1] = create_fingerprint(line1[k].c_str());
    if (fingerprints[i][i1]->n_body_type!=nwords-1) {error->one(filename,linenum,"invalid fingerprint for element combination");}
    k++;
    fingerprints[i][i1]->init(atomtypes,utils::inumeric(filename,linenum,line1[k++].c_str(),1,lmp));
    fingerprintcount[i]++;
  }
  delete[] atomtypes;
}

void PairRANN::read_fingerprint_constants(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i,j,k,i1,*atomtypes;
  bool found;
  int nwords = line.size();
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
  int n_body_type = nwords-4;
  atomtypes = new int[n_body_type];
  for (i=1;i<=n_body_type;i++) {
    found = false;
    for (j=0;j<nelementsp;j++) {
      if (line[i].compare(elementsp[j])==0) {
        atomtypes[i-1]=j;
        found = true;
        break;
      }
    }
    if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
  }
  i = atomtypes[0];
  found = false;
  for (k=0;k<fingerprintperelement[i];k++) {
    if (fingerprints[i][k]->empty) {continue;}
    if (n_body_type!=fingerprints[i][k]->n_body_type) {continue;}
    for (j=0;j<n_body_type;j++) {
      if (fingerprints[i][k]->atomtypes[j]!=atomtypes[j]) {break;}
      if (j==n_body_type-1) {
        if (line[nwords-3].compare(fingerprints[i][k]->style)==0 && utils::inumeric(filename,linenum,line[nwords-2].c_str(),1,lmp)==fingerprints[i][k]->id) {
          found=true;
          i1 = k;
          break;
        }
      }
    }
    if (found) {break;}
  }
  if (!found) {error->one(filename,linenum-1,"cannot define constants for unknown fingerprint");}
  fingerprints[i][i1]->fullydefined=fingerprints[i][i1]->parse_values(line[nwords-1],line1);
  delete[] atomtypes;
}

void PairRANN::read_network_layers(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i,j;
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before network layers in potential file.");
  for (i=0;i<nelements;i++) {
    if (line[1].compare(elements[i])==0) {
      net[i].layers = utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
      if (net[i].layers < 1)error->one(filename,linenum,"invalid number of network layers");
      delete[] net[i].dimensions;
      weightdefined[i] = new bool [net[i].layers];
      biasdefined[i] = new bool [net[i].layers];
      net[i].dimensions = new int[net[i].layers];
      net[i].Weights = new double * [net[i].layers-1];
      net[i].Biases = new double * [net[i].layers-1];
      for (j=0;j<net[i].layers;j++) {
        net[i].dimensions[j]=0;
        weightdefined[i][j] = false;
        biasdefined[i][j] = false;
      }
      activation[i]=new RANN::Activation* [net[i].layers-1];
      for (int j=0;j<net[i].layers-1;j++) {
        activation[i][j]= new RANN::Activation(this);
      }
      return;
    }
  }
  error->one(filename,linenum-1,"network layers element not found in atom types");
}

void PairRANN::read_layer_size(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i;
  for (i=0;i<nelements;i++) {
    if (line[1].compare(elements[i])==0) {
      if (net[i].layers==0)error->one(filename,linenum-1,"networklayers for each atom type must be defined before the corresponding layer sizes.");
      int j = utils::inumeric(filename,linenum,line[2].c_str(),1,lmp);
      if (j>=net[i].layers || j<0) {error->one(filename,linenum,"invalid layer in layer size definition");};
      net[i].dimensions[j]= utils::inumeric(filename,linenum,line1[0].c_str(),1,lmp);
      return;
    }
  }
  error->one(filename,linenum-1,"layer size element not found in atom types");
}

void PairRANN::read_weight(std::vector<std::string> line,std::vector<std::string> line1,FILE* fp,char *filename,int *linenum) {
  int i,j,k,l,nwords;
  char *ptr;
  const int longline = 4096;
  char linetemp [longline];
  for (l=0;l<nelements;l++) {
    if (line[1].compare(elements[l])==0) {
      if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before weights.");
      i=utils::inumeric(filename,*linenum,line[2].c_str(),1,lmp);
      if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid weight layer");
      if (net[l].dimensions[i]==0 || net[l].dimensions[i+1]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding weight");
      net[l].Weights[i] = new double[net[l].dimensions[i]*net[l].dimensions[i+1]];
      weightdefined[l][i] = true;
      nwords = line1.size();
      if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
      for (k=0;k<net[l].dimensions[i];k++) {
        net[l].Weights[i][k] = utils::numeric(filename,*linenum,line1[k].c_str(),1,lmp);
      }
      for (j=1;j<net[l].dimensions[i+1];j++) {
        ptr = fgets(linetemp,longline,fp);
        (*linenum)++;
        Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
        line1 = values1.as_vector();
        if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
        nwords = line1.size();
        if (nwords != net[l].dimensions[i])error->one(filename,*linenum,"invalid weights per line");
        for (k=0;k<net[l].dimensions[i];k++) {
          net[l].Weights[i][j*net[l].dimensions[i]+k] = utils::numeric(filename,*linenum,line1[k].c_str(),1,lmp);
        }
      }
      return;
    }
  }
  error->one(filename,*linenum-1,"weight element not found in atom types");
}

void PairRANN::read_bias(std::vector<std::string> line,std::vector<std::string> line1,FILE* fp,char *filename,int *linenum) {
  int i,j,l;
  char linetemp[MAXLINE],*ptr;
  for (l=0;l<nelements;l++) {
    if (line[1].compare(elements[l])==0) {
      if (net[l].layers==0)error->one(filename,*linenum-1,"networklayers must be defined before biases.");
      i=utils::inumeric(filename,*linenum,line[2].c_str(),1,lmp);
      if (i>=net[l].layers || i<0)error->one(filename,*linenum-1,"invalid bias layer");
      if (net[l].dimensions[i]==0) error->one(filename,*linenum-1,"network layer sizes must be defined before corresponding bias");
      biasdefined[l][i] = true;
      net[l].Biases[i] = new double[net[l].dimensions[i+1]];
      net[l].Biases[i][0] = utils::numeric(filename,*linenum,line1[0].c_str(),1,lmp);
      for (j=1;j<net[l].dimensions[i+1];j++) {
        ptr=fgets(linetemp,MAXLINE,fp);
        if (ptr==nullptr)error->one(filename,*linenum,"unexpected end of potential file!");
        (*linenum)++;
        Tokenizer values1 = Tokenizer(linetemp,": ,\t_\n");
        line1 = values1.as_vector();
        net[l].Biases[i][j] = utils::numeric(filename,*linenum,line1[0].c_str(),1,lmp);
      }
      return;
    }
  }
  error->one(filename,*linenum-1,"bias element not found in atom types");
}

void PairRANN::read_activation_functions(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i,l;
  for (l=0;l<nelements;l++) {
    if (line[1].compare(elements[l])==0) {
      if (net[l].layers==0)error->one(filename,linenum-1,"networklayers must be defined before activation functions.");
      i = strtol(line[2].c_str(),nullptr,10);
      if (i>=net[l].layers || i<0)error->one(filename,linenum-1,"invalid activation layer");
      delete activation[l][i];
      activation[l][i]=create_activation(line1[0].c_str());
      return;
    }
  }
  error->one(filename,linenum-1,"activation function element not found in atom types");
}

void PairRANN::read_screening(std::vector<std::string> line,std::vector<std::string> line1,char *filename,int linenum) {
  int i,j,k,*atomtypes;
  bool found;
  int nwords = line.size();
  if (nelements == -1)error->one(filename,linenum-1,"atom types must be defined before fingerprints in potential file.");
  if (nwords!=5)error->one(filename,linenum-1,"invalid screening command");
  int n_body_type = 3;
  atomtypes = new int[n_body_type];
  for (i=1;i<=n_body_type;i++) {
    found = false;
    for (j=0;j<nelementsp;j++) {
      if (line[i].compare(elementsp[j])==0) {
        atomtypes[i-1]=j;
        found = true;
        break;
      }
    }
    if (!found) {error->one(filename,linenum-1,"fingerprint element not found in atom types");}
  }
  i = atomtypes[0];
  j = atomtypes[1];
  k = atomtypes[2];
  int index = i*nelements*nelements+j*nelements+k;
  if (line[4].compare("Cmin")==0)  {
    screening_min[index] = utils::numeric(filename,linenum,line1[0].c_str(),1,lmp);
  }
  else if (line[4].compare("Cmax")==0) {
    screening_max[index] = utils::numeric(filename,linenum,line1[0].c_str(),1,lmp);
  }
  else error->one(filename,linenum-1,"unrecognized screening keyword");
  delete[] atomtypes;
}

//Called after finishing reading the potential file to make sure it is complete. True is bad.
//also does the rest of the memory allocation.
bool PairRANN::check_potential() {
  int i,j,k,l;
  if (nelements==-1) {return true;}
  for (i=0;i<=nelements;i++) {
    if (i<nelements) {
      if (mass[i]<0)return true;//uninitialized mass
    }
    if (net[i].layers==0)break;//no definitions for this starting element, not considered an error.
    net[i].maxlayer=0;
    for (j=0;j<net[i].layers;j++) {
      if (net[i].dimensions[j]==0)return true;//incomplete network definition
      if (net[i].dimensions[j]>net[i].maxlayer)net[i].maxlayer = net[i].dimensions[j];
    }
    if (net[i].maxlayer>fnmax) {fnmax = net[i].maxlayer;}
    if (net[i].dimensions[net[i].layers-1]!=1)return true;//output layer must have single neuron (the energy)
    if (net[i].dimensions[0]>fmax)fmax=net[i].dimensions[0];
    for (j=0;j<net[i].layers-1;j++) {
      if (!weightdefined[i][j])return true;//undefined weights
      if (!biasdefined[i][j])return true;//undefined biases
      if (activation[i][j]->empty)return true;//undefined activations
      for (k=0;k<net[i].dimensions[j+1];k++) {
        for (l=0;l<net[i].dimensions[j];l++) {
          if (net[i].Weights[j][k*net[i].dimensions[j]+l]==0)return true;//undefined weights
        }
        if (net[i].Biases[j][k]==0)return true;//undefined biases
      }
    }
    for (j=0;j<fingerprintperelement[i];j++) {
      if (fingerprints[i][j]->fullydefined==false)return true;
      fingerprints[i][j]->startingneuron = fingerprintlength[i];
      fingerprintlength[i] +=fingerprints[i][j]->get_length();
      if (fingerprints[i][j]->rc>cutmax) {cutmax = fingerprints[i][j]->rc;}
    }
    if (net[i].dimensions[0]!=fingerprintlength[i])return true;
  }
  memory->create(xn,nmax1,"pair:xn");
  memory->create(yn,nmax1,"pair:yn");
  memory->create(zn,nmax1,"pair:zn");
  memory->create(tn,nmax1,"pair:tn");
  memory->create(jl,nmax1,"pair:jl");
  memory->create(features,fmax,"pair:features");
  memory->create(dfeaturesx,fmax*nmax2,"pair:dfeaturesx");
  memory->create(dfeaturesy,fmax*nmax2,"pair:dfeaturesy");
  memory->create(dfeaturesz,fmax*nmax2,"pair:dfeaturesz");
  memory->create(layer,fnmax,"pair:layer");
  memory->create(sum,fnmax,"pair:sum");
  memory->create(sum1,fnmax,"pair:sum1");
  memory->create(dlayerx,nmax2,fnmax,"pair:dlayerx");
  memory->create(dlayery,nmax2,fnmax,"pair:dlayery");
  memory->create(dlayerz,nmax2,fnmax,"pair:dlayerz");
  memory->create(dlayersumx,nmax2,fnmax,"pair:dlayersumx");
  memory->create(dlayersumy,nmax2,fnmax,"pair:dlayersumy");
  memory->create(dlayersumz,nmax2,fnmax,"pair:dlayersumz");
  if (doscreen) {
    memory->create(Sik,nmax2,"pair:Sik");
    memory->create(Bij,nmax2,"pair:Bij");
    memory->create(dSikx,nmax2,"pair:dSikx");
    memory->create(dSiky,nmax2,"pair:dSiky");
    memory->create(dSikz,nmax2,"pair:dSikz");
    memory->create(dSijkx,nmax2*nmax2,"pair:dSijkx");
    memory->create(dSijky,nmax2*nmax2,"pair:dSijky");
    memory->create(dSijkz,nmax2*nmax2,"pair:dSijkz");
    memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc");
    memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc");
    memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc");
  }
  if (dospin) {
    memory->create(sx,fmax*nmax2,"pair:sx");
    memory->create(sy,fmax*nmax2,"pair:sy");
    memory->create(sz,fmax*nmax2,"pair:sz");
    memory->create(dsx,nmax2,fnmax,"pair:dsx");
    memory->create(dsy,nmax2,fnmax,"pair:dsy");
    memory->create(dsz,nmax2,fnmax,"pair:dsz");
    memory->create(dssumx,nmax2,fnmax,"pair:dssumx");
    memory->create(dssumy,nmax2,fnmax,"pair:dssumy");
    memory->create(dssumz,nmax2,fnmax,"pair:dssumz");
  }
  return false;//everything looks good
}

void PairRANN::compute(int eflag, int vflag)
{
  //perform force/energy computation_
  if (dospin) {
    if (strcmp(update->unit_style,"metal") != 0)
      error->one(FLERR,"Spin pair styles require metal units");
    if (!atom->sp_flag)
        error->one(FLERR,"Spin pair styles requires atom/spin style");
  }
  ev_init(eflag,vflag);

  // only global virial via fdotr is supported by this pair style

  if (vflag_atom)
    error->all(FLERR,"Pair style rann does not support computing per-atom stress");
  if (vflag && !vflag_fdotr)
    error->all(FLERR,"Pair style rann does not support 'pair_modify nofdotr'");

  int ii,i,j;
  int nn = 0;
  sims = new Simulation[1];
  sims->inum = listfull->inum;
  sims->ilist=listfull->ilist;
  sims->id = listfull->ilist;
  sims->type = atom->type;
  sims->x = atom->x;
  sims->numneigh = listfull->numneigh;
  sims->firstneigh = listfull->firstneigh;
  if (dospin) {
    sims->s = atom->sp;
  }
  int itype,f,jnum,len;
  if (eflag_global) {eng_vdwl=0;eng_coul=0;}
  double energy=0;
  double **force = atom->f;
  double **fm = atom->fm;

  //loop over atoms
  for (ii=0;ii<sims->inum;ii++) {
      i = sims->ilist[ii];
      itype = map[sims->type[i]];
      f = net[itype].dimensions[0];
      jnum = sims->numneigh[i];
      if (jnum>nmax1) {
        nmax1 = jnum;
        memory->grow(xn,nmax1,"pair:xn");
        memory->grow(yn,nmax1,"pair:yn");
        memory->grow(zn,nmax1,"pair:zn");
        memory->grow(tn,nmax1,"pair:tn");
        memory->grow(jl,nmax1,"pair:jl");
      }
      cull_neighbor_list(&jnum,i,0);
      if (jnum>nmax2) {
        nmax2=jnum;
        memory->grow(dfeaturesx,fmax*nmax2,"pair:dfeaturesx");
        memory->grow(dfeaturesy,fmax*nmax2,"pair:dfeaturesy");
        memory->grow(dfeaturesz,fmax*nmax2,"pair:dfeaturesz");
        memory->grow(layer,fnmax,"pair:layer");
        memory->grow(sum,fnmax,"pair:sum");
        memory->grow(sum1,fnmax,"pair:sum1");
        memory->grow(dlayerx,nmax2,fnmax,"pair:dlayerx");
        memory->grow(dlayery,nmax2,fnmax,"pair:dlayery");
        memory->grow(dlayerz,nmax2,fnmax,"pair:dlayerz");
        memory->grow(dlayersumx,nmax2,fnmax,"pair:dlayersumx");
        memory->grow(dlayersumy,nmax2,fnmax,"pair:dlayersumy");
        memory->grow(dlayersumz,nmax2,fnmax,"pair:dlayersumz");
        if (doscreen) {
          memory->grow(Sik,nmax2,"pair:Sik");
          memory->grow(Bij,nmax2,"pair:Bij");
          memory->grow(dSikx,nmax2,"pair:dSikx");
          memory->grow(dSiky,nmax2,"pair:dSiky");
          memory->grow(dSikz,nmax2,"pair:dSikz");
          memory->grow(dSijkx,nmax2*nmax2,"pair:dSijkx");
          memory->grow(dSijky,nmax2*nmax2,"pair:dSijky");
          memory->grow(dSijkz,nmax2*nmax2,"pair:dSijkz");
          memory->destroy(dSijkxc);
          memory->destroy(dSijkyc);
          memory->destroy(dSijkzc);
          memory->create(dSijkxc,nmax2,nmax2,"pair:dSijkxc");
          memory->create(dSijkyc,nmax2,nmax2,"pair:dSijkyc");
          memory->create(dSijkzc,nmax2,nmax2,"pair:dSijkzc");
        }
        if (dospin) {
          memory->grow(sx,fmax*nmax2,"pair:sx");
          memory->grow(sy,fmax*nmax2,"pair:sy");
          memory->grow(sz,fmax*nmax2,"pair:sz");
          memory->grow(dsx,nmax2,fnmax,"pair:dsx");
          memory->grow(dsy,nmax2,fnmax,"pair:dsy");
          memory->grow(dsz,nmax2,fnmax,"pair:dsz");
          memory->grow(dssumx,nmax2,fnmax,"pair:dssumx");
          memory->grow(dssumy,nmax2,fnmax,"pair:dssumy");
          memory->grow(dssumz,nmax2,fnmax,"pair:dssumz");
        }
      }
      for (j=0;j<f;j++) {
        features[j]=0;
      }
      for (j=0;j<f*jnum;j++) {
        dfeaturesx[j]=dfeaturesy[j]=dfeaturesz[j]=0;
      }
      //screening is calculated once for all atoms if any fingerprint uses it.
      if (dospin) {
        for (j=0;j<f*jnum;j++) {
          sx[j]=sy[j]=sz[j]=0;
        }
      }
      if (doscreen) screening(ii,0,jnum-1);
      if (allscreen) screen_neighbor_list(&jnum);
      //do fingerprints for atom type
      len = fingerprintperelement[itype];
      for (j=0;j<len;j++) {
        if      (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==true  && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==true  && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
      }
      itype = nelements;
      //do fingerprints for type "all"
      len = fingerprintperelement[itype];
      for (j=0;j<len;j++) {
        if      (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==false && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==true  && fingerprints[itype][j]->screen==false)fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,ii,nn,xn,yn,zn,tn,jnum-1,jl);
        else if (fingerprints[itype][j]->spin==true  && fingerprints[itype][j]->screen==true) fingerprints[itype][j]->compute_fingerprint(features,dfeaturesx,dfeaturesy,dfeaturesz,sx,sy,sz,Sik,dSikx,dSiky,dSikz,dSijkx,dSijky,dSijkz,Bij,ii,nn,xn,yn,zn,tn,jnum-1,jl);
      }
      //run fingerprints through network
      if (dospin) {
        propagateforwardspin(&energy,force,fm,ii,jnum);
      } else {
        propagateforward(&energy,force,ii,jnum);
      }
  }
  if (vflag_fdotr) virial_fdotr_compute();
  delete[] sims;
}

void PairRANN::cull_neighbor_list(int* jnum,int i,int sn) {
  int *jlist,j,count,jj,*type,jtype;
  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
  double **x = sims[sn].x;
  xtmp = x[i][0];
  ytmp = x[i][1];
  ztmp = x[i][2];
  type = sims[sn].type;
  jlist = sims[sn].firstneigh[i];
  count = 0;
  for (jj=0;jj<jnum[0];jj++) {
    j = jlist[jj];
    j &= NEIGHMASK;
    jtype = map[type[j]];
    delx = xtmp - x[j][0];
    dely = ytmp - x[j][1];
    delz = ztmp - x[j][2];
    rsq = delx*delx + dely*dely + delz*delz;
    if (rsq>cutmax*cutmax) {
      continue;
    }
    xn[count]=delx;
    yn[count]=dely;
    zn[count]=delz;
    tn[count]=jtype;
    jl[count]=j;
    count++;
  }
  jnum[0]=count+1;
}

void PairRANN::screen_neighbor_list(int *jnum) {
  int jj,kk,count,count1;
  count = 0;
  for (jj=0;jj<jnum[0]-1;jj++) {
    if (Bij[jj]) {
      count1 = 0;
      if (jj!=count) {
      xn[count]=xn[jj];
      yn[count]=yn[jj];
      zn[count]=zn[jj];
      tn[count]=tn[jj];
      jl[count]=jl[jj];
      Sik[count]=Sik[jj];
      dSikx[count]=dSikx[jj];
      dSiky[count]=dSiky[jj];
      dSikz[count]=dSikz[jj];
      }
      for (kk=0;kk<jnum[0]-1;kk++) {
        if (Bij[kk]) {
          dSijkxc[count][count1] = dSijkx[jj*(jnum[0]-1)+kk];
          dSijkyc[count][count1] = dSijky[jj*(jnum[0]-1)+kk];
          dSijkzc[count][count1] = dSijkz[jj*(jnum[0]-1)+kk];
          count1++;
        }
      }
      count++;
    }
  }
  jnum[0]=count+1;
  for (jj=0;jj<count;jj++) {
    Bij[jj]=true;
    for (kk=0;kk<count;kk++) {
      dSijkx[jj*count+kk] = dSijkxc[jj][kk];
      dSijky[jj*count+kk] = dSijkyc[jj][kk];
      dSijkz[jj*count+kk] = dSijkzc[jj][kk];
    }
  }
}


void PairRANN::screening(int ii,int sid,int jnum)
{
  //see Baskes, Materials Chemistry and Physics 50 (1997) 152-1.58
  int i,jj,kk,itype,jtype,ktype;
  double Sijk,Cijk,Cn,Cd,Dij,Dik,Djk,C,dfc,dC;
  PairRANN::Simulation *sim = &sims[sid];
  double delx,dely,delz,rij,delx2,dely2,delz2,rik,delx3,dely3,delz3,rjk;
  i = sim->ilist[ii];
  itype = map[sim->type[i]];
  for (int jj=0;jj<jnum;jj++) {
    Sik[jj]=1;
    Bij[jj]=true;
    dSikx[jj]=0;
    dSiky[jj]=0;
    dSikz[jj]=0;
    for (kk=0;kk<jnum;kk++) {
      dSijkx[jj*jnum+kk]=0;
      dSijky[jj*jnum+kk]=0;
      dSijkz[jj*jnum+kk]=0;
    }
  }
  for (kk=0;kk<jnum;kk++) {//outer sum over k in accordance with source, some others reorder to outer sum over jj
    if (Bij[kk]==false) {continue;}
    ktype = tn[kk];
    delx2 = xn[kk];
    dely2 = yn[kk];
    delz2 = zn[kk];
    rik = delx2*delx2+dely2*dely2+delz2*delz2;
    if (rik>cutmax*cutmax) {
      Bij[kk]= false;
      continue;
    }
    for (jj=0;jj<jnum;jj++) {
      if (jj==kk) {continue;}
      if (Bij[jj]==false) {continue;}
      jtype = tn[jj];
      delx = xn[jj];
      dely = yn[jj];
      delz = zn[jj];
      rij = delx*delx+dely*dely+delz*delz;
      if (rij>cutmax*cutmax) {
        Bij[jj] = false;
        continue;
      }
      delx3 = delx2-delx;
      dely3 = dely2-dely;
      delz3 = delz2-delz;
      rjk = delx3*delx3+dely3*dely3+delz3*delz3;
      if (rik+rjk<=rij) {continue;}//bond angle > 90 degrees
      if (rik+rij<=rjk) {continue;}//bond angle > 90 degrees
      double Cmax = screening_max[itype*nelements*nelements+jtype*nelements+ktype];
      double Cmin = screening_min[itype*nelements*nelements+jtype*nelements+ktype];
      double temp1 = rij-rik+rjk;
      Cn = temp1*temp1-4*rij*rjk;
      //a few expanded computations provided for clarity:
      //Cn = (rij-rik+rjk)*(rij-rik+rjk)-4*rij*rjk;
      //Cd = (rij-rjk)*(rij-rjk)-rik*rik;
      //Cijk = 1+2*(rik*rij+rik*rjk-rik*rik)/(rik*rik-(rij-rjk)*(rij-rjk));
      temp1 = rij-rjk;
      Cd = temp1*temp1-rik*rik;
      Cijk = Cn/Cd;
      C = (Cijk-Cmin)/(Cmax-Cmin);
      if (C>=1) {continue;}
      else if (C<=0) {
        Bij[kk]=false;
        break;
      }
      dC = Cmax-Cmin;
      dC *= dC;
      dC *= dC;
      temp1 = 1-C;
      temp1 *= temp1;
      temp1 *= temp1;
      Sijk = 1-temp1;
      Sijk *= Sijk;
      Dij = 4*rik*(Cn+4*rjk*(rij+rik-rjk))/Cd/Cd;
      Dik = -4*(rij*Cn+rjk*Cn+8*rij*rik*rjk)/Cd/Cd;
      Djk = 4*rik*(Cn+4*rij*(rik-rij+rjk))/Cd/Cd;
      temp1 = Cijk-Cmax;
      double temp2 = temp1*temp1;
      dfc = 8*temp1*temp2/(temp2*temp2-dC);
      Sik[kk] *= Sijk;
      dSijkx[kk*jnum+jj] = dfc*(delx*Dij-delx3*Djk);
      dSikx[kk] += dfc*(delx2*Dik+delx3*Djk);
      dSijky[kk*jnum+jj] = dfc*(dely*Dij-dely3*Djk);
      dSiky[kk] += dfc*(dely2*Dik+dely3*Djk);
      dSijkz[kk*jnum+jj] = dfc*(delz*Dij-delz3*Djk);
      dSikz[kk] += dfc*(delz2*Dik+delz3*Djk);
    }
  }
}


//Called by getproperties. Propagate features and dfeatures through network. Updates force and energy
void PairRANN::propagateforward(double *energy,double **force,int ii,int jnum) {
  int i,j,k,jj,j1,itype,i1;
  int *ilist;
  ilist = listfull->ilist;
  int *type = atom->type;
  i1=ilist[ii];
  itype = map[type[i1]];
  NNarchitecture net1 = net[itype];
  int L = net1.layers-1;
  //energy output with forces from analytical derivatives
  double dsum1;
  int f = net1.dimensions[0];
  for (i=0;i<net1.layers-1;i++) {
    for (j=0;j<net1.dimensions[i+1];j++) {
      //energy forward propagation
      sum[j]=0;
      for (k=0;k<net1.dimensions[i];k++) {
        if (i==0&&j==0) {
          layer[k]=features[k];
        }
        sum[j] += net1.Weights[i][j*net1.dimensions[i]+k]*layer[k];
      }
      sum[j] += net1.Biases[i][j];
      dsum1 = activation[itype][i]->dactivation_function(sum[j]);
      sum[j] = activation[itype][i]->activation_function(sum[j]);
      if (i==L-1) {
        energy[j] = sum[j];
        if (eflag_atom) eatom[i1]=sum[j];
        if (eflag_global) eng_vdwl +=sum[j];
      }
      //force propagation
      for (jj=0;jj<jnum;jj++) {
        dlayersumx[jj][j]=0;
        dlayersumy[jj][j]=0;
        dlayersumz[jj][j]=0;
        for (k=0;k<net1.dimensions[i];k++) {
          if (i==0&&j==0) {
            dlayerx[jj][k]=dfeaturesx[jj*f+k];
            dlayery[jj][k]=dfeaturesy[jj*f+k];
            dlayerz[jj][k]=dfeaturesz[jj*f+k];
          }
          double w1 = net1.Weights[i][j*net1.dimensions[i]+k];
          dlayersumx[jj][j] += w1*dlayerx[jj][k];
          dlayersumy[jj][j] += w1*dlayery[jj][k];
          dlayersumz[jj][j] += w1*dlayerz[jj][k];
        }
        dlayersumx[jj][j] = dsum1*dlayersumx[jj][j];
        dlayersumy[jj][j] = dsum1*dlayersumy[jj][j];
        dlayersumz[jj][j] = dsum1*dlayersumz[jj][j];
        if (i==L-1 && jj < (jnum-1)) {
          j1 = jl[jj];
          force[j1][0]+=dlayersumx[jj][j];
          force[j1][1]+=dlayersumy[jj][j];
          force[j1][2]+=dlayersumz[jj][j];
        }
      }
      if (i==L-1) {
        j1 = ilist[ii];
        jj = jnum-1;
        force[j1][0]+=dlayersumx[jj][j];
        force[j1][1]+=dlayersumy[jj][j];
        force[j1][2]+=dlayersumz[jj][j];
      }
    }
    //update values for next iteration
    for (j=0;j<net1.dimensions[i+1];j++) {
      layer[j]=sum[j];
      for (jj=0;jj<jnum;jj++) {
        dlayerx[jj][j] = dlayersumx[jj][j];
        dlayery[jj][j] = dlayersumy[jj][j];
        dlayerz[jj][j] = dlayersumz[jj][j];
      }
    }
  }
}

//Called by getproperties. Propagate features and dfeatures through network. Updates force and energy
void PairRANN::propagateforwardspin(double * energy,double **force,double **fm,int ii,int jnum) {
  int i,j,k,jj,j1,itype,i1;
  int *ilist;
  ilist = listfull->ilist;
  int *type = atom->type;
  i1=ilist[ii];
  itype = map[type[i1]];
  NNarchitecture net1 = net[itype];
  int L = net1.layers-1;
  //energy output with forces from analytical derivatives
  double dsum1;
  int f = net1.dimensions[0];
  for (i=0;i<net1.layers-1;i++) {
    for (j=0;j<net1.dimensions[i+1];j++) {
      //energy forward propagation
      sum[j]=0;
      for (k=0;k<net1.dimensions[i];k++) {
        if (i==0&&j==0) {
          layer[k]=features[k];
        }
        sum[j] += net1.Weights[i][j*net1.dimensions[i]+k]*layer[k];
      }
      sum[j] += net1.Biases[i][j];
      dsum1 = activation[itype][i]->dactivation_function(sum[j]);
      sum[j] = activation[itype][i]->activation_function(sum[j]);
      if (i==L-1) {
        energy[j] = sum[j];
        if (eflag_atom) eatom[i1]=sum[j];
        if (eflag_global) eng_vdwl +=sum[j];
      }
      //force propagation
      for (jj=0;jj<jnum;jj++) {
        dlayersumx[jj][j]=0;
        dlayersumy[jj][j]=0;
        dlayersumz[jj][j]=0;
        dssumx[jj][j]=0;
        dssumy[jj][j]=0;
        dssumz[jj][j]=0;
        for (k=0;k<net1.dimensions[i];k++) {
          if (i==0&&j==0) {
            dlayerx[jj][k]=dfeaturesx[jj*f+k];
            dlayery[jj][k]=dfeaturesy[jj*f+k];
            dlayerz[jj][k]=dfeaturesz[jj*f+k];
            dsx[jj][k]=-sx[jj*f+k];
            dsy[jj][k]=-sy[jj*f+k];
            dsz[jj][k]=-sz[jj*f+k];
          }
          double w1 = net1.Weights[i][j*net1.dimensions[i]+k];
          dlayersumx[jj][j] += w1*dlayerx[jj][k];
          dlayersumy[jj][j] += w1*dlayery[jj][k];
          dlayersumz[jj][j] += w1*dlayerz[jj][k];
          dssumx[jj][j] += w1*dsx[jj][k];
          dssumy[jj][j] += w1*dsy[jj][k];
          dssumz[jj][j] += w1*dsz[jj][k];
        }
        dlayersumx[jj][j] = dsum1*dlayersumx[jj][j];
        dlayersumy[jj][j] = dsum1*dlayersumy[jj][j];
        dlayersumz[jj][j] = dsum1*dlayersumz[jj][j];
        dssumx[jj][j] *= dsum1;
        dssumy[jj][j] *= dsum1;
        dssumz[jj][j] *= dsum1;
        if (i==L-1 && jj < (jnum-1)) {
          j1 = jl[jj];
          force[j1][0]+=dlayersumx[jj][j];
          force[j1][1]+=dlayersumy[jj][j];
          force[j1][2]+=dlayersumz[jj][j];
          fm[j1][0]+=dssumx[jj][j];
          fm[j1][1]+=dssumy[jj][j];
          fm[j1][2]+=dssumz[jj][j];
        }
      }
      if (i==L-1) {
        j1 = ilist[ii];
        jj = jnum-1;
        force[j1][0]+=dlayersumx[jj][j];
        force[j1][1]+=dlayersumy[jj][j];
        force[j1][2]+=dlayersumz[jj][j];
        fm[j1][0]+=dssumx[jj][j];
        fm[j1][1]+=dssumy[jj][j];
        fm[j1][2]+=dssumz[jj][j];
      }
    }
    //update values for next iteration
    for (j=0;j<net1.dimensions[i+1];j++) {
      layer[j]=sum[j];
      for (jj=0;jj<jnum;jj++) {
        dlayerx[jj][j] = dlayersumx[jj][j];
        dlayery[jj][j] = dlayersumy[jj][j];
        dlayerz[jj][j] = dlayersumz[jj][j];
        dsx[jj][j] = dssumx[jj][j];
        dsy[jj][j] = dssumy[jj][j];
        dsz[jj][j] = dssumz[jj][j];
      }
    }
  }
}

void PairRANN::init_list(int /*which*/, NeighList *ptr)
{
  listfull = ptr;
}

void PairRANN::init_style()
{
    int irequest_full = neighbor->request(this,instance_me);
    neighbor->requests[irequest_full]->id = 1;
    neighbor->requests[irequest_full]->half = 0;
    neighbor->requests[irequest_full]->full = 1;
}


/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

double PairRANN::init_one(int /*i*/, int /*j*/)
{
  return cutmax;
}

void PairRANN::errorf(const char *file, int line, const char * message) {
  error->one(file,line,message);
}

int PairRANN::factorial(int n) {
        return round(MathSpecial::factorial(n));
}

RANN::Fingerprint *PairRANN::create_fingerprint(const char *style)
{
  if (strcmp(style,"radial")==0) {
          return new RANN::Fingerprint_radial(this);
  }
  else if (strcmp(style,"radialscreened")==0) {
          return new RANN::Fingerprint_radialscreened(this);
  }
  else if (strcmp(style,"radialscreenedspin")==0) {
          return new RANN::Fingerprint_radialscreenedspin(this);
  }
  else if (strcmp(style,"radialspin")==0) {
          return new RANN::Fingerprint_radialspin(this);
  }
  else if (strcmp(style,"bond")==0) {
          return new RANN::Fingerprint_bond(this);
  }
  else if (strcmp(style,"bondscreened")==0) {
          return new RANN::Fingerprint_bondscreened(this);
  }
  else if (strcmp(style,"bondscreenedspin")==0) {
          return new RANN::Fingerprint_bondscreenedspin(this);
  }
  else if (strcmp(style,"bondspin")==0) {
          return new RANN::Fingerprint_bondspin(this);
  }
  error->one(FLERR,"Unknown fingerprint style {}",style);
  return nullptr;
}


RANN::Activation *PairRANN::create_activation(const char *style)
{
  if (strcmp(style,"linear")==0) {
          return new RANN::Activation_linear(this);
  }
  else if (strcmp(style,"sigI")==0) {
          return new RANN::Activation_sigI(this);
  }
  error->one(FLERR,"Unknown activation style {}",style);
  return nullptr;
}