diff --git a/PG-PuReMD/src/box.c b/PG-PuReMD/src/box.c
index 525f24e5cc43f1b50f7bcdaee912b634e2ecfc71..397d5e0daa3b896a8d273f933719747484482264 100644
--- a/PG-PuReMD/src/box.c
+++ b/PG-PuReMD/src/box.c
@@ -204,7 +204,7 @@ void Setup_My_Box( reax_system *system, control_params *control )
     for ( d = 0; d < 3; ++d )
     {
         my_box->min[d] = big_box->box_norms[d] * system->my_coords[d] /
-                         control->procs_by_dim[d];
+            control->procs_by_dim[d];
         my_box->box[d][d] = big_box->box_norms[d] / control->procs_by_dim[d];
         //my_box->max[d] = big_box->box_norms[d] * (system->my_coords[d] + 1) /
         //control->procs_by_dim[d];
diff --git a/PG-PuReMD/src/control.c b/PG-PuReMD/src/control.c
index 5cb8ce6f1099fde86343d81d8c7b147d41bc4864..4469443cdd8e5f626ee604c7a749f1a6831b8b9f 100644
--- a/PG-PuReMD/src/control.c
+++ b/PG-PuReMD/src/control.c
@@ -161,6 +161,12 @@ char Read_Control_File( char *control_file, control_params* control,
         }
         else if ( strcmp(tmp[0], "proc_by_dim") == 0 )
         {
+            if ( c < 4 )
+            {
+                fprintf( stderr, "invalid number of control file parameters (procs_by_dim). terminating!\n" );
+                MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+            }
+
             ival = atoi(tmp[1]);
             control->procs_by_dim[0] = ival;
             ival = atoi(tmp[2]);
diff --git a/PG-PuReMD/src/geo_tools.c b/PG-PuReMD/src/geo_tools.c
index dff292e74019a18c78daa34519bd15b5246d0224..6ac42e6be6afa355ae85d2431131d69e255d4f5c 100644
--- a/PG-PuReMD/src/geo_tools.c
+++ b/PG-PuReMD/src/geo_tools.c
@@ -162,10 +162,10 @@ char Read_Geo( char* geo_file, reax_system* system, control_params *control,
 int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
 {
     char *cryst;
-    char  line[MAX_LINE + 1];
-    char  descriptor[9];
-    char  s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
-    char  s_group[12], s_zValue[12];
+    char line[MAX_LINE + 1];
+    char descriptor[9];
+    char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
+    char s_group[12], s_zValue[12];
 
     /* initialize variables */
     fseek( geo, 0, SEEK_SET ); // set the pointer to the beginning of the file
@@ -181,23 +181,24 @@ int Read_Box_Info( reax_system *system, FILE *geo, int geo_format )
 
     /* locate the cryst line in the geo file, read it and
        initialize the big box */
-    while ( !feof( geo ) )
+    while ( fgets( line, MAX_LINE, geo ) )
     {
-        fgets( line, MAX_LINE, geo );
-
         if ( strncmp( line, cryst, 6 ) == 0 )
         {
             if ( geo_format == PDB )
+            {
                 sscanf( line, PDB_CRYST1_FORMAT,
                         &descriptor[0],
                         &s_a[0], &s_b[0], &s_c[0],
                         &s_alpha[0], &s_beta[0], &s_gamma[0],
                         &s_group[0], &s_zValue[0] );
+            }
 
             /* compute full volume tensor from the angles */
             Setup_Big_Box( atof(s_a),  atof(s_b), atof(s_c),
-                           atof(s_alpha), atof(s_beta), atof(s_gamma),
-                           &(system->big_box) );
+                    atof(s_alpha), atof(s_beta), atof(s_gamma),
+                    &(system->big_box) );
+
             return SUCCESS;
         }
     }
@@ -215,13 +216,13 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
 
     /* initialize variables */
     fseek( geo, 0, SEEK_SET ); /* set the pointer to the beginning of the file */
-    system->bigN = system->n = system->N = 0;
+    system->bigN = 0;
+    system->n = 0;
+    system->N = 0;
 
     /* increment number of atoms for each line denoting an atom desc */
-    while ( !feof( geo ) )
+    while ( fgets( line, MAX_LINE, geo ) )
     {
-        fgets( line, MAX_LINE, geo );
-
         if ( strncmp( line, "ATOM", 4 ) == 0 ||
                 strncmp( line, "HETATM", 6 ) == 0 )
         {
@@ -239,9 +240,13 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
 
             /* if the point is inside my_box, add it to my lists */
             if ( is_Inside_Box(&(system->my_box), x) )
+            {
                 ++system->n;
-            //if( is_Inside_Box(&(system->my_ext_box), x) )
-            //  ++system->N;
+            }
+//            if ( is_Inside_Box(&(system->my_ext_box), x) )
+//            {
+//              ++system->N;
+//            }
         }
     }
 
@@ -307,18 +312,19 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
     fprintf( stderr, "p%d: starting to read the pdb file\n", system->my_rank );
 #endif
     fseek( pdb, 0, SEEK_SET );
-    c  = 0;
+    c = 0;
     c1 = 0;
     top = 0;
-    while ( !feof( pdb ) )
+
+    s[0] = 0;
+    for ( i = 0; i < c1; ++i )
     {
-        /* clear previous input line */
-        s[0] = 0;
-        for ( i = 0; i < c1; ++i )
-            tmp[i][0] = 0;
+        tmp[i][0] = 0;
+    }
 
+    while ( fgets( s, MAX_LINE, pdb ) )
+    {
         /* read new line and tokenize it */
-        fgets( s, MAX_LINE, pdb );
         strncpy( s1, s, MAX_LINE - 1 );
         c1 = Tokenize( s, &tmp );
 
@@ -400,7 +406,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
 
             Fit_to_Periodic_Box( &(system->big_box), &x );
 
-            if ( is_Inside_Box( &(system->my_box), x ) )
+            if ( is_Inside_Box( &(system->my_box), x ) == TRUE )
             {
                 /* store orig_id, type, name and coord info of the new atom */
                 atom = &(system->my_atoms[top]);
@@ -474,6 +480,13 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
                 // fprintf( stderr, "\n" );
             }
         }
+
+        /* clear previous input line */
+        s[0] = 0;
+        for ( i = 0; i < c1; ++i )
+        {
+            tmp[i][0] = 0;
+        }
     }
 
     fclose( pdb );
diff --git a/PG-PuReMD/src/init_md.c b/PG-PuReMD/src/init_md.c
index d071f0dd3e5175ca55bc00d36a266344b41bd059..63e70dc40c5a3cd5f235fdcb502201a2a3a38335 100644
--- a/PG-PuReMD/src/init_md.c
+++ b/PG-PuReMD/src/init_md.c
@@ -427,7 +427,7 @@ void Init_Workspace( reax_system *system, control_params *control,
 
 
 /************** setup communication data structures  **************/
-int Init_MPI_Datatypes( reax_system *system, storage *workspace,
+void Init_MPI_Datatypes( reax_system *system, storage *workspace,
         mpi_datatypes *mpi_data, char *msg )
 {
     int i, block[11];
@@ -442,12 +442,19 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     /* setup the world */
     mpi_data->world = MPI_COMM_WORLD;
 
-    /* mpi_atom - [orig_id, imprt_id, type, num_bonds, num_hbonds, name,
-                   x, v, f_old, s, t] */
-    block[0] = block[1] = block[2] = block[3] = block[4] = 1;
+    /* mpi_atom: orig_id, imprt_id, type, num_bonds, num_hbonds, name,
+     * x, v, f_old, s, t */
+    block[0] = 1;
+    block[1] = 1;
+    block[2] = 1;
+    block[3] = 1;
+    block[4] = 1;
     block[5] = MAX_ATOM_NAME_LEN;
-    block[6] = block[7] = block[8] = 3;
-    block[9] = block[10] = 4;
+    block[6] = 3;
+    block[7] = 3;
+    block[8] = 3;
+    block[9] = 4;
+    block[10] = 4;
 
 //    MPI_Get_address( &sample, &base );
 //    MPI_Get_address( &(sample.orig_id), disp + 0 );
@@ -477,15 +484,27 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     disp[9] = offsetof( mpi_atom, s );
     disp[10] = offsetof( mpi_atom, t );
 
-    type[0] = type[1] = type[2] = type[3] = type[4] = MPI_INT;
+    type[0] = MPI_INT;
+    type[1] = MPI_INT;
+    type[2] = MPI_INT;
+    type[3] = MPI_INT;
+    type[4] = MPI_INT;
     type[5] = MPI_CHAR;
-    type[6] = type[7] = type[8] = type[9] = type[10] = MPI_DOUBLE;
+    type[6] = MPI_DOUBLE;
+    type[7] = MPI_DOUBLE;
+    type[8] = MPI_DOUBLE;
+    type[9] = MPI_DOUBLE;
+    type[10] = MPI_DOUBLE;
 
     MPI_Type_create_struct( 11, block, disp, type, &(mpi_data->mpi_atom_type) );
     MPI_Type_commit( &(mpi_data->mpi_atom_type) );
 
     /* boundary_atom - [orig_id, imprt_id, type, num_bonds, num_hbonds, x] */
-    block[0] = block[1] = block[2] = block[3] = block[4] = 1;
+    block[0] = 1;
+    block[1] = 1;
+    block[2] = 1;
+    block[3] = 1;
+    block[4] = 1;
     block[5] = 3;
 
 //    MPI_Get_address( &b_sample, &base );
@@ -506,7 +525,11 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     disp[4] = offsetof( boundary_atom, num_hbonds );
     disp[5] = offsetof( boundary_atom, x );
 
-    type[0] = type[1] = type[2] = type[3] = type[4] = MPI_INT;
+    type[0] = MPI_INT;
+    type[1] = MPI_INT;
+    type[2] = MPI_INT;
+    type[3] = MPI_INT;
+    type[4] = MPI_INT;
     type[5] = MPI_DOUBLE;
 
     MPI_Type_create_struct( 6, block, disp, type, &(mpi_data->boundary_atom_type) );
@@ -545,9 +568,11 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     MPI_Type_commit( &(mpi_data->mpi_rvec2) );
 
     /* restart_atom - [orig_id, type, name, x, v] */
-    block[0] = block[1] = 1 ;
+    block[0] = 1;
+    block[1] = 1 ;
     block[2] = MAX_ATOM_NAME_LEN;
-    block[3] = block[4] = 3;
+    block[3] = 3;
+    block[4] = 3;
 
 //    MPI_Get_address( &r_sample, &base );
 //    MPI_Get_address( &(r_sample.orig_id), disp + 0 );
@@ -565,14 +590,14 @@ int Init_MPI_Datatypes( reax_system *system, storage *workspace,
     disp[3] = offsetof( restart_atom, x );
     disp[4] = offsetof( restart_atom, v );
 
-    type[0] = type[1] = MPI_INT;
+    type[0] = MPI_INT;
+    type[1] = MPI_INT;
     type[2] = MPI_CHAR;
-    type[3] = type[4] = MPI_DOUBLE;
+    type[3] = MPI_DOUBLE;
+    type[4] = MPI_DOUBLE;
 
     MPI_Type_create_struct( 5, block, disp, type, &(mpi_data->restart_atom_type) );
     MPI_Type_commit( &(mpi_data->restart_atom_type) );
-
-    return SUCCESS;
 }
 
 
@@ -703,14 +728,7 @@ void Initialize( reax_system *system, control_params *control,
 
     char msg[MAX_STR];
 
-    if ( Init_MPI_Datatypes( system, workspace, mpi_data, msg ) == FAILURE )
-    {
-        fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
-                 system->my_rank );
-        fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
-                 system->my_rank );
-        MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
-    }
+    Init_MPI_Datatypes( system, workspace, mpi_data, msg );
 
 #if defined(DEBUG)
     fprintf( stderr, "p%d: initialized mpi datatypes\n", system->my_rank );
@@ -865,14 +883,7 @@ void Initialize( reax_system *system, control_params *control,
     fprintf( stderr, "p%d: initialized workspace\n", system->my_rank );
 #endif
 
-    if ( Init_MPI_Datatypes( system, workspace, mpi_data, msg ) == FAILURE )
-    {
-        fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
-                 system->my_rank );
-        fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
-                 system->my_rank );
-        MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
-    }
+    Init_MPI_Datatypes( system, workspace, mpi_data, msg );
 
 #if defined(DEBUG)
     fprintf( stderr, "p%d: initialized mpi datatypes\n", system->my_rank );
diff --git a/PG-PuReMD/src/init_md.h b/PG-PuReMD/src/init_md.h
index c5222cbd3d2b2a4d2ddef5fe3b976e7d0644ffab..18ec71111056a5e4837680fca0ad03cdc247e622 100644
--- a/PG-PuReMD/src/init_md.h
+++ b/PG-PuReMD/src/init_md.h
@@ -31,7 +31,7 @@ extern "C" {
 
 void Generate_Initial_Velocities( reax_system *, real );
 
-int Init_MPI_Datatypes( reax_system *, storage *, mpi_datatypes *, char * );
+void Init_MPI_Datatypes( reax_system *, storage *, mpi_datatypes *, char * );
 
 void Initialize( reax_system*, control_params*, simulation_data*,
         storage*, reax_list**, output_controls*, mpi_datatypes* );
diff --git a/PG-PuReMD/src/io_tools.c b/PG-PuReMD/src/io_tools.c
index 3cce1777cfea5efff1793db2aa1f7019471b9a04..77b222da1bbfba0d55cc147f4aa6d3ae5bbe7bca 100644
--- a/PG-PuReMD/src/io_tools.c
+++ b/PG-PuReMD/src/io_tools.c
@@ -58,7 +58,9 @@ int Init_Output_Files( reax_system *system, control_params *control,
     {
         ret = Init_Traj( system, control, out_control, mpi_data, msg );
         if ( ret == FAILURE )
+        {
             return ret;
+        }
     }
 
     if ( system->my_rank == MASTER_NODE )
@@ -470,30 +472,58 @@ int Init_Output_Files( reax_system *system, control_params *control,
 
 
 /************************ close output files ************************/
-int Close_Output_Files( reax_system *system, control_params *control,
+void Close_Output_Files( reax_system *system, control_params *control,
         output_controls *out_control, mpi_datatypes *mpi_data )
 {
     if ( out_control->write_steps > 0 )
+    {
         End_Traj( system->my_rank, out_control );
+    }
 
     if ( system->my_rank == MASTER_NODE )
     {
         if ( out_control->energy_update_freq > 0 )
         {
-            fclose( out_control->out );
-            fclose( out_control->pot );
+            if ( out_control->out )
+            {
+                fclose( out_control->out );
+            }
+            if ( out_control->pot )
+            {
+                fclose( out_control->pot );
+            }
+
 #if defined(LOG_PERFORMANCE)
-            fclose( out_control->log );
+            if ( out_control->log )
+            {
+                fclose( out_control->log );
+            }
 #endif
         }
 
         if ( control->ensemble == NPT || control->ensemble == iNPT ||
                 control->ensemble == sNPT )
-            fclose( out_control->prs );
+        {
+            if ( out_control->prs )
+            {
+                fclose( out_control->prs );
+            }
+        }
 
-        if ( control->dipole_anal ) fclose( out_control->dpl );
-        if ( control->diffusion_coef ) fclose( out_control->drft );
-        if ( control->molecular_analysis ) fclose( out_control->mol );
+        if ( control->dipole_anal )
+        {
+            fclose( out_control->dpl );
+        }
+
+        if ( control->diffusion_coef )
+        {
+            fclose( out_control->drft );
+        }
+
+        if ( control->molecular_analysis )
+        {
+            fclose( out_control->mol );
+        }
     }
 
 #ifdef TEST_ENERGY
@@ -541,8 +571,6 @@ int Close_Output_Files( reax_system *system, control_params *control,
     fclose( out_control->nlist );
 #endif
 #endif
-
-    return SUCCESS;
 }
 
 
@@ -673,8 +701,11 @@ void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs )
 
     /* loop over neighbor processes */
     for ( r[0] = -1; r[0] <= 1; ++r[0])
+    {
         for ( r[1] = -1; r[1] <= 1; ++r[1] )
+        {
             for ( r[2] = -1; r[2] <= 1; ++r[2] )
+            {
                 if ( (nbr = Relative_Coord_Encoding( r )) != MYSELF )
                 {
                     nbr_pr = &(my_nbrs[nbr]);
@@ -696,8 +727,11 @@ void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs )
                              nbr_pr->end_recv[0], nbr_pr->end_recv[1],
                              nbr_pr->end_recv[2] );
                 }
+            }
+        }
+    }
 
-    fclose(f);
+    fclose( f );
 }
 
 
@@ -719,7 +753,9 @@ void Print_Native_GCells( reax_system *system )
     g = &(system->my_grid);
 
     for ( i = g->native_str[0]; i < g->native_end[0]; i++ )
+    {
         for ( j = g->native_str[1]; j < g->native_end[1]; j++ )
+        {
             for ( k = g->native_str[2]; k < g->native_end[2]; k++ )
             {
                 gc = &( g->cells[ index_grid_3d(i, j, k, g) ] );
@@ -733,20 +769,24 @@ void Print_Native_GCells( reax_system *system )
 
                 //for( l = gc->str; l < gc->end; ++l )
                 for ( l = g->str[index_grid_3d(i, j, k, g)]; l < g->end[index_grid_3d(i, j, k, g)]; ++l )
+                {
                     fprintf( f, "%5d", system->my_atoms[l].orig_id );
+                }
                 fprintf( f, "\n" );
             }
+        }
+    }
 
-    fclose(f);
+    fclose( f );
 }
 
 
 void Print_All_GCells( reax_system *system )
 {
-    int        i, j, k, l;
-    char       fname[100];
-    FILE      *f;
-    grid      *g;
+    int i, j, k, l;
+    char fname[100];
+    FILE *f;
+    grid *g;
     grid_cell *gc;
     char gcell_type_text[10][12] =
     {
@@ -759,7 +799,9 @@ void Print_All_GCells( reax_system *system )
     g = &(system->my_grid);
 
     for ( i = 0; i < g->ncells[0]; i++ )
+    {
         for ( j = 0; j < g->ncells[1]; j++ )
+        {
             for ( k = 0; k < g->ncells[2]; k++ )
             {
                 gc = &( g->cells[ index_grid_3d(i, j, k, g) ] );
@@ -773,11 +815,15 @@ void Print_All_GCells( reax_system *system )
 
                 //for( l = gc->str; l < gc->end; ++l )
                 for ( l = g->str[index_grid_3d(i, j, k, g)]; l < g->end[index_grid_3d(i, j, k, g)]; ++l )
+                {
                     fprintf( f, "%5d", system->my_atoms[l].orig_id );
+                }
                 fprintf( f, "\n" );
             }
+        }
+    }
 
-    fclose(f);
+    fclose( f );
 }
 
 
@@ -798,12 +844,14 @@ void Print_My_Atoms( reax_system *system )
     //   system->my_rank, system->n );
 
     for ( i = 0; i < system->n; ++i )
+    {
         fprintf( fh, "p%-2d %-5d %2d %24.15e%24.15e%24.15e\n",
                  system->my_rank,
                  system->my_atoms[i].orig_id, system->my_atoms[i].type,
                  system->my_atoms[i].x[0],
                  system->my_atoms[i].x[1],
                  system->my_atoms[i].x[2] );
+    }
 
     fclose( fh );
 }
@@ -826,12 +874,14 @@ void Print_My_Ext_Atoms( reax_system *system )
     //   system->my_rank, system->n );
 
     for ( i = 0; i < system->N; ++i )
+    {
         fprintf( fh, "p%-2d %-5d imprt%-5d %2d %24.15e%24.15e%24.15e\n",
                  system->my_rank, system->my_atoms[i].orig_id,
                  system->my_atoms[i].imprt_id, system->my_atoms[i].type,
                  system->my_atoms[i].x[0],
                  system->my_atoms[i].x[1],
                  system->my_atoms[i].x[2] );
+    }
 
     fclose( fh );
 }
@@ -840,13 +890,13 @@ void Print_My_Ext_Atoms( reax_system *system )
 void Print_Far_Neighbors( reax_system *system, reax_list **lists,
         control_params *control )
 {
-    char  fname[100];
-    int   i, j, id_i, id_j, nbr, natoms;
+    char fname[100];
+    int i, j, id_i, id_j, nbr, natoms;
     FILE *fout;
     reax_list *far_nbrs;
 
     sprintf( fname, "%s.far_nbrs.%d", control->sim_name, system->my_rank );
-    fout      = fopen( fname, "w" );
+    fout = fopen( fname, "w" );
     far_nbrs = (*lists) + FAR_NBRS;
     natoms = system->N;
 
@@ -996,14 +1046,16 @@ void Print_LinSys_Soln( reax_system *system, real *x, real *b_prm, real *b )
 {
     int i;
     char fname[100];
-    FILE  *fout;
+    FILE *fout;
 
     sprintf( fname, "qeq.%d.out", system->my_rank );
     fout = fopen( fname, "w" );
 
     for ( i = 0; i < system->n; ++i )
+    {
         fprintf( fout, "%6d%10.4f%10.4f%10.4f\n",
                  system->my_atoms[i].orig_id, x[i], b_prm[i], b[i] );
+    }
 
     fclose( fout );
 }
@@ -1011,19 +1063,21 @@ void Print_LinSys_Soln( reax_system *system, real *x, real *b_prm, real *b )
 
 void Print_Charges( reax_system *system )
 {
-    int    i;
-    char   fname[100];
-    FILE  *fout;
+    int i;
+    char fname[100];
+    FILE *fout;
 
     sprintf( fname, "q.%d.out", system->my_rank );
     fout = fopen( fname, "w" );
 
     for ( i = 0; i < system->n; ++i )
+    {
         fprintf( fout, "%6d %10.7f %10.7f %10.7f\n",
                  system->my_atoms[i].orig_id,
                  system->my_atoms[i].s[0],
                  system->my_atoms[i].t[0],
                  system->my_atoms[i].q );
+    }
 
     fclose( fout );
 }
@@ -1134,12 +1188,16 @@ void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
             nbr = bonds->select.bond_list[pj].nbr;
             id_j = system->my_atoms[nbr].orig_id;
             if ( id_i < id_j )
+            {
                 temp[num++] = id_j;
+            }
         }
 
-        qsort(&temp, num, sizeof(int), fn_qsort_intcmp);
-        for (j = 0; j < num; j++)
+        qsort( &temp, num, sizeof(int), fn_qsort_intcmp );
+        for ( j = 0; j < num; j++ )
+        {
             fprintf(f, "%6d", temp[j] );
+        }
         fprintf(f, "\n");
     }
 }
@@ -1255,7 +1313,7 @@ void Output_Results( reax_system *system, control_params *control,
                          data->tot_press[0], data->tot_press[1], data->tot_press[2],
                          system->big_box.V );
 
-                fflush( out_control->prs);
+                fflush( out_control->prs );
             }
 
             fflush( out_control->out );
diff --git a/PG-PuReMD/src/io_tools.h b/PG-PuReMD/src/io_tools.h
index f83c9686400a03a70e35a06f777b1a9965e0f02e..cb69ad0501a020bccf6bdd21f4f5f8c1c0054c94 100644
--- a/PG-PuReMD/src/io_tools.h
+++ b/PG-PuReMD/src/io_tools.h
@@ -32,7 +32,7 @@ extern "C" {
 int Init_Output_Files( reax_system*, control_params*,
         output_controls*, mpi_datatypes*, char* );
 
-int Close_Output_Files( reax_system*, control_params*,
+void Close_Output_Files( reax_system*, control_params*,
         output_controls*, mpi_datatypes* );
 
 void Print_Box( simulation_box*, char*, FILE* );
diff --git a/PG-PuReMD/src/parallelreax.c b/PG-PuReMD/src/parallelreax.c
index da9074127fec719403a9ea0e892961314c7a7b29..3b3757f186b619540d9cd18cc117ffca0738d028 100644
--- a/PG-PuReMD/src/parallelreax.c
+++ b/PG-PuReMD/src/parallelreax.c
@@ -183,10 +183,12 @@ int main( int argc, char* argv[] )
     real t_begin, t_end;
 #endif
 
+    MPI_Init( &argc, &argv );
+
     if ( argc != 4 )
     {
         usage( argv );
-        exit( INVALID_INPUT );
+        MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
     }
 
 #ifdef HAVE_CUDA
@@ -228,7 +230,6 @@ int main( int argc, char* argv[] )
     }
 
     /* setup MPI environment */
-    MPI_Init( &argc, &argv );
     MPI_Comm_size( MPI_COMM_WORLD, &(control->nprocs) );
     MPI_Comm_rank( MPI_COMM_WORLD, &(system->my_rank) );
 
@@ -443,7 +444,6 @@ int main( int argc, char* argv[] )
     }
 
     /* setup MPI environment */
-    MPI_Init( &argc, &argv );
     MPI_Comm_size( MPI_COMM_WORLD, &(control->nprocs) );
     MPI_Comm_rank( MPI_COMM_WORLD, &(system->my_rank) );
 
@@ -561,17 +561,6 @@ int main( int argc, char* argv[] )
 //    Write_PDB( &system, &(lists[BOND]), &out_control );
     Close_Output_Files( system, control, out_control, mpi_data );
 
-    MPI_Finalize( );
-
-    /* deallocate data structures */
-    sfree( system, "main::system" );
-    sfree( control, "main::control" );
-    sfree( data, "main::data" );
-    sfree( workspace, "main::workspace" );
-    sfree( lists, "main::lists" );
-    sfree( out_control, "main::out_control" );
-    sfree( mpi_data, "main::mpi_data" );
-
 #if defined(TEST_ENERGY) || defined(TEST_FORCES)
 //    Integrate_Results(control);
 #endif
@@ -580,5 +569,24 @@ int main( int argc, char* argv[] )
     fprintf( stderr, "p%d has reached the END\n", system->my_rank );
 #endif
 
+    MPI_Finalized( &ret );
+    if ( !ret )
+    { 
+        MPI_Finalize( );
+    }
+
+    /* deallocate data structures */
+    sfree( mpi_data, "main::mpi_data" );
+    sfree( out_control, "main::out_control" );
+    for ( i = 0; i < LIST_N; ++i )
+    {
+        sfree( lists[i], "main::lists[i]" );
+    }
+    sfree( lists, "main::lists" );
+    sfree( workspace, "main::workspace" );
+    sfree( data, "main::data" );
+    sfree( control, "main::control" );
+    sfree( system, "main::system" );
+
     return 0;
 }