added support for USER-AWPMD package

src/USER-AWPMD/atom_vec_wavepacket.cpp

@@ -61,6 +61,20 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
   setup_fields();
 }
 
+/* ----------------------------------------------------------------------
+   set local copies of all grow ptrs used by this class, except defaults
+   needed in replicate when 2 atom classes exist and it calls pack_restart()
+------------------------------------------------------------------------- */
+
+void AtomVecWavepacket::grow_pointers()
+{
+  q = atom->q;
+  spin = atom->spin;
+  eradius = atom->eradius;
+  ervel = atom->ervel;
+  erforce = atom->erforce;
+}
+
 /* ----------------------------------------------------------------------
    clear extra forces starting at atom N
    nbytes = # of bytes to clear for a per-atom vector
@@ -68,7 +82,7 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
 
 void AtomVecWavepacket::force_clear(int n, size_t nbytes)
 {
-  memset(&atom->erforce[n],0,nbytes);
+  memset(&erforce[n],0,nbytes);
 }
 
 /* ----------------------------------------------------------------------
@@ -78,7 +92,7 @@ void AtomVecWavepacket::force_clear(int n, size_t nbytes)
 
 void AtomVecWavepacket::create_atom_post(int ilocal)
 {
-  atom->q[ilocal] = 1.0;
+  q[ilocal] = 1.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -88,7 +102,7 @@ void AtomVecWavepacket::create_atom_post(int ilocal)
 
 void AtomVecWavepacket::data_atom_post(int ilocal)
 {
-  atom->ervel[ilocal] = 0.0;
+  ervel[ilocal] = 0.0;
 }
 
 /* ----------------------------------------------------------------------
@@ -113,33 +127,28 @@ int AtomVecWavepacket::property_atom(char *name)
 void AtomVecWavepacket::pack_property_atom(int index, double *buf,
                                            int nvalues, int groupbit)
 {
-  int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
   int n = 0;
   if (index == 0) {
-    int *spin = atom->spin;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = spin[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 1) {
-    double *eradius = atom->eradius;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = eradius[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 2) {
-    double *ervel = atom->ervel;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = ervel[i];
       else buf[n] = 0.0;
       n += nvalues;
     }
   } else if (index == 3) {
-    double *erforce = atom->erforce;
     for (int i = 0; i < nlocal; i++) {
       if (mask[i] & groupbit) buf[n] = erforce[i];
       else buf[n] = 0.0;

src/USER-AWPMD/atom_vec_wavepacket.h

@@ -27,11 +27,17 @@ namespace LAMMPS_NS {
 class AtomVecWavepacket : public AtomVec {
 public:
   AtomVecWavepacket(class LAMMPS *);
+
+  void grow_pointers();
   void force_clear(int, size_t);
   void create_atom_post(int);
   void data_atom_post(int);
   int property_atom(char *);
   void pack_property_atom(int, double *, int, int);
+
+ private:
+  int *spin;
+  double *q,*eradius,*ervel,*erforce;
 };
 
 }

src/USER-AWPMD/pair_awpmd_cut.cpp

@@ -259,7 +259,7 @@ void PairAWPMDCut::compute(int eflag, int vflag)
       Vector_3 xx=Vector_3(x[i][0],x[i][1],x[i][2]);
       Vector_3 rv=m*Vector_3(v[i][0],v[i][1],v[i][2]);
       double pv=ermscale*m*atom->ervel[i];
-      Vector_2 cc=Vector_2(atom->cs[2*i],atom->cs[2*i+1]);
+      Vector_2 cc=Vector_2(atom->cs[i][0],atom->cs[i][1]);
       gmap[i]=wpmd->add_split(xx,rv,atom->eradius[i],pv,cc,1.,atom->q[i],i<nlocal ? atom->tag[i] : -atom->tag[i]);
       // resetting for the case constraints were applied
       v[i][0]=rv[0]/m;
@@ -284,7 +284,7 @@ void PairAWPMDCut::compute(int eflag, int vflag)
     } else { // electron
       int iel=gmap[i];
       int s=spin[i] >0 ? 0 : 1;
-      wpmd->get_wp_force(s,iel,(Vector_3 *)f[i],(Vector_3 *)(atom->vforce+3*i),atom->erforce+i,atom->ervelforce+i,(Vector_2 *)(atom->csforce+2*i));
+      wpmd->get_wp_force(s,iel,(Vector_3 *)f[i],(Vector_3 *)(atom->vforce[i]),atom->erforce+i,atom->ervelforce+i,(Vector_2 *)(atom->csforce[i]));
     }
   }
 
@@ -671,11 +671,11 @@ void PairAWPMDCut::min_xf_get(int /* ignore */)
   double *eradius = atom->eradius;
   double *erforce = atom->erforce;
   double **v=atom->v;
-  double *vforce=atom->vforce;
+  double **vforce=atom->vforce;
   double *ervel=atom->ervel;
   double *ervelforce=atom->ervelforce;
-  double *cs=atom->cs;
+  double **cs=atom->cs;
-  double *csforce=atom->csforce;
+  double **csforce=atom->csforce;
 
   int *spin = atom->spin;
   int nlocal = atom->nlocal;
@@ -686,14 +686,14 @@ void PairAWPMDCut::min_xf_get(int /* ignore */)
       min_varforce[7*i] = eradius[i]*erforce[i];
       for(int j=0;j<3;j++){
         min_var[7*i+1+3*j] = v[i][j];
-        min_varforce[7*i+1+3*j] = vforce[3*i+j];
+        min_varforce[7*i+1+3*j] = vforce[i][j];
       }
       min_var[7*i+4] = ervel[i];
       min_varforce[7*i+4] = ervelforce[i];
-      min_var[7*i+5] = cs[2*i];
+      min_var[7*i+5] = cs[i][0];
-      min_varforce[7*i+5] = csforce[2*i];
+      min_varforce[7*i+5] = csforce[i][0];
-      min_var[7*i+6] = cs[2*i+1];
+      min_var[7*i+6] = cs[i][1];
-      min_varforce[7*i+6] = csforce[2*i+1];
+      min_varforce[7*i+6] = csforce[i][1];
 
     } else {
       for(int j=0;j<7;j++)
@@ -710,7 +710,7 @@ void PairAWPMDCut::min_x_set(int /* ignore */)
   double *eradius = atom->eradius;
   double **v=atom->v;
   double *ervel=atom->ervel;
-  double *cs=atom->cs;
+  double **cs=atom->cs;
 
   int *spin = atom->spin;
   int nlocal = atom->nlocal;
@@ -721,8 +721,8 @@ void PairAWPMDCut::min_x_set(int /* ignore */)
       for(int j=0;j<3;j++)
         v[i][j]=min_var[7*i+1+3*j];
       ervel[i]=min_var[7*i+4];
-      cs[2*i]=min_var[7*i+5];
+      cs[i][0]=min_var[7*i+5];
-      cs[2*i+1]=min_var[7*i+6];
+      cs[i][1]=min_var[7*i+6];
     }
   }
 }

src/USER-SMD/atom_vec_smd.cpp

@@ -90,7 +90,7 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp)
     "x0 vest vfrac rmass radius contact_radius molecule e "
     "eff_plastic_strain eff_plastic_strain_rate smd_data_9 smd_stress damage";
   fields_data_atom = (char *) 
-    "id type molecule vfrac rmass radius contact_radius x";
+    "id type molecule vfrac rmass radius contact_radius x0 x";
   fields_data_vel = (char *) "id v";
 
   // set these array sizes based on defines

src/atom.cpp

@@ -134,8 +134,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 
   spin = NULL;
   eradius = ervel = erforce = NULL;
-  ervelforce = cs = csforce = NULL;
+  ervelforce = NULL;
-  vforce = NULL;
+  cs = csforce = vforce = NULL;
   etag = NULL;
 
   // USER-DPD package
@@ -512,8 +512,8 @@ void Atom::peratom_create()
 
   // USER-AWPMD package
 
-  add_peratom("cs",&cs,DOUBLE,0);
+  add_peratom("cs",&cs,DOUBLE,2);
-  add_peratom("csforce",&csforce,DOUBLE,0);
+  add_peratom("csforce",&csforce,DOUBLE,2);
   add_peratom("vforce",&vforce,DOUBLE,3);
   add_peratom("ervelforce",&ervelforce,DOUBLE,0);
   add_peratom("etag",&etag,INT,0);

src/atom.h

@@ -109,8 +109,8 @@ class Atom : protected Pointers {
 
   int *spin;
   double *eradius,*ervel,*erforce;
-  double *ervelforce,*cs,*csforce;
+  double *ervelforce;
-  double **vforce;
+  double **cs,**csforce,**vforce;
   int *etag;
 
   // USER-DPD package
@@ -322,7 +322,7 @@ class Atom : protected Pointers {
   inline int get_map_size() {return map_tag_max+1;};
   inline int get_map_maxarray() {return map_maxarray+1;};
 
-  // NOTE: placeholder method until AtomVec is refactored
+  // NOTE: placeholder method until KOKKOS/AtomVec is refactored
   int memcheck(const char *) {return 1;}
 
   bigint memory_usage();

src/atom_vec.cpp

@@ -1368,7 +1368,7 @@ int AtomVec::size_restart()
   int i,nn,cols,collength,ncols;
   void *plength;
 
-  // NOTE: need to worry about overflow of returned int ??
+  // NOTE: need to worry about overflow of returned int N
 
   int nlocal = atom->nlocal;
 
@@ -1738,9 +1738,9 @@ void AtomVec::data_atom(double *coord, imageint imagetmp, char **values)
 
   // error checks applicable to all styles
 
-  if (atom->tag[nlocal] <= 0)
+  if (tag[nlocal] <= 0)
     error->one(FLERR,"Invalid atom ID in Atoms section of data file");
-  if (atom->type[nlocal] <= 0 || atom->type[nlocal] > atom->ntypes)
+  if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
 
   // if needed, modify unpacked values or initialize other peratom values
@@ -2388,7 +2388,6 @@ void AtomVec::setup_fields()
   }
 
   // set style-specific sizes
-  // NOTE: check for others vars in atom_vec.cpp/h ??
 
   comm_x_only = 1;
   if (ncomm) comm_x_only = 0;

src/atom_vec_body.cpp

@@ -132,15 +132,12 @@ void AtomVecBody::process_args(int narg, char **arg)
 }
 
 /* ----------------------------------------------------------------------
-   grow atom arrays
+   set local copies of all grow ptrs used by this class, except defaults
-   must set local copy of body ptr
+   needed in replicate when 2 atom classes exist and it calls pack_restart()
-   needed in replicate when 2 atom classes exist and pack_restart() is called
 ------------------------------------------------------------------------- */
 
-void AtomVecBody::grow(int n)
+void AtomVecBody::grow_pointers()
 {
-  AtomVec::grow(n);
-
   body = atom->body;
   rmass = atom->rmass;
   radius = atom->radius;

src/atom_vec_body.h

@@ -43,7 +43,7 @@ class AtomVecBody : public AtomVec {
   ~AtomVecBody();
   void process_args(int, char **);
 
-  void grow(int);
+  void grow_pointers();
   void copy_bonus(int, int, int);
   void clear_bonus();
   int pack_comm_bonus(int, int *, double *);

src/atom_vec_hybrid.cpp

@@ -31,12 +31,10 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp) : AtomVec(lmp)
   keywords = NULL;
   fieldstrings = NULL;
 
+  bonus_flag = 0;
   nstyles_bonus = 0;
   styles_bonus = NULL;
 
-  // NOTE: set bonus_flag if any substyle does
-  //       set nstyles_bonus, styles_bonus
-
   // these strings will be concatenated from sub-style strings
   // fields_data_atom & fields_data_vel start with fields common to all styles
 
@@ -45,6 +43,8 @@ AtomVecHybrid::AtomVecHybrid(LAMMPS *lmp) : AtomVec(lmp)
   fields_exchange = fields_restart = fields_create = (char *) "";
   fields_data_atom = (char *) "id type x";
   fields_data_vel = (char *) "id v";
+
+  fields_allocated = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -56,7 +56,10 @@ AtomVecHybrid::~AtomVecHybrid()
   for (int k = 0; k < nstyles; k++) delete [] keywords[k];
   delete [] keywords;
 
-  // NOTE: need to check these have actually been allocated
+  for (int k = 0; k < nstyles_bonus; k++) delete styles_bonus[k];
+  delete [] styles_bonus;
+
+  if (!fields_allocated) return;
 
   delete [] fields_grow;
   delete [] fields_copy;
@@ -198,6 +201,8 @@ void AtomVecHybrid::process_args(int narg, char **arg)
   fields_data_atom = merge_fields(10,fields_data_atom,0,null);
   fields_data_vel = merge_fields(11,fields_data_vel,0,null);
 
+  fields_allocated = 1;
+
   // check concat_grow for multiple special-case fields
   // may cause issues with style-specific create_atom() and data_atom() methods
   // issue warnings if appear in multiple sub-styles
@@ -219,6 +224,23 @@ void AtomVecHybrid::process_args(int narg, char **arg)
 
   delete [] concat_grow;
 
+  // set bonus_flag if any substyle has bonus data
+  // set nstyles_bonus & styles_bonus
+
+  nstyles_bonus = 0;
+  for (int k = 0; k < nstyles; k++)
+    if (styles[k]->bonus_flag) nstyles_bonus++;
+
+  if (nstyles_bonus) {
+    bonus_flag = 1;
+    styles_bonus = new AtomVec*[nstyles_bonus];
+    nstyles_bonus = 0;
+    for (int k = 0; k < nstyles; k++) {
+      if (styles[k]->bonus_flag)
+        styles_bonus[nstyles_bonus++] = styles[k];
+    }
+  }
+
   // parent AtomVec can now operate on merged fields
 
   setup_fields();
@@ -326,7 +348,13 @@ void AtomVecHybrid::copy_bonus(int i, int j, int delflag)
     styles_bonus[k]->copy_bonus(i,j,delflag);
 }
 
-// NOTE: need a clear_bonus() ?
+/* ---------------------------------------------------------------------- */
+
+void AtomVecHybrid::clear_bonus()
+{
+  for (int k = 0; k < nstyles_bonus; k++)
+    styles_bonus[k]->clear_bonus();
+}
 
 /* ---------------------------------------------------------------------- */
 

src/atom_vec_hybrid.h

@@ -38,7 +38,7 @@ class AtomVecHybrid : public AtomVec {
   void grow_pointers();
   void force_clear(int, size_t);
   void copy_bonus(int, int, int);
-  void clear_bonus() {}
+  void clear_bonus();
   int pack_comm_bonus(int, int *, double *);
   void unpack_comm_bonus(int, int, double *);
   int pack_border_bonus(int, int *, double *);
@@ -64,6 +64,7 @@ class AtomVecHybrid : public AtomVec {
  private:
   int nallstyles;
   char **allstyles;
+  int fields_allocated;
 
   struct FieldStrings {
     char **fstr;