diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index f88753297e..c799816fd3 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -21,6 +21,7 @@
 #include "stdlib.h"
 #include "math.h"
 #include "pppm.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "comm.h"
 #include "neighbor.h"
@@ -854,7 +855,7 @@ void PPPM::set_grid()
   acons[7][6] = 4887769399.0 / 37838389248.0;
 
   double q2 = qsqsum / force->dielectric;
-  double natoms = atom->natoms;
+  bigint natoms = atom->natoms;
 
   // use xprd,yprd,zprd even if triclinic so grid size is the same
   // adjust z dimension for 2d slab PPPM
@@ -1006,7 +1007,7 @@ int PPPM::factorable(int n)
    compute RMS precision for a dimension
 ------------------------------------------------------------------------- */
 
-double PPPM::rms(double h, double prd, double natoms,
+double PPPM::rms(double h, double prd, bigint natoms,
 		 double q2, double **acons)
 {
   double sum = 0.0;
@@ -1027,7 +1028,7 @@ double PPPM::diffpr(double hx, double hy, double hz, double q2, double **acons)
   double xprd = domain->xprd;
   double yprd = domain->yprd;
   double zprd = domain->zprd;
-  double natoms = atom->natoms;
+  bigint natoms = atom->natoms;
 
   lprx = rms(hx,xprd,natoms,q2,acons);
   lpry = rms(hy,yprd,natoms,q2,acons);
diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h
index 5778bda703..9ca04a7540 100644
--- a/src/KSPACE/pppm.h
+++ b/src/KSPACE/pppm.h
@@ -21,6 +21,7 @@ KSpaceStyle(pppm,PPPM)
 #define LMP_PPPM_H
 
 #include "kspace.h"
+#include "lmptype.h"
 
 namespace LAMMPS_NS {
 
@@ -83,7 +84,7 @@ class PPPM : public KSpace {
   void allocate();
   void deallocate();
   int factorable(int);
-  double rms(double, double, double, double, double **);
+  double rms(double, double, bigint, double, double **);
   double diffpr(double, double, double, double, double **);
   void compute_gf_denom();
   double gf_denom(double, double, double);
diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp
index 8dfecad9f7..8e7cb00497 100644
--- a/src/REPLICA/fix_neb.cpp
+++ b/src/REPLICA/fix_neb.cpp
@@ -110,10 +110,10 @@ void FixNEB::init()
   memory->destroy_2d_double_array(xprev);
   memory->destroy_2d_double_array(xnext);
   memory->destroy_2d_double_array(tangent);
-  natoms = atom->nlocal;
-  xprev = memory->create_2d_double_array(natoms,3,"neb:xprev");
-  xnext = memory->create_2d_double_array(natoms,3,"neb:xnext");
-  tangent = memory->create_2d_double_array(natoms,3,"neb:tangent");
+  nebatoms = atom->nlocal;
+  xprev = memory->create_2d_double_array(nebatoms,3,"neb:xprev");
+  xnext = memory->create_2d_double_array(nebatoms,3,"neb:xnext");
+  tangent = memory->create_2d_double_array(nebatoms,3,"neb:tangent");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -167,7 +167,7 @@ void FixNEB::min_post_force(int vflag)
   double **x = atom->x;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
-  if (nlocal != natoms) error->one("Atom count changed in fix neb");
+  if (nlocal != nebatoms) error->one("Atom count changed in fix neb");
 
   if (ireplica == 0) 
     MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);
diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h
index 066259b166..06b2c24ca6 100644
--- a/src/REPLICA/fix_neb.h
+++ b/src/REPLICA/fix_neb.h
@@ -46,7 +46,7 @@ class FixNEB : public Fix {
   char *id_pe;
   class Compute *pe;
 
-  int natoms;
+  int nebatoms;
   double **xprev,**xnext;
   double **tangent;
 };
diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp
index 96539381cc..7864eb28a9 100644
--- a/src/REPLICA/prd.cpp
+++ b/src/REPLICA/prd.cpp
@@ -20,6 +20,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "prd.h"
+#include "lmptype.h"
 #include "universe.h"
 #include "update.h"
 #include "atom.h"
@@ -47,8 +48,6 @@
 
 using namespace LAMMPS_NS;
 
-#define MAXINT 0x7FFFFFFF
-
 /* ---------------------------------------------------------------------- */
 
 PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {}
@@ -117,7 +116,7 @@ void PRD::command(int narg, char **arg)
 
   // workspace for inter-replica communication via gathers
 
-  natoms = static_cast<int> (atom->natoms);
+  natoms = atom->natoms;
 
   displacements = NULL;
   tagall = NULL;
@@ -432,7 +431,7 @@ void PRD::dephase()
   timer->barrier_start(TIME_LOOP);
 
   for (int i = 0; i < n_dephase; i++) {
-    int seed = static_cast<int> (random_dephase->uniform() * MAXINT);
+    int seed = static_cast<int> (random_dephase->uniform() * MAXINT32);
     if (seed == 0) seed = 1;
     velocity->create(temp_dephase,seed);
     update->integrate->run(t_dephase);
diff --git a/src/atom.cpp b/src/atom.cpp
index 16f21748f0..6224725143 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -45,7 +45,8 @@ using namespace LAMMPS_NS;
 
 Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
 {
-  natoms = nlocal = nghost = nmax = 0;
+  natoms = 0;
+  nlocal = nghost = nmax = 0;
   ntypes = 0;
   nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
   nbonds = nangles = ndihedrals = nimpropers = 0;
diff --git a/src/atom.h b/src/atom.h
index b12b337062..6c3f84d969 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -15,6 +15,7 @@
 #define LMP_ATOM_H
 
 #include "pointers.h"
+#include "lmptype.h"
 
 namespace LAMMPS_NS {
 
@@ -25,7 +26,7 @@ class Atom : protected Pointers {
 
   // atom counts
 
-  double natoms;                // total # of atoms in system, could be 0
+  bigint natoms;                // total # of atoms in system, could be 0
   int nlocal,nghost;            // # of owned and ghost atoms on this proc
   int nmax;                     // max # of owned+ghost in arrays on this proc
   int tag_enable;               // 0/1 if atom ID tags are defined
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index 3d7ac95a73..5f54d8fe80 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -15,6 +15,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "create_atoms.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "comm.h"
@@ -28,7 +29,6 @@
 
 using namespace LAMMPS_NS;
 
-#define MAXATOMS 0x7FFFFFFF
 #define BIG      1.0e30
 #define EPSILON  1.0e-6
 
@@ -145,7 +145,7 @@ void CreateAtoms::command(int narg, char **arg)
 
   // add atoms
 
-  double natoms_previous = atom->natoms;
+  bigint natoms_previous = atom->natoms;
   int nlocal_previous = atom->nlocal;
 
   if (style == SINGLE) add_single();
@@ -168,16 +168,16 @@ void CreateAtoms::command(int narg, char **arg)
 
   // new total # of atoms
 
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&atom->natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
 
   // print status
 
   if (comm->me == 0) {
     if (screen)
-      fprintf(screen,"Created %.15g atoms\n",atom->natoms-natoms_previous);
+      fprintf(screen,"Created %lu atoms\n",atom->natoms-natoms_previous);
     if (logfile)
-      fprintf(logfile,"Created %.15g atoms\n",atom->natoms-natoms_previous);
+      fprintf(logfile,"Created %lu atoms\n",atom->natoms-natoms_previous);
   }
 
   // reset simulation now that more atoms are defined
@@ -185,8 +185,8 @@ void CreateAtoms::command(int narg, char **arg)
   // if global map exists, reset it
   // if a molecular system, set nspecial to 0 for new atoms
 
-  if (atom->natoms > MAXATOMS) atom->tag_enable = 0;
-  if (atom->natoms <= MAXATOMS) atom->tag_extend();
+  if (atom->natoms > MAXINT32) atom->tag_enable = 0;
+  if (atom->natoms <= MAXINT32) atom->tag_extend();
 
   if (atom->map_style) {
     atom->nghost = 0;
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 8ccd3d0287..b2fec4e2a7 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -14,6 +14,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "delete_atoms.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "comm.h"
@@ -49,7 +50,7 @@ void DeleteAtoms::command(int narg, char **arg)
 
   // store state before delete
 
-  int natoms_previous = static_cast<int> (atom->natoms);
+  bigint natoms_previous = atom->natoms;
 
   // delete the atoms
 
@@ -91,8 +92,8 @@ void DeleteAtoms::command(int narg, char **arg)
   // reset atom->map if it exists
   // set nghost to 0 so old ghosts of deleted atoms won't be mapped
 
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&atom->natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
   if (atom->map_style) {
     atom->nghost = 0;
     atom->map_init();
@@ -101,12 +102,12 @@ void DeleteAtoms::command(int narg, char **arg)
 
   // print before and after atom count
 
-  int ndelete = static_cast<int> (natoms_previous - atom->natoms);
+  bigint ndelete = natoms_previous - atom->natoms;
 
   if (comm->me == 0) {
-    if (screen) fprintf(screen,"Deleted %d atoms, new total = %.15g\n",
+    if (screen) fprintf(screen,"Deleted %lu atoms, new total = %lu\n",
 			ndelete,atom->natoms);
-    if (logfile) fprintf(logfile,"Deleted %d atoms, new total = %.15g\n",
+    if (logfile) fprintf(logfile,"Deleted %lu atoms, new total = %lu\n",
 			 ndelete,atom->natoms);
   }
 }
diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp
index 29ffe56c4f..bd9171106c 100644
--- a/src/displace_atoms.cpp
+++ b/src/displace_atoms.cpp
@@ -15,6 +15,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "displace_atoms.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "modify.h"
 #include "domain.h"
@@ -210,12 +211,12 @@ void DisplaceAtoms::command(int narg, char **arg)
 
   // check if any atoms were lost
 
-  double natoms;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint natoms;
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
   if (natoms != atom->natoms) {
     char str[128];
-    sprintf(str,"Lost atoms via displace_atoms: original %.15g current %.15g",
+    sprintf(str,"Lost atoms via displace_atoms: original %lu current %lu",
 	    atom->natoms,natoms);
     error->all(str);
   }
diff --git a/src/displace_box.cpp b/src/displace_box.cpp
index 58962b795d..a8b5b351f8 100644
--- a/src/displace_box.cpp
+++ b/src/displace_box.cpp
@@ -16,6 +16,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "displace_box.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "modify.h"
 #include "domain.h"
@@ -376,12 +377,12 @@ void DisplaceBox::command(int narg, char **arg)
 
   // check if any atoms were lost
 
-  double natoms;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint natoms;
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
   if (natoms != atom->natoms) {
     char str[128];
-    sprintf(str,"Lost atoms via displace_box: original %.15g current %.15g",
+    sprintf(str,"Lost atoms via displace_box: original %lu current %lu",
 	    atom->natoms,natoms);
     error->all(str);
   }
diff --git a/src/dump.cpp b/src/dump.cpp
index a5d64b5452..02e316caf2 100644
--- a/src/dump.cpp
+++ b/src/dump.cpp
@@ -16,6 +16,7 @@
 #include "string.h"
 #include "stdio.h"
 #include "dump.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "irregular.h"
 #include "update.h"
@@ -31,7 +32,6 @@ using namespace LAMMPS_NS;
 
 Dump *Dump::dumpptr;
 
-#define MAXATOMS 0x7FFFFFFF
 #define BIG 1.0e20
 #define IBIG 2147483647
 #define EPSILON 1.0e-6
@@ -174,7 +174,7 @@ void Dump::init()
       irregular = new Irregular(lmp);
 
     double size = group->count(igroup);
-    if (size > MAXATOMS) error->all("Too many atoms to dump sort");
+    if (size > MAXINT32) error->all("Too many atoms to dump sort");
 
     // set reorderflag = 1 if can simply reorder local atoms rather than sort
     // criteria: sorting by ID, atom IDs are consecutive from 1 to Natoms
diff --git a/src/dump_dcd.cpp b/src/dump_dcd.cpp
index 6794cfbc64..0ae2b707ea 100644
--- a/src/dump_dcd.cpp
+++ b/src/dump_dcd.cpp
@@ -1,366 +1,366 @@
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under 
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-/* ----------------------------------------------------------------------
-   Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
-                        Axel Kohlmeyer (Temple U), support for groups
-------------------------------------------------------------------------- */
-
-#include "math.h"
-#include "inttypes.h"
-#include "stdio.h"
-#include "time.h"
-#include "string.h"
-#include "dump_dcd.h"
-#include "domain.h"
-#include "atom.h"
-#include "update.h"
-#include "output.h"
-#include "group.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-#define NFILE_POS 8L
-#define NSTEP_POS 20L
-
-// necessary to set SEEK params b/c MPI-2 messes with these settings
-
-#ifndef SEEK_SET
-#define SEEK_SET        0
-#define SEEK_CUR        1
-#define SEEK_END        2
-#endif
-
-/* ---------------------------------------------------------------------- */
-
-static inline void fwrite_int32(FILE* fd, uint32_t i)
-{
-  fwrite(&i,sizeof(uint32_t),1,fd);
-}
-
-/* ---------------------------------------------------------------------- */
-
-DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
-{
-  if (narg != 5) error->all("Illegal dump dcd command");
-  if (binary || compressed || multifile || multiproc)
-    error->all("Invalid dump dcd filename");
-
-  size_one = 3;
-  sort_flag = 1;
-  sortcol = 0;
-
-  unwrap_flag = 0;
-  format_default = NULL;
-    
-  // allocate global array for atom coords
-
-  if (igroup == 0) natoms = static_cast<int> (atom->natoms);
-  else natoms = static_cast<int> (group->count(igroup));
-  if (natoms <= 0) error->all("Invalid natoms for dump dcd");
-
-  coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
-  xf = &coords[0*natoms];
-  yf = &coords[1*natoms];
-  zf = &coords[2*natoms];
-
-  openfile();
-  headerflag = 0;
-  nevery_save = 0;
-  ntotal = 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-DumpDCD::~DumpDCD()
-{
-  memory->sfree(coords);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::init_style()
-{
-  if (sort_flag == 0 || sortcol != 0)
-    error->all("Dump dcd requires sorting by atom ID");
-
-  // check that dump frequency has not changed and is not a variable
-
-  int idump;
-  for (idump = 0; idump < output->ndump; idump++)
-    if (strcmp(id,output->dump[idump]->id) == 0) break;
-  if (output->every_dump[idump] == 0)
-    error->all("Cannot use variable every setting for dump dcd");
-
-  if (nevery_save == 0) nevery_save = output->every_dump[idump];
-  else if (nevery_save != output->every_dump[idump])
-    error->all("Cannot change dump_modify every for dump dcd");
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::openfile()
-{
-  if (me == 0) {
-    fp = fopen(filename,"wb");
-    if (fp == NULL) error->one("Cannot open dump file");
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::write_header(int n)
-{
-  if (n != natoms) error->all("Dump dcd of non-matching # of atoms");
-
-  // first time, write header for entire file
-
-  if (headerflag == 0) {
-    if (me == 0) write_dcd_header("Written by LAMMPS");
-    headerflag = 1;
-    nframes = 0;
-  }
-
-  // dim[] = size and angle cosines of orthogonal or triclinic box
-  // dim[0] = a = length of unit cell vector along x-axis
-  // dim[1] = gamma = cosine of angle between a and b
-  // dim[2] = b = length of unit cell vector in xy-plane
-  // dim[3] = beta = cosine of angle between a and c
-  // dim[4] = alpha = cosine of angle between b and c
-  // dim[5] = c = length of final unit cell vector
-  // 48 = 6 doubles
-
-  double dim[6];
-  if (domain->triclinic) {
-    double *h = domain->h;
-    double alen = h[0];
-    double blen = sqrt(h[5]*h[5] + h[1]*h[1]);
-    double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]);
-    dim[0] = alen;
-    dim[2] = blen;
-    dim[5] = clen;
-    dim[4] = (h[5]*h[4] + h[1]*h[3]) / blen/clen; // alpha
-    dim[3] = (h[0]*h[4]) / alen/clen;             // beta
-    dim[1] = (h[0]*h[5]) / alen/blen;             // gamma
-  } else {
-    dim[0] = domain->xprd;
-    dim[2] = domain->yprd;
-    dim[5] = domain->zprd;
-    dim[1] = dim[3] = dim[4] = 0.0;
-  }
-  
-  if (me == 0) {
-    uint32_t out_integer = 48;
-    fwrite_int32(fp,out_integer);
-    fwrite(dim,out_integer,1,fp); 
-    fwrite_int32(fp,out_integer);
-    if (flush_flag) fflush(fp);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int DumpDCD::count()
-{
-  if (igroup == 0) return atom->nlocal;
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  int m = 0;
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) m++;
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::pack(int *ids)
-{
-  int m,n;
-
-  int *tag = atom->tag;
-  double **x = atom->x;
-  int *image = atom->image;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  m = n = 0;
-  if (unwrap_flag) {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
-    double xy = domain->xy;
-    double xz = domain->xz;
-    double yz = domain->yz;
-
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-	int ix = (image[i] & 1023) - 512;
-	int iy = (image[i] >> 10 & 1023) - 512;
-	int iz = (image[i] >> 20) - 512;
-
-	if (domain->triclinic) {
-	  buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz;
-	  buf[m++] = x[i][1] + iy * yprd + iz * yz;
-	  buf[m++] = x[i][2] + iz * zprd;
-	} else {
-	  buf[m++] = x[i][0] + ix * xprd;
-	  buf[m++] = x[i][1] + iy * yprd;
-	  buf[m++] = x[i][2] + iz * zprd;
-	}
-	ids[n++] = tag[i];
-      }
-    }
-
-  } else {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	buf[m++] = x[i][0];
-	buf[m++] = x[i][1];
-	buf[m++] = x[i][2];
-	ids[n++] = tag[i];
-      }
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::write_data(int n, double *mybuf)
-{
-  // copy buf atom coords into 3 global arrays
-
-  int m = 0;
-  for (int i = 0; i < n; i++) {
-    xf[ntotal] = mybuf[m++];
-    yf[ntotal] = mybuf[m++];
-    zf[ntotal] = mybuf[m++];
-    ntotal++;
-  }
-
-  // if last chunk of atoms in this snapshot, write global arrays to file
-
-  if (ntotal == natoms) {
-    write_frame();
-    ntotal = 0;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int DumpDCD::modify_param(int narg, char **arg)
-{
-  if (strcmp(arg[0],"unwrap") == 0) {
-    if (narg < 2) error->all("Illegal dump_modify command");
-    if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
-    else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
-    else error->all("Illegal dump_modify command");
-    return 2;
-  }
-  return 0;
-}
-
-/* ----------------------------------------------------------------------
-   return # of bytes of allocated memory in buf and global coords array
-------------------------------------------------------------------------- */
-
-double DumpDCD::memory_usage()
-{
-  double bytes = Dump::memory_usage();
-  bytes += 3*natoms * sizeof(float);
-  return bytes;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::write_frame()
-{
-  // write coords
-
-  uint32_t out_integer = natoms*sizeof(float);
-  fwrite_int32(fp,out_integer);
-  fwrite(xf,out_integer,1,fp);
-  fwrite_int32(fp,out_integer);
-  fwrite_int32(fp,out_integer);
-  fwrite(yf,out_integer,1,fp);
-  fwrite_int32(fp,out_integer);
-  fwrite_int32(fp,out_integer);
-  fwrite(zf,out_integer,1,fp);
-  fwrite_int32(fp,out_integer);
-  
-  // update NFILE and NSTEP fields in DCD header
-
-  nframes++;
-  out_integer = nframes;
-  fseek(fp,NFILE_POS,SEEK_SET);
-  fwrite_int32(fp,out_integer);
-  out_integer = update->ntimestep;
-  fseek(fp,NSTEP_POS,SEEK_SET);
-  fwrite_int32(fp,out_integer);
-  fseek(fp,0,SEEK_END);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpDCD::write_dcd_header(const char *remarks)
-{
-  uint32_t out_integer;
-  float out_float;
-  char title_string[200];
-  time_t cur_time;
-  struct tm *tmbuf;
-
-  out_integer = 84;
-  fwrite_int32(fp,out_integer);
-  strcpy(title_string,"CORD");
-  fwrite(title_string,4,1,fp);
-  fwrite_int32(fp,0);                    // NFILE = # of snapshots in file
-  fwrite_int32(fp,update->ntimestep);    // START = timestep of first snapshot
-  fwrite_int32(fp,nevery_save);          // SKIP = interval between snapshots
-  fwrite_int32(fp,update->ntimestep);    // NSTEP = timestep of last snapshot
-  fwrite_int32(fp,0);			 // NAMD writes NSTEP or ISTART
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  out_float = update->dt;
-  fwrite(&out_float,sizeof(float),1,fp);
-  fwrite_int32(fp,1);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,0);
-  fwrite_int32(fp,24);                   // pretend to be Charmm version 24
-  fwrite_int32(fp,84);
-  fwrite_int32(fp,164);
-  fwrite_int32(fp,2);
-  strncpy(title_string,remarks,80);
-  title_string[79] = '\0';
-  fwrite(title_string,80,1,fp);
-  cur_time=time(NULL);
-  tmbuf=localtime(&cur_time);
-  memset(title_string,' ',81);
-  strftime(title_string,80,"REMARKS Created %d %B,%Y at %R",tmbuf);
-  fwrite(title_string,80,1,fp);
-  fwrite_int32(fp,164);
-  fwrite_int32(fp,4);
-  fwrite_int32(fp,natoms);                // number of atoms in each snapshot
-  fwrite_int32(fp,4);
-  if (flush_flag) fflush(fp);
-}
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
+                        Axel Kohlmeyer (Temple U), support for groups
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "inttypes.h"
+#include "stdio.h"
+#include "time.h"
+#include "string.h"
+#include "dump_dcd.h"
+#include "domain.h"
+#include "atom.h"
+#include "update.h"
+#include "output.h"
+#include "group.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define NFILE_POS 8L
+#define NSTEP_POS 20L
+
+// necessary to set SEEK params b/c MPI-2 messes with these settings
+
+#ifndef SEEK_SET
+#define SEEK_SET        0
+#define SEEK_CUR        1
+#define SEEK_END        2
+#endif
+
+/* ---------------------------------------------------------------------- */
+
+static inline void fwrite_int32(FILE* fd, uint32_t i)
+{
+  fwrite(&i,sizeof(uint32_t),1,fd);
+}
+
+/* ---------------------------------------------------------------------- */
+
+DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
+{
+  if (narg != 5) error->all("Illegal dump dcd command");
+  if (binary || compressed || multifile || multiproc)
+    error->all("Invalid dump dcd filename");
+
+  size_one = 3;
+  sort_flag = 1;
+  sortcol = 0;
+
+  unwrap_flag = 0;
+  format_default = NULL;
+    
+  // allocate global array for atom coords
+
+  if (igroup == 0) natoms = static_cast<int> (atom->natoms);
+  else natoms = static_cast<int> (group->count(igroup));
+  if (natoms <= 0) error->all("Invalid natoms for dump dcd");
+
+  coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
+  xf = &coords[0*natoms];
+  yf = &coords[1*natoms];
+  zf = &coords[2*natoms];
+
+  openfile();
+  headerflag = 0;
+  nevery_save = 0;
+  ntotal = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+DumpDCD::~DumpDCD()
+{
+  memory->sfree(coords);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::init_style()
+{
+  if (sort_flag == 0 || sortcol != 0)
+    error->all("Dump dcd requires sorting by atom ID");
+
+  // check that dump frequency has not changed and is not a variable
+
+  int idump;
+  for (idump = 0; idump < output->ndump; idump++)
+    if (strcmp(id,output->dump[idump]->id) == 0) break;
+  if (output->every_dump[idump] == 0)
+    error->all("Cannot use variable every setting for dump dcd");
+
+  if (nevery_save == 0) nevery_save = output->every_dump[idump];
+  else if (nevery_save != output->every_dump[idump])
+    error->all("Cannot change dump_modify every for dump dcd");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::openfile()
+{
+  if (me == 0) {
+    fp = fopen(filename,"wb");
+    if (fp == NULL) error->one("Cannot open dump file");
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::write_header(int n)
+{
+  if (n != natoms) error->all("Dump dcd of non-matching # of atoms");
+
+  // first time, write header for entire file
+
+  if (headerflag == 0) {
+    if (me == 0) write_dcd_header("Written by LAMMPS");
+    headerflag = 1;
+    nframes = 0;
+  }
+
+  // dim[] = size and angle cosines of orthogonal or triclinic box
+  // dim[0] = a = length of unit cell vector along x-axis
+  // dim[1] = gamma = cosine of angle between a and b
+  // dim[2] = b = length of unit cell vector in xy-plane
+  // dim[3] = beta = cosine of angle between a and c
+  // dim[4] = alpha = cosine of angle between b and c
+  // dim[5] = c = length of final unit cell vector
+  // 48 = 6 doubles
+
+  double dim[6];
+  if (domain->triclinic) {
+    double *h = domain->h;
+    double alen = h[0];
+    double blen = sqrt(h[5]*h[5] + h[1]*h[1]);
+    double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]);
+    dim[0] = alen;
+    dim[2] = blen;
+    dim[5] = clen;
+    dim[4] = (h[5]*h[4] + h[1]*h[3]) / blen/clen; // alpha
+    dim[3] = (h[0]*h[4]) / alen/clen;             // beta
+    dim[1] = (h[0]*h[5]) / alen/blen;             // gamma
+  } else {
+    dim[0] = domain->xprd;
+    dim[2] = domain->yprd;
+    dim[5] = domain->zprd;
+    dim[1] = dim[3] = dim[4] = 0.0;
+  }
+  
+  if (me == 0) {
+    uint32_t out_integer = 48;
+    fwrite_int32(fp,out_integer);
+    fwrite(dim,out_integer,1,fp); 
+    fwrite_int32(fp,out_integer);
+    if (flush_flag) fflush(fp);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int DumpDCD::count()
+{
+  if (igroup == 0) return atom->nlocal;
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int m = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) m++;
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::pack(int *ids)
+{
+  int m,n;
+
+  int *tag = atom->tag;
+  double **x = atom->x;
+  int *image = atom->image;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  m = n = 0;
+  if (unwrap_flag) {
+    double xprd = domain->xprd;
+    double yprd = domain->yprd;
+    double zprd = domain->zprd;
+    double xy = domain->xy;
+    double xz = domain->xz;
+    double yz = domain->yz;
+
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	int ix = (image[i] & 1023) - 512;
+	int iy = (image[i] >> 10 & 1023) - 512;
+	int iz = (image[i] >> 20) - 512;
+
+	if (domain->triclinic) {
+	  buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz;
+	  buf[m++] = x[i][1] + iy * yprd + iz * yz;
+	  buf[m++] = x[i][2] + iz * zprd;
+	} else {
+	  buf[m++] = x[i][0] + ix * xprd;
+	  buf[m++] = x[i][1] + iy * yprd;
+	  buf[m++] = x[i][2] + iz * zprd;
+	}
+	ids[n++] = tag[i];
+      }
+    }
+
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	buf[m++] = x[i][0];
+	buf[m++] = x[i][1];
+	buf[m++] = x[i][2];
+	ids[n++] = tag[i];
+      }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::write_data(int n, double *mybuf)
+{
+  // copy buf atom coords into 3 global arrays
+
+  int m = 0;
+  for (int i = 0; i < n; i++) {
+    xf[ntotal] = mybuf[m++];
+    yf[ntotal] = mybuf[m++];
+    zf[ntotal] = mybuf[m++];
+    ntotal++;
+  }
+
+  // if last chunk of atoms in this snapshot, write global arrays to file
+
+  if (ntotal == natoms) {
+    write_frame();
+    ntotal = 0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int DumpDCD::modify_param(int narg, char **arg)
+{
+  if (strcmp(arg[0],"unwrap") == 0) {
+    if (narg < 2) error->all("Illegal dump_modify command");
+    if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
+    else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
+    else error->all("Illegal dump_modify command");
+    return 2;
+  }
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   return # of bytes of allocated memory in buf and global coords array
+------------------------------------------------------------------------- */
+
+double DumpDCD::memory_usage()
+{
+  double bytes = Dump::memory_usage();
+  bytes += 3*natoms * sizeof(float);
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::write_frame()
+{
+  // write coords
+
+  uint32_t out_integer = natoms*sizeof(float);
+  fwrite_int32(fp,out_integer);
+  fwrite(xf,out_integer,1,fp);
+  fwrite_int32(fp,out_integer);
+  fwrite_int32(fp,out_integer);
+  fwrite(yf,out_integer,1,fp);
+  fwrite_int32(fp,out_integer);
+  fwrite_int32(fp,out_integer);
+  fwrite(zf,out_integer,1,fp);
+  fwrite_int32(fp,out_integer);
+  
+  // update NFILE and NSTEP fields in DCD header
+
+  nframes++;
+  out_integer = nframes;
+  fseek(fp,NFILE_POS,SEEK_SET);
+  fwrite_int32(fp,out_integer);
+  out_integer = update->ntimestep;
+  fseek(fp,NSTEP_POS,SEEK_SET);
+  fwrite_int32(fp,out_integer);
+  fseek(fp,0,SEEK_END);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpDCD::write_dcd_header(const char *remarks)
+{
+  uint32_t out_integer;
+  float out_float;
+  char title_string[200];
+  time_t cur_time;
+  struct tm *tmbuf;
+
+  out_integer = 84;
+  fwrite_int32(fp,out_integer);
+  strcpy(title_string,"CORD");
+  fwrite(title_string,4,1,fp);
+  fwrite_int32(fp,0);                    // NFILE = # of snapshots in file
+  fwrite_int32(fp,update->ntimestep);    // START = timestep of first snapshot
+  fwrite_int32(fp,nevery_save);          // SKIP = interval between snapshots
+  fwrite_int32(fp,update->ntimestep);    // NSTEP = timestep of last snapshot
+  fwrite_int32(fp,0);			 // NAMD writes NSTEP or ISTART
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  out_float = update->dt;
+  fwrite(&out_float,sizeof(float),1,fp);
+  fwrite_int32(fp,1);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,0);
+  fwrite_int32(fp,24);                   // pretend to be Charmm version 24
+  fwrite_int32(fp,84);
+  fwrite_int32(fp,164);
+  fwrite_int32(fp,2);
+  strncpy(title_string,remarks,80);
+  title_string[79] = '\0';
+  fwrite(title_string,80,1,fp);
+  cur_time=time(NULL);
+  tmbuf=localtime(&cur_time);
+  memset(title_string,' ',81);
+  strftime(title_string,80,"REMARKS Created %d %B,%Y at %R",tmbuf);
+  fwrite(title_string,80,1,fp);
+  fwrite_int32(fp,164);
+  fwrite_int32(fp,4);
+  fwrite_int32(fp,natoms);                // number of atoms in each snapshot
+  fwrite_int32(fp,4);
+  if (flush_flag) fflush(fp);
+}
diff --git a/src/finish.cpp b/src/finish.cpp
index ac27ad2bdf..4512a3f97c 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -44,9 +44,10 @@ void Finish::end(int flag)
   int i,m,nneigh,nneighfull;
   int histo[10];
   int loopflag,minflag,prdflag,timeflag,fftflag,histoflag,neighflag;
-  double time,tmp,ave,max,min,natoms;
+  double time,tmp,ave,max,min;
   double time_loop,time_other;
-  
+  bigint natoms;
+
   int me,nprocs;
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
@@ -85,17 +86,17 @@ void Finish::end(int flag)
     // overall loop time
     // use actual natoms, in case atoms were lost
 
-    double rlocal = atom->nlocal;
-    MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+    bigint nblocal = atom->nlocal;
+    MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
     
     if (me == 0) {
       if (screen) 
 	fprintf(screen,
-		"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
+		"Loop time of %g on %d procs for %d steps with %lu atoms\n",
 		time_loop,nprocs,update->nsteps,natoms);
       if (logfile)
 	fprintf(logfile,
-		"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
+		"Loop time of %g on %d procs for %d steps with %lu atoms\n",
 		time_loop,nprocs,update->nsteps,natoms);
     }
 
diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp
index 696ad1e520..5a7ea44709 100644
--- a/src/fix_orient_fcc.cpp
+++ b/src/fix_orient_fcc.cpp
@@ -404,13 +404,13 @@ void FixOrientFCC::post_force(int vflag)
 
     if (me == 0) { 
       if (screen) {
-	fprintf(screen,"orient step %d: %g atoms have %d neighbors\n",
+	fprintf(screen,"orient step %d: %lu atoms have %d neighbors\n",
 		update->ntimestep,atom->natoms,total);
 	fprintf(screen,"  neighs: min = %d, max = %d, ave = %g\n",
 		min,max,ave);
       }
       if (logfile) {
-	fprintf(logfile,"orient step %d: %g atoms have %d neighbors\n",
+	fprintf(logfile,"orient step %d: %lu atoms have %d neighbors\n",
 		update->ntimestep,atom->natoms,total);
 	fprintf(logfile,"  neighs: min = %d, max = %d, ave = %g\n",
 		min,max,ave);
diff --git a/src/library.cpp b/src/library.cpp
index 55d63a74d1..546ad4f3f7 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -18,6 +18,7 @@
 #include "string.h"
 #include "stdlib.h"
 #include "library.h"
+#include "lmptype.h"
 #include "lammps.h"
 #include "input.h"
 #include "atom.h"
@@ -347,6 +348,7 @@ void *lammps_extract_variable(void *ptr, char *name, char *group)
 int lammps_get_natoms(void *ptr)
 {
   LAMMPS *lmp = (LAMMPS *) ptr;
+  if (lmp->atom->natoms > MAXINT32) return 0;
   int natoms = static_cast<int> (lmp->atom->natoms);
   return natoms;
 }
@@ -360,6 +362,7 @@ void lammps_get_coords(void *ptr, double *coords)
   // error if tags are not defined or not consecutive
 
   if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) return;
+  if (lmp->atom->natoms > MAXINT32) return;
 
   int natoms = static_cast<int> (lmp->atom->natoms);
   double *copy = new double[3*natoms];
@@ -391,6 +394,7 @@ void lammps_put_coords(void *ptr, double *coords)
   // error if no map defined by LAMMPS
 
   if (lmp->atom->map_style == 0) return;
+  if (lmp->atom->natoms > MAXINT32) return;
 
   int natoms = static_cast<int> (lmp->atom->natoms);
   double **x = lmp->atom->x;
diff --git a/src/memory.cpp b/src/memory.cpp
index 4e46b87b57..c4eb9eba2a 100644
--- a/src/memory.cpp
+++ b/src/memory.cpp
@@ -16,13 +16,25 @@
 #include "stdlib.h"
 #include "string.h"
 #include "memory.h"
+#include "lmptype.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-Memory::Memory(LAMMPS *lmp) : Pointers(lmp) {}
+Memory::Memory(LAMMPS *lmp) : Pointers(lmp)
+{
+  // check datatype sizes
+
+  if (sizeof(bigint) != 8)
+    error->all("No support for 8-byte unsigned integers");
+
+  int mpisize;
+  MPI_Type_size(MPI_UNSIGNED_LONG,&mpisize);
+  if (mpisize != 8)
+    error->all("MPI_UNSIGNED_LONG is not 8-byte data type");
+}
 
 /* ----------------------------------------------------------------------
    safe malloc 
diff --git a/src/min_cg.cpp b/src/min_cg.cpp
index 17c5bf446f..63ed35fddd 100644
--- a/src/min_cg.cpp
+++ b/src/min_cg.cpp
@@ -11,10 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "mpi.h"
 #include "math.h"
 #include "string.h"
-#include "mpi.h"
 #include "min_cg.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "update.h"
 #include "output.h"
@@ -22,8 +23,6 @@
 
 using namespace LAMMPS_NS;
 
-#define MAXATOMS 0x7FFFFFFF
-
 // EPS_ENERGY = minimum normalization for energy tolerance
 
 #define EPS_ENERGY 1.0e-8
@@ -52,7 +51,7 @@ int MinCG::iterate(int maxiter)
   // nlimit = max # of CG iterations before restarting
   // set to ndoftotal unless too big
 
-  int nlimit = static_cast<int> (MIN(MAXATOMS,ndoftotal));
+  int nlimit = static_cast<int> (MIN(MAXINT32,ndoftotal));
 
   // initialize working vectors
 
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 1e769eda3c..fb4d3efe90 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -16,6 +16,7 @@
 #include "string.h"
 #include "stdlib.h"
 #include "read_data.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "comm.h"
@@ -335,7 +336,7 @@ void ReadData::header(int flag)
 
     // search line for header keyword and set corresponding variable
 
-    if (strstr(line,"atoms")) sscanf(line,"%lg",&atom->natoms);
+    if (strstr(line,"atoms")) sscanf(line,"%lu",&atom->natoms);
     else if (strstr(line,"bonds")) sscanf(line,"%d",&atom->nbonds);
     else if (strstr(line,"angles")) sscanf(line,"%d",&atom->nangles);
     else if (strstr(line,"dihedrals")) sscanf(line,"%d",&atom->ndihedrals);
@@ -402,19 +403,18 @@ void ReadData::header(int flag)
 
 /* ----------------------------------------------------------------------
    read all atoms
-   accumulate nread in double precision to allow natoms > 2^31
 ------------------------------------------------------------------------- */
 
 void ReadData::atoms()
 {
   int i,m,nchunk;
   
-  double nread = 0.0;
-  double natoms = atom->natoms;
+  bigint nread = 0;
+  bigint natoms = atom->natoms;
 
   while (nread < natoms) {
     if (natoms-nread > CHUNK) nchunk = CHUNK;
-    else nchunk = static_cast<int> (natoms - nread);
+    else nchunk = natoms-nread;
     if (me == 0) {
       char *eof;
       m = 0;
@@ -449,24 +449,27 @@ void ReadData::atoms()
   // if any atom ID > 0, error if any atom ID == 0
   // not checking if atom IDs > natoms or are unique
   
+  int nlocal = atom->nlocal;
+  int *tag = atom->tag;
+
   int flag = 0;
-  for (int i = 0; i < atom->nlocal; i++)
-    if (atom->tag[i] < 0) flag = 1;
+  for (int i = 0; i < nlocal; i++)
+    if (tag[i] < 0) flag = 1;
   int flag_all;
   MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
   if (flag_all)
     error->all("Invalid atom ID in Atoms section of data file");
 
   flag = 0;
-  for (int i = 0; i < atom->nlocal; i++)
-    if (atom->tag[i] > 0) flag = 1;
+  for (int i = 0; i < nlocal; i++)
+    if (tag[i] > 0) flag = 1;
   MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
   if (flag_all == 0) atom->tag_enable = 0;
 
   if (atom->tag_enable) {
     flag = 0;
-    for (int i = 0; i < atom->nlocal; i++)
-      if (atom->tag[i] == 0) flag = 1;
+    for (int i = 0; i < nlocal; i++)
+      if (tag[i] == 0) flag = 1;
     MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
     if (flag_all)
       error->all("Invalid atom ID in Atoms section of data file");
@@ -483,7 +486,6 @@ void ReadData::atoms()
 /* ----------------------------------------------------------------------
    read all velocities
    to find atoms, must build atom map if not a molecular system 
-   accumulate nread in double precision to allow natoms > 2^31
 ------------------------------------------------------------------------- */
 
 void ReadData::velocities()
@@ -498,12 +500,12 @@ void ReadData::velocities()
     atom->map_set();
   }
 
-  double nread = 0.0;
-  double natoms = atom->natoms;
+  bigint nread = 0;
+  bigint natoms = atom->natoms;
 
   while (nread < natoms) {
     if (natoms-nread > CHUNK) nchunk = CHUNK;
-    else nchunk = static_cast<int> (natoms - nread);
+    else nchunk = natoms-nread;
     if (me == 0) {
       char *eof;
       m = 0;
@@ -560,9 +562,11 @@ void ReadData::bonds()
 
   // check that bonds were assigned correctly
 
+  int nlocal = atom->nlocal;
+
   int sum;
   int n = 0;
-  for (i = 0; i < atom->nlocal; i++) n += atom->num_bond[i];
+  for (i = 0; i < nlocal; i++) n += atom->num_bond[i];
   MPI_Allreduce(&n,&sum,1,MPI_INT,MPI_SUM,world);
   int factor = 1;
   if (!force->newton_bond) factor = 2;
@@ -602,9 +606,11 @@ void ReadData::angles()
 
   // check that angles were assigned correctly
 
+  int nlocal = atom->nlocal;
+
   int sum;
   int n = 0;
-  for (i = 0; i < atom->nlocal; i++) n += atom->num_angle[i];
+  for (i = 0; i < nlocal; i++) n += atom->num_angle[i];
   MPI_Allreduce(&n,&sum,1,MPI_INT,MPI_SUM,world);
   int factor = 1;
   if (!force->newton_bond) factor = 3;
@@ -644,9 +650,11 @@ void ReadData::dihedrals()
 
   // check that dihedrals were assigned correctly
 
+  int nlocal = atom->nlocal;
+
   int sum;
   int n = 0;
-  for (i = 0; i < atom->nlocal; i++) n += atom->num_dihedral[i];
+  for (i = 0; i < nlocal; i++) n += atom->num_dihedral[i];
   MPI_Allreduce(&n,&sum,1,MPI_INT,MPI_SUM,world);
   int factor = 1;
   if (!force->newton_bond) factor = 4;
@@ -687,9 +695,11 @@ void ReadData::impropers()
 
   // check that impropers were assigned correctly
 
+  int nlocal = atom->nlocal;
+
   int sum;
   int n = 0;
-  for (i = 0; i < atom->nlocal; i++) n += atom->num_improper[i];
+  for (i = 0; i < nlocal; i++) n += atom->num_improper[i];
   MPI_Allreduce(&n,&sum,1,MPI_INT,MPI_SUM,world);
   int factor = 1;
   if (!force->newton_bond) factor = 4;
diff --git a/src/read_restart.cpp b/src/read_restart.cpp
index 2b8f59e3bd..1b759d8840 100644
--- a/src/read_restart.cpp
+++ b/src/read_restart.cpp
@@ -17,6 +17,7 @@
 #include "sys/types.h"
 #include "dirent.h"
 #include "read_restart.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "domain.h"
@@ -290,13 +291,13 @@ void ReadRestart::command(int narg, char **arg)
 
   // check that all atoms were assigned to procs
 
-  double natoms;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint natoms;
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
 
   if (me == 0) {
-    if (screen) fprintf(screen,"  %.15g atoms\n",natoms);
-    if (logfile) fprintf(logfile,"  %.15g atoms\n",natoms);
+    if (screen) fprintf(screen,"  %lu atoms\n",natoms);
+    if (logfile) fprintf(logfile,"  %lu atoms\n",natoms);
   }
 
   if (natoms != atom->natoms) error->all("Did not assign all atoms correctly");
@@ -594,7 +595,7 @@ void ReadRestart::header()
       delete [] style;
 
     } else if (flag == NATOMS) {
-      atom->natoms = read_double();
+      atom->natoms = read_bigint();
     } else if (flag == NTYPES) {
       atom->ntypes = read_int();
     } else if (flag == NBONDS) {
@@ -809,3 +810,15 @@ char *ReadRestart::read_char()
   MPI_Bcast(value,n,MPI_CHAR,0,world);
   return value;
 }
+
+/* ----------------------------------------------------------------------
+   read a bigint from restart file and bcast it
+------------------------------------------------------------------------- */
+
+bigint ReadRestart::read_bigint()
+{
+  bigint value;
+  if (me == 0) fread(&value,sizeof(bigint),1,fp);
+  MPI_Bcast(&value,1,MPI_UNSIGNED_LONG,0,world);
+  return value;
+}
diff --git a/src/read_restart.h b/src/read_restart.h
index fd35f8c97b..d5364a1447 100644
--- a/src/read_restart.h
+++ b/src/read_restart.h
@@ -22,6 +22,7 @@ CommandStyle(read_restart,ReadRestart)
 
 #include "stdio.h"
 #include "pointers.h"
+#include "lmptype.h"
 
 namespace LAMMPS_NS {
 
@@ -43,6 +44,7 @@ class ReadRestart : protected Pointers {
   int read_int();
   double read_double();
   char *read_char();
+  bigint read_bigint();
 };
 
 }
diff --git a/src/replicate.cpp b/src/replicate.cpp
index 5b6951ab56..6e8d1bffe4 100644
--- a/src/replicate.cpp
+++ b/src/replicate.cpp
@@ -13,7 +13,9 @@
 
 #include "stdlib.h"
 #include "string.h"
+#include "stdint.h"
 #include "replicate.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "atom_vec_hybrid.h"
@@ -27,7 +29,6 @@
 using namespace LAMMPS_NS;
 
 #define LB_FACTOR 1.1
-#define MAXATOMS  0x7FFFFFFF
 #define EPSILON   1.0e-6
 
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
@@ -129,16 +130,16 @@ void Replicate::command(int narg, char **arg)
   // if molecular, N/Nbonds/etc cannot be > 2^31 else tags/counts invalid
 
   double rep = nrep;
-  if (rep*old->natoms > MAXATOMS) atom->tag_enable = 0;
+  if (rep*old->natoms > MAXINT32) atom->tag_enable = 0;
 
   if (atom->tag_enable == 0)
     for (int i = 0; i < atom->nlocal; i++)
       atom->tag[i] = 0;
 
   if (atom->molecular) {
-    if (rep*old->natoms > MAXATOMS || rep*old->nbonds > MAXATOMS ||
-	rep*old->nangles > MAXATOMS || rep*old->ndihedrals > MAXATOMS ||
-	rep*old->nimpropers > MAXATOMS)
+    if (rep*old->natoms > MAXINT32 || rep*old->nbonds > MAXINT32 ||
+	rep*old->nangles > MAXINT32 || rep*old->ndihedrals > MAXINT32 ||
+	rep*old->nimpropers > MAXINT32)
       error->all("Too big a problem to replicate with molecular atom style");
   }
 
@@ -363,13 +364,13 @@ void Replicate::command(int narg, char **arg)
 
   // check that all atoms were assigned to procs
 
-  double natoms;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint natoms;
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
 
   if (me == 0) {
-    if (screen) fprintf(screen,"  %.15g atoms\n",natoms);
-    if (logfile) fprintf(logfile,"  %.15g atoms\n",natoms);
+    if (screen) fprintf(screen,"  %lu atoms\n",natoms);
+    if (logfile) fprintf(logfile,"  %lu atoms\n",natoms);
   }
 
   if (natoms != atom->natoms)
diff --git a/src/thermo.cpp b/src/thermo.cpp
index e31fd8a8e2..60997baf29 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -56,7 +56,7 @@ using namespace LAMMPS_NS;
 
 enum{IGNORE,WARN,ERROR};           // same as write_restart.cpp
 enum{ONELINE,MULTILINE};
-enum{INT,FLOAT};
+enum{INT,FLOAT,BIGINT};
 enum{SCALAR,VECTOR,ARRAY};
 
 #define INVOKED_SCALAR 1
@@ -135,12 +135,15 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
   format_multi = (char *) "---------------- Step %8d ----- "
                           "CPU = %11.4f (sec) ----------------";
+  format_float_one_def = (char *) "%12.8g";
+  format_float_multi_def = (char *) "%14.4f";
   format_int_one_def = (char *) "%8d";
   format_int_multi_def = (char *) "%14d";
-  format_g_def = (char *) "%12.8g";
-  format_f_def = (char *) "%14.4f";
-  format_int_user = NULL;
+  format_bigint_one_def = (char *) "%8lu";
+  format_bigint_multi_def = (char *) "%14lu";
   format_float_user = NULL;
+  format_int_user = NULL;
+  format_bigint_user = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -154,8 +157,9 @@ Thermo::~Thermo()
 
   // format strings
 
-  delete [] format_int_user;
   delete [] format_float_user;
+  delete [] format_int_user;
+  delete [] format_bigint_user;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -188,13 +192,19 @@ void Thermo::init()
     if (lineflag == MULTILINE && i % 3 == 0) strcat(format[i],"\n");
 
     if (format_user[i]) ptr = format_user[i];
-    else if (vtype[i] == INT && format_int_user) ptr = format_int_user;
-    else if (vtype[i] == INT && lineflag == ONELINE) ptr = format_int_one_def;
-    else if (vtype[i] == INT && lineflag == MULTILINE) 
-      ptr = format_int_multi_def;
-    else if (vtype[i] == FLOAT && format_float_user) ptr = format_float_user;
-    else if (lineflag == ONELINE) ptr = format_g_def;
-    else if (lineflag == MULTILINE) ptr = format_f_def;
+    else if (vtype[i] == FLOAT) {
+      if (format_float_user) ptr = format_float_user;
+      else if (lineflag == ONELINE) ptr = format_float_one_def;
+      else if (lineflag == MULTILINE) ptr = format_float_multi_def;
+    } else if (vtype[i] == INT) {
+      if (format_int_user) ptr = format_int_user;
+      else if (lineflag == ONELINE) ptr = format_int_one_def;
+      else if (lineflag == MULTILINE) ptr = format_int_multi_def;
+    } else if (vtype[i] == BIGINT) {
+      if (format_bigint_user) ptr = format_bigint_user;
+      else if (lineflag == ONELINE) ptr = format_bigint_one_def;
+      else if (lineflag == MULTILINE) ptr = format_bigint_multi_def;
+    }
 
     n = strlen(format[i]);
     if (lineflag == ONELINE) sprintf(&format[i][n],"%s ",ptr);
@@ -310,8 +320,13 @@ void Thermo::compute(int flag)
 
   for (ifield = 0; ifield < nfield; ifield++) {
     (this->*vfunc[ifield])();
-    if (vtype[ifield] == INT) loc += sprintf(&line[loc],format[ifield],ivalue);
-    else loc += sprintf(&line[loc],format[ifield],dvalue);
+    if (vtype[ifield] == FLOAT)
+      loc += sprintf(&line[loc],format[ifield],dvalue);
+    else if (vtype[ifield] == INT) 
+      loc += sprintf(&line[loc],format[ifield],ivalue);
+    else if (vtype[ifield] == BIGINT) {
+      loc += sprintf(&line[loc],format[ifield],bivalue);
+    }
   }
 
   // kludge for RedStorm timing issue
@@ -332,13 +347,13 @@ void Thermo::compute(int flag)
    check for lost atoms, return current number of atoms
 ------------------------------------------------------------------------- */
 
-double Thermo::lost_check()
+bigint Thermo::lost_check()
 {
   // ntotal = current # of atoms
 
-  double ntotal;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&ntotal,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint ntotal;
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&ntotal,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
   if (ntotal == atom->natoms) return ntotal;
 
   // if not checking or already warned, just return
@@ -350,15 +365,14 @@ double Thermo::lost_check()
 
   if (lostflag == ERROR) {
     char str[128];
-    sprintf(str,"Lost atoms: original %.15g current %.15g",
-	    atom->natoms,ntotal);
+    sprintf(str,"Lost atoms: original %lu current %lu",atom->natoms,ntotal);
     error->all(str);
   }
 
   // warning message
 
   char str[128];
-  sprintf(str,"Lost atoms: original %.15g current %.15g",atom->natoms,ntotal);
+  sprintf(str,"Lost atoms: original %lu current %lu",atom->natoms,ntotal);
   if (me == 0) error->warning(str,0);
   lostbefore = 1;
   return ntotal;
@@ -484,6 +498,15 @@ void Thermo::modify_params(int narg, char **arg)
 	int n = strlen(arg[iarg+2]) + 1;
 	format_int_user = new char[n];
 	strcpy(format_int_user,arg[iarg+2]);
+	if (format_bigint_user) delete [] format_bigint_user;
+	n = strlen(format_int_user) + 2;
+	format_bigint_user = new char[n];
+	char *ptr = strchr(format_int_user,'d');
+	if (ptr == NULL) 
+	  error->all("Thermo_modify int format does not contain d character");
+	*ptr = '\0';
+	sprintf(format_bigint_user,"%s%s%s",format_int_user,"lu",ptr+1);
+	*ptr = 'd';
       } else if (strcmp(arg[iarg+1],"float") == 0) {
 	if (format_float_user) delete [] format_float_user;
 	int n = strlen(arg[iarg+2]) + 1;
@@ -610,7 +633,7 @@ void Thermo::parse_fields(char *str)
       addfield("S/CPU",&Thermo::compute_spcpu,FLOAT);
 
     } else if (strcmp(word,"atoms") == 0) {
-      addfield("Atoms",&Thermo::compute_atoms,INT);
+      addfield("Atoms",&Thermo::compute_atoms,BIGINT);
     } else if (strcmp(word,"temp") == 0) {
       addfield("Temp",&Thermo::compute_temp,FLOAT);
       index_temp = add_compute(id_temp,SCALAR);
@@ -942,7 +965,6 @@ int Thermo::evaluate_keyword(char *word, double *answer)
 
   } else if (strcmp(word,"atoms") == 0) {
     compute_atoms();
-    dvalue = ivalue;
 
   } else if (strcmp(word,"temp") == 0) {
     if (!temperature)
@@ -1409,7 +1431,7 @@ void Thermo::compute_spcpu()
 
 void Thermo::compute_atoms()
 {
-  ivalue = static_cast<int> (natoms);
+  bivalue = natoms;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/thermo.h b/src/thermo.h
index c35b91e168..526fcfa89f 100644
--- a/src/thermo.h
+++ b/src/thermo.h
@@ -15,6 +15,7 @@
 #define LMP_THERMO_H
 
 #include "pointers.h"
+#include "lmptype.h"
 
 namespace LAMMPS_NS {
 
@@ -30,7 +31,7 @@ class Thermo : protected Pointers {
   Thermo(class LAMMPS *, int, char **);
   ~Thermo();
   void init();
-  double lost_check();
+  bigint lost_check();
   void modify_params(int, char **);
   void header();
   void compute(int);
@@ -44,10 +45,12 @@ class Thermo : protected Pointers {
   char **keyword;
   int *vtype;
 
-  char *format_multi,*format_int_one_def,*format_int_multi_def;
-  char *format_g_def,*format_f_def;
-  char *format_int_user,*format_float_user;
   char **format,**format_user;
+  char *format_multi;
+  char *format_float_one_def,*format_float_multi_def;
+  char *format_int_one_def,*format_int_multi_def;
+  char *format_bigint_one_def,*format_bigint_multi_def;
+  char *format_float_user,*format_int_user,*format_bigint_user;
 
   int normvalue;         // use this for normflag unless natoms = 0
   int normuserflag;      // 0 if user has not set, 1 if has
@@ -61,9 +64,12 @@ class Thermo : protected Pointers {
   double last_time;
   int last_step;
 
+  bigint natoms;
+
                          // data used by routines that compute single values
   int ivalue;            // integer value to print
-  double dvalue,natoms;  // dvalue = double value to print
+  double dvalue;         // dvalue = double value to print
+  bigint bivalue;        // big integer value to print
   int ifield;            // which field in thermo output is being computed
   int *field2index;      // which compute,fix,variable calcs this field
   int *argindex1;        // indices into compute,fix scalar,vector
diff --git a/src/velocity.cpp b/src/velocity.cpp
index 506e21f288..a3f84537b6 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -17,6 +17,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "velocity.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "update.h"
 #include "domain.h"
@@ -222,6 +223,8 @@ void Velocity::create(double t_desired, int seed)
 
     // error check
 
+    if (atom->natoms > MAXINT32)
+      error->all("Too big a problem to use velocity create loop all");
     if (atom->tag_enable == 0)
       error->all("Cannot use velocity create loop all unless atoms have IDs");
     if (atom->tag_consecutive() == 0)
diff --git a/src/write_restart.cpp b/src/write_restart.cpp
index fc9748fde2..c67e88271e 100644
--- a/src/write_restart.cpp
+++ b/src/write_restart.cpp
@@ -14,6 +14,7 @@
 #include "mpi.h"
 #include "string.h"
 #include "write_restart.h"
+#include "lmptype.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "atom_vec_hybrid.h"
@@ -120,8 +121,8 @@ void WriteRestart::write(char *file)
   // natoms = sum of nlocal = value to write into restart file
   // if unequal and thermo lostflag is "error", don't write restart file
 
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+  bigint nblocal = atom->nlocal;
+  MPI_Allreduce(&nblocal,&natoms,1,MPI_UNSIGNED_LONG,MPI_SUM,world);
   if (natoms != atom->natoms && output->thermo->lostflag == ERROR) 
     error->all("Atom count is inconsistent, cannot write restart file");
 
@@ -331,7 +332,7 @@ void WriteRestart::header()
     }
   }
 
-  write_double(NATOMS,natoms);
+  write_bigint(NATOMS,natoms);
   write_int(NTYPES,atom->ntypes);
   write_int(NBONDS,atom->nbonds);
   write_int(NBONDTYPES,atom->nbondtypes);
@@ -484,3 +485,14 @@ void WriteRestart::write_char(int flag, char *value)
   fwrite(&n,sizeof(int),1,fp);
   fwrite(value,sizeof(char),n,fp);
 }
+
+/* ----------------------------------------------------------------------
+   write a flag and a bigint into restart file 
+------------------------------------------------------------------------- */
+
+void WriteRestart::write_bigint(int flag, bigint value)
+{
+  fwrite(&flag,sizeof(int),1,fp);
+  fwrite(&value,sizeof(bigint),1,fp);
+}
+
diff --git a/src/write_restart.h b/src/write_restart.h
index c73af2becc..db729ec010 100644
--- a/src/write_restart.h
+++ b/src/write_restart.h
@@ -22,6 +22,7 @@ CommandStyle(write_restart,WriteRestart)
 
 #include "stdio.h"
 #include "pointers.h"
+#include "lmptype.h"
 
 namespace LAMMPS_NS {
 
@@ -34,7 +35,7 @@ class WriteRestart : protected Pointers {
  private:
   int me,nprocs;
   FILE *fp;
-  double natoms;         // natoms (sum of nlocal) to write into file
+  bigint natoms;         // natoms (sum of nlocal) to write into file
 
   void header();
   void type_arrays();
@@ -43,6 +44,7 @@ class WriteRestart : protected Pointers {
   void write_int(int, int);
   void write_double(int, double);
   void write_char(int, char *);
+  void write_bigint(int, bigint);
 };
 
 }