diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
index 5da0f33e3a..af443ede92 100644
--- a/doc/src/fix_bond_react.txt
+++ b/doc/src/fix_bond_react.txt
@@ -385,6 +385,10 @@ No parameter of this fix can be used with the {start/stop} keywords of
 the "run"_run.html command.  This fix is not invoked during "energy
 minimization"_minimize.html.
 
+When fix bond/react is 'unfixed,' all internally-created groups are
+deleted. Therefore, fix bond/react can only be unfixed after unfixing
+all other fixes that use any group created by fix bond/react.
+
 [Restrictions:]
 
 This fix is part of the USER-MISC package.  It is only enabled if
diff --git a/src/MOLECULE/dihedral_helix.cpp b/src/MOLECULE/dihedral_helix.cpp
index ae23b77951..0998d654af 100644
--- a/src/MOLECULE/dihedral_helix.cpp
+++ b/src/MOLECULE/dihedral_helix.cpp
@@ -39,7 +39,10 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralHelix::DihedralHelix(LAMMPS *lmp) : Dihedral(lmp) {}
+DihedralHelix::DihedralHelix(LAMMPS *lmp) : Dihedral(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/MOLECULE/dihedral_multi_harmonic.cpp b/src/MOLECULE/dihedral_multi_harmonic.cpp
index f6461abb6e..09c90f8d81 100644
--- a/src/MOLECULE/dihedral_multi_harmonic.cpp
+++ b/src/MOLECULE/dihedral_multi_harmonic.cpp
@@ -34,7 +34,10 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *lmp) : Dihedral(lmp) {}
+DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *lmp) : Dihedral(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/USER-MISC/dihedral_cosine_shift_exp.cpp b/src/USER-MISC/dihedral_cosine_shift_exp.cpp
index 82da173f8e..6a04228226 100644
--- a/src/USER-MISC/dihedral_cosine_shift_exp.cpp
+++ b/src/USER-MISC/dihedral_cosine_shift_exp.cpp
@@ -36,7 +36,10 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralCosineShiftExp::DihedralCosineShiftExp(LAMMPS *lmp) : Dihedral(lmp) {}
+DihedralCosineShiftExp::DihedralCosineShiftExp(LAMMPS *lmp) : Dihedral(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/USER-MISC/dihedral_nharmonic.cpp b/src/USER-MISC/dihedral_nharmonic.cpp
index f8e8850680..aa3f6a6f3e 100644
--- a/src/USER-MISC/dihedral_nharmonic.cpp
+++ b/src/USER-MISC/dihedral_nharmonic.cpp
@@ -35,7 +35,8 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp) {
+DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp)
+{
   writedata = 1;
 }
 
diff --git a/src/USER-MISC/dihedral_quadratic.cpp b/src/USER-MISC/dihedral_quadratic.cpp
index 1b64b52faf..c9bce4f02c 100644
--- a/src/USER-MISC/dihedral_quadratic.cpp
+++ b/src/USER-MISC/dihedral_quadratic.cpp
@@ -38,7 +38,10 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralQuadratic::DihedralQuadratic(LAMMPS *lmp) : Dihedral(lmp) {}
+DihedralQuadratic::DihedralQuadratic(LAMMPS *lmp) : Dihedral(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/USER-MISC/dihedral_spherical.cpp b/src/USER-MISC/dihedral_spherical.cpp
index 77fa885b7a..8ce58243ac 100644
--- a/src/USER-MISC/dihedral_spherical.cpp
+++ b/src/USER-MISC/dihedral_spherical.cpp
@@ -41,7 +41,10 @@ using namespace MathExtra;
 
 /* ---------------------------------------------------------------------- */
 
-DihedralSpherical::DihedralSpherical(LAMMPS *lmp) : Dihedral(lmp) {}
+DihedralSpherical::DihedralSpherical(LAMMPS *lmp) : Dihedral(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -817,10 +820,11 @@ void DihedralSpherical::write_data(FILE *fp)
   for (int i = 1; i <= atom->ndihedraltypes; i++) {
     fprintf(fp,"%d %d ", i , nterms[i]);
     for (int j = 0; j < nterms[i]; j++) {
-      fprintf(fp, "%g %g %g %g %g %g %g %g %g ",
-              phi_mult[i][j],    phi_shift[i][j],    phi_offset[i][j],
-              theta1_mult[i][j], theta1_shift[i][j], theta1_offset[i][j],
-              theta2_mult[i][j], theta2_shift[i][j], theta2_offset[i][j]);
+      fprintf(fp, "%g %g %g %g %g %g %g %g %g %g ", Ccoeff[i][j],
+              phi_mult[i][j], phi_shift[i][j]*180.0/MY_PI, phi_offset[i][j],
+              theta1_mult[i][j], theta1_shift[i][j]*180.0/MY_PI,
+              theta1_offset[i][j], theta2_mult[i][j],
+              theta2_shift[i][j]*180.0/MY_PI, theta2_offset[i][j]);
     }
     fprintf(fp,"\n");
   }
diff --git a/src/USER-MISC/dihedral_table_cut.cpp b/src/USER-MISC/dihedral_table_cut.cpp
index 6ebe094e50..e9596df6c8 100644
--- a/src/USER-MISC/dihedral_table_cut.cpp
+++ b/src/USER-MISC/dihedral_table_cut.cpp
@@ -425,17 +425,6 @@ DihedralTableCut::DihedralTableCut(LAMMPS *lmp) : Dihedral(lmp)
 DihedralTableCut::~DihedralTableCut()
 {
   if (allocated) {
-    memory->destroy(setflag);
-    memory->destroy(setflag_d);
-    memory->destroy(setflag_aat);
-
-    memory->destroy(k1);
-    memory->destroy(k2);
-    memory->destroy(k3);
-    memory->destroy(phi1);
-    memory->destroy(phi2);
-    memory->destroy(phi3);
-
     memory->destroy(aat_k);
     memory->destroy(aat_theta0_1);
     memory->destroy(aat_theta0_2);
@@ -761,27 +750,15 @@ void DihedralTableCut::allocate()
   allocated = 1;
   int n = atom->ndihedraltypes;
 
-  memory->create(k1,n+1,"dihedral:k1");
-  memory->create(k2,n+1,"dihedral:k2");
-  memory->create(k3,n+1,"dihedral:k3");
-  memory->create(phi1,n+1,"dihedral:phi1");
-  memory->create(phi2,n+1,"dihedral:phi2");
-  memory->create(phi3,n+1,"dihedral:phi3");
-
   memory->create(aat_k,n+1,"dihedral:aat_k");
   memory->create(aat_theta0_1,n+1,"dihedral:aat_theta0_1");
   memory->create(aat_theta0_2,n+1,"dihedral:aat_theta0_2");
 
-  memory->create(setflag,n+1,"dihedral:setflag");
-  memory->create(setflag_d,n+1,"dihedral:setflag_d");
-  memory->create(setflag_aat,n+1,"dihedral:setflag_aat");
-
   memory->create(tabindex,n+1,"dihedral:tabindex");
-  //memory->create(phi0,n+1,"dihedral:phi0"); <-equilibrium angles not supported
   memory->create(setflag,n+1,"dihedral:setflag");
 
   for (int i = 1; i <= n; i++)
-    setflag[i] = setflag_d[i] = setflag_aat[i] = 0;
+    setflag[i] = 0;
 }
 
 void DihedralTableCut::settings(int narg, char **arg)
@@ -824,9 +801,6 @@ void DihedralTableCut::coeff(int narg, char **arg)
   int ilo,ihi;
   force->bounds(FLERR,arg[0],atom->ndihedraltypes,ilo,ihi);
 
-  int count = 0;
-
-
   double k_one = force->numeric(FLERR,arg[2]);
   double theta0_1_one = force->numeric(FLERR,arg[3]);
   double theta0_2_one = force->numeric(FLERR,arg[4]);
@@ -837,8 +811,6 @@ void DihedralTableCut::coeff(int narg, char **arg)
     aat_k[i] = k_one;
     aat_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
     aat_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
-    setflag_aat[i] = 1;
-    count++;
   }
 
   int me;
@@ -998,8 +970,7 @@ void DihedralTableCut::coeff(int narg, char **arg)
             //  To be nice and report something, I do the same thing here.)
             cyc_splintD(tb->phi, tb->e, tb->e2, tablength, MY_2PI,phi);
           f = -dU_dphi;
-        }
-        else
+        } else
           // Otherwise we calculated the tb->f[] array.  Report its contents.
           f = tb->f[i];
         if (tb->use_degrees) {
@@ -1015,9 +986,8 @@ void DihedralTableCut::coeff(int narg, char **arg)
   } // if (me == 0)
 
   // store ptr to table in tabindex
-  count = 0;
-  for (int i = ilo; i <= ihi; i++)
-  {
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
     tabindex[i] = ntables;
     //phi0[i] = tb->phi0; <- equilibrium dihedral angles not supported
     setflag[i] = 1;
@@ -1025,12 +995,7 @@ void DihedralTableCut::coeff(int narg, char **arg)
   }
   ntables++;
 
-
   if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
-
-  for (int i = ilo; i <= ihi; i++)
-    if (setflag_d[i] == 1 && setflag_aat[i] == 1 )
-      setflag[i] = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -1041,16 +1006,6 @@ void DihedralTableCut::write_restart(FILE *fp)
 {
   fwrite(&tabstyle,sizeof(int),1,fp);
   fwrite(&tablength,sizeof(int),1,fp);
-  fwrite(&k1[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&k2[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&k3[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&phi1[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&phi2[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&phi3[1],sizeof(double),atom->ndihedraltypes,fp);
-
-  fwrite(&aat_k[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);
-  fwrite(&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -1064,33 +1019,10 @@ void DihedralTableCut::read_restart(FILE *fp)
   if (comm->me == 0) {
     fread(&tabstyle,sizeof(int),1,fp);
     fread(&tablength,sizeof(int),1,fp);
-    fread(&k1[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&k2[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&k3[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&phi1[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&phi2[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&phi3[1],sizeof(double),atom->ndihedraltypes,fp);
-
-    fread(&aat_k[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);
-    fread(&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);
   }
 
-  MPI_Bcast(&k1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&k2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&k3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&phi1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&phi2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&phi3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-
-  MPI_Bcast(&aat_k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&aat_theta0_1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
-  MPI_Bcast(&aat_theta0_2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
   MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
   MPI_Bcast(&tablength,1,MPI_INT,0,world);
-
-  allocate();
-  for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1474,22 +1406,3 @@ void DihedralTableCut::bcast_table(Table *tb)
 }
 
 
-
-/* ----------------------------------------------------------------------
-   proc 0 writes to data file
-------------------------------------------------------------------------- */
-
-void DihedralTableCut::write_data(FILE *fp)
-{
-  for (int i = 1; i <= atom->ndihedraltypes; i++)
-    fprintf(fp,"%d %g %g %g %g %g %g\n",i,
-            k1[i],phi1[i]*180.0/MY_PI,
-            k2[i],phi2[i]*180.0/MY_PI,
-            k3[i],phi3[i]*180.0/MY_PI);
-
-  fprintf(fp,"\nAngleAngleTorsion Coeffs\n\n");
-  for (int i = 1; i <= atom->ndihedraltypes; i++)
-    fprintf(fp,"%d %g %g %g\n",i,aat_k[i],
-            aat_theta0_1[i]*180.0/MY_PI,aat_theta0_2[i]*180.0/MY_PI);
-
-}
diff --git a/src/USER-MISC/dihedral_table_cut.h b/src/USER-MISC/dihedral_table_cut.h
index d01903c88b..dd645bedda 100644
--- a/src/USER-MISC/dihedral_table_cut.h
+++ b/src/USER-MISC/dihedral_table_cut.h
@@ -20,7 +20,6 @@ DihedralStyle(table/cut,DihedralTableCut)
 #ifndef LMP_DIHEDRAL_TABLE_CUT_H
 #define LMP_DIHEDRAL_TABLE_CUT_H
 
-#include <stdio.h>
 #include "dihedral.h"
 
 namespace LAMMPS_NS {
@@ -34,26 +33,18 @@ class DihedralTableCut : public Dihedral {
   void coeff(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void write_data(FILE *);
-  double single(int type, int i1, int i2, int i3, int i4);
 
  protected:
-  double *k1,*k2,*k3;
-  double *phi1,*phi2,*phi3;
   double *aat_k,*aat_theta0_1,*aat_theta0_2;
-  int *setflag_d;
-  int *setflag_aat;
 
   void allocate();
 
   int tabstyle,tablength;
-  // double *phi0;       <- equilibrium angles not supported
   char *checkU_fname;
   char *checkF_fname;
 
   struct Table {
     int ninput;
-    //double phi0;      <-equilibrium angles not supported
     int f_unspecified; // boolean (but MPI does not like type "bool")
     int use_degrees;   // boolean (but MPI does not like type "bool")
     double *phifile,*efile,*ffile;
diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index c5df3ed438..5d624ad7f4 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -457,8 +457,6 @@ FixBondReact::~FixBondReact()
   memory->destroy(global_mega_glove);
 
   if (stabilization_flag == 1) {
-    delete [] exclude_group;
-
     // check nfix in case all fixes have already been deleted
     if (id_fix1 && modify->nfix) modify->delete_fix(id_fix1);
     delete [] id_fix1;
@@ -473,6 +471,18 @@ FixBondReact::~FixBondReact()
   delete [] statted_id;
   delete [] guess_branch;
   delete [] pioneer_count;
+
+  char **newarg;
+  newarg = new char*[2];
+  newarg[0] = master_group;
+  newarg[1] = (char *) "delete";
+  group->assign(2,newarg);
+  if (stabilization_flag == 1) {
+    newarg[0] = exclude_group;
+    group->assign(2,newarg);
+    delete [] exclude_group;
+  }
+  delete [] newarg;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-PHONON/dynamical_matrix.cpp b/src/USER-PHONON/dynamical_matrix.cpp
index 372f4e4e31..d94bd11a80 100644
--- a/src/USER-PHONON/dynamical_matrix.cpp
+++ b/src/USER-PHONON/dynamical_matrix.cpp
@@ -138,6 +138,9 @@ void DynamicalMatrix::command(int narg, char **arg)
     else if (style == ESKM) options(narg-3,&arg[3]); //COME BACK
     else if (comm->me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
 
+    if (atom->map_style == 0)
+      error->all(FLERR,"Dynamical_matrix command requires an atom map, see atom_modify");
+
     // move atoms by 3-vector or specified variable(s)
 
     if (style == REGULAR) {
@@ -240,7 +243,7 @@ void DynamicalMatrix::calculateMatrix()
     int nlocal = atom->nlocal;
     bigint natoms = atom->natoms;
     int *type = atom->type;
-    int *gm = groupmap;
+    bigint *gm = groupmap;
     double imass; // dynamical matrix element
     double *m = atom->mass;
     double **f = atom->f;
@@ -256,20 +259,30 @@ void DynamicalMatrix::calculateMatrix()
     //initialize dynmat to all zeros
     dynmat_clear(dynmat);
 
-    if (comm->me == 0 && screen) fprintf(screen,"Calculating Dynamical Matrix...\n");
+    if (comm->me == 0 && screen) {
+        fprintf(screen,"Calculating Dynamical Matrix ...\n");
+        fprintf(screen,"  Total # of atoms = " BIGINT_FORMAT "\n", natoms);
+        fprintf(screen,"  Atoms in group = " BIGINT_FORMAT "\n", gcount);
+        fprintf(screen,"  Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
+    }
+    
+    // emit dynlen rows of dimalpha*dynlen*dimbeta elements
 
     update->nsteps = 0;
+    int prog = 0;
     for (bigint i=1; i<=natoms; i++){
         local_idx = atom->map(i);
+        if (gm[i-1] < 0)
+            continue;
         for (bigint alpha=0; alpha<3; alpha++){
             displace_atom(local_idx, alpha, 1);
             update_force();
             for (bigint j=1; j<=natoms; j++){
                 local_jdx = atom->map(j);
                 if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
-                    && gm[i-1] >= 0 && gm[j-1] >= 0){
+                    && gm[j-1] >= 0){
                     for (int beta=0; beta<3; beta++){
-                        dynmat[alpha][(gm[j-1])*3+beta] -= f[local_jdx][beta];
+                        dynmat[alpha][gm[j-1]*3+beta] -= f[local_jdx][beta];
                     }
                 }
             }
@@ -278,15 +291,15 @@ void DynamicalMatrix::calculateMatrix()
             for (bigint j=1; j<=natoms; j++){
                 local_jdx = atom->map(j);
                 if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
-                    && gm[i-1] >= 0 && gm[j-1] >= 0){
+                    && gm[j-1] >= 0){
                     for (bigint beta=0; beta<3; beta++){
                         if (atom->rmass_flag == 1)
                             imass = sqrt(m[local_idx] * m[local_jdx]);
                         else
                             imass = sqrt(m[type[local_idx]] * m[type[local_jdx]]);
-                        dynmat[alpha][(gm[j-1])*3+beta] -= -f[local_jdx][beta];
-                        dynmat[alpha][(gm[j-1])*3+beta] /= (2 * del * imass);
-                        dynmat[alpha][(gm[j-1])*3+beta] *= conversion;
+                        dynmat[alpha][gm[j-1]*3+beta] -= -f[local_jdx][beta];
+                        dynmat[alpha][gm[j-1]*3+beta] /= (2 * del * imass);
+                        dynmat[alpha][gm[j-1]*3+beta] *= conversion;
                     }
                 }
             }
@@ -297,7 +310,16 @@ void DynamicalMatrix::calculateMatrix()
         if (me == 0)
             writeMatrix(fdynmat);
         dynmat_clear(dynmat);
+        if (comm->me == 0 && screen) {
+            int p = 10 * gm[i-1] / gcount;
+            if (p > prog) {
+              prog = p;
+              fprintf(screen," %d%%",p*10);
+              fflush(screen);
+            }
+        }
     }
+    if (comm->me == 0 && screen) fprintf(screen,"\n");
 
     for (int i=0; i < 3; i++)
         delete [] dynmat[i];
@@ -316,27 +338,24 @@ void DynamicalMatrix::calculateMatrix()
 
 void DynamicalMatrix::writeMatrix(double **dynmat)
 {
-    if (me != 0 || fp == NULL) return;
+    if (me != 0 || !fp)
+        return;
 
-    // print file comment lines
-
-    if (!binaryflag && fp) {
-        clearerr(fp);
+    clearerr(fp);
+    if (binaryflag) {
+        for (int i=0; i<3; i++)
+            fwrite(dynmat[i], sizeof(double), dynlen, fp);
+        if (ferror(fp))
+            error->one(FLERR, "Error writing to binary file");
+    } else {
         for (int i = 0; i < 3; i++) {
             for (bigint j = 0; j < dynlen; j++) {
                 if ((j+1)%3==0) fprintf(fp, "%4.8f\n", dynmat[i][j]);
                 else fprintf(fp, "%4.8f ",dynmat[i][j]);
             }
         }
-    }
-    if (ferror(fp))
-        error->one(FLERR,"Error writing to file");
-
-    if (binaryflag && fp) {
-        clearerr(fp);
-        fwrite(&dynmat[0], sizeof(double), 3 * dynlen, fp);
         if (ferror(fp))
-            error->one(FLERR, "Error writing to binary file");
+            error->one(FLERR,"Error writing to file");
     }
 }
 
@@ -490,6 +509,7 @@ void DynamicalMatrix::convert_units(const char *style)
 void DynamicalMatrix::create_groupmap()
 {
     //Create a group map which maps atom order onto group
+    // groupmap[global atom index-1] = output column/row
 
     int local_idx; // local index
     int gid = 0; //group index
@@ -498,7 +518,7 @@ void DynamicalMatrix::create_groupmap()
     bigint natoms = atom->natoms;
     int *recv = new int[comm->nprocs];
     int *displs = new int[comm->nprocs];
-    int *temp_groupmap = new int[natoms];
+    bigint *temp_groupmap = new bigint[natoms];
 
     //find number of local atoms in the group (final_gid)
     for (bigint i=1; i<=natoms; i++){
@@ -507,7 +527,7 @@ void DynamicalMatrix::create_groupmap()
             gid += 1; // gid at the end of loop is final_Gid
     }
     //create an array of length final_gid
-    int *sub_groupmap = new int[gid];
+    bigint *sub_groupmap = new bigint[gid];
 
     gid = 0;
     //create a map between global atom id and group atom id for each proc
@@ -524,7 +544,7 @@ void DynamicalMatrix::create_groupmap()
         recv[i] = 0;
     }
     recv[comm->me] = gid;
-    MPI_Allreduce(recv,displs,4,MPI_INT,MPI_SUM,world);
+    MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
     for (int i=0; i<comm->nprocs; i++){
         recv[i]=displs[i];
         if (i>0) displs[i] = displs[i-1]+recv[i-1];
@@ -532,15 +552,17 @@ void DynamicalMatrix::create_groupmap()
     }
 
     //combine subgroup maps into total temporary groupmap
-    MPI_Allgatherv(sub_groupmap,gid,MPI_INT,temp_groupmap,recv,displs,MPI_INT,world);
+    MPI_Allgatherv(sub_groupmap,gid,MPI_LMP_BIGINT,temp_groupmap,recv,displs,MPI_LMP_BIGINT,world);
     std::sort(temp_groupmap,temp_groupmap+gcount);
 
     //populate member groupmap based on temp groupmap
-    for (bigint i=0; i<natoms; i++){
-        if (i==temp_groupmap[i]-1)
-            groupmap[i] = temp_groupmap[i]-1;
+    bigint j = 0;
+    for (bigint i=1; i<=natoms; i++){
+        // flag groupmap contents that are in temp_groupmap
+        if (j < gcount && i == temp_groupmap[j])
+            groupmap[i-1] = j++;
         else
-            groupmap[i] = -1;
+            groupmap[i-1] = -1;
     }
 
     //free that memory!
diff --git a/src/USER-PHONON/dynamical_matrix.h b/src/USER-PHONON/dynamical_matrix.h
index 56de814e97..8ff11044ea 100644
--- a/src/USER-PHONON/dynamical_matrix.h
+++ b/src/USER-PHONON/dynamical_matrix.h
@@ -58,7 +58,7 @@ namespace LAMMPS_NS {
         bigint dynlen;             // rank of dynamical matrix
         int scaleflag;
         int me;
-        int *groupmap;
+        bigint *groupmap;
 
         int compressed;            // 1 if dump file is written compressed, 0 no
         int binaryflag;            // 1 if dump file is written binary, 0 no