diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp
index 4a98cca61..751560def 100644
--- a/src/pair_hybrid.cpp
+++ b/src/pair_hybrid.cpp
@@ -1,971 +1,974 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, 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.
 ------------------------------------------------------------------------- */
 
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
 #include <ctype.h>
 #include "pair_hybrid.h"
 #include "atom.h"
 #include "force.h"
 #include "pair.h"
 #include "neighbor.h"
 #include "neigh_request.h"
 #include "update.h"
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
 #include "respa.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp),
   styles(NULL), keywords(NULL), multiple(NULL), nmap(NULL),
   map(NULL), special_lj(NULL), special_coul(NULL), compute_tally(NULL)
 {
   nstyles = 0;
 
   outerflag = 0;
   respaflag = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 PairHybrid::~PairHybrid()
 {
   if (nstyles) {
     for (int m = 0; m < nstyles; m++) {
       delete styles[m];
       delete [] keywords[m];
       if (special_lj[m])   delete [] special_lj[m];
       if (special_coul[m]) delete [] special_coul[m];
     }
   }
   delete [] styles;
   delete [] keywords;
   delete [] multiple;
 
   delete [] special_lj;
   delete [] special_coul;
   delete [] compute_tally;
 
   delete [] svector;
 
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
     memory->destroy(cutghost);
     memory->destroy(nmap);
     memory->destroy(map);
   }
 }
 
 /* ----------------------------------------------------------------------
   call each sub-style's compute() or compute_outer() function
   accumulate sub-style global/peratom energy/virial in hybrid
   for global vflag = 1:
     each sub-style computes own virial[6]
     sum sub-style virial[6] to hybrid's virial[6]
   for global vflag = 2:
     call sub-style with adjusted vflag to prevent it calling
       virial_fdotr_compute()
     hybrid calls virial_fdotr_compute() on final accumulated f
 ------------------------------------------------------------------------- */
 
 void PairHybrid::compute(int eflag, int vflag)
 {
   int i,j,m,n;
 
   // if no_virial_fdotr_compute is set and global component of
   //   incoming vflag = 2, then
   // reset vflag as if global component were 1
   // necessary since one or more sub-styles cannot compute virial as F dot r
 
   if (no_virial_fdotr_compute && vflag % 4 == 2) vflag = 1 + vflag/4 * 4;
 
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = eflag_global = vflag_global =
          eflag_atom = vflag_atom = 0;
 
   // check if global component of incoming vflag = 2
   // if so, reset vflag passed to substyle as if it were 0
   // necessary so substyle will not invoke virial_fdotr_compute()
 
   int vflag_substyle;
   if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4;
   else vflag_substyle = vflag;
 
   double *saved_special = save_special();
 
   // check if we are running with r-RESPA using the hybrid keyword
 
   Respa *respa = NULL;
   respaflag = 0;
   if (strstr(update->integrate_style,"respa")) {
     respa = (Respa *) update->integrate;
     if (respa->nhybrid_styles > 0) respaflag = 1;
   }
 
   for (m = 0; m < nstyles; m++) {
 
     set_special(m);
 
     if (!respaflag || (respaflag && respa->hybrid_compute[m])) {
 
       // invoke compute() unless compute flag is turned off or
       // outerflag is set and sub-style has a compute_outer() method
 
       if (styles[m]->compute_flag == 0) continue;
       if (outerflag && styles[m]->respa_enable)
         styles[m]->compute_outer(eflag,vflag_substyle);
       else styles[m]->compute(eflag,vflag_substyle);
     }
 
     restore_special(saved_special);
 
     // jump to next sub-style if r-RESPA does not want global accumulated data
 
     if (respaflag && !respa->tally_global) continue;
 
     if (eflag_global) {
       eng_vdwl += styles[m]->eng_vdwl;
       eng_coul += styles[m]->eng_coul;
     }
     if (vflag_global) {
       for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
     }
     if (eflag_atom) {
       n = atom->nlocal;
       if (force->newton_pair) n += atom->nghost;
       double *eatom_substyle = styles[m]->eatom;
       for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
     }
     if (vflag_atom) {
       n = atom->nlocal;
       if (force->newton_pair) n += atom->nghost;
       double **vatom_substyle = styles[m]->vatom;
       for (i = 0; i < n; i++)
         for (j = 0; j < 6; j++)
           vatom[i][j] += vatom_substyle[i][j];
     }
   }
 
   delete [] saved_special;
 
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::add_tally_callback(Compute *ptr)
 {
   for (int m = 0; m < nstyles; m++)
     if (compute_tally[m]) styles[m]->add_tally_callback(ptr);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::del_tally_callback(Compute *ptr)
 {
   for (int m = 0; m < nstyles; m++)
     if (compute_tally[m]) styles[m]->del_tally_callback(ptr);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::compute_inner()
 {
   for (int m = 0; m < nstyles; m++)
     if (styles[m]->respa_enable) styles[m]->compute_inner();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::compute_middle()
 {
   for (int m = 0; m < nstyles; m++)
     if (styles[m]->respa_enable) styles[m]->compute_middle();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::compute_outer(int eflag, int vflag)
 {
   outerflag = 1;
   compute(eflag,vflag);
   outerflag = 0;
 }
 
 /* ----------------------------------------------------------------------
    allocate all arrays
 ------------------------------------------------------------------------- */
 
 void PairHybrid::allocate()
 {
   allocated = 1;
   int n = atom->ntypes;
 
   memory->create(setflag,n+1,n+1,"pair:setflag");
   for (int i = 1; i <= n; i++)
     for (int j = i; j <= n; j++)
       setflag[i][j] = 0;
 
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
   memory->create(cutghost,n+1,n+1,"pair:cutghost");
 
   memory->create(nmap,n+1,n+1,"pair:nmap");
   memory->create(map,n+1,n+1,nstyles,"pair:map");
   for (int i = 1; i <= n; i++)
     for (int j = i; j <= n; j++)
       nmap[i][j] = 0;
 }
 
 /* ----------------------------------------------------------------------
    create one pair style for each arg in list
 ------------------------------------------------------------------------- */
 
 void PairHybrid::settings(int narg, char **arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal pair_style command");
 
   // delete old lists, since cannot just change settings
 
   if (nstyles) {
     for (int m = 0; m < nstyles; m++) delete styles[m];
     delete [] styles;
     for (int m = 0; m < nstyles; m++) delete [] keywords[m];
     delete [] keywords;
   }
 
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
     memory->destroy(cutghost);
     memory->destroy(nmap);
     memory->destroy(map);
   }
   allocated = 0;
 
   // allocate list of sub-styles as big as possibly needed if no extra args
 
   styles = new Pair*[narg];
   keywords = new char*[narg];
   multiple = new int[narg];
 
   special_lj = new double*[narg];
   special_coul = new double*[narg];
 
   compute_tally = new int[narg];
 
   // allocate each sub-style
   // allocate uses suffix, but don't store suffix version in keywords,
   //   else syntax in coeff() will not match
   // call settings() with set of args that are not pair style names
   // use force->pair_map to determine which args these are
 
   int iarg,jarg,dummy;
 
   iarg = 0;
   nstyles = 0;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"hybrid") == 0)
       error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
     if (strcmp(arg[iarg],"none") == 0)
       error->all(FLERR,"Pair style hybrid cannot have none as an argument");
 
     styles[nstyles] = force->new_pair(arg[iarg],1,dummy);
     force->store_style(keywords[nstyles],arg[iarg],0);
     special_lj[nstyles] = special_coul[nstyles] = NULL;
     compute_tally[nstyles] = 1;
 
     jarg = iarg + 1;
     while (jarg < narg && !force->pair_map->count(arg[jarg])) jarg++;
     styles[nstyles]->settings(jarg-iarg-1,&arg[iarg+1]);
     iarg = jarg;
     nstyles++;
   }
 
   // multiple[i] = 1 to M if sub-style used multiple times, else 0
 
   for (int i = 0; i < nstyles; i++) {
     int count = 0;
     for (int j = 0; j < nstyles; j++) {
       if (strcmp(keywords[j],keywords[i]) == 0) count++;
       if (j == i) multiple[i] = count;
     }
     if (count == 1) multiple[i] = 0;
   }
 
   // set pair flags from sub-style flags
 
   flags();
 }
 
 /* ----------------------------------------------------------------------
    set top-level pair flags from sub-style flags
 ------------------------------------------------------------------------- */
 
 void PairHybrid::flags()
 {
   int m;
 
   // set comm_forward, comm_reverse, comm_reverse_off to max of any sub-style
 
   for (m = 0; m < nstyles; m++) {
     if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
     if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
     if (styles[m]) comm_reverse_off = MAX(comm_reverse_off,
                                           styles[m]->comm_reverse_off);
   }
 
   // single_enable = 1 if any sub-style is set
   // respa_enable = 1 if any sub-style is set
   // manybody_flag = 1 if any sub-style is set
   // no_virial_fdotr_compute = 1 if any sub-style is set
   // ghostneigh = 1 if any sub-style is set
   // ewaldflag, pppmflag, msmflag, dipoleflag, dispersionflag, tip4pflag = 1
   //   if any sub-style is set
   // compute_flag = 1 if any sub-style is set
 
   single_enable = 0;
   compute_flag = 0;
   for (m = 0; m < nstyles; m++) {
     if (styles[m]->single_enable) single_enable = 1;
     if (styles[m]->respa_enable) respa_enable = 1;
     if (styles[m]->manybody_flag) manybody_flag = 1;
     if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1;
     if (styles[m]->ghostneigh) ghostneigh = 1;
     if (styles[m]->ewaldflag) ewaldflag = 1;
     if (styles[m]->pppmflag) pppmflag = 1;
     if (styles[m]->msmflag) msmflag = 1;
     if (styles[m]->dipoleflag) dipoleflag = 1;
     if (styles[m]->dispersionflag) dispersionflag = 1;
     if (styles[m]->tip4pflag) tip4pflag = 1;
     if (styles[m]->compute_flag) compute_flag = 1;
   }
 
   // single_extra = min of all sub-style single_extra
   // allocate svector
 
   single_extra = styles[0]->single_extra;
   for (m = 1; m < nstyles; m++)
     single_extra = MIN(single_extra,styles[m]->single_extra);
 
   if (single_extra) {
     delete [] svector;
     svector = new double[single_extra];
   }
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */
 
 void PairHybrid::coeff(int narg, char **arg)
 {
   if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
+  if (lmp->kokkos)
+    error->all(FLERR,"Cannot yet use pair hybrid with Kokkos");
+
   int ilo,ihi,jlo,jhi;
   force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
   force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
 
   // 3rd arg = pair sub-style name
   // 4th arg = pair sub-style index if name used multiple times
   // allow for "none" as valid sub-style name
 
   int multflag;
   int m;
 
   for (m = 0; m < nstyles; m++) {
     multflag = 0;
     if (strcmp(arg[2],keywords[m]) == 0) {
       if (multiple[m]) {
         multflag = 1;
         if (narg < 4) error->all(FLERR,"Incorrect args for pair coefficients");
         if (!isdigit(arg[3][0]))
           error->all(FLERR,"Incorrect args for pair coefficients");
         int index = force->inumeric(FLERR,arg[3]);
         if (index == multiple[m]) break;
         else continue;
       } else break;
     }
   }
 
   int none = 0;
   if (m == nstyles) {
     if (strcmp(arg[2],"none") == 0) none = 1;
     else error->all(FLERR,"Pair coeff for hybrid has invalid style");
   }
 
   // move 1st/2nd args to 2nd/3rd args
   // if multflag: move 1st/2nd args to 3rd/4th args
   // just copy ptrs, since arg[] points into original input line
 
   arg[2+multflag] = arg[1];
   arg[1+multflag] = arg[0];
 
   // invoke sub-style coeff() starting with 1st remaining arg
 
   if (!none) styles[m]->coeff(narg-1-multflag,&arg[1+multflag]);
 
   // if sub-style only allows one pair coeff call (with * * and type mapping)
   // then unset setflag/map assigned to that style before setting it below
   // in case pair coeff for this sub-style is being called for 2nd time
 
   if (!none && styles[m]->one_coeff)
     for (int i = 1; i <= atom->ntypes; i++)
       for (int j = i; j <= atom->ntypes; j++)
         if (nmap[i][j] && map[i][j][0] == m) {
           setflag[i][j] = 0;
           nmap[i][j] = 0;
         }
 
   // set setflag and which type pairs map to which sub-style
   // if sub-style is none: set hybrid setflag, wipe out map
   // else: set hybrid setflag & map only if substyle setflag is set
   //       previous mappings are wiped out
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       if (none) {
         setflag[i][j] = 1;
         nmap[i][j] = 0;
         count++;
       } else if (styles[m]->setflag[i][j]) {
         setflag[i][j] = 1;
         nmap[i][j] = 1;
         map[i][j][0] = m;
         count++;
       }
     }
   }
 
   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
 }
 
 /* ----------------------------------------------------------------------
    init specific to this pair style
 ------------------------------------------------------------------------- */
 
 void PairHybrid::init_style()
 {
   int i,m,itype,jtype,used,istyle,skip;
 
   // error if a sub-style is not used
 
   int ntypes = atom->ntypes;
 
   for (istyle = 0; istyle < nstyles; istyle++) {
     used = 0;
     for (itype = 1; itype <= ntypes; itype++)
       for (jtype = itype; jtype <= ntypes; jtype++)
         for (m = 0; m < nmap[itype][jtype]; m++)
           if (map[itype][jtype][m] == istyle) used = 1;
     if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used");
   }
 
   // check if special_lj/special_coul overrides are compatible
 
   for (istyle = 0; istyle < nstyles; istyle++) {
     if (special_lj[istyle]) {
       for (i = 1; i < 4; ++i) {
         if (((force->special_lj[i] == 0.0) || (force->special_lj[i] == 1.0))
             && (force->special_lj[i] != special_lj[istyle][i]))
           error->all(FLERR,"Pair_modify special setting for pair hybrid "
                      "incompatible with global special_bonds setting");
       }
     }
 
     if (special_coul[istyle]) {
       for (i = 1; i < 4; ++i) {
         if (((force->special_coul[i] == 0.0)
              || (force->special_coul[i] == 1.0))
             && (force->special_coul[i] != special_coul[istyle][i]))
           error->all(FLERR,"Pair_modify special setting for pair hybrid "
                      "incompatible with global special_bonds setting");
       }
     }
   }
 
   // each sub-style makes its neighbor list request(s)
 
   for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();
 
   // create skip lists inside each pair neigh request
   // any kind of list can have its skip flag set in this loop
 
   for (i = 0; i < neighbor->nrequest; i++) {
     if (!neighbor->requests[i]->pair) continue;
 
     // istyle = associated sub-style for the request
 
     for (istyle = 0; istyle < nstyles; istyle++)
       if (styles[istyle] == neighbor->requests[i]->requestor) break;
 
     // allocate iskip and ijskip
     // initialize so as to skip all pair types
     // set ijskip = 0 if type pair matches any entry in sub-style map
     // set ijskip = 0 if mixing will assign type pair to this sub-style
     //   will occur if type pair is currently unassigned
     //   and both I,I and J,J are assigned to single sub-style
     //   and sub-style for both I,I and J,J match istyle
     // set iskip = 1 only if all ijskip for itype are 1
 
     int *iskip = new int[ntypes+1];
     int **ijskip;
     memory->create(ijskip,ntypes+1,ntypes+1,"pair_hybrid:ijskip");
 
     for (itype = 1; itype <= ntypes; itype++)
       for (jtype = 1; jtype <= ntypes; jtype++)
         ijskip[itype][jtype] = 1;
 
     for (itype = 1; itype <= ntypes; itype++)
       for (jtype = itype; jtype <= ntypes; jtype++) {
         for (m = 0; m < nmap[itype][jtype]; m++)
           if (map[itype][jtype][m] == istyle)
             ijskip[itype][jtype] = ijskip[jtype][itype] = 0;
         if (nmap[itype][jtype] == 0 &&
             nmap[itype][itype] == 1 && map[itype][itype][0] == istyle &&
             nmap[jtype][jtype] == 1 && map[jtype][jtype][0] == istyle)
           ijskip[itype][jtype] = ijskip[jtype][itype] = 0;
       }
 
     for (itype = 1; itype <= ntypes; itype++) {
       iskip[itype] = 1;
       for (jtype = 1; jtype <= ntypes; jtype++)
         if (ijskip[itype][jtype] == 0) iskip[itype] = 0;
     }
 
     // if any skipping occurs
     // set request->skip and copy iskip and ijskip into request
     // else delete iskip and ijskip
     // no skipping if pair style assigned to all type pairs
 
     skip = 0;
     for (itype = 1; itype <= ntypes; itype++)
       for (jtype = 1; jtype <= ntypes; jtype++)
         if (ijskip[itype][jtype] == 1) skip = 1;
 
     if (skip) {
       neighbor->requests[i]->skip = 1;
       neighbor->requests[i]->iskip = iskip;
       neighbor->requests[i]->ijskip = ijskip;
     } else {
       delete [] iskip;
       memory->destroy(ijskip);
     }
   }
 }
 
 /* ----------------------------------------------------------------------
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
 
 double PairHybrid::init_one(int i, int j)
 {
   // if I,J is not set explicitly:
   // perform mixing only if I,I sub-style = J,J sub-style
   // also require I,I and J,J are both assigned to single sub-style
 
   if (setflag[i][j] == 0) {
     if (nmap[i][i] != 1 || nmap[j][j] != 1 || map[i][i][0] != map[j][j][0])
       error->one(FLERR,"All pair coeffs are not set");
     nmap[i][j] = 1;
     map[i][j][0] = map[i][i][0];
   }
 
   // call init/mixing for all sub-styles of I,J
   // set cutsq in sub-style just as Pair::init() does via call to init_one()
   // set cutghost for I,J and J,I just as sub-style does
   // sum tail corrections for I,J
   // return max cutoff of all sub-styles assigned to I,J
   // if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned
 
   double cutmax = 0.0;
   cutghost[i][j] = cutghost[j][i] = 0.0;
   if (tail_flag) etail_ij = ptail_ij = 0.0;
 
   nmap[j][i] = nmap[i][j];
 
   for (int k = 0; k < nmap[i][j]; k++) {
     map[j][i][k] = map[i][j][k];
     double cut = styles[map[i][j][k]]->init_one(i,j);
     styles[map[i][j][k]]->cutsq[i][j] =
       styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
     if (styles[map[i][j][k]]->ghostneigh)
       cutghost[i][j] = cutghost[j][i] =
         MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
     if (tail_flag) {
       etail_ij += styles[map[i][j][k]]->etail_ij;
       ptail_ij += styles[map[i][j][k]]->ptail_ij;
     }
     cutmax = MAX(cutmax,cut);
   }
 
   return cutmax;
 }
 
 /* ----------------------------------------------------------------------
    invoke setup for each sub-style
 ------------------------------------------------------------------------- */
 
 void PairHybrid::setup()
 {
   for (int m = 0; m < nstyles; m++) styles[m]->setup();
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes to restart file
 ------------------------------------------------------------------------- */
 
 void PairHybrid::write_restart(FILE *fp)
 {
   fwrite(&nstyles,sizeof(int),1,fp);
 
   // each sub-style writes its settings, but no coeff info
 
   int n;
   for (int m = 0; m < nstyles; m++) {
     n = strlen(keywords[m]) + 1;
     fwrite(&n,sizeof(int),1,fp);
     fwrite(keywords[m],sizeof(char),n,fp);
     styles[m]->write_restart_settings(fp);
     // write out per style special settings, if present
     n = (special_lj[m] == NULL) ? 0 : 1;
     fwrite(&n,sizeof(int),1,fp);
     if (n) fwrite(special_lj[m],sizeof(double),4,fp);
     n = (special_coul[m] == NULL) ? 0 : 1;
     fwrite(&n,sizeof(int),1,fp);
     if (n) fwrite(special_coul[m],sizeof(double),4,fp);
   }
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads from restart file, bcasts
 ------------------------------------------------------------------------- */
 
 void PairHybrid::read_restart(FILE *fp)
 {
   int me = comm->me;
   if (me == 0) fread(&nstyles,sizeof(int),1,fp);
   MPI_Bcast(&nstyles,1,MPI_INT,0,world);
 
   // allocate list of sub-styles
 
   styles = new Pair*[nstyles];
   keywords = new char*[nstyles];
   multiple = new int[nstyles];
 
   special_lj = new double*[nstyles];
   special_coul = new double*[nstyles];
 
   // each sub-style is created via new_pair()
   // each reads its settings, but no coeff info
 
   int n,dummy;
   for (int m = 0; m < nstyles; m++) {
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     keywords[m] = new char[n];
     if (me == 0) fread(keywords[m],sizeof(char),n,fp);
     MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
     styles[m] = force->new_pair(keywords[m],0,dummy);
     styles[m]->read_restart_settings(fp);
     // read back per style special settings, if present
     special_lj[m] = special_coul[m] = NULL;
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     if (n > 0 ) {
       special_lj[m] = new double[4];
       if (me == 0) fread(special_lj[m],sizeof(double),4,fp);
       MPI_Bcast(special_lj[m],4,MPI_DOUBLE,0,world);
     }
     if (me == 0) fread(&n,sizeof(int),1,fp);
     MPI_Bcast(&n,1,MPI_INT,0,world);
     if (n > 0 ) {
       special_coul[m] = new double[4];
       if (me == 0) fread(special_coul[m],sizeof(double),4,fp);
       MPI_Bcast(special_coul[m],4,MPI_DOUBLE,0,world);
     }
   }
 
   // multiple[i] = 1 to M if sub-style used multiple times, else 0
 
   for (int i = 0; i < nstyles; i++) {
     int count = 0;
     for (int j = 0; j < nstyles; j++) {
       if (strcmp(keywords[j],keywords[i]) == 0) count++;
       if (j == i) multiple[i] = count;
     }
     if (count == 1) multiple[i] = 0;
   }
 
   // set pair flags from sub-style flags
 
   flags();
 }
 
 /* ----------------------------------------------------------------------
    call sub-style to compute single interaction
    error if sub-style does not support single() call
    since overlay could have multiple sub-styles, sum results explicitly
 ------------------------------------------------------------------------- */
 
 double PairHybrid::single(int i, int j, int itype, int jtype,
                           double rsq, double factor_coul, double factor_lj,
                           double &fforce)
 {
   if (nmap[itype][jtype] == 0)
     error->one(FLERR,"Invoked pair single on pair style none");
 
   double fone;
   fforce = 0.0;
   double esum = 0.0;
 
   for (int m = 0; m < nmap[itype][jtype]; m++) {
     if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) {
       if (styles[map[itype][jtype][m]]->single_enable == 0)
         error->one(FLERR,"Pair hybrid sub-style does not support single call");
 
       if ((special_lj[map[itype][jtype][m]] != NULL) ||
           (special_coul[map[itype][jtype][m]] != NULL))
         error->one(FLERR,"Pair hybrid single calls do not support"
                    " per sub-style special bond values");
 
       esum += styles[map[itype][jtype][m]]->
         single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone);
       fforce += fone;
 
       // copy substyle extra values into hybrid's svector
 
       if (single_extra && styles[map[itype][jtype][m]]->single_extra)
         for (m = 0; m < single_extra; m++)
           svector[m] = styles[map[itype][jtype][m]]->svector[m];
     }
   }
 
   return esum;
 }
 
 /* ----------------------------------------------------------------------
    modify parameters of the pair style and its sub-styles
 ------------------------------------------------------------------------- */
 
 void PairHybrid::modify_params(int narg, char **arg)
 {
   if (narg == 0) error->all(FLERR,"Illegal pair_modify command");
 
   // if 1st keyword is pair, apply other keywords to one sub-style
 
   if (strcmp(arg[0],"pair") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal pair_modify command");
     int m;
     for (m = 0; m < nstyles; m++)
       if (strcmp(arg[1],keywords[m]) == 0) break;
     if (m == nstyles) error->all(FLERR,"Unknown pair_modify hybrid sub-style");
     int iarg = 2;
 
     if (multiple[m]) {
       if (narg < 3) error->all(FLERR,"Illegal pair_modify command");
       int multiflag = force->inumeric(FLERR,arg[2]);
       for (m = 0; m < nstyles; m++)
         if (strcmp(arg[1],keywords[m]) == 0 && multiflag == multiple[m]) break;
       if (m == nstyles)
         error->all(FLERR,"Unknown pair_modify hybrid sub-style");
       iarg = 3;
     }
 
     // if 2nd keyword (after pair) is special:
     // invoke modify_special() for the sub-style
 
     if (iarg < narg && strcmp(arg[iarg],"special") == 0) {
       if (narg < iarg+5)
         error->all(FLERR,"Illegal pair_modify special command");
       modify_special(m,narg-iarg,&arg[iarg+1]);
       iarg += 5;
     }
 
     // if 2nd keyword (after pair) is compute/tally:
     // set flag to register USER-TALLY computes accordingly
 
     if (iarg < narg && strcmp(arg[iarg],"compute/tally") == 0) {
       if (narg < iarg+2)
         error->all(FLERR,"Illegal pair_modify compute/tally command");
       if (strcmp(arg[iarg+1],"yes") == 0) {
         compute_tally[m] = 1;
       } else if (strcmp(arg[iarg+1],"no") == 0) {
         compute_tally[m] = 0;
       } else error->all(FLERR,"Illegal pair_modify compute/tally command");
       iarg += 2;
     }
 
     // apply the remaining keywords to the base pair style itself and the
     // sub-style except for "pair" and "special".
     // the former is important for some keywords like "tail" or "compute"
 
     if (narg-iarg > 0) {
       Pair::modify_params(narg-iarg,&arg[iarg]);
       styles[m]->modify_params(narg-iarg,&arg[iarg]);
     }
 
   // apply all keywords to pair hybrid itself and every sub-style
 
   } else {
     Pair::modify_params(narg,arg);
     for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
   }
 
   // reset global compute_flag since there may have been changes
   // to any of the substyles
   compute_flag = 0;
   for (int m = 0; m < nstyles; m++)
     if (styles[m]->compute_flag) compute_flag = 1;
 }
 
 /* ----------------------------------------------------------------------
    store a local per pair style override for special_lj and special_coul
 ------------------------------------------------------------------------- */
 
 void PairHybrid::modify_special(int m, int narg, char **arg)
 {
   double special[4];
   int i;
 
   special[0] = 1.0;
   special[1] = force->numeric(FLERR,arg[1]);
   special[2] = force->numeric(FLERR,arg[2]);
   special[3] = force->numeric(FLERR,arg[3]);
 
   if (strcmp(arg[0],"lj/coul") == 0) {
     if (!special_lj[m]) special_lj[m] = new double[4];
     if (!special_coul[m]) special_coul[m] = new double[4];
     for (i = 0; i < 4; ++i)
       special_lj[m][i] = special_coul[m][i] = special[i];
 
   } else if (strcmp(arg[0],"lj") == 0) {
     if (!special_lj[m]) special_lj[m] = new double[4];
     for (i = 0; i < 4; ++i)
       special_lj[m][i] = special[i];
 
   } else if (strcmp(arg[0],"coul") == 0) {
     if (!special_coul[m]) special_coul[m] = new double[4];
     for (i = 0; i < 4; ++i)
       special_coul[m][i] = special[i];
 
   } else error->all(FLERR,"Illegal pair_modify special command");
 }
 
 /* ----------------------------------------------------------------------
    override global special bonds settings with per substyle values
 ------------------------------------------------------------------------- */
 
 void PairHybrid::set_special(int m)
 {
   int i;
   if (special_lj[m])
     for (i = 0; i < 4; ++i) force->special_lj[i] = special_lj[m][i];
   if (special_coul[m])
     for (i = 0; i < 4; ++i) force->special_coul[i] = special_coul[m][i];
 }
 
 /* ----------------------------------------------------------------------
    store global special settings
 ------------------------------------------------------------------------- */
 
 double * PairHybrid::save_special()
 {
   double *saved = new double[8];
 
   for (int i = 0; i < 4; ++i) {
     saved[i] = force->special_lj[i];
     saved[i+4] = force->special_coul[i];
   }
   return saved;
 }
 
 /* ----------------------------------------------------------------------
    restore global special settings from saved data
 ------------------------------------------------------------------------- */
 
 void PairHybrid::restore_special(double *saved)
 {
   for (int i = 0; i < 4; ++i) {
     force->special_lj[i] = saved[i];
     force->special_coul[i] = saved[i+4];
   }
 }
 
 /* ----------------------------------------------------------------------
    extract a ptr to a particular quantity stored by pair
    pass request thru to sub-styles
    return first non-NULL result except for cut_coul request
    for cut_coul, insure all non-NULL results are equal since required by Kspace
 ------------------------------------------------------------------------- */
 
 void *PairHybrid::extract(const char *str, int &dim)
 {
   void *cutptr = NULL;
   void *ptr;
   double cutvalue = 0.0;
 
   for (int m = 0; m < nstyles; m++) {
     ptr = styles[m]->extract(str,dim);
     if (ptr && strcmp(str,"cut_coul") == 0) {
       double *p_newvalue = (double *) ptr;
       double newvalue = *p_newvalue;
       if (cutptr && newvalue != cutvalue)
         error->all(FLERR,
                    "Coulomb cutoffs of pair hybrid sub-styles do not match");
       cutptr = ptr;
       cutvalue = newvalue;
     } else if (ptr) return ptr;
   }
 
   if (strcmp(str,"cut_coul") == 0) return cutptr;
   return NULL;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairHybrid::reset_dt()
 {
   for (int m = 0; m < nstyles; m++) styles[m]->reset_dt();
 }
 
 /* ----------------------------------------------------------------------
    check if itype,jtype maps to sub-style
 ------------------------------------------------------------------------- */
 
 int PairHybrid::check_ijtype(int itype, int jtype, char *substyle)
 {
   for (int m = 0; m < nmap[itype][jtype]; m++)
     if (strcmp(keywords[map[itype][jtype][m]],substyle) == 0) return 1;
   return 0;
 }
 
 /* ----------------------------------------------------------------------
    memory usage of each sub-style
 ------------------------------------------------------------------------- */
 
 double PairHybrid::memory_usage()
 {
   double bytes = maxeatom * sizeof(double);
   bytes += maxvatom*6 * sizeof(double);
   for (int m = 0; m < nstyles; m++) bytes += styles[m]->memory_usage();
   return bytes;
 }
diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h
index 2364b16f4..463ae00ec 100644
--- a/src/pair_hybrid.h
+++ b/src/pair_hybrid.h
@@ -1,154 +1,158 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, 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.
 ------------------------------------------------------------------------- */
 
 #ifdef PAIR_CLASS
 
 PairStyle(hybrid,PairHybrid)
 
 #else
 
 #ifndef LMP_PAIR_HYBRID_H
 #define LMP_PAIR_HYBRID_H
 
 #include <stdio.h>
 #include "pair.h"
 
 namespace LAMMPS_NS {
 
 class PairHybrid : public Pair {
   friend class FixGPU;
   friend class FixIntel;
   friend class FixOMP;
   friend class Force;
   friend class Respa;
   friend class Info;
  public:
   PairHybrid(class LAMMPS *);
   virtual ~PairHybrid();
   virtual void compute(int, int);
   void settings(int, char **);
   virtual void coeff(int, char **);
   void init_style();
   double init_one(int, int);
   void setup();
   void write_restart(FILE *);
   void read_restart(FILE *);
   double single(int, int, int, int, double, double, double, double &);
   void modify_params(int narg, char **arg);
   double memory_usage();
 
   void compute_inner();
   void compute_middle();
   void compute_outer(int, int);
   void *extract(const char *, int &);
   void reset_dt();
 
   int check_ijtype(int, int, char *);
 
   virtual void add_tally_callback(class Compute *);
   virtual void del_tally_callback(class Compute *);
 
  protected:
   int nstyles;                  // # of sub-styles
   Pair **styles;                // list of Pair style classes
   char **keywords;              // style name of each Pair style
   int *multiple;                // 0 if style used once, else Mth instance
 
   int outerflag;                // toggle compute() when invoked by outer()
   int respaflag;                // 1 if different substyles are assigned to
                                 // different r-RESPA levels
 
   int **nmap;                   // # of sub-styles itype,jtype points to
   int ***map;                   // list of sub-styles itype,jtype points to
   double **special_lj;          // list of per style LJ exclusion factors
   double **special_coul;        // list of per style Coulomb exclusion factors
   int *compute_tally;           // list of on/off flags for tally computes
 
   void allocate();
   void flags();
 
   void modify_special(int, int, char**);
   double *save_special();
   void set_special(int);
   void restore_special(double *);
 };
 
 }
 
 #endif
 #endif
 
 /* ERROR/WARNING messages:
 
+E: Cannot yet use pair hybrid with Kokkos
+
+This feature is not yet supported.
+
 E: Illegal ... command
 
 Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
 E: Pair style hybrid cannot have hybrid as an argument
 
 Self-explanatory.
 
 E: Pair style hybrid cannot have none as an argument
 
 Self-explanatory.
 
 E: Incorrect args for pair coefficients
 
 Self-explanatory.  Check the input script or data file.
 
 E: Pair coeff for hybrid has invalid style
 
 Style in pair coeff must have been listed in pair_style command.
 
 E: Pair hybrid sub-style is not used
 
 No pair_coeff command used a sub-style specified in the pair_style
 command.
 
 E: Pair_modify special setting for pair hybrid incompatible with global special_bonds setting
 
 Cannot override a setting of 0.0 or 1.0 or change a setting between
 0.0 and 1.0.
 
 E: All pair coeffs are not set
 
 All pair coefficients must be set in the data file or by the
 pair_coeff command before running a simulation.
 
 E: Invoked pair single on pair style none
 
 A command (e.g. a dump) attempted to invoke the single() function on a
 pair style none, which is illegal.  You are probably attempting to
 compute per-atom quantities with an undefined pair style.
 
 E: Pair hybrid sub-style does not support single call
 
 You are attempting to invoke a single() call on a pair style
 that doesn't support it.
 
 E: Pair hybrid single calls do not support per sub-style special bond values
 
 Self-explanatory.
 
 E: Unknown pair_modify hybrid sub-style
 
 The choice of sub-style is unknown.
 
 E: Coulomb cutoffs of pair hybrid sub-styles do not match
 
 If using a Kspace solver, all Coulomb cutoffs of long pair styles must
 be the same.
 
 */