diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index 588aa34009750fde01dcbdc69f40e86135600246..7a20fa9fd825806a107f77cf38a3850ba3cb2496 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -230,7 +230,7 @@ void Update_Box_Semi_Isotropic( simulation_box *box, rvec mu )
  * is within the simulation box
  *
  * Assumption: (0, 0, 0) is the minimum box coordinate,
- * the the maximum coordinate is strictly positive in each dimension
+ * the maximum coordinate is strictly positive in each dimension
  *
  * Inputs:
  *  x: current atom position 
@@ -253,30 +253,22 @@ void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_bo
         /* new atomic position after update along dimension i */
         tmp = x[i] + dx[i];
 
-        /* outside of simulation box boundary [0.0, d_i],
+        /* outside of simulation box boundary [0.0, d_i),
          * where d_i is the simulation box length along dimesion i */
-        if ( tmp < 0.0 || tmp > box->box_norms[i] )
+        if ( tmp < 0.0 || tmp >= box->box_norms[i] )
         {
             remapped = FALSE;
 
-            /* re-map the position to be in the range [-d_i, d_i] */
-            if ( tmp <= -box->box_norms[i] || tmp >= box->box_norms[i] )
+            /* re-map the position to be in the range (-d_i, d_i) */
+            if ( tmp <= -1.0 * box->box_norms[i] || tmp >= box->box_norms[i] )
             {
-                if ( tmp >= 0.0 )
-                {
-                    rel_map[i] += ((int) CEIL( tmp / box->box_norms[i] )) - 1;
-                }
-                else
-                {
-                    rel_map[i] += (int) FLOOR( tmp / box->box_norms[i] );
-                }
-
+                rel_map[i] += (int) (tmp / box->box_norms[i]);
                 tmp = FMOD( tmp, box->box_norms[i] );
                 remapped = TRUE;
             }
 
             /* if the new position after re-mapping is negative,
-             * translate back to the range [0, d_i] (i.e., inside the box) */
+             * translate back to the range [0, d_i) (i.e., inside the box) */
             if ( tmp < 0.0 )
             {
                 tmp += box->box_norms[i];
@@ -288,7 +280,7 @@ void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_bo
             }
         }
 
-        assert( tmp >= 0.0 && tmp <= box->box_norms[i] );
+        assert( tmp >= 0.0 && tmp < box->box_norms[i] );
 
         x[i] = tmp;
     }
@@ -327,7 +319,7 @@ void Update_Atom_Position_Non_Periodic( rvec x, rvec dx, ivec rel_map, simulatio
         }
         else if ( tmp > box->box_norms[i] )
         {
-            tmp = box->box_norms[i];
+            tmp = box->box_norms[i] - 1.0E-10;
         }
 
         x[i] = tmp;
@@ -665,17 +657,17 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
     kmax = (int) CEIL( vlist_cut / box->box_norms[2] );
 
     /* assumption: orthogonal coordinates */
-    for ( i = -imax; i <= imax; ++i )
+    for ( i = -1 * imax; i <= imax; ++i )
     {
         d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
         sqr_d_i = SQR( d_i );
 
-        for ( j = -jmax; j <= jmax; ++j )
+        for ( j = -1 * jmax; j <= jmax; ++j )
         {
             d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
             sqr_d_j = SQR( d_j );
 
-            for ( k = -kmax; k <= kmax; ++k )
+            for ( k = -1 * kmax; k <= kmax; ++k )
             {
                 d_k = (x2[2] + k * box->box_norms[2]) - x1[2];
                 sqr_d_k = SQR( d_k );
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index ad4e4fb6feff108d44825fb9bd4f86c905c06cae..f5bfde6686432c6d760dea90b34893b6fbe913b6 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -26,6 +26,11 @@
 #include "vector.h"
 
 
+/* Counts the number of atoms binned to each cell within the grid
+ *
+ * NOTE: this assumes the atom positions are within the 
+ * simulation box boundaries [0, d_i) where d_i is the simulation
+ * box length along a particular dimension */
 static int Estimate_GCell_Population( reax_system * const system )
 {
     int i, j, k, l;
@@ -38,8 +43,10 @@ static int Estimate_GCell_Population( reax_system * const system )
 
     for ( l = 0; l < system->N; l++ )
     {
-        /* NOTE: this assumes the atom positions
-         * are within the simulation box boundaries */
+        assert( system->atoms[l].x[0] >= 0.0 && system->atoms[l].x[0] < system->box.box_norms[0] );
+        assert( system->atoms[l].x[1] >= 0.0 && system->atoms[l].x[1] < system->box.box_norms[1] );
+        assert( system->atoms[l].x[2] >= 0.0 && system->atoms[l].x[2] < system->box.box_norms[2] );
+
         i = (int) (system->atoms[l].x[0] * g->inv_len[0]);
         j = (int) (system->atoms[l].x[1] * g->inv_len[1]);
         k = (int) (system->atoms[l].x[2] * g->inv_len[2]);
@@ -406,15 +413,15 @@ static void Find_Neighbor_Grid_Cells( grid * const g )
                 stack_top = 0;
 
                 /* choose an unmarked neighbor cell */
-                for ( di = -g->spread[0]; di <= g->spread[0]; di++ )
+                for ( di = -1 * g->spread[0]; di <= g->spread[0]; di++ )
                 {
                     x = Shift( i, di, 0, g );
 
-                    for ( dj = -g->spread[1]; dj <= g->spread[1]; dj++ )
+                    for ( dj = -1 * g->spread[1]; dj <= g->spread[1]; dj++ )
                     {
                         y = Shift( j, dj, 1, g );
 
-                        for ( dk = -g->spread[2]; dk <= g->spread[2]; dk++ )
+                        for ( dk = -1 * g->spread[2]; dk <= g->spread[2]; dk++ )
                         {
                             z = Shift( k, dk, 2, g );
 
@@ -614,6 +621,11 @@ void Update_Grid( reax_system * const system )
 }
 
 
+/* Utilize a cell method approach to bin atoms to each cell within the grid
+ *
+ * NOTE: this assumes the atom positions are within the 
+ * simulation box boundaries [0, d_i) where d_i is the simulation
+ * box length along a particular dimension */
 void Bin_Atoms( reax_system * const system, static_storage * const workspace )
 {
     int i, j, k, l;
@@ -626,9 +638,13 @@ void Bin_Atoms( reax_system * const system, static_storage * const workspace )
 
     for ( l = 0; l < system->N; l++ )
     {
-        i = (int)(system->atoms[l].x[0] * g->inv_len[0]);
-        j = (int)(system->atoms[l].x[1] * g->inv_len[1]);
-        k = (int)(system->atoms[l].x[2] * g->inv_len[2]);
+        assert( system->atoms[l].x[0] >= 0.0 && system->atoms[l].x[0] < system->box.box_norms[0] );
+        assert( system->atoms[l].x[1] >= 0.0 && system->atoms[l].x[1] < system->box.box_norms[1] );
+        assert( system->atoms[l].x[2] >= 0.0 && system->atoms[l].x[2] < system->box.box_norms[2] );
+
+        i = (int) (system->atoms[l].x[0] * g->inv_len[0]);
+        j = (int) (system->atoms[l].x[1] * g->inv_len[1]);
+        k = (int) (system->atoms[l].x[2] * g->inv_len[2]);
         /* atom index in grid cell => atom number */
         g->atoms[i][j][k][g->top[i][j][k]] = l;
         /* count of current atoms in this grid cell */
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 7c9fc2e97302986609bd231604c77be11d58dd72..c97368cb6e67a068bbf09184ce6ab26d0f444be0 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -193,24 +193,24 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
     switch ( control->ensemble )
     {
     case NVE:
-        data->N_f = 3 * system->N;
+        data->N_f = 3.0 * system->N;
         *Evolve = &Velocity_Verlet_NVE;
         break;
 
     case bNVT:
-        data->N_f = 3 * system->N;
+        data->N_f = 3.0 * system->N;
         *Evolve = &Velocity_Verlet_Berendsen_NVT;
         break;
 
     case nhNVT:
-        data->N_f = 3 * system->N + 1;
+        data->N_f = 3.0 * system->N + 1.0;
         *Evolve = &Velocity_Verlet_Nose_Hoover_NVT_Klein;
 
         if ( control->restart == FALSE
                 || (control->restart == TRUE && control->random_vel == TRUE) )
         {
             data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin
-                    - data->N_f * K_B * control->T);
+                    - data->N_f * K_B / F_CONV * control->T);
             data->therm.v_xi = data->therm.G_xi * control->dt;
             data->therm.v_xi_old = 0.0;
             data->therm.xi = 0.0;
@@ -222,13 +222,13 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
         fprintf( stderr, "[ERROR] THIS OPTION IS NOT YET IMPLEMENTED! TERMINATING...\n" );
         exit( UNKNOWN_OPTION );
 
-        data->N_f = 3 * system->N + 9;
+        data->N_f = 3.0 * system->N + 9.0;
         *Evolve = &Velocity_Verlet_Berendsen_Isotropic_NPT;
 
         if ( control->restart == FALSE )
         {
             data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
-                    data->N_f * K_B * control->T);
+                    data->N_f * K_B / F_CONV * control->T);
             data->therm.v_xi = data->therm.G_xi * control->dt;
             data->iso_bar.eps = 1.0 / 3.0 * LOG( system->box.volume );
 //            data->inv_W = 1.0 / (data->N_f * K_B * control->T * SQR(control->Tau_P));
@@ -238,13 +238,13 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
 
     /* semi-isotropic NPT */
     case sNPT:
-        data->N_f = 3 * system->N + 4;
+        data->N_f = 3.0 * system->N + 4.0;
         *Evolve = &Velocity_Verlet_Berendsen_Semi_Isotropic_NPT;
         break;
 
     /* isotropic NPT */
     case iNPT:
-        data->N_f = 3 * system->N + 2;
+        data->N_f = 3.0 * system->N + 2.0;
         *Evolve = &Velocity_Verlet_Berendsen_Isotropic_NPT;
         break;
 
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 23e725ba35e7bcae684bfdafa3bc5f24495e7cb2..fc1a95f14865d73a8bc43386922fe170c273ace2 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -32,7 +32,8 @@
 #include "vector.h"
 
 
-/* Velocity Verlet integrator for microcanonical ensemble. */
+/* Perform simulation using the microcanonical ensemble with
+ * the Velocity Verlet integrator. */
 void Velocity_Verlet_NVE( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
@@ -83,8 +84,8 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
 }
 
 
-/* Velocity Verlet integrator for constant volume and temperature
- *  with Berendsen thermostat. */
+/* Perform simulation using constant volume and temperature (NVT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator. */
 void Velocity_Verlet_Berendsen_NVT( reax_system* system,
         control_params* control, simulation_data *data,
         static_storage *workspace, reax_list **lists,
@@ -162,8 +163,8 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
 }
 
 
-/* Velocity Verlet integrator for constant volume and constant temperature
- *  with Nose-Hoover thermostat.
+/* Perform simulation using constant volume and temperature (NVT ensemble)
+ * with a Nose-Hoover thermostat and the Velocity Verlet integrator.
  *
  * Reference: Understanding Molecular Simulation, Frenkel and Smit
  *  Academic Press Inc. San Diego, 1996 p. 388-391 */
@@ -258,7 +259,8 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params*
 }
 
 
-/* Velocity Verlet integrator for constant pressure and constant temperature.
+/* Perform simulation using constant pressure and temperature (NPT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator.
  *
  * NOTE: All box dimensions are scaled by the same amount, and
  * there is no change in the angles between axes. */
@@ -369,7 +371,8 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
 }
 
 
-/* Velocity Verlet integrator for constant pressure and constant temperature. */
+/* Perform simulation using constant pressure and temperature (NPT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator. */
 void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
         control_params* control, simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
@@ -510,7 +513,7 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
     /* Compute zeta(t + dt/2), E_Kininetic(t + dt/2)
      * IMPORTANT: What will be the initial value of zeta? and what is g? */
     data->therm.xi += 0.5 * dt * control->Tau_T
-        * (2.0 * data->E_Kin - data->N_f * K_B * control->T);
+        * (2.0 * data->E_Kin - data->N_f * K_B / F_CONV * control->T);
 
     Reset( system, control, data, workspace );
     fprintf(out_control->log, "reset-");
@@ -544,10 +547,10 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
 
     /* Compute zeta(t + dt) */
     data->therm.xi += 0.5 * dt * control->Tau_T
-        * (2.0 * data->E_Kin - data->N_f * K_B * control->T);
+        * (2.0 * data->E_Kin - data->N_f * K_B / F_CONV * control->T);
 
     fprintf( out_control->log, "Xi: %8.3f %8.3f %8.3f\n",
-             data->therm.xi, data->E_Kin, data->N_f * K_B * control->T );
+             data->therm.xi, data->E_Kin, data->N_f * K_B / F_CONV * control->T );
     fflush( out_control->log );
 }
 
diff --git a/sPuReMD/src/io_tools.c b/sPuReMD/src/io_tools.c
index 819d76781b6da386c26152914d4d4bdd5015beaa..24fda64895773228647cffb54e4b562eeeb2420b 100644
--- a/sPuReMD/src/io_tools.c
+++ b/sPuReMD/src/io_tools.c
@@ -460,8 +460,10 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
 
 
 /* far nbrs contain only i-j nbrhood info, no j-i. */
-void Print_Far_Neighbors( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, reax_list **lists )
+void Print_Far_Neighbors( reax_system const * const system,
+        control_params const * const control,
+        simulation_data const * const data,
+        static_storage const * const workspace, reax_list ** const lists )
 {
     int i, pj, id_i, id_j;
     char fname[MAX_STR];
@@ -482,7 +484,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
             id_j = workspace->orig_id[far_nbrs->far_nbr_list[pj].nbr];
 
             fprintf( fout, "%6d%6d%3d%3d%3d %12.7f\n",
-                     id_i, id_j,
+                     id_j, id_i,
                      far_nbrs->far_nbr_list[pj].rel_box[0],
                      far_nbrs->far_nbr_list[pj].rel_box[1],
                      far_nbrs->far_nbr_list[pj].rel_box[2],
@@ -595,7 +597,7 @@ void Output_Results( reax_system *system, control_params *control,
 #if defined(TEST_ENERGY) || defined(DEBUG) || defined(DEBUG_FOCUS)
         fprintf( out_control->out,
                  "%-6d%24.15f%24.15f%24.15f%13.5f%13.5f%16.5f%13.5f%13.5f\n",
-                 data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
+                 data->step, data->E_Tot, data->E_Pot, data->E_Kin,
                  data->therm.T, control->T, system->box.volume, data->iso_bar.P,
                  (control->P[0] + control->P[1] + control->P[2]) / 3.0 );
 
@@ -610,7 +612,7 @@ void Output_Results( reax_system *system, control_params *control,
 #else
         fprintf( out_control->out,
                  "%-6d%16.2f%16.2f%16.2f%11.2f%11.2f%13.2f%13.5f%13.5f\n",
-                 data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
+                 data->step, data->E_Tot, data->E_Pot, data->E_Kin,
                  data->therm.T, control->T, system->box.volume, data->iso_bar.P,
                  (control->P[0] + control->P[1] + control->P[2]) / 3.0 );
 
diff --git a/sPuReMD/src/io_tools.h b/sPuReMD/src/io_tools.h
index 742a070c1c4a1b24ad4367ed310279f408ef4753..7b37aeb2e03cb351ba102886c31fe5b68d9a7560 100644
--- a/sPuReMD/src/io_tools.h
+++ b/sPuReMD/src/io_tools.h
@@ -32,8 +32,9 @@ char *Get_Atom_Name( reax_system*, int );
 void Print_Near_Neighbors( reax_system*, control_params*, static_storage*,
         reax_list** );
 
-void Print_Far_Neighbors( reax_system*, control_params*, simulation_data*,
-        static_storage*, reax_list** );
+void Print_Far_Neighbors( reax_system const * const, control_params const * const,
+        simulation_data const * const, static_storage const * const,
+        reax_list ** const );
 
 void Print_Total_Force( reax_system*, control_params*, simulation_data*,
         static_storage*, reax_list**, output_controls* );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 4f006123ac87d560ac717a355fe45917764c12ea..e95616c38880cfdd3f5c936ca6e51933b034e5ca 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -73,8 +73,8 @@
   /* unit conversion from ??? to kcal / mol */
   #define C_ELE (332.0638)
   /* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
-  #define K_B (0.831687)
-//  #define K_B (0.8314510)
+//  #define K_B (0.831687)
+  #define K_B (0.8314510)
   /* unit conversion for atomic force to AMU * A / ps^2 */
   #define F_CONV (4.184e2)
   /* energy conversion constant from electron volts to kilo-calories per mole */
@@ -122,7 +122,8 @@
 /* unit conversion for atomic energy */
 #if !defined(E_CONV)
   /* AMU * Angstroms^2 / ps^2 --> kcal / mol */
-  #define E_CONV (0.002391)
+//  #define E_CONV (0.002391)
+  #define E_CONV (1.0 / 418.40)
 #endif
 /* energy conversion constant from electron volts to kilo-calories per mole */
 #if !defined(EV_to_KCALpMOL)
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 05f820003cea38039d2375b91a4e220e04c1ae4a..57a060a80b07512ed1b8987d87dd320af0863e59 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -196,10 +196,10 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
         rvec_Scale( p, m, system->atoms[i].v );
         data->E_Kin += rvec_Dot( p, system->atoms[i].v );
     }
-    data->E_Kin *= 0.5;
+    data->E_Kin *= 0.5 * E_CONV;
 
     /* temperature scalar */
-    data->therm.T = (2.0 * data->E_Kin) / (data->N_f * K_B);
+    data->therm.T = (2.0 * data->E_Kin) / (data->N_f * K_B / F_CONV);
 
     /* avoid T being an absolute zero! */
     if ( FABS( data->therm.T ) < ALMOST_ZERO )
@@ -209,17 +209,18 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
 }
 
 
+/* Compute potential and total energies */
 void Compute_Total_Energy( simulation_data* data )
 {
     data->E_Pot = data->E_BE + data->E_Ov + data->E_Un  + data->E_Lp +
         data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
         data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
 
-
-    data->E_Tot = data->E_Pot + data->E_Kin * E_CONV;
+    data->E_Tot = data->E_Pot + data->E_Kin;
 }
 
 
+/* Check for numeric breakdowns in the energies */
 void Check_Energy( simulation_data* data )
 {
     if ( IS_NAN_REAL(data->E_Pol) )
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 45411f1bcb3a1582698196a00ebcbf32fcb940eb..0483e048cec24801e1bf05a1e0f4851818ce11dc 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -101,7 +101,10 @@ void Transform_to_UnitBox( rvec x1, simulation_box *box, int flag, rvec x2 )
 }
 
 
-/* determine whether point p is inside the box */
+/* Check and remap (if necessary) an atom position to fall
+ * within the boundaries of a periodic simulation box where
+ * the boundaries are [0, d_i) with d_i being the length
+ * of the simulation box in a particular dimension */
 void Fit_to_Periodic_Box( simulation_box *box, rvec p )
 {
     int i;
@@ -116,10 +119,10 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec p )
                 p[i] += box->box_norms[i];
             }
         }
-        else if ( p[i] > box->max[i] )
+        else if ( p[i] >= box->max[i] )
         {
             /* handle higher coords */
-            while ( p[i] > box->max[i] )
+            while ( p[i] >= box->max[i] )
             {
                 p[i] -= box->box_norms[i];
             }
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index a19d4a89450dbf94c4a59c3b35e7d4296fdeabd2..896c040016dfe704fe6c6a2e0872509df3b0700b 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -371,7 +371,7 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
     /* calculate total frame length*/
     snprintf( buffer, SIZE, FRAME_GLOBALS,
              data->step, data->time,
-             data->E_Tot, data->E_Pot, E_CONV * data->E_Kin, data->therm.T,
+             data->E_Tot, data->E_Pot, data->E_Kin, data->therm.T,
              P, system->box.volume,
              system->box.box_norms[0],
              system->box.box_norms[1],
@@ -620,7 +620,7 @@ int Append_xyz_Frame( reax_system *system, control_params *control,
 
     out_control->write( out_control->trj, "%d\t%8.3f\t%8.3f\t%8.3f\t%8.3f\n",
             data->step, data->E_Tot, data->E_Pot,
-            E_CONV * data->E_Kin, data->therm.T );
+            data->E_Kin, data->therm.T );
 
     for ( i = 0; i < system->N; ++i )
     {