From 17f8e293d535ad54eac79608d7b65c4904929918 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Mon, 26 Feb 2018 16:07:50 -0500
Subject: [PATCH] sPuReMD: general code clean-up. Restrict function visibility
 to file scope as applicable.

---
 sPuReMD/src/allocate.c                 |  93 +++++-------
 sPuReMD/src/allocate.h                 |  10 +-
 sPuReMD/src/analyze.c                  |  50 +++---
 sPuReMD/src/analyze.h                  |   2 +
 sPuReMD/src/bond_orders.c              |   5 +-
 sPuReMD/src/bond_orders.h              |   3 +
 sPuReMD/src/box.c                      |   1 +
 sPuReMD/src/box.h                      |  15 +-
 sPuReMD/src/charges.c                  | 190 ++++++-----------------
 sPuReMD/src/charges.h                  |   2 +
 sPuReMD/src/control.c                  |   7 +-
 sPuReMD/src/control.h                  |   3 +-
 sPuReMD/src/ffield.c                   |  21 +--
 sPuReMD/src/ffield.h                   |   5 +-
 sPuReMD/src/forces.c                   | 102 +++----------
 sPuReMD/src/forces.h                   |   9 +-
 sPuReMD/src/four_body_interactions.c   |  24 ++-
 sPuReMD/src/four_body_interactions.h   |   4 +-
 sPuReMD/src/geo_tools.c                | 202 +++++++++++--------------
 sPuReMD/src/geo_tools.h                |  15 +-
 sPuReMD/src/grid.c                     |  20 +--
 sPuReMD/src/grid.h                     |  10 +-
 sPuReMD/src/init_md.c                  |  40 +++--
 sPuReMD/src/integrate.c                |  33 ++--
 sPuReMD/src/lin_alg.c                  | 100 +++---------
 sPuReMD/src/list.c                     |  40 +----
 sPuReMD/src/lookup.c                   |  11 +-
 sPuReMD/src/lookup.h                   |   3 -
 sPuReMD/src/mytypes.h                  |  16 +-
 sPuReMD/src/neighbors.c                |   5 -
 sPuReMD/src/neighbors.h                |   7 +-
 sPuReMD/src/print_utils.c              |   3 +-
 sPuReMD/src/print_utils.h              |  14 +-
 sPuReMD/src/random.c                   |  14 +-
 sPuReMD/src/random.h                   |  10 +-
 sPuReMD/src/reset_utils.c              |   8 +-
 sPuReMD/src/reset_utils.h              |   6 +-
 sPuReMD/src/single_body_interactions.c |  31 ++--
 sPuReMD/src/single_body_interactions.h |   5 +-
 sPuReMD/src/system_props.c             |   1 +
 sPuReMD/src/system_props.h             |   8 +-
 sPuReMD/src/three_body_interactions.c  |   1 +
 sPuReMD/src/three_body_interactions.h  |  12 +-
 sPuReMD/src/traj.h                     |  75 ++++-----
 sPuReMD/src/two_body_interactions.h    |  14 +-
 45 files changed, 490 insertions(+), 760 deletions(-)

diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 6c07fe4f..907fb1a0 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -24,38 +24,37 @@
 #include "list.h"
 #include "tool_box.h"
 
+
 /* allocate space for atoms */
-int PreAllocate_Space( reax_system *system, control_params *control,
+void PreAllocate_Space( reax_system *system, control_params *control,
         static_storage *workspace )
 {
     int i;
 
     system->atoms = (reax_atom*) scalloc( system->N,
-            sizeof(reax_atom), "atoms" );
+            sizeof(reax_atom), "PreAllocate_Space::system->atoms" );
     workspace->orig_id = (int*) scalloc( system->N,
-            sizeof(int), "orid_id" );
+            sizeof(int), "PreAllocate_Space::workspace->orid_id" );
 
     /* space for keeping restriction info, if any */
     if ( control->restrict_bonds )
     {
         workspace->restricted = (int*) scalloc( system->N,
-                sizeof(int), "restricted_atoms" );
+                sizeof(int), "PreAllocate_Space::workspace->restricted_atoms" );
 
         workspace->restricted_list = (int**) scalloc( system->N,
-                sizeof(int*), "restricted_list" );
+                sizeof(int*), "PreAllocate_Space::workspace->restricted_list" );
 
         for ( i = 0; i < system->N; ++i )
         {
             workspace->restricted_list[i] = (int*) scalloc( MAX_RESTRICT,
-                    sizeof(int), "restricted_list[i]" );
+                    sizeof(int), "PreAllocate_Space::workspace->restricted_list[i]" );
         }
     }
-
-    return SUCCESS;
 }
 
 
-void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
+static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
 {
     Delete_List( TYP_FAR_NEIGHBOR, far_nbrs );
 
@@ -70,36 +69,36 @@ void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
 }
 
 
-/* dynamic allocation of memory for matrix in CSR format */
-int Allocate_Matrix( sparse_matrix **pH, int n, int m )
+/* Dynamic allocation of memory for matrix in CSR format
+ *
+ * pH (output): pointer to sparse matrix for which to allocate
+ * n: dimension of the matrix
+ * m: number of nonzeros to allocate space for in matrix
+ * */
+void Allocate_Matrix( sparse_matrix **pH, int n, int m )
 {
     sparse_matrix *H;
 
-    if ( (*pH = (sparse_matrix*) smalloc( sizeof(sparse_matrix),
-                    "Allocate_Matrix::pH" )) == NULL )
-    {
-        return FAILURE;
-    }
+    *pH = (sparse_matrix*) smalloc( sizeof(sparse_matrix),
+            "Allocate_Matrix::pH" );
 
     H = *pH;
     H->n = n;
     H->m = m;
 
-    if ( (H->start = (unsigned int*) smalloc( sizeof(unsigned int) * (n + 1),
-                    "Allocate_Matrix::H->start" )) == NULL
-            || (H->j = (unsigned int*) smalloc( sizeof(unsigned int) * m,
-                    "Allocate_Matrix::H->j" )) == NULL
-            || (H->val = (real*) smalloc( sizeof(real) * m,
-                    "Allocate_Matrix::H->val" )) == NULL )
-    {
-        return FAILURE;
-    }
-
-    return SUCCESS;
+    H->start = (unsigned int*) smalloc( sizeof(unsigned int) * (n + 1),
+            "Allocate_Matrix::H->start" );
+    H->j = (unsigned int*) smalloc( sizeof(unsigned int) * m,
+            "Allocate_Matrix::H->j" );
+    H->val = (real*) smalloc( sizeof(real) * m,
+            "Allocate_Matrix::H->val" );
 }
 
 
-/* deallocate memory for matrix in CSR format */
+/* Deallocate memory for matrix in CSR format
+ *
+ * H (output): pointer to sparse matrix for which to allocate
+ * */
 void Deallocate_Matrix( sparse_matrix *H )
 {
     sfree( H->start, "Deallocate_Matrix::H->start" );
@@ -109,28 +108,15 @@ void Deallocate_Matrix( sparse_matrix *H )
 }
 
 
-int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name )
+static void Reallocate_Matrix( sparse_matrix **H, int n, int m )
 {
     Deallocate_Matrix( *H );
 
-    if ( Allocate_Matrix( H, n, m ) == FAILURE )
-    {
-        fprintf(stderr, "not enough space for %s matrix. terminating!\n", name);
-        exit( INSUFFICIENT_MEMORY );
-    }
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "reallocating %s matrix, n = %d, m = %d\n",
-             name, n, m );
-    fprintf( stderr, "memory allocated: %s = %ldMB\n",
-             name, m * sizeof(sparse_matrix_entry) / (1024 * 1024) );
-#endif
-
-    return SUCCESS;
+    Allocate_Matrix( H, n, m );
 }
 
 
-int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
+void Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
         reax_list *hbonds )
 {
     int i, num_hbonds;
@@ -164,19 +150,14 @@ int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
     fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
              num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
 #endif
-    return SUCCESS;
 }
 
 
-int Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
+static void Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
 {
     int i;
     int *hb_top;
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "reallocating hbonds\n" );
-#endif
-
     hb_top = (int *) scalloc( n, sizeof(int), "Reallocate_HBonds_List::hb_top" );
     for ( i = 0; i < n; ++i )
     {
@@ -191,12 +172,10 @@ int Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
     Allocate_HBond_List( n, num_h, h_index, hb_top, hbonds );
 
     sfree( hb_top, "Reallocate_HBonds_List::hb_top" );
-
-    return SUCCESS;
 }
 
 
-int Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
+void Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
 {
     int i, num_bonds;
 
@@ -223,11 +202,10 @@ int Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
     fprintf( stderr, "memory allocated: bonds = %ldMB\n",
              num_bonds * sizeof(bond_data) / (1024 * 1024) );
 #endif
-    return SUCCESS;
 }
 
 
-int Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
+static void Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
 {
     int i;
     int *bond_top;
@@ -251,8 +229,6 @@ int Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3bo
     *num_bonds = bond_top[n - 1];
 
     sfree( bond_top, "Reallocate_Bonds_List::bond_top" );
-
-    return SUCCESS;
 }
 
 
@@ -276,7 +252,8 @@ void Reallocate( reax_system *system, control_params *control, static_storage *w
 
     if ( realloc->Htop > 0 )
     {
-        Reallocate_Matrix(&(workspace->H), system->N_cm, realloc->Htop * SAFE_ZONE, "H");
+        Reallocate_Matrix( &(workspace->H), system->N_cm,
+                realloc->Htop * SAFE_ZONE );
         realloc->Htop = -1;
 
         Deallocate_Matrix( workspace->L );
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index b961521e..10a401de 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -24,16 +24,18 @@
 
 #include "mytypes.h"
 
-int PreAllocate_Space( reax_system*, control_params*, static_storage* );
+
+void PreAllocate_Space( reax_system*, control_params*, static_storage* );
 
 void Reallocate( reax_system*, control_params*, static_storage*, reax_list**, int );
 
-int Allocate_Matrix( sparse_matrix**, int, int );
+void Allocate_Matrix( sparse_matrix**, int, int );
 
 void Deallocate_Matrix( sparse_matrix* );
 
-int Allocate_HBond_List( int, int, int*, int*, reax_list* );
+void Allocate_HBond_List( int, int, int*, int*, reax_list* );
+
+void Allocate_Bond_List( int, int*, reax_list* );
 
-int Allocate_Bond_List( int, int*, reax_list* );
 
 #endif
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index 1b260eea..3f63f09d 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -60,7 +60,7 @@ typedef struct
 
 
 // copy bond list into old bond list
-void Copy_Bond_List( reax_system *system, control_params *control,
+static void Copy_Bond_List( reax_system *system, control_params *control,
                      reax_list **lists )
 {
     int i, j, top_old;
@@ -91,7 +91,7 @@ void Copy_Bond_List( reax_system *system, control_params *control,
 
 
 // ASSUMPTION: Bond lists are sorted
-int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
+static int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
 {
     int oldp, newp;
     reax_list *new_bonds = &(*lists)[BONDS];
@@ -169,7 +169,7 @@ int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
 }
 
 
-void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
+static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
                    control_params *control, reax_list *bonds, int print,
                    FILE *fout )
 {
@@ -193,7 +193,7 @@ void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
 }
 
 
-void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
+static void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
        size_t size )
 {
     int j, atom;
@@ -220,7 +220,7 @@ void Print_Molecule( reax_system *system, molecule *m, int mode, char *s,
 }
 
 
-void Analyze_Molecules( reax_system *system, control_params *control,
+static void Analyze_Molecules( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, FILE *fout )
 {
@@ -318,7 +318,7 @@ void Analyze_Molecules( reax_system *system, control_params *control,
 }
 
 
-void Report_Bond_Change( reax_system *system, control_params *control,
+static void Report_Bond_Change( reax_system *system, control_params *control,
                          static_storage *workspace,  reax_list *old_bonds,
                          reax_list *new_bonds, int a1, int a2, int flag,
                          FILE *fout )
@@ -389,7 +389,7 @@ void Report_Bond_Change( reax_system *system, control_params *control,
 
 
 /* ASSUMPTION: Bond lists are sorted */
-void Compare_Bonding( int atom, reax_system *system, control_params *control,
+static void Compare_Bonding( int atom, reax_system *system, control_params *control,
                       static_storage *workspace, reax_list *old_bonds,
                       reax_list *new_bonds, FILE *fout )
 {
@@ -507,7 +507,7 @@ void Compare_Bonding( int atom, reax_system *system, control_params *control,
 }
 
 
-void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
+static void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
                   control_params *control, reax_list *bonds, int ignore )
 {
     int i, t, start, end, nbr;
@@ -531,8 +531,7 @@ void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
 }
 
 
-
-void Analyze_Fragments( reax_system *system, control_params *control,
+static void Analyze_Fragments( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, FILE *fout, int ignore )
 {
@@ -601,7 +600,7 @@ void Analyze_Fragments( reax_system *system, control_params *control,
 }
 
 
-void Analyze_Silica( reax_system *system, control_params *control,
+static void Analyze_Silica( reax_system *system, control_params *control,
                      simulation_data *data, static_storage *workspace,
                      reax_list **lists, FILE *fout )
 {
@@ -724,7 +723,7 @@ void Analyze_Silica( reax_system *system, control_params *control,
 }
 
 
-int Get_Type_of_Molecule( molecule *m )
+static int Get_Type_of_Molecule( molecule *m )
 {
     if ( m->atom_count == 3 && m->mtypes[1] == 2 && m->mtypes[2] == 1 )
     {
@@ -735,7 +734,7 @@ int Get_Type_of_Molecule( molecule *m )
 }
 
 
-void Calculate_Dipole_Moment( reax_system *system, control_params *control,
+static void Calculate_Dipole_Moment( reax_system *system, control_params *control,
                               simulation_data *data, static_storage *workspace,
                               reax_list *bonds, FILE *fout )
 {
@@ -781,7 +780,7 @@ void Calculate_Dipole_Moment( reax_system *system, control_params *control,
 }
 
 
-void Copy_Positions( reax_system *system, static_storage *workspace )
+static void Copy_Positions( reax_system *system, static_storage *workspace )
 {
     int i;
 
@@ -790,7 +789,7 @@ void Copy_Positions( reax_system *system, static_storage *workspace )
 }
 
 
-void Calculate_Drift( reax_system *system, control_params *control,
+static void Calculate_Drift( reax_system *system, control_params *control,
                       simulation_data *data, static_storage *workspace,
                       FILE *fout )
 {
@@ -846,7 +845,7 @@ void Calculate_Drift( reax_system *system, control_params *control,
 }
 
 
-void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
+static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
                                FILE *fout )
 {
     int i, j, k;
@@ -913,7 +912,7 @@ void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
 }
 
 
-void Calculate_Density_Slice( reax_system *system, simulation_data *data,
+static void Calculate_Density_Slice( reax_system *system, simulation_data *data,
                               FILE *fout )
 {
     real slice_thickness = 0.5;
@@ -957,14 +956,17 @@ void Analysis( reax_system *system, control_params *control,
     int steps;
 
     steps = data->step - data->prev_steps;
-    // fprintf( stderr, "prev_steps: %d\n", data->prev_steps );
 
     if ( steps == 1 )
     {
         if ( control->molec_anal == REACTIONS )
+        {
             Copy_Bond_List( system, control, lists );
+        }
         if ( control->diffusion_coef )
+        {
             Copy_Positions( system, workspace );
+        }
     }
 
     /****** Molecular Analysis ******/
@@ -977,25 +979,29 @@ void Analysis( reax_system *system, control_params *control,
                                out_control->mol, 0 );
             /* discover fragments without the ignored atoms */
             if ( control->num_ignored )
+            {
                 Analyze_Fragments( system, control, data, workspace, lists,
                                    out_control->ign, 1 );
+            }
         }
         else if ( control->molec_anal == REACTIONS )
+        {
             /* discover molecular changes - reactions */
             Analyze_Molecules( system, control, data, workspace,
                                lists, out_control->mol );
+        }
     }
 
     /****** Electric Dipole Moment ******/
     if ( control->dipole_anal && steps % control->freq_dipole_anal == 0 )
+    {
         Calculate_Dipole_Moment( system, control, data, workspace,
                 &(*lists)[BONDS], out_control->dpl );
+    }
 
     /****** Drift ******/
     if ( control->diffusion_coef && steps % control->freq_diffusion_coef == 0 )
+    {
         Calculate_Drift( system, control, data, workspace, out_control->drft );
-
-#if defined(DEBUG)
-    fprintf( stderr, "analysis... done\n" );
-#endif
+    }
 }
diff --git a/sPuReMD/src/analyze.h b/sPuReMD/src/analyze.h
index af2b056f..c8328240 100644
--- a/sPuReMD/src/analyze.h
+++ b/sPuReMD/src/analyze.h
@@ -24,6 +24,7 @@
 
 #include "mytypes.h"
 
+
 void Analysis( reax_system*, control_params*, simulation_data*,
                static_storage*, reax_list**, output_controls* );
 
@@ -50,4 +51,5 @@ void Analysis( reax_system*, control_params*, simulation_data*,
 
 //void Calculate_Density_Slice( reax_system*, simulation_data*, FILE* );
 
+
 #endif
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index d8f71ee0..7ee01ac1 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -20,6 +20,7 @@
   ----------------------------------------------------------------------*/
 
 #include "bond_orders.h"
+
 #include "list.h"
 #include "lookup.h"
 #include "print_utils.h"
@@ -671,7 +672,7 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
 /* Locate j on i's list.
    This function assumes that j is there for sure!
    And this is the case given our method of neighbor generation*/
-int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
+static int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
 {
     int start = Start_Index(i, bonds);
     int end = End_Index(i, bonds);
@@ -722,7 +723,7 @@ static inline void Copy_Bond_Order_Data( bond_order_data *dest, bond_order_data
 }
 
 
-int compare_bonds( const void *p1, const void *p2 )
+static int compare_bonds( const void *p1, const void *p2 )
 {
     return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
 }
diff --git a/sPuReMD/src/bond_orders.h b/sPuReMD/src/bond_orders.h
index 67592ed8..37cd6242 100644
--- a/sPuReMD/src/bond_orders.h
+++ b/sPuReMD/src/bond_orders.h
@@ -24,6 +24,7 @@
 
 #include "mytypes.h"
 
+
 typedef struct
 {
     real C1dbo, C2dbo, C3dbo;
@@ -32,6 +33,7 @@ typedef struct
     real C1dDelta, C2dDelta, C3dDelta;
 } dbond_coefficients;
 
+
 #ifdef TEST_FORCES
 void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* );
 void Get_dBOpinpi2( reax_system*, reax_list**, int, int, real, real, rvec*, rvec* );
@@ -52,4 +54,5 @@ void Add_dBond_to_Forces_NPT( int, int, reax_system*, simulation_data*,
                               static_storage*, reax_list** );
 void Calculate_Bond_Orders( reax_system*, control_params*, simulation_data*,
                             static_storage*, reax_list**, output_controls* );
+
 #endif
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index f3646ae3..fba8bd32 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -20,6 +20,7 @@
   ----------------------------------------------------------------------*/
 
 #include "box.h"
+
 #include "tool_box.h"
 #include "vector.h"
 
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index 4902a340..8282c14c 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -24,16 +24,17 @@
 
 #include "mytypes.h"
 
-void Setup_Box( real, real, real, real, real, real, simulation_box* );
-
-/* Initializes box from box rtensor */
-void Update_Box(rtensor, simulation_box* /*, int*/);
-void Update_Box_Isotropic(simulation_box*, real /*, int*/);
-void Update_Box_SemiIsotropic( simulation_box*, rvec /*, int*/ );
 
 /* Computes all the transformations,
    metric and other quantities from box rtensor */
-void Make_Consistent(simulation_box*/*, int*/ );
+void Make_Consistent( simulation_box* );
+
+void Setup_Box( real, real, real, real, real, real, simulation_box* );
+
+/* Initializes box from box rtensor */
+void Update_Box( rtensor, simulation_box* );
+void Update_Box_Isotropic( simulation_box*, real );
+void Update_Box_SemiIsotropic( simulation_box*, rvec );
 
 int Are_Far_Neighbors( rvec, rvec, simulation_box*, real, far_neighbor_data* );
 void Get_NonPeriodic_Far_Neighbors( rvec, rvec, simulation_box*,
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index a554cd81..5d2faa75 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -790,25 +790,16 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
 
             if ( workspace->L == NULL )
             {
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
-
+                Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+                Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
             }
             else if ( workspace->L->m < fillin )
             {
                 Deallocate_Matrix( workspace->L );
                 Deallocate_Matrix( workspace->U );
 
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+                Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
             }
             break;
 
@@ -816,12 +807,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -829,12 +816,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -845,12 +828,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
             {
                 /* TODO: safest storage estimate is ILU(0)
                  * (same as lower triangular portion of H), could improve later */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -859,12 +838,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
 
                 /* TODO: safest storage estimate is ILU(0)
                  * (same as lower triangular portion of H), could improve later */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -872,12 +847,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -885,12 +856,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -964,25 +931,16 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
 
             if ( workspace->L == NULL )
             {
-                if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
-
+                Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm );
+                Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm );
             }
             else if ( workspace->L->m < fillin + system->N_cm )
             {
                 Deallocate_Matrix( workspace->L );
                 Deallocate_Matrix( workspace->U );
 
-                if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm );
+                Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm );
             }
             break;
 
@@ -990,12 +948,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1003,12 +957,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -1019,12 +969,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
             {
                 /* TODO: safest storage estimate is ILU(0)
                  * (same as lower triangular portion of H), could improve later */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1033,12 +979,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
 
                 /* TODO: safest storage estimate is ILU(0)
                  * (same as lower triangular portion of H), could improve later */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -1046,12 +988,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1059,8 +997,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
                 {
                     fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
                     exit( INSUFFICIENT_MEMORY );
@@ -1146,12 +1084,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
 
             if ( workspace->L == NULL )
             {
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+                Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
             }
             else if ( workspace->L->m < fillin )
             {
@@ -1159,12 +1093,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, fillin );
+                Allocate_Matrix( &(workspace->U), Hptr->n, fillin );
             }
             break;
 
@@ -1172,12 +1102,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1185,12 +1111,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -1201,12 +1123,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
             {
                 /* TODO: safest storage estimate is ILU(0)
                  * (same as lower triangular portion of H), could improve later */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1214,12 +1132,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
@@ -1227,12 +1141,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
             if ( workspace->L == NULL )
             {
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             else if ( workspace->L->m < Hptr->m )
             {
@@ -1240,12 +1150,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Deallocate_Matrix( workspace->U );
 
                 /* factors have sparsity pattern as H */
-                if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
-                        Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
-                {
-                    fprintf( stderr, "[ERROR] not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m );
+                Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m );
             }
             break;
 
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index e0520fc6..48028b0c 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -24,8 +24,10 @@
 
 #include "mytypes.h"
 
+
 void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
           static_storage* const, const reax_list* const,
           const output_controls* const );
 
+
 #endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 4f5616ab..c48a579d 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -19,14 +19,15 @@
   <http://www.gnu.org/licenses/>.
   ----------------------------------------------------------------------*/
 
+#include "control.h"
+
 #include <ctype.h>
 
-#include "control.h"
 #include "traj.h"
 #include "tool_box.h"
 
 
-char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
+void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
         output_controls *out_control )
 {
     char *s, **tmp;
@@ -599,6 +600,4 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
 
     fprintf(stderr, "control file read\n" );
 #endif
-
-    return SUCCESS;
 }
diff --git a/sPuReMD/src/control.h b/sPuReMD/src/control.h
index 66d0dde7..da79fe52 100644
--- a/sPuReMD/src/control.h
+++ b/sPuReMD/src/control.h
@@ -24,6 +24,7 @@
 
 #include "mytypes.h"
 
-char Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
+
+void Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
 
 #endif
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 1fb05c45..f8b39c8e 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -19,13 +19,14 @@
   <http://www.gnu.org/licenses/>.
   ----------------------------------------------------------------------*/
 
+#include "ffield.h"
+
 #include <ctype.h>
 
-#include "ffield.h"
 #include "tool_box.h"
 
 
-char Read_Force_Field( FILE* fp, reax_interaction* reax )
+void Read_Force_Field( FILE* fp, reax_interaction* reax )
 {
     char *s;
     char **tmp;
@@ -46,7 +47,6 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
     /* reading first header comment */
     fgets( s, MAX_LINE, fp );
 
-
     /* line 2 is number of global parameters */
     fgets( s, MAX_LINE, fp );
     Tokenize( s, &tmp );
@@ -55,8 +55,8 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
     n = atoi(tmp[0]);
     if (n < 1)
     {
-        fprintf( stderr, "WARNING: number of globals in ffield file is 0!\n" );
-        return 1;
+        fprintf( stderr, "[WARNING] number of globals in ffield file is 0!\n" );
+        return;
     }
 
     reax->gp.n_global = n;
@@ -74,19 +74,16 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
         reax->gp.l[i] = val;
     }
 
-
     /* next line is number of atom types and some comments */
     fgets( s, MAX_LINE, fp );
     Tokenize( s, &tmp );
     reax->num_atom_types = atoi(tmp[0]);
 
-
     /* 3 lines of comments */
     fgets(s, MAX_LINE, fp);
     fgets(s, MAX_LINE, fp);
     fgets(s, MAX_LINE, fp);
 
-
     /* Allocating structures in reax_interaction */
     reax->sbp = (single_body_parameters*)
                 scalloc( reax->num_atom_types, sizeof(single_body_parameters),
@@ -743,7 +740,7 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
         k = atoi(tmp[1]) - 1;
         m = atoi(tmp[2]) - 1;
 
-        if (j < reax->num_atom_types && m < reax->num_atom_types)
+        if ( j < reax->num_atom_types && m < reax->num_atom_types )
         {
             val = atof(tmp[3]);
             reax->hbp[j][k][m].r0_hb = val;
@@ -785,10 +782,4 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
     }
 
     sfree( tor_flag, "Read_Force_Field::tor_flag" );
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "force field read\n" );
-#endif
-
-    return SUCCESS;
 }
diff --git a/sPuReMD/src/ffield.h b/sPuReMD/src/ffield.h
index 4aaf32a6..e2e88a10 100644
--- a/sPuReMD/src/ffield.h
+++ b/sPuReMD/src/ffield.h
@@ -23,6 +23,9 @@
 #define __FFIELD_H_
 
 #include "mytypes.h"
-char Read_Force_Field( FILE*, reax_interaction* );
+
+
+void Read_Force_Field( FILE*, reax_interaction* );
+
 
 #endif
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 24e5c291..78783798 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -41,7 +41,7 @@ typedef enum
 } MATRIX_ENTRY_POSITION;
 
 
-void Dummy_Interaction( reax_system *system, control_params *control,
+static void Dummy_Interaction( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
 {
@@ -72,13 +72,11 @@ void Init_Bonded_Force_Functions( control_params *control,
 }
 
 
-void Compute_Bonded_Forces( reax_system *system, control_params *control,
+static void Compute_Bonded_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control, interaction_function *interaction_functions )
 {
-
     int i;
-    //real t_start, t_end, t_elapsed;
 
 #ifdef TEST_ENERGY
     /* Mark beginning of a new timestep in each energy file */
@@ -116,10 +114,6 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
         (interaction_functions[i])( system, control, data, workspace,
                 lists, out_control );
 
-#if defined(DEBUG_FOCUS)
-        fprintf( stderr, "f%d-", i );
-#endif
-
 #ifdef TEST_FORCES
         (Print_Interactions[i])(system, control, data, workspace,
                                 lists, out_control);
@@ -128,7 +122,7 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
 }
 
 
-void Compute_NonBonded_Forces( reax_system *system, control_params *control,
+static void Compute_NonBonded_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list** lists, output_controls *out_control )
 {
@@ -146,10 +140,6 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
     t_elapsed = Get_Timing_Info( t_start );
     data->timing.cm += t_elapsed;
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "qeq - " );
-#endif
-
     if ( control->tabulate <= 0 )
     {
         vdW_Coulomb_Energy( system, control, data, workspace, lists, out_control );
@@ -160,10 +150,6 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
                 lists, out_control );
     }
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "nonb forces - " );
-#endif
-
 #ifdef TEST_FORCES
     Print_vdW_Coulomb_Forces( system, control, data, workspace,
             lists, out_control );
@@ -173,7 +159,7 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
 
 /* This version of Compute_Total_Force computes forces from coefficients
    accumulated by all interaction functions. Saves enormous time & space! */
-void Compute_Total_Force( reax_system *system, control_params *control,
+static void Compute_Total_Force( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace, reax_list **lists )
 {
     int i;
@@ -212,8 +198,6 @@ void Compute_Total_Force( reax_system *system, control_params *control,
         }
 
 #ifdef _OPENMP
-        #pragma omp barrier
-
         #pragma omp for schedule(static)
         for ( i = 0; i < system->N; ++i )
         {
@@ -227,7 +211,7 @@ void Compute_Total_Force( reax_system *system, control_params *control,
 }
 
 
-void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
+static void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
         int Hmax, int Htop, int num_bonds, int num_hbonds )
 {
     int i, flag;
@@ -243,7 +227,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
         if ( Htop > Hmax )
         {
             fprintf( stderr,
-                     "step%d - ran out of space on H matrix: Htop=%d, max = %d",
+                     "[ERROR] step%d - ran out of space on H matrix: Htop=%d, max = %d",
                      step, Htop, Hmax );
             exit( INSUFFICIENT_MEMORY );
         }
@@ -266,7 +250,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
 
     if ( flag > -1 )
     {
-        fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
+        fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
                  step, flag, End_Index(flag, bonds), Start_Index(flag + 1, bonds) );
         exit( INSUFFICIENT_MEMORY );
     }
@@ -277,7 +261,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
 
         if ( End_Index(i, bonds) > bonds->total_intrs )
         {
-            fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
+            fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
                      step, flag, End_Index(i, bonds), bonds->total_intrs );
             exit( INSUFFICIENT_MEMORY );
         }
@@ -304,7 +288,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
 
         if ( flag > -1 )
         {
-            fprintf( stderr, "step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
+            fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
                      step, flag, End_Index(flag, hbonds), Start_Index(flag + 1, hbonds) );
             exit( INSUFFICIENT_MEMORY );
         }
@@ -316,7 +300,7 @@ void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int
 
             if ( End_Index(i, hbonds) > hbonds->total_intrs )
             {
-                fprintf( stderr, "step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
+                fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
                          step, flag, End_Index(i, hbonds), hbonds->total_intrs );
                 exit( INSUFFICIENT_MEMORY );
             }
@@ -366,14 +350,14 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
             break;
 
             default:
-                fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+                fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
                 exit( INVALID_INPUT );
             break;
         }
         break;
 
     default:
-        fprintf( stderr, "Invalid charge method. Terminating...\n" );
+        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
         exit( INVALID_INPUT );
         break;
     }
@@ -419,14 +403,14 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
             break;
 
             default:
-                fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+                fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
                 exit( INVALID_INPUT );
             break;
         }
         break;
 
     default:
-        fprintf( stderr, "Invalid charge method. Terminating...\n" );
+        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
         exit( INVALID_INPUT );
         break;
     }
@@ -602,7 +586,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
 }
 
 
-void Init_Forces( reax_system *system, control_params *control,
+static void Init_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
 {
@@ -835,8 +819,8 @@ void Init_Forces( reax_system *system, control_params *control,
                         bo_ij->BO -= control->bo_cut;
                         bo_ji->BO_s -= control->bo_cut;
                         bo_ji->BO -= control->bo_cut;
-                        workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
-                        workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+                        workspace->total_bond_order[i] += bo_ij->BO;
+                        workspace->total_bond_order[j] += bo_ji->BO;
                         bo_ij->Cdbo = 0.0;
                         bo_ij->Cdbopi = 0.0;
                         bo_ij->Cdbopi2 = 0.0;
@@ -844,35 +828,6 @@ void Init_Forces( reax_system *system, control_params *control,
                         bo_ji->Cdbopi = 0.0;
                         bo_ji->Cdbopi2 = 0.0;
 
-                        /*fprintf( stderr, "%d %d %g %g %g\n",
-                          i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/
-
-                        /*fprintf( stderr, "Cln_BOp_s: %f, pbo2: %f, C12:%f\n",
-                          Cln_BOp_s, twbp->p_bo2, C12 );
-                          fprintf( stderr, "Cln_BOp_pi: %f, pbo4: %f, C34:%f\n",
-                          Cln_BOp_pi, twbp->p_bo4, C34 );
-                          fprintf( stderr, "Cln_BOp_pi2: %f, pbo6: %f, C56:%f\n",
-                          Cln_BOp_pi2, twbp->p_bo6, C56 );*/
-                        /*fprintf(stderr, "pbo1: %f, pbo2:%f\n", twbp->p_bo1, twbp->p_bo2);
-                          fprintf(stderr, "pbo3: %f, pbo4:%f\n", twbp->p_bo3, twbp->p_bo4);
-                          fprintf(stderr, "pbo5: %f, pbo6:%f\n", twbp->p_bo5, twbp->p_bo6);
-                          fprintf( stderr, "r_s: %f, r_p: %f, r_pp: %f\n",
-                          twbp->r_s, twbp->r_p, twbp->r_pp );
-                          fprintf( stderr, "C12: %g, C34:%g, C56:%g\n", C12, C34, C56 );*/
-
-                        /*fprintf( stderr, "\tfactors: %g %g %g\n",
-                          -(bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi +
-                          bo_ij->BO_pi2 * Cln_BOp_pp),
-                          -bo_ij->BO_pi * Cln_BOp_pi, -bo_ij->BO_pi2 * Cln_BOp_pi2 );*/
-                        /*fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
-                          bo_ij->dBOp[0], bo_ij->dBOp[1], bo_ij->dBOp[2] );
-                          fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
-                          bo_ij->dln_BOp_pi[0], bo_ij->dln_BOp_pi[1],
-                          bo_ij->dln_BOp_pi[2] );
-                          fprintf( stderr, "dBOpi2:\t[%g, %g, %g]\n\n",
-                          bo_ij->dln_BOp_pi2[0], bo_ij->dln_BOp_pi2[1],
-                          bo_ij->dln_BOp_pi2[2] );*/
-
                         Set_End_Index( j, btop_j + 1, bonds );
                     }
                 }
@@ -915,7 +870,7 @@ void Init_Forces( reax_system *system, control_params *control,
 }
 
 
-void Init_Forces_Tab( reax_system *system, control_params *control,
+static void Init_Forces_Tab( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
 {
@@ -1148,8 +1103,8 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
                         bo_ij->BO -= control->bo_cut;
                         bo_ji->BO_s -= control->bo_cut;
                         bo_ji->BO -= control->bo_cut;
-                        workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
-                        workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
+                        workspace->total_bond_order[i] += bo_ij->BO;
+                        workspace->total_bond_order[j] += bo_ji->BO;
                         bo_ij->Cdbo = 0.0;
                         bo_ij->Cdbopi = 0.0;
                         bo_ij->Cdbopi2 = 0.0;
@@ -1322,34 +1277,21 @@ void Compute_Forces( reax_system *system, control_params *control,
     t_elapsed = Get_Timing_Info( t_start );
     data->timing.init_forces += t_elapsed;
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "init_forces - ");
-#endif
-
     t_start = Get_Time( );
     Compute_Bonded_Forces( system, control, data, workspace, lists, out_control,
            interaction_functions );
     t_elapsed = Get_Timing_Info( t_start );
     data->timing.bonded += t_elapsed;
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "bonded_forces - ");
-#endif
-
     t_start = Get_Time( );
     Compute_NonBonded_Forces( system, control, data, workspace,
             lists, out_control );
     t_elapsed = Get_Timing_Info( t_start );
     data->timing.nonb += t_elapsed;
 
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "nonbondeds - ");
-#endif
-
     Compute_Total_Force( system, control, data, workspace, lists );
 
 #if defined(DEBUG_FOCUS)
-    fprintf( stderr, "totalforces - ");
     //Print_Total_Force( system, control, data, workspace, lists, out_control );
 #endif
 
@@ -1357,8 +1299,4 @@ void Compute_Forces( reax_system *system, control_params *control,
     Print_Total_Force( system, control, data, workspace, lists, out_control );
     Compare_Total_Forces( system, control, data, workspace, lists, out_control );
 #endif
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "forces - ");
-#endif
 }
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index eeccac22..01c789b9 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -24,13 +24,16 @@
 
 #include "mytypes.h"
 
+
 void Init_Bonded_Force_Functions( control_params*,
        interaction_function* );
 
 void Compute_Forces( reax_system*, control_params*, simulation_data*,
-                     static_storage*, reax_list**, output_controls*,
-                     interaction_function* );
+        static_storage*, reax_list**, output_controls*,
+        interaction_function* );
 
 void Estimate_Storage_Sizes( reax_system*, control_params*, reax_list**,
-                             int*, int*, int*, int* );
+        int*, int*, int*, int* );
+
+
 #endif
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index 22d10476..c46e9bc6 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -20,17 +20,17 @@
   ----------------------------------------------------------------------*/
 
 #include "four_body_interactions.h"
+
 #include "bond_orders.h"
 #include "box.h"
 #include "list.h"
 #include "lookup.h"
 #include "vector.h"
-#include "math.h"
 
 #define MIN_SINE 1e-10
 
 
-real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
+static real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
         rvec dvec_kl, real r_kl, rvec dvec_li, real r_li,
         three_body_interaction_data *p_ijk,
         three_body_interaction_data *p_jkl,
@@ -118,10 +118,22 @@ real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
        -p_jkl->dcos_dk[1]/sin_jkl,
        -p_jkl->dcos_dk[2]/sin_jkl );*/
 
-    if ( sin_ijk >= 0 && sin_ijk <= MIN_SINE ) sin_ijk = MIN_SINE;
-    else if ( sin_ijk <= 0 && sin_ijk >= -MIN_SINE ) sin_ijk = -MIN_SINE;
-    if ( sin_jkl >= 0 && sin_jkl <= MIN_SINE ) sin_jkl = MIN_SINE;
-    else if ( sin_jkl <= 0 && sin_jkl >= -MIN_SINE ) sin_jkl = -MIN_SINE;
+    if ( sin_ijk >= 0 && sin_ijk <= MIN_SINE )
+    {
+        sin_ijk = MIN_SINE;
+    }
+    else if ( sin_ijk <= 0 && sin_ijk >= -MIN_SINE )
+    {
+        sin_ijk = -MIN_SINE;
+    }
+    if ( sin_jkl >= 0 && sin_jkl <= MIN_SINE )
+    {
+        sin_jkl = MIN_SINE;
+    }
+    else if ( sin_jkl <= 0 && sin_jkl >= -MIN_SINE )
+    {
+        sin_jkl = -MIN_SINE;
+    }
 
     // dcos_omega_di
     rvec_ScaledSum( dcos_omega_di, (htra - arg * hnra) / r_ij, dvec_ij, -1., dvec_li );
diff --git a/sPuReMD/src/four_body_interactions.h b/sPuReMD/src/four_body_interactions.h
index 88261c51..e7a8e64f 100644
--- a/sPuReMD/src/four_body_interactions.h
+++ b/sPuReMD/src/four_body_interactions.h
@@ -24,7 +24,9 @@
 
 #include "mytypes.h"
 
+
 void Four_Body_Interactions( reax_system*, control_params*, simulation_data*,
-                             static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
+
 
 #endif
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 549d29ff..05ece8c8 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -27,8 +27,7 @@
 #include "vector.h"
 
 
-/********************* geo format routines ******************/
-void Count_Geo_Atoms( FILE *geo, reax_system *system )
+static void Count_Geo_Atoms( FILE *geo, reax_system *system )
 {
     int i, serial;
     rvec x;
@@ -45,7 +44,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
         Fit_to_Periodic_Box( &(system->box), &x );
     }
 
-    fseek( geo, 0, SEEK_SET ); // set the pointer to the beginning of the file
+    fseek( geo, 0, SEEK_SET );
     fgets( line, MAX_LINE, geo );
     fgets( line, MAX_LINE, geo );
 
@@ -55,75 +54,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
 }
 
 
-char Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
-        simulation_data *data, static_storage *workspace )
-{
-
-    FILE *geo;
-    char descriptor[9];
-    int i, serial, top;
-    real box_x, box_y, box_z, alpha, beta, gamma;
-    rvec x;
-    char element[3], name[9];
-    reax_atom *atom;
-
-    /* open the geometry file */
-    if ( (geo = fopen(geo_file, "r")) == NULL )
-    {
-        fprintf( stderr, "Error opening the geo file! terminating...\n" );
-        exit( FILE_NOT_FOUND );
-    }
-
-    /* read box information */
-    fscanf( geo, CUSTOM_BOXGEO_FORMAT,
-            descriptor, &box_x, &box_y, &box_z, &alpha, &beta, &gamma );
-    /* initialize the box */
-    Setup_Box( box_x, box_y, box_z, alpha, beta, gamma, &(system->box) );
-
-    /* count my atoms & allocate storage */
-    Count_Geo_Atoms( geo, system );
-    if ( PreAllocate_Space( system, control, workspace ) == FAILURE )
-    {
-        fprintf( stderr, "PreAllocate_Space: not enough memory!" );
-        fprintf( stderr, "terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
-
-    /* read in my atom info */
-    top = 0;
-    for ( i = 0; i < system->N; ++i )
-    {
-        fscanf( geo, CUSTOM_ATOM_FORMAT,
-                &serial, element, name, &x[0], &x[1], &x[2] );
-        Fit_to_Periodic_Box( &(system->box), &x );
-#if defined(DEBUG)
-        fprintf( stderr, "atom%d: %s %s %f %f %f\n",
-                 serial, element, name, x[0], x[1], x[2] );
-#endif
-
-        atom = &(system->atoms[top]);
-        workspace->orig_id[i] = serial;
-        atom->type = Get_Atom_Type( &(system->reaxprm), element );
-        strcpy( atom->name, name );
-        rvec_Copy( atom->x, x );
-        rvec_MakeZero( atom->v );
-        rvec_MakeZero( atom->f );
-        atom->q = 0.;
-
-        top++;
-    }
-
-    fclose( geo );
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "finished reading the geo file\n" );
-#endif
-
-    return SUCCESS;
-}
-
-
-int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
+static int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
 {
     int ret;
     char *cryst;
@@ -181,7 +112,64 @@ int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
 }
 
 
-void Count_PDB_Atoms( FILE *geo, reax_system *system )
+void Read_Geo( const char * const geo_file, reax_system* system, control_params *control,
+        simulation_data *data, static_storage *workspace )
+{
+
+    FILE *geo;
+    char descriptor[9];
+    int i, serial, top;
+    real box_x, box_y, box_z, alpha, beta, gamma;
+    rvec x;
+    char element[3], name[9];
+    reax_atom *atom;
+
+    /* open the geometry file */
+    if ( (geo = fopen(geo_file, "r")) == NULL )
+    {
+        fprintf( stderr, "Error opening the geo file! terminating...\n" );
+        exit( FILE_NOT_FOUND );
+    }
+
+    /* read box information */
+    fscanf( geo, CUSTOM_BOXGEO_FORMAT,
+            descriptor, &box_x, &box_y, &box_z, &alpha, &beta, &gamma );
+    /* initialize the box */
+    Setup_Box( box_x, box_y, box_z, alpha, beta, gamma, &(system->box) );
+
+    /* count my atoms & allocate storage */
+    Count_Geo_Atoms( geo, system );
+    PreAllocate_Space( system, control, workspace );
+
+    /* read in my atom info */
+    top = 0;
+    for ( i = 0; i < system->N; ++i )
+    {
+        fscanf( geo, CUSTOM_ATOM_FORMAT,
+                &serial, element, name, &x[0], &x[1], &x[2] );
+        Fit_to_Periodic_Box( &(system->box), &x );
+#if defined(DEBUG)
+        fprintf( stderr, "atom%d: %s %s %f %f %f\n",
+                 serial, element, name, x[0], x[1], x[2] );
+#endif
+
+        atom = &(system->atoms[top]);
+        workspace->orig_id[i] = serial;
+        atom->type = Get_Atom_Type( &(system->reaxprm), element );
+        strcpy( atom->name, name );
+        rvec_Copy( atom->x, x );
+        rvec_MakeZero( atom->v );
+        rvec_MakeZero( atom->f );
+        atom->q = 0.;
+
+        top++;
+    }
+
+    fclose( geo );
+}
+
+
+static void Count_PDB_Atoms( FILE *geo, reax_system *system )
 {
     char *endptr = NULL;
     char line[MAX_LINE + 1];
@@ -221,7 +209,7 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
 }
 
 
-char Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
+void Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
         simulation_data *data, static_storage *workspace )
 {
 
@@ -258,18 +246,9 @@ char Read_PDB( const char * const pdb_file, reax_system* system, control_params
     }
 
     Count_PDB_Atoms( pdb, system );
-    if ( PreAllocate_Space( system, control, workspace ) == FAILURE )
-    {
-        fprintf( stderr, "PreAllocate_Space: not enough memory!" );
-        fprintf( stderr, "terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    PreAllocate_Space( system, control, workspace );
 
     /* start reading and processing the pdb file */
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "starting to read the pdb file\n" );
-#endif
-
     fseek( pdb, 0, SEEK_SET );
     c = 0;
     c1 = 0;
@@ -420,26 +399,21 @@ char Read_PDB( const char * const pdb_file, reax_system* system, control_params
     }
     if ( ferror( pdb ) )
     {
-        return FAILURE;
+        fprintf( stderr, "[ERROR] Unable to read PDB file. Terminating...\n" );
+        exit( INVALID_INPUT );
     }
 
     fclose( pdb );
 
     Deallocate_Tokenizer_Space( &s, &s1, &tmp );
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "finished reading the pdb file\n" );
-#endif
-
-    return SUCCESS;
 } 
 
 
 /* PDB serials are written without regard to the order, we'll see if this
-   cause trouble, if so we'll have to rethink this approach
-   Also, we do not write connect lines yet.
-*/
-char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
+ * cause trouble, if so we'll have to rethink this approach
+ * Also, we do not write connect lines yet.
+ * */
+void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
         control_params *control, static_storage *workspace, output_controls *out_control )
 {
     int i, buffer_req, buffer_len;
@@ -476,7 +450,7 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
                    system->box.box[2][2] * system->box.box[1][2]) /
                   (system->box.box_norms[2] * system->box.box_norms[1]) );
 
-    /*open pdb and write header*/
+    /* open pdb and write header */
     snprintf( fname, MAX_STR + 9, "%s-%d.pdb", control->sim_name, data->step );
     pdb = fopen(fname, "w");
     fprintf( pdb, PDB_CRYST1_FORMAT_O,
@@ -487,7 +461,7 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
     fprintf( out_control->log, "Box written\n" );
     fflush( out_control->log );
 
-    /*write atom lines to buffer*/
+    /* write atom lines to buffer */
     for ( i = 0; i < system->N; i++)
     {
         p_atom = &(system->atoms[i]);
@@ -508,6 +482,12 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
 
     buffer_len = system->N * PDB_ATOM_FORMAT_O_LENGTH;
     buffer[buffer_len] = 0;
+    
+    if ( ferror( pdb ) )
+    {
+        fprintf( stderr, "[ERROR] Unable to write PDB file. Terminating...\n" );
+        exit( INVALID_INPUT );
+    }
 
     fprintf( pdb, "%s", buffer );
     fclose( pdb );
@@ -533,12 +513,10 @@ char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
 
     sfree( buffer, "Write_PDB::buffer" );
     sfree( line, "Write_PDB::line" );
-
-    return SUCCESS;
 }
 
 
-char Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
+void Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
         simulation_data *data, static_storage *workspace )
 {
     FILE *bgf;
@@ -584,17 +562,22 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
         token_cnt = Tokenize( line, &tokens );
 
         if ( !strcmp( tokens[0], "ATOM" ) || !strcmp( tokens[0], "HETATM" ) )
-            (system->N)++;
+        {
+            ++system->N;
+        }
 
         line[0] = 0;
     }
-    if ( ferror ( bgf ) )
+    if ( ferror( bgf ) )
     {
-        return FAILURE;
+        fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+        exit( INVALID_INPUT );
     }
+
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "system->N: %d\n", system->N );
 #endif
+
     fclose( bgf );
 
     /* memory allocations for atoms, atom maps, bond restrictions */
@@ -745,7 +728,7 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
         {
             /* check number of restrictions */
             Check_Input_Range( token_cnt - 2, 0, MAX_RESTRICT,
-                               "CONECT line exceeds max restrictions allowed.\n" );
+                    "CONECT line exceeds max restrictions allowed.\n" );
 
             /* read bond restrictions */
             if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[1]) ) )
@@ -777,16 +760,11 @@ char Read_BGF( const char * const bgf_file, reax_system* system, control_params
             tokens[i][0] = 0;
         }
     }
-    if ( ferror ( bgf ) )
+    if ( ferror( bgf ) )
     {
-        return FAILURE;
+        fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+        exit( INVALID_INPUT );
     }
 
     fclose( bgf );
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "bgf file read\n" );
-#endif
-
-    return SUCCESS;
 }
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index 5e70721d..61cec37a 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -24,14 +24,12 @@
 
 #include "mytypes.h"
 
+
 // CUSTOM_BOXGEO: BOXGEO box_x box_y box_z  angle1 angle2 angle3
 #define CUSTOM_BOXGEO_FORMAT " %s %lf %lf %lf %lf %lf %lf"
 // CUSTOM ATOM: serial element name x y z
 #define CUSTOM_ATOM_FORMAT " %d %s %s %lf %lf %lf"
 
-char Read_Geo( const char * const, reax_system*, control_params*,
-        simulation_data*, static_storage* );
-
 /* PDB format :
 http://www.rcsb.org/pdb/file_formats/pdb/pdbguide2.2/guide2.2_frame.html
 
@@ -117,13 +115,18 @@ COLUMNS       DATA TYPE       FIELD         DEFINITION
 
 #define BGF_CRYSTX_FORMAT "%8s%11s%11s%11s%11s%11s%11s"
 
-char Read_PDB( const char * const, reax_system*, control_params*,
+
+void Read_Geo( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
-char Read_BGF( const char * const, reax_system*, control_params*,
+void Read_PDB( const char * const, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
-char Write_PDB( reax_system*, reax_list*, simulation_data*,
+void Read_BGF( const char * const, reax_system*, control_params*,
+        simulation_data*, static_storage* );
+
+void Write_PDB( reax_system*, reax_list*, simulation_data*,
         control_params*, static_storage*, output_controls* );
 
+
 #endif
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 9c04c3f2..b0749b06 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -26,7 +26,7 @@
 #include "vector.h"
 
 
-int Estimate_GCell_Population( reax_system* system )
+static int Estimate_GCell_Population( reax_system* system )
 {
     int i, j, k, l;
     int max_atoms;
@@ -66,7 +66,7 @@ int Estimate_GCell_Population( reax_system* system )
 }
 
 
-void Allocate_Space_for_Grid( reax_system *system )
+static void Allocate_Space_for_Grid( reax_system *system )
 {
     int i, j, k, l;
     grid *g;
@@ -172,7 +172,7 @@ void Allocate_Space_for_Grid( reax_system *system )
 }
 
 
-void Deallocate_Grid_Space( grid *g )
+static void Deallocate_Grid_Space( grid *g )
 {
     int i, j, k;
 
@@ -216,7 +216,7 @@ void Deallocate_Grid_Space( grid *g )
 }
 
 
-int Shift(int p, int dp, int dim, grid *g )
+int Shift( int p, int dp, int dim, grid *g )
 {
     int dim_len = 0;
     int newp = p + dp;
@@ -248,8 +248,8 @@ int Shift(int p, int dp, int dim, grid *g )
 
 /* finds the closest point between two grid cells denoted by c1 and c2.
    periodic boundary conditions are taken into consideration as well. */
-void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
-                         int c2x, int c2y, int c2z, rvec closest_point )
+static void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
+        int c2x, int c2y, int c2z, rvec closest_point )
 {
     int  i, d;
     ivec c1 = { c1x, c1y, c1z };
@@ -294,7 +294,7 @@ void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
 }
 
 
-void Find_Neighbor_GridCells( grid *g )
+static void Find_Neighbor_GridCells( grid *g )
 {
     int i, j, k;
     int di, dj, dk;
@@ -543,7 +543,7 @@ static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
 }
 
 
-void Copy_Storage( reax_system *system, static_storage *workspace,
+static void Copy_Storage( reax_system *system, static_storage *workspace,
         control_params *control, int top, int old_id, int old_type, int *num_H,
         real **v, real **s, real **t, int *orig_id, rvec *f_old )
 {
@@ -578,7 +578,7 @@ void Copy_Storage( reax_system *system, static_storage *workspace,
 }
 
 
-void Free_Storage( static_storage *workspace, control_params * control )
+static void Free_Storage( static_storage *workspace, control_params * control )
 {
     int i;
 
@@ -600,7 +600,7 @@ void Free_Storage( static_storage *workspace, control_params * control )
 }
 
 
-void Assign_New_Storage( static_storage *workspace,
+static void Assign_New_Storage( static_storage *workspace,
         real **v, real **s, real **t, int *orig_id, rvec *f_old )
 {
     workspace->v = v;
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index 6e83b740..af1cf0ac 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -25,19 +25,17 @@
 #include "mytypes.h"
 
 
+int Shift( int, int, int, grid* );
+
 void Setup_Grid( reax_system* );
 
 void Update_Grid( reax_system* );
 
-void Finalize_Grid( reax_system* );
+void Bin_Atoms( reax_system*, static_storage* );
 
-int  Shift( int, int, int, grid* );
+void Finalize_Grid( reax_system* );
 
 void Cluster_Atoms( reax_system *, static_storage *, control_params * );
 
-void Bin_Atoms( reax_system*, static_storage* );
-
-void Reset_Marks( grid*, ivec*, int );
-
 
 #endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index c809655f..30e9ef4c 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -35,7 +35,7 @@
 #include "vector.h"
 
 
-void Generate_Initial_Velocities( reax_system *system, real T )
+static void Generate_Initial_Velocities( reax_system *system, real T )
 {
     int i;
     real scale, norm;
@@ -73,7 +73,7 @@ void Generate_Initial_Velocities( reax_system *system, real T )
 }
 
 
-void Init_System( reax_system *system, control_params *control,
+static void Init_System( reax_system *system, control_params *control,
         simulation_data *data )
 {
     int i;
@@ -128,7 +128,7 @@ void Init_System( reax_system *system, control_params *control,
 }
 
 
-void Init_Simulation_Data( reax_system *system, control_params *control,
+static void Init_Simulation_Data( reax_system *system, control_params *control,
         simulation_data *data, output_controls *out_control,
         evolve_function *Evolve )
 {
@@ -230,7 +230,7 @@ void Init_Simulation_Data( reax_system *system, control_params *control,
 
 
 /* Initialize Taper params */
-void Init_Taper( control_params *control )
+static void Init_Taper( control_params *control )
 {
     real d1, d7;
     real swa, swa2, swa3;
@@ -273,7 +273,7 @@ void Init_Taper( control_params *control )
 }
 
 
-void Init_Workspace( reax_system *system, control_params *control,
+static void Init_Workspace( reax_system *system, control_params *control,
         static_storage *workspace )
 {
     int i;
@@ -682,7 +682,7 @@ void Init_Workspace( reax_system *system, control_params *control,
 }
 
 
-void Init_Lists( reax_system *system, control_params *control,
+static void Init_Lists( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
 {
@@ -728,20 +728,12 @@ void Init_Lists( reax_system *system, control_params *control,
             break;
     }
 
-    if ( Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz ) == FAILURE )
-    {
-        fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz );
     /* TODO: better estimate for H_sp?
      *   If so, need to refactor Estimate_Storage_Sizes
      *   to use various cut-off distances as parameters
      *   (non-bonded, hydrogen, 3body, etc.) */
-    if ( Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz ) == FAILURE )
-    {
-        fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz );
 
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
@@ -814,7 +806,7 @@ void Init_Lists( reax_system *system, control_params *control,
 }
 
 
-void Init_Out_Controls( reax_system *system, control_params *control,
+static void Init_Out_Controls( reax_system *system, control_params *control,
         static_storage *workspace, output_controls *out_control )
 {
 #define TEMP_SIZE (1000)
@@ -1105,7 +1097,7 @@ void Initialize( reax_system *system, control_params *control,
 }
 
 
-void Finalize_System( reax_system *system, control_params *control,
+static void Finalize_System( reax_system *system, control_params *control,
         simulation_data *data )
 {
     int i, j, k;
@@ -1147,13 +1139,13 @@ void Finalize_System( reax_system *system, control_params *control,
 }
 
 
-void Finalize_Simulation_Data( reax_system *system, control_params *control,
+static void Finalize_Simulation_Data( reax_system *system, control_params *control,
         simulation_data *data, output_controls *out_control )
 {
 }
 
 
-void Finalize_Workspace( reax_system *system, control_params *control,
+static void Finalize_Workspace( reax_system *system, control_params *control,
         static_storage *workspace )
 {
     int i;
@@ -1369,7 +1361,7 @@ void Finalize_Workspace( reax_system *system, control_params *control,
 }
 
 
-void Finalize_Lists( control_params *control, reax_list **lists )
+static void Finalize_Lists( control_params *control, reax_list **lists )
 {
     Delete_List( TYP_FAR_NEIGHBOR, &(*lists)[FAR_NBRS] );
     if ( control->hb_cut > 0.0 )
@@ -1386,7 +1378,7 @@ void Finalize_Lists( control_params *control, reax_list **lists )
 }
 
 
-void Finalize_Out_Controls( reax_system *system, control_params *control,
+static void Finalize_Out_Controls( reax_system *system, control_params *control,
         static_storage *workspace, output_controls *out_control )
 {
     /* close trajectory file */
@@ -1523,7 +1515,6 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
     }
 #endif
 
-
 #ifdef TEST_FORCES
     /* close bond orders file */
     if ( out_control->fbo )
@@ -1594,6 +1585,9 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
 }
 
 
+/* Deallocate top-level data structures, close file handles, etc.
+ *
+ */
 void Finalize( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control, const int output_enabled )
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index aff852c3..fbcc1667 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -92,7 +92,6 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
 }
 
 
-
 void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params* control,
         simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control, interaction_function *interaction_functions )
@@ -204,7 +203,6 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
 }
 
 
-
 /* uses Berendsen-type coupling for both T and P.
    All box dimensions are scaled by the same amount,
    there is no change in the angles between axes. */
@@ -220,6 +218,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
     dt = control->dt;
     steps = data->step - data->prev_steps;
     renbr = (steps % control->reneighbor == 0);
+
 #if defined(DEBUG_FOCUS)
     //fprintf( out_control->prs,
     //         "tau_t: %g  tau_p: %g  dt/tau_t: %g  dt/tau_p: %g\n",
@@ -291,9 +290,13 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
     /* temperature scaler */
     lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
     if ( lambda < MIN_dT )
+    {
         lambda = MIN_dT;
+    }
     else if (lambda > MAX_dT )
+    {
         lambda = MAX_dT;
+    }
     lambda = SQRT( lambda );
 
     /* Scale velocities and positions at t+dt */
@@ -446,13 +449,10 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
 /************************************************/
 
 #ifdef ANISOTROPIC
-
-void Velocity_Verlet_Nose_Hoover_NVT(reax_system* system,
-                                     control_params* control,
-                                     simulation_data *data,
-                                     static_storage *workspace,
-                                     reax_list **lists,
-                                     output_controls *out_control )
+void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
+        control_params* control, simulation_data *data,
+        static_storage *workspace, reax_list **lists,
+        output_controls *out_control,
         interaction_function *interaction_functions )
 {
     int i;
@@ -522,7 +522,6 @@ void Velocity_Verlet_Nose_Hoover_NVT(reax_system* system,
 }
 
 
-
 void Velocity_Verlet_Isotropic_NPT( reax_system* system,
         control_params* control, simulation_data *data,
         static_storage *workspace, reax_list **lists,
@@ -585,7 +584,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
     //Update_Box_Isotropic( EXP( 3.0 * iso_bar->eps ), &(system->box) );
     Update_Box_Isotropic( &(system->box), EXP( 3.0 * iso_bar->eps ) );
 
-
     // Calculate new forces, f(t + dt)
     Reset( system, control, data, workspace );
     fprintf(out_control->log, "reset-");
@@ -604,7 +602,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
     fprintf(out_control->log, "forces\n");
     fflush( out_control->log );
 
-
     // Compute iteration constants for each atom's velocity and for P_internal
     // Compute kinetic energy for initial velocities of the iteration
     P_int_const = E_kin = 0;
@@ -624,7 +621,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
                   rvec_Dot( system->atoms[i].v, system->atoms[i].v ) );
     }
 
-
     // Compute initial p_int
     inv_3V = 1.0 / (3.0 * system->box.volume);
     P_int = inv_3V * ( 2.0 * E_kin + P_int_const );
@@ -640,7 +636,6 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
         v_xi_old = v_xi_new;
         v_eps_old = v_eps_new;
 
-
         for ( i = 0; i < system->N; ++i )
         {
             coef_v = 1.0 / (1.0 + 0.5 * dt * exp_deps *
@@ -655,13 +650,11 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
         v_eps_new = coef_v_eps * ( iso_bar->v_eps +
                                    0.5 * dt * ( iso_bar->a_eps + a_eps_new ) );
 
-
         G_xi_new = control->Tau_T * ( 2.0 * E_kin +
                                       SQR( v_eps_old ) / control->Tau_P -
                                       (data->N_f + 1) * K_B * control->T );
         v_xi_new = therm->v_xi + 0.5 * dt * ( therm->G_xi + G_xi_new );
 
-
         E_kin = 0;
         for ( i = 0; i < system->N; ++i )
             E_kin += (0.5 * system->reaxprm.sbp[system->atoms[i].type].mass *
@@ -669,14 +662,12 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
 
         P_int = inv_3V * ( 2.0 * E_kin + P_int_const );
 
-
         fprintf( out_control->log,
                  "itr %d E_kin: %8.3f veps_n:%8.3f veps_o:%8.3f vxi_n:%8.3f vxi_o: %8.3f\n",
                  itr, E_kin, v_eps_new, v_eps_old, v_xi_new, v_xi_old );
     }
     while ( FABS(v_eps_new - v_eps_old) + FABS(v_xi_new - v_xi_old) > 2e-3 );
 
-
     therm->v_xi_old = therm->v_xi;
     therm->v_xi = v_xi_new;
     therm->G_xi = G_xi_new;
@@ -696,6 +687,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
 
 #endif
 
+
 /* uses Berendsen-type coupling for both T and P.
    All box dimensions are scaled by the same amount,
    there is no change in the angles between axes. */
@@ -710,7 +702,6 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
     rvec dx;
     reax_atom *atom;
 
-
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "step%d\n", data->step );
 #endif
@@ -764,9 +755,13 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
     Compute_Kinetic_Energy( system, data );
     lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
     if ( lambda < MIN_dT )
+    {
         lambda = MIN_dT;
+    }
     else if (lambda > MAX_dT )
+    {
         lambda = MAX_dT;
+    }
     lambda = SQRT( lambda );
 
     fprintf (stderr, "step:%d lambda -> %f \n", data->step, lambda);
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index bb24dfc2..1f17ed77 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -56,11 +56,7 @@ static sparse_matrix * create_test_mat( void )
     unsigned int i, n;
     sparse_matrix *H_test;
 
-    if ( Allocate_Matrix( &H_test, 3, 6 ) == FAILURE )
-    {
-        fprintf( stderr, "not enough memory for test matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &H_test, 3, 6 );
 
     //3x3, SPD, store lower half
     i = 0;
@@ -173,27 +169,15 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
 
     if ( *A_full == NULL )
     {
-        if ( Allocate_Matrix( A_full, A->n, 2 * A->m - A->n ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] not enough memory for full A. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
     }
     else if ( (*A_full)->m < 2 * A->m - A->n )
     {
         Deallocate_Matrix( *A_full );
-        if ( Allocate_Matrix( A_full, A->n, 2 * A->m - A->n ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] not enough memory for full A. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_full, A->n, 2 * A->m - A->n );
     }
 
-    if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "[ERROR] not enough memory for full A. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &A_t, A->n, A->m );
 
     /* Set up the sparse matrix data structure for A. */
     Transpose( A, A_t );
@@ -250,20 +234,12 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
 
     if ( *A_spar_patt == NULL )
     {
-        if ( Allocate_Matrix( A_spar_patt, A->n, A->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] Not enough memory for preconditioning matrices. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_spar_patt, A->n, A->m );
     }
     else if ( ((*A_spar_patt)->m) < (A->m) )
     {
         Deallocate_Matrix( *A_spar_patt );
-        if ( Allocate_Matrix( A_spar_patt, A->n, A->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] Not enough memory for preconditioning matrices. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_spar_patt, A->n, A->m );
     }
 
 
@@ -370,11 +346,7 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
     {
         /* A_app_inv has the same sparsity pattern
          * * as A_spar_patt_full (omit non-zero values) */
-        if ( Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] Not enough memory for approximate inverse matrix. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
     }
     else if ( ((*A_app_inv)->m) < (A->m) )
     {
@@ -382,11 +354,7 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A, sparse_matrix *
 
         /* A_app_inv has the same sparsity pattern
          * * as A_spar_patt_full (omit non-zero values) */
-        if ( Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] Not enough memory for approximate inverse matrix. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( A_app_inv, (*A_spar_patt_full)->n, (*A_spar_patt_full)->m );
     }
 }
 
@@ -617,11 +585,7 @@ real SuperLU_Factorize( const sparse_matrix * const A,
                                     "SuperLU_Factorize::Ltop" );
     Utop = (unsigned int*) smalloc( (A->n + 1) * sizeof(unsigned int),
                                     "SuperLU_Factorize::Utop" );
-    if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &A_t, A->n, A->m );
 
     /* Set up the sparse matrix data structure for A. */
     Transpose( A, A_t );
@@ -1017,11 +981,7 @@ real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
     D = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D" );
     D_inv = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D_inv" );
     Utop = (int*) smalloc( (A->n + 1) * sizeof(int), "ICHOL_PAR::Utop" );
-    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "not enough memory for ICHOL_PAR preconditioning matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &DAD, A->n, A->m );
 
 #ifdef _OPENMP
     #pragma omp parallel for schedule(static) \
@@ -1194,11 +1154,7 @@ real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
 
     D = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D" );
     D_inv = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D_inv" );
-    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "[ILU_PAR] Not enough memory for preconditioning matrices. Terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &DAD, A->n, A->m );
 
 #ifdef _OPENMP
     #pragma omp parallel for schedule(static) \
@@ -1407,13 +1363,9 @@ real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
 
     start = Get_Time( );
 
-    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE ||
-            Allocate_Matrix( &L_temp, A->n, A->m ) == FAILURE ||
-            Allocate_Matrix( &U_temp, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "not enough memory for ILUT_PAR preconditioning matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &DAD, A->n, A->m );
+    Allocate_Matrix( &L_temp, A->n, A->m );
+    Allocate_Matrix( &U_temp, A->n, A->m );
 
     D = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D" );
     D_inv = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D_inv" );
@@ -1976,11 +1928,7 @@ void Transpose_I( sparse_matrix * const A )
 {
     sparse_matrix * A_t;
 
-    if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
-    {
-        fprintf( stderr, "not enough memory for transposing matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &A_t, A->n, A->m );
 
     Transpose( A, A_t );
 
@@ -2505,11 +2453,7 @@ void permute_matrix( const static_storage * const workspace,
     int i, pj, nr, nc;
     sparse_matrix *LUtemp;
 
-    if ( Allocate_Matrix( &LUtemp, LU->n, LU->m ) == FAILURE )
-    {
-        fprintf( stderr, "Not enough space for graph coloring (factor permutation). Terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Matrix( &LUtemp, LU->n, LU->m );
 
     /* count nonzeros in each row of permuted factor (re-use color_top for counting) */
     memset( workspace->color_top, 0, sizeof(unsigned int) * (LU->n + 1) );
@@ -2644,20 +2588,12 @@ void setup_graph_coloring( const control_params * const control,
 {
     if ( *H_p == NULL )
     {
-        if ( Allocate_Matrix( H_p, H->n, H->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] not enough memory for graph coloring. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( H_p, H->n, H->m );
     }
     else if ( (*H_p)->m < H->m )
     {
         Deallocate_Matrix( *H_p );
-        if ( Allocate_Matrix( H_p, H->n, H->m ) == FAILURE )
-        {
-            fprintf( stderr, "[ERROR] not enough memory for graph coloring. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        Allocate_Matrix( H_p, H->n, H->m );
     }
 
     compute_full_sparse_matrix( H, H_full );
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index e6fbd2f4..fb7ddaed 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -97,52 +97,28 @@ void Delete_List( int type, reax_list* l )
     switch ( type )
     {
     case TYP_VOID:
-        if ( l->select.v != NULL )
-        {
-            sfree( l->select.v, "Delete_List::l->select.v" );
-        }
+        sfree( l->select.v, "Delete_List::l->select.v" );
         break;
     case TYP_THREE_BODY:
-        if ( l->select.three_body_list != NULL )
-        {
-            sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
-        }
+        sfree( l->select.three_body_list, "Delete_List::l->select.three_body_list" );
         break;
     case TYP_BOND:
-        if ( l->select.bond_list != NULL )
-        {
-            sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
-        }
+        sfree( l->select.bond_list, "Delete_List::l->select.bond_list" );
         break;
     case TYP_DBO:
-        if ( l->select.dbo_list != NULL )
-        {
-            sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
-        }
+        sfree( l->select.dbo_list, "Delete_List::l->select.dbo_list" );
         break;
     case TYP_DDELTA:
-        if ( l->select.dDelta_list != NULL )
-        {
-            sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
-        }
+        sfree( l->select.dDelta_list, "Delete_List::l->select.dDelta_list" );
         break;
     case TYP_FAR_NEIGHBOR:
-        if ( l->select.far_nbr_list != NULL )
-        {
-            sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
-        }
+        sfree( l->select.far_nbr_list, "Delete_List::l->select.far_nbr_list" );
         break;
     case TYP_NEAR_NEIGHBOR:
-        if ( l->select.near_nbr_list != NULL )
-        {
-            sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
-        }
+        sfree( l->select.near_nbr_list, "Delete_List::l->select.near_nbr_list" );
         break;
     case TYP_HBOND:
-        if ( l->select.hbond_list != NULL )
-        {
-            sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
-        }
+        sfree( l->select.hbond_list, "Delete_List::l->select.hbond_list" );
         break;
 
     default:
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index c1093143..aba004d8 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -52,7 +52,7 @@ void Make_Lookup_Table( real xmin, real xmax, int n,
 
 
 /* Fills solution into x. Warning: will modify c and d! */
-void Tridiagonal_Solve( const real *a, const real *b,
+static void Tridiagonal_Solve( const real *a, const real *b,
         real *c, real *d, real *x, unsigned int n)
 {
     int i;
@@ -77,7 +77,7 @@ void Tridiagonal_Solve( const real *a, const real *b,
 }
 
 
-void Natural_Cubic_Spline( const real *h, const real *f,
+static void Natural_Cubic_Spline( const real *h, const real *f,
         cubic_spline_coef *coef, unsigned int n )
 {
     int i;
@@ -151,8 +151,7 @@ void Natural_Cubic_Spline( const real *h, const real *f,
 }
 
 
-
-void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
+static void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
         cubic_spline_coef *coef, unsigned int n )
 {
     int i;
@@ -219,7 +218,7 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
 }
 
 
-void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
+static void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
 {
     int i;
     real base, dif;
@@ -502,7 +501,7 @@ void Finalize_LR_Lookup_Table( reax_system *system, control_params *control,
 }
 
 
-int Lookup_Index_Of( real x, lookup_table* t )
+static int Lookup_Index_Of( real x, lookup_table* t )
 {
     return (int)( t->a * ( x - t->xmin ) );
 }
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index 6fbd6ffa..804226c8 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -25,14 +25,11 @@
 #include "mytypes.h"
 
 
-/* Function pointer definitions */
 typedef real (*lookup_function)(real);
 
 
 void Make_Lookup_Table( real, real, int, lookup_function, lookup_table* );
 
-int Lookup_Index_Of( real, lookup_table* );
-
 real Lookup( real, lookup_table* );
 
 void Make_LR_Lookup_Table( reax_system*, control_params*,
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 61de424a..00d7d098 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -27,14 +27,14 @@
   #include "config.h"
 #endif
 
-#include "math.h"
-#include "random.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "sys/time.h"
-#include "time.h"
-#include "zlib.h"
+#include <math.h>
+#include <random.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/time.h>
+#include <time.h>
+#include <zlib.h>
 
 #ifdef _OPENMP
   #include <omp.h>
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index b63a7fca..1b9e38a6 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -259,7 +259,6 @@ int Estimate_NumNeighbors( reax_system *system, control_params *control,
 
 
 #if defined DONE
-
 void Choose_Neighbor_Finder( reax_system *system, control_params *control,
                              get_far_neighbors_function *Get_Far_Neighbors )
 {
@@ -620,7 +619,6 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
 }
 
 
-
 void Generate_Neighbor_Lists( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
         reax_list **lists, output_controls *out_control )
@@ -745,7 +743,4 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
              system->N * far_nbrs->intrs_per_unit );
 #endif
 }
-
-
-
 #endif
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index fde27ebe..4239433c 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -24,9 +24,12 @@
 
 #include "mytypes.h"
 
+
 void Generate_Neighbor_Lists( reax_system*, control_params*, simulation_data*,
-                              static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
 
 int Estimate_NumNeighbors( reax_system*, control_params*,
-                           static_storage*, reax_list** );
+        static_storage*, reax_list** );
+
+
 #endif
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 8c3ea309..2ddd9eac 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -506,7 +506,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
 }
 
 
-int fn_qsort_intcmp( const void *a, const void *b )
+static int fn_qsort_intcmp( const void *a, const void *b )
 {
     return ( *(int *)a - * (int *)b);
 }
@@ -745,7 +745,6 @@ void Output_Results( reax_system *system, control_params *control,
 }
 
 
-
 void Print_Linear_System( reax_system *system, control_params *control,
         static_storage *workspace, int step )
 {
diff --git a/sPuReMD/src/print_utils.h b/sPuReMD/src/print_utils.h
index fe75d198..de0fa4ab 100644
--- a/sPuReMD/src/print_utils.h
+++ b/sPuReMD/src/print_utils.h
@@ -24,9 +24,6 @@
 
 #include "mytypes.h"
 
-typedef void (*print_interaction)(reax_system*, control_params*, simulation_data*,
-                                  static_storage*, reax_list**, output_controls*);
-print_interaction Print_Interactions[NO_OF_INTERACTIONS];
 
 char *Get_Element( reax_system*, int );
 
@@ -60,28 +57,39 @@ void Print_Sparse_Matrix2( sparse_matrix*, char*, char* );
 void Print_Sparse_Matrix_Binary( sparse_matrix*, char* );
 
 void Print_Bonds( reax_system*, reax_list*, char* );
+
 void Print_Bond_List2( reax_system*, reax_list*, char* );
 
 #ifdef TEST_FORCES
 void Dummy_Printer( reax_system*, control_params*, simulation_data*,
                     static_storage*, reax_list**, output_controls* );
+
 void Print_Bond_Forces( reax_system*, control_params*, simulation_data*,
                         static_storage*, reax_list**, output_controls* );
+
 void Print_LonePair_Forces( reax_system*, control_params*, simulation_data*,
                             static_storage*, reax_list**, output_controls* );
+
 void Print_OverUnderCoor_Forces(reax_system*, control_params*, simulation_data*,
                                 static_storage*, reax_list**, output_controls*);
+
 void Print_Three_Body_Forces( reax_system*, control_params*, simulation_data*,
                               static_storage*, reax_list**, output_controls* );
+
 void Print_Hydrogen_Bond_Forces(reax_system*, control_params*, simulation_data*,
                                 static_storage*, reax_list**, output_controls*);
+
 void Print_Four_Body_Forces( reax_system*, control_params*, simulation_data*,
                              static_storage*, reax_list**, output_controls* );
+
 void Print_vdW_Coulomb_Forces( reax_system*, control_params*, simulation_data*,
                                static_storage*, reax_list**, output_controls* );
+
 void Compare_Total_Forces( reax_system*, control_params*, simulation_data*,
                            static_storage*, reax_list**, output_controls* );
+
 void Init_Force_Test_Functions( );
 #endif
 
+
 #endif
diff --git a/sPuReMD/src/random.c b/sPuReMD/src/random.c
index 9b09e752..4afe6aea 100644
--- a/sPuReMD/src/random.c
+++ b/sPuReMD/src/random.c
@@ -23,24 +23,26 @@
 
 
 /* System random number generator used linear congruance method with
-   large periodicity for generation of pseudo random number. function
-   Random returns this random number appropriately scaled so that
-   0 <= Random(range) < range */
+ * large periodicity for generation of pseudo random number. function
+ * Random returns this random number appropriately scaled so that
+ * 0 <= Random(range) < range */
 double Random(double range)
 {
     return (random() * range) / 2147483647L;
 }
 
+
 /* This function seeds the system pseudo random number generator with
-   current time. Use this function once in the begining to initialize
-   the system */
+ * current time. Use this function once in the begining to initialize
+ * the system */
 void Randomize()
 {
     srandom(time(NULL));
 }
 
+
 /* GRandom return random number with gaussian distribution with mean
-   and standard deviation "sigma" */
+ * and standard deviation "sigma" */
 double GRandom(double mean, double sigma)
 {
     double v1 = Random(2.0) - 1.0;
diff --git a/sPuReMD/src/random.h b/sPuReMD/src/random.h
index a90b181d..1fd1642f 100644
--- a/sPuReMD/src/random.h
+++ b/sPuReMD/src/random.h
@@ -22,21 +22,23 @@
 #ifndef __RANDOM_H_
 #define __RANDOM_H_
 
-#include <mytypes.h>
+#include "mytypes.h"
+
 
 /* System random number generator used linear congruance method with
    large periodicity for generation of pseudo random number. function
    Random returns this random number appropriately scaled so that
    0 <= Random(range) < range */
-double Random(double);
+double Random( double );
 
 /* This function seeds the system pseudo random number generator with
    current time. Use this function once in the begining to initialize
    the system */
-void Randomize();
+void Randomize( );
 
 /* GRandom return random number with gaussian distribution with mean
    and standard deviation "sigma" */
-double GRandom(double, double);
+double GRandom( double, double );
+
 
 #endif
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index 3edfe907..89f35972 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -36,7 +36,7 @@ void Reset_Atoms( reax_system* system )
 }
 
 
-void Reset_Pressures( simulation_data *data )
+static void Reset_Pressures( simulation_data *data )
 {
     rtensor_MakeZero( data->flex_bar.P );
     data->iso_bar.P = 0.0;
@@ -115,7 +115,7 @@ void Reset_Workspace( reax_system *system, static_storage *workspace )
 
 
 void Reset_Neighbor_Lists( reax_system *system, control_params *control,
-                           static_storage *workspace, reax_list **lists )
+        static_storage *workspace, reax_list **lists )
 {
     int i, tmp;
     reax_list *bonds;
@@ -159,10 +159,6 @@ void Reset( reax_system *system, control_params *control,
     Reset_Workspace( system, workspace );
 
     Reset_Neighbor_Lists( system, control, workspace, lists );
-
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "reset - ");
-#endif
 }
 
 
diff --git a/sPuReMD/src/reset_utils.h b/sPuReMD/src/reset_utils.h
index 4355460b..97598adb 100644
--- a/sPuReMD/src/reset_utils.h
+++ b/sPuReMD/src/reset_utils.h
@@ -24,9 +24,8 @@
 
 #include "mytypes.h"
 
-void Reset_Atoms( reax_system* );
 
-void Reset_Pressures( simulation_data* );
+void Reset_Atoms( reax_system* );
 
 void Reset_Simulation_Data( simulation_data* );
 
@@ -42,10 +41,9 @@ void Reset_Neighbor_Lists( reax_system*, control_params*,
 void Reset( reax_system*, control_params*, simulation_data*,
             static_storage*, reax_list** );
 
-//void Reset_Neighbor_Lists( reax_system*, static_storage*, reax_list** );
-
 void Reset_Grid( grid* );
 
 void Reset_Marks( grid*, ivec*, int );
 
+
 #endif
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 4a9ff310..cf22fb6a 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -92,7 +92,9 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
         /* correction for C2 */
         if ( system->reaxprm.gp.l[5] > 0.001 &&
                 !strcmp( system->reaxprm.sbp[type_i].name, "C" ) )
+        {
             for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
+            {
                 if ( i < bonds->select.bond_list[pj].nbr )
                 {
                     j = bonds->select.bond_list[pj].nbr;
@@ -128,9 +130,10 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
                     }
 
                 }
+            }
+        }
     }
 
-
     for ( i = 0; i < system->N; ++i )
     {
         type_i = system->atoms[i].type;
@@ -138,8 +141,13 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
 
         /* over-coordination energy */
         if ( sbp_i->mass > 21.0 )
+        {
             dfvl = 0.0;
-        else dfvl = 1.0; // only for 1st-row elements
+        }
+        else
+        {
+            dfvl = 1.0; // only for 1st-row elements
+        }
 
         p_ovun2 = sbp_i->p_ovun2;
         sum_ovun1 = 0;
@@ -178,7 +186,6 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
         CEover4 = CEover2 * (dfvl * workspace->Delta_lp_temp[i]) *
                   p_ovun4 * exp_ovun1 * SQR(inv_exp_ovun1);
 
-
         /* under-coordination potential */
         p_ovun2 = sbp_i->p_ovun2;
         p_ovun5 = sbp_i->p_ovun5;
@@ -190,7 +197,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
         inv_exp_ovun8 = 1.0 / (1.0 + exp_ovun8);
 
         data->E_Un += e_un =
-                          -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
+            -p_ovun5 * (1.0 - exp_ovun6) * inv_exp_ovun2n * inv_exp_ovun8;
 
         CEunder1 = inv_exp_ovun2n * ( p_ovun5 * p_ovun6 * exp_ovun6 * inv_exp_ovun8 +
                                       p_ovun2 * e_un * exp_ovun2n);
@@ -211,7 +218,6 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
         Add_dDelta( system, lists, i, CEunder3, workspace->f_un ); // UnCoor - 1st
 #endif
 
-
         for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
         {
             pbond = &(bonds->select.bond_list[pj]);
@@ -220,23 +226,20 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
             bo_ij = &(pbond->bo_data);
             twbp  = &(system->reaxprm.tbp[ type_i ][ type_j ]);
 
-
             bo_ij->Cdbo += CEover1 * twbp->p_ovun1 * twbp->De_s; // OvCoor - 1st
             workspace->CdDelta[j] += CEover4 * (1.0 - dfvl * workspace->dDelta_lp[j]) *
-                                     (bo_ij->BO_pi + bo_ij->BO_pi2); // OvCoor - 3a
+                (bo_ij->BO_pi + bo_ij->BO_pi2); // OvCoor - 3a
             bo_ij->Cdbopi += CEover4 *
-                             (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //OvCoor-3b
+                (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //OvCoor-3b
             bo_ij->Cdbopi2 += CEover4 *
-                              (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //OvCoor-3b
-
+                (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //OvCoor-3b
 
             workspace->CdDelta[j] += CEunder4 * (1.0 - dfvl * workspace->dDelta_lp[j]) *
-                                     (bo_ij->BO_pi + bo_ij->BO_pi2);   // UnCoor - 2a
+                (bo_ij->BO_pi + bo_ij->BO_pi2);   // UnCoor - 2a
             bo_ij->Cdbopi += CEunder4 *
-                             (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
+                (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
             bo_ij->Cdbopi2 += CEunder4 *
-                              (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
-
+                (workspace->Delta[j] - dfvl * workspace->Delta_lp_temp[j]); //UnCoor-2b
 
 #ifdef TEST_ENERGY
             /* fprintf( out_control->eov, "%6d%23.15e%23.15e"
diff --git a/sPuReMD/src/single_body_interactions.h b/sPuReMD/src/single_body_interactions.h
index 8283c624..f575453b 100644
--- a/sPuReMD/src/single_body_interactions.h
+++ b/sPuReMD/src/single_body_interactions.h
@@ -22,9 +22,12 @@
 #ifndef __SINGLE_BODY_INTERACTIONS_H_
 #define __SINGLE_BODY_INTERACTIONS_H_
 
-#include <mytypes.h>
+#include "mytypes.h"
+
 
 void LonePair_OverUnder_Coordination_Energy( reax_system*, control_params*,
         simulation_data*, static_storage*,
         reax_list**, output_controls* );
+
+
 #endif
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index fad302b4..ad98d171 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -20,6 +20,7 @@
   ----------------------------------------------------------------------*/
 
 #include "system_props.h"
+
 #include "tool_box.h"
 #include "vector.h"
 
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index bf8f36e0..996e6a45 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -22,7 +22,7 @@
 #ifndef __SYSTEM_PROP_H_
 #define __SYSTEM_PROP_H_
 
-#include <mytypes.h>
+#include "mytypes.h"
 
 
 void Temperature_Control( control_params*, simulation_data*, output_controls* );
@@ -33,9 +33,11 @@ void Compute_Center_of_Mass( reax_system*, simulation_data*, FILE* );
 
 void Compute_Kinetic_Energy( reax_system*, simulation_data* );
 
-void Compute_Pressure( reax_system*, simulation_data*, static_storage* );
-
 void Compute_Pressure_Isotropic( reax_system*, control_params*, simulation_data*, output_controls* );
 
+void Compute_Pressure_Isotropic_Klein( reax_system*, simulation_data* );
+
+void Compute_Pressure( reax_system*, simulation_data*, static_storage* );
+
 
 #endif
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 485be5fd..1ea18a87 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -20,6 +20,7 @@
   ----------------------------------------------------------------------*/
 
 #include "three_body_interactions.h"
+
 #include "bond_orders.h"
 #include "list.h"
 #include "lookup.h"
diff --git a/sPuReMD/src/three_body_interactions.h b/sPuReMD/src/three_body_interactions.h
index c654751e..30fe5c1c 100644
--- a/sPuReMD/src/three_body_interactions.h
+++ b/sPuReMD/src/three_body_interactions.h
@@ -24,14 +24,16 @@
 
 #include "mytypes.h"
 
-void Three_Body_Interactions( reax_system*, control_params*, simulation_data*,
-                              static_storage*, reax_list**, output_controls* );
-
-void Hydrogen_Bonds( reax_system*, control_params*, simulation_data*,
-                     static_storage*, reax_list**, output_controls* );
 
 void Calculate_Theta( rvec, real, rvec, real, real*, real* );
 
 void Calculate_dCos_Theta( rvec, real, rvec, real, rvec*, rvec*, rvec* );
 
+void Three_Body_Interactions( reax_system*, control_params*, simulation_data*,
+        static_storage*, reax_list**, output_controls* );
+
+void Hydrogen_Bonds( reax_system*, control_params*, simulation_data*,
+        static_storage*, reax_list**, output_controls* );
+
+
 #endif
diff --git a/sPuReMD/src/traj.h b/sPuReMD/src/traj.h
index e45921df..738bc193 100644
--- a/sPuReMD/src/traj.h
+++ b/sPuReMD/src/traj.h
@@ -23,7 +23,8 @@
 #define __TRAJ_H__
 
 #include "mytypes.h"
-#include "zlib.h"
+
+#include <zlib.h>
 
 
 #define BLOCK_MARK "REAX_BLOCK_MARK "
@@ -74,24 +75,28 @@
 #define SIZE_INFO_LINE3 "%-10d %-10d %-10d\n"
 #define SIZE_INFO_LEN3 33
 
-enum ATOM_LINE_OPTS {OPT_NOATOM = 0, OPT_ATOM_BASIC = 4, OPT_ATOM_wF = 5,
-                     OPT_ATOM_wV = 6, OPT_ATOM_FULL = 7
-                    };
-enum BOND_LINE_OPTS {OPT_NOBOND, OPT_BOND_BASIC, OPT_BOND_FULL};
-enum ANGLE_LINE_OPTS {OPT_NOANGLE, OPT_ANGLE_BASIC};
 
-typedef struct
+enum ATOM_LINE_OPTS
 {
-    int no_of_sub_blocks;
-    int size;
-    char* buffer;
-    struct block** sub_blocks;
-} block;
-
+    OPT_NOATOM = 0,
+    OPT_ATOM_BASIC = 4,
+    OPT_ATOM_wF = 5,
+    OPT_ATOM_wV = 6,
+    OPT_ATOM_FULL = 7,
+};
+
+enum BOND_LINE_OPTS
+{
+    OPT_NOBOND = 0,
+    OPT_BOND_BASIC = 1,
+    OPT_BOND_FULL = 2,
+};
 
-int Write_Block( gzFile, block* );
-int Read_Next_Block( gzFile, block*, int* );
-int Skip_Next_Block( gzFile, int*);
+enum ANGLE_LINE_OPTS
+{
+    OPT_NOANGLE = 0,
+    OPT_ANGLE_BASIC = 1,
+};
 
 
 /*
@@ -141,12 +146,11 @@ int Skip_Next_Block( gzFile, int*);
   No. of torsion entries (int)
   Torsion info lines as per torsion format.
 */
-
-
 int Write_Custom_Header( reax_system*, control_params*,
-                         static_storage*, output_controls* );
-int Write_xyz_Header   ( reax_system*, control_params*,
-                         static_storage*, output_controls* );
+        static_storage*, output_controls* );
+
+int Write_xyz_Header( reax_system*, control_params*,
+        static_storage*, output_controls* );
 
 /*
   Write_Traj_Header( gzfile file,
@@ -157,34 +161,13 @@ int Write_xyz_Header   ( reax_system*, control_params*,
  */
 char Write_Traj_Header( FILE*, int, char**, char**, control_params* );
 
-
-/*
-  Push_Traj_Frame(gzfile file,
-                  reax_system* system,
-          control_params* control,
-          simulation_data* data,
-          static_storage* workspace,
-          reax_list** lists,
-          char** various flags);
-*/
-int Push_Traj_Frame( /*gzfile*/ FILE*, reax_system*, control_params*,
-                                simulation_data*, static_storage*, reax_list**, char** );
-
-/*
-  Append_Traj_Frame( gzfile file,
-                        reax_system* system,
-                        control_params* control,
-                simulation_data* data,
-                static_storage* workspace,
-                reax_list** lists,
-                char** various flags);
-*/
 int Append_Custom_Frame( reax_system*, control_params*, simulation_data*,
-                         static_storage*, reax_list**, output_controls* );
-int Append_xyz_Frame   ( reax_system*, control_params*, simulation_data*,
-                         static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
 
+int Append_xyz_Frame( reax_system*, control_params*, simulation_data*,
+        static_storage*, reax_list**, output_controls* );
 
 void Read_Traj( output_controls*, char * );
 
+
 #endif
diff --git a/sPuReMD/src/two_body_interactions.h b/sPuReMD/src/two_body_interactions.h
index b176b5eb..042ee246 100644
--- a/sPuReMD/src/two_body_interactions.h
+++ b/sPuReMD/src/two_body_interactions.h
@@ -22,13 +22,19 @@
 #ifndef __TWO_BODY_INTERACTIONS_H_
 #define __TWO_BODY_INTERACTIONS_H_
 
-#include <mytypes.h>
+#include "mytypes.h"
+
 
 void Bond_Energy( reax_system*, control_params*, simulation_data*,
-                  static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
+
 void vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
-                         static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
+
 void LR_vdW_Coulomb( reax_system*, control_params*, int, int, real, LR_data* );
+
 void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
-                                   static_storage*, reax_list**, output_controls* );
+        static_storage*, reax_list**, output_controls* );
+
+
 #endif
-- 
GitLab