diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 9a439d87d6b987a3f9db5e4b93a1ded74b8d7338..4dc0c5cbefa6e32d0f4ff971fe9c4a0db86eb5c9 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -666,7 +666,7 @@ void FixPour::pre_exchange()
 void FixPour::find_maxid()
 {
   tagint *tag = atom->tag;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
 
   tagint max = 0;
@@ -676,7 +676,7 @@ void FixPour::find_maxid()
   if (mode == MOLECULE && molecule) {
     max = 0;
     for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
-    MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world);
+    MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
   }
 }
 
diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h
index 34d13da9855fb3aa68d5e1176b6004aa38baec15..e57379906e0a00ea6f5b96de2a779dd0bb9b14c9 100644
--- a/src/GRANULAR/fix_pour.h
+++ b/src/GRANULAR/fix_pour.h
@@ -62,8 +62,7 @@ class FixPour : public Fix {
   int *recvcounts,*displs;
   int nfreq,nfirst,ninserted,nper;
   double lo_current,hi_current;
-  tagint maxtag_all;
-  int maxmol_all;
+  tagint maxtag_all,maxmol_all;
   class RanPark *random,*random2;
 
   void find_maxid();
diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp
index 31bdb75682a3ea094373357a06010fefe286cb7c..79b298391ae9a512dc90db9483653de765ca59c8 100644
--- a/src/MC/fix_bond_swap.cpp
+++ b/src/MC/fix_bond_swap.cpp
@@ -126,7 +126,8 @@ void FixBondSwap::init()
   // require an atom style with molecule IDs
 
   if (atom->molecule == NULL)
-    error->all(FLERR,"Must use atom style with molecule IDs with fix bond/swap");
+    error->all(FLERR,
+               "Must use atom style with molecule IDs with fix bond/swap");
 
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0)
@@ -193,7 +194,7 @@ void FixBondSwap::pre_neighbor()
 
   tagint *tag = atom->tag;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *num_bond = atom->num_bond;
   tagint **bond_atom = atom->bond_atom;
   int **bond_type = atom->bond_type;
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index ab85e6e21b619ee800480fcbf9c1d45183a873f4..90202ad7eb546fde02c92858a316f4e9d1fc5da2 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -252,7 +252,7 @@ void FixGCMC::init()
   // if molflag not set, warn if any deletable atom has a mol ID
 
   if (molflag == 0 && atom->molecule_flag) {
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int *mask = atom->mask;
     int flag = 0;
     for (int i = 0; i < atom->nlocal; i++)
@@ -268,7 +268,7 @@ void FixGCMC::init()
   // if molflag set, check for unset mol IDs
 
   if (molflag == 1) {
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int *mask = atom->mask;
     int flag = 0;
     for (int i = 0; i < atom->nlocal; i++)
@@ -578,7 +578,7 @@ void FixGCMC::attempt_molecule_translation()
 
   if (ngas == 0) return;
 
-  int translation_molecule = pick_random_gas_molecule();
+  tagint translation_molecule = pick_random_gas_molecule();
   if (translation_molecule == -1) return;
 
   double energy_before_sum = molecule_energy(translation_molecule);
@@ -637,7 +637,7 @@ void FixGCMC::attempt_molecule_rotation()
 
   if (ngas == 0) return;
 
-  int rotation_molecule = pick_random_gas_molecule();
+  tagint rotation_molecule = pick_random_gas_molecule();
   if (rotation_molecule == -1) return;
   
   double energy_before_sum = molecule_energy(rotation_molecule);
@@ -719,7 +719,7 @@ void FixGCMC::attempt_molecule_deletion()
 
   if (ngas == 0) return;
   
-  int deletion_molecule = pick_random_gas_molecule();
+  tagint deletion_molecule = pick_random_gas_molecule();
   if (deletion_molecule == -1) return;
 
   double deletion_energy_sum = molecule_energy(deletion_molecule);
@@ -756,13 +756,20 @@ void FixGCMC::attempt_molecule_insertion()
   double com_coord[3];
   if (regionflag) {
     int region_attempt = 0;
-    com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
-    com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
-    com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
-    while (domain->regions[iregion]->match(com_coord[0],com_coord[1],com_coord[2]) == 0) {
-      com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
-      com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
-      com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
+    com_coord[0] = region_xlo + random_equal->uniform() * 
+      (region_xhi-region_xlo);
+    com_coord[1] = region_ylo + random_equal->uniform() * 
+      (region_yhi-region_ylo);
+    com_coord[2] = region_zlo + random_equal->uniform() * 
+      (region_zhi-region_zlo);
+    while (domain->regions[iregion]->match(com_coord[0],com_coord[1],
+                                           com_coord[2]) == 0) {
+      com_coord[0] = region_xlo + random_equal->uniform() * 
+        (region_xhi-region_xlo);
+      com_coord[1] = region_ylo + random_equal->uniform() * 
+        (region_yhi-region_ylo);
+      com_coord[2] = region_zlo + random_equal->uniform() * 
+        (region_zhi-region_zlo);
       region_attempt++;
       if (region_attempt >= max_region_attempts) return;
     }
@@ -779,9 +786,12 @@ void FixGCMC::attempt_molecule_insertion()
   double insertion_energy = 0.0;
   bool procflag[natoms_per_molecule];
   for (int i = 0; i < natoms_per_molecule; i++) {
-    atom_coord[i][0] = rot[0]*model_x[i][0] + rot[1]*model_x[i][1] + rot[2]*model_x[i][2] + com_coord[0];
-    atom_coord[i][1] = rot[3]*model_x[i][0] + rot[4]*model_x[i][1] + rot[5]*model_x[i][2] + com_coord[1];
-    atom_coord[i][2] = rot[6]*model_x[i][0] + rot[7]*model_x[i][1] + rot[8]*model_x[i][2] + com_coord[2];
+    atom_coord[i][0] = rot[0]*model_x[i][0] + 
+      rot[1]*model_x[i][1] + rot[2]*model_x[i][2] + com_coord[0];
+    atom_coord[i][1] = rot[3]*model_x[i][0] + 
+      rot[4]*model_x[i][1] + rot[5]*model_x[i][2] + com_coord[1];
+    atom_coord[i][2] = rot[6]*model_x[i][0] + 
+      rot[7]*model_x[i][1] + rot[8]*model_x[i][2] + com_coord[2];
 
     double xtmp[3];
     xtmp[0] = atom_coord[i][0];
@@ -799,11 +809,13 @@ void FixGCMC::attempt_molecule_insertion()
   }
 
   double insertion_energy_sum = 0.0;
-  MPI_Allreduce(&insertion_energy,&insertion_energy_sum,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&insertion_energy,&insertion_energy_sum,1,
+                MPI_DOUBLE,MPI_SUM,world);
 
-  if (random_equal->uniform() < zz*volume*natoms_per_molecule*exp(-beta*insertion_energy_sum)/(ngas+1)) {  
+  if (random_equal->uniform() < zz*volume*natoms_per_molecule*
+      exp(-beta*insertion_energy_sum)/(ngas+1)) {  
     maxmol++;
-    if (maxmol >= MAXSMALLINT) 
+    if (maxmol >= MAXTAGINT) 
       error->all(FLERR,"Fix gcmc ran out of available molecule IDs");
 
     tagint maxtag = 0;
@@ -883,13 +895,13 @@ void FixGCMC::attempt_molecule_insertion()
    compute particle's interaction energy with the rest of the system
 ------------------------------------------------------------------------- */
 
-double FixGCMC::energy(int i, int itype, int imolecule, double *coord)
+double FixGCMC::energy(int i, int itype, tagint imolecule, double *coord)
 {
   double delx,dely,delz,rsq;
 
   double **x = atom->x;
   int *type = atom->type;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
   pair = force->pair;
   cutsq = force->pair->cutsq;
@@ -941,7 +953,7 @@ int FixGCMC::pick_random_gas_atom()
 int FixGCMC::pick_random_gas_molecule()
 {
   int iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
-  int gas_molecule_id = 0;
+  tagint gas_molecule_id = 0;
   if ((iwhichglobal >= ngas_before) &&
       (iwhichglobal < ngas_before + ngas_local)) {
     int iwhichlocal = iwhichglobal - ngas_before;
@@ -949,8 +961,9 @@ int FixGCMC::pick_random_gas_molecule()
     gas_molecule_id = atom->molecule[i];
   }
 
-  int gas_molecule_id_all = 0;
-  MPI_Allreduce(&gas_molecule_id,&gas_molecule_id_all,1,MPI_INT,MPI_MAX,world);
+  tagint gas_molecule_id_all = 0;
+  MPI_Allreduce(&gas_molecule_id,&gas_molecule_id_all,1,
+                MPI_LMP_TAGINT,MPI_MAX,world);
   
   return gas_molecule_id_all;
 }
@@ -960,7 +973,7 @@ int FixGCMC::pick_random_gas_molecule()
    sum across all procs that own atoms of the given molecule
 ------------------------------------------------------------------------- */
 
-double FixGCMC::molecule_energy(int gas_molecule_id)
+double FixGCMC::molecule_energy(tagint gas_molecule_id)
 {
   double mol_energy = 0.0;
   for (int i = 0; i < atom->nlocal; i++)
@@ -1040,8 +1053,8 @@ void FixGCMC::get_model_molecule()
   if (atom->molecular) {
     for (int i = 0; i < atom->nlocal; i++) 
       maxmol = MAX(atom->molecule[i],maxmol);
-    int maxmol_all;
-    MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_INT,MPI_MAX,world);
+    tagint maxmol_all;
+    MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
     maxmol = maxmol_all;
   }
 
diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h
index 7fdd1f70c6df84c59e72f842c80d6702f36ba02e..ee532820b8ad1a39f15ec3a7a30141477c647d43 100644
--- a/src/MC/fix_gcmc.h
+++ b/src/MC/fix_gcmc.h
@@ -41,7 +41,7 @@ class FixGCMC : public Fix {
   void attempt_molecule_insertion();
   double energy(int, int, int, double *);
   int pick_random_gas_atom();
-  int pick_random_gas_molecule();
+  tagint pick_random_gas_molecule();
   double molecule_energy(int);
   void get_rotation_matrix(double, double *);
   void get_model_molecule();
@@ -63,9 +63,10 @@ class FixGCMC : public Fix {
   int regionflag;           // 0 = anywhere in box, 1 = specific region
   int iregion;              // GCMC region
   char *idregion;           // GCMC region id
-  bool pressure_flag;       // true if user specified reservoir pressure, false otherwise
+  bool pressure_flag;       // true if user specified reservoir pressure
+                            // else false
 
-  int maxmol;               // largest molecule tag across all existing atoms
+  tagint maxmol;            // largest molecule tag across all existing atoms
   int natoms_per_molecule;  // number of atoms in each gas molecule
 
   double ntranslation_attempts;
diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp
index ea8fb37be9a4578ff80258fd93910e92d93e1d38..f83af6d7d6cda172ffa0680556b3af3190579a27 100644
--- a/src/MISC/fix_deposit.cpp
+++ b/src/MISC/fix_deposit.cpp
@@ -523,7 +523,7 @@ void FixDeposit::pre_exchange()
 void FixDeposit::find_maxid()
 {
   tagint *tag = atom->tag;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
 
   tagint max = 0;
@@ -533,7 +533,7 @@ void FixDeposit::find_maxid()
   if (mode == MOLECULE && molecule) {
     max = 0;
     for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
-    MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world);
+    MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
   }
 }
 
diff --git a/src/MISC/fix_deposit.h b/src/MISC/fix_deposit.h
index 56bf34c6ddaf0cedd1fc2ae84eea3c656932ea1c..0c43a423d9fcd60ca53f19df49e4e265ee3db6df 100644
--- a/src/MISC/fix_deposit.h
+++ b/src/MISC/fix_deposit.h
@@ -55,8 +55,7 @@ class FixDeposit : public Fix {
   class Fix *fixrigid,*fixshake;
 
   int nfirst,ninserted;
-  tagint maxtag_all;
-  int maxmol_all;
+  tagint maxtag_all,maxmol_all;
   class RanPark *random;
 
   void find_maxid();
diff --git a/src/MISC/fix_evaporate.cpp b/src/MISC/fix_evaporate.cpp
index ab153f805e63c5ab7ba9fbfe2fac87d720df23d0..421c22e66638eea68d099030dc99899632018c9c 100644
--- a/src/MISC/fix_evaporate.cpp
+++ b/src/MISC/fix_evaporate.cpp
@@ -137,7 +137,7 @@ void FixEvaporate::init()
   // if molflag not set, warn if any deletable atom has a mol ID
 
   if (molflag == 0 && atom->molecule_flag) {
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
     int flag = 0;
@@ -234,8 +234,9 @@ void FixEvaporate::pre_exchange()
   // keep ndel,ndeltopo,ncount,nall,nbefore current after each mol deletion
 
   } else {
-    int me,proc,iatom,imolecule,ndelone,ndelall;
-    int *molecule = atom->molecule;
+    int me,proc,iatom,ndelone,ndelall;
+    tagint imol;
+    tagint *molecule = atom->molecule;
 
     ndeltopo[0] = ndeltopo[1] = ndeltopo[2] = ndeltopo[3] = 0;
 
@@ -247,7 +248,7 @@ void FixEvaporate::pre_exchange()
       if (iwhichglobal >= nbefore && iwhichglobal < nbefore + ncount) {
         iwhichlocal = iwhichglobal - nbefore;
         iatom = list[iwhichlocal];
-        imolecule = molecule[iatom];
+        imol = molecule[iatom];
         me = comm->me;
       } else me = -1;
 
@@ -257,10 +258,10 @@ void FixEvaporate::pre_exchange()
       // be careful to delete correct # of bond,angle,etc for newton on or off
 
       MPI_Allreduce(&me,&proc,1,MPI_INT,MPI_MAX,world);
-      MPI_Bcast(&imolecule,1,MPI_INT,proc,world);
+      MPI_Bcast(&imol,1,MPI_LMP_TAGINT,proc,world);
       ndelone = 0;
       for (i = 0; i < nlocal; i++) {
-        if (imolecule && molecule[i] == imolecule) {
+        if (imol && molecule[i] == imol) {
           mark[i] = 1;
           ndelone++;
 
diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp
index 9b303451c3bdd94fde16efe9cb16cbf50b109832..bafac2391f2b619f172b09711e4a3e6dbc001e34 100644
--- a/src/MOLECULE/atom_vec_angle.cpp
+++ b/src/MOLECULE/atom_vec_angle.cpp
@@ -478,7 +478,7 @@ void AtomVecAngle::unpack_border(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   }
 
   if (atom->nextra_border)
@@ -503,7 +503,7 @@ void AtomVecAngle::unpack_border_vel(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
@@ -524,7 +524,7 @@ int AtomVecAngle::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++)
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   return m;
 }
 
@@ -599,7 +599,7 @@ int AtomVecAngle::unpack_exchange(double *buf)
   mask[nlocal] = (int) ubuf(buf[m++]).i;
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -725,7 +725,7 @@ int AtomVecAngle::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -794,7 +794,7 @@ void AtomVecAngle::data_atom(double *coord, imageint imagetmp, char **values)
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -822,7 +822,7 @@ void AtomVecAngle::data_atom(double *coord, imageint imagetmp, char **values)
 
 int AtomVecAngle::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
 
   num_bond[nlocal] = 0;
   num_angle[nlocal] = 0;
@@ -867,8 +867,9 @@ int AtomVecAngle::pack_data_hybrid(int i, double *buf)
 void AtomVecAngle::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT " %d %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT 
+            " %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,
             buf[i][3],buf[i][4],buf[i][5],
             (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
@@ -881,7 +882,7 @@ void AtomVecAngle::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecAngle::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d",(int) ubuf(buf[0]).i);
+  fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[0]).i);
   return 1;
 }
 
diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp
index 7e5e63991895a1a898d73972e61782cb99bc9434..030a86ac862dce88993489c3e01abda42ee74455 100644
--- a/src/MOLECULE/atom_vec_bond.cpp
+++ b/src/MOLECULE/atom_vec_bond.cpp
@@ -436,7 +436,7 @@ int AtomVecBond::pack_border_hybrid(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    buf[m++] = molecule[j];
+    buf[m++] = ubuf(molecule[j]).d;
   }
   return m;
 }
@@ -457,7 +457,7 @@ void AtomVecBond::unpack_border(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   }
 
   if (atom->nextra_border)
@@ -482,7 +482,7 @@ void AtomVecBond::unpack_border_vel(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
@@ -503,7 +503,7 @@ int AtomVecBond::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++)
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   return m;
 }
 
@@ -570,7 +570,7 @@ int AtomVecBond::unpack_exchange(double *buf)
   mask[nlocal] = (int) ubuf(buf[m++]).i;
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -680,7 +680,7 @@ int AtomVecBond::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -740,7 +740,7 @@ void AtomVecBond::data_atom(double *coord, imageint imagetmp, char **values)
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -767,7 +767,7 @@ void AtomVecBond::data_atom(double *coord, imageint imagetmp, char **values)
 
 int AtomVecBond::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
 
   num_bond[nlocal] = 0;
 
@@ -811,8 +811,9 @@ int AtomVecBond::pack_data_hybrid(int i, double *buf)
 void AtomVecBond::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT " %d %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
+            " %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,
             buf[i][3],buf[i][4],buf[i][5],
             (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
@@ -825,7 +826,7 @@ void AtomVecBond::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecBond::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d",(int) ubuf(buf[0]).i);
+  fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[0]).i);
   return 1;
 }
 
diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp
index 09175dd0346a2d66018db849de628bf7ef8a0035..f132fdc59902e0dcc0e55a780af20fe73bd88fff 100644
--- a/src/MOLECULE/atom_vec_full.cpp
+++ b/src/MOLECULE/atom_vec_full.cpp
@@ -473,7 +473,7 @@ int AtomVecFull::pack_border_vel(int n, int *list, double *buf,
         buf[m++] = ubuf(type[j]).d;
         buf[m++] = ubuf(mask[j]).d;
         buf[m++] = q[j];
-        buf[m++] = molecule[j];
+        buf[m++] = ubuf(molecule[j]).d;
         buf[m++] = v[j][0];
         buf[m++] = v[j][1];
         buf[m++] = v[j][2];
@@ -544,7 +544,7 @@ void AtomVecFull::unpack_border(int n, int first, double *buf)
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
     q[i] = buf[m++];
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   }
 
   if (atom->nextra_border)
@@ -570,7 +570,7 @@ void AtomVecFull::unpack_border_vel(int n, int first, double *buf)
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
     q[i] = buf[m++];
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
@@ -592,7 +592,7 @@ int AtomVecFull::unpack_border_hybrid(int n, int first, double *buf)
   last = first + n;
   for (i = first; i < last; i++) {
     q[i] = buf[m++];
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   }
   return m;
 }
@@ -688,7 +688,7 @@ int AtomVecFull::unpack_exchange(double *buf)
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
   q[nlocal] = buf[m++];
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -853,7 +853,7 @@ int AtomVecFull::unpack_restart(double *buf)
   v[nlocal][2] = buf[m++];
 
   q[nlocal] = buf[m++];
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -943,7 +943,7 @@ void AtomVecFull::data_atom(double *coord, imageint imagetmp, char **values)
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -975,7 +975,7 @@ void AtomVecFull::data_atom(double *coord, imageint imagetmp, char **values)
 
 int AtomVecFull::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
   q[nlocal] = atof(values[1]);
 
   num_bond[nlocal] = 0;
@@ -1025,9 +1025,9 @@ int AtomVecFull::pack_data_hybrid(int i, double *buf)
 void AtomVecFull::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT 
-            " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
+            " %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,
             buf[i][3],buf[i][4],buf[i][5],buf[i][6],
             (int) ubuf(buf[i][7]).i,(int) ubuf(buf[i][8]).i,
@@ -1040,7 +1040,7 @@ void AtomVecFull::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecFull::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d %-1.16e",(int) ubuf(buf[0]).i,buf[1]);
+  fprintf(fp," " TAGINT_FORMAT " %-1.16e",(tagint) ubuf(buf[0]).i,buf[1]);
   return 2;
 }
 
diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp
index 2f2265b5df43cae19d40284281b85014ea016137..766da9e35443ed7eb924b8551dd4c2e4357b0e85 100644
--- a/src/MOLECULE/atom_vec_molecular.cpp
+++ b/src/MOLECULE/atom_vec_molecular.cpp
@@ -535,7 +535,7 @@ void AtomVecMolecular::unpack_border(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   }
 
   if (atom->nextra_border)
@@ -560,7 +560,7 @@ void AtomVecMolecular::unpack_border_vel(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     v[i][0] = buf[m++];
     v[i][1] = buf[m++];
     v[i][2] = buf[m++];
@@ -581,7 +581,7 @@ int AtomVecMolecular::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++)
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
   return m;
 }
 
@@ -674,7 +674,7 @@ int AtomVecMolecular::unpack_exchange(double *buf)
   mask[nlocal] = (int) ubuf(buf[m++]).i;
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -837,7 +837,7 @@ int AtomVecMolecular::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
 
   num_bond[nlocal] = (int) ubuf(buf[m++]).i;
   for (k = 0; k < num_bond[nlocal]; k++) {
@@ -927,7 +927,7 @@ void AtomVecMolecular::data_atom(double *coord, imageint imagetmp,
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -957,7 +957,7 @@ void AtomVecMolecular::data_atom(double *coord, imageint imagetmp,
 
 int AtomVecMolecular::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
 
   num_bond[nlocal] = 0;
   num_angle[nlocal] = 0;
@@ -1004,8 +1004,9 @@ int AtomVecMolecular::pack_data_hybrid(int i, double *buf)
 void AtomVecMolecular::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT " %d %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT 
+            " %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,
             buf[i][3],buf[i][4],buf[i][5],
             (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
@@ -1018,7 +1019,7 @@ void AtomVecMolecular::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecMolecular::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d",(int) ubuf(buf[0]).i);
+  fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[0]).i);
   return 1;
 }
 
diff --git a/src/MOLECULE/atom_vec_molecular.h b/src/MOLECULE/atom_vec_molecular.h
index 58aa039e1e3c2f2883b9d0fc8f966d1520ea3e26..293b44a8264e646183aebcf799e0e6e12301bea4 100644
--- a/src/MOLECULE/atom_vec_molecular.h
+++ b/src/MOLECULE/atom_vec_molecular.h
@@ -61,7 +61,7 @@ class AtomVecMolecular : public AtomVec {
   int *type,*mask;
   imageint *image;
   double **x,**v,**f;
-  int *molecule;
+  tagint *molecule;
   int **nspecial;
   tagint **special;
   int *num_bond;
diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp
index 1aea5795a1105b3cb6d62ba540e84162dbdd290d..cbada23c5a416b7cf118fdfb0679f3377d2387d3 100644
--- a/src/POEMS/fix_poems.cpp
+++ b/src/POEMS/fix_poems.cpp
@@ -144,7 +144,7 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
       if (!(mask[i] & groupbit)) natom2body[i] = 0;
 
   // each molecule in fix group is a rigid body
-  // maxmol = largest molecule #
+  // maxmol = largest molecule ID
   // ncount = # of atoms in each molecule (have to sum across procs)
   // nbody = # of non-zero ncount values
   // use nall as incremented ptr to set atom2body[] values for each atom
@@ -152,31 +152,36 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
   } else if (strcmp(arg[3],"molecule") == 0) {
     if (narg != 4) error->all(FLERR,"Illegal fix poems command");
     if (atom->molecular == 0)
-      error->all(FLERR,"Must use a molecular atom style with fix poems molecule");
+      error->all(FLERR,
+                 "Must use a molecular atom style with fix poems molecule");
 
     int *mask = atom->mask;
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int nlocal = atom->nlocal;
 
-    int maxmol = -1;
+    tagint maxmol_tag = -1;
     for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);
+      if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]);
 
-    int itmp;
-    MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world);
-    maxmol = itmp + 1;
+    tagint itmp;
+    MPI_Allreduce(&maxmol_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
+    if (itmp+1 > MAXSMALLINT) 
+      error->all(FLERR,"Too many molecules for fix poems");
+    maxmol = (int) itmp;
 
-    int *ncount = new int[maxmol];
-    for (i = 0; i < maxmol; i++) ncount[i] = 0;
+    int *ncount;
+    memory->create(ncount,maxmol+1,"rigid:ncount");
+    for (i = 0; i <= maxmol; i++) ncount[i] = 0;
 
     for (i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) ncount[molecule[i]]++;
 
-    int *nall = new int[maxmol];
-    MPI_Allreduce(ncount,nall,maxmol,MPI_INT,MPI_SUM,world);
+    int *nall;
+    memory->create(nall,maxmol+1,"rigid:ncount");
+    MPI_Allreduce(ncount,nall,maxmol+1,MPI_INT,MPI_SUM,world);
 
     nbody = 0;
-    for (i = 0; i < maxmol; i++)
+    for (i = 0; i <= maxmol; i++)
       if (nall[i]) nall[i] = nbody++;
       else nall[i] = -1;
 
@@ -188,8 +193,8 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
       }
     }
 
-    delete [] ncount;
-    delete [] nall;
+    memory->destroy(ncount);
+    memory->destroy(nall);
 
   } else error->all(FLERR,"Illegal fix poems command");
 
@@ -203,7 +208,8 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
     if (natom2body[i] > MAXBODY) flag = 1;
   int flagall;
   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
-  if (flagall) error->all(FLERR,"Atom in too many rigid bodies - boost MAXBODY");
+  if (flagall) 
+    error->all(FLERR,"Atom in too many rigid bodies - boost MAXBODY");
 
   // create all nbody-length arrays
 
diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index eee93e95c633a8ca675a27f9aa5a5060b0e861b5..c79fa2359e8475ebadcd40721239c86c255b7158 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -109,7 +109,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     }
 
   // each molecule in fix group is a rigid body
-  // maxmol = largest molecule #
+  // maxmol = largest molecule ID
   // ncount = # of atoms in each molecule (have to sum across procs)
   // nbody = # of non-zero ncount values
   // use nall as incremented ptr to set body[] values for each atom
@@ -121,16 +121,18 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
       error->all(FLERR,"Fix rigid molecule requires atom attribute molecule");
 
     int *mask = atom->mask;
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int nlocal = atom->nlocal;
 
-    maxmol = -1;
+    tagint maxmol_tag = -1;
     for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);
+      if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]);
 
-    int itmp;
-    MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world);
-    maxmol = itmp;
+    tagint itmp;
+    MPI_Allreduce(&maxmol_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
+    if (itmp+1 > MAXSMALLINT) 
+      error->all(FLERR,"Too many molecules for fix rigid");
+    maxmol = (int) itmp;
 
     int *ncount;
     memory->create(ncount,maxmol+1,"rigid:ncount");
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 46809d1b344c4f0c32daaae3064dec81b1122db6..78ef0ed0f28a758ca544973e2c262ac69b384126 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -107,15 +107,15 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   // maxmol = largest molecule #
 
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
 
   maxmol = -1;
   for (i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);
 
-  int itmp;
-  MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world);
+  tagint itmp;
+  MPI_Allreduce(&maxmol,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
   maxmol = itmp;
 
   // parse optional args
@@ -2147,8 +2147,8 @@ void FixRigidSmall::setup_bodies_dynamic()
 
 void FixRigidSmall::readfile(int which, double **array, int *inbody)
 {
-  int i,j,m,nchunk,id,eofflag;
-  int nlines;
+  int i,j,m,nchunk,eofflag,nlines;
+  tagint id;
   FILE *fp;
   char *eof,*start,*next,*buf;
   char line[MAXLINE];
@@ -2157,10 +2157,10 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody)
   // key = mol ID of bodies my atoms own
   // value = index into local body array
 
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
 
-  hash = new std::map<int,int>();
+  hash = new std::map<tagint,int>();
   for (int i = 0; i < nlocal; i++)
     if (bodyown[i] >= 0) (*hash)[atom->molecule[i]] = bodyown[i];
 
@@ -2219,28 +2219,28 @@ void FixRigidSmall::readfile(int which, double **array, int *inbody)
       for (j = 1; j < nwords; j++)
         values[j] = strtok(NULL," \t\n\r\f");
 
-      id = atoi(values[0]);
+      id = ATOTAGINT(values[0]);
       if (id <= 0 || id > maxmol) 
         error->all(FLERR,"Invalid rigid body ID in fix rigid/small file");
       if (hash->find(id) == hash->end()) {
         buf = next + 1;
         continue;
       }
-      id = (*hash)[id];
-      inbody[id] = 1;
+      m = (*hash)[id];
+      inbody[m] = 1;
 
       if (which == 0) {
-        body[id].mass = atof(values[1]);
-        body[id].xcm[0] = atof(values[2]);
-        body[id].xcm[1] = atof(values[3]);
-        body[id].xcm[2] = atof(values[4]);
+        body[m].mass = atof(values[1]);
+        body[m].xcm[0] = atof(values[2]);
+        body[m].xcm[1] = atof(values[3]);
+        body[m].xcm[2] = atof(values[4]);
       } else {
-        array[id][0] = atof(values[5]);
-        array[id][1] = atof(values[6]);
-        array[id][2] = atof(values[7]);
-        array[id][3] = atof(values[10]);
-        array[id][4] = atof(values[9]);
-        array[id][5] = atof(values[8]);
+        array[m][0] = atof(values[5]);
+        array[m][1] = atof(values[6]);
+        array[m][2] = atof(values[7]);
+        array[m][3] = atof(values[10]);
+        array[m][4] = atof(values[9]);
+        array[m][5] = atof(values[8]);
       }
 
       buf = next + 1;
diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h
index 6f73a9c03a0ed4a360a940f0831fc1fa2fd3305a..392a928bee587a234f2b3189f60fcfd45e8f9a23 100644
--- a/src/RIGID/fix_rigid_small.h
+++ b/src/RIGID/fix_rigid_small.h
@@ -81,7 +81,7 @@ class FixRigidSmall : public Fix {
   int firstflag;            // 1 for first-time setup of rigid bodies
   int commflag;             // various modes of forward/reverse comm
   int nbody;                // total # of rigid bodies
-  int maxmol;               // max mol-ID
+  tagint maxmol;            // max mol-ID
   double maxextent;         // furthest distance from body owner to body atom
 
   struct Body {
diff --git a/src/USER-LB/fix_lb_rigid_pc_sphere.cpp b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
index a0e3016698319224de5d8a7cf09e25f520bca8e9..6a1e4e531eb7f3abb3eb14e5faea79b9aa88f6c2 100755
--- a/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
+++ b/src/USER-LB/fix_lb_rigid_pc_sphere.cpp
@@ -95,31 +95,36 @@ FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) :
   } else if (strcmp(arg[3],"molecule") == 0) {
     iarg = 4;
     if (atom->molecular == 0)
-      error->all(FLERR,"Must use a molecular atom style with fix lb/rigid/pc/sphere molecule");
+      error->all(FLERR,"Must use a molecular atom style with "
+                 "fix lb/rigid/pc/sphere molecule");
 
     int *mask = atom->mask;
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     int nlocal = atom->nlocal;
 
-    int maxmol = -1;
+    tagint maxmol_tag = -1;
     for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) maxmol = MAX(maxmol,molecule[i]);
+      if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]);
 
-    int itmp;
-    MPI_Allreduce(&maxmol,&itmp,1,MPI_INT,MPI_MAX,world);
-    maxmol = itmp + 1;
+    tagint itmp;
+    MPI_Allreduce(&maxmol_tag,&itmp,1,MPI__LMP_TAGINT,MPI_MAX,world);
+    if (itmp+1 > MAXSMALLINT) 
+      error->all(FLERR,"Too many molecules for fix lb/rigid/pc/sphere");
+    maxmol = (int) itmp;
 
-    int *ncount = new int[maxmol];
-    for (i = 0; i < maxmol; i++) ncount[i] = 0;
+    int *ncount;
+    memory->create(ncount,maxmol+1,"rigid:ncount");
+    for (i = 0; i <= maxmol; i++) ncount[i] = 0;
 
     for (i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) ncount[molecule[i]]++;
 
-    int *nall = new int[maxmol];
-    MPI_Allreduce(ncount,nall,maxmol,MPI_INT,MPI_SUM,world);
+    int *nall;
+    memory->create(nall,maxmol+1,"rigid:ncount");
+    MPI_Allreduce(ncount,nall,maxmol+1,MPI_LMP_TAGINT,MPI_SUM,world);
 
     nbody = 0;
-    for (i = 0; i < maxmol; i++)
+    for (i = 0; i <= maxmol; i++)
       if (nall[i]) nall[i] = nbody++;
       else nall[i] = -1;
 
@@ -128,8 +133,8 @@ FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) :
       if (mask[i] & groupbit) body[i] = nall[molecule[i]];
     }
 
-    delete [] ncount;
-    delete [] nall;
+    memory->destroy(ncount);
+    memory->destroy(nall);
 
   // each listed group is a rigid body
   // check if all listed groups exist
@@ -140,7 +145,8 @@ FixLbRigidPCSphere::FixLbRigidPCSphere(LAMMPS *lmp, int narg, char **arg) :
     if (narg < 5) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
     nbody = atoi(arg[4]);
     if (nbody <= 0) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
-    if (narg < 5+nbody) error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
+    if (narg < 5+nbody) 
+      error->all(FLERR,"Illegal fix lb/rigid/pc/sphere command");
     iarg = 5 + nbody;
 
     int *igroups = new int[nbody];
diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp
index e609627897c6f993643e6c6703614521b42997e6..d1230555ce01b318a8065ab3ae1490ff86087a69 100644
--- a/src/USER-OMP/neigh_full_omp.cpp
+++ b/src/USER-OMP/neigh_full_omp.cpp
@@ -50,7 +50,7 @@ void Neighbor::full_nsq_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
 
@@ -138,7 +138,7 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -242,7 +242,7 @@ void Neighbor::full_bin_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -341,7 +341,7 @@ void Neighbor::full_bin_ghost_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -469,7 +469,7 @@ void Neighbor::full_multi_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
diff --git a/src/USER-OMP/neigh_gran_omp.cpp b/src/USER-OMP/neigh_gran_omp.cpp
index 5f0402f36397e230670e84b409b9a6b71c1492a5..e293102e9a2e628beee72f071afc68b1836a9e4e 100644
--- a/src/USER-OMP/neigh_gran_omp.cpp
+++ b/src/USER-OMP/neigh_gran_omp.cpp
@@ -63,7 +63,7 @@ void Neighbor::granular_nsq_no_newton_omp(NeighList *list)
   int *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
 
   int *ilist = list->ilist;
@@ -192,7 +192,7 @@ void Neighbor::granular_nsq_newton_omp(NeighList *list)
   int *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
 
   int *ilist = list->ilist;
@@ -304,7 +304,7 @@ void Neighbor::granular_bin_no_newton_omp(NeighList *list)
   int *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -443,7 +443,7 @@ void Neighbor::granular_bin_newton_omp(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -552,7 +552,7 @@ void Neighbor::granular_bin_newton_tri_omp(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp
index f501684a19642d7a79782c644da6426f02f0e2a0..76f36ea965d3ab35a0bca1594e0e97c21c649667 100644
--- a/src/USER-OMP/neigh_half_bin_omp.cpp
+++ b/src/USER-OMP/neigh_half_bin_omp.cpp
@@ -56,7 +56,7 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -159,7 +159,7 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -288,7 +288,7 @@ void Neighbor::half_bin_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -414,7 +414,7 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
diff --git a/src/USER-OMP/neigh_half_multi_omp.cpp b/src/USER-OMP/neigh_half_multi_omp.cpp
index b46594870fa38100b833d3d457ee4e605e030abe..a804e8008f519e65ac996967acd797441d9c8eeb 100644
--- a/src/USER-OMP/neigh_half_multi_omp.cpp
+++ b/src/USER-OMP/neigh_half_multi_omp.cpp
@@ -58,7 +58,7 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -164,7 +164,7 @@ void Neighbor::half_multi_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -298,7 +298,7 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp
index a98c4721f3d5dfe2a91be852874e191a618d5ac2..fa90945d5567f94ec792d5abaa0871e619f247bf 100644
--- a/src/USER-OMP/neigh_half_nsq_omp.cpp
+++ b/src/USER-OMP/neigh_half_nsq_omp.cpp
@@ -52,7 +52,7 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -141,7 +141,7 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -250,7 +250,7 @@ void Neighbor::half_nsq_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
 
diff --git a/src/USER-OMP/neigh_respa_omp.cpp b/src/USER-OMP/neigh_respa_omp.cpp
index 11c7fd9ae307a5c00ccd853461fdfbf8ec349213..458b2bbac4f2b873ea08a86bbfba0b6d757bf95c 100644
--- a/src/USER-OMP/neigh_respa_omp.cpp
+++ b/src/USER-OMP/neigh_respa_omp.cpp
@@ -59,7 +59,7 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
 
@@ -210,7 +210,7 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nall = atom->nlocal + atom->nghost;
   int molecular = atom->molecular;
 
@@ -382,7 +382,7 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -545,7 +545,7 @@ void Neighbor::respa_bin_newton_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
@@ -749,7 +749,7 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int molecular = atom->molecular;
 
   int *ilist = list->ilist;
diff --git a/src/USER-REAXC/fix_reaxc_bonds.cpp b/src/USER-REAXC/fix_reaxc_bonds.cpp
index 3ced20a2891765f8ea66b31fa424a6e7ec0da715..8d8818a49b327cbd24cce2a2755323e5be5579b6 100644
--- a/src/USER-REAXC/fix_reaxc_bonds.cpp
+++ b/src/USER-REAXC/fix_reaxc_bonds.cpp
@@ -290,7 +290,7 @@ void FixReaxCBonds::RecvBuffer(double *buf, int nbuf, int nbuf_local,
         avqtmp = buf[j+3];
         numbonds = nint(buf[j+4]);
 
-        fprintf(fp," %d %d %d",itag,itype,numbonds);
+        fprintf(fp," " TAGINT_FORMAT " %d %d",itag,itype,numbonds);
 
         for (k = 5; k < 5+numbonds; k++) {
           jtag = static_cast<tagint> (buf[j+k]);
@@ -298,7 +298,7 @@ void FixReaxCBonds::RecvBuffer(double *buf, int nbuf, int nbuf_local,
         }
         j += (5+numbonds);
 
-        fprintf(fp," %d",nint(buf[j]));
+        fprintf(fp," " TAGINT_FORMAT,static_cast<tagint> (buf[j]));
         j ++;
 
         for (k = 0; k < numbonds; k++) {
diff --git a/src/atom.h b/src/atom.h
index 8a5a0b643f524e4e57a244f068a691bf37072041..c552cad6c239a49e4a9a8a485e73031206d6647a 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -50,7 +50,8 @@ class Atom : protected Pointers {
   imageint *image;
   double **x,**v,**f;
 
-  int *molecule,*molindex,*molatom;
+  tagint *molecule;
+  int *molindex,*molatom;
   double *q,**mu;
   double **omega,**angmom,**torque;
   double *radius,*rmass,*vfrac,*s0;
diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp
index 5c89d237fc7ec51e68927b1c654fe3c86dcec1c3..bd5a5ea2b3aa9c82ddc85ad650cc134dcab65d3f 100644
--- a/src/atom_vec_line.cpp
+++ b/src/atom_vec_line.cpp
@@ -642,7 +642,7 @@ int AtomVecLine::pack_border_hybrid(int n, int *list, double *buf)
   m = 0;
   for (i = 0; i < n; i++) {
     j = list[i];
-    buf[m++] = molecule[j];
+    buf[m++] = ubuf(molecule[j]).d;
     if (line[j] < 0) buf[m++] = ubuf(0).d;
     else {
       buf[m++] = ubuf(1).d;
@@ -669,7 +669,7 @@ void AtomVecLine::unpack_border(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     line[i] = (int) ubuf(buf[m++]).i;
     if (line[i] == 0) line[i] = -1;
     else {
@@ -705,7 +705,7 @@ void AtomVecLine::unpack_border_vel(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     line[i] = (int) ubuf(buf[m++]).i;
     if (line[i] == 0) line[i] = -1;
     else {
@@ -740,7 +740,7 @@ int AtomVecLine::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     line[i] = (int) ubuf(buf[m++]).i;
     if (line[i] == 0) line[i] = -1;
     else {
@@ -816,7 +816,7 @@ int AtomVecLine::unpack_exchange(double *buf)
   mask[nlocal] = (int) ubuf(buf[m++]).i;
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
   rmass[nlocal] = buf[m++];
   omega[nlocal][0] = buf[m++];
   omega[nlocal][1] = buf[m++];
@@ -931,7 +931,7 @@ int AtomVecLine::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
   rmass[nlocal] = buf[m++];
   omega[nlocal][0] = buf[m++];
   omega[nlocal][1] = buf[m++];
@@ -1000,7 +1000,7 @@ void AtomVecLine::data_atom(double *coord, imageint imagetmp, char **values)
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -1038,7 +1038,7 @@ void AtomVecLine::data_atom(double *coord, imageint imagetmp, char **values)
 
 int AtomVecLine::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
 
   line[nlocal] = atoi(values[1]);
   if (line[nlocal] == 0) line[nlocal] = -1;
@@ -1166,9 +1166,9 @@ int AtomVecLine::pack_data_hybrid(int i, double *buf)
 void AtomVecLine::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT 
-            " %d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
+            " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,(int) ubuf(buf[i][3]).i,
             buf[i][4],buf[i][5],buf[i][6],buf[i][7],
             (int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i,
@@ -1181,7 +1181,8 @@ void AtomVecLine::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecLine::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d %d %-1.16e",(int) ubuf(buf[0]).i,(int) ubuf(buf[1]).i,buf[2]);
+  fprintf(fp," " TAGINT_FORMAT " %d %-1.16e",
+          (tagint) ubuf(buf[0]).i,(int) ubuf(buf[1]).i,buf[2]);
   return 3;
 }
 
diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp
index d170c8b5d0efa1ff28881b1bd836a9c3f7d2c2e9..20528ff0d6230c12084f3a322201930e5e8da15f 100644
--- a/src/atom_vec_tri.cpp
+++ b/src/atom_vec_tri.cpp
@@ -885,7 +885,7 @@ void AtomVecTri::unpack_border(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     tri[i] = (int) ubuf(buf[m++]).i;
     if (tri[i] == 0) tri[i] = -1;
     else {
@@ -941,7 +941,7 @@ void AtomVecTri::unpack_border_vel(int n, int first, double *buf)
     tag[i] = (tagint) ubuf(buf[m++]).i;
     type[i] = (int) ubuf(buf[m++]).i;
     mask[i] = (int) ubuf(buf[m++]).i;
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     tri[i] = (int) ubuf(buf[m++]).i;
     if (tri[i] == 0) tri[i] = -1;
     else {
@@ -996,7 +996,7 @@ int AtomVecTri::unpack_border_hybrid(int n, int first, double *buf)
   m = 0;
   last = first + n;
   for (i = first; i < last; i++) {
-    molecule[i] = (int) ubuf(buf[m++]).i;
+    molecule[i] = (tagint) ubuf(buf[m++]).i;
     tri[i] = (int) ubuf(buf[m++]).i;
     if (tri[i] == 0) tri[i] = -1;
     else {
@@ -1110,7 +1110,7 @@ int AtomVecTri::unpack_exchange(double *buf)
   mask[nlocal] = (int) ubuf(buf[m++]).i;
   image[nlocal] = (imageint) ubuf(buf[m++]).i;
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
   rmass[nlocal] = buf[m++];
   angmom[nlocal][0] = buf[m++];
   angmom[nlocal][1] = buf[m++];
@@ -1263,7 +1263,7 @@ int AtomVecTri::unpack_restart(double *buf)
   v[nlocal][1] = buf[m++];
   v[nlocal][2] = buf[m++];
 
-  molecule[nlocal] = (int) ubuf(buf[m++]).i;
+  molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
   rmass[nlocal] = buf[m++];
   angmom[nlocal][0] = buf[m++];
   angmom[nlocal][1] = buf[m++];
@@ -1351,7 +1351,7 @@ void AtomVecTri::data_atom(double *coord, imageint imagetmp, char **values)
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = ATOTAGINT(values[0]);
-  molecule[nlocal] = atoi(values[1]);
+  molecule[nlocal] = ATOTAGINT(values[1]);
   type[nlocal] = atoi(values[2]);
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one(FLERR,"Invalid atom type in Atoms section of data file");
@@ -1389,7 +1389,7 @@ void AtomVecTri::data_atom(double *coord, imageint imagetmp, char **values)
 
 int AtomVecTri::data_atom_hybrid(int nlocal, char **values)
 {
-  molecule[nlocal] = atoi(values[0]);
+  molecule[nlocal] = ATOTAGINT(values[0]);
 
   tri[nlocal] = atoi(values[1]);
   if (tri[nlocal] == 0) tri[nlocal] = -1;
@@ -1614,9 +1614,9 @@ int AtomVecTri::pack_data_hybrid(int i, double *buf)
 void AtomVecTri::write_data(FILE *fp, int n, double **buf)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,TAGINT_FORMAT 
-            " %d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
-            (tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
+    fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
+            " %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
+            (tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (int) ubuf(buf[i][2]).i,(int) ubuf(buf[i][3]).i,
             buf[i][4],buf[i][5],buf[i][6],buf[i][7],
             (int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i,
@@ -1629,7 +1629,8 @@ void AtomVecTri::write_data(FILE *fp, int n, double **buf)
 
 int AtomVecTri::write_data_hybrid(FILE *fp, double *buf)
 {
-  fprintf(fp," %d %d %-1.16e",(int) ubuf(buf[0]).i,(int) ubuf(buf[1]).i,buf[2]);
+  fprintf(fp," " TAGINT_FORMAT " %d %-1.16e",
+          (tagint) ubuf(buf[0]).i,(int) ubuf(buf[1]).i,buf[2]);
   return 3;
 }
 
diff --git a/src/compute.cpp b/src/compute.cpp
index 9be14029dc81de5f691a5540da74e7acbf01fb70..f9a95e8d10b97b8cdcf3f1bd60bb6046b1a34780 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -29,7 +29,7 @@
 using namespace LAMMPS_NS;
 
 #define DELTA 4
-#define BIG 2000000000
+#define BIG MAXTAGINT
 
 /* ---------------------------------------------------------------------- */
 
@@ -215,7 +215,7 @@ void Compute::clearstep()
          return idlo and idhi
 ------------------------------------------------------------------------- */
 
-int Compute::molecules_in_group(int &idlo, int &idhi)
+int Compute::molecules_in_group(tagint &idlo, tagint &idhi)
 {
   int i;
 
@@ -225,12 +225,12 @@ int Compute::molecules_in_group(int &idlo, int &idhi)
   // find lo/hi molecule ID for any atom in group
   // warn if atom in group has ID = 0
 
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  int lo = BIG;
-  int hi = -BIG;
+  tagint lo = BIG;
+  tagint hi = -BIG;
   int flag = 0;
   for (i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -245,14 +245,18 @@ int Compute::molecules_in_group(int &idlo, int &idhi)
     error->warning(FLERR,"Atom with molecule ID = 0 included in "
                    "compute molecule group");
 
-  MPI_Allreduce(&lo,&idlo,1,MPI_INT,MPI_MIN,world);
-  MPI_Allreduce(&hi,&idhi,1,MPI_INT,MPI_MAX,world);
+  MPI_Allreduce(&lo,&idlo,1,MPI_LMP_TAGINT,MPI_MIN,world);
+  MPI_Allreduce(&hi,&idhi,1,MPI_LMP_TAGINT,MPI_MAX,world);
   if (idlo == BIG) return 0;
 
   // molmap = vector of length nlen
   // set to 1 for IDs that appear in group across all procs, else 0
 
-  int nlen = idhi-idlo+1;
+  tagint nlen_tag = idhi-idlo+1;
+  if (nlen_tag > MAXSMALLINT) 
+    error->all(FLERR,"Too many molecules for compute");
+  int nlen = (int) nlen_tag;
+
   memory->create(molmap,nlen,"compute:molmap");
   for (i = 0; i < nlen; i++) molmap[i] = 0;
 
diff --git a/src/compute.h b/src/compute.h
index e7cb31bae16967f2293bd5970c6671cf59edfdc1..3ae91e5b1aa72c407228312a4d7a80aac5217235 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -130,7 +130,7 @@ class Compute : protected Pointers {
 
   int *molmap;                 // convert molecule ID to local index
 
-  int molecules_in_group(int &, int &);
+  int molecules_in_group(tagint &, tagint &);
 
   inline int sbmask(int j) {
     return j >> SBBITS & 3;
diff --git a/src/compute_atom_molecule.h b/src/compute_atom_molecule.h
index 769f58d5ec5c538b7095e15a2aeb4f829f0730e1..cd135e47d8ded98b4f82d4b9a3010efe681ce8b0 100644
--- a/src/compute_atom_molecule.h
+++ b/src/compute_atom_molecule.h
@@ -35,7 +35,7 @@ class ComputeAtomMolecule : public Compute {
 
  private:
   int nvalues,nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
 
   int *which,*argindex,*flavor,*value2index;
   char **ids;
diff --git a/src/compute_com_molecule.cpp b/src/compute_com_molecule.cpp
index 130072cc630e52d0c7850f10c2f46e534bdaabc9..74b2cb5ac3689d0a9ec9242ec1084bf770453866 100644
--- a/src/compute_com_molecule.cpp
+++ b/src/compute_com_molecule.cpp
@@ -48,18 +48,18 @@ ComputeCOMMolecule::ComputeCOMMolecule(LAMMPS *lmp, int narg, char **arg) :
   // compute masstotal for each molecule
 
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  int i,imol;
+  tagint imol;
   double massone;
 
-  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
@@ -95,18 +95,18 @@ void ComputeCOMMolecule::init()
 
 void ComputeCOMMolecule::compute_array()
 {
-  int i,imol;
+  tagint imol;
   double massone;
   double unwrap[3];
 
   invoked_array = update->ntimestep;
 
-  for (i = 0; i < nmolecules; i++)
+  for (int i = 0; i < nmolecules; i++)
     com[i][0] = com[i][1] = com[i][2] = 0.0;
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
@@ -128,7 +128,7 @@ void ComputeCOMMolecule::compute_array()
 
   MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
                 MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < nmolecules; i++) {
+  for (int i = 0; i < nmolecules; i++) {
     comall[i][0] /= masstotal[i];
     comall[i][1] /= masstotal[i];
     comall[i][2] /= masstotal[i];
@@ -141,8 +141,8 @@ void ComputeCOMMolecule::compute_array()
 
 double ComputeCOMMolecule::memory_usage()
 {
-  double bytes = 2*nmolecules * sizeof(double);
+  double bytes = (bigint) nmolecules * 2 * sizeof(double);
   if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
-  bytes += 2*nmolecules*3 * sizeof(double);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
   return bytes;
 }
diff --git a/src/compute_com_molecule.h b/src/compute_com_molecule.h
index a4f941c52aa8e563503d46a8c1cea7dbd15ed25f..3e6fdedd99e8b27edfe681d039daac6509ad4bdb 100644
--- a/src/compute_com_molecule.h
+++ b/src/compute_com_molecule.h
@@ -34,7 +34,7 @@ class ComputeCOMMolecule : public Compute {
 
  private:
   int nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
 
   double *massproc,*masstotal;
   double **com,**comall;
diff --git a/src/compute_gyration_molecule.cpp b/src/compute_gyration_molecule.cpp
index bf4a70a296df54e48930852af5b62a3aa07cfbda..f81a373830e59a88d31b8e39ff0abaef786b0509 100644
--- a/src/compute_gyration_molecule.cpp
+++ b/src/compute_gyration_molecule.cpp
@@ -73,18 +73,18 @@ ComputeGyrationMolecule::ComputeGyrationMolecule(LAMMPS *lmp,
   // compute masstotal for each molecule
 
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  int i,imol;
+  tagint imol;
   double massone;
 
-  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
@@ -122,7 +122,7 @@ void ComputeGyrationMolecule::init()
 
 void ComputeGyrationMolecule::compute_vector()
 {
-  int i,imol;
+  tagint imol;
   double dx,dy,dz,massone;
   double unwrap[3];
 
@@ -130,18 +130,18 @@ void ComputeGyrationMolecule::compute_vector()
 
   molcom();
 
-  for (i = 0; i < nmolecules; i++) rg[i] = 0.0;
+  for (int i = 0; i < nmolecules; i++) rg[i] = 0.0;
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       imol = molecule[i];
       if (molmap) imol = molmap[imol-idlo];
@@ -157,14 +157,16 @@ void ComputeGyrationMolecule::compute_vector()
 
   MPI_Allreduce(rg,vector,nmolecules,MPI_DOUBLE,MPI_SUM,world);
 
-  for (i = 0; i < nmolecules; i++) vector[i] = sqrt(vector[i]/masstotal[i]);
+  for (int i = 0; i < nmolecules; i++)
+    vector[i] = sqrt(vector[i]/masstotal[i]);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeGyrationMolecule::compute_array()
 {
-  int i,j,imol;
+  int i,j;
+  tagint imol;
   double dx,dy,dz,massone;
   double unwrap[3];
 
@@ -178,7 +180,7 @@ void ComputeGyrationMolecule::compute_array()
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
@@ -220,23 +222,23 @@ void ComputeGyrationMolecule::compute_array()
 
 void ComputeGyrationMolecule::molcom()
 {
-  int i,imol;
+  tagint imol;
   double dx,dy,dz,massone;
   double unwrap[3];
 
-  for (i = 0; i < nmolecules; i++)
+  for (int i = 0; i < nmolecules; i++)
     com[i][0] = com[i][1] = com[i][2] = 0.0;
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       imol = molecule[i];
       if (molmap) imol = molmap[imol-idlo];
@@ -251,7 +253,7 @@ void ComputeGyrationMolecule::molcom()
 
   MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
                 MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < nmolecules; i++) {
+  for (int i = 0; i < nmolecules; i++) {
     comall[i][0] /= masstotal[i];
     comall[i][1] /= masstotal[i];
     comall[i][2] /= masstotal[i];
@@ -264,10 +266,10 @@ void ComputeGyrationMolecule::molcom()
 
 double ComputeGyrationMolecule::memory_usage()
 {
-  double bytes = 2*nmolecules * sizeof(double);
+  double bytes = (bigint) nmolecules * 2 * sizeof(double);
   if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
-  bytes += 2*nmolecules*3 * sizeof(double);
-  if (tensor) bytes += 2*6*nmolecules * sizeof(double);
-  else bytes += 2*nmolecules * sizeof(double);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+  if (tensor) bytes += (bigint) nmolecules * 2*6 * sizeof(double);
+  else bytes += (bigint) nmolecules * 2 * sizeof(double);
   return bytes;
 }
diff --git a/src/compute_gyration_molecule.h b/src/compute_gyration_molecule.h
index 99910782fae8eda1d98b28c23adcf28afbe1c6d2..78b49231a2fe2030109aad0af4ee374b76ddab40 100644
--- a/src/compute_gyration_molecule.h
+++ b/src/compute_gyration_molecule.h
@@ -36,7 +36,7 @@ class ComputeGyrationMolecule : public Compute {
  private:
   int tensor;
   int nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
 
   double *massproc,*masstotal;
   double **com,**comall;
diff --git a/src/compute_inertia_molecule.cpp b/src/compute_inertia_molecule.cpp
index 8cdb3cb64cdc7b60396c881144202b4c683269ed..df9958a04a82f78635b3631923e1f0abec808fd6 100644
--- a/src/compute_inertia_molecule.cpp
+++ b/src/compute_inertia_molecule.cpp
@@ -51,18 +51,18 @@ ComputeInertiaMolecule(LAMMPS *lmp, int narg, char **arg) :
   // compute masstotal for each molecule
 
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  int i,imol;
+  tagint imol;
   double massone;
 
-  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
@@ -100,7 +100,8 @@ void ComputeInertiaMolecule::init()
 
 void ComputeInertiaMolecule::compute_array()
 {
-  int i,j,imol;
+  int i,j;
+  tagint imol;
   double dx,dy,dz,massone;
   double unwrap[3];
 
@@ -108,7 +109,7 @@ void ComputeInertiaMolecule::compute_array()
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
@@ -176,9 +177,9 @@ void ComputeInertiaMolecule::compute_array()
 
 double ComputeInertiaMolecule::memory_usage()
 {
-  double bytes = 2*nmolecules * sizeof(double);
+  double bytes = (bigint) nmolecules * 2 * sizeof(double);
   if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
-  bytes += 2*nmolecules*3 * sizeof(double);
-  bytes += 2*nmolecules*6 * sizeof(double);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+  bytes += (bigint) nmolecules * 2*6 * sizeof(double);
   return bytes;
 }
diff --git a/src/compute_inertia_molecule.h b/src/compute_inertia_molecule.h
index 80657108e7c7a982a9866c33f2c7eafce4129c44..a081dd1ad54f064054074ea09e2c85396db9a1bc 100644
--- a/src/compute_inertia_molecule.h
+++ b/src/compute_inertia_molecule.h
@@ -34,7 +34,7 @@ class ComputeInertiaMolecule : public Compute {
 
  private:
   int nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
 
   double *massproc,*masstotal;
   double **com,**comall;
diff --git a/src/compute_msd_molecule.cpp b/src/compute_msd_molecule.cpp
index 736ac9cfc6de75370e7ae87467fe4e38175a7cf2..515436ba3a0bdb5e3b0abcd1d67b1534a19eccfc 100644
--- a/src/compute_msd_molecule.cpp
+++ b/src/compute_msd_molecule.cpp
@@ -50,18 +50,18 @@ ComputeMSDMolecule::ComputeMSDMolecule(LAMMPS *lmp, int narg, char **arg) :
   // compute masstotal for each molecule
 
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
   int nlocal = atom->nlocal;
 
-  int i,imol;
+  tagint imol;
   double massone;
 
-  for (i = 0; i < nmolecules; i++) massproc[i] = 0.0;
+  for (int i = 0; i < nmolecules; i++) massproc[i] = 0.0;
 
-  for (i = 0; i < nlocal; i++)
+  for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
       if (rmass) massone = rmass[i];
       else massone = mass[type[i]];
@@ -77,7 +77,7 @@ ComputeMSDMolecule::ComputeMSDMolecule(LAMMPS *lmp, int narg, char **arg) :
 
   firstflag = 1;
   compute_array();
-  for (i = 0; i < nmolecules; i++) {
+  for (int i = 0; i < nmolecules; i++) {
     cominit[i][0] = comall[i][0];
     cominit[i][1] = comall[i][1];
     cominit[i][2] = comall[i][2];
@@ -110,7 +110,7 @@ void ComputeMSDMolecule::init()
 
 void ComputeMSDMolecule::compute_array()
 {
-  int i,imol;
+  tagint imol;
   double dx,dy,dz,massone;
   double unwrap[3];
 
@@ -118,12 +118,12 @@ void ComputeMSDMolecule::compute_array()
 
   // compute current COM positions
 
-  for (i = 0; i < nmolecules; i++)
+  for (int i = 0; i < nmolecules; i++)
     com[i][0] = com[i][1] = com[i][2] = 0.0;
 
   double **x = atom->x;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *type = atom->type;
   imageint *image = atom->image;
   double *mass = atom->mass;
@@ -145,7 +145,7 @@ void ComputeMSDMolecule::compute_array()
 
   MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,
                 MPI_DOUBLE,MPI_SUM,world);
-  for (i = 0; i < nmolecules; i++) {
+  for (int i = 0; i < nmolecules; i++) {
     comall[i][0] /= masstotal[i];
     comall[i][1] /= masstotal[i];
     comall[i][2] /= masstotal[i];
@@ -156,7 +156,7 @@ void ComputeMSDMolecule::compute_array()
 
   if (firstflag) return;
 
-  for (i = 0; i < nmolecules; i++) {
+  for (int i = 0; i < nmolecules; i++) {
     dx = comall[i][0] - cominit[i][0];
     dy = comall[i][1] - cominit[i][1];
     dz = comall[i][2] - cominit[i][2];
@@ -173,9 +173,9 @@ void ComputeMSDMolecule::compute_array()
 
 double ComputeMSDMolecule::memory_usage()
 {
-  double bytes = 2*nmolecules * sizeof(double);
+  double bytes = (bigint) nmolecules * 2 * sizeof(double);
   if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
-  bytes += 2*nmolecules*3 * sizeof(double);
-  bytes += nmolecules*4 * sizeof(double);
+  bytes += (bigint) nmolecules * 2*3 * sizeof(double);
+  bytes += (bigint) nmolecules * 4 * sizeof(double);
   return bytes;
 }
diff --git a/src/compute_msd_molecule.h b/src/compute_msd_molecule.h
index b4c5bd53a19e3c6e410783836da5ab4aba96ca3b..839698cabd638a07feaadb81f63295b5cd421c18 100644
--- a/src/compute_msd_molecule.h
+++ b/src/compute_msd_molecule.h
@@ -34,7 +34,7 @@ class ComputeMSDMolecule : public Compute {
 
  private:
   int nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
   int firstflag;
 
   double *massproc,*masstotal;
diff --git a/src/compute_property_molecule.cpp b/src/compute_property_molecule.cpp
index c0bcafc9ea6a2660080495916c0c3c873b5379e6..dafe2db68098072f88a4997dd7f4c92decfcb28f 100644
--- a/src/compute_property_molecule.cpp
+++ b/src/compute_property_molecule.cpp
@@ -43,7 +43,8 @@ ComputePropertyMolecule(LAMMPS *lmp, int narg, char **arg) :
       pack_choice[i] = &ComputePropertyMolecule::pack_mol;
     else if (strcmp(arg[iarg],"count") == 0)
       pack_choice[i] = &ComputePropertyMolecule::pack_count;
-    else error->all(FLERR,"Invalid keyword in compute property/molecule command");
+    else error->all(FLERR,
+                    "Invalid keyword in compute property/molecule command");
   }
 
   // setup molecule-based data
@@ -116,7 +117,7 @@ void ComputePropertyMolecule::compute_array()
 
 double ComputePropertyMolecule::memory_usage()
 {
-  double bytes = nmolecules*nvalues * sizeof(double);
+  double bytes = (bigint) nmolecules * nvalues * sizeof(double);
   if (molmap) bytes += (idhi-idlo+1) * sizeof(int);
   return bytes;
 }
@@ -129,7 +130,7 @@ double ComputePropertyMolecule::memory_usage()
 
 void ComputePropertyMolecule::pack_mol(int n)
 {
-  for (int m = idlo; m <= idhi; m++)
+  for (tagint m = idlo; m <= idhi; m++)
     if (molmap == NULL || molmap[m-idlo] >= 0) {
       buf[n] = m;
       n += nvalues;
@@ -140,12 +141,13 @@ void ComputePropertyMolecule::pack_mol(int n)
 
 void ComputePropertyMolecule::pack_count(int n)
 {
-  int i,m,imol;
+  int i,m;
+  tagint imol;
 
   int *count_one = new int[nmolecules];
   for (m = 0; m < nmolecules; m++) count_one[m] = 0;
 
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
diff --git a/src/compute_property_molecule.h b/src/compute_property_molecule.h
index 0c5ae68df08c2b515af08bdfbb208e39da165ee7..27cebe8fe18cd00d901f58bb0cf02613773619e8 100644
--- a/src/compute_property_molecule.h
+++ b/src/compute_property_molecule.h
@@ -35,7 +35,7 @@ class ComputePropertyMolecule : public Compute {
 
  private:
   int nvalues,nmolecules;
-  int idlo,idhi;
+  tagint idlo,idhi;
 
   double *buf;
 
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index 1e1f5716da8962526a8d2617f2f6e44756a87788..40c7fe4ec4ed83574e3986fb1a32d0dfea753bb7 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -292,13 +292,13 @@ void CreateAtoms::command(int narg, char **arg)
     // moloffset = max molecule ID for all molecules owned by previous procs
     //             including molecules existing before this creation
 
-    int moloffset;
-    int *molecule = atom->molecule;
+    tagint moloffset;
+    tagint *molecule = atom->molecule;
     if (molecule) {
-      int max = 0;
+      tagint max = 0;
       for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]);
-      int maxmol;
-      MPI_Allreduce(&max,&maxmol,1,MPI_INT,MPI_MAX,world);
+      tagint maxmol;
+      MPI_Allreduce(&max,&maxmol,1,MPI_LMP_TAGINT,MPI_MAX,world);
       MPI_Scan(&molcreate,&moloffset,1,MPI_INT,MPI_SUM,world);
       moloffset = moloffset - molcreate + maxmol;
     }
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index ac1c354d2907c89f06ce6941d6214c1130952424..5c353c7fe561032d10add5c414b4f87e7bb11232 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -172,15 +172,15 @@ void DeleteAtoms::delete_region(int narg, char **arg)
   // store list of molecule IDs I delete atoms from in list
   // pass list from proc to proc via ring communication
 
-  hash = new std::map<int,int>();
+  hash = new std::map<tagint,int>();
 
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   for (int i = 0; i < nlocal; i++)
     if (dlist[i] && hash->find(molecule[i]) == hash->end())
       (*hash)[molecule[i]] = 1;
 
   int n = hash->size();
-  int *list;
+  tagint *list;
   memory->create(list,n,"delete_atoms:list");
 
   n = 0;
@@ -188,7 +188,7 @@ void DeleteAtoms::delete_region(int narg, char **arg)
   for (pos = hash->begin(); pos != hash->end(); ++pos) list[n++] = pos->first;
 
   cptr = this;
-  comm->ring(n,sizeof(int),list,1,molring,NULL);
+  comm->ring(n,sizeof(tagint),list,1,molring,NULL);
 
   delete hash;
   memory->destroy(list);
@@ -197,16 +197,16 @@ void DeleteAtoms::delete_region(int narg, char **arg)
 /* ----------------------------------------------------------------------
    callback from comm->ring()
    cbuf = list of N molecule IDs, put them in hash
-   loop over my atoms, if matches moleculed ID in hash, delete that atom
+   loop over my atoms, if matches molecule ID in hash, delete that atom
 ------------------------------------------------------------------------- */
 
 void DeleteAtoms::molring(int n, char *cbuf)
 {
-  int *list = (int *) cbuf;
+  tagint *list = (int *) cbuf;
   int *dlist = cptr->dlist;
-  std::map<int,int> *hash = cptr->hash;
+  std::map<tagint,int> *hash = cptr->hash;
   int nlocal = cptr->atom->nlocal;
-  int *molecule = cptr->atom->molecule;
+  tagint *molecule = cptr->atom->molecule;
 
   hash->clear();
   for (int i = 0; i < n; i++) (*hash)[list[i]] = 1;
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index f0c0371044e5fee6b5c689336772712dac174507..9619214b94ba9c08945178c016e1583fc6bb259f 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -446,7 +446,7 @@ int DumpCustom::count()
         if (!atom->molecule_flag)
           error->all(FLERR,
                      "Threshhold for an atom property that isn't allocated");
-        int *molecule = atom->molecule;
+        tagint *molecule = atom->molecule;
         for (i = 0; i < nlocal; i++) dchoose[i] = molecule[i];
         ptr = dchoose;
         nstride = 1;
@@ -1725,7 +1725,7 @@ void DumpCustom::pack_id(int n)
 
 void DumpCustom::pack_molecule(int n)
 {
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
 
   for (int i = 0; i < nchoose; i++) {
     buf[n] = molecule[clist[i]];
diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp
index a1be2b556bb5c6ca30a00da16763d7a4fb77d951..115750e271720c5863945bc4015d47dc8a99cf79 100644
--- a/src/fix_property_atom.cpp
+++ b/src/fix_property_atom.cpp
@@ -222,7 +222,7 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf)
 
     if ((m = atom->map(itag)) >= 0) {
       for (j = 0; j < nvalue; j++) {
-        if (style[j] == MOLECULE) atom->molecule[m] = atoi(values[j+1]);
+        if (style[j] == MOLECULE) atom->molecule[m] = ATOTAGINT(values[j+1]);
         else if (style[j] == CHARGE) atom->q[m] = atof(values[j+1]);
         else if (style[j] == INTEGER)
           atom->ivector[index[j]][m] = atoi(values[j+1]);
@@ -284,7 +284,7 @@ void FixPropertyAtom::write_data_section_pack(int mth, double **buf)
   for (int m = 0; m < nvalue; m++) {
     int mp1 = m+1;
     if (style[m] == MOLECULE) {
-      int *molecule = atom->molecule;
+      tagint *molecule = atom->molecule;
       for (i = 0; i < nlocal; i++) buf[i][mp1] = ubuf(molecule[i]).d;
     } else if (style[m] == CHARGE) {
       double *q = atom->q;
@@ -327,7 +327,9 @@ void FixPropertyAtom::write_data_section(int mth, FILE *fp,
   for (int i = 0; i < n; i++) {
     fprintf(fp,TAGINT_FORMAT,(tagint) ubuf(buf[i][0]).i);
     for (m = 0; m < nvalue; m++) {
-      if (style[m] == MOLECULE || style[m] == INTEGER)
+      if (style[m] == MOLECULE)
+        fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[i][m+1]).i);
+      else if (style[m] == INTEGER)
         fprintf(fp," %d",(int) ubuf(buf[i][m+1]).i);
       else fprintf(fp," %g",buf[i][m+1]);
     }
@@ -343,7 +345,7 @@ double FixPropertyAtom::memory_usage()
 {
   double bytes = 0.0;
   for (int m = 0; m < nvalue; m++) {
-    if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(int);
+    if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(tagint);
     else if (style[m] == CHARGE) bytes = atom->nmax * sizeof(double);
     else if (style[m] == INTEGER) bytes = atom->nmax * sizeof(int);
     else if (style[m] == DOUBLE) bytes = atom->nmax * sizeof(double);
@@ -363,7 +365,7 @@ void FixPropertyAtom::grow_arrays(int nmax)
   for (int m = 0; m < nvalue; m++) {
     if (style[m] == MOLECULE) {
       memory->grow(atom->molecule,nmax,"atom:molecule");
-      size_t nbytes = (nmax-nmax_old) * sizeof(int);
+      size_t nbytes = (nmax-nmax_old) * sizeof(tagint);
       memset(&atom->molecule[nmax_old],0,nbytes);
     } else if (style[m] == CHARGE) {
       memory->grow(atom->q,nmax,"atom:q");
@@ -412,7 +414,7 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf)
   int m = 0;
   for (k = 0; k < nvalue; k++) {
     if (style[k] == MOLECULE) {
-      int *molecule = atom->molecule;
+      tagint *molecule = atom->molecule;
       for (i = 0; i < n; i++) {
         j = list[i];
         buf[m++] = ubuf(molecule[j]).d;
@@ -452,10 +454,10 @@ int FixPropertyAtom::unpack_border(int n, int first, double *buf)
   int m = 0;
   for (k = 0; k < nvalue; k++) {
     if (style[k] == MOLECULE) {
-      int *molecule = atom->molecule;
+      tagint *molecule = atom->molecule;
       last = first + n;
       for (i = first; i < last; i++)
-        molecule[i] = (int) ubuf(buf[m++]).i;
+        molecule[i] = (tagint) ubuf(buf[m++]).i;
     } else if (style[k] == CHARGE) {
       double *q = atom->q;
       last = first + n;
@@ -500,7 +502,7 @@ int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
 {
   for (int m = 0; m < nvalue; m++) {
     if (style[m] == MOLECULE) 
-      atom->molecule[nlocal] = (int) ubuf(buf[m]).i;
+      atom->molecule[nlocal] = (tagint) ubuf(buf[m]).i;
     else if (style[m] == CHARGE)
       atom->q[nlocal] = buf[m];
     else if (style[m] == INTEGER) 
@@ -546,7 +548,7 @@ void FixPropertyAtom::unpack_restart(int nlocal, int nth)
 
   for (int i = 0; i < nvalue; i++) {
     if (style[i] == MOLECULE) 
-      atom->molecule[nlocal] = (int) ubuf(extra[nlocal][m++]).i;
+      atom->molecule[nlocal] = (tagint) ubuf(extra[nlocal][m++]).i;
     else if (style[i] == CHARGE)
       atom->q[nlocal] = extra[nlocal][m++];
     else if (style[i] == INTEGER) 
diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp
index addcffcab477cb5332a9cfd657ca1938a588b12c..6863576884d3ae909dff8ecbb4f7298faccb01c0 100644
--- a/src/fix_store_state.cpp
+++ b/src/fix_store_state.cpp
@@ -603,7 +603,7 @@ void FixStoreState::pack_id(int n)
 
 void FixStoreState::pack_molecule(int n)
 {
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
diff --git a/src/group.cpp b/src/group.cpp
index 22b9035cc4eedfecc21438b1d0d5296a38229169..ce99834a6e569986d4243b34bb54d48883c126fa 100644
--- a/src/group.cpp
+++ b/src/group.cpp
@@ -199,11 +199,11 @@ void Group::assign(int narg, char **arg)
           bound2 = force->bnumeric(FLERR,arg[4]);
       } else if (narg != 4) error->all(FLERR,"Illegal group command");
 
-      tagint *tag = atom->tag;
-      int *attribute;
+      int *attribute = NULL;
+      tagint *tattribute = NULL;
       if (category == TYPE) attribute = atom->type;
-      else if (category == MOLECULE) attribute = atom->molecule;
-      else if (category == ID) attribute = NULL;
+      else if (category == MOLECULE) tattribute = atom->molecule;
+      else if (category == ID) tattribute = atom->tag;
 
       // add to group if meets condition
 
@@ -234,25 +234,25 @@ void Group::assign(int narg, char **arg)
       } else {
         if (condition == LT) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] < bound1) mask[i] |= bit;
+            if (tattribute[i] < bound1) mask[i] |= bit;
         } else if (condition == LE) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] <= bound1) mask[i] |= bit;
+            if (tattribute[i] <= bound1) mask[i] |= bit;
         } else if (condition == GT) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] > bound1) mask[i] |= bit;
+            if (tattribute[i] > bound1) mask[i] |= bit;
         } else if (condition == GE) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] >= bound1) mask[i] |= bit;
+            if (tattribute[i] >= bound1) mask[i] |= bit;
         } else if (condition == EQ) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] == bound1) mask[i] |= bit;
+            if (tattribute[i] == bound1) mask[i] |= bit;
         } else if (condition == NEQ) {
           for (i = 0; i < nlocal; i++) 
-            if (tag[i] != bound1) mask[i] |= bit;
+            if (tattribute[i] != bound1) mask[i] |= bit;
         } else if (condition == BETWEEN) {
           for (i = 0; i < nlocal; i++)
-            if (tag[i] >= bound1 && tag[i] <= bound2)
+            if (tattribute[i] >= bound1 && tattribute[i] <= bound2)
               mask[i] |= bit;
         }
       }
@@ -260,38 +260,25 @@ void Group::assign(int narg, char **arg)
     // args = list of values
       
     } else {
-      tagint *tag = atom->tag;
-      int *attribute;
+      int *attribute = NULL;
+      tagint *tattribute = NULL;
       if (category == TYPE) attribute = atom->type;
-      else if (category == MOLECULE) attribute = atom->molecule;
-      else if (category == ID) attribute = NULL;
+      else if (category == MOLECULE) tattribute = atom->molecule;
+      else if (category == ID) tattribute = atom->tag;
                                  
       char *ptr;
       tagint start,stop,delta;
 
       for (int iarg = 2; iarg < narg; iarg++) {
-        if (attribute) {
-          if (strchr(arg[iarg],':')) {
-            start = atoi(strtok(arg[iarg],":")); 
-            stop = atoi(strtok(NULL,":"));
-            ptr = strtok(NULL,":");
-            if (ptr) delta = atoi(ptr);
-            else delta = 1;
-          } else {
-            start = stop = atoi(arg[iarg]);
-            delta = 1;
-          }
+        if (strchr(arg[iarg],':')) {
+          start = ATOTAGINT(strtok(arg[iarg],":")); 
+          stop = ATOTAGINT(strtok(NULL,":"));
+          ptr = strtok(NULL,":");
+          if (ptr) delta = ATOTAGINT(ptr);
+          else delta = 1;
         } else {
-          if (strchr(arg[iarg],':')) {
-            start = ATOTAGINT(strtok(arg[iarg],":")); 
-            stop = ATOTAGINT(strtok(NULL,":"));
-            ptr = strtok(NULL,":");
-            if (ptr) delta = ATOTAGINT(ptr);
-            else delta = 1;
-          } else {
-            start = stop = ATOTAGINT(arg[iarg]);
-            delta = 1;
-          }
+          start = stop = ATOTAGINT(arg[iarg]);
+          delta = 1;
         }
 
         // add to group if attribute matches value or sequence
@@ -302,8 +289,8 @@ void Group::assign(int narg, char **arg)
                 (attribute[i]-start) % delta == 0) mask[i] |= bit;
         } else {
           for (i = 0; i < nlocal; i++)
-            if (tag[i] >= start && tag[i] <= stop &&
-                (tag[i]-start) % delta == 0) mask[i] |= bit;
+            if (tattribute[i] >= start && tattribute[i] <= stop &&
+                (tattribute[i]-start) % delta == 0) mask[i] |= bit;
         }
       }
     }
diff --git a/src/lammps.cpp b/src/lammps.cpp
index a896c752d2564f521a6a129dad9d125ca2d28f77..5bd15ace425d640b52682a9bfef3ba0bad5f0a60 100644
--- a/src/lammps.cpp
+++ b/src/lammps.cpp
@@ -366,7 +366,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
     }
   }
 
-  // check datatype settings in lmptype.h
+  // check consistency of datatype settings in lmptype.h
 
   if (sizeof(smallint) != sizeof(int))
     error->all(FLERR,"Smallint setting in lmptype.h is invalid");
@@ -378,6 +378,10 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
     error->all(FLERR,"Bigint setting in lmptype.h is invalid");
 
   int mpisize;
+  MPI_Type_size(MPI_LMP_TAGINT,&mpisize);
+  if (mpisize != sizeof(tagint))
+      error->all(FLERR,"MPI_LMP_TAGINT and tagint in "
+                 "lmptype.h are not compatible");
   MPI_Type_size(MPI_LMP_BIGINT,&mpisize);
   if (mpisize != sizeof(bigint))
       error->all(FLERR,"MPI_LMP_BIGINT and bigint in "
diff --git a/src/lmptype.h b/src/lmptype.h
index a23aa9113b13cfdd46110ee6ec633bb6c43e5952..29301058b8cb44585d85c533753efd495a40002e 100644
--- a/src/lmptype.h
+++ b/src/lmptype.h
@@ -15,13 +15,13 @@
 
 // smallint = variables for on-procesor system (nlocal, nmax, etc)
 // imageint = variables for atom image flags (image)
-// tagint = variables for atom IDs (tag)
+// tagint = variables for atom IDs and molecule IDs (tag,molecule)
 // bigint = variables for total system (natoms, ntimestep, etc)
 
 // smallint must be an int, as defined by C compiler
 // imageint can be 32-bit or 64-bit int, must be >= smallint
 // tagint can be 32-bit or 64-bit int, must be >= smallint
-// bigint can be 32-bit or 64-bit int, must be >= tagint
+// bigint can be 32-bit or 64-bit int, must be >= imageint,tagint
 
 // MPI_LMP_BIGINT = MPI data type corresponding to a bigint
 
diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp
index 7b536f469991dde548abf30be86acce65653aa60..6df4bfb5129959bd6a6e4103d5d50bee649984d3 100644
--- a/src/neigh_full.cpp
+++ b/src/neigh_full.cpp
@@ -39,7 +39,7 @@ void Neighbor::full_nsq(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -122,7 +122,7 @@ void Neighbor::full_nsq_ghost(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -219,7 +219,7 @@ void Neighbor::full_bin(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
@@ -310,7 +310,7 @@ void Neighbor::full_bin_ghost(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -431,7 +431,7 @@ void Neighbor::full_multi(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
diff --git a/src/neigh_gran.cpp b/src/neigh_gran.cpp
index 465e6af863fdf20b07eed6860e4df5ea2cf30b6d..e6ef93c5ebba2c878d29b7c68fabffd4f22851f2 100644
--- a/src/neigh_gran.cpp
+++ b/src/neigh_gran.cpp
@@ -50,7 +50,7 @@ void Neighbor::granular_nsq_no_newton(NeighList *list)
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   if (includegroup) {
@@ -178,7 +178,7 @@ void Neighbor::granular_nsq_newton(NeighList *list)
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   if (includegroup) {
@@ -284,7 +284,7 @@ void Neighbor::granular_bin_no_newton(NeighList *list)
   tagint *tag = atom->tag;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
@@ -420,7 +420,7 @@ void Neighbor::granular_bin_newton(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
@@ -522,7 +522,7 @@ void Neighbor::granular_bin_newton_tri(NeighList *list)
   double *radius = atom->radius;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
 
diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp
index 75736e97d28ce31b300440ba1a7c46ef8eb09d33..68f7832fa2a641ce9f73e1055d64ff35df6d7689 100644
--- a/src/neigh_half_bin.cpp
+++ b/src/neigh_half_bin.cpp
@@ -46,7 +46,7 @@ void Neighbor::half_bin_no_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
   int molecular = atom->molecular;
@@ -141,7 +141,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -263,7 +263,7 @@ void Neighbor::half_bin_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
   int molecular = atom->molecular;
@@ -382,7 +382,7 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   if (includegroup) nlocal = atom->nfirst;
   int molecular = atom->molecular;
diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp
index fddc683fb7c8d9fbd7e04dd06fea4b86cb3475ec..6dc03a43efaef34672416472c5f9a80709bec5d4 100644
--- a/src/neigh_half_multi.cpp
+++ b/src/neigh_half_multi.cpp
@@ -48,7 +48,7 @@ void Neighbor::half_multi_no_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
@@ -147,7 +147,7 @@ void Neighbor::half_multi_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
@@ -274,7 +274,7 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp
index f8fad972bbd96d85db547cc8eeb4d8b716149f77..efdd26a311efd90a1580300372b50f328769dcea 100644
--- a/src/neigh_half_nsq.cpp
+++ b/src/neigh_half_nsq.cpp
@@ -39,7 +39,7 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -123,7 +123,7 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -228,7 +228,7 @@ void Neighbor::half_nsq_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp
index 738426abdb2e2ac5f8be83760d8347f9288bfcc6..b2d19b11e07f14d007324ed6b15ef1c1058227a3 100644
--- a/src/neigh_respa.cpp
+++ b/src/neigh_respa.cpp
@@ -43,7 +43,7 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -186,7 +186,7 @@ void Neighbor::respa_nsq_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int molecular = atom->molecular;
@@ -351,7 +351,7 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
@@ -504,7 +504,7 @@ void Neighbor::respa_bin_newton(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
@@ -698,7 +698,7 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
   double **x = atom->x;
   int *type = atom->type;
   int *mask = atom->mask;
-  int *molecule = atom->molecule;
+  tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
   if (includegroup) nlocal = atom->nfirst;
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 10f7ce0a105f9de6712807020392e13dd1bdc6a4..15ac28378b4407139a9519eecf8b4fec1dc7128c 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1929,7 +1929,7 @@ int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz)
 ------------------------------------------------------------------------- */
 
 int Neighbor::exclusion(int i, int j, int itype, int jtype,
-                        int *mask, int *molecule) const {
+                        int *mask, tagint *molecule) const {
   int m;
 
   if (nex_type && ex_type[itype][jtype]) return 1;
diff --git a/src/neighbor.h b/src/neighbor.h
index 751e15523165ec6d435a17db248b2166ca6609be..2bf33d22e15b3f71af10757d417cedc25b95958a 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -171,7 +171,7 @@ class Neighbor : protected Pointers {
   int coord2bin(double *, int &, int &, int&); // ditto
 
   int exclusion(int, int, int,
-                int, int *, int *) const;  // test for pair exclusion
+                int, int *, tagint *) const;    // test for pair exclusion
 
   virtual void choose_build(int, class NeighRequest *);
   void choose_stencil(int, class NeighRequest *);
diff --git a/src/replicate.cpp b/src/replicate.cpp
index 9b4e9bcf4852e062b75a4eeeb0e75187586421f5..b31506d6c1be0075eab029e433af9268cad72672 100644
--- a/src/replicate.cpp
+++ b/src/replicate.cpp
@@ -83,11 +83,11 @@ void Replicate::command(int narg, char **arg)
 
   // maxmol = largest molecule tag across all existing atoms
 
-  int maxmol = 0;
+  tagint maxmol = 0;
   if (atom->molecule_flag) {
     for (i = 0; i < atom->nlocal; i++) maxmol = MAX(atom->molecule[i],maxmol);
-    int maxmol_all;
-    MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_INT,MPI_MAX,world);
+    tagint maxmol_all;
+    MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
     maxmol = maxmol_all;
   }
 
@@ -254,8 +254,8 @@ void Replicate::command(int narg, char **arg)
   AtomVec *old_avec = old->avec;
   AtomVec *avec = atom->avec;
 
-  int ix,iy,iz,mol_offset;
-  tagint atom_offset;
+  int ix,iy,iz;
+  tagint atom_offset,mol_offset;
   imageint image;
   double x[3],lamda[3];
   double *coord;
diff --git a/src/set.cpp b/src/set.cpp
index 562b8f94cedb0b5a0688175379f0abf37beb3074..108aecc2391df785bb38498d2f9a36edd0e2ee1e 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -451,11 +451,12 @@ void Set::selection(int n)
   } else if (style == MOL_SELECT) {
     if (atom->molecule_flag == 0)
       error->all(FLERR,"Cannot use set mol with no molecule IDs defined");
-    else force->bounds(id,BIG,nlo,nhi,0);
+    bigint nlobig,nhibig;
+    force->boundsbig(id,MAXTAGINT,nlobig,nhibig);
 
-    int *molecule = atom->molecule;
+    tagint *molecule = atom->molecule;
     for (int i = 0; i < n; i++)
-      if (molecule[i] >= nlo && molecule[i] <= nhi) select[i] = 1;
+      if (molecule[i] >= nlobig && molecule[i] <= nhibig) select[i] = 1;
       else select[i] = 0;
 
   } else if (style == TYPE_SELECT) {
diff --git a/src/variable.cpp b/src/variable.cpp
index 87fed075c1839b79125a9499781a0f5477460275..c866a02ad79a2a27cfd378206aa97bd6e00be8cc 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -3521,6 +3521,7 @@ void Variable::atom_vector(char *word, Tree **tree,
       newtree->barray = (bigint *) atom->tag;
     }
     newtree->nstride = 1;
+
   } else if (strcmp(word,"mass") == 0) {
     if (atom->rmass) {
       newtree->nstride = 1;
@@ -3529,17 +3530,25 @@ void Variable::atom_vector(char *word, Tree **tree,
       newtree->type = TYPEARRAY;
       newtree->array = atom->mass;
     }
+
   } else if (strcmp(word,"type") == 0) {
     newtree->type = INTARRAY;
     newtree->nstride = 1;
     newtree->iarray = atom->type;
+
   } else if (strcmp(word,"mol") == 0) {
     if (!atom->molecule_flag) 
       error->one(FLERR,"Variable uses atom property that isn't allocated");
-    newtree->type = INTARRAY;
+    if (sizeof(tagint) == sizeof(smallint)) {
+      newtree->type = INTARRAY;
+      newtree->iarray = atom->molecule;
+    } else {
+      newtree->type = BIGINTARRAY;
+      newtree->barray = (bigint *) atom->molecule;
+    }
     newtree->nstride = 1;
-    newtree->iarray = atom->molecule;
   }
+
   else if (strcmp(word,"x") == 0) newtree->array = &atom->x[0][0];
   else if (strcmp(word,"y") == 0) newtree->array = &atom->x[0][1];
   else if (strcmp(word,"z") == 0) newtree->array = &atom->x[0][2];