From 127eb167c1162552294dfc60b60138615bb79dea Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 22 Jan 2014 17:44:32 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11293
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/GRANULAR/fix_pour.cpp              |  4 +-
 src/GRANULAR/fix_pour.h                |  3 +-
 src/MC/fix_bond_swap.cpp               |  5 +-
 src/MC/fix_gcmc.cpp                    | 65 +++++++++++++++-----------
 src/MC/fix_gcmc.h                      |  7 +--
 src/MISC/fix_deposit.cpp               |  4 +-
 src/MISC/fix_deposit.h                 |  3 +-
 src/MISC/fix_evaporate.cpp             | 13 +++---
 src/MOLECULE/atom_vec_angle.cpp        | 21 +++++----
 src/MOLECULE/atom_vec_bond.cpp         | 23 ++++-----
 src/MOLECULE/atom_vec_full.cpp         | 24 +++++-----
 src/MOLECULE/atom_vec_molecular.cpp    | 21 +++++----
 src/MOLECULE/atom_vec_molecular.h      |  2 +-
 src/POEMS/fix_poems.cpp                | 38 ++++++++-------
 src/RIGID/fix_rigid.cpp                | 16 ++++---
 src/RIGID/fix_rigid_small.cpp          | 40 ++++++++--------
 src/RIGID/fix_rigid_small.h            |  2 +-
 src/USER-LB/fix_lb_rigid_pc_sphere.cpp | 36 ++++++++------
 src/USER-OMP/neigh_full_omp.cpp        | 10 ++--
 src/USER-OMP/neigh_gran_omp.cpp        | 10 ++--
 src/USER-OMP/neigh_half_bin_omp.cpp    |  8 ++--
 src/USER-OMP/neigh_half_multi_omp.cpp  |  6 +--
 src/USER-OMP/neigh_half_nsq_omp.cpp    |  6 +--
 src/USER-OMP/neigh_respa_omp.cpp       | 10 ++--
 src/USER-REAXC/fix_reaxc_bonds.cpp     |  4 +-
 src/atom.h                             |  3 +-
 src/atom_vec_line.cpp                  | 25 +++++-----
 src/atom_vec_tri.cpp                   | 23 ++++-----
 src/compute.cpp                        | 20 ++++----
 src/compute.h                          |  2 +-
 src/compute_atom_molecule.h            |  2 +-
 src/compute_com_molecule.cpp           | 20 ++++----
 src/compute_com_molecule.h             |  2 +-
 src/compute_gyration_molecule.cpp      | 42 +++++++++--------
 src/compute_gyration_molecule.h        |  2 +-
 src/compute_inertia_molecule.cpp       | 19 ++++----
 src/compute_inertia_molecule.h         |  2 +-
 src/compute_msd_molecule.cpp           | 26 +++++------
 src/compute_msd_molecule.h             |  2 +-
 src/compute_property_molecule.cpp      | 12 +++--
 src/compute_property_molecule.h        |  2 +-
 src/create_atoms.cpp                   | 10 ++--
 src/delete_atoms.cpp                   | 16 +++----
 src/dump_custom.cpp                    |  4 +-
 src/fix_property_atom.cpp              | 22 +++++----
 src/fix_store_state.cpp                |  2 +-
 src/group.cpp                          | 63 ++++++++++---------------
 src/lammps.cpp                         |  6 ++-
 src/lmptype.h                          |  4 +-
 src/neigh_full.cpp                     | 10 ++--
 src/neigh_gran.cpp                     | 10 ++--
 src/neigh_half_bin.cpp                 |  8 ++--
 src/neigh_half_multi.cpp               |  6 +--
 src/neigh_half_nsq.cpp                 |  6 +--
 src/neigh_respa.cpp                    | 10 ++--
 src/neighbor.cpp                       |  2 +-
 src/neighbor.h                         |  2 +-
 src/replicate.cpp                      | 10 ++--
 src/set.cpp                            |  7 +--
 src/variable.cpp                       | 13 +++++-
 60 files changed, 421 insertions(+), 375 deletions(-)

diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 9a439d87d6..4dc0c5cbef 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 34d13da985..e57379906e 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 31bdb75682..79b298391a 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 ab85e6e21b..90202ad7eb 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 7fdd1f70c6..ee532820b8 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 ea8fb37be9..f83af6d7d6 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 56bf34c6dd..0c43a423d9 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 ab153f805e..421c22e666 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 9b303451c3..bafac2391f 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 7e5e639918..030a86ac86 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 09175dd034..f132fdc599 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 2f2265b5df..766da9e354 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 58aa039e1e..293b44a826 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 1aea5795a1..cbada23c5a 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 eee93e95c6..c79fa2359e 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 46809d1b34..78ef0ed0f2 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 6f73a9c03a..392a928bee 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 a0e3016698..6a1e4e531e 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 e609627897..d1230555ce 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 5f0402f363..e293102e9a 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 f501684a19..76f36ea965 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 b46594870f..a804e8008f 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 a98c4721f3..fa90945d55 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 11c7fd9ae3..458b2bbac4 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 3ced20a289..8d8818a49b 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 8a5a0b643f..c552cad6c2 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 5c89d237fc..bd5a5ea2b3 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 d170c8b5d0..20528ff0d6 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 9be14029dc..f9a95e8d10 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 e7cb31bae1..3ae91e5b1a 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 769f58d5ec..cd135e47d8 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 130072cc63..74b2cb5ac3 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 a4f941c52a..3e6fdedd99 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 bf4a70a296..f81a373830 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 99910782fa..78b49231a2 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 8cdb3cb64c..df9958a04a 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 80657108e7..a081dd1ad5 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 736ac9cfc6..515436ba3a 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 b4c5bd53a1..839698cabd 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 c0bcafc9ea..dafe2db680 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 0c5ae68df0..27cebe8fe1 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 1e1f5716da..40c7fe4ec4 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 ac1c354d29..5c353c7fe5 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 f0c0371044..9619214b94 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 a1be2b556b..115750e271 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 addcffcab4..6863576884 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 22b9035cc4..ce99834a6e 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 a896c752d2..5bd15ace42 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 a23aa9113b..29301058b8 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 7b536f4699..6df4bfb512 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 465e6af863..e6ef93c5eb 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 75736e97d2..68f7832fa2 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 fddc683fb7..6dc03a43ef 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 f8fad972bb..efdd26a311 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 738426abdb..b2d19b11e0 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 10f7ce0a10..15ac28378b 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 751e155231..2bf33d22e1 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 9b4e9bcf48..b31506d6c1 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 562b8f94ce..108aecc239 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 87fed075c1..c866a02ad7 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];
-- 
GitLab