diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index 47821b1056..4f42028f41 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -205,6 +205,9 @@ void PairReaxC::settings(int narg, char **arg)
 
   qeqflag = 1;
   control->lgflag = 0;
+  system->mincap = MIN_CAP;
+  system->safezone = SAFE_ZONE;
+  system->saferzone = SAFER_ZONE;
 
   // process optional keywords
 
@@ -223,6 +226,19 @@ void PairReaxC::settings(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"no") == 0) control->lgflag = 0;
       else error->all(FLERR,"Illegal pair_style reax/c command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"safezone") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style reax/c command");
+      system->safezone = atof(arg[iarg+1]);
+      if (system->safezone < 0.0) 
+	error->all(FLERR,"Illegal pair_style reax/c safezone command");
+      system->saferzone = system->safezone + 0.2;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"mincap") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style reax/c command");
+      system->mincap = atoi(arg[iarg+1]);
+      if (system->mincap < 0) 
+	error->all(FLERR,"Illegal pair_style reax/c mincap command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal pair_style reax/c command");
   }
 
@@ -339,6 +355,8 @@ void PairReaxC::init_style( )
 void PairReaxC::setup( )
 {
   int oldN;
+  int mincap = system->mincap;
+  double safezone = system->safezone;
 
   system->n = atom->nlocal; // my atoms
   system->N = atom->nlocal + atom->nghost; // mine + ghosts
@@ -356,8 +374,8 @@ void PairReaxC::setup( )
 
     // determine the local and total capacity
 
-    system->local_cap = MAX( (int)(system->n * SAFE_ZONE), MIN_CAP );
-    system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
+    system->local_cap = MAX( (int)(system->n * safezone), mincap );
+    system->total_cap = MAX( (int)(system->N * safezone), mincap );
 
     // initialize my data structures
 
@@ -532,6 +550,9 @@ void PairReaxC::write_reax_atoms()
   int *num_bonds = fix_reax->num_bonds;
   int *num_hbonds = fix_reax->num_hbonds;
 
+  if (system->N > system->total_cap)
+    error->all(FLERR,"Too many ghost atoms");
+
   for( int i = 0; i < system->N; ++i ){
     system->my_atoms[i].orig_id = atom->tag[i];
     system->my_atoms[i].type = map[atom->type[i]];
@@ -578,6 +599,9 @@ int PairReaxC::estimate_reax_lists()
   reax_list *far_nbrs;
   far_neighbor_data *far_list;
 
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+
   x = atom->x;
   nlocal = atom->nlocal;
   nghost = atom->nghost;
@@ -616,7 +640,7 @@ int PairReaxC::estimate_reax_lists()
   free( marked );
   free( dist );
 
-  return static_cast<int> (MAX( num_nbrs*SAFE_ZONE, MIN_CAP*MIN_NBRS ));
+  return static_cast<int> (MAX( num_nbrs*safezone, mincap*MIN_NBRS ));
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-REAXC/pair_reax_c.h b/src/USER-REAXC/pair_reax_c.h
index 111074ac88..30b1d4e14a 100644
--- a/src/USER-REAXC/pair_reax_c.h
+++ b/src/USER-REAXC/pair_reax_c.h
@@ -79,3 +79,13 @@ class PairReaxC : public Pair {
 
 #endif
 #endif
+
+/* ERROR/WARNING messages:
+
+E: Too many ghost atoms
+
+Number of ghost atoms has increased too much during simulation and has exceeded 
+the size of reax/c arrays.  Increase safe_zone and min_cap in pair_style reax/c 
+command
+
+*/
diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp
index 039b2a00cb..1fdf444c51 100644
--- a/src/USER-REAXC/reaxc_allocate.cpp
+++ b/src/USER-REAXC/reaxc_allocate.cpp
@@ -49,9 +49,14 @@ int PreAllocate_Space( reax_system *system, control_params *control,
 {
   int  i;
 
-  /* determine the local and total capacity */
-  system->local_cap = MAX( (int)(system->n * SAFE_ZONE), MIN_CAP );
-  system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+
+  // determine the local and total capacity
+
+  system->local_cap = MAX( (int)(system->n * safezone), mincap );
+  system->total_cap = MAX( (int)(system->N * safezone), mincap );
+
 #if defined(DEBUG)
   fprintf( stderr, "p%d: local_cap=%d total_cap=%d\n",
            system->my_rank, system->local_cap, system->total_cap );
@@ -490,6 +495,9 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
 {
   int i, id, total_hbonds;
 
+  int mincap = system->mincap;
+  double saferzone = system->saferzone;
+
   total_hbonds = 0;
   for( i = 0; i < system->n; ++i )
     if( (id = system->my_atoms[i].Hindex) >= 0 ) {
@@ -498,7 +506,7 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
       //                                   MIN_HBONDS);
       total_hbonds += system->my_atoms[i].num_hbonds;
     }
-  total_hbonds = (int)(MAX( total_hbonds*SAFER_ZONE, MIN_CAP*MIN_HBONDS ));
+  total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
 
   Delete_List( hbonds, comm );
   if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds, comm ) ) {
@@ -515,6 +523,9 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
 {
   int i;
 
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+
   *total_bonds = 0;
   *est_3body = 0;
   for( i = 0; i < system->N; ++i ){
@@ -523,7 +534,7 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
     // system->my_atoms[i].num_bonds = MAX( Num_Entries(i,bonds)*2, MIN_BONDS );
     *total_bonds += system->my_atoms[i].num_bonds;
   }
-  *total_bonds = (int)(MAX( *total_bonds * SAFE_ZONE, MIN_CAP*MIN_BONDS ));
+  *total_bonds = (int)(MAX( *total_bonds * safezone, mincap*MIN_BONDS ));
 
   Delete_List( bonds, comm );
   if(!Make_List(system->total_cap, *total_bonds, TYP_BOND, bonds, comm)) {
@@ -760,6 +771,10 @@ void ReAllocate( reax_system *system, control_params *control,
   MPI_Comm comm;
   char msg[200];
 
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+  double saferzone = system->saferzone;
+
   realloc = &(workspace->realloc);
   g = &(system->my_grid);
   comm = mpi_data->world;
@@ -787,14 +802,14 @@ void ReAllocate( reax_system *system, control_params *control,
   if( system->n >= DANGER_ZONE * system->local_cap ||
       (0 && system->n <= LOOSE_ZONE * system->local_cap) ) {
     nflag = 1;
-    system->local_cap = MAX( (int)(system->n * SAFE_ZONE), MIN_CAP );
+    system->local_cap = MAX( (int)(system->n * safezone), mincap );
   }
 
   Nflag = 0;
   if( system->N >= DANGER_ZONE * system->total_cap ||
       (0 && system->N <= LOOSE_ZONE * system->total_cap) ) {
     Nflag = 1;
-    system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
+    system->total_cap = MAX( (int)(system->N * safezone), mincap );
   }
 
   if( Nflag ) {
@@ -839,7 +854,7 @@ void ReAllocate( reax_system *system, control_params *control,
       }
 
       newsize = static_cast<int>
-        (MAX( realloc->num_far*SAFE_ZONE, MIN_CAP*MIN_NBRS ));
+        (MAX( realloc->num_far*safezone, mincap*MIN_NBRS ));
 #if defined(DEBUG_FOCUS)
       fprintf( stderr, "p%d: reallocating far_nbrs: num_fars=%d, space=%dMB\n",
                system->my_rank, (int)(realloc->num_far*SAFE_ZONE),
@@ -869,7 +884,7 @@ void ReAllocate( reax_system *system, control_params *control,
 #endif
 
     newsize = static_cast<int>
-        (MAX( realloc->Htop*SAFE_ZONE, MIN_CAP*MIN_NBRS ));
+        (MAX( realloc->Htop*safezone, mincap*MIN_NBRS ));
     Reallocate_Matrix( &(workspace->H), system->local_cap,
                        newsize, "H", comm );
     //Deallocate_Matrix( workspace->L );
@@ -886,7 +901,7 @@ void ReAllocate( reax_system *system, control_params *control,
     if( system->numH >= DANGER_ZONE * system->Hcap ||
         (0 && system->numH <= LOOSE_ZONE * system->Hcap) ) {
       Hflag = 1;
-      system->Hcap = int(MAX( system->numH * SAFER_ZONE, MIN_CAP ));
+      system->Hcap = int(MAX( system->numH * saferzone, mincap ));
     }
 
     if( Hflag || realloc->hbonds ) {
@@ -926,7 +941,7 @@ void ReAllocate( reax_system *system, control_params *control,
     if( num_bonds == -1 )
       num_bonds = ((*lists)+BONDS)->num_intrs;
 
-    realloc->num_3body = (int)(MAX(realloc->num_3body*SAFE_ZONE, MIN_3BODIES));
+    realloc->num_3body = (int)(MAX(realloc->num_3body*safezone, MIN_3BODIES));
 
     if( !Make_List( num_bonds, realloc->num_3body, TYP_THREE_BODY,
                     (*lists)+THREE_BODIES, comm ) ) {
diff --git a/src/USER-REAXC/reaxc_defs.h b/src/USER-REAXC/reaxc_defs.h
index 33eb8dc170..16fdfdaead 100644
--- a/src/USER-REAXC/reaxc_defs.h
+++ b/src/USER-REAXC/reaxc_defs.h
@@ -1,3 +1,5 @@
+
+
 /*----------------------------------------------------------------------
   PuReMD - Purdue ReaxFF Molecular Dynamics Program
 
@@ -88,7 +90,7 @@
 #define MIN_NBRS       100
 #define MIN_HENTRIES   100
 #define MAX_BONDS      30
-#define MIN_BONDS      15
+#define MIN_BONDS      25
 #define MIN_HBONDS     25
 #define MIN_3BODIES    1000
 #define MIN_GCELL_POPL 50
diff --git a/src/USER-REAXC/reaxc_forces.cpp b/src/USER-REAXC/reaxc_forces.cpp
index 578f25e895..9fac507f78 100644
--- a/src/USER-REAXC/reaxc_forces.cpp
+++ b/src/USER-REAXC/reaxc_forces.cpp
@@ -194,6 +194,8 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
   reallocate_data *realloc;
   realloc = &(workspace->realloc);
 
+  double saferzone = system->saferzone;
+
   /* bond list */
   if( N > 0 ) {
     bonds = *lists + BONDS;
@@ -226,7 +228,7 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
       Hindex = system->my_atoms[i].Hindex;
       if( Hindex > -1 ) {
         system->my_atoms[i].num_hbonds =
-          (int)(MAX( Num_Entries(Hindex, hbonds)*SAFER_ZONE, MIN_HBONDS ));
+          (int)(MAX( Num_Entries(Hindex, hbonds)*saferzone, MIN_HBONDS ));
 
         //if( Num_Entries(i, hbonds) >=
         //(Start_Index(i+1,hbonds)-Start_Index(i,hbonds))*0.90/*DANGER_ZONE*/){
@@ -721,6 +723,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
   far_neighbor_data *nbr_pj;
   reax_atom *atom_i, *atom_j;
 
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+  double saferzone = system->saferzone;
+
   far_nbrs = *lists + FAR_NBRS;
   *Htop = 0;
   memset( hb_top, 0, sizeof(int) * system->local_cap );
@@ -806,9 +812,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
     }
   }
 
-  *Htop = (int)(MAX( *Htop * SAFE_ZONE, MIN_CAP * MIN_HENTRIES ));
+  *Htop = (int)(MAX( *Htop * safezone, mincap * MIN_HENTRIES ));
   for( i = 0; i < system->n; ++i )
-    hb_top[i] = (int)(MAX( hb_top[i] * SAFER_ZONE, MIN_HBONDS ));
+    hb_top[i] = (int)(MAX( hb_top[i] * saferzone, MIN_HBONDS ));
 
   for( i = 0; i < system->N; ++i ) {
     *num_3body += SQR(bond_top[i]);
diff --git a/src/USER-REAXC/reaxc_init_md.cpp b/src/USER-REAXC/reaxc_init_md.cpp
index 14e0a575ad..0b21bdcab1 100644
--- a/src/USER-REAXC/reaxc_init_md.cpp
+++ b/src/USER-REAXC/reaxc_init_md.cpp
@@ -280,9 +280,14 @@ int Init_System( reax_system *system, control_params *control, char *msg )
   int i;
   reax_atom *atom;
 
-  /* determine the local and total capacity */
-  system->local_cap = MAX( (int)(system->n * SAFE_ZONE), MIN_CAP );
-  system->total_cap = MAX( (int)(system->N * SAFE_ZONE), MIN_CAP );
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+  double saferzone = system->saferzone;
+
+  // determine the local and total capacity
+
+  system->local_cap = MAX( (int)(system->n * safezone), mincap);
+  system->total_cap = MAX( (int)(system->N * safezone), mincap);
 
   /* estimate numH and Hcap */
   system->numH = 0;
@@ -293,7 +298,7 @@ int Init_System( reax_system *system, control_params *control, char *msg )
         atom->Hindex = system->numH++;
       else atom->Hindex = -1;
     }
-  system->Hcap = (int)(MAX( system->numH * SAFER_ZONE, MIN_CAP ));
+  system->Hcap = (int)(MAX( system->numH * saferzone, mincap ));
 
 #if defined(DEBUG_FOCUS)
   fprintf( stderr, "p%d: n=%d local_cap=%d\n",
@@ -639,6 +644,10 @@ int  Init_Lists( reax_system *system, control_params *control,
   int nrecv[MAX_NBRS];
   MPI_Comm comm;
 
+  int mincap = system->mincap;
+  double safezone = system->safezone;
+  double saferzone = system->saferzone;
+
   comm = mpi_data->world;
   bond_top = (int*) calloc( system->total_cap, sizeof(int) );
   hb_top = (int*) calloc( system->local_cap, sizeof(int) );
@@ -652,7 +661,7 @@ int  Init_Lists( reax_system *system, control_params *control,
       system->my_atoms[i].num_hbonds = hb_top[i];
       total_hbonds += hb_top[i];
     }
-    total_hbonds = (int)(MAX( total_hbonds*SAFER_ZONE, MIN_CAP*MIN_HBONDS ));
+    total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
 
     if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
                     *lists+HBONDS, comm ) ) {
@@ -674,7 +683,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     system->my_atoms[i].num_bonds = bond_top[i];
     total_bonds += bond_top[i];
   }
-  bond_cap = (int)(MAX( total_bonds*SAFE_ZONE, MIN_CAP*MIN_BONDS ));
+  bond_cap = (int)(MAX( total_bonds*safezone, mincap*MIN_BONDS ));
 
   if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
                   *lists+BONDS, comm ) ) {
@@ -688,7 +697,7 @@ int  Init_Lists( reax_system *system, control_params *control,
 #endif
 
   /* 3bodies list */
-  cap_3body = (int)(MAX( num_3body*SAFE_ZONE, MIN_3BODIES ));
+  cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
   if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
                   *lists+THREE_BODIES, comm ) ){
     fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h
index 90feca2dd4..91c35352a5 100644
--- a/src/USER-REAXC/reaxc_types.h
+++ b/src/USER-REAXC/reaxc_types.h
@@ -488,6 +488,8 @@ typedef struct
 
   class Pair *pair_ptr;
   int my_bonds;
+  int mincap;
+  real safezone, saferzone;
 
 } reax_system;