diff --git a/src/USER-REAXC/fix_reaxc_species.cpp b/src/USER-REAXC/fix_reaxc_species.cpp
index 56dd17f828..b370814933 100644
--- a/src/USER-REAXC/fix_reaxc_species.cpp
+++ b/src/USER-REAXC/fix_reaxc_species.cpp
@@ -12,7 +12,8 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Ray Shan (Sandia, tnshan@sandia.gov)
+   Contributing authors: Ray Shan (Sandia, tnshan@sandia.gov)
+   			 Oleg Sergeev (VNIIA, sergeev@vniia.ru)
 ------------------------------------------------------------------------- */
 
 #include "lmptype.h"
@@ -20,9 +21,10 @@
 #include "math.h"
 #include "atom.h"
 #include "string.h"
+#include "fix_ave_atom.h"
 #include "fix_reaxc_species.h"
-#include "update.h"
 #include "domain.h"
+#include "update.h"
 #include "pair_reax_c.h"
 #include "modify.h"
 #include "neighbor.h"
@@ -36,8 +38,6 @@
 #include "memory.h"
 #include "error.h"
 #include "reaxc_list.h"
-#include "reaxc_types.h"
-#include "reaxc_defs.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -54,32 +54,59 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 2;
 
+  peratom_flag = 1;
+  size_peratom_cols = 0;
+  peratom_freq = 1;
+
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
   ntypes = atom->ntypes;
 
-  nevery = force->inumeric(FLERR,arg[3]);
-  nrepeat = force->inumeric(FLERR,arg[4]);
-  global_freq = nfreq = force->inumeric(FLERR,arg[5]);
+  nevery = atoi(arg[3]);
+  nrepeat = atoi(arg[4]);
+  global_freq = nfreq = atoi(arg[5]);
 
   comm_forward = 1;
   
-  if (nevery == 0 || nrepeat <= 0 || nfreq <= 0)
+  if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
     error->all(FLERR,"Illegal fix reax/c/species command");
   if (nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all(FLERR,"Illegal fix reax/c/species command");
+  
+  // Neighbor lists must stay unchanged during averaging of bonds, 
+  // but may be updated when no averaging is performed.
+  
+  int rene_flag = 0;
+  if (nfreq % neighbor->every != 0 || neighbor->every < nevery * nrepeat) {
+    int newneighborevery = nevery * nrepeat;
+    while (nfreq % newneighborevery != 0 && newneighborevery <= nfreq / 2)
+      newneighborevery++;
 
-  if (neighbor->every != nfreq || neighbor->delay != 0 || neighbor->dist_check != 0){
-    if (me == 0) {
-      char str[128];
-      sprintf(str,"Resetting reneighboring criteria for fix reax/c/species");
-      error->warning(FLERR,str);
-    }
-    neighbor->every = nfreq;
+    if (nfreq % newneighborevery != 0)
+      newneighborevery = nfreq;
+
+    neighbor->every = newneighborevery;
+    rene_flag = 1;
+  }
+
+  if (neighbor->delay != 0 || neighbor->dist_check != 0) {
     neighbor->delay = 0;
     neighbor->dist_check = 0;
+    rene_flag = 1;
   }
 
+  if (me == 0 && rene_flag) {
+    char str[128];
+    sprintf(str,"Resetting reneighboring criteria for fix reax/c/species");
+    error->warning(FLERR,str);
+  }
+
+  tmparg = NULL;
+  memory->create(tmparg,4,4,"reax/c/species:tmparg");
+  strcpy(tmparg[0],arg[3]);
+  strcpy(tmparg[1],arg[4]);
+  strcpy(tmparg[2],arg[5]);
+
   if (me == 0) {
     fp = fopen(arg[6],"w");
     if (fp == NULL) {
@@ -89,12 +116,17 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
-  tmpq = NULL;
-  tmpx = NULL;
-  abo = NULL;
-  BOCut = NULL;
+  x0 = NULL;
+  PBCconnected = NULL;
   clusterID = NULL;
 
+  int ntmp = 1;
+  memory->create(x0,ntmp,"reax/c/species:x0");
+  memory->create(PBCconnected,ntmp,"reax/c/species:PBCconnected");
+  memory->create(clusterID,ntmp,"reax/c/species:clusterID");
+  vector_atom = clusterID;
+
+  BOCut = NULL;
   Name = NULL;
   MolName = NULL;
   MolType = NULL;
@@ -102,43 +134,48 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
   nd = NULL;
   molmap = NULL;
 
+  nmax = 0;
+  setupflag = 0;
+
+  // set default bond order cutoff
+  int n, i, j, itype, jtype;
+  double bo_cut;
+  bg_cut = 0.30;
+  n = ntypes+1;
+  memory->create(BOCut,n,n,"reax/c/species:BOCut");
+  for (i = 1; i < n; i ++)
+    for (j = 1; j < n; j ++)
+      BOCut[i][j] = bg_cut;
+
+  // optional args
+  eletype = NULL;
+  ele = filepos = NULL;
+  eleflag = posflag = padflag = 0;
+
   singlepos_opened = multipos_opened = 0;
   multipos = 0;
   posfreq = 0;
 
-  // initialize bond order cutoff
-  int n, i, j;
-  n = ntypes+1;
-  memory->create(BOCut,n,n,"reaxc/c/species:BOCut");
-  for (i = 1; i < n; i ++)
-    for (j = 1; j < n; j ++)
-      BOCut[i][j] = 0.0;
-
-  // optional args
-  eletype = NULL;
-  ele = posspec = filepos = NULL;
-  eleflag = posflag = padflag = 0;
-
   int iarg = 7;
-  int itype, jtype;
-  double bo_cut;
-
   while (iarg < narg) {
+
     // set BO cutoff
     if (strcmp(arg[iarg],"cutoff") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix reax/c/species command");
-      itype = force->inumeric(FLERR,arg[iarg+1]);
-      jtype = force->inumeric(FLERR,arg[iarg+2]);
-      bo_cut = force->numeric(FLERR,arg[iarg+3]);
+      itype = atoi(arg[iarg+1]);
+      jtype = atoi(arg[iarg+2]);
+      bo_cut = atof(arg[iarg+3]);
       if (itype > ntypes || jtype > ntypes) 
-	error->all(FLERR,"Illegal fix reax/c/species command");
+      	error->all(FLERR,"Illegal fix reax/c/species command");
       if (itype <= 0 || jtype <= 0) 
-	error->all(FLERR,"Illegal fix reax/c/species command");
+      	error->all(FLERR,"Illegal fix reax/c/species command");
       if (bo_cut > 1.0 || bo_cut < 0.0)
-	error->all(FLERR,"Illegal fix reax/c/species command");
+      	error->all(FLERR,"Illegal fix reax/c/species command");
+
       BOCut[itype][jtype] = bo_cut; 
       BOCut[jtype][itype] = bo_cut;
       iarg += 4;
+
     // modify element type names
     } else if (strcmp(arg[iarg],"element") == 0) {
       if (iarg+ntypes+1 > narg) error->all(FLERR,"Illegal fix reax/c/species command");
@@ -146,26 +183,30 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
       eletype = (char**) malloc(ntypes*sizeof(char*));
       for (int i = 0; i < ntypes; i ++) {
         eletype[i] = (char*) malloc(2*sizeof(char));
-	strcpy(eletype[i],arg[iarg+1+i]);
+      	strcpy(eletype[i],arg[iarg+1+i]);
       }
       eleflag = 1;
       iarg += ntypes + 1;
+
     // position of molecules
     } else if (strcmp(arg[iarg],"position") == 0) {
-      if (iarg+3 > narg) error->all(FLERR,"Illegal fix species command");
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix reax/c/species command");
       posflag = 1;
-      posfreq = force->inumeric(FLERR,arg[iarg+1]);
-      filepos = new char[n];
+      posfreq = atoi(arg[iarg+1]);
+      if (posfreq < nfreq || (posfreq%nfreq != 0))
+	error->all(FLERR,"Illegal fix reax/c/species command");
+
+      filepos = new char[255];
       strcpy(filepos,arg[iarg+2]);
       if (strchr(filepos,'*')) {
-	multipos = 1;
+        multipos = 1;
       } else {
-	if (me == 0) {
-	  pos = fopen(filepos, "w");
-	  if (pos == NULL) error->one(FLERR,"Cannot open fix species position file");
-	}
-	singlepos_opened = 1;
-	multipos = 0;
+        if (me == 0) {
+          pos = fopen(filepos, "w");
+          if (pos == NULL) error->one(FLERR,"Cannot open fix reax/c/species position file");
+        }
+      	singlepos_opened = 1;
+      	multipos = 0;
       }
       iarg += 3;
     } else error->all(FLERR,"Illegal fix reax/c/species command");
@@ -174,18 +215,16 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
   if (!eleflag) {
     memory->create(ele,ntypes+1,"reax/c/species:ele");
     ele[0]='C';
-    ele[1]='H';
-    ele[2]='O';
-    ele[3]='N';
+    if (ntypes > 1)
+      ele[1]='H';
+    if (ntypes > 2)
+      ele[2]='O';
+    if (ntypes > 3)
+      ele[3]='N';
   }
 
-  nmax = 0;
-  vector_nmole = vector_nspec = 0;
-
-  irepeat = 0;
-  nvalid = nextvalid();
-
-  memory->create(Name,ntypes,"reax/c/species:Name");
+  vector_nmole = 0;
+  vector_nspec = 0;
 
 }
 
@@ -193,15 +232,27 @@ FixReaxCSpecies::FixReaxCSpecies(LAMMPS *lmp, int narg, char **arg) :
 
 FixReaxCSpecies::~FixReaxCSpecies()
 {
-  memory->destroy(tmpq);
-  memory->destroy(tmpx);
-  memory->destroy(abo);
-  memory->destroy(Name);
+  memory->destroy(ele);
   memory->destroy(BOCut);
   memory->destroy(clusterID);
+  memory->destroy(PBCconnected);
+  memory->destroy(x0);
+
+  memory->destroy(nd);
+  memory->destroy(Name);
+  memory->destroy(NMol);
+  memory->destroy(MolType);
+  memory->destroy(MolName);
+  memory->destroy(tmparg);
+
+  if (filepos)
+    delete [] filepos;
 
   if (me == 0) fclose(fp);
   if (me == 0 && posflag && multipos_opened) fclose(pos);
+
+  modify->delete_compute("SPECATOM");
+  modify->delete_fix("SPECBOND");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -218,6 +269,7 @@ int FixReaxCSpecies::setmask()
 void FixReaxCSpecies::setup(int vflag)
 {
   ntotal = static_cast<int> (atom->natoms);
+  memory->create(Name,ntypes,"reax/c/species:Name");
 
   post_integrate();
 }
@@ -227,30 +279,130 @@ void FixReaxCSpecies::setup(int vflag)
 void FixReaxCSpecies::init()
 {
   if (atom->tag_enable == 0)
-    error->all(FLERR,"Cannot use fix reax/c/specis unless atoms have IDs");
+    error->all(FLERR,"Cannot use fix reax/c/species unless atoms have IDs");
 
   reaxc = (PairReaxC *) force->pair_match("reax/c",1);
-  if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/specis without "
+  if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/species without "
 		  "pair_style reax/c");
 
-  if (nvalid < update->ntimestep) {
-    irepeat = 0;
-    nvalid = nextvalid();
-  }
+  reaxc->fixspecies_flag = 1;
+  nvalid = update->ntimestep+nfreq;
 
+  // check if this fix has been called twice
   int count = 0;
   for (int i = 0; i < modify->nfix; i++)
     if (strcmp(modify->fix[i]->style,"reax/c/species") == 0) count++;
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one fix reax/c/species");
 
-  // set default bond order cutoff
-  int n, i, j;
-  bg_cut = reaxc->control->bg_cut;
-  n = ntypes+1;
-  for (i = 1; i < n; i ++)
-    for (j = 1; j < n; j ++)
-      if (BOCut[i][j] == 0.0) BOCut[i][j] = bg_cut;
+  if (!setupflag) {
+    // create a compute to store properties
+    create_compute();
+
+    // create a fix to point to fix_ave_atom for averaging stored properties
+    create_fix();
+
+    setupflag = 1;
+  }
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpecies::create_compute()
+{
+  int narg;
+  char **args;
+
+  narg = 34;
+  args = new char*[narg];
+  args[0]  = (char *) "SPECATOM";
+  args[1]  = (char *) "all";
+  args[2]  = (char *) "spec/atom";
+  args[3]  = (char *) "q";
+  args[4]  = (char *) "x";
+  args[5]  = (char *) "y";
+  args[6]  = (char *) "z";
+  args[7]  = (char *) "vx";
+  args[8]  = (char *) "vy";
+  args[9]  = (char *) "vz";
+  args[10] = (char *) "abo01";
+  args[11] = (char *) "abo02";
+  args[12] = (char *) "abo03";
+  args[13] = (char *) "abo04";
+  args[14] = (char *) "abo05";
+  args[15] = (char *) "abo06";
+  args[16] = (char *) "abo07";
+  args[17] = (char *) "abo08";
+  args[18] = (char *) "abo09";
+  args[19] = (char *) "abo10";
+  args[20] = (char *) "abo11";
+  args[21] = (char *) "abo12";
+  args[22] = (char *) "abo13";
+  args[23] = (char *) "abo14";
+  args[24] = (char *) "abo15";
+  args[25] = (char *) "abo16";
+  args[26] = (char *) "abo17";
+  args[27] = (char *) "abo18";
+  args[28] = (char *) "abo19";
+  args[29] = (char *) "abo20";
+  args[30] = (char *) "abo21";
+  args[31] = (char *) "abo22";
+  args[32] = (char *) "abo23";
+  args[33] = (char *) "abo24";
+  modify->add_compute(narg,args);
+  delete [] args;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpecies::create_fix()
+{
+  int narg;
+  char **args;
+
+  narg = 37;
+  args = new char*[narg];
+  args[0]  = (char *) "SPECBOND";
+  args[1]  = (char *) "all";
+  args[2]  = (char *) "ave/atom";
+  args[3]  = tmparg[0];
+  args[4]  = tmparg[1];
+  args[5]  = tmparg[2];
+  args[6]  = (char *) "c_SPECATOM[1]";	 // q, array_atoms[i][0]
+  args[7]  = (char *) "c_SPECATOM[2]";	 // x, 1
+  args[8]  = (char *) "c_SPECATOM[3]";	 // y, 2
+  args[9]  = (char *) "c_SPECATOM[4]";	 // z, 3
+  args[10] = (char *) "c_SPECATOM[5]";	 // vx, 4
+  args[11] = (char *) "c_SPECATOM[6]";	 // vy, 5
+  args[12] = (char *) "c_SPECATOM[7]";	 // vz, 6
+  args[13] = (char *) "c_SPECATOM[8]";	 // abo01, 7
+  args[14] = (char *) "c_SPECATOM[9]";	
+  args[15] = (char *) "c_SPECATOM[10]";	
+  args[16] = (char *) "c_SPECATOM[11]"; 
+  args[17] = (char *) "c_SPECATOM[12]";
+  args[18] = (char *) "c_SPECATOM[13]";
+  args[19] = (char *) "c_SPECATOM[14]"; 
+  args[20] = (char *) "c_SPECATOM[15]";
+  args[21] = (char *) "c_SPECATOM[16]";
+  args[22] = (char *) "c_SPECATOM[17]";
+  args[23] = (char *) "c_SPECATOM[18]";
+  args[24] = (char *) "c_SPECATOM[19]"; // abo12, 18
+  args[25] = (char *) "c_SPECATOM[20]";	
+  args[26] = (char *) "c_SPECATOM[21]"; 
+  args[27] = (char *) "c_SPECATOM[22]";
+  args[28] = (char *) "c_SPECATOM[23]";
+  args[29] = (char *) "c_SPECATOM[24]"; 
+  args[30] = (char *) "c_SPECATOM[25]";
+  args[31] = (char *) "c_SPECATOM[26]";
+  args[32] = (char *) "c_SPECATOM[27]";
+  args[33] = (char *) "c_SPECATOM[28]";
+  args[34] = (char *) "c_SPECATOM[29]"; 
+  args[35] = (char *) "c_SPECATOM[30]";
+  args[36] = (char *) "c_SPECATOM[31]";
+  modify->add_fix(narg,args);
+  f_SPECBOND = (FixAveAtom *) modify->fix[modify->nfix-1];
+  delete [] args;
 
 }
 
@@ -265,125 +417,81 @@ void FixReaxCSpecies::init_list(int id, NeighList *ptr)
 
 void FixReaxCSpecies::post_integrate()
 {
-  bigint ntimestep = update->ntimestep;
-  if (ntimestep != nvalid) return;
-
   Output_ReaxC_Bonds(update->ntimestep,fp);
   if (me == 0) fflush(fp);
-
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixReaxCSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
+
 {
-  int i, j, ii,jj;
   int Nmole, Nspec;
 
-  MPI_Barrier(world);
+  // point to fix_ave_atom
+  f_SPECBOND->end_of_step();
+
+  if (ntimestep != nvalid) return;
+
   nlocal = atom->nlocal;
 
   if (atom->nmax > nmax) {
     nmax = atom->nmax;
-    memory->destroy(tmpq);
-    memory->destroy(tmpx);
-    memory->destroy(abo);
+    memory->destroy(x0);
+    memory->destroy(PBCconnected);
     memory->destroy(clusterID);
-    memory->create(tmpq,nmax,"reax/c/species:tmpq");
-    memory->create(tmpx,nmax,3,"reax/c/species:tmpx");
-    memory->create(abo,nmax,nmax,"reax/c/species:abo");
+    memory->create(x0,nmax,"reax/c/species:x0");
+    memory->create(PBCconnected,nmax,"reax/c/species:PBCconnected");
     memory->create(clusterID,nmax,"reax/c/species:clusterID");
+    vector_atom = clusterID;
+  }
+
+  for (int i = 0; i < nmax; i++) {
+    PBCconnected[i] = 0;
+    x0[i].x = x0[i].y = x0[i].z = 0.0;
   }
 
-  repeat = nrepeat;
   Nmole = Nspec = 0;
-  vector_nmole = vector_nspec = 0;
-
-  if (irepeat == 0)
-    for (i = 0; i < nmax; i++) {
-      tmpq[i] = 0.0;
-      for (j = 0; j < nmax; j++)
-        abo[i][j] = 0.0;
-      for (j = 0; j < 3; j++)
-        tmpx[i][j] = 0.0;
-    }
-
-  GatherBondOrder(lists);
-
-  irepeat++;
-  if (irepeat < nrepeat) {
-    nvalid += nevery;
-    return;
-  }
-  irepeat = 0;
-  nvalid = ntimestep + nfreq - (nrepeat-1)*nevery;
 
   FindMolecule();
 
-  SortMolecule( Nmole);
+  SortMolecule (Nmole);
 
   FindSpecies(Nmole, Nspec);
 
   vector_nmole = Nmole;
   vector_nspec = Nspec;
 
-  if (me == 0) WriteFormulae( Nmole, Nspec);
+  if (me == 0 && ntimestep >= 0)
+    WriteFormulas (Nmole, Nspec);
 
-  if (posflag && (ntimestep%posfreq==0)) WritePos(Nmole, Nspec);
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixReaxCSpecies::GatherBondOrder(struct _reax_list *lists)
-{
-  int i, ii, j, jj, nj, rj, inum, jnum;
-  int *ilist, *jlist, *numneigh, **firstneigh;
-  double bo_tmp;
-  bond_data *bo_ij;
-
-  inum = reaxc->list->inum;
-  ilist = reaxc->list->ilist;
-  numneigh = reaxc->list->numneigh;
-  firstneigh = reaxc->list->firstneigh;
-
-  for (ii = 0; ii < inum; ii++) {
-    i = ilist[ii];
-
-    tmpq[i] += atom->q[i]/repeat;
-    for (jj = 0; jj < 3; jj++)
-      tmpx[i][jj] += atom->x[i][jj]/repeat;
-    
-    jnum = numneigh[i];
-    jlist = firstneigh[i];
-    
-    for (jj = 0; jj < jnum; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
-
-      for( nj = Start_Index(i, reaxc->lists); nj < End_Index(i, reaxc->lists); ++nj ) {
-        bo_ij = &( reaxc->lists->select.bond_list[nj] );
-        rj = bo_ij->nbr;
-
-	if (atom->tag[j] == atom->tag[rj]) {
-          bo_tmp = bo_ij->bo_data.BO;
-	  abo[i][j] += bo_tmp/repeat;
-	}
-      }
-    }
+  if (posflag && ((ntimestep)%posfreq==0)) {
+    WritePos(Nmole, Nspec);
+    if (me == 0) fflush(pos);
   }
 
+  nvalid += nfreq;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixReaxCSpecies::FindMolecule()
+AtomCoord chAnchor(AtomCoord in1, AtomCoord in2)
 {
-  int i,j,ii,jj,inum,jnum,n,itype,jtype;
-  int change, done, anychange, loop, looptot;
+  if (in1.x < in2.x) 
+    return in1;
+  return in2;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpecies::FindMolecule ()
+{
+  int i,j,ii,jj,inum,jnum,n,itype,jtype,itag,jtag,loop,looptot;
+  int change,done,anychange;
   int *mask = atom->mask;
   int *ilist, *jlist, *numneigh, **firstneigh;
-  double bo_tmp, bo_cut;
+  double bo_tmp,bo_cut;
+  double **spec_atom = f_SPECBOND->array_atom;
 
   inum = reaxc->list->inum;
   ilist = reaxc->list->ilist;
@@ -392,7 +500,12 @@ void FixReaxCSpecies::FindMolecule()
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
-    if (mask[i] & groupbit) clusterID[i] = atom->tag[i];
+    if (mask[i] & groupbit) {
+      clusterID[i] = atom->tag[i];
+      x0[i].x = spec_atom[i][1];
+      x0[i].y = spec_atom[i][2];
+      x0[i].z = spec_atom[i][3];
+    }
     else clusterID[i] = 0.0;
   }
 
@@ -406,29 +519,35 @@ void FixReaxCSpecies::FindMolecule()
       done = 1;
 
       for (ii = 0; ii < inum; ii++) {
-	i = ilist[ii];
-	if (!(mask[i] & groupbit)) continue;
-	
-	itype = atom->type[i];
-	jlist = firstneigh[i];
-	jnum = numneigh[i];
+      	i = ilist[ii];
+      	if (!(mask[i] & groupbit)) continue;
 
-	for (jj = 0; jj < jnum; jj++) {
-	  j = jlist[jj];
+      	itype = atom->type[i];
 
-	  j &= NEIGHMASK;
-	  if (!(mask[j] & groupbit)) continue;
-	  if (clusterID[i] == clusterID[j]) continue;
+        for (jj = 0; jj < MAXSPECBOND; jj++) {
+      	  j = reaxc->tmpid[i][jj];
 
-	  jtype = atom->type[j];
-	  bo_cut = BOCut[itype][jtype];
-	  bo_tmp = abo[i][j];
+      	  if (j < i) continue;
+      	  if (!(mask[j] & groupbit)) continue;
 
-	  if (bo_tmp > bo_cut) {
-	    clusterID[i] = clusterID[j] = MIN(clusterID[i],clusterID[j]);
-	    done = 0;
-	  }
-	}
+      	  if (clusterID[i] == clusterID[j] && PBCconnected[i] == PBCconnected[j] 
+	    && x0[i].x == x0[j].x && x0[i].y == x0[j].y && x0[i].z == x0[j].z) continue;
+
+          jtype = atom->type[j];
+      	  bo_cut = BOCut[itype][jtype];
+      	  bo_tmp = spec_atom[i][jj+7];
+
+      	  if (bo_tmp > bo_cut) {
+            clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]);
+            PBCconnected[i] = PBCconnected[j] = MAX(PBCconnected[i], PBCconnected[j]);
+            x0[i] = x0[j] = chAnchor(x0[i], x0[j]);
+            if ((fabs(spec_atom[i][1] - spec_atom[j][1]) > reaxc->control->bond_cut)
+             || (fabs(spec_atom[i][2] - spec_atom[j][2]) > reaxc->control->bond_cut) 
+             || (fabs(spec_atom[i][3] - spec_atom[j][3]) > reaxc->control->bond_cut))
+              PBCconnected[i] = PBCconnected[j] = 1;
+      	    done = 0;
+      	  }
+      	}
       }
       if (!done) change = 1;
       if (done) break;
@@ -437,7 +556,8 @@ void FixReaxCSpecies::FindMolecule()
     if (!anychange) break;
 
     MPI_Allreduce(&loop,&looptot,1,MPI_INT,MPI_SUM,world);
-    if (looptot >= 200*nprocs) break;
+    if (looptot >= 400*nprocs) break;
+
   }
 }
 
@@ -500,12 +620,13 @@ void FixReaxCSpecies::SortMolecule(int &Nmole)
 
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
   if (flagall && comm->me == 0)
-    error->warning(FLERR,"One or more reax/c cluster has atoms not in group");
+    error->warning(FLERR,"One or more cluster has atoms not in group");
 
   for (n = 0; n < nlocal; n++) {
     if (!(mask[n] & groupbit)) continue;
-    clusterID[n] = molmap[nint(clusterID[n])-idlo]+1;
+    clusterID[n] = molmap[nint(clusterID[n])-idlo] + 1;
   }
+
   memory->destroy(molmap);
   molmap = NULL;
 
@@ -549,7 +670,7 @@ void FixReaxCSpecies::FindSpecies(int Nmole, int &Nspec)
     MPI_Allreduce(&flag_mol,&flag_tmp,1,MPI_INT,MPI_MAX,world);
     flag_mol = flag_tmp;
 
-    MPI_Reduce(Name,Nameall,ntypes,MPI_INT,MPI_SUM,0,world);
+    MPI_Allreduce(Name,Nameall,ntypes,MPI_INT,MPI_SUM,world);
     for (n = 0; n < ntypes; n++) Name[n] = Nameall[n];
 
     if (flag_mol == 1) {
@@ -603,7 +724,7 @@ int FixReaxCSpecies::CheckExistence(int id, int ntypes)
 
 /* ---------------------------------------------------------------------- */
 
-void FixReaxCSpecies::WriteFormulae(int Nmole, int Nspec)
+void FixReaxCSpecies::WriteFormulas(int Nmole, int Nspec)
 {
   int i, j, itemp;
   bigint ntimestep = update->ntimestep;
@@ -638,6 +759,7 @@ void FixReaxCSpecies::WriteFormulae(int Nmole, int Nspec)
 }
 
 /* ---------------------------------------------------------------------- */
+
 void FixReaxCSpecies::OpenPos()
 {
   char *filecurrent;
@@ -659,31 +781,35 @@ void FixReaxCSpecies::OpenPos()
 
   if (me == 0) {
     pos = fopen(filecurrent, "w");
-    if (pos == NULL) error->one(FLERR,"Cannot open fix species position file");
+    if (pos == NULL) error->one(FLERR,"Cannot open fix reax/c/species position file");
   } else pos = NULL;
   multipos_opened = 1;
 
-  delete [] filecurrent;
+  free(filecurrent);
 }
 
 /* ---------------------------------------------------------------------- */
+
 void FixReaxCSpecies::WritePos(int Nmole, int Nspec)
 {
-  int i,itype,cid;
+  int i, itype, cid;
   int count, count_tmp, m, n, k;
-  int *ilist, *jlist, *numneigh, **firstneigh;
-  double avq, avq_tmp, avx[3], avx_tmp, box[3];
-  int *mask =atom->mask;
   int *Nameall;
+  int *mask =atom->mask;
+  double avq, avq_tmp, avx[3], avx_tmp, box[3], halfbox[3];
+  double **spec_atom = f_SPECBOND->array_atom;
 
   if (multipos) OpenPos();
-
+  
   box[0] = domain->boxhi[0] - domain->boxlo[0];
   box[1] = domain->boxhi[1] - domain->boxlo[1];
   box[2] = domain->boxhi[2] - domain->boxlo[2];
 
+  for (int j = 0; j < 3; j++)
+    halfbox[j] = box[j] / 2;
+
   if (me == 0) {
-    fprintf(pos,"Timestep" BIGINT_FORMAT "NMole %d  NSpec %d  xlo %f  "
+    fprintf(pos,"Timestep "BIGINT_FORMAT " NMole %d  NSpec %d  xlo %f  "
 		"xhi %f  ylo %f  yhi %f  zlo %f  zhi %f\n",
 	        update->ntimestep,Nmole, Nspec,
 		domain->boxlo[0],domain->boxhi[0],
@@ -694,14 +820,16 @@ void FixReaxCSpecies::WritePos(int Nmole, int Nspec)
   }
 
   Nameall = NULL;
-  memory->create(Nameall,ntypes,"species:Nameall");
+  memory->create(Nameall,ntypes,"reax/c/species:Nameall");
 
   for (m = 1; m <= Nmole; m ++) {
 
     count = 0;
     avq = 0.0;
-    for (n = 0; n < 3; n++) avx[n] = 0.0;
-    for (n = 0; n < ntypes; n ++) Name[n] = 0;
+    for (n = 0; n < 3; n++) 
+      avx[n] = 0.0;
+    for (n = 0; n < ntypes; n ++) 
+      Name[n] = 0;
 
     for (i = 0; i < nlocal; i ++) {
       if (!(mask[i] & groupbit)) continue;
@@ -709,9 +837,24 @@ void FixReaxCSpecies::WritePos(int Nmole, int Nspec)
       if (cid == m) {
         itype = atom->type[i]-1;
         Name[itype] ++;
-	count ++;
-	avq += tmpq[i];
-	for (n = 0; n < 3; n++) avx[n] += tmpx[i][n];
+      	count ++;
+      	avq += spec_atom[i][0];
+        if (PBCconnected[i]) {
+          if ((x0[i].x - spec_atom[i][1]) > halfbox[0])
+            spec_atom[i][1] += box[0];
+          if ((spec_atom[i][1] - x0[i].x) > halfbox[0])
+             spec_atom[i][1] -= box[0];
+          if ((x0[i].y - spec_atom[i][2]) > halfbox[1])
+              spec_atom[i][2] += box[1];
+          if ((spec_atom[i][2] - x0[i].y) > halfbox[1])
+             spec_atom[i][2] -= box[1];
+          if ((x0[i].z - spec_atom[i][3]) > halfbox[2])
+            spec_atom[i][3] += box[2];
+          if ((spec_atom[i][3] - x0[i].z) > halfbox[2])
+           spec_atom[i][3] -= box[2];
+        }
+        for (n = 0; n < 3; n++)
+          avx[n] += spec_atom[i][n+1];
       }
     }
 
@@ -732,21 +875,26 @@ void FixReaxCSpecies::WritePos(int Nmole, int Nspec)
     for (n = 0; n < ntypes; n++) Name[n] = Nameall[n];
 
     if (me == 0) { 
-      fprintf(pos,"%d\t%d\t\t",m,count);
+      fprintf(pos,"%d\t%d\t",m,count);
       for (n = 0; n < ntypes; n++) { 
         if (Name[n] != 0) {
           if (eletype) fprintf(pos,"%s",eletype[n]);
           else fprintf(pos,"%c",ele[n]);
           if (Name[n] != 1) fprintf(pos,"%d",Name[n]);
-	}
+        }
       }
       if (count > 0) {
         avq /= count;
         for (k = 0; k < 3; k++) { 
-	  avx[k] /= count;
-	  avx[k] -= domain->boxlo[k];
-	  avx[k] /= box[k];
-	}
+      	  avx[k] /= count;
+          if (avx[k] >= domain->boxhi[k])
+            avx[k] -= box[k];
+          if (avx[k] < domain->boxlo[k])
+            avx[k] += box[k];
+
+      	  avx[k] -= domain->boxlo[k];
+      	  avx[k] /= box[k];
+        }
         fprintf(pos,"\t%.8f \t%.8f \t%.8f \t%.8f",
 		    avq,avx[0],avx[1],avx[2]);
       }
@@ -761,12 +909,12 @@ void FixReaxCSpecies::WritePos(int Nmole, int Nspec)
 
 double FixReaxCSpecies::compute_vector(int n)
 {
-  if (n == 0) {
+  if (n == 0) 
     return vector_nmole;
-  } else if (n == 1) {
+  if (n == 1) 
     return vector_nspec;
-  }
   return 0.0;
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -779,23 +927,6 @@ int FixReaxCSpecies::nint(const double &r)
   return i;
 }
 
-/* ----------------------------------------------------------------------
-   calculate nvalid = next step on which end_of_step does something
-   can be this timestep if multiple of nfreq and nrepeat = 1
-   else backup from next multiple of nfreq
-------------------------------------------------------------------------- */
-
-bigint FixReaxCSpecies::nextvalid()
-{
-  bigint nvalid = (update->ntimestep/nfreq)*nfreq + nfreq;
-  if (nvalid-nfreq == update->ntimestep && nrepeat == 1)
-    nvalid = update->ntimestep;
-  else
-    nvalid -= (nrepeat-1)*nevery;
-  if (nvalid < update->ntimestep) nvalid += nfreq;
-  return nvalid-0;
-}
-
 /* ---------------------------------------------------------------------- */
 
 int FixReaxCSpecies::pack_comm(int n, int *list, double *buf, 
@@ -806,9 +937,14 @@ int FixReaxCSpecies::pack_comm(int n, int *list, double *buf,
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    buf[m++] = clusterID[j];
+    buf[m] = clusterID[j];
+    buf[m+1] = (double)PBCconnected[j];
+    buf[m+2] = x0[j].x;
+    buf[m+3] = x0[j].y;
+    buf[m+4] = x0[j].z;
+    m += 5;
   }
-  return 1;
+  return 5;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -819,7 +955,14 @@ void FixReaxCSpecies::unpack_comm(int n, int first, double *buf)
 
   m = 0;
   last = first + n;
-  for (i = first; i < last; i++) clusterID[i] = buf[m++];
+  for (i = first; i < last; i++) {
+    clusterID[i] = buf[m];
+    PBCconnected[i] = (int)buf[m+1];
+    x0[i].x = buf[m+2];
+    x0[i].y = buf[m+3];
+    x0[i].z = buf[m+4];
+    m += 5;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -828,8 +971,9 @@ double FixReaxCSpecies::memory_usage()
 {
   double bytes;
 
-  bytes = 2.0*nmax*sizeof(double);
-  bytes += nmax*nmax*sizeof(double);
+  bytes = 5*nmax*sizeof(double);  // clusterID + PBCconnected + x0
 
   return bytes;
 }
+
+/* ---------------------------------------------------------------------- */
diff --git a/src/USER-REAXC/fix_reaxc_species.h b/src/USER-REAXC/fix_reaxc_species.h
index 517e3f7edb..12c0aa98c5 100644
--- a/src/USER-REAXC/fix_reaxc_species.h
+++ b/src/USER-REAXC/fix_reaxc_species.h
@@ -20,12 +20,21 @@ FixStyle(reax/c/species,FixReaxCSpecies)
 #ifndef LMP_FIX_REAXC_SPECIES_H
 #define LMP_FIX_REAXC_SPECIES_H
 
-#include "stdio.h"
 #include "fix.h"
 #include "pointers.h"
 
+#include "pair_reax_c.h"
+#include "reaxc_types.h"
+#include "reaxc_defs.h"
+
+#define BUFLEN 1000
+
 namespace LAMMPS_NS {
 
+typedef struct {
+  double x, y, z;
+} AtomCoord;
+
 class FixReaxCSpecies : public Fix {
  public:
   FixReaxCSpecies(class LAMMPS *, int, char **);
@@ -38,42 +47,45 @@ class FixReaxCSpecies : public Fix {
   double compute_vector(int);
 
  private:
-  class PairReaxC *reaxc;
-  class NeighList *list;
-
   int me, nprocs, nmax, nlocal, ntypes, ntotal;
-  int nrepeat, irepeat, repeat, nfreq, posfreq;
+  bigint natoms;
+  int nrepeat, nfreq, posfreq;
   int Nmoltype, vector_nmole, vector_nspec;
   int *Name, *MolName, *NMol, *nd, *MolType, *molmap;
   double *clusterID;
+  int *PBCconnected;
+  AtomCoord *x0;
 
   double bg_cut;
-  double *tmpq, **tmpx;
-  double **BOCut, **abo;
+  double **BOCut;
+  char **tmparg;
 
   FILE *fp, *pos;
-  int eleflag, posflag, padflag, multipos;
+  int eleflag, posflag, multipos, padflag, setupflag;
   int singlepos_opened, multipos_opened;
-  char *ele, **eletype, *posspec, *filepos;
+  char *ele, **eletype, *filepos;
 
   void Output_ReaxC_Bonds(bigint, FILE *);
-  void GatherBondOrder(struct _reax_list*);
+  void create_compute();
+  void create_fix();
   void FindMolecule();
   void SortMolecule(int &);
   void FindSpecies(int, int &);
-  void WriteFormulae(int, int);
-  void OpenPos();
-  void WritePos(int, int);
-
+  void WriteFormulas(int, int);
   int CheckExistence(int, int);
-  int nint(const double &);
 
+  int nint(const double &);
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
+  void OpenPos();
+  void WritePos(int, int);
   double memory_usage();
 
-  bigint nvalid, nextvalid();
-  struct _reax_list *lists;
+  bigint nvalid;
+
+  class NeighList *list;
+  class FixAveAtom *f_SPECBOND;
+  class PairReaxC *reaxc;
 
 };
 }
diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index dd7c0e558e..57212a0c4d 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -108,11 +108,14 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
   system->pair_ptr = this;
 
   fix_reax = NULL;
+  tmpid = NULL;
+  tmpbo = NULL;
 
   nextra = 14;
   pvector = new double[nextra];
 
   setup_flag = 0;
+  fixspecies_flag = 0;
 
   nmax = 0;
 }
@@ -159,6 +162,9 @@ PairReaxC::~PairReaxC()
     delete [] gamma;
   }
 
+  memory->destroy(tmpid);
+  memory->destroy(tmpbo);
+
   delete [] pvector;
 
 }
@@ -534,6 +540,26 @@ void PairReaxC::compute(int eflag, int vflag)
 
   Output_Results( system, control, data, &lists, out_control, mpi_data );
 
+  // populate tmpid and tmpbo arrays for fix reax/c/species
+  int i, j;
+
+  if(fixspecies_flag) {
+    if (system->N > nmax) {
+      memory->destroy(tmpid);
+      memory->destroy(tmpbo);
+      nmax = system->N;
+      memory->create(tmpid,nmax,MAXSPECBOND,"pair:tmpid");
+      memory->create(tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");
+    }
+   
+    for (i = 0; i < system->N; i ++)
+      for (j = 0; j < MAXSPECBOND; j ++) {
+        tmpbo[i][j] = 0.0;
+	tmpid[i][j] = 0;
+      }
+    FindBond();
+  }
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -736,3 +762,65 @@ void *PairReaxC::extract(const char *str, int &dim)
 }
 
 /* ---------------------------------------------------------------------- */
+
+double PairReaxC::memory_usage()
+{
+  double bytes = 0.0;
+  
+  // From pair_reax_c
+  bytes += 1.0 * system->N * sizeof(int);
+  bytes += 1.0 * system->N * sizeof(double);
+
+  // From reaxc_allocate: BO
+  bytes += 1.0 * system->total_cap * sizeof(reax_atom);
+  bytes += 19.0 * system->total_cap * sizeof(real);
+  bytes += 3.0 * system->total_cap * sizeof(int);
+
+  double mem1 = bytes;
+
+  // From reaxc_lists
+  bytes += 2.0 * lists->n * sizeof(int);
+  bytes += lists->num_intrs * sizeof(three_body_interaction_data);
+  bytes += lists->num_intrs * sizeof(bond_data);
+  bytes += lists->num_intrs * sizeof(dbond_data);
+  bytes += lists->num_intrs * sizeof(dDelta_data);
+  bytes += lists->num_intrs * sizeof(far_neighbor_data);
+  bytes += lists->num_intrs * sizeof(hbond_data);
+
+  if(fixspecies_flag)
+    bytes += 2 * nmax * MAXSPECBOND * sizeof(double);
+
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairReaxC::FindBond()
+{
+  int i, ii, j, pj, jtag, nj, jtmp, jj;
+  double bo_tmp, bo_cut, rij, rsq, r_tmp;
+
+  bond_data *bo_ij;
+  bo_cut = 0.10;
+
+  for (i = 0; i < system->n; i++) {
+    nj = 0;
+    for( pj = Start_Index(i, lists); pj < End_Index(i, lists); ++pj ) {
+      bo_ij = &( lists->select.bond_list[pj] );
+      j = bo_ij->nbr;
+      if (j < i) continue;
+
+      bo_tmp = bo_ij->bo_data.BO;
+      r_tmp = bo_ij->d;
+
+      if (bo_tmp >= bo_cut ) {
+	tmpid[i][nj] = j;
+	tmpbo[i][nj] = bo_tmp;
+	nj ++;
+	if (nj > MAXSPECBOND) error->all(FLERR,"Increase MAXSPECBOND in fix_reaxc_species.h");
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
diff --git a/src/USER-REAXC/pair_reax_c.h b/src/USER-REAXC/pair_reax_c.h
index 060ef2ce8e..5f89a95d87 100644
--- a/src/USER-REAXC/pair_reax_c.h
+++ b/src/USER-REAXC/pair_reax_c.h
@@ -78,6 +78,7 @@ class PairReaxC : public Pair {
 
   int nmax;
   void FindBond();
+  double memory_usage();
 };
 
 }