Added a lot of boilerplate for inner data, not working yet

src/ML-IAP/mliap_descriptor_snap.cpp

@@ -46,7 +46,8 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename):
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   ndescriptors = snaptr->ncoeff;
 }
@@ -162,10 +163,12 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data)
       int j = snaptr->inside[jj];
       if (chemflag)
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
       else
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, 0);
+			       snaptr->rcutij[jj], jj, 0,
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
 
       snaptr->compute_deidrj(fij);
 
@@ -236,10 +239,12 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data)
 
       if (chemflag)
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
       else
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, 0);
+			       snaptr->rcutij[jj], jj, 0,
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
 
       snaptr->compute_dbidrj();
 
@@ -309,10 +314,12 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data)
     for (int jj = 0; jj < ninside; jj++) {
       if (chemflag)
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
       else
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, 0);
+			       snaptr->rcutij[jj], jj, 0,
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
 
       snaptr->compute_dbidrj();
 
@@ -361,6 +368,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
 
   for (int i = 0; i < nelements; i++) delete[] elements[i];
   delete[] elements;

src/ML-IAP/mliap_descriptor_snap.h

@@ -39,6 +39,7 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor {
 
   int twojmax, switchflag, bzeroflag;
   int chemflag, bnormflag, wselfallflag;
+  int switchinnerflag;
   double rfac0, rmin0;
 };
 

src/ML-SNAP/compute_sna_atom.cpp

@@ -54,6 +54,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
   nelements = 1;
 
   // offset by 1 to match up with types
@@ -138,7 +139,8 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   ncoeff = snaptr->ncoeff;
   size_peratom_cols = ncoeff;

src/ML-SNAP/compute_sna_atom.h

@@ -44,6 +44,7 @@ class ComputeSNAAtom : public Compute {
   double *wjelem;
   int *map;    // map types to [0,nelements)
   int nelements, chemflag;
+  int switchinnerflag;
   class SNA *snaptr;
   double cutmax;
   int quadraticflag;

src/ML-SNAP/compute_snad_atom.cpp

@@ -54,6 +54,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
   nelements = 1;
 
   // process required arguments
@@ -136,7 +137,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   ncoeff = snaptr->ncoeff;
   nperdim = ncoeff;
@@ -300,7 +302,8 @@ void ComputeSNADAtom::compute_peratom()
       for (int jj = 0; jj < ninside; jj++) {
         const int j = snaptr->inside[jj];
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                                    snaptr->rcutij[jj], jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
         snaptr->compute_dbidrj();
 
         // Accumulate -dBi/dRi, -dBi/dRj

src/ML-SNAP/compute_snad_atom.h

@@ -46,6 +46,7 @@ class ComputeSNADAtom : public Compute {
   double *wjelem;
   int *map;    // map types to [0,nelements)
   int nelements, chemflag;
+  int switchinnerflag;
   class SNA *snaptr;
   double cutmax;
   int quadraticflag;

src/ML-SNAP/compute_snap.cpp

@@ -61,6 +61,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
   nelements = 1;
 
   // process required arguments
@@ -148,7 +149,8 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   ncoeff = snaptr->ncoeff;
   nperdim = ncoeff;
@@ -354,7 +356,8 @@ void ComputeSnap::compute_array()
       for (int jj = 0; jj < ninside; jj++) {
         const int j = snaptr->inside[jj];
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                                    snaptr->rcutij[jj], jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
         snaptr->compute_dbidrj();
 
         // Accumulate dBi/dRi, -dBi/dRj

src/ML-SNAP/compute_snap.h

@@ -46,6 +46,7 @@ class ComputeSnap : public Compute {
   double *wjelem;
   int *map;    // map types to [0,nelements)
   int nelements, chemflag;
+  int switchinnerflag;
   class SNA *snaptr;
   double cutmax;
   int quadraticflag;

src/ML-SNAP/compute_snav_atom.cpp

@@ -53,6 +53,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
   nelements = 1;
 
   // process required arguments
@@ -131,7 +132,8 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   ncoeff = snaptr->ncoeff;
   nperdim = ncoeff;
@@ -294,7 +296,8 @@ void ComputeSNAVAtom::compute_peratom()
         const int j = snaptr->inside[jj];
 
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                                    snaptr->rcutij[jj], jj, snaptr->element[jj]);
+			       snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
         snaptr->compute_dbidrj();
 
         // Accumulate -dBi/dRi*Ri, -dBi/dRj*Rj

src/ML-SNAP/compute_snav_atom.h

@@ -46,6 +46,7 @@ class ComputeSNAVAtom : public Compute {
   double *wjelem;
   int *map;    // map types to [0,nelements)
   int nelements, chemflag;
+  int switchinnerflag;
   class SNA *snaptr;
   int quadraticflag;
 };

src/ML-SNAP/pair_snap.cpp

@@ -174,10 +174,13 @@ void PairSNAP::compute(int eflag, int vflag)
       int j = snaptr->inside[jj];
       if (chemflag)
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, snaptr->element[jj]);
+                               snaptr->rcutij[jj], jj, snaptr->element[jj],
+			       snaptr->rinnerij[jj], snaptr->drinnerij[jj]);
       else
         snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
-                               snaptr->rcutij[jj],jj, 0);
+                               snaptr->rcutij[jj], jj, 0,
+			       //			       snaptr->rinnerij[jj], snaptr->drinnerij[jj],jj);
+			       0.0, 1.0);
 
       snaptr->compute_deidrj(fij);
 
@@ -403,7 +406,8 @@ void PairSNAP::coeff(int narg, char **arg)
 
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag, nelements);
+                   chemflag, bnormflag, wselfallflag,
+		   nelements, switchinnerflag);
 
   if (ncoeff != snaptr->ncoeff) {
     if (comm->me == 0)
@@ -626,6 +630,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
+  switchinnerflag = 1;
   chunksize = 32768;
   parallel_thresh = 8192;
 

src/ML-SNAP/pair_snap.h

@@ -59,6 +59,7 @@ class PairSNAP : public Pair {
   double **scale;         // for thermodynamic integration
   int twojmax, switchflag, bzeroflag, bnormflag;
   int chemflag, wselfallflag;
+  int switchinnerflag;
   int chunksize,parallel_thresh;
   double rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag;    // flags for required parameters

src/ML-SNAP/sna.cpp

@@ -103,7 +103,7 @@ using namespace MathSpecial;
    j = |j1-j2|, |j1-j2|+2,...,j1+j2-2,j1+j2
 
    [1] Albert Bartok-Partay, "Gaussian Approximation..."
-   Doctoral Thesis, Cambrindge University, (2009)
+   Doctoral Thesis, Cambridge University, (2009)
 
    [2] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,
    "Quantum Theory of Angular Momentum," World Scientific (1988)
@@ -112,14 +112,15 @@ using namespace MathSpecial;
 
 SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in,
          double rmin0_in, int switch_flag_in, int bzero_flag_in,
-         int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in) : Pointers(lmp)
+         int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in,
+	 int nelements_in, int switch_inner_flag_in) : Pointers(lmp)
 {
   wself = 1.0;
 
   rfac0 = rfac0_in;
   rmin0 = rmin0_in;
   switch_flag = switch_flag_in;
-  switch_inner_flag = 0;
+  switch_inner_flag = switch_inner_flag_in;
   bzero_flag = bzero_flag_in;
   chem_flag = chem_flag_in;
   bnorm_flag = bnorm_flag_in;
@@ -142,6 +143,8 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in,
   inside = nullptr;
   wj = nullptr;
   rcutij = nullptr;
+  rinnerij = nullptr;
+  drinnerij = nullptr;
   element = nullptr;
   nmax = 0;
   idxz = nullptr;
@@ -171,7 +174,11 @@ SNA::~SNA()
   memory->destroy(inside);
   memory->destroy(wj);
   memory->destroy(rcutij);
-  memory->destroy(element);
+  if (switch_inner_flag) {
+    memory->destroy(rinnerij);
+    memory->destroy(drinnerij);
+  }
+if (chem_flag) memory->destroy(element);
   memory->destroy(ulist_r_ij);
   memory->destroy(ulist_i_ij);
   delete[] idxz;
@@ -308,7 +315,6 @@ void SNA::init()
   init_rootpqarray();
 }
 
-
 void SNA::grow_rij(int newnmax)
 {
   if (newnmax <= nmax) return;
@@ -319,14 +325,22 @@ void SNA::grow_rij(int newnmax)
   memory->destroy(inside);
   memory->destroy(wj);
   memory->destroy(rcutij);
-  memory->destroy(element);
+  if (switch_inner_flag) {
+    memory->destroy(rinnerij);
+    memory->destroy(drinnerij);
+  }
+  if (chem_flag) memory->destroy(element);
   memory->destroy(ulist_r_ij);
   memory->destroy(ulist_i_ij);
   memory->create(rij, nmax, 3, "pair:rij");
   memory->create(inside, nmax, "pair:inside");
   memory->create(wj, nmax, "pair:wj");
   memory->create(rcutij, nmax, "pair:rcutij");
-  memory->create(element, nmax, "sna:element");
+  if (switch_inner_flag) {
+    memory->create(rinnerij, nmax, "pair:rinnerij");
+    memory->create(drinnerij, nmax, "pair:drinnerij");
+  }
+  if (chem_flag) memory->create(element, nmax, "sna:element");
   memory->create(ulist_r_ij, nmax, idxu_max, "sna:ulist_ij");
   memory->create(ulist_i_ij, nmax, idxu_max, "sna:ulist_ij");
 }
@@ -360,9 +374,9 @@ void SNA::compute_ui(int jnum, int ielem)
 
     compute_uarray(x, y, z, z0, r, j);
     if (chem_flag)
-      add_uarraytot(r, wj[j], rcutij[j], j, element[j]);
+      add_uarraytot(r, wj[j], rcutij[j], j, element[j], rinnerij[j], drinnerij[j]);
     else
-      add_uarraytot(r, wj[j], rcutij[j], j, 0);
+      add_uarraytot(r, wj[j], rcutij[j], j, 0, rinnerij[j], drinnerij[j]);
   }
 
 }
@@ -952,7 +966,8 @@ void SNA::compute_dbidrj()
    calculate derivative of Ui w.r.t. atom j
 ------------------------------------------------------------------------- */
 
-void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, int jelem)
+void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj,
+			 int jelem, double rinner, double drinner)
 {
   double rsq, r, x, y, z, z0, theta0, cs, sn;
   double dz0dr;
@@ -970,7 +985,7 @@ void SNA::compute_duidrj(double* rij, double wj, double rcut, int jj, int jelem)
   dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq;
 
   elem_duarray = jelem;
-  compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut, jj);
+  compute_duarray(x, y, z, z0, r, dz0dr, wj, rcut, jj, rinner, drinner);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -999,18 +1014,17 @@ void SNA::zero_uarraytot(int ielem)
    add Wigner U-functions for one neighbor to the total
 ------------------------------------------------------------------------- */
 
-void SNA::add_uarraytot(double r, double wj, double rcut, int jj, int jelem)
+void SNA::add_uarraytot(double r, double wj, double rcut, int jj, int jelem,
+			double rinner, double drinner)
 {
   double sfac;
-  double rinner = 0.0;
-  double drinner = 1.0;
 
   sfac = compute_sfac(r, rcut);
   sfac *= wj;
 
   if (switch_inner_flag)
     sfac *= compute_sfac_inner(r, rinner, drinner);
-
+  
   double* ulist_r = ulist_r_ij[jj];
   double* ulist_i = ulist_i_ij[jj];
 
@@ -1123,7 +1137,8 @@ void SNA::compute_uarray(double x, double y, double z,
 
 void SNA::compute_duarray(double x, double y, double z,
                           double z0, double r, double dz0dr,
-                          double wj, double rcut, int jj)
+                          double wj, double rcut, int jj,
+			  double rinner, double drinner)
 {
   double r0inv;
   double a_r, a_i, b_r, b_i;
@@ -1254,11 +1269,9 @@ void SNA::compute_duarray(double x, double y, double z,
   double dsfac = compute_dsfac(r, rcut);
 
   if (switch_inner_flag) {
-    double rinner = 0.0;
+    double sfac_inner = compute_sfac_inner(r, rinner, drinner);
-    double drinner = 1.0;
+    dsfac = dsfac*sfac_inner + sfac*compute_dsfac_inner(r, rinner, drinner);
-  
+    sfac *= sfac_inner;
-    sfac *= compute_sfac_inner(r, rinner, drinner);
-    dsfac *= compute_dsfac_inner(r, rinner, drinner);
   }
   
   sfac *= wj;
@@ -1323,7 +1336,12 @@ double SNA::memory_usage()
   bytes += (double)nmax * sizeof(int);                           // inside
   bytes += (double)nmax * sizeof(double);                        // wj
   bytes += (double)nmax * sizeof(double);                        // rcutij
-
+  if (switch_inner_flag) {
+    bytes += (double)nmax * sizeof(double);                      // rinnerij
+    bytes += (double)nmax * sizeof(double);                      // drinnerij
+  }
+  if (chem_flag) bytes += (double)nmax * sizeof(int);            // element
+  
   return bytes;
 }
 /* ---------------------------------------------------------------------- */

src/ML-SNAP/sna.h

@@ -33,7 +33,7 @@ struct SNA_BINDICES {
 class SNA : protected Pointers {
 
  public:
-  SNA(LAMMPS *, double, int, double, int, int, int, int, int, int);
+  SNA(LAMMPS *, double, int, double, int, int, int, int, int, int, int);
 
   SNA(LAMMPS *lmp) : Pointers(lmp){};
   ~SNA() override;
@@ -53,7 +53,7 @@ class SNA : protected Pointers {
 
   // functions for derivatives
 
-  void compute_duidrj(double *, double, double, int, int);
+  void compute_duidrj(double *, double, double, int, int, double, double);
   void compute_dbidrj();
   void compute_deidrj(double *);
   double compute_sfac(double, double);
@@ -64,12 +64,14 @@ class SNA : protected Pointers {
   
   double *blist;
   double **dblist;
-  double **rij;
+  double **rij;     // short rij list
-  int *inside;
+  int *inside;      // short neighbor list
-  double *wj;
+  double *wj;       // short weight list
-  double *rcutij;
+  double *rcutij;   // short cutoff list 
-  int *element;    // index on [0,nelements)
+  double *rinnerij; // short inner cutoff list
-  int nmax;
+  double *drinnerij;// short inner range list
+  int *element;     // short element list [0,nelements)
+  int nmax;         // allocated size of short lists
 
   void grow_rij(int);
 
@@ -90,7 +92,7 @@ class SNA : protected Pointers {
   int ***idxcg_block;
 
   double *ulisttot_r, *ulisttot_i;
-  double **ulist_r_ij, **ulist_i_ij;
+  double **ulist_r_ij, **ulist_i_ij; // short u list
   int *idxu_block;
 
   double *zlist_r, *zlist_i;
@@ -107,11 +109,12 @@ class SNA : protected Pointers {
   void print_clebsch_gordan();
   void init_rootpqarray();
   void zero_uarraytot(int);
-  void add_uarraytot(double, double, double, int, int);
+  void add_uarraytot(double, double, double, int, int, double, double);
   void compute_uarray(double, double, double, double, double, int);
   double deltacg(int, int, int);
   void compute_ncoeff();
-  void compute_duarray(double, double, double, double, double, double, double, double, int);
+  void compute_duarray(double, double, double, double, double, double,
+		       double, double, int, double, double);
 
   // Sets the style for the switching function
   // 0 = none