diff --git a/PG-PuReMD/src/cuda/cuda_init_md.cu b/PG-PuReMD/src/cuda/cuda_init_md.cu
index 81b9b3e13f82f5eb0755fc4c3980738e99be324b..5e91eaee32fc15ce63731fb7131a0c800c1ea4fe 100644
--- a/PG-PuReMD/src/cuda/cuda_init_md.cu
+++ b/PG-PuReMD/src/cuda/cuda_init_md.cu
@@ -301,14 +301,7 @@ void Cuda_Initialize( reax_system *system, control_params *control,
     Cuda_Init_ScratchArea( );
 
     /* MPI_DATATYPES */
-    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 );
 
     /* SYSTEM */
     if ( Cuda_Init_System( system, control, data, workspace, mpi_data, msg ) == FAILURE )
diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index 7d7af3e972125b4bf93be4de28365ce04dc969eb..0753ea1090ec9b2f4c60007047c86f7c4a91728e 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -814,7 +814,7 @@ void ReAllocate( reax_system *system, control_params *control,
     /* far neighbors */
     if ( renbr )
     {
-        far_nbrs = *lists + FAR_NBRS;
+        far_nbrs = lists[FAR_NBRS];
 
         if ( Nflag || realloc->num_far >= far_nbrs->num_intrs * DANGER_ZONE )
         {
@@ -878,7 +878,7 @@ void ReAllocate( reax_system *system, control_params *control,
 
         if ( Hflag || realloc->hbonds )
         {
-            ret = Reallocate_HBonds_List( system, (*lists) + HBONDS, comm );
+            ret = Reallocate_HBonds_List( system, lists[HBONDS], comm );
             realloc->hbonds = 0;
 #if defined(DEBUG_FOCUS)
             fprintf(stderr, "p%d: reallocating hbonds: total_hbonds=%d space=%dMB\n",
@@ -891,7 +891,7 @@ void ReAllocate( reax_system *system, control_params *control,
     num_bonds = est_3body = -1;
     if ( Nflag || realloc->bonds )
     {
-        Reallocate_Bonds_List( system, (*lists) + BONDS, &num_bonds,
+        Reallocate_Bonds_List( system, lists[BONDS], &num_bonds,
                                &est_3body, comm );
         realloc->bonds = 0;
         realloc->num_3body = MAX( realloc->num_3body, est_3body );
@@ -911,15 +911,15 @@ void ReAllocate( reax_system *system, control_params *control,
                  (int)(realloc->num_3body * sizeof(three_body_interaction_data) /
                        (1024 * 1024)) );
 #endif
-        Delete_List( (*lists) + THREE_BODIES, comm );
+        Delete_List( lists[THREE_BODIES], comm );
 
         if ( num_bonds == -1 )
-            num_bonds = ((*lists) + BONDS)->num_intrs;
+            num_bonds = lists[BONDS]->num_intrs;
 
         realloc->num_3body = (int)(MAX(realloc->num_3body * SAFE_ZONE, MIN_3BODIES));
 
         if ( !Make_List( num_bonds, realloc->num_3body, TYP_THREE_BODY,
-                         (*lists) + THREE_BODIES, comm ) )
+                         lists[THREE_BODIES], comm ) )
         {
             fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
             MPI_Abort( comm, CANNOT_INITIALIZE );
diff --git a/PuReMD/src/analyze.c b/PuReMD/src/analyze.c
index 228e7a6f25cc44f3fab360e653c09e1c1166d56b..5307d97feb32a6574b64de424008bba4822d2a6c 100644
--- a/PuReMD/src/analyze.c
+++ b/PuReMD/src/analyze.c
@@ -89,8 +89,8 @@ void Analyze_Fragments( reax_system *system, control_params *control,
     char fragments[MAX_FRAGMENT_TYPES][MAX_ATOM_TYPES];
     int  fragment_count[MAX_FRAGMENT_TYPES];
     molecule m;
-    reax_list *new_bonds = (*lists) + BONDS;
-    //list *old_bonds = (*lists) + OLD_BONDS;
+    reax_list *new_bonds = lists[BONDS];
+    //list *old_bonds = lists[OLD_BONDS];
 
     /* fragment analysis */
     fprintf( fout, "step%d fragments\n", data->step );
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index 71dd1f44ecf387efcba9313275187055c655c8b5..ca54db6f77f6701bdc05b851920c9f9cde5f4471 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -35,8 +35,8 @@
 void Get_dBO( reax_system *system, reax_list **lists,
               int i, int pj, real C, rvec *v )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     int start_pj, end_pj, k;
 
     pj = bonds->bond_list[pj].dbond_index;
@@ -52,8 +52,8 @@ void Get_dBO( reax_system *system, reax_list **lists,
 void Get_dBOpinpi2( reax_system *system, reax_list **lists,
                     int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -73,8 +73,8 @@ void Get_dBOpinpi2( reax_system *system, reax_list **lists,
 void Add_dBO( reax_system *system, reax_list **lists,
               int i, int pj, real C, rvec *v )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     int start_pj, end_pj, k;
 
     pj = bonds->bond_list[pj].dbond_index;
@@ -92,8 +92,8 @@ void Add_dBO( reax_system *system, reax_list **lists,
 void Add_dBOpinpi2( reax_system *system, reax_list **lists,
                     int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -113,8 +113,8 @@ void Add_dBOpinpi2( reax_system *system, reax_list **lists,
 void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
                         int i, int pj, real C )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     int start_pj, end_pj, k;
 
     pj = bonds->bond_list[pj].dbond_index;
@@ -130,8 +130,8 @@ void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
 void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
                               int i, int pj, real Cpi, real Cpi2 )
 {
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -150,7 +150,7 @@ void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
 
 void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v )
 {
-    reax_list *dDeltas = &((*lists)[DDELTAS]);
+    reax_list *dDeltas = lists[DDELTAS];
     int start = Start_Index(i, dDeltas);
     int end = End_Index(i, dDeltas);
     int k;
@@ -164,7 +164,7 @@ void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v
 void Add_dDelta_to_Forces( reax_system *system, reax_list **lists,
                            int i, real C )
 {
-    reax_list *dDeltas = &((*lists)[DDELTAS]);
+    reax_list *dDeltas = lists[DDELTAS];
     int start = Start_Index(i, dDeltas);
     int end = End_Index(i, dDeltas);
     int k;
@@ -194,8 +194,8 @@ void Calculate_dBO( int i, int pj,
     //  memset(due_j_pi, 0, sizeof(rvec)*1000 );
     //  memset(due_i_pi, 0, sizeof(rvec)*1000 );
 
-    bonds = *lists + BONDS;
-    dBOs  = *lists + DBOS;
+    bonds = lists[BONDS];
+    dBOs = lists[DBOS];
 
     start_i = Start_Index(i, bonds);
     end_i = End_Index(i, bonds);
@@ -395,7 +395,7 @@ void Calculate_dBO( int i, int pj,
 void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
                               storage *workspace, reax_list **lists )
 {
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
     bond_data *nbr_j, *nbr_k;
     bond_order_data *bo_ij, *bo_ji;
     dbond_coefficients coef;
@@ -551,7 +551,7 @@ void Add_dBond_to_Forces_NPT( int i, int pj, simulation_data *data,
 void Add_dBond_to_Forces( int i, int pj,
                           storage *workspace, reax_list **lists )
 {
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
     bond_data *nbr_j, *nbr_k;
     bond_order_data *bo_ij, *bo_ji;
     dbond_coefficients coef;
@@ -822,7 +822,7 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
     single_body_parameters *sbp_i, *sbp_j;
     two_body_parameters *twbp;
     bond_order_data *bo_ij, *bo_ji;
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
 #ifdef TEST_FORCES
     int  k, pk, start_j, end_j;
     int  top_dbo, top_dDelta;
@@ -831,8 +831,8 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
     reax_list *dDeltas, *dBOs;
     top_dbo = 0;
     top_dDelta = 0;
-    dDeltas = (*lists) + DDELTAS;
-    dBOs = (*lists) + DBOS;
+    dDeltas = lists[DDELTAS];
+    dBOs = lists[DBOS];
 
     //for( i = 0; i < system->N; ++i )
     //  qsort( &(bonds->bond_list[Start_Index(i, bonds)]),
diff --git a/PuReMD/src/bonds.c b/PuReMD/src/bonds.c
index 5b0b31f6df68b824837b252fc37534504d482e21..dec58cdb6daa0be7f01d33523db428ad7e39cd50 100644
--- a/PuReMD/src/bonds.c
+++ b/PuReMD/src/bonds.c
@@ -51,7 +51,7 @@ void Bonds( reax_system *system, control_params *control,
     bond_order_data *bo_ij;
     reax_list *bonds;
 
-    bonds = (*lists) + BONDS;
+    bonds = lists[BONDS];
     gp3 = system->reax_param.gp.l[3];
     gp4 = system->reax_param.gp.l[4];
     gp7 = system->reax_param.gp.l[7];
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 93ac4dee2b14297b7246f15149454d1a0c770411..8419e3efe44377e7d1680fe09b2bc06ef5756ca6 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -245,7 +245,7 @@ void Pack_Boundary_Atom( boundary_atom *matm, reax_atom *ratm, int i )
 {
     matm->orig_id  = ratm->orig_id;
     matm->imprt_id = i;
-    matm->type     = ratm->type;
+    matm->type = ratm->type;
     matm->num_bonds = ratm->num_bonds;
     matm->num_hbonds = ratm->num_hbonds;
     rvec_Copy( matm->x, ratm->x );
@@ -278,12 +278,15 @@ void Sort_Boundary_Atoms( reax_system *system, int start, int end,
     fprintf( stderr, "p%d sort_exchange: start=%d end=%d dim=%d starting...\n",
              system->my_rank, start, end, dim );
 #endif
+
     atoms = system->my_atoms;
     my_box = &( system->my_box );
 
     /* place each atom into the appropriate outgoing list */
     for ( i = start; i < end; ++i )
+    {
         for ( d = dim; d < 3; ++d )
+        {
             for ( p = 2 * d; p < 2 * d + 2; ++p )
             {
                 nbr_pr = &( system->my_nbrs[p] );
@@ -296,12 +299,16 @@ void Sort_Boundary_Atoms( reax_system *system, int start, int end,
                     Pack_Boundary_Atom( out_buf + out_cnt, atoms + i, i );
                 }
             }
+        }
+    }
 
 #if defined(DEBUG_FOCUS)
     for ( i = 2 * dim; i < 2 * dim + 2; ++i )
+    {
         fprintf( stderr, "p%d to p%d(nbr%d) # of exchanges to send = %d\n",
                  system->my_rank, system->my_nbrs[i].rank, i,
                  out_bufs[i].cnt );
+    }
     fprintf( stderr, "p%d sort_exchange: start=%d end=%d dim=%d done!\n",
              system->my_rank, start, end, dim );
 #endif
@@ -356,6 +363,7 @@ void Estimate_Boundary_Atoms( reax_system *system, int start, int end,
 
     /* sort the atoms to their outgoing buffers */
     for ( i = 0; i < end; ++i )
+    {
         for ( p = 2 * d; p < 2 * d + 2; ++p )
         {
             nbr_pr = &( system->my_nbrs[p] );
@@ -368,6 +376,7 @@ void Estimate_Boundary_Atoms( reax_system *system, int start, int end,
                 Pack_Boundary_Atom( out_buf + out_cnt, atoms + i, i );
             }
         }
+    }
 
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d estimate_exchange: end=%d dim=%d done!\n",
@@ -525,7 +534,11 @@ int SendRecv( reax_system* system, mpi_datatypes *mpi_data, MPI_Datatype type,
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d sendrecv: entered\n", system->my_rank );
 #endif
-    if ( clr ) Reset_Out_Buffers( mpi_data->out_buffers, system->num_nbrs );
+
+    if ( clr )
+    {
+        Reset_Out_Buffers( mpi_data->out_buffers, system->num_nbrs );
+    }
     comm = mpi_data->comm_mesh3D;
     in1 = mpi_data->in1_buffer;
     in2 = mpi_data->in2_buffer;
@@ -545,8 +558,10 @@ int SendRecv( reax_system* system, mpi_datatypes *mpi_data, MPI_Datatype type,
 
         /* for estimates in1_buffer & in2_buffer will be NULL */
         if ( est_flag )
+        {
             Estimate_Init_Storage( system->my_rank, nbr1, nbr2, d,
-                                   &max, nrecv, &in1, &in2, comm );
+                    &max, nrecv, &in1, &in2, comm );
+        }
 
         /* initiate recvs */
         MPI_Irecv( in1, nrecv[2 * d], type, nbr1->rank, 2 * d + 1, comm, &req1 );
@@ -554,9 +569,9 @@ int SendRecv( reax_system* system, mpi_datatypes *mpi_data, MPI_Datatype type,
 
         /* send both messages in dimension d */
         MPI_Send( out_bufs[2 * d].out_atoms, out_bufs[2 * d].cnt, type,
-                  nbr1->rank, 2 * d, comm );
+                nbr1->rank, 2 * d, comm );
         MPI_Send( out_bufs[2 * d + 1].out_atoms, out_bufs[2 * d + 1].cnt, type,
-                  nbr2->rank, 2 * d + 1, comm );
+                nbr2->rank, 2 * d + 1, comm );
 
         /* recv and unpack atoms from nbr1 in dimension d */
         MPI_Wait( &req1, &stat1 );
@@ -593,42 +608,62 @@ void Comm_Atoms( reax_system *system, control_params *control,
     real t_start = 0, t_elapsed = 0;
 
     if ( system->my_rank == MASTER_NODE )
+    {
         t_start = Get_Time( );
+    }
 #endif
 
     if ( renbr )
     {
-        for ( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = system->est_trans;
+        for ( i = 0; i < MAX_NBRS; ++i )
+        {
+            nrecv[i] = system->est_trans;
+        }
         system->n = SendRecv( system, mpi_data, mpi_data->mpi_atom_type, nrecv,
-                              Sort_Transfer_Atoms, Unpack_Transfer_Message, 1 );
+                Sort_Transfer_Atoms, Unpack_Transfer_Message, 1 );
         Bin_My_Atoms( system, &(workspace->realloc) );
         Reorder_My_Atoms( system, workspace );
+
 #if defined(DEBUG_FOCUS)
         fprintf( stderr, "p%d updated local atoms, n=%d\n",
                  system->my_rank, system->n );
         MPI_Barrier( mpi_data->world );
 #endif
 
-        for ( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = system->my_nbrs[i].est_recv;
-        system->N = SendRecv(system, mpi_data, mpi_data->boundary_atom_type, nrecv,
-                             Sort_Boundary_Atoms, Unpack_Exchange_Message, 1);
+        for ( i = 0; i < MAX_NBRS; ++i )
+        {
+            nrecv[i] = system->my_nbrs[i].est_recv;
+        }
+        system->N = SendRecv( system, mpi_data, mpi_data->boundary_atom_type, nrecv,
+                Sort_Boundary_Atoms, Unpack_Exchange_Message, 1 );
+
 #if defined(DEBUG_FOCUS)
         fprintf( stderr, "p%d: exchanged boundary atoms, N=%d\n",
-                 system->my_rank, system->N );
+                system->my_rank, system->N );
+
         for ( i = 0; i < MAX_NBRS; ++i )
+        {
             fprintf( stderr, "p%d: nbr%d(p%d) str=%d cnt=%d end=%d\n",
                      system->my_rank, i, system->my_nbrs[i].rank,
                      system->my_nbrs[i].atoms_str,  system->my_nbrs[i].atoms_cnt,
                      system->my_nbrs[i].atoms_str + system->my_nbrs[i].atoms_cnt );
+        }
+
         MPI_Barrier( mpi_data->world );
 #endif
+
         Bin_Boundary_Atoms( system );
     }
     else
     {
-        for ( i = 0; i < MAX_NBRS; ++i ) nrecv[i] = system->my_nbrs[i].atoms_cnt;
+        for ( i = 0; i < MAX_NBRS; ++i )
+        {
+            nrecv[i] = system->my_nbrs[i].atoms_cnt;
+        }
+
         SendRecv( system, mpi_data, mpi_data->mpi_rvec, nrecv,
-                  Sort_Position_Updates, Unpack_Position_Updates, 0 );
+                Sort_Position_Updates, Unpack_Position_Updates, 0 );
+
 #if defined(DEBUG_FOCUS)
         fprintf( stderr, "p%d: updated positions\n", system->my_rank );
         MPI_Barrier( mpi_data->world );
@@ -642,6 +677,7 @@ void Comm_Atoms( reax_system *system, control_params *control,
         data->timing.comm += t_elapsed;
     }
 #endif
+
 #if defined(DEBUG_FOCUS)
     fprintf(stderr, "p%d @ renbr=%d: comm_atoms done\n", system->my_rank, renbr);
     fprintf( stderr, "p%d: system->n = %d, system->N = %d\n",
diff --git a/PuReMD/src/control.c b/PuReMD/src/control.c
index 07b46d3a4d096f10d390ad203d13b14d151e36f9..cd95c2465dc9a17edcf951295babb9527506b0c4 100644
--- a/PuReMD/src/control.c
+++ b/PuReMD/src/control.c
@@ -114,9 +114,8 @@ char Read_Control_File( char *control_file, control_params* control,
         tmp[i] = (char*) malloc(sizeof(char) * MAX_LINE);
 
     /* read control parameters file */
-    while (!feof(fp))
+    while ( fgets( s, MAX_LINE, fp ) )
     {
-        fgets( s, MAX_LINE, fp );
         c = Tokenize( s, &tmp );
         //fprintf( stderr, "%s\n", s );
 
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 888c6acae9db372721326582558605acc73c464e..d438a0eba65d50ae11a8f6bd03fa5fe79492474b 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -142,7 +142,7 @@ void Compute_Total_Force( reax_system *system, control_params *control,
                           reax_list **lists, mpi_datatypes *mpi_data )
 {
     int i, pj;
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
 
     for ( i = 0; i < system->N; ++i )
         for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
@@ -194,7 +194,7 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
     /* bond list */
     if ( N > 0 )
     {
-        bonds = *lists + BONDS;
+        bonds = lists[BONDS];
 
         for ( i = 0; i < N; ++i )
         {
@@ -221,7 +221,7 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
     /* hbonds list */
     if ( numH > 0 )
     {
-        hbonds = *lists + HBONDS;
+        hbonds = lists[HBONDS];
 
         for ( i = 0; i < n; ++i )
         {
@@ -337,9 +337,9 @@ void Init_Forces( reax_system *system, control_params *control,
     far_neighbor_data *nbr_pj;
     reax_atom *atom_i, *atom_j;
 
-    far_nbrs = *lists + FAR_NBRS;
-    bonds = *lists + BONDS;
-    hbonds = *lists + HBONDS;
+    far_nbrs = lists[FAR_NBRS];
+    bonds = lists[BONDS];
+    hbonds = lists[HBONDS];
 
     for ( i = 0; i < system->n; ++i )
         workspace->bond_mark[i] = 0;
@@ -605,9 +605,9 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
     far_neighbor_data *nbr_pj;
     reax_atom *atom_i, *atom_j;
 
-    far_nbrs = *lists + FAR_NBRS;
-    bonds = *lists + BONDS;
-    hbonds = *lists + HBONDS;
+    far_nbrs = lists[FAR_NBRS];
+    bonds = lists[BONDS];
+    hbonds = lists[HBONDS];
 
     for ( i = 0; i < system->n; ++i )
         workspace->bond_mark[i] = 0;
@@ -791,7 +791,7 @@ void Estimate_Storages( reax_system *system, control_params *control,
     far_neighbor_data *nbr_pj;
     reax_atom *atom_i, *atom_j;
 
-    far_nbrs = *lists + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     *Htop = 0;
     memset( hb_top, 0, sizeof(int) * system->local_cap );
     memset( bond_top, 0, sizeof(int) * system->total_cap );
diff --git a/PuReMD/src/geo_tools.c b/PuReMD/src/geo_tools.c
index dd5115c07649d3c7571ae36fc79cd221eb9ad0f1..c1e3549fedf2039f96af84aa263cd114c7dee3cb 100644
--- a/PuReMD/src/geo_tools.c
+++ b/PuReMD/src/geo_tools.c
@@ -212,10 +212,8 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
     system->bigN = system->n = 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 )
         {
diff --git a/PuReMD/src/hydrogen_bonds.c b/PuReMD/src/hydrogen_bonds.c
index 6aea423234e68ba53729b8841d6cae15155bbb9e..cba267fffa71ef7a11ac23def3c72cae92838653 100644
--- a/PuReMD/src/hydrogen_bonds.c
+++ b/PuReMD/src/hydrogen_bonds.c
@@ -59,9 +59,9 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
     bond_data *bond_list;
     hbond_data *hbond_list;
 
-    bonds = (*lists) + BONDS;
+    bonds = lists[BONDS];
     bond_list = bonds->bond_list;
-    hbonds = (*lists) + HBONDS;
+    hbonds = lists[HBONDS];
     hbond_list = hbonds->hbond_list;
 
     /* loops below discover the Hydrogen bonds between i-j-k triplets.
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index b35a46887e638424b5d570841cd12e63c00462f9..caad5420453b4e9329cc0b4ee9382462c0d9e862 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -553,8 +553,8 @@ int  Init_Lists( reax_system *system, control_params *control,
              system->my_rank, system->local_cap, system->total_cap );
 #endif
 
-    if (!Make_List( system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
-                    *lists + FAR_NBRS, comm ))
+    if ( !Make_List( system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
+                lists[FAR_NBRS], comm ) )
     {
         fprintf(stderr, "Problem in initializing far nbrs list. Terminating!\n");
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -605,7 +605,7 @@ int  Init_Lists( reax_system *system, control_params *control,
         total_hbonds = MAX( total_hbonds * SAFER_ZONE, MIN_CAP * MIN_HBONDS );
 
         if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
-                         *lists + HBONDS, comm ) )
+                         lists[HBONDS], comm ) )
         {
             fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
             MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -618,7 +618,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     }
 
     /* bonds list */
-    //Allocate_Bond_List( system->N, bond_top, (*lists)+BONDS );
+    //Allocate_Bond_List( system->N, bond_top, lists[BONDS] );
     //num_bonds = bond_top[system->N-1];
     total_bonds = 0;
     for ( i = 0; i < system->N; ++i )
@@ -629,7 +629,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     bond_cap = MAX( total_bonds * SAFE_ZONE, MIN_CAP * MIN_BONDS );
 
     if ( !Make_List( system->total_cap, bond_cap, TYP_BOND,
-                     *lists + BONDS, comm ) )
+                     lists[BONDS], comm ) )
     {
         fprintf( stderr, "not enough space for bonds list. terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -643,7 +643,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     /* 3bodies list */
     cap_3body = MAX( num_3body * SAFE_ZONE, MIN_3BODIES );
     if ( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
-                     *lists + THREE_BODIES, comm ) )
+                     lists[THREE_BODIES], comm ) )
     {
         fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -656,7 +656,7 @@ int  Init_Lists( reax_system *system, control_params *control,
 
 #if defined(TEST_FORCES)
     if ( !Make_List( system->total_cap, bond_cap * 8, TYP_DDELTA,
-                     *lists + DDELTAS, comm ) )
+                     lists[DDELTAS], comm ) )
     {
         fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -665,7 +665,7 @@ int  Init_Lists( reax_system *system, control_params *control,
              system->my_rank, bond_cap * 30,
              bond_cap * 8 * sizeof(dDelta_data) / (1024 * 1024) );
 
-    if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, *lists + DBOS, comm ) )
+    if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, lists[DBOS], comm ) )
     {
         fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -709,7 +709,7 @@ int  Init_Lists( reax_system *system, control_params *control,
         total_hbonds = (int)(MAX( total_hbonds * SAFER_ZONE, MIN_CAP * MIN_HBONDS ));
 
         if ( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
-                         *lists + HBONDS, comm ) )
+                         lists[HBONDS], comm ) )
         {
             fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
             MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -722,7 +722,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     }
 
     /* bonds list */
-    //Allocate_Bond_List( system->N, bond_top, (*lists)+BONDS );
+    //Allocate_Bond_List( system->N, bond_top, lists[BONDS] );
     //num_bonds = bond_top[system->N-1];
     total_bonds = 0;
     for ( i = 0; i < system->N; ++i )
@@ -733,7 +733,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     bond_cap = (int)(MAX( total_bonds * SAFE_ZONE, MIN_CAP * MIN_BONDS ));
 
     if ( !Make_List( system->total_cap, bond_cap, TYP_BOND,
-                     *lists + BONDS, comm ) )
+                     lists[BONDS], comm ) )
     {
         fprintf( stderr, "not enough space for bonds list. terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -747,7 +747,7 @@ int  Init_Lists( reax_system *system, control_params *control,
     /* 3bodies list */
     cap_3body = (int)(MAX( num_3body * SAFE_ZONE, MIN_3BODIES ));
     if ( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
-                     *lists + THREE_BODIES, comm ) )
+                     lists[THREE_BODIES], comm ) )
     {
         fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -760,7 +760,7 @@ int  Init_Lists( reax_system *system, control_params *control,
 
 #if defined(TEST_FORCES)
     if ( !Make_List( system->total_cap, bond_cap * 8, TYP_DDELTA,
-                     *lists + DDELTAS, comm ) )
+                     lists[DDELTAS], comm ) )
     {
         fprintf( stderr, "Problem in initializing dDelta list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
@@ -769,7 +769,7 @@ int  Init_Lists( reax_system *system, control_params *control,
              system->my_rank, bond_cap * 30,
              bond_cap * 8 * sizeof(dDelta_data) / (1024 * 1024) );
 
-    if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, (*lists) + DBOS, comm ) )
+    if ( !Make_List( bond_cap, bond_cap * 50, TYP_DBO, lists[DBOS], comm ) )
     {
         fprintf( stderr, "Problem in initializing dBO list. Terminating!\n" );
         MPI_Abort( comm, INSUFFICIENT_MEMORY );
diff --git a/PuReMD/src/io_tools.c b/PuReMD/src/io_tools.c
index 5a5e6f2a08026ba5c9b6f785aa7f607360300c07..6c0cc5c8e755ff0b23fdf8234e7530a3abceed5d 100644
--- a/PuReMD/src/io_tools.c
+++ b/PuReMD/src/io_tools.c
@@ -844,7 +844,7 @@ void Print_Far_Neighbors( reax_system *system, reax_list **lists,
 
     sprintf( fname, "%s.far_nbrs.%d", control->sim_name, system->my_rank );
     fout      = fopen( fname, "w" );
-    far_nbrs = (*lists) + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     natoms = system->N;
 
     for ( i = 0; i < natoms; ++i )
@@ -1258,8 +1258,8 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
     int i, pj, pk;
     bond_order_data *bo_ij;
     dbond_data *dbo_k;
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *dBOs  = (*lists) + DBOS;
+    reax_list *bonds = lists[BONDS];
+    reax_list *dBOs = lists[DBOS];
 
     /* bond orders */
     fprintf( out_control->fbo, "step: %d\n", data->step );
@@ -1421,7 +1421,7 @@ void Print_Far_Neighbors_List( reax_system *system, reax_list **lists,
     int temp[500];
     reax_list *far_nbrs;
 
-    far_nbrs = (*lists) + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     fprintf( out_control->flist, "step: %d\n", data->step );
     fprintf( out_control->flist, "%6s\t%-38s\n", "atom", "Far_nbrs_list");
 
@@ -1452,7 +1452,7 @@ void Print_Bond_List( reax_system *system, control_params *control,
                       output_controls *out_control)
 {
     int i, j, id_i, id_j, nbr, pj;
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
 
     int temp[500];
     int num = 0;
diff --git a/PuReMD/src/list.c b/PuReMD/src/list.c
index 19f233db605687076fea2e03786fc0b597b6a373..4ccb03ed94174c6a5bb184ae48687a6e40c73515 100644
--- a/PuReMD/src/list.c
+++ b/PuReMD/src/list.c
@@ -41,11 +41,12 @@ int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
     l->end_index = (int*) smalloc( n * sizeof(int), "list:end_index", comm );
 
     l->type = type;
+
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "list: n=%d num_intrs=%d type=%d\n", l->n, l->num_intrs, l->type );
 #endif
 
-    switch (l->type)
+    switch ( l->type )
     {
     case TYP_VOID:
         l->v = (void*) smalloc(l->num_intrs * sizeof(void*), "list:v", comm);
@@ -129,31 +130,3 @@ void Delete_List( reax_list *l, MPI_Comm comm )
         MPI_Abort( comm, INVALID_INPUT );
     }
 }
-
-
-#if defined(PURE_REAX)
-inline int Num_Entries( int i, reax_list *l )
-{
-    return l->end_index[i] - l->index[i];
-}
-
-inline int Start_Index( int i, reax_list *l )
-{
-    return l->index[i];
-}
-
-inline int End_Index( int i, reax_list *l )
-{
-    return l->end_index[i];
-}
-
-inline void Set_Start_Index( int i, int val, reax_list *l )
-{
-    l->index[i] = val;
-}
-
-inline void Set_End_Index( int i, int val, reax_list *l )
-{
-    l->end_index[i] = val;
-}
-#endif
diff --git a/PuReMD/src/list.h b/PuReMD/src/list.h
index 5cc3a62630fdb29f3ac07ffb8b30469de49f34ea..da400b76e1303bf5e587e949f664bc93225de3dc 100644
--- a/PuReMD/src/list.h
+++ b/PuReMD/src/list.h
@@ -24,40 +24,40 @@
 
 #include "reax_types.h"
 
-int  Make_List( int, int, int, reax_list*, MPI_Comm );
+
+int Make_List( int, int, int, reax_list*, MPI_Comm );
+
 void Delete_List( reax_list*, MPI_Comm );
 
-inline int  Num_Entries(int, reax_list*);
-inline int  Start_Index( int, reax_list* );
-inline int  End_Index( int, reax_list* );
-inline void Set_Start_Index(int, int, reax_list*);
-inline void Set_End_Index(int, int, reax_list*);
 
-#if defined(LAMMPS_REAX)
-inline int Num_Entries( int i, reax_list *l )
+static inline int Num_Entries( int i, reax_list *l )
 {
     return l->end_index[i] - l->index[i];
 }
 
-inline int Start_Index( int i, reax_list *l )
+
+static inline int Start_Index( int i, reax_list *l )
 {
     return l->index[i];
 }
 
-inline int End_Index( int i, reax_list *l )
+
+static inline int End_Index( int i, reax_list *l )
 {
     return l->end_index[i];
 }
 
-inline void Set_Start_Index( int i, int val, reax_list *l )
+
+static inline void Set_Start_Index( int i, int val, reax_list *l )
 {
     l->index[i] = val;
 }
 
-inline void Set_End_Index( int i, int val, reax_list *l )
+
+static inline void Set_End_Index( int i, int val, reax_list *l )
 {
     l->end_index[i] = val;
 }
-#endif // LAMMPS_REAX
+
 
 #endif
diff --git a/PuReMD/src/multi_body.c b/PuReMD/src/multi_body.c
index 1c80cd86f0b24b9af718d6c8c877a347cca6e98c..7ef1f12d1103c5d3ef73fe15487329627b26caee 100644
--- a/PuReMD/src/multi_body.c
+++ b/PuReMD/src/multi_body.c
@@ -53,7 +53,7 @@ void Atom_Energy( reax_system *system, control_params *control,
     two_body_parameters *twbp;
     bond_data *pbond;
     bond_order_data *bo_ij;
-    reax_list *bonds = (*lists) + BONDS;
+    reax_list *bonds = lists[BONDS];
 
     /* Initialize parameters */
     p_lp1 = system->reax_param.gp.l[15];
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index 8f96541e2dd347a4565622f542a7b784352769a8..5be0016eda64a6e4ec56f02d77695cd212fcdd55 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -86,7 +86,7 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
 
     // fprintf( stderr, "\n\tentered nbrs - " );
     g = &( system->my_grid );
-    far_nbrs = (*lists) + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     num_far = 0;
 
     /* first pick up a cell in the grid */
diff --git a/PuReMD/src/nonbonded.c b/PuReMD/src/nonbonded.c
index 1bdd2ab6e6a5bd2b932a126aedb3b2ca3abc9814..ab25b807d5e0a0263a4d8f89ecc475c8615b8bdd 100644
--- a/PuReMD/src/nonbonded.c
+++ b/PuReMD/src/nonbonded.c
@@ -52,7 +52,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
     // rtensor temp_rtensor, total_rtensor;
 
     natoms = system->n;
-    far_nbrs = (*lists) + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     p_vdW1 = system->reax_param.gp.l[28];
     p_vdW1i = 1.0 / p_vdW1;
     e_core = 0;
@@ -230,7 +230,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
     LR_lookup_table *t;
 
     natoms = system->n;
-    far_nbrs = (*lists) + FAR_NBRS;
+    far_nbrs = lists[FAR_NBRS];
     steps = data->step - data->prev_steps;
     update_freq = out_control->energy_update_freq;
     update_energies = update_freq > 0 && steps % update_freq == 0;
diff --git a/PuReMD/src/parallelreax.c b/PuReMD/src/parallelreax.c
index 6d0955c70c4dfd5813a2a29ae3be191996270dea..4b401ad603db53f5b982a2c7fc24861a9b9a6495 100644
--- a/PuReMD/src/parallelreax.c
+++ b/PuReMD/src/parallelreax.c
@@ -125,46 +125,60 @@ int main( int argc, char* argv[] )
     int i;
     real t_start = 0, t_elapsed;
 
+    MPI_Init( &argc, &argv );
+
     if ( argc != 4 )
     {
-        usage(argv);
-        exit( INVALID_INPUT );
+        usage( argv );
+        MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
     }
 
     /* allocated main datastructures */
     system = (reax_system *)
-             smalloc( sizeof(reax_system), "system", MPI_COMM_WORLD );
+        smalloc( sizeof(reax_system), "system", MPI_COMM_WORLD );
     control = (control_params *)
-              smalloc( sizeof(control_params), "control", MPI_COMM_WORLD );
+        smalloc( sizeof(control_params), "control", MPI_COMM_WORLD );
     data = (simulation_data *)
-           smalloc( sizeof(simulation_data), "data", MPI_COMM_WORLD );
+        smalloc( sizeof(simulation_data), "data", MPI_COMM_WORLD );
 
     workspace = (storage *)
-                smalloc( sizeof(storage), "workspace", MPI_COMM_WORLD );
+        smalloc( sizeof(storage), "workspace", MPI_COMM_WORLD );
     lists = (reax_list **)
-            smalloc( LIST_N * sizeof(reax_list*), "lists", MPI_COMM_WORLD );
+        smalloc( LIST_N * sizeof(reax_list*), "lists", MPI_COMM_WORLD );
     for ( i = 0; i < LIST_N; ++i )
     {
         lists[i] = (reax_list *)
-                   smalloc( sizeof(reax_list), "lists[i]", MPI_COMM_WORLD );
+            smalloc( sizeof(reax_list), "lists[i]", MPI_COMM_WORLD );
         lists[i]->allocated = 0;
+        lists[i]->n = 0;
+        lists[i]->num_intrs = 0;
+        lists[i]->index = NULL;
+        lists[i]->end_index = NULL;
+        lists[i]->type = 0;
+        lists[i]->v = NULL;
+        lists[i]->three_body_list = NULL;
+        lists[i]->bond_list = NULL;
+        lists[i]->dbo_list = NULL;
+        lists[i]->dDelta_list = NULL;
+        lists[i]->far_nbr_list = NULL;
+        lists[i]->hbond_list = NULL;
     }
     out_control = (output_controls *)
-                  smalloc( sizeof(output_controls), "out_control", MPI_COMM_WORLD );
+        smalloc( sizeof(output_controls), "out_control", MPI_COMM_WORLD );
     mpi_data = (mpi_datatypes *)
-               smalloc( sizeof(mpi_datatypes), "mpi_data", MPI_COMM_WORLD );
+        smalloc( sizeof(mpi_datatypes), "mpi_data", MPI_COMM_WORLD );
 
     /* setup the parallel environment */
-    MPI_Init( &argc, &argv );
     MPI_Comm_size( MPI_COMM_WORLD, &(control->nprocs) );
     MPI_Comm_rank( MPI_COMM_WORLD, &(system->my_rank) );
     system->wsize = control->nprocs;
     system->global_offset = (int*)
-                            scalloc( system->wsize + 1, sizeof(int), "global_offset", MPI_COMM_WORLD );
+        scalloc( system->wsize + 1, sizeof(int), "global_offset", MPI_COMM_WORLD );
 
     /* read system description files */
     Read_System( argv[1], argv[2], argv[3], system, control,
-                 data, workspace, out_control, mpi_data );
+            data, workspace, out_control, mpi_data );
+
 #if defined(DEBUG)
     fprintf( stderr, "p%d: read simulation info\n", system->my_rank );
     MPI_Barrier( MPI_COMM_WORLD );
@@ -176,6 +190,7 @@ int main( int argc, char* argv[] )
 
     /* initialize datastructures */
     Initialize( system, control, data, workspace, lists, out_control, mpi_data );
+
 #if defined(DEBUG)
     fprintf( stderr, "p%d: initializated data structures\n", system->my_rank );
     MPI_Barrier( mpi_data->world );
@@ -186,9 +201,10 @@ int main( int argc, char* argv[] )
     Reset( system, control, data, workspace, lists, MPI_COMM_WORLD );
     Generate_Neighbor_Lists( system, data, workspace, lists );
     Compute_Forces( system, control, data, workspace,
-                    lists, out_control, mpi_data );
+            lists, out_control, mpi_data );
     Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
     Output_Results( system, control, data, lists, out_control, mpi_data );
+
 #if defined(DEBUG)
     fprintf( stderr, "p%d: computed forces at t0\n", system->my_rank );
     MPI_Barrier( mpi_data->world );
@@ -214,6 +230,7 @@ int main( int argc, char* argv[] )
             else if ( out_control->restart_format == WRITE_BINARY )
                 Write_Binary_Restart( system, control, data, out_control, mpi_data );
         }
+
 #if defined(DEBUG)
         fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
         MPI_Barrier( mpi_data->world );
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 073918c4a8309130207e3fdce05a7d7a0be0f0ca..3ebcd73cb15c89910d87df80714eb73dcec5604b 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -22,15 +22,15 @@
 #ifndef __REAX_TYPES_H_
 #define __REAX_TYPES_H_
 
-#include "ctype.h"
-#include "math.h"
-#include "mpi.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "sys/time.h"
-#include "time.h"
-#include "zlib.h"
+#include <ctype.h>
+#include <math.h>
+#include <mpi.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/time.h>
+#include <time.h>
+#include <zlib.h>
 
 /************* SOME DEFS - crucial for reax_types.h *********/
 
@@ -678,8 +678,11 @@ typedef struct
 {
     int thb;
     int pthb; // pointer to the third body on the central atom's nbrlist
-    real theta, cos_theta;
-    rvec dcos_di, dcos_dj, dcos_dk;
+    real theta;
+    real cos_theta;
+    rvec dcos_di;
+    rvec dcos_dj;
+    rvec dcos_dk;
 } three_body_interaction_data;
 
 
@@ -715,12 +718,28 @@ typedef struct
 
 typedef struct
 {
-    real BO, BO_s, BO_pi, BO_pi2;
-    real Cdbo, Cdbopi, Cdbopi2;
-    real C1dbo, C2dbo, C3dbo;
-    real C1dbopi, C2dbopi, C3dbopi, C4dbopi;
-    real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
-    rvec dBOp, dln_BOp_s, dln_BOp_pi, dln_BOp_pi2;
+    real BO;
+    real BO_s;
+    real BO_pi;
+    real BO_pi2;
+    real Cdbo;
+    real Cdbopi;
+    real Cdbopi2;
+    real C1dbo;
+    real C2dbo;
+    real C3dbo;
+    real C1dbopi;
+    real C2dbopi;
+    real C3dbopi;
+    real C4dbopi;
+    real C1dbopi2;
+    real C2dbopi2;
+    real C3dbopi2;
+    real C4dbopi2;
+    rvec dBOp;
+    rvec dln_BOp_s;
+    rvec dln_BOp_pi;
+    rvec dln_BOp_pi2;
 } bond_order_data;
 
 typedef struct
@@ -729,7 +748,7 @@ typedef struct
     int sym_index;
     int dbond_index;
     ivec rel_box;
-    //  rvec ext_factor;
+//    rvec ext_factor;
     real d;
     rvec dvec;
     bond_order_data bo_data;
@@ -867,15 +886,15 @@ typedef struct
     int *end_index;
 
     int type;
-    list_type select;
+//    list_type select;
 
     void *v;
     three_body_interaction_data *three_body_list;
-    bond_data          *bond_list;
-    dbond_data         *dbo_list;
-    dDelta_data        *dDelta_list;
-    far_neighbor_data  *far_nbr_list;
-    hbond_data         *hbond_list;
+    bond_data *bond_list;
+    dbond_data *dbo_list;
+    dDelta_data *dDelta_list;
+    far_neighbor_data *far_nbr_list;
+    hbond_data *hbond_list;
 } reax_list;
 
 
diff --git a/PuReMD/src/reset_tools.c b/PuReMD/src/reset_tools.c
index af869e980d8e5b1705bde5e19d156bbeef364b3f..eb7c64652a34f4825121305d256843bb749ace26 100644
--- a/PuReMD/src/reset_tools.c
+++ b/PuReMD/src/reset_tools.c
@@ -74,7 +74,11 @@ void Reset_Energies( energy_data *en )
 
 void Reset_Temperatures( simulation_data *data )
 {
-    data->therm.T = 0;
+    data->therm.T = 0.0;
+    data->therm.xi = 0.0;
+    data->therm.v_xi = 0.0;
+    data->therm.v_xi_old = 0.0;
+    data->therm.G_xi = 0.0;
 }
 
 
@@ -181,7 +185,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
     /* bonds list */
     if ( system->N > 0 )
     {
-        bonds = (*lists) + BONDS;
+        bonds = lists[BONDS];
         total_bonds = 0;
 
         /* reset start-end indexes */
@@ -212,7 +216,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
     /* hbonds list */
     if ( control->hbond_cut > 0 && system->numH > 0 )
     {
-        hbonds = (*lists) + HBONDS;
+        hbonds = lists[HBONDS];
         total_hbonds = 0;
 
         /* reset start-end indexes */
diff --git a/PuReMD/src/system_props.c b/PuReMD/src/system_props.c
index 50e328a45f00c986a75cbbc87587363892c7e9cd..05b74323b902a625f15a1a8681c92ee7e3446f47 100644
--- a/PuReMD/src/system_props.c
+++ b/PuReMD/src/system_props.c
@@ -56,7 +56,7 @@ void Temperature_Control( control_params *control, simulation_data *data )
 
 
 void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
-                             MPI_Comm comm )
+        MPI_Comm comm )
 {
     int i;
     rvec p;
@@ -64,7 +64,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
 
     data->my_en.e_kin = 0.0;
     data->sys_en.e_kin = 0.0;
-    data->therm.T = 0;
+    data->therm.T = 0.0;
 
     for ( i = 0; i < system->n; i++ )
     {
@@ -75,13 +75,15 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
     }
 
     MPI_Allreduce( &data->my_en.e_kin,  &data->sys_en.e_kin,
-                   1, MPI_DOUBLE, MPI_SUM, comm );
+            1, MPI_DOUBLE, MPI_SUM, comm );
 
-    data->therm.T = (2. * data->sys_en.e_kin) / (data->N_f * K_B);
+    data->therm.T = (2.0 * data->sys_en.e_kin) / (data->N_f * K_B);
 
     // avoid T being an absolute zero, might cause F.P.E!
-    if ( fabs(data->therm.T) < ALMOST_ZERO )
+    if ( fabs( data->therm.T ) < ALMOST_ZERO )
+    {
         data->therm.T = ALMOST_ZERO;
+    }
 }
 
 
diff --git a/PuReMD/src/torsion_angles.c b/PuReMD/src/torsion_angles.c
index d0947d3efa669056ab06060418a0aa088399fa24..f915fdb8325de91f5c5ea72f47e36ae0ab306992 100644
--- a/PuReMD/src/torsion_angles.c
+++ b/PuReMD/src/torsion_angles.c
@@ -191,8 +191,8 @@ void Torsion_Angles( reax_system *system, control_params *control,
     real p_tor3 = system->reax_param.gp.l[24];
     real p_tor4 = system->reax_param.gp.l[25];
     real p_cot2 = system->reax_param.gp.l[27];
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *thb_intrs = (*lists) + THREE_BODIES;
+    reax_list *bonds = lists[BONDS];
+    reax_list *thb_intrs = lists[THREE_BODIES];
     // char  fname[100];
     // FILE *ftor;
 
diff --git a/PuReMD/src/traj.c b/PuReMD/src/traj.c
index e455d74b3accf3d925a7df81a2c387e0983d288f..40df8bbd1860eff00141ed97013e596430daf69c 100644
--- a/PuReMD/src/traj.c
+++ b/PuReMD/src/traj.c
@@ -1097,10 +1097,10 @@ int Append_Frame( reax_system *system, control_params *control,
         Write_Atoms( system, control, out_control, mpi_data );
 
     if ( out_control->write_bonds )
-        Write_Bonds( system, control, (*lists + BONDS), out_control, mpi_data );
+        Write_Bonds( system, control, lists[BONDS], out_control, mpi_data );
 
     if ( out_control->write_angles )
-        Write_Angles( system, control, (*lists + BONDS), (*lists + THREE_BODIES),
+        Write_Angles( system, control, lists[BONDS], lists[THREE_BODIES],
                       out_control, mpi_data );
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d: appended frame %d\n", system->my_rank, data->step );
diff --git a/PuReMD/src/valence_angles.c b/PuReMD/src/valence_angles.c
index db823abd7a65444d1c2951d93701437b54aa39e7..3e90162b4bb95ef31568626ba44707e3eece5c56 100644
--- a/PuReMD/src/valence_angles.c
+++ b/PuReMD/src/valence_angles.c
@@ -107,8 +107,8 @@ void Valence_Angles( reax_system *system, control_params *control,
     three_body_interaction_data *p_ijk, *p_kji;
     bond_data *pbond_ij, *pbond_jk, *pbond_jt;
     bond_order_data *bo_ij, *bo_jk, *bo_jt;
-    reax_list *bonds = (*lists) + BONDS;
-    reax_list *thb_intrs =  (*lists) + THREE_BODIES;
+    reax_list *bonds = lists[BONDS];
+    reax_list *thb_intrs =  lists[THREE_BODIES];
 
 
     /* global parameters used in these calculations */
diff --git a/sPuReMD/src/testmd.c b/sPuReMD/src/testmd.c
index 620b8a010feaac83a66de5985076826cee7b5b6a..9ef277465c6ff5ca5235087a91cc75a29bfcd9e4 100644
--- a/sPuReMD/src/testmd.c
+++ b/sPuReMD/src/testmd.c
@@ -223,5 +223,5 @@ int main( int argc, char* argv[] )
 
     sfree( lists, "main::lists" );
 
-    return SUCCESS;
+    return 0;
 }