diff --git a/src/MISC/fix_efield.h b/src/MISC/fix_efield.h
index 4b0e9d63bc..1a79d34e8f 100644
--- a/src/MISC/fix_efield.h
+++ b/src/MISC/fix_efield.h
@@ -47,7 +47,6 @@ class FixEfield : public Fix {
   int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
   int nlevels_respa;
   double qe2f;
-  double fdotx;
   int qflag,muflag;
 
   int maxatom;
diff --git a/src/accelerator_kokkos.h b/src/accelerator_kokkos.h
index 5d7124e15a..2b2ba1a957 100644
--- a/src/accelerator_kokkos.h
+++ b/src/accelerator_kokkos.h
@@ -20,16 +20,15 @@
 #ifdef LMP_KOKKOS
 
 #include "kokkos.h"
+#include "atom_kokkos.h"
+#include "comm_kokkos.h"
+#include "modify_kokkos.h"
 
 #else
 
 // dummy interface to KOKKOS
 // needed for compiling when KOKKOS is not installed
 
-//#include "comm.h"
-//#include "modify.h"
-//#include "verlet.h"
-
 namespace LAMMPS_NS {
 
 class KokkosLMP {
@@ -41,18 +40,18 @@ class KokkosLMP {
   void accelerator(int, char **) {}
 };
 
+class AtomKokkos : public Atom {
+ public:
+ AtomKokkos(class LAMMPS *lmp) : Atom(lmp) {}
+  ~AtomKokkos() {}
+};
+
 class CommKokkos : public Comm {
  public:
  CommKokkos(class LAMMPS *lmp) : Comm(lmp) {}
   ~CommKokkos() {}
 };
 
-class DomainKokkos : public Domain {
- public:
- DomainKokkos(class LAMMPS *lmp) : Domain(lmp) {}
-  ~DomainKokkos() {}
-};
-
 class NeighborKokkos : public Neighbor {
  public:
  NeighborKokkos(class LAMMPS *lmp) : Neighbor(lmp) {}
@@ -65,12 +64,6 @@ class ModifyKokkos : public Modify {
   ~ModifyKokkos() {}
 };
 
-class VerletKokkos : public Verlet {
- public:
- VerletKokkos(class LAMMPS *lmp, int narg, char **arg) : Verlet(lmp,narg,arg) {}
-  ~VerletKokkos() {}
-};
-
 }
 
 #endif
diff --git a/src/atom.h b/src/atom.h
index fae87653cd..0056ecee55 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -237,7 +237,7 @@ class Atom : protected Pointers {
   void map_delete();
   int map_find_hash(tagint);
 
- private:
+ protected:
 
   // global to local ID mapping
 
diff --git a/src/atom_masks.h b/src/atom_masks.h
index 21c5e0bb97..82e460c301 100644
--- a/src/atom_masks.h
+++ b/src/atom_masks.h
@@ -18,7 +18,7 @@
 
 #define EMPTY_MASK     0x00000000
 #define ALL_MASK       0xffffffff
-#define SAMETAG_MASK   0x01000000
+#define SAMETAG_MASK   0x40000000
 #define EXTENDED_MASK  0x80000000
 
 // standard
@@ -39,49 +39,51 @@
 #define IMPROPER_MASK  0x00002000
 #define SPECIAL_MASK   0x00004000
 #define MAP_MASK       0x00008000
+#define ENERGY_MASK    0x00010000
+#define VIRIAL_MASK    0x00020000
 
 // granular
 
-#define RADIUS_MASK    0x00010000
-#define DENSITY_MASK   0x00020000
-#define OMEGA_MASK     0x00040000
-#define TORQUE_MASK    0x00080000
-#define ANGMOM_MASK    0x00100000
-#define GRANULAR_MASK  0x001f0000
+#define RADIUS_MASK    0x00100000
+#define DENSITY_MASK   0x00200000
+#define OMEGA_MASK     0x00400000
+#define TORQUE_MASK    0x00800000
+#define ANGMOM_MASK    0x01000000
+#define GRANULAR_MASK  0x01f00000
 
 // peridynamics
 
-#define VFRAC_MASK     0x00200000
-#define S0_MASK        0x00400000
-#define X0_MASK        0x00800000
-#define PERI_MASK      0x00e00000
+#define VFRAC_MASK     0x00000001
+#define S0_MASK        0x00000002
+#define X0_MASK        0x00000004
+#define PERI_MASK      0x00000007
 
-#define ELLIPSOID_MASK 0x00000001
-#define LINE_MASK      0x00000002
-#define TRI_MASK       0x00000004
+#define ELLIPSOID_MASK 0x00000008
+#define LINE_MASK      0x00000010
+#define TRI_MASK       0x00000020
 
 // electron
 
-#define SPIN_MASK      0x00000010
-#define ERADIUS_MASK   0x00000020
-#define ERVEL_MASK     0x00000040
-#define ERFORCE_MASK   0x00000080
-#define ERVELFORCE_MASK 0x00000100
+#define SPIN_MASK      0x00000100
+#define ERADIUS_MASK   0x00000200
+#define ERVEL_MASK     0x00000400
+#define ERFORCE_MASK   0x00000800
+#define ERVELFORCE_MASK 0x00001000
 
-#define CS_MASK        0x00000200
-#define CSFORCE_MASK   0x00000400
-#define VFORCE_MASK    0x00000800
+#define CS_MASK        0x00002000
+#define CSFORCE_MASK   0x00004000
+#define VFORCE_MASK    0x00008000
 
-#define ELECTRON_MASK  0x00000f00
+#define ELECTRON_MASK  0x0000ff00
 
 // SPH
 
-#define ETAG_MASK      0x00001000
-#define RHO_MASK       0x00002000
-#define DRHO_MASK      0x00004000
-#define E_MASK         0x00008000
-#define DE_MASK        0x00010000
-#define VEST_MASK      0x00020000
-#define CV_MASK        0x00040000
+#define ETAG_MASK      0x00010000
+#define RHO_MASK       0x00020000
+#define DRHO_MASK      0x00040000
+#define E_MASK         0x00080000
+#define DE_MASK        0x00100000
+#define VEST_MASK      0x00200000
+#define CV_MASK        0x00400000
 
 #endif
diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp
index b381845e46..24e9e2d713 100644
--- a/src/atom_vec.cpp
+++ b/src/atom_vec.cpp
@@ -29,7 +29,7 @@ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
   bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 0;
   mass_type = dipole_type = 0;
   size_data_bonus = 0;
-  cudable = false;
+  cudable = kokkosable = 0;
 
   nargcopy = 0;
   argcopy = NULL;
diff --git a/src/atom_vec.h b/src/atom_vec.h
index 224feb1438..241e84fadb 100644
--- a/src/atom_vec.h
+++ b/src/atom_vec.h
@@ -43,6 +43,7 @@ class AtomVec : protected Pointers {
   int nset;                            // # of molecules in list
 
   int cudable;                         // 1 if atom style is CUDA-enabled
+  int kokkosable;                      // 1 if atom style is KOKKOS-enabled
   int *maxsend;                        // CUDA-specific variable
 
   int nargcopy;          // copy of command-line args for atom_style command
diff --git a/src/fix.cpp b/src/fix.cpp
index 982a6e0472..7bfd612ebd 100644
--- a/src/fix.cpp
+++ b/src/fix.cpp
@@ -75,8 +75,13 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   maxvatom = 0;
   vatom = NULL;
 
+  // CUDA and KOKKOS per-fix data masks
+
   datamask = ALL_MASK;
   datamask_ext = ALL_MASK;
+
+  datamask_read = datamask_read_ext = ALL_MASK;
+  datamask_modify = datamask_modify_ext = ALL_MASK;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/fix.h b/src/fix.h
index 6f11272be6..bf2aa7e46a 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -83,6 +83,15 @@ class Fix : protected Pointers {
   double **vatom;                // accumulated per-atom virial
 
   int restart_reset;             // 1 if restart just re-initialized fix
+
+  // KOKKOS host/device flag and per-fix data masks
+
+  ExecutionSpace execution_space;
+  unsigned int datamask_read, datamask_read_ext;
+  unsigned int datamask_modify, datamask_modify_ext;
+
+  // USER-CUDA per-fix data masks
+
   unsigned int datamask;
   unsigned int datamask_ext;
 
diff --git a/src/lammps.cpp b/src/lammps.cpp
index 1aaa0365f8..07d983d7a4 100644
--- a/src/lammps.cpp
+++ b/src/lammps.cpp
@@ -425,6 +425,11 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
     error->all(FLERR,"Small to big integers are not sized correctly");
 #endif
 
+  // error check on accelerator packages
+
+  if (cudaflag == 1 && kokkosflag == 1) 
+    error->all(FLERR,"Cannot use -cuda on and -kokkos on");
+
   // create Cuda class if USER-CUDA installed, unless explicitly switched off
   // instantiation creates dummy Cuda class if USER-CUDA is not installed
 
@@ -558,6 +563,7 @@ void LAMMPS::create()
   // so that nthreads is defined when create_avec invokes grow()
 
   if (cuda) comm = new CommCuda(this);
+  else if (kokkos) comm = new CommKokkos(this);
   else comm = new Comm(this);
 
   if (cuda) neighbor = new NeighborCuda(this);
@@ -570,13 +576,15 @@ void LAMMPS::create()
   else domain = new Domain(this);
 #endif
 
-  atom = new Atom(this);
+  if (kokkos) atom = new AtomKokkos(this);
+  else atom = new Atom(this);
   atom->create_avec("atomic",0,NULL,suffix);
 
   group = new Group(this);
   force = new Force(this);    // must be after group, to create temperature
 
   if (cuda) modify = new ModifyCuda(this);
+  else if (kokkos) modify = new ModifyKokkos(this);
   else modify = new Modify(this);
 
   output = new Output(this);  // must be after group, so "all" exists
diff --git a/src/lmptype.h b/src/lmptype.h
index 29301058b8..1390d32477 100644
--- a/src/lmptype.h
+++ b/src/lmptype.h
@@ -48,6 +48,10 @@
 
 namespace LAMMPS_NS {
 
+// enum used for KOKKOS host/device flag
+
+enum ExecutionSpace{Host,Device};
+
 // reserve 2 hi bits in molecular system neigh list for special bonds flag
 // max local + ghost atoms per processor = 2^30 - 1
 
diff --git a/src/memory.h b/src/memory.h
index e2d1eecc47..b14d0590c2 100644
--- a/src/memory.h
+++ b/src/memory.h
@@ -16,6 +16,9 @@
 
 #include "lmptype.h"
 #include "pointers.h"
+#ifdef LMP_KOKKOS
+#include "kokkos_type.h"
+#endif
 
 namespace LAMMPS_NS {
 
@@ -28,10 +31,17 @@ class Memory : protected Pointers {
   void sfree(void *);
   void fail(const char *);
 
+  // Kokkos memory allocation functions
+
+#ifdef LMP_KOKKOS
+#include "memory_kokkos.h"
+#endif
+
 /* ----------------------------------------------------------------------
    create/grow/destroy vecs and multidim arrays with contiguous memory blocks
    only use with primitive data types, e.g. 1d vec of ints, 2d array of doubles
-   cannot use with pointers, e.g. 1d vec of int*, due to mismatched destroy
+   fail() prevents use with pointers, 
+     e.g. 1d vec of int*, due to mismatched destroy
    avoid use with non-primitive data types to avoid code bloat
    for these other cases, use smalloc/srealloc/sfree directly
 ------------------------------------------------------------------------- */
@@ -41,42 +51,42 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *create(TYPE *&array, int n, const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n;
-      array = (TYPE *) smalloc(nbytes,name);
-      return array;
-    }
+  TYPE *create(TYPE *&array, int n, const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n;
+    array = (TYPE *) smalloc(nbytes,name);
+    return array;
+  }
 
   template <typename TYPE>
-    TYPE **create(TYPE **&array, int n, const char *name) {fail(name);}
+  TYPE **create(TYPE **&array, int n, const char *name) 
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *grow(TYPE *&array, int n, const char *name)
-    {
-      if (array == NULL) return create(array,n,name);
+  TYPE *grow(TYPE *&array, int n, const char *name)
+  {
+    if (array == NULL) return create(array,n,name);
 
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n;
-      array = (TYPE *) srealloc(array,nbytes,name);
-      return array;
-    }
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n;
+    array = (TYPE *) srealloc(array,nbytes,name);
+    return array;
+  }
 
   template <typename TYPE>
-    TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);}
+  TYPE **grow(TYPE **&array, int n, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    destroy a 1d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE *array)
-    {
-      sfree(array);
-    }
+  void destroy(TYPE *array)
+  {sfree(array);}
 
 /* ----------------------------------------------------------------------
    create a 1d array with index from nlo to nhi inclusive
@@ -84,51 +94,51 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * (nhi-nlo+1);
-      array = (TYPE *) smalloc(nbytes,name);
-      array -= nlo;
-      return array;
-    }
+  TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * (nhi-nlo+1);
+    array = (TYPE *) smalloc(nbytes,name);
+    array -= nlo;
+    return array;
+  }
 
   template <typename TYPE>
-    TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name)
-    {fail(name);}
+  TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    destroy a 1d array with index offset
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy1d_offset(TYPE *array, int offset)
-    {
-      if (array) sfree(&array[offset]);
-    }
+  void destroy1d_offset(TYPE *array, int offset)
+  {
+    if (array) sfree(&array[offset]);
+  }
 
 /* ----------------------------------------------------------------------
    create a 2d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **create(TYPE **&array, int n1, int n2, const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2;
-      TYPE *data = (TYPE *) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1;
-      array = (TYPE **) smalloc(nbytes,name);
-
-      bigint n = 0;
-      for (int i = 0; i < n1; i++) {
-        array[i] = &data[n];
-        n += n2;
-      }
-      return array;
+  TYPE **create(TYPE **&array, int n1, int n2, const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2;
+    TYPE *data = (TYPE *) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1;
+    array = (TYPE **) smalloc(nbytes,name);
+    
+    bigint n = 0;
+    for (int i = 0; i < n1; i++) {
+      array[i] = &data[n];
+      n += n2;
     }
+    return array;
+  }
 
   template <typename TYPE>
-    TYPE ***create(TYPE ***&array, int n1, int n2, const char *name)
-    {fail(name);}
+  TYPE ***create(TYPE ***&array, int n1, int n2, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1st dim of a 2d array
@@ -136,128 +146,128 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **grow(TYPE **&array, int n1, int n2, const char *name)
-    {
-      if (array == NULL) return create(array,n1,n2,name);
-
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2;
-      TYPE *data = (TYPE *) srealloc(array[0],nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1;
-      array = (TYPE **) srealloc(array,nbytes,name);
-
-      bigint n = 0;
-      for (int i = 0; i < n1; i++) {
-        array[i] = &data[n];
-        n += n2;
-      }
-      return array;
+  TYPE **grow(TYPE **&array, int n1, int n2, const char *name)
+  {
+    if (array == NULL) return create(array,n1,n2,name);
+    
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2;
+    TYPE *data = (TYPE *) srealloc(array[0],nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1;
+    array = (TYPE **) srealloc(array,nbytes,name);
+    
+    bigint n = 0;
+    for (int i = 0; i < n1; i++) {
+      array[i] = &data[n];
+      n += n2;
     }
-
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name)
-    {fail(name);}
-
+  TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name)
+  {fail(name); return NULL;}
+  
 /* ----------------------------------------------------------------------
    destroy a 2d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE **array)
-    {
-      if (array == NULL) return;
-      sfree(array[0]);
-      sfree(array);
-    }
+  void destroy(TYPE **array)
+  {
+    if (array == NULL) return;
+    sfree(array[0]);
+    sfree(array);
+  }
 
 /* ----------------------------------------------------------------------
    create a 2d array with a ragged 2nd dimension
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **create_ragged(TYPE **&array, int n1, int *n2, const char *name)
-    {
-      bigint n2sum = 0;
-      for (int i = 0; i < n1; i++) n2sum += n2[i];
-
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n2sum;
-      TYPE *data = (TYPE *) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1;
-      array = (TYPE **) smalloc(nbytes,name);
-
-      bigint n = 0;
-      for (int i = 0; i < n1; i++) {
-        array[i] = &data[n];
-        n += n2[i];
-      }
-      return array;
+  TYPE **create_ragged(TYPE **&array, int n1, int *n2, const char *name)
+  {
+    bigint n2sum = 0;
+    for (int i = 0; i < n1; i++) n2sum += n2[i];
+    
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n2sum;
+    TYPE *data = (TYPE *) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1;
+    array = (TYPE **) smalloc(nbytes,name);
+    
+    bigint n = 0;
+    for (int i = 0; i < n1; i++) {
+      array[i] = &data[n];
+      n += n2[i];
     }
-
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ***create_ragged(TYPE ***&array, int n1, int *n2, const char *name)
-    {fail(name);}
-
+  TYPE ***create_ragged(TYPE ***&array, int n1, int *n2, const char *name)
+  {fail(name); return NULL;}
+  
 /* ----------------------------------------------------------------------
    create a 2d array with 2nd index from n2lo to n2hi inclusive
    cannot grow it
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi,
-                           const char *name)
-    {
-      int n2 = n2hi - n2lo + 1;
-      create(array,n1,n2,name);
-      for (int i = 0; i < n1; i++) array[i] -= n2lo;
-      return array;
-    }
-
+  TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi,
+                         const char *name)
+  {
+    int n2 = n2hi - n2lo + 1;
+    create(array,n1,n2,name);
+    for (int i = 0; i < n1; i++) array[i] -= n2lo;
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi,
-                            const char *name) {fail(name);}
-
+  TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi,
+                          const char *name) {fail(name); return NULL;}
+  
 /* ----------------------------------------------------------------------
    destroy a 2d array with 2nd index offset
 ------------------------------------------------------------------------- */
-
+  
   template <typename TYPE>
-    void destroy2d_offset(TYPE **array, int offset)
-    {
-      if (array == NULL) return;
-      sfree(&array[0][offset]);
-      sfree(array);
-    }
-
+  void destroy2d_offset(TYPE **array, int offset)
+  {
+    if (array == NULL) return;
+    sfree(&array[0][offset]);
+    sfree(array);
+  }
+  
 /* ----------------------------------------------------------------------
    create a 3d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
-      TYPE *data = (TYPE *) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1*n2;
-      TYPE **plane = (TYPE **) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE **)) * n1;
-      array = (TYPE ***) smalloc(nbytes,name);
-
-      int i,j;
-      bigint m;
-      bigint n = 0;
-      for (i = 0; i < n1; i++) {
-        m = ((bigint) i) * n2;
-        array[i] = &plane[m];
-        for (j = 0; j < n2; j++) {
-          plane[m+j] = &data[n];
-          n += n3;
-        }
+  TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
+    TYPE *data = (TYPE *) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1*n2;
+    TYPE **plane = (TYPE **) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE **)) * n1;
+    array = (TYPE ***) smalloc(nbytes,name);
+    
+    int i,j;
+    bigint m;
+    bigint n = 0;
+    for (i = 0; i < n1; i++) {
+      m = ((bigint) i) * n2;
+      array[i] = &plane[m];
+      for (j = 0; j < n2; j++) {
+        plane[m+j] = &data[n];
+        n += n3;
       }
-      return array;
     }
+    return array;
+  }
 
   template <typename TYPE>
-    TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name)
-    {fail(name);}
+  TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1st dim of a 3d array
@@ -265,47 +275,47 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name)
-    {
-      if (array == NULL) return create(array,n1,n2,n3,name);
-
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
-      TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1*n2;
-      TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE **)) * n1;
-      array = (TYPE ***) srealloc(array,nbytes,name);
-
-      int i,j;
-      bigint m;
-      bigint n = 0;
-      for (i = 0; i < n1; i++) {
-        m = ((bigint) i) * n2;
-        array[i] = &plane[m];
-        for (j = 0; j < n2; j++) {
-          plane[m+j] = &data[n];
-          n += n3;
-        }
+  TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name)
+  {
+    if (array == NULL) return create(array,n1,n2,n3,name);
+    
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
+    TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1*n2;
+    TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE **)) * n1;
+    array = (TYPE ***) srealloc(array,nbytes,name);
+    
+    int i,j;
+    bigint m;
+    bigint n = 0;
+    for (i = 0; i < n1; i++) {
+      m = ((bigint) i) * n2;
+      array[i] = &plane[m];
+      for (j = 0; j < n2; j++) {
+        plane[m+j] = &data[n];
+        n += n3;
       }
-      return array;
     }
-
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name)
-    {fail(name);}
+  TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    destroy a 3d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE ***array)
-    {
-      if (array == NULL) return;
-      sfree(array[0][0]);
-      sfree(array[0]);
-      sfree(array);
-    }
+  void destroy(TYPE ***array)
+  {
+    if (array == NULL) return;
+    sfree(array[0][0]);
+    sfree(array[0]);
+    sfree(array);
+  }
 
 /* ----------------------------------------------------------------------
    create a 3d array with 1st index from n1lo to n1hi inclusive
@@ -313,29 +323,29 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi,
-                            int n2, int n3, const char *name)
-    {
-      int n1 = n1hi - n1lo + 1;
-      create(array,n1,n2,n3,name);
-      array -= n1lo;
-      return array;
-    }
-
+  TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi,
+                          int n2, int n3, const char *name)
+  {
+    int n1 = n1hi - n1lo + 1;
+    create(array,n1,n2,n3,name);
+    array -= n1lo;
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi,
-                             int n2, int n3, const char *name)
-    {fail(name);}
+  TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi,
+                           int n2, int n3, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    free a 3d array with 1st index offset
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy3d_offset(TYPE ***array, int offset)
-    {
-      if (array) destroy(&array[offset]);
-    }
+  void destroy3d_offset(TYPE ***array, int offset)
+  {
+    if (array) destroy(&array[offset]);
+  }
 
 /* ----------------------------------------------------------------------
    create a 3d array with
@@ -346,97 +356,97 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi,
-                            int n2lo, int n2hi, int n3lo, int n3hi,
-                            const char *name)
-    {
-      int n1 = n1hi - n1lo + 1;
-      int n2 = n2hi - n2lo + 1;
-      int n3 = n3hi - n3lo + 1;
-      create(array,n1,n2,n3,name);
-
-      bigint m = ((bigint) n1) * n2;
-      for (bigint i = 0; i < m; i++) array[0][i] -= n3lo;
-      for (int i = 0; i < n1; i++) array[i] -= n2lo;
-      array -= n1lo;
-      return array;
-    }
-
+  TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi,
+                          int n2lo, int n2hi, int n3lo, int n3hi,
+                          const char *name)
+  {
+    int n1 = n1hi - n1lo + 1;
+    int n2 = n2hi - n2lo + 1;
+    int n3 = n3hi - n3lo + 1;
+    create(array,n1,n2,n3,name);
+    
+    bigint m = ((bigint) n1) * n2;
+    for (bigint i = 0; i < m; i++) array[0][i] -= n3lo;
+    for (int i = 0; i < n1; i++) array[i] -= n2lo;
+    array -= n1lo;
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi,
-                             int n2lo, int n2hi, int n3lo, int n3hi,
-                             const char *name)
-    {fail(name);}
+  TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi,
+                           int n2lo, int n2hi, int n3lo, int n3hi,
+                           const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    free a 3d array with all 3 indices offset
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy3d_offset(TYPE ***array,
-                          int n1_offset, int n2_offset, int n3_offset)
-    {
-      if (array == NULL) return;
-      sfree(&array[n1_offset][n2_offset][n3_offset]);
-      sfree(&array[n1_offset][n2_offset]);
-      sfree(&array[n1_offset]);
-    }
-
+  void destroy3d_offset(TYPE ***array,
+                        int n1_offset, int n2_offset, int n3_offset)
+  {
+    if (array == NULL) return;
+    sfree(&array[n1_offset][n2_offset][n3_offset]);
+    sfree(&array[n1_offset][n2_offset]);
+    sfree(&array[n1_offset]);
+  }
+  
 /* ----------------------------------------------------------------------
    create a 4d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4,
-                    const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4;
-      TYPE *data = (TYPE *) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1*n2*n3;
-      TYPE **cube = (TYPE **) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE **)) * n1*n2;
-      TYPE ***plane = (TYPE ***) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE ***)) * n1;
-      array = (TYPE ****) smalloc(nbytes,name);
+  TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4,
+                  const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4;
+    TYPE *data = (TYPE *) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1*n2*n3;
+    TYPE **cube = (TYPE **) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE **)) * n1*n2;
+    TYPE ***plane = (TYPE ***) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE ***)) * n1;
+    array = (TYPE ****) smalloc(nbytes,name);
 
-      int i,j,k;
-      bigint m1,m2;
-      bigint n = 0;
-      for (i = 0; i < n1; i++) {
-        m2 = ((bigint) i) * n2;
-        array[i] = &plane[m2];
-        for (j = 0; j < n2; j++) {
-          m1 = ((bigint) i) * n2 + j;
-          m2 = ((bigint) i) * n2*n3 + j*n3;
-          plane[m1] = &cube[m2];
-          for (k = 0; k < n3; k++) {
-            m1 = ((bigint) i) * n2*n3 + j*n3 + k;
-            cube[m1] = &data[n];
-            n += n4;
-          }
+    int i,j,k;
+    bigint m1,m2;
+    bigint n = 0;
+    for (i = 0; i < n1; i++) {
+      m2 = ((bigint) i) * n2;
+      array[i] = &plane[m2];
+      for (j = 0; j < n2; j++) {
+        m1 = ((bigint) i) * n2 + j;
+        m2 = ((bigint) i) * n2*n3 + j*n3;
+        plane[m1] = &cube[m2];
+        for (k = 0; k < n3; k++) {
+          m1 = ((bigint) i) * n2*n3 + j*n3 + k;
+          cube[m1] = &data[n];
+          n += n4;
         }
       }
-      return array;
     }
-
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
-                     const char *name)
-    {fail(name);}
+  TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
+                   const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    destroy a 4d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE ****array)
-    {
-      if (array == NULL) return;
-      sfree(array[0][0][0]);
-      sfree(array[0][0]);
-      sfree(array[0]);
-      sfree(array);
-    }
+  void destroy(TYPE ****array)
+  {
+    if (array == NULL) return;
+    sfree(array[0][0][0]);
+    sfree(array[0][0]);
+    sfree(array[0]);
+    sfree(array);
+  }
 
 /* ----------------------------------------------------------------------
    create a 4d array with indices
@@ -447,145 +457,147 @@ class Memory : protected Pointers {
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ****create4d_offset(TYPE ****&array, int n1, int n2lo, int n2hi,
-                             int n3lo, int n3hi, int n4lo, int n4hi,
-                             const char *name)
-    {
-      int n2 = n2hi - n2lo + 1;
-      int n3 = n3hi - n3lo + 1;
-      int n4 = n4hi - n4lo + 1;
-      create(array,n1,n2,n3,n4,name);
-      
-      bigint m = ((bigint) n1) * n2 * n3;
-      for (bigint i = 0; i < m; i++) array[0][0][i] -= n4lo;
-      m = ((bigint) n1) * n2;
-      for (bigint i = 0; i < m; i++) array [0][i] -= n3lo;
-      for (int i = 0; i < n1; i++) array[i] -= n2lo;
-      return array;
-    }
-
+  TYPE ****create4d_offset(TYPE ****&array, int n1, int n2lo, int n2hi,
+                           int n3lo, int n3hi, int n4lo, int n4hi,
+                           const char *name)
+  {
+    int n2 = n2hi - n2lo + 1;
+    int n3 = n3hi - n3lo + 1;
+    int n4 = n4hi - n4lo + 1;
+    create(array,n1,n2,n3,n4,name);
+    
+    bigint m = ((bigint) n1) * n2 * n3;
+    for (bigint i = 0; i < m; i++) array[0][0][i] -= n4lo;
+    m = ((bigint) n1) * n2;
+    for (bigint i = 0; i < m; i++) array [0][i] -= n3lo;
+    for (int i = 0; i < n1; i++) array[i] -= n2lo;
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ****create4d_offset(TYPE *****&array, int n1, int n2lo, int n2hi,
-                             int n3lo, int n3hi, int n4lo, int n4hi,
-                             const char *name)
-    {fail(name);}
-
+  TYPE ****create4d_offset(TYPE *****&array, int n1, int n2lo, int n2hi,
+                           int n3lo, int n3hi, int n4lo, int n4hi,
+                           const char *name)
+  {fail(name); return NULL;}
+  
 /* ----------------------------------------------------------------------
    free a 4d array with indices 2,3, and 4 offset
 ------------------------------------------------------------------------- */
+
   template <typename TYPE>
-    void destroy4d_offset(TYPE ****array, int n2_offset, int n3_offset, int n4_offset)
-    {
-      if (array == NULL) return;
-      sfree(&array[0][n2_offset][n3_offset][n4_offset]);
-      sfree(&array[0][n2_offset][n3_offset]);
-      sfree(&array[0][n2_offset]);
-      sfree(array);
-    }
+  void destroy4d_offset(TYPE ****array, 
+                        int n2_offset, int n3_offset, int n4_offset)
+  {
+    if (array == NULL) return;
+    sfree(&array[0][n2_offset][n3_offset][n4_offset]);
+    sfree(&array[0][n2_offset][n3_offset]);
+    sfree(&array[0][n2_offset]);
+    sfree(array);
+  }
 
 /* ----------------------------------------------------------------------
    create a 5d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
-                     int n5, const char *name)
-    {
-      bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4*n5;
-      TYPE *data = (TYPE *) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE *)) * n1*n2*n3*n4;
-      TYPE **level4 = (TYPE **) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE **)) * n1*n2*n3;
-      TYPE ***level3 = (TYPE ***) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE ***)) * n1*n2;
-      TYPE ****level2 = (TYPE ****) smalloc(nbytes,name);
-      nbytes = ((bigint) sizeof(TYPE ****)) * n1;
-      array = (TYPE *****) smalloc(nbytes,name);
-
-      int i,j,k,l;
-      bigint m1,m2;
-      bigint n = 0;
-      for (i = 0; i < n1; i++) {
-        m2 = ((bigint) i) * n2;
-        array[i] = &level2[m2];
-        for (j = 0; j < n2; j++) {
-          m1 = ((bigint) i) * n2 + j;
-          m2 = ((bigint) i) * n2*n3 +  ((bigint) j) * n3;
-          level2[m1] = &level3[m2];
-          for (k = 0; k < n3; k++) {
-            m1 = ((bigint) i) * n2*n3 +  ((bigint) j) * n3 + k;
-            m2 = ((bigint) i) * n2*n3*n4 +
-              ((bigint) j) * n3*n4 + ((bigint) k) * n4;
-            level3[m1] = &level4[m2];
-            for (l = 0; l < n4; l++) {
-              m1 = ((bigint) i) * n2*n3*n4 +
-                ((bigint) j) * n3*n4 + ((bigint) k) * n4 + l;
-              level4[m1] = &data[n];
-              n += n5;
-            }
+  TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
+                   int n5, const char *name)
+  {
+    bigint nbytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4*n5;
+    TYPE *data = (TYPE *) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE *)) * n1*n2*n3*n4;
+    TYPE **level4 = (TYPE **) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE **)) * n1*n2*n3;
+    TYPE ***level3 = (TYPE ***) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE ***)) * n1*n2;
+    TYPE ****level2 = (TYPE ****) smalloc(nbytes,name);
+    nbytes = ((bigint) sizeof(TYPE ****)) * n1;
+    array = (TYPE *****) smalloc(nbytes,name);
+    
+    int i,j,k,l;
+    bigint m1,m2;
+    bigint n = 0;
+    for (i = 0; i < n1; i++) {
+      m2 = ((bigint) i) * n2;
+      array[i] = &level2[m2];
+      for (j = 0; j < n2; j++) {
+        m1 = ((bigint) i) * n2 + j;
+        m2 = ((bigint) i) * n2*n3 +  ((bigint) j) * n3;
+        level2[m1] = &level3[m2];
+        for (k = 0; k < n3; k++) {
+          m1 = ((bigint) i) * n2*n3 +  ((bigint) j) * n3 + k;
+          m2 = ((bigint) i) * n2*n3*n4 +
+            ((bigint) j) * n3*n4 + ((bigint) k) * n4;
+          level3[m1] = &level4[m2];
+          for (l = 0; l < n4; l++) {
+            m1 = ((bigint) i) * n2*n3*n4 +
+              ((bigint) j) * n3*n4 + ((bigint) k) * n4 + l;
+            level4[m1] = &data[n];
+            n += n5;
           }
         }
       }
-      return array;
     }
-
+    return array;
+  }
+  
   template <typename TYPE>
-    TYPE ******create(TYPE ******&array, int n1, int n2, int n3, int n4,
-                      int n5, const char *name)
-    {fail(name);}
+  TYPE ******create(TYPE ******&array, int n1, int n2, int n3, int n4,
+                    int n5, const char *name)
+  {fail(name); return NULL;}
 
 /* ----------------------------------------------------------------------
    destroy a 5d array
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE *****array)
-    {
-      if (array == NULL) return;
-      sfree(array[0][0][0][0]);
-      sfree(array[0][0][0]);
-      sfree(array[0][0]);
-      sfree(array[0]);
-      sfree(array);
-    }
-
+  void destroy(TYPE *****array)
+  {
+    if (array == NULL) return;
+    sfree(array[0][0][0][0]);
+    sfree(array[0][0][0]);
+    sfree(array[0][0]);
+    sfree(array[0]);
+    sfree(array);
+  }
+  
 /* ----------------------------------------------------------------------
    memory usage of arrays, including pointers
 ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    bigint usage(TYPE *array, int n)
-    {
-      bigint bytes = ((bigint) sizeof(TYPE)) * n;
-      return bytes;
-    }
-
+  bigint usage(TYPE *array, int n)
+  {
+    bigint bytes = ((bigint) sizeof(TYPE)) * n;
+    return bytes;
+  }
+  
   template <typename TYPE>
-    bigint usage(TYPE **array, int n1, int n2)
-    {
-      bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2;
-      bytes += ((bigint) sizeof(TYPE *)) * n1;
-      return bytes;
-    }
-
+  bigint usage(TYPE **array, int n1, int n2)
+  {
+    bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2;
+    bytes += ((bigint) sizeof(TYPE *)) * n1;
+    return bytes;
+  }
+  
   template <typename TYPE>
-    bigint usage(TYPE ***array, int n1, int n2, int n3)
-    {
-      bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
-      bytes += ((bigint) sizeof(TYPE *)) * n1*n2;
-      bytes += ((bigint) sizeof(TYPE **)) * n1;
-      return bytes;
-    }
-
+  bigint usage(TYPE ***array, int n1, int n2, int n3)
+  {
+    bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2*n3;
+    bytes += ((bigint) sizeof(TYPE *)) * n1*n2;
+    bytes += ((bigint) sizeof(TYPE **)) * n1;
+    return bytes;
+  }
+  
   template <typename TYPE>
-    bigint usage(TYPE ****array, int n1, int n2, int n3, int n4)
-    {
-      bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4;
-      bytes += ((bigint) sizeof(TYPE *)) * n1*n2*n3;
-      bytes += ((bigint) sizeof(TYPE **)) * n1*n2;
-      bytes += ((bigint) sizeof(TYPE ***)) * n1;
-      return bytes;
-    }
+  bigint usage(TYPE ****array, int n1, int n2, int n3, int n4)
+  {
+    bigint bytes = ((bigint) sizeof(TYPE)) * n1*n2*n3*n4;
+    bytes += ((bigint) sizeof(TYPE *)) * n1*n2*n3;
+    bytes += ((bigint) sizeof(TYPE **)) * n1*n2;
+    bytes += ((bigint) sizeof(TYPE ***)) * n1;
+    return bytes;
+  }
 };
 
 }
diff --git a/src/modify.h b/src/modify.h
index b632c36147..422c77d2fe 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -51,7 +51,6 @@ class Modify : protected Pointers {
   virtual void setup_pre_force(int);
   virtual void initial_integrate(int);
   virtual void post_integrate();
-  void pre_decide();
   virtual void pre_exchange();
   virtual void pre_neighbor();
   virtual void pre_force(int);
@@ -61,27 +60,27 @@ class Modify : protected Pointers {
   virtual double thermo_energy();
   virtual void post_run();
 
-  void setup_pre_force_respa(int, int);
-  void initial_integrate_respa(int, int, int);
-  void post_integrate_respa(int, int);
-  void pre_force_respa(int, int, int);
-  void post_force_respa(int, int, int);
-  void final_integrate_respa(int, int);
+  virtual void setup_pre_force_respa(int, int);
+  virtual void initial_integrate_respa(int, int, int);
+  virtual void post_integrate_respa(int, int);
+  virtual void pre_force_respa(int, int, int);
+  virtual void post_force_respa(int, int, int);
+  virtual void final_integrate_respa(int, int);
 
-  void min_pre_exchange();
-  void min_pre_neighbor();
-  void min_pre_force(int);
-  void min_post_force(int);
+  virtual void min_pre_exchange();
+  virtual void min_pre_neighbor();
+  virtual void min_pre_force(int);
+  virtual void min_post_force(int);
 
-  double min_energy(double *);
-  void min_store();
-  void min_step(double, double *);
-  void min_clearstore();
-  void min_pushstore();
-  void min_popstore();
-  int min_reset_ref();
-  double max_alpha(double *);
-  int min_dof();
+  virtual double min_energy(double *);
+  virtual void min_store();
+  virtual void min_step(double, double *);
+  virtual void min_clearstore();
+  virtual void min_pushstore();
+  virtual void min_popstore();
+  virtual double max_alpha(double *);
+  virtual int min_dof();
+  virtual int min_reset_ref();
 
   void add_fix(int, char **, char *suffix = NULL);
   void modify_fix(int, char **);
diff --git a/src/read_data.cpp b/src/read_data.cpp
index f5a6542280..c563f6b620 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -218,6 +218,15 @@ void ReadData::command(int narg, char **arg)
       domain->box_exist = 1;
       update->ntimestep = 0;
     
+      // insure extra settings are applied before grow(),
+      //   even if no topology in file
+      // if topology is in file, realloc and another grow() is done below
+
+      atom->bond_per_atom = atom->extra_bond_per_atom;
+      atom->angle_per_atom = atom->extra_angle_per_atom;
+      atom->dihedral_per_atom = atom->extra_dihedral_per_atom;
+      atom->improper_per_atom = atom->extra_improper_per_atom;
+
       int n;
       if (comm->nprocs == 1) n = static_cast<int> (atom->natoms);
       else n = static_cast<int> (LB_FACTOR * atom->natoms / comm->nprocs);
diff --git a/src/verlet.h b/src/verlet.h
index f3da03d6c3..7a19880f18 100644
--- a/src/verlet.h
+++ b/src/verlet.h
@@ -39,7 +39,7 @@ class Verlet : public Integrate {
   int torqueflag,erforceflag;
   int e_flag,rho_flag;
 
-  void force_clear();
+  virtual void force_clear();
 };
 
 }