Merge pull request #3142 from charlessievers/chem_snap

SNAP bik flag

doc/src/compute_sna_atom.rst

@@ -33,7 +33,7 @@ Syntax
 * R_1, R_2,... = list of cutoff radii, one for each type (distance units)
 * w_1, w_2,... = list of neighbor weights, one for each type
 * zero or more keyword/value pairs may be appended
-* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag*
+* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag*
 
   .. parsed-literal::
 
@@ -56,6 +56,9 @@ Syntax
        *wselfallflag* value = *0* or *1*
           *0* = self-contribution only for element of central atom
           *1* = self-contribution for all elements
+       *bikflag* value = *0* or *1* (only implemented for compute snap)
+          *0* = per-atom bispectrum descriptors are summed over atoms
+          *1* = per-atom bispectrum descriptors are not summed over atoms
 
 Examples
 """"""""
@@ -296,6 +299,15 @@ This option is typically used in conjunction with the *chem* keyword,
 and LAMMPS will generate a warning if both *chem* and *bnormflag*
 are not both set or not both unset.
 
+The keyword *bikflag* determines whether or not to expand the bispectrum
+rows of the global array returned by compute snap. If *bikflag* is set
+to *1* then the bispectrum row, which is typically the per-atom bispectrum
+descriptors :math:`B_{i,k}` summed over all atoms *i* to produce
+:math:`B_k`, becomes bispectrum rows equal to the number of atoms. Thus,
+the resulting bispectrum rows are :math:`B_{i,k}` instead of just
+:math:`B_k`. In this case, the entries in the final column for these rows
+are set to zero.
+
 .. note::
 
    If you have a bonded system, then the settings of :doc:`special_bonds

src/ML-SNAP/compute_snap.cpp

@@ -57,6 +57,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   switchflag = 1;
   bzeroflag = 1;
   quadraticflag = 0;
+  bikflag = 0;
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
@@ -137,6 +138,11 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Illegal compute snap command");
       wselfallflag = atoi(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"bikflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute snap command");
+      bikflag = atoi(arg[iarg+1]);
+      iarg += 2;
     } else error->all(FLERR,"Illegal compute snap command");
   }
 
@@ -152,7 +158,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   yoffset = nperdim;
   zoffset = 2*nperdim;
   natoms = atom->natoms;
-  size_array_rows = 1+ndims_force*natoms+ndims_virial;
+  bik_rows = 1;
+  if (bikflag) bik_rows = natoms;
+  size_array_rows = bik_rows+ndims_force*natoms+ndims_virial;
   size_array_cols = nperdim*atom->ntypes+1;
   lastcol = size_array_cols-1;
 
@@ -287,6 +295,8 @@ void ComputeSnap::compute_array()
   const int* const mask = atom->mask;
 
   for (int ii = 0; ii < inum; ii++) {
+    int irow = 0;
+    if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1;
     const int i = ilist[ii];
     if (mask[i] & groupbit) {
 
@@ -417,17 +427,17 @@ void ComputeSnap::compute_array()
 
       int k = typeoffset_global;
       for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-        snap[0][k++] += snaptr->blist[icoeff];
+        snap[irow][k++] += snaptr->blist[icoeff];
 
       // quadratic contributions
 
       if (quadraticflag) {
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
           double bveci = snaptr->blist[icoeff];
-          snap[0][k++] += 0.5*bveci*bveci;
+          snap[irow][k++] += 0.5*bveci*bveci;
           for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
             double bvecj = snaptr->blist[jcoeff];
-            snap[0][k++] += bveci*bvecj;
+            snap[irow][k++] += bveci*bvecj;
           }
         }
       }
@@ -443,10 +453,10 @@ void ComputeSnap::compute_array()
       for (int i = 0; i < ntotal; i++) {
         double *snadi = snap_peratom[i]+typeoffset_local;
         int iglobal = atom->tag[i];
-        int irow = 3*(iglobal-1)+1;
-        snap[irow][icoeff+typeoffset_global] += snadi[icoeff];
-        snap[irow+1][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
-        snap[irow+2][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
+        int irow = 3*(iglobal-1)+bik_rows;
+        snap[irow++][icoeff+typeoffset_global] += snadi[icoeff];
+        snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
+        snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
       }
     }
   }
@@ -455,10 +465,10 @@ void ComputeSnap::compute_array()
 
   for (int i = 0; i < atom->nlocal; i++) {
     int iglobal = atom->tag[i];
-    int irow = 3*(iglobal-1)+1;
-    snap[irow][lastcol] = atom->f[i][0];
-    snap[irow+1][lastcol] = atom->f[i][1];
-    snap[irow+2][lastcol] = atom->f[i][2];
+    int irow = 3*(iglobal-1)+bik_rows;
+    snap[irow++][lastcol] = atom->f[i][0];
+    snap[irow++][lastcol] = atom->f[i][1];
+    snap[irow][lastcol] = atom->f[i][2];
   }
 
   // accumulate bispectrum virial contributions to global array
@@ -470,22 +480,22 @@ void ComputeSnap::compute_array()
   MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
 
   // assign energy to last column
-
+  for (int i = 0; i < bik_rows; i++) snapall[i][lastcol] = 0;
   int irow = 0;
   double reference_energy = c_pe->compute_scalar();
-  snapall[irow++][lastcol] = reference_energy;
+  snapall[irow][lastcol] = reference_energy;
 
   // assign virial stress to last column
   // switch to Voigt notation
 
   c_virial->compute_vector();
-  irow += 3*natoms;
+  irow += 3*natoms+bik_rows;
   snapall[irow++][lastcol] = c_virial->vector[0];
   snapall[irow++][lastcol] = c_virial->vector[1];
   snapall[irow++][lastcol] = c_virial->vector[2];
   snapall[irow++][lastcol] = c_virial->vector[5];
   snapall[irow++][lastcol] = c_virial->vector[4];
-  snapall[irow++][lastcol] = c_virial->vector[3];
+  snapall[irow][lastcol] = c_virial->vector[3];
 
 }
 
@@ -497,7 +507,7 @@ void ComputeSnap::compute_array()
 void ComputeSnap::dbdotr_compute()
 {
   double **x = atom->x;
-  int irow0 = 1+ndims_force*natoms;
+  int irow0 = bik_rows+ndims_force*natoms;
 
   // sum over bispectrum contributions to forces
   // on all particles including ghosts
@@ -518,7 +528,7 @@ void ComputeSnap::dbdotr_compute()
         snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][2];
         snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][1];
         snap[irow++][icoeff+typeoffset_global] += dbdz*x[i][0];
-        snap[irow++][icoeff+typeoffset_global] += dbdy*x[i][0];
+        snap[irow][icoeff+typeoffset_global] += dbdy*x[i][0];
       }
     }
 }

src/ML-SNAP/compute_snap.h

@@ -49,6 +49,8 @@ class ComputeSnap : public Compute {
   class SNA *snaptr;
   double cutmax;
   int quadraticflag;
+  int bikflag;
+  int bik_rows;
 
   Compute *c_pe;
   Compute *c_virial;