diff --git a/PG-PuReMD/src/charges.c b/PG-PuReMD/src/charges.c
index 9c87b6ea712772eb4f09c143c42455e47c4d5028..31f6889dc2019daa339e6c431e02696f90702721 100644
--- a/PG-PuReMD/src/charges.c
+++ b/PG-PuReMD/src/charges.c
@@ -30,9 +30,9 @@
 #include "tool_box.h"
 
 
-int compare_matrix_entry(const void *v1, const void *v2)
+int compare_matrix_entry( const void *v1, const void *v2 )
 {
-    return ((sparse_matrix_entry *)v1)->j - ((sparse_matrix_entry *)v2)->j;
+    return ((sparse_matrix_entry *) v1)->j - ((sparse_matrix_entry *) v2)->j;
 }
 
 
@@ -44,487 +44,891 @@ void Sort_Matrix_Rows( sparse_matrix *A )
     {
         si = A->start[i];
         ei = A->end[i];
-        qsort( &(A->entries[si]), ei - si,
+        qsort( &A->entries[si], ei - si,
                 sizeof(sparse_matrix_entry), compare_matrix_entry );
     }
 }
 
 
-void Calculate_Droptol( sparse_matrix *A, real *droptol, real dtol )
+static void Init_Linear_Solver( reax_system *system, simulation_data *data,
+        control_params *control, storage *workspace, mpi_datatypes *mpi_data )
 {
-    int i, j, k;
-    real val;
+    int i;
+    reax_atom *atom;
 
-    /* init droptol to 0 - not necessary for an upper-triangular A */
-    for ( i = 0; i < A->n; ++i )
+    /* initialize solution vectors for linear solves in charge method */
+    switch ( control->charge_method )
     {
-        droptol[i] = 0;
+        case QEQ_CM:
+            for ( i = 0; i < system->n; ++i )
+            {
+                atom = &system->my_atoms[i];
+
+                workspace->b_s[i] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+                workspace->b_t[i] = -1.0;
+                workspace->b[i][0] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+                workspace->b[i][1] = -1.0;
+            }
+            break;
+
+        case EE_CM:
+            for ( i = 0; i < system->n; ++i )
+            {
+                atom = &system->my_atoms[i];
+
+                workspace->b_s[i] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+
+                //TODO: check if unused (redundant)
+                workspace->b[i][0] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+            }
+
+            if ( system->my_rank == 0 )
+            {
+                workspace->b_s[system->n] = control->cm_q_net;
+                workspace->b[system->n][0] = control->cm_q_net;
+            }
+            break;
+
+        case ACKS2_CM:
+            for ( i = 0; i < system->n; ++i )
+            {
+                atom = &system->my_atoms[i];
+
+                workspace->b_s[i] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+
+                //TODO: check if unused (redundant)
+                workspace->b[i][0] = -1.0 * system->reax_param.sbp[ atom->type ].chi;
+            }
+
+            if ( system->my_rank == 0 )
+            {
+                workspace->b_s[system->n] = control->cm_q_net;
+                workspace->b[system->n][0] = control->cm_q_net;
+            }
+
+            for ( i = system->n + 1; i < system->N_cm; ++i )
+            {
+                atom = &system->my_atoms[i];
+
+                workspace->b_s[i] = 0.0;
+
+                //TODO: check if unused (redundant)
+                workspace->b[i][0] = 0.0;
+            }
+
+            if ( system->my_rank == 0 )
+            {
+                workspace->b_s[system->n] = control->cm_q_net;
+                workspace->b[system->n][0] = control->cm_q_net;
+            }
+            break;
+
+        default:
+            fprintf( stderr, "[ERROR] Unknown charge method type. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
     }
+}
 
-    /* calculate sqaure of the norm of each row */
-    for ( i = 0; i < A->n; ++i )
-    {
-        val = A->entries[A->start[i]].val; // diagonal entry
-        droptol[i] += val * val;
 
-        // only within my block
-        for ( k = A->start[i] + 1; A->entries[k].j < A->n; ++k )
+static void Extrapolate_Charges_QEq( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
+{
+    int i;
+    real s_tmp, t_tmp;
+
+    /* spline extrapolation for s & t */
+    //TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
+#ifdef _OPENMP
+    #pragma omp parallel for schedule(static) \
+        default(none) private(i, s_tmp, t_tmp)
+#endif
+    for ( i = 0; i < system->N_cm; ++i )
+    {
+        /* no extrapolation, previous solution as initial guess */
+        if ( control->cm_init_guess_extrap1 == 0 )
+        {
+            s_tmp = system->my_atoms[i].s[0];
+        }
+        /* linear */
+        else if ( control->cm_init_guess_extrap1 == 1 )
         {
-            j = A->entries[k].j;
-            val = A->entries[k].val;
+            s_tmp = 2.0 * system->my_atoms[i].s[0] - system->my_atoms[i].s[1];
+        }
+        /* quadratic */
+        else if ( control->cm_init_guess_extrap1 == 2 )
+        {
+            s_tmp = system->my_atoms[i].s[2] + 3.0 * (system->my_atoms[i].s[0] - system->my_atoms[i].s[1]);
+        }
+        /* cubic */
+        else if ( control->cm_init_guess_extrap1 == 3 )
+        {
+            s_tmp = 4.0 * (system->my_atoms[i].s[0] + system->my_atoms[i].s[2]) -
+                    (6.0 * system->my_atoms[i].s[1] + system->my_atoms[i].s[3]);
+        }
+        /* 4th order */
+        else if ( control->cm_init_guess_extrap1 == 4 )
+        {
+            s_tmp = 5.0 * (system->my_atoms[i].s[0] - system->my_atoms[i].s[3]) +
+                10.0 * (-1.0 * system->my_atoms[i].s[1] + system->my_atoms[i].s[2]) + system->my_atoms[i].s[4];
+        }
+        else
+        {
+            s_tmp = 0.0;
+        }
 
-            droptol[i] += val * val;
-            droptol[j] += val * val;
+        /* no extrapolation, previous solution as initial guess */
+        if ( control->cm_init_guess_extrap2 == 0 )
+        {
+            t_tmp = system->my_atoms[i].t[0];
+        }
+        /* linear */
+        else if ( control->cm_init_guess_extrap2 == 1 )
+        {
+            t_tmp = 2.0 * system->my_atoms[i].t[0] - system->my_atoms[i].t[1];
+        }
+        /* quadratic */
+        else if ( control->cm_init_guess_extrap2 == 2 )
+        {
+            t_tmp = system->my_atoms[i].t[2] + 3.0 * (system->my_atoms[i].t[0] - system->my_atoms[i].t[1]);
+        }
+        /* cubic */
+        else if ( control->cm_init_guess_extrap2 == 3 )
+        {
+            t_tmp = 4.0 * (system->my_atoms[i].t[0] + system->my_atoms[i].t[2]) -
+                (6.0 * system->my_atoms[i].t[1] + system->my_atoms[i].t[3]);
+        }
+        /* 4th order */
+        else if ( control->cm_init_guess_extrap2 == 4 )
+        {
+            t_tmp = 5.0 * (system->my_atoms[i].t[0] - system->my_atoms[i].t[3]) +
+                10.0 * (-1.0 * system->my_atoms[i].t[1] + system->my_atoms[i].t[2]) + system->my_atoms[i].t[4];
+        }
+        else
+        {
+            t_tmp = 0.0;
         }
-    }
 
-    /* calculate local droptol for each row */
-    //fprintf( stderr, "droptol: " );
-    for ( i = 0; i < A->n; ++i )
-    {
-        droptol[i] = SQRT( droptol[i] ) * dtol;
-        //fprintf( stderr, "%f\n", droptol[i] );
+        system->my_atoms[i].s[4] = system->my_atoms[i].s[3];
+        system->my_atoms[i].s[3] = system->my_atoms[i].s[2];
+        system->my_atoms[i].s[2] = system->my_atoms[i].s[1];
+        system->my_atoms[i].s[1] = system->my_atoms[i].s[0];
+        system->my_atoms[i].s[0] = s_tmp;
+        /* x is used as initial guess to solver */
+        workspace->x[i][0] = s_tmp;
+
+        system->my_atoms[i].t[4] = system->my_atoms[i].t[3];
+        system->my_atoms[i].t[3] = system->my_atoms[i].t[2];
+        system->my_atoms[i].t[2] = system->my_atoms[i].t[1];
+        system->my_atoms[i].t[1] = system->my_atoms[i].t[0];
+        system->my_atoms[i].t[0] = t_tmp;
+        /* x is used as initial guess to solver */
+        workspace->x[i][1] = t_tmp;
     }
-    //fprintf( stderr, "\n" );
 }
 
 
-int Estimate_LU_Fill( sparse_matrix *A, real *droptol )
+static void Extrapolate_Charges_EE( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
 {
-    int i, j, pj;
-    int fillin;
-    real val;
+    int i;
+    real s_tmp;
 
-    fillin = 0;
-    for ( i = 0; i < A->n; ++i )
+    /* spline extrapolation for s */
+    //TODO: good candidate for vectorization, avoid moving data with head pointer and circular buffer
+#ifdef _OPENMP
+    #pragma omp parallel for schedule(static) \
+        default(none) private(i, s_tmp)
+#endif
+    for ( i = 0; i < system->N_cm; ++i )
     {
-        for ( pj = A->start[i] + 1; A->entries[pj].j < A->n; ++pj )
+        /* no extrapolation, previous solution as initial guess */
+        if ( control->cm_init_guess_extrap1 == 0 )
         {
-            j = A->entries[pj].j;
-            val = A->entries[pj].val;
-            if ( FABS(val) > droptol[i] )
-                ++fillin;
+            s_tmp = system->my_atoms[i].s[0];
+        }
+        /* linear */
+        else if ( control->cm_init_guess_extrap1 == 1 )
+        {
+            s_tmp = 2.0 * system->my_atoms[i].s[0] - system->my_atoms[i].s[1];
+        }
+        /* quadratic */
+        else if ( control->cm_init_guess_extrap1 == 2 )
+        {
+            s_tmp = system->my_atoms[i].s[2] + 3.0 * (system->my_atoms[i].s[0] - system->my_atoms[i].s[1]);
+        }
+        /* cubic */
+        else if ( control->cm_init_guess_extrap1 == 3 )
+        {
+            s_tmp = 4.0 * (system->my_atoms[i].s[0] + system->my_atoms[i].s[2]) -
+                    (6.0 * system->my_atoms[i].s[1] + system->my_atoms[i].s[3]);
+        }
+        /* 4th order */
+        else if ( control->cm_init_guess_extrap1 == 4 )
+        {
+            s_tmp = 5.0 * (system->my_atoms[i].s[0] - system->my_atoms[i].s[3]) +
+                10.0 * (-1.0 * system->my_atoms[i].s[1] + system->my_atoms[i].s[2]) + system->my_atoms[i].s[4];
+        }
+        else
+        {
+            s_tmp = 0.0;
         }
-    }
 
-    return fillin + A->n;
+        system->my_atoms[i].s[4] = system->my_atoms[i].s[3];
+        system->my_atoms[i].s[3] = system->my_atoms[i].s[2];
+        system->my_atoms[i].s[2] = system->my_atoms[i].s[1];
+        system->my_atoms[i].s[1] = system->my_atoms[i].s[0];
+        system->my_atoms[i].s[0] = s_tmp;
+        /* x is used as initial guess to solver */
+        workspace->x[i][0] = s_tmp;
+    }
 }
 
 
-void ICHOLT( sparse_matrix *A, real *droptol,
-        sparse_matrix *L, sparse_matrix *U )
+/* Compute preconditioner for QEq
+ */
+static void Compute_Preconditioner_QEq( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
 {
-    sparse_matrix_entry tmp[1000];
-    int i, j, pj, k1, k2, tmptop, Utop;
-    real val, dval;
-    int *Ltop;
-
-    Ltop = (int*) smalloc( A->n * sizeof(int), "ICHOLT::Ltop" );
+    sparse_matrix *Hptr;
 
-    // clear data structures
-    Utop = 0;
-    tmptop = 0;
-    for ( i = 0; i < A->n; ++i )
+    if ( control->cm_domain_sparsify_enabled == TRUE )
     {
-        U->start[i] = L->start[i] = L->end[i] = Ltop[i] = 0;
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
     }
 
-    for ( i = A->n - 1; i >= 0; --i )
+    switch ( control->cm_solver_pre_comp_type )
     {
-        U->start[i] = Utop;
-        tmptop = 0;
+        case NONE_PC:
+            break;
 
-        for ( pj = A->end[i] - 1; A->entries[pj].j >= A->n; --pj ); // skip ghosts
-        for ( ; pj > A->start[i]; --pj )
-        {
-            j = A->entries[pj].j;
-            val = A->entries[pj].val;
-            fprintf( stderr, "i: %d, j: %d  val=%f ", i, j, val );
-            //fprintf( stdout, "%d %d %24.16f\n", 6540-i, 6540-j, val );
-            //fprintf( stdout, "%d %d %24.16f\n", 6540-j, 6540-i, val );
+        case DIAG_PC:
+            data->timing.cm_solver_pre_comp +=
+                diag_pre_comp( system, workspace->Hdia_inv );
+//                diag_pre_comp( Hptr, workspace->Hdia_inv );
+            break;
 
-            if ( FABS(val) > droptol[i] )
-            {
-                k1 = tmptop - 1;
-                k2 = U->start[j] + 1;
-                while ( k1 >= 0 && k2 < U->end[j] )
-                {
-                    if ( tmp[k1].j < U->entries[k2].j )
-                    {
-                        k1--;
-                    }
-                    else if ( tmp[k1].j > U->entries[k2].j )
-                    {
-                        k2++;
-                    }
-                    else
-                    {
-                        val -= (tmp[k1--].val * U->entries[k2++].val);
-                    }
-                }
-
-                // U matrix is upper triangular
-                val /= U->entries[U->start[j]].val;
-                fprintf( stderr, " newval=%f", val );
-                tmp[tmptop].j = j;
-                tmp[tmptop].val = val;
-                tmptop++;
-            }
-            fprintf( stderr, "\n" );
-        }
-        //fprintf( stderr, "i = %d - tmptop = %d\n", i, tmptop );
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
+            exit( INVALID_INPUT );
+            break;
 
-        // compute the ith diagonal in U
-        dval = A->entries[A->start[i]].val;
-        //fprintf( stdout, "%d %d %24.16f\n", 6540-i, 6540-i, dval );
-        for ( k1 = 0; k1 < tmptop; ++k1 )
-        {
-            //if( FABS(tmp[k1].val) > droptol[i] )
-            dval -= SQR(tmp[k1].val);
-        }
-        dval = SQRT(dval);
-        // keep the diagonal in any case
-        U->entries[Utop].j = i;
-        U->entries[Utop].val = dval;
-        Utop++;
-
-        fprintf(stderr, "row%d: droptol=%.15f val=%.15f\n", i, droptol[i], dval);
-        for ( k1 = tmptop - 1; k1 >= 0; --k1 )
-        {
-            // apply the dropping rule once again
-            if ( FABS(tmp[k1].val) > droptol[i] / dval )
-            {
-                U->entries[Utop].j = tmp[k1].j;
-                U->entries[Utop].val = tmp[k1].val;
-                Utop++;
-                Ltop[tmp[k1].j]++;
-                fprintf( stderr, "%d(%.15f)\n", tmp[k1].j, tmp[k1].val );
-            }
-        }
+        case SAI_PC:
+#if defined(HAVE_LAPACKE) || defined(HAVE_LAPACKE_MKL)
+            //TODO: implement
+//            data->timing.cm_solver_pre_comp +=
+//                sparse_approx_inverse( workspace->H_full, workspace->H_spar_patt_full,
+//                        &workspace->H_app_inv );
+#else
+            fprintf( stderr, "[ERROR] LAPACKE support disabled. Re-compile before enabling. Terminating...\n" );
+            exit( INVALID_INPUT );
+#endif
+            break;
 
-        U->end[i] = Utop;
-        //fprintf( stderr, "i = %d - Utop = %d\n", i, Utop );
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
     }
+}
+
 
-#if defined(DEBUG)
-    // print matrix U
-    fprintf( stderr, "nnz(U): %d\n", Utop );
-    for ( i = 0; i < U->n; ++i )
+/* Compute preconditioner for EE
+ */
+static void Compute_Preconditioner_EE( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
+{
+    sparse_matrix *Hptr;
+
+    if ( control->cm_domain_sparsify_enabled == TRUE )
     {
-        fprintf( stderr, "row%d: ", i );
-        for ( pj = U->start[i]; pj < U->end[i]; ++pj )
-        {
-            fprintf( stderr, "%d ", U->entries[pj].j );
-        }
-        fprintf( stderr, "\n" );
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
     }
+
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 1.0;
+    }
+    
+    switch ( control->cm_solver_pre_comp_type )
+    {
+        case NONE_PC:
+            break;
+
+        case DIAG_PC:
+            data->timing.cm_solver_pre_comp +=
+                diag_pre_comp( system, workspace->Hdia_inv );
+//                diag_pre_comp( Hptr, workspace->Hdia_inv );
+            break;
+
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
+            exit( INVALID_INPUT );
+            break;
+
+        case SAI_PC:
+#if defined(HAVE_LAPACKE) || defined(HAVE_LAPACKE_MKL)
+            //TODO: implement
+//            data->timing.cm_solver_pre_comp +=
+//                sparse_approx_inverse( workspace->H_full, workspace->H_spar_patt_full,
+//                        &workspace->H_app_inv );
+#else
+            fprintf( stderr, "[ERROR] LAPACKE support disabled. Re-compile before enabling. Terminating...\n" );
+            exit( INVALID_INPUT );
 #endif
+            break;
+
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
+    }
 
-    // transpose matrix U into L
-    L->start[0] = L->end[0] = 0;
-    for ( i = 1; i < L->n; ++i )
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
     {
-        L->start[i] = L->end[i] = L->start[i - 1] + Ltop[i - 1] + 1;
-        //fprintf( stderr, "i=%d  L->start[i]=%d\n", i, L->start[i] );
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 0.0;
     }
+}
+
 
-    for ( i = 0; i < U->n; ++i )
+/* Compute preconditioner for ACKS2
+ */
+static void Compute_Preconditioner_ACKS2( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
+{
+    sparse_matrix *Hptr;
+
+    if ( control->cm_domain_sparsify_enabled == TRUE )
     {
-        for ( pj = U->start[i]; pj < U->end[i]; ++pj )
-        {
-            j = U->entries[pj].j;
-            L->entries[L->end[j]].j = i;
-            L->entries[L->end[j]].val = U->entries[pj].val;
-            L->end[j]++;
-        }
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
     }
 
-#if defined(DEBUG)
-    // print matrix L
-    fprintf( stderr, "nnz(L): %d\n", L->end[L->n - 1] );
-    for ( i = 0; i < L->n; ++i )
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
     {
-        fprintf( stderr, "row%d: ", i );
-        for ( pj = L->start[i]; pj < L->end[i]; ++pj )
-        {
-            fprintf( stderr, "%d ", L->entries[pj].j );
-        }
-        fprintf( stderr, "\n" );
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 1.0;
+        Hptr->entries[Hptr->start[system->N_cm] - 1].val = 1.0;
     }
+    
+    switch ( control->cm_solver_pre_comp_type )
+    {
+        case NONE_PC:
+            break;
+
+        case DIAG_PC:
+            data->timing.cm_solver_pre_comp +=
+                diag_pre_comp( system, workspace->Hdia_inv );
+//                diag_pre_comp( Hptr, workspace->Hdia_inv );
+            break;
+
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
+            exit( INVALID_INPUT );
+            break;
+
+        case SAI_PC:
+#if defined(HAVE_LAPACKE) || defined(HAVE_LAPACKE_MKL)
+            //TODO: implement
+//            data->timing.cm_solver_pre_comp +=
+//                sparse_approx_inverse( workspace->H_full, workspace->H_spar_patt_full,
+//                        &workspace->H_app_inv );
+#else
+            fprintf( stderr, "[ERROR] LAPACKE support disabled. Re-compile before enabling. Terminating...\n" );
+            exit( INVALID_INPUT );
 #endif
+            break;
 
-    sfree( Ltop, "Ltop" );
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
+    }
+
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 0.0;
+        Hptr->entries[Hptr->start[system->N_cm] - 1].val = 0.0;
+    }
 }
 
 
-void Init_MatVec( reax_system *system, simulation_data *data,
-        control_params *control, storage *workspace, mpi_datatypes *mpi_data )
+static void Setup_Preconditioner_QEq( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
 {
-    int i; //, fillin;
-    reax_atom *atom;
+    real time;
+    sparse_matrix *Hptr;
 
-//    if( (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0 ||
-//            workspace->L == NULL )
-//    {
-////        Print_Linear_System( system, control, workspace, data->step );
-//        Sort_Matrix_Rows( workspace->H );
-//        fprintf( stderr, "H matrix sorted\n" );
-//
-//        Calculate_Droptol( workspace->H, workspace->droptol, control->droptol );
-//        fprintf( stderr, "drop tolerances calculated\n" );
-//
-//        if( workspace->L == NULL )
-//        {
-//            fillin = Estimate_LU_Fill( workspace->H, workspace->droptol );
-//
-//            if( Allocate_Matrix( &(workspace->L), workspace->H->cap, fillin ) == 0 ||
-//                    Allocate_Matrix( &(workspace->U), workspace->H->cap, fillin ) == 0 )
-//            {
-//                fprintf( stderr, "not enough memory for LU matrices. terminating.\n" );
-//                MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
-//            }
-//
-//            workspace->L->n = workspace->H->n;
-//            workspace->U->n = workspace->H->n;
-//
-//#if defined(DEBUG_FOCUS)
-//            fprintf( stderr, "p%d: n=%d, fillin = %d\n",x
-//                    system->my_rank, workspace->L->n, fillin );
-//            fprintf( stderr, "p%d: allocated memory: L = U = %ldMB\n",
-//                    system->my_rank,fillin*sizeof(sparse_matrix_entry)/(1024*1024) );
-//#endif
-//        }
-//
-//      ICHOLT( workspace->H, workspace->droptol, workspace->L, workspace->U );
-//#if defined(DEBUG_FOCUS)
-//    fprintf( stderr, "p%d: icholt finished\n", system->my_rank );
-////    sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
-////    Print_Sparse_Matrix2( workspace->L, fname );
-////    Print_Sparse_Matrix( U );
-//#endif
-//    }
+    if ( control->cm_domain_sparsify_enabled == TRUE )
+    {
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
+    }
 
-    for ( i = 0; i < system->n; ++i )
+    /* sort H needed for SpMV's in linear solver, H or H_sp needed for preconditioning */
+    time = Get_Time( );
+    Sort_Matrix_Rows( &workspace->H );
+    if ( control->cm_domain_sparsify_enabled == TRUE )
     {
-        atom = &( system->my_atoms[i] );
+        Sort_Matrix_Rows( &workspace->H_sp );
+    }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
-        /* initialize diagonal inverse preconditioner vectors */
-        workspace->Hdia_inv[i] = 1. / system->reax_param.sbp[ atom->type ].eta;
+    switch ( control->cm_solver_pre_comp_type )
+    {
+        case NONE_PC:
+            break;
 
-        /* linear extrapolation for s and for t */
-        // newQEq: no extrapolation!
-        //workspace->s[i] = 2 * atom->s[0] - atom->s[1]; //0;
-        //workspace->t[i] = 2 * atom->t[0] - atom->t[1]; //0;
-        //workspace->x[i][0] = 2 * atom->s[0] - atom->s[1]; //0;
-        //workspace->x[i][1] = 2 * atom->t[0] - atom->t[1]; //0;
+        case DIAG_PC:
+            if ( workspace->Hdia_inv == NULL )
+            {
+//                workspace->Hdia_inv = scalloc( Hptr->n, sizeof( real ),
+//                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
+            }
+            break;
 
-        /* quadratic extrapolation for s and t */
-        // workspace->s[i] = atom->s[2] + 3 * ( atom->s[0] - atom->s[1] );
-        // workspace->t[i] = atom->t[2] + 3 * ( atom->t[0] - atom->t[1] );
-        //workspace->x[i][0] = atom->s[2] + 3 * ( atom->s[0] - atom->s[1] );
-        workspace->x[i][1] = atom->t[2] + 3 * ( atom->t[0] - atom->t[1] );
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
+            exit( INVALID_INPUT );
+            break;
 
-        /* cubic extrapolation for s and t */
-        workspace->x[i][0] = 4 * (atom->s[0] + atom->s[2]) - (6 * atom->s[1] + atom->s[3]);
-        //workspace->x[i][1] = 4*(atom->t[0]+atom->t[2])-(6*atom->t[1]+atom->t[3]);
+        case SAI_PC:
+            //TODO: implement
+//            setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
+//                    &workspace->H_spar_patt_full, &workspace->H_app_inv,
+//                    control->cm_solver_pre_comp_sai_thres );
+            break;
 
-//        fprintf(stderr, "i=%d s=%f t=%f\n", i, workspace->s[i], workspace->t[i]);
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
     }
+}
 
-    /* initialize solution vectors for linear solves in charge method */
-    switch ( control->charge_method )
+
+/* Setup routines before computing the preconditioner for EE
+ */
+static void Setup_Preconditioner_EE( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
+{
+    real time;
+    sparse_matrix *Hptr;
+
+    if ( control->cm_domain_sparsify_enabled == TRUE )
     {
-        case QEQ_CM:
-            for ( i = 0; i < system->n; ++i )
-            {
-                atom = &( system->my_atoms[i] );
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
+    }
 
-                workspace->b_s[i] = -system->reax_param.sbp[ atom->type ].chi;
-                workspace->b_t[i] = -1.0;
-                workspace->b[i][0] = -system->reax_param.sbp[ atom->type ].chi;
-                workspace->b[i][1] = -1.0;
-            }
+    /* sorted H needed for SpMV's in linear solver, H or H_sp needed for preconditioning */
+    time = Get_Time( );
+    Sort_Matrix_Rows( &workspace->H );
+    if ( control->cm_domain_sparsify_enabled == TRUE )
+    {
+        Sort_Matrix_Rows( &workspace->H_sp );
+    }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
+
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 1.0;
+    }
+
+    switch ( control->cm_solver_pre_comp_type )
+    {
+        case NONE_PC:
             break;
 
-        case EE_CM:
-            for ( i = 0; i < system->n; ++i )
+        case DIAG_PC:
+            if ( workspace->Hdia_inv == NULL )
             {
-                atom = &( system->my_atoms[i] );
+//                workspace->Hdia_inv = scalloc( system->N_cm, sizeof( real ),
+//                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
+            }
+            break;
 
-                workspace->b_s[i] = -system->reax_param.sbp[ atom->type ].chi;
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
+            exit( INVALID_INPUT );
+            break;
 
-                //TODO: check if unused (redundant)
-                workspace->b[i][0] = -system->reax_param.sbp[ atom->type ].chi;
-            }
+        case SAI_PC:
+            //TODO: implement
+//            setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
+//                    &workspace->H_spar_patt_full, &workspace->H_app_inv,
+//                    control->cm_solver_pre_comp_sai_thres );
+            break;
 
-            if ( system->my_rank == 0 )
-            {
-                workspace->b_s[system->n] = control->cm_q_net;
-                workspace->b[system->n][0] = control->cm_q_net;
-            }
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
             break;
+    }
 
-        case ACKS2_CM:
-            for ( i = 0; i < system->n; ++i )
-            {
-                atom = &( system->my_atoms[i] );
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 0.0;
+    }
+}
 
-                workspace->b_s[i] = -system->reax_param.sbp[ atom->type ].chi;
 
-                //TODO: check if unused (redundant)
-                workspace->b[i][0] = -system->reax_param.sbp[ atom->type ].chi;
-            }
+/* Setup routines before computing the preconditioner for ACKS2
+ */
+static void Setup_Preconditioner_ACKS2( const reax_system * const system,
+        const control_params * const control,
+        simulation_data * const data, storage * const workspace )
+{
+    real time;
+    sparse_matrix *Hptr;
 
-            if ( system->my_rank == 0 )
-            {
-                workspace->b_s[system->n] = control->cm_q_net;
-                workspace->b[system->n][0] = control->cm_q_net;
-            }
+    if ( control->cm_domain_sparsify_enabled == TRUE )
+    {
+        Hptr = &workspace->H_sp;
+    }
+    else
+    {
+        Hptr = &workspace->H;
+    }
 
-            for ( i = system->n + 1; i < system->N_cm; ++i )
-            {
-                atom = &( system->my_atoms[i] );
+    /* sort H needed for SpMV's in linear solver, H or H_sp needed for preconditioning */
+    time = Get_Time( );
+    Sort_Matrix_Rows( &workspace->H );
+    if ( control->cm_domain_sparsify_enabled == TRUE )
+    {
+        Sort_Matrix_Rows( &workspace->H_sp );
+    }
+    data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
 
-                workspace->b_s[i] = 0.0;
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 1.0;
+        Hptr->entries[Hptr->start[system->N_cm] - 1].val = 1.0;
+    }
 
-                //TODO: check if unused (redundant)
-                workspace->b[i][0] = 0.0;
-            }
+    switch ( control->cm_solver_pre_comp_type )
+    {
+        case NONE_PC:
+            break;
 
-            if ( system->my_rank == 0 )
+        case DIAG_PC:
+            if ( workspace->Hdia_inv == NULL )
             {
-                workspace->b_s[system->n] = control->cm_q_net;
-                workspace->b[system->n][0] = control->cm_q_net;
+//                workspace->Hdia_inv = scalloc( Hptr->n, sizeof( real ),
+//                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
             }
             break;
 
-        default:
-            fprintf( stderr, "Unknown charge method type. Terminating...\n" );
+        case ICHOLT_PC:
+        case ILU_PAR_PC:
+        case ILUT_PAR_PC:
+        case ILU_SUPERLU_MT_PC:
+            fprintf( stderr, "[ERROR] Unsupported preconditioner computation method. Terminating...\n" );
             exit( INVALID_INPUT );
             break;
+
+        case SAI_PC:
+            //TODO: implement
+//            setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
+//                    &workspace->H_spar_patt_full, &workspace->H_app_inv,
+//                    control->cm_solver_pre_comp_sai_thres );
+            break;
+
+        default:
+            fprintf( stderr, "[ERROR] Unrecognized preconditioner computation method. Terminating...\n" );
+            exit( UNKNOWN_OPTION );
+            break;
+    }
+
+    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
+            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
+            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
+    {
+        Hptr->entries[Hptr->start[system->N + 1] - 1].val = 0.0;
+        Hptr->entries[Hptr->start[system->N_cm] - 1].val = 0.0;
     }
 }
 
 
-void Calculate_Charges( reax_system *system, storage *workspace,
-        mpi_datatypes *mpi_data )
+/* Combine ficticious charges s and t to get atomic charge q for QEq method
+ */
+static void Calculate_Charges_QEq( const reax_system * const system,
+        storage * const workspace, const mpi_datatypes * const mpi_data )
 {
     int i;
-    real u;//, s_sum, t_sum;
+    real u;
     rvec2 my_sum, all_sum;
-    reax_atom *atom;
     real *q;
 
-    q = (real*) smalloc( system->N * sizeof(real), "Calculate_Charges::q" );
+    q = smalloc( sizeof(real) * system->N, "Calculate_Charges_QEq::q" );
 
-    //s_sum = Parallel_Vector_Acc(workspace->s, system->n, mpi_data->world);
-    //t_sum = Parallel_Vector_Acc(workspace->t, system->n, mpi_data->world);
-    my_sum[0] = my_sum[1] = 0;
+    my_sum[0] = 0.0;
+    my_sum[1] = 0.0;
     for ( i = 0; i < system->n; ++i )
     {
         my_sum[0] += workspace->x[i][0];
         my_sum[1] += workspace->x[i][1];
     }
 
-#if defined(DEBUG)
-    fprintf( stderr, "Host : my_sum[0]: %f and %f \n", my_sum[0], my_sum[1] );
-#endif
-
     MPI_Allreduce( &my_sum, &all_sum, 2, MPI_DOUBLE, MPI_SUM, mpi_data->world );
 
     u = all_sum[0] / all_sum[1];
-
-#if defined(DEBUG)
-    fprintf( stderr, "Host : u: %f \n", u );
-#endif
-
     for ( i = 0; i < system->n; ++i )
     {
-        atom = &( system->my_atoms[i] );
-
-        /* compute charge based on s & t */
-        //atom->q = workspace->s[i] - u * workspace->t[i];
-        q[i] = atom->q = workspace->x[i][0] - u * workspace->x[i][1];
-
-        /* backup s & t */
-        atom->s[3] = atom->s[2];
-        atom->s[2] = atom->s[1];
-        atom->s[1] = atom->s[0];
-        //atom->s[0] = workspace->s[i];
-        atom->s[0] = workspace->x[i][0];
-
-        atom->t[3] = atom->t[2];
-        atom->t[2] = atom->t[1];
-        atom->t[1] = atom->t[0];
-        //atom->t[0] = workspace->t[i];
-        atom->t[0] = workspace->x[i][1];
+        q[i] = workspace->x[i][0] - u * workspace->x[i][1];
+        system->my_atoms[i].q = q[i];
     }
 
     Dist( system, mpi_data, q, REAL_PTR_TYPE, MPI_DOUBLE, real_packer );
+
     for ( i = system->n; i < system->N; ++i )
     {
         system->my_atoms[i].q = q[i];
     }
 
-    sfree( q, "Calculate_Charges::q" );
+    sfree( q, "Calculate_Charges_QEq::q" );
 }
 
 
-void QEq( reax_system *system, control_params *control, simulation_data *data,
-        storage *workspace, output_controls *out_control,
-        mpi_datatypes *mpi_data )
+/* Get atomic charge q for EE method
+ */
+static void Calculate_Charges_EE( const reax_system * const system,
+        storage * const workspace, const mpi_datatypes * const mpi_data )
 {
-    int s_matvecs, t_matvecs;
-
-    Init_MatVec( system, data, control, workspace, mpi_data );
-
-    //if( data->step == 50010 ) {
-    //  Print_Linear_System( system, control, workspace, data->step );
-    // }
+    int i;
 
-#if defined(DEBUG)
-    fprintf( stderr, "p%d: initialized qEq\n", system->my_rank );
-    //Print_Linear_System( system, control, workspace, data->step );
-#endif
+    for ( i = 0; i < system->N; ++i )
+    {
+        system->my_atoms[i].q = workspace->s[i];
+    }
+}
 
-    //MATRIX CHANGES
-    s_matvecs = dual_CG( system, workspace, &workspace->H, workspace->b,
-            control->cm_solver_q_err, workspace->x, mpi_data, out_control->log, data );
-    t_matvecs = 0;
-    //fprintf (stderr, "Host: First CG complated with iterations: %d \n", s_matvecs);
 
-    //s_matvecs = CG(system, workspace, workspace->H, workspace->b_s, //newQEq sCG
-    // control->cm_solver_q_err, workspace->s, mpi_data, out_control->log );
-    //s_matvecs = PCG( system, workspace, workspace->H, workspace->b_s,
-    //   control->cm_solver_q_err, workspace->L, workspace->U, workspace->s,
-    //   mpi_data, out_control->log );
+/* Main driver method for QEq kernel
+ *  1) init / setup routines for preconditioning of linear solver
+ *  2) compute preconditioner
+ *  3) extrapolate charges
+ *  4) perform 2 linear solves
+ *  5) compute atomic charges based on output of (4)
+ */
+static void QEq( reax_system * const system, control_params * const control,
+        simulation_data * const data, storage * const workspace,
+        const output_controls * const out_control,
+        const mpi_datatypes * const mpi_data )
+{
+    int iters;
 
-#if defined(DEBUG)
-    fprintf( stderr, "p%d: first CG completed\n", system->my_rank );
-#endif
+    Init_Linear_Solver( system, data, control, workspace, mpi_data );
 
-    //t_matvecs = CG(system, workspace, workspace->H, workspace->b_t, //newQEq sCG
-    // control->cm_solver_q_err, workspace->t, mpi_data, out_control->log );
-    //t_matvecs = PCG( system, workspace, workspace->H, workspace->b_t,
-    //   control->cm_solver_q_err, workspace->L, workspace->U, workspace->t,
-    //   mpi_data, out_control->log );
+    if ( control->cm_solver_pre_comp_refactor > 0 &&
+            ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
+        
+    {
+        Setup_Preconditioner_QEq( system, control, data, workspace );
 
-#if defined(DEBUG)
-    fprintf( stderr, "p%d: second CG completed\n", system->my_rank );
-#endif
+        Compute_Preconditioner_QEq( system, control, data, workspace );
+    }
 
-    Calculate_Charges( system, workspace, mpi_data );
+    Extrapolate_Charges_QEq( system, control, data, workspace );
 
-#if defined(DEBUG)
-    fprintf( stderr, "p%d: computed charges\n", system->my_rank );
-    //Print_Charges( system );
-#endif
+    switch ( control->cm_solver_type )
+    {
+    case CG_S:
+        iters = dual_CG( system, workspace, &workspace->H, workspace->b,
+                control->cm_solver_q_err, workspace->x, mpi_data, out_control->log, data );
+
+//        iters = CG( system, workspace, workspace->H, workspace->b_s, //newQEq sCG
+//                control->cm_solver_q_err, workspace->s, mpi_data, out_control->log );
+//        iters += PCG( system, workspace, workspace->H, workspace->b_s,
+//                control->cm_solver_q_err, workspace->L, workspace->U, workspace->s,
+//                mpi_data, out_control->log );
+        break;
+
+    case GMRES_S:
+    case GMRES_H_S:
+    case SDM_S:
+    case BiCGStab_S:
+    default:
+        fprintf( stderr, "[ERROR] Unrecognized QEq solver selection. Terminating...\n" );
+        exit( INVALID_INPUT );
+        break;
+    }
 
 #if defined(LOG_PERFORMANCE)
     if ( system->my_rank == MASTER_NODE )
     {
-        data->timing.s_matvecs += s_matvecs;
-        data->timing.t_matvecs += t_matvecs;
+        data->timing.cm_solver_iters += iters;
     }
 #endif
+
+    Calculate_Charges_QEq( system, workspace, mpi_data );
+}
+
+
+/* Main driver method for EE kernel
+ *  1) init / setup routines for preconditioning of linear solver
+ *  2) compute preconditioner
+ *  3) extrapolate charges
+ *  4) perform 1 linear solve
+ *  5) compute atomic charges based on output of (4)
+ */
+static void EE( reax_system * const system, control_params * const control,
+        simulation_data * const data, storage * const workspace,
+        const output_controls * const out_control,
+        const mpi_datatypes * const mpi_data )
+{
+    int iters;
+
+    Init_Linear_Solver( system, data, control, workspace, mpi_data );
+
+    if ( control->cm_solver_pre_comp_refactor > 0 &&
+            ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
+    {
+        Setup_Preconditioner_EE( system, control, data, workspace );
+
+        Compute_Preconditioner_EE( system, control, data, workspace );
+    }
+
+    Extrapolate_Charges_EE( system, control, data, workspace );
+
+    switch ( control->cm_solver_type )
+    {
+    case GMRES_S:
+    case GMRES_H_S:
+    case CG_S:
+    case SDM_S:
+    case BiCGStab_S:
+    default:
+        fprintf( stderr, "[ERROR] Unrecognized EE solver selection. Terminating...\n" );
+        exit( INVALID_INPUT );
+        break;
+    }
+
+    data->timing.cm_solver_iters += iters;
+
+    Calculate_Charges_EE( system, workspace, mpi_data );
+}
+
+
+/* Main driver method for ACKS2 kernel
+ *  1) init / setup routines for preconditioning of linear solver
+ *  2) compute preconditioner
+ *  3) extrapolate charges
+ *  4) perform 1 linear solve
+ *  5) compute atomic charges based on output of (4)
+ */
+static void ACKS2( reax_system * const system, control_params * const control,
+        simulation_data * const data, storage * const workspace,
+        const output_controls * const out_control,
+        const mpi_datatypes * const mpi_data )
+{
+    int iters;
+
+    Init_Linear_Solver( system, data, control, workspace, mpi_data );
+
+    if ( control->cm_solver_pre_comp_refactor > 0 &&
+            ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
+    {
+        Setup_Preconditioner_ACKS2( system, control, data, workspace );
+
+        Compute_Preconditioner_ACKS2( system, control, data, workspace );
+    }
+
+    Extrapolate_Charges_EE( system, control, data, workspace );
+
+    switch ( control->cm_solver_type )
+    {
+    case GMRES_S:
+    case GMRES_H_S:
+    case CG_S:
+    case SDM_S:
+    case BiCGStab_S:
+    default:
+        fprintf( stderr, "[ERROR] Unrecognized ACKS2 solver selection. Terminating...\n" );
+        exit( INVALID_INPUT );
+        break;
+    }
+
+    data->timing.cm_solver_iters += iters;
+
+    Calculate_Charges_EE( system, workspace, mpi_data );
+}
+
+
+void Compute_Charges( reax_system * const system, control_params * const control,
+        simulation_data * const data, storage * const workspace,
+        const output_controls * const out_control,
+        const mpi_datatypes * const mpi_data )
+{
+    switch ( control->charge_method )
+    {
+    case QEQ_CM:
+        QEq( system, control, data, workspace, out_control, mpi_data );
+        break;
+
+    case EE_CM:
+        EE( system, control, data, workspace, out_control, mpi_data );
+        break;
+
+    case ACKS2_CM:
+        ACKS2( system, control, data, workspace, out_control, mpi_data );
+        break;
+
+    default:
+        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
+        exit( UNKNOWN_OPTION );
+        break;
+    }
 }
diff --git a/PG-PuReMD/src/charges.h b/PG-PuReMD/src/charges.h
index 08af5641406e9cb7d57d260701dbf6df702e4e47..83ffe17d31f9092c6b2eac2e26445a0942b966b6 100644
--- a/PG-PuReMD/src/charges.h
+++ b/PG-PuReMD/src/charges.h
@@ -29,8 +29,9 @@
 extern "C" {
 #endif
 
-void QEq( reax_system*, control_params*, simulation_data*,
-        storage*, output_controls*, mpi_datatypes* );
+void Compute_Charges( reax_system * const, control_params * const,
+        simulation_data * const, storage * const,
+        const output_controls * const, const mpi_datatypes * const );
 
 #ifdef __cplusplus
 }
diff --git a/PG-PuReMD/src/control.c b/PG-PuReMD/src/control.c
index 863032b5a7cb75b1681a9e9cf8d5294c874d8292..de09e775ce21a11e6eae9656f0646eea148e49f7 100644
--- a/PG-PuReMD/src/control.c
+++ b/PG-PuReMD/src/control.c
@@ -85,6 +85,8 @@ void Read_Control_File( char *control_file, control_params* control,
     control->cm_solver_restart = 50;
     control->cm_solver_q_err = 0.000001;
     control->cm_domain_sparsify_enabled = FALSE;
+    control->cm_init_guess_extrap1 = 3;
+    control->cm_init_guess_extrap2 = 2;
     control->cm_domain_sparsity = 1.0;
     control->cm_solver_pre_comp_type = DIAG_PC;
     control->cm_solver_pre_comp_sweeps = 3;
@@ -309,6 +311,16 @@ void Read_Control_File( char *control_file, control_params* control,
                     control->cm_domain_sparsify_enabled = TRUE;
                 }
             }
+            else if ( strncmp(tmp[0], "cm_init_guess_extrap1", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_init_guess_extrap1 = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_init_guess_extrap2", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_init_guess_extrap2 = ival;
+            }
             else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
             {
                 ival = atoi( tmp[1] );
@@ -329,6 +341,11 @@ void Read_Control_File( char *control_file, control_params* control,
                 ival = atoi( tmp[1] );
                 control->cm_solver_pre_comp_sweeps = ival;
             }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_sai_thres", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_solver_pre_comp_sai_thres = val;
+            }
             else if ( strcmp(tmp[0], "cm_solver_pre_app_type") == 0 )
             {
                 ival = atoi( tmp[1] );
diff --git a/PG-PuReMD/src/cuda/cuda_charges.cu b/PG-PuReMD/src/cuda/cuda_charges.cu
index c22058e9e81af68be9e7da236ca27c8e62eae10c..bf67d73410d29cb6b57cdcf2e1f8a55ce25b8a61 100644
--- a/PG-PuReMD/src/cuda/cuda_charges.cu
+++ b/PG-PuReMD/src/cuda/cuda_charges.cu
@@ -250,7 +250,7 @@ void Cuda_QEq( reax_system *system, control_params *control, simulation_data
         *data, storage *workspace, output_controls *out_control, mpi_datatypes
         *mpi_data )
 {
-    int s_matvecs, t_matvecs;
+    int iters;
 
     Cuda_Init_MatVec( system, workspace );
 
@@ -271,10 +271,9 @@ void Cuda_QEq( reax_system *system, control_params *control, simulation_data
         break;
 
     case CG_S:
-        s_matvecs = Cuda_dual_CG( system, control, workspace, &dev_workspace->H,
+        iters = Cuda_dual_CG( system, control, workspace, &dev_workspace->H,
                 dev_workspace->b, control->cm_solver_q_err, dev_workspace->x, mpi_data,
                 out_control->log, data );
-        t_matvecs = 0;
         break;
 
 
@@ -289,8 +288,7 @@ void Cuda_QEq( reax_system *system, control_params *control, simulation_data
 #if defined(LOG_PERFORMANCE)
     if ( system->my_rank == MASTER_NODE )
     {
-        data->timing.s_matvecs += s_matvecs;
-        data->timing.t_matvecs += t_matvecs;
+        data->timing.cm_solver_iters += iters;
     }
 #endif
 }
@@ -300,7 +298,7 @@ void Cuda_EE( reax_system *system, control_params *control, simulation_data
         *data, storage *workspace, output_controls *out_control, mpi_datatypes
         *mpi_data )
 {
-    int s_matvecs, t_matvecs;
+    int iters;
 
     Cuda_Init_MatVec( system, workspace );
 
@@ -314,9 +312,8 @@ void Cuda_EE( reax_system *system, control_params *control, simulation_data
         break;
 
     case CG_S:
-        s_matvecs = Cuda_CG( system, control, workspace, &dev_workspace->H,
+        iters = Cuda_CG( system, control, workspace, &dev_workspace->H,
                 dev_workspace->b_s, control->cm_solver_q_err, dev_workspace->s, mpi_data );
-        t_matvecs = 0;
         break;
 
 
@@ -331,8 +328,7 @@ void Cuda_EE( reax_system *system, control_params *control, simulation_data
 #if defined(LOG_PERFORMANCE)
     if ( system->my_rank == MASTER_NODE )
     {
-        data->timing.s_matvecs += s_matvecs;
-        data->timing.t_matvecs += t_matvecs;
+        data->timing.cm_solver_iters += iters;
     }
 #endif
 }
@@ -342,7 +338,7 @@ void Cuda_ACKS2( reax_system *system, control_params *control, simulation_data
         *data, storage *workspace, output_controls *out_control, mpi_datatypes
         *mpi_data )
 {
-    int s_matvecs, t_matvecs;
+    int iters;
 
     Cuda_Init_MatVec( system, workspace );
 
@@ -356,14 +352,13 @@ void Cuda_ACKS2( reax_system *system, control_params *control, simulation_data
         break;
 
     case CG_S:
-        s_matvecs = Cuda_CG( system, control, workspace, &dev_workspace->H,
+        iters = Cuda_CG( system, control, workspace, &dev_workspace->H,
                 dev_workspace->b_s, control->cm_solver_q_err, dev_workspace->s, mpi_data );
-        t_matvecs = 0;
         break;
 
 
     default:
-        fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
+        fprintf( stderr, "[ERROR] Unrecognized QEq solver selection. Terminating...\n" );
         exit( INVALID_INPUT );
         break;
     }
@@ -373,8 +368,7 @@ void Cuda_ACKS2( reax_system *system, control_params *control, simulation_data
 #if defined(LOG_PERFORMANCE)
     if ( system->my_rank == MASTER_NODE )
     {
-        data->timing.s_matvecs += s_matvecs;
-        data->timing.t_matvecs += t_matvecs;
+        data->timing.cm_solver_iters += iters;
     }
 #endif
 }
diff --git a/PG-PuReMD/src/forces.c b/PG-PuReMD/src/forces.c
index f8b6e7f34e0272393bd50e3773cc83c5d5d76c29..491f888eff7921fab6168ffd2251d9dbf5d892dd 100644
--- a/PG-PuReMD/src/forces.c
+++ b/PG-PuReMD/src/forces.c
@@ -966,7 +966,7 @@ int Compute_Forces( reax_system *system, control_params *control,
 #if defined(PURE_REAX)
         if ( charge_flag == TRUE )
         {
-            QEq( system, control, data, workspace, out_control, mpi_data );
+            Compute_Charges( system, control, data, workspace, out_control, mpi_data );
         }
 
 #if defined(LOG_PERFORMANCE)
diff --git a/PG-PuReMD/src/io_tools.c b/PG-PuReMD/src/io_tools.c
index e0e763afc3d3782b065f804a4128367c7974a0d3..7f704aeee1c572a5797ddd78bd383327dc5da505 100644
--- a/PG-PuReMD/src/io_tools.c
+++ b/PG-PuReMD/src/io_tools.c
@@ -1343,7 +1343,7 @@ void Output_Results( reax_system *system, control_params *control,
                     data->timing.nbrs * denom, data->timing.init_forces * denom,
                     data->timing.bonded * denom, data->timing.nonb * denom,
                     data->timing.cm * denom,
-                    (int)((data->timing.s_matvecs + data->timing.t_matvecs) * denom),
+                    (int)(data->timing.cm_solver_iters * denom),
                     data->timing.num_retries );
 
             Reset_Timing( &(data->timing) );
diff --git a/PG-PuReMD/src/lin_alg.c b/PG-PuReMD/src/lin_alg.c
index f16d9e235e52f61a34dc23b2368ec4b3a1878c71..7686e2d3dceddab5baa5dea84c4f2611f4295a0d 100644
--- a/PG-PuReMD/src/lin_alg.c
+++ b/PG-PuReMD/src/lin_alg.c
@@ -73,6 +73,31 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
 }
 
 
+/* Diagonal (Jacobi) preconditioner computation */
+real diag_pre_comp( const reax_system * const system, real * const Hdia_inv )
+{
+    unsigned int i;
+    real start;
+
+    start = Get_Time( );
+
+    for ( i = 0; i < system->n; ++i )
+    {
+//        if ( H->entries[H->start[i + 1] - 1].val != 0.0 )
+//        {
+//            Hdia_inv[i] = 1.0 / H->entries[H->start[i + 1] - 1].val;
+            Hdia_inv[i] = 1.0 / system->reax_param.sbp[ system->my_atoms[i].type ].eta;
+//        }
+//        else
+//        {
+//            Hdia_inv[i] = 1.0;
+//        }
+    }
+
+    return Get_Timing_Info( start );
+}
+
+
 int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H, rvec2
         *b, real tol, rvec2 *x, mpi_datatypes* mpi_data, FILE *fout,
         simulation_data *data )
diff --git a/PG-PuReMD/src/lin_alg.h b/PG-PuReMD/src/lin_alg.h
index 4b7ba2f0368e81d498f93d91e33db074ebddd9ca..1047686d4e3bd60fdf7acd69d1333297dcf7146f 100644
--- a/PG-PuReMD/src/lin_alg.h
+++ b/PG-PuReMD/src/lin_alg.h
@@ -29,6 +29,9 @@
 extern "C" {
 #endif
 
+//real diag_pre_comp( const sparse_matrix * const, real * const );
+real diag_pre_comp( const reax_system * const, real * const );
+
 int GMRES( reax_system*, storage*, sparse_matrix*,
         real*, real, real*, mpi_datatypes*, FILE* );
 
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index e0ffa7bb37bcf18d73e2383074d6b0374108a328..f984eb3963c238169d4e80e6b4a37eadac426776 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -315,7 +315,6 @@
 #endif
 
 
-/******************* ENUMERATIONS *************************/
 /* ensemble type */
 enum ensemble
 {
@@ -331,11 +330,11 @@ enum ensemble
 /* interaction list type */
 enum lists
 {
-    BONDS = 0,
-    OLD_BONDS = 1,
-    THREE_BODIES = 2,
-    HBONDS = 3,
-    FAR_NBRS = 4,
+    FAR_NBRS = 0,
+    BONDS = 1,
+    OLD_BONDS = 2,
+    THREE_BODIES = 3,
+    HBONDS = 4,
     DBOS = 5,
     DDELTAS = 6,
     LIST_N = 7,
@@ -422,6 +421,7 @@ enum solver
     GMRES_H_S = 1,
     CG_S = 2,
     SDM_S = 3,
+    BiCGStab_S = 4,
 };
 
 enum pre_comp
@@ -432,6 +432,7 @@ enum pre_comp
     ILU_PAR_PC = 3,
     ILUT_PAR_PC = 4,
     ILU_SUPERLU_MT_PC = 5,
+    SAI_PC = 6,
 };
 
 enum pre_app
@@ -480,13 +481,6 @@ enum traj_methods
     TF_N = 2,
 };
 
-/* ??? */
-enum molecules
-{
-    UNKNOWN = 0,
-    WATER = 1,
-};
-
 
 /********************** TYPE DEFINITIONS ********************/
 /* 3D vector, integer values */
@@ -1390,27 +1384,34 @@ typedef struct
     /* flag to control if force computations are tablulated */
     int tabulate;
 
-    /**/
+    /* method for computing atomic charges */
     unsigned int charge_method;
     /* frequency (in terms of simulation time steps) at which to
      * re-compute atomic charge distribution */
     int charge_freq;
-    /**/
+    /* iterative linear solver type */
     unsigned int cm_solver_type;
-    /**/
+    /* system net charge */
     real cm_q_net;
-    /**/
+    /* max. iterations for linear solver */
     unsigned int cm_solver_max_iters;
-    /**/
+    /* max. iterations before restarting in specific solvers, e.g., GMRES(k) */
     unsigned int cm_solver_restart;
     /* error tolerance of solution produced by charge distribution
      * sparse iterative linear solver */
     real cm_solver_q_err;
-    /**/
+    /* ratio used in computing sparser charge matrix,
+     * between 0.0 and 1.0 */
     real cm_domain_sparsity;
-    /**/
+    /* TRUE if enabled, FALSE otherwise */
     unsigned int cm_domain_sparsify_enabled;
-    /**/
+    /* order of spline extrapolation used for computing initial guess
+     * to linear solver */
+    unsigned int cm_init_guess_extrap1;
+    /* order of spline extrapolation used for computing initial guess
+     * to linear solver */
+    unsigned int cm_init_guess_extrap2;
+    /* preconditioner type for linear solver */
     unsigned int cm_solver_pre_comp_type;
     /* frequency (in terms of simulation time steps) at which to recompute
      * incomplete factorizations */
@@ -1418,11 +1419,17 @@ typedef struct
     /* drop tolerance of incomplete factorization schemes (ILUT, ICHOLT, etc.)
      * used for preconditioning the iterative linear solver used in charge distribution */
     real cm_solver_pre_comp_droptol;
-    /**/
+    /* num. of sweeps for computing preconditioner factors
+     * in fine-grained iterative methods (FG-ICHOL, FG-ILU) */
     unsigned int cm_solver_pre_comp_sweeps;
-    /**/
+    /* relative num. of non-zeros to charge matrix used to
+     * compute the sparse approximate inverse preconditioner,
+     * between 0.0 and 1.0 */
+    real cm_solver_pre_comp_sai_thres;
+    /* preconditioner application type */
     unsigned int cm_solver_pre_app_type;
-    /**/
+    /* num. of iterations used to apply preconditioner via
+     * Jacobi relaxation scheme (truncated Neumann series) */
     unsigned int cm_solver_pre_app_jacobi_iters;
 
     /* initial temperature of simulation, in Kelvin */
@@ -1544,12 +1551,11 @@ typedef struct
     real end;
     /* total elapsed time of event */
     real elapsed;
-
     /* total simulation time */
     real total;
     /* communication time */
     real comm;
-    /* neighbor (i.e., Verlet) list generation time */
+    /* neighbor list generation time */
     real nbrs;
     /* force initialization time */
     real init_forces;
@@ -1559,10 +1565,22 @@ typedef struct
     real nonb;
     /* atomic charge distribution calculation time */
     real cm;
-    /* num. of steps in iterative linear solver for charge distribution (QEq, first solve) */
-    int s_matvecs;
-    /* num. of steps in iterative linear solver for charge distribution (QEq, second solve) */
-    int t_matvecs;
+    /**/
+    real cm_sort_mat_rows;
+    /**/
+    real cm_solver_pre_comp;
+    /**/
+    real cm_solver_pre_app;
+    /* num. of steps in iterative linear solver for charge distribution */
+    int cm_solver_iters;
+    /**/
+    real cm_solver_spmv;
+    /**/
+    real cm_solver_vector_ops;
+    /**/
+    real cm_solver_orthog;
+    /**/
+    real cm_solver_tri_solve;
     /* num. of retries in main sim. loop */
     int num_retries;
 } reax_timing;
@@ -1948,15 +1966,31 @@ typedef struct
     /**/
     int *bond_mark;
 
-    /* charge matrix storage */
+    /* charge method storage */
     /* charge matrix */
     sparse_matrix H;
-    /* preconditioner */
+    /* charge matrix (full) */
+    sparse_matrix H_full;
+    /* sparser charge matrix */
+    sparse_matrix H_sp;
+    /* permuted charge matrix (graph coloring) */
+    sparse_matrix H_p;
+    /* sparsity pattern of charge matrix, used in
+     * computing a sparse approximate inverse preconditioner */
+    sparse_matrix H_spar_patt;
+    /* sparsity pattern of charge matrix (full), used in
+     * computing a sparse approximate inverse preconditioner */
+    sparse_matrix H_spar_patt_full;
+    /* sparse approximate inverse preconditioner */
+    sparse_matrix H_app_inv;
+    /* incomplete Cholesky or LU preconditioner */
     sparse_matrix L;
-    /* preconditioner */
+    /* incomplete Cholesky or LU preconditioner */
     sparse_matrix U;
-    /* preconditioner */
+    /* Jacobi preconditioner */
     real *Hdia_inv;
+    /* row drop tolerences for incomplete Cholesky preconditioner */
+    real *droptol;
     /**/
     real *b_s;
     /**/
@@ -1970,8 +2004,6 @@ typedef struct
     /**/
     real *t;
     /**/
-    real *droptol;
-    /**/
     rvec2 *b;
     /**/
     rvec2 *x;
diff --git a/PG-PuReMD/src/reset_tools.c b/PG-PuReMD/src/reset_tools.c
index 75863438e29d91414656f89ff2b48604675933bf..79386be592d1e47db9230e3006494ecb72e2e46a 100644
--- a/PG-PuReMD/src/reset_tools.c
+++ b/PG-PuReMD/src/reset_tools.c
@@ -118,8 +118,7 @@ void Reset_Timing( reax_timing *rt )
     rt->bonded = 0.0;
     rt->nonb = 0.0;
     rt->cm = 0.0;
-    rt->s_matvecs = 0;
-    rt->t_matvecs = 0;
+    rt->cm_solver_iters = 0;
     rt->num_retries = 0;
 }
 
diff --git a/PG-PuReMD/src/tool_box.c b/PG-PuReMD/src/tool_box.c
index 1093a99e385dd12bee96077bef8bc027ba6ae732..36022bcd9f8b91d79f48b6bbf8e72984dbb636d4 100644
--- a/PG-PuReMD/src/tool_box.c
+++ b/PG-PuReMD/src/tool_box.c
@@ -223,13 +223,13 @@ int Get_Atom_Type( reax_interaction *reax_param, char *s )
 
 char *Get_Element( reax_system *system, int i )
 {
-    return &( system->reax_param.sbp[system->my_atoms[i].type].name[0] );
+    return &system->reax_param.sbp[system->my_atoms[i].type].name[0];
 }
 
 
 char *Get_Atom_Name( reax_system *system, int i )
 {
-    return &(system->my_atoms[i].name[0]);
+    return &system->my_atoms[i].name[0];
 }
 
 
@@ -254,12 +254,12 @@ int Tokenize( const char* s, char*** tok )
 {
     char test[MAX_LINE];
     char *sep = "\t \n!=";
-    char *word;
+    char *word, *saveptr;
     int count = 0;
 
     strncpy( test, s, MAX_LINE );
 
-    for ( word = strtok(test, sep); word; word = strtok(NULL, sep) )
+    for ( word = strtok_r(test, sep, &saveptr); word; word = strtok_r(NULL, sep, &saveptr) )
     {
         strncpy( (*tok)[count], word, MAX_LINE );
         count++;
@@ -302,7 +302,7 @@ void * smalloc( size_t n, const char *name )
     }
 
 #if defined(DEBUG)
-    fprintf( stderr, "[INFO] granted memory at address: %p\n", (void *) ptr );
+    fprintf( stderr, "[INFO] granted memory at address: %p\n", ptr );
     fflush( stderr );
 #endif
 
@@ -335,7 +335,7 @@ void * srealloc( void *ptr, size_t n, const char *name )
                 n, name );
     }
 
-    fprintf( stderr, "[INFO] requesting memory for %s\n", name );
+    fprintf( stderr, "[INFO] requesting %zu bytes for %s\n", n, name );
     fflush( stderr );
 #endif
 
@@ -351,7 +351,7 @@ void * srealloc( void *ptr, size_t n, const char *name )
     }
 
 #if defined(DEBUG)
-    fprintf( stderr, "[INFO] address: %p [SREALLOC]\n", (void *) new_ptr );
+    fprintf( stderr, "[INFO] granted memory at address: %p\n", new_ptr );
     fflush( stderr );
 #endif
 
@@ -379,7 +379,7 @@ void * scalloc( size_t n, size_t size, const char *name )
     }
 
 #if defined(DEBUG)
-    fprintf( stderr, "[INFO] requesting memory for %s\n", name );
+    fprintf( stderr, "[INFO] requesting %zu bytes for %s\n", n * size, name );
     fflush( stderr );
 #endif
 
@@ -393,7 +393,7 @@ void * scalloc( size_t n, size_t size, const char *name )
     }
 
 #if defined(DEBUG)
-    fprintf( stderr, "[INFO] address: %p [SCALLOC]\n", (void *) ptr );
+    fprintf( stderr, "[INFO] granted memory at address: %p\n", ptr );
     fflush( stderr );
 #endif
 
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index 1d0b9c36185f1b9b98c8a356c0c6eec59f845e61..c2a6a6ba1e6bb24cf25f2882c3a0f1e5b9331f5b 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -34,14 +34,9 @@
 
 enum atoms
 {
-    C_ATOM = 0,
-    H_ATOM = 1,
-    O_ATOM = 2,
-    N_ATOM = 3,
-    S_ATOM = 4,
-    SI_ATOM = 5,
-    GE_ATOM = 6,
-    X_ATOM = 7,
+    O_ATOM = 0,
+    N_ATOM = 1,
+    SI_ATOM = 2,
 };
 
 enum molecule_type
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 4ff5a3426ba83949b084c2e078ebe974311fee90..875ce97e916d0b836eb5a2f8a07638fdb2d12914 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -55,24 +55,24 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
         /* linear */
         else if ( control->cm_init_guess_extrap1 == 1 )
         {
-            s_tmp = 2 * workspace->s[0][i] - workspace->s[1][i];
+            s_tmp = 2.0 * workspace->s[0][i] - workspace->s[1][i];
         }
         /* quadratic */
         else if ( control->cm_init_guess_extrap1 == 2 )
         {
-            s_tmp = workspace->s[2][i] + 3 * (workspace->s[0][i]-workspace->s[1][i]);
+            s_tmp = workspace->s[2][i] + 3.0 * (workspace->s[0][i]-workspace->s[1][i]);
         }
         /* cubic */
         else if ( control->cm_init_guess_extrap1 == 3 )
         {
-            s_tmp = 4 * (workspace->s[0][i] + workspace->s[2][i]) -
-                    (6 * workspace->s[1][i] + workspace->s[3][i] );
+            s_tmp = 4.0 * (workspace->s[0][i] + workspace->s[2][i]) -
+                    (6.0 * workspace->s[1][i] + workspace->s[3][i]);
         }
         /* 4th order */
         else if ( control->cm_init_guess_extrap1 == 4 )
         {
-            s_tmp = 5 * (workspace->s[0][i] - workspace->s[3][i]) +
-                10 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
+            s_tmp = 5.0 * (workspace->s[0][i] - workspace->s[3][i]) +
+                10.0 * (-1.0 * workspace->s[1][i] + workspace->s[2][i]) + workspace->s[4][i];
         }
         else
         {
@@ -92,19 +92,19 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
         /* quadratic */
         else if ( control->cm_init_guess_extrap2 == 2 )
         {
-            t_tmp = workspace->t[2][i] + 3 * (workspace->t[0][i] - workspace->t[1][i]);
+            t_tmp = workspace->t[2][i] + 3.0 * (workspace->t[0][i] - workspace->t[1][i]);
         }
         /* cubic */
         else if ( control->cm_init_guess_extrap2 == 3 )
         {
             t_tmp = 4.0 * (workspace->t[0][i] + workspace->t[2][i]) -
-                (6.0 * workspace->t[1][i] + workspace->t[3][i] );
+                (6.0 * workspace->t[1][i] + workspace->t[3][i]);
         }
         /* 4th order */
         else if ( control->cm_init_guess_extrap2 == 4 )
         {
             t_tmp = 5.0 * (workspace->t[0][i] - workspace->t[3][i]) +
-                10.0 * (-workspace->t[1][i] + workspace->t[2][i] ) + workspace->t[4][i];
+                10.0 * (-1.0 * workspace->t[1][i] + workspace->t[2][i]) + workspace->t[4][i];
         }
         else
         {
@@ -186,8 +186,7 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
  */
 static void Compute_Preconditioner_QEq( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     real time;
     sparse_matrix *Hptr;
@@ -292,8 +291,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
  */
 //static void Compute_Preconditioner_EE( const reax_system * const system,
 //        const control_params * const control,
-//        simulation_data * const data, static_storage * const workspace,
-//        const reax_list * const far_nbrs )
+//        simulation_data * const data, static_storage * const workspace )
 //{
 //    int i, top;
 //    static real * ones = NULL, * x = NULL, * y = NULL;
@@ -536,8 +534,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
  */
 static void Compute_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     real time;
     sparse_matrix *Hptr;
@@ -668,8 +665,7 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
  */
 static void Compute_Preconditioner_ACKS2( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     real time;
     sparse_matrix *Hptr;
@@ -800,8 +796,7 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
 
 static void Setup_Preconditioner_QEq( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     int fillin;
     real time;
@@ -934,8 +929,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
  */
 static void Setup_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     int fillin;
     real time;
@@ -1086,8 +1080,7 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
  */
 static void Setup_Preconditioner_ACKS2( const reax_system * const system,
         const control_params * const control,
-        simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs )
+        simulation_data * const data, static_storage * const workspace )
 {
     int fillin;
     real time;
@@ -1291,7 +1284,7 @@ static void Calculate_Charges_EE( const reax_system * const system,
  */
 static void QEq( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs, const output_controls * const out_control )
+        const output_controls * const out_control )
 {
     int iters;
 
@@ -1299,9 +1292,9 @@ static void QEq( reax_system * const system, control_params * const control,
             ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
         
     {
-        Setup_Preconditioner_QEq( system, control, data, workspace, far_nbrs );
+        Setup_Preconditioner_QEq( system, control, data, workspace );
 
-        Compute_Preconditioner_QEq( system, control, data, workspace, far_nbrs );
+        Compute_Preconditioner_QEq( system, control, data, workspace );
     }
 
     Extrapolate_Charges_QEq( system, control, data, workspace );
@@ -1382,16 +1375,16 @@ static void QEq( reax_system * const system, control_params * const control,
  */
 static void EE( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs, const output_controls * const out_control )
+        const output_controls * const out_control )
 {
     int iters;
 
     if ( control->cm_solver_pre_comp_refactor > 0 &&
             ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
     {
-        Setup_Preconditioner_EE( system, control, data, workspace, far_nbrs );
+        Setup_Preconditioner_EE( system, control, data, workspace );
 
-        Compute_Preconditioner_EE( system, control, data, workspace, far_nbrs );
+        Compute_Preconditioner_EE( system, control, data, workspace );
     }
 
     Extrapolate_Charges_EE( system, control, data, workspace );
@@ -1454,16 +1447,16 @@ static void EE( reax_system * const system, control_params * const control,
  */
 static void ACKS2( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs, const output_controls * const out_control )
+        const output_controls * const out_control )
 {
     int iters;
 
     if ( control->cm_solver_pre_comp_refactor > 0 &&
             ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) )
     {
-        Setup_Preconditioner_ACKS2( system, control, data, workspace, far_nbrs );
+        Setup_Preconditioner_ACKS2( system, control, data, workspace );
 
-        Compute_Preconditioner_ACKS2( system, control, data, workspace, far_nbrs );
+        Compute_Preconditioner_ACKS2( system, control, data, workspace );
     }
 
 //   Print_Linear_System( system, control, workspace, data->step );
@@ -1533,7 +1526,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
 
 void Compute_Charges( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const reax_list * const far_nbrs, const output_controls * const out_control )
+        const output_controls * const out_control )
 {
 #if defined(DEBUG_FOCUS)
 #define SIZE (200)
@@ -1562,20 +1555,20 @@ void Compute_Charges( reax_system * const system, control_params * const control
     switch ( control->charge_method )
     {
     case QEQ_CM:
-        QEq( system, control, data, workspace, far_nbrs, out_control );
+        QEq( system, control, data, workspace, out_control );
         break;
 
     case EE_CM:
-        EE( system, control, data, workspace, far_nbrs, out_control );
+        EE( system, control, data, workspace, out_control );
         break;
 
     case ACKS2_CM:
-        ACKS2( system, control, data, workspace, far_nbrs, out_control );
+        ACKS2( system, control, data, workspace, out_control );
         break;
 
     default:
-        fprintf( stderr, "Invalid charge method. Terminating...\n" );
-        exit( INVALID_INPUT );
+        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
+        exit( UNKNOWN_OPTION );
         break;
     }
 }
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index 48028b0c8fe9fa355c7dcd6f658d00e08d0db74b..85ac44d32fca2a17955f1286a5546fd1a32a5074 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -26,8 +26,7 @@
 
 
 void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
-          static_storage* const, const reax_list* const,
-          const output_controls* const );
+          static_storage* const, const output_controls* const );
 
 
 #endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index bd95e4e9ed6de618f0a5100f80c6f71b946ec00d..baacc07b006e15a77ce32fc003c5a889f105c1ce 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -32,9 +32,8 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
         output_controls *out_control )
 {
     char *s, **tmp;
-    int i;
+    int c, i, ival;
     real val;
-    int ival;
 
     /* assign default values */
     strncpy( control->sim_name, "default.sim", MAX_STR );
@@ -149,410 +148,413 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
     /* read control parameters file */
     while ( fgets( s, MAX_LINE, fp ) )
     {
-        Tokenize( s, &tmp );
+        c = Tokenize( s, &tmp );
 
-        if ( strncmp(tmp[0], "simulation_name", MAX_LINE) == 0 )
+        if ( c > 0 )
         {
-            strncpy( control->sim_name, tmp[1], MAX_STR );
-        }
-        //else if( strncmp(tmp[0], "restart", MAX_LINE) == 0 ) {
-        //  ival = atoi(tmp[1]);
-        //  control->restart = ival;
-        //}
-        else if ( strncmp(tmp[0], "restart_format", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->restart_format = ival;
-        }
-        else if ( strncmp(tmp[0], "restart_freq", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->restart_freq = ival;
-        }
-        else if ( strncmp(tmp[0], "random_vel", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->random_vel = ival;
-        }
-        else if ( strncmp(tmp[0], "reposition_atoms", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->reposition_atoms = ival;
-        }
-        else if ( strncmp(tmp[0], "ensemble_type", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->ensemble = ival;
-        }
-        else if ( strncmp(tmp[0], "nsteps", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->nsteps = ival;
-        }
-        else if ( strncmp(tmp[0], "dt", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->dt = val * 1.e-3;  // convert dt from fs to ps!
-        }
-        else if ( strncmp(tmp[0], "periodic_boundaries", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->periodic_boundaries = ival;
-        }
-        else if ( strncmp(tmp[0], "geo_format", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->geo_format = ival;
-        }
-        else if ( strncmp(tmp[0], "restrict_bonds", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->restrict_bonds = ival;
-        }
-        else if ( strncmp(tmp[0], "tabulate_long_range", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->tabulate = ival;
-        }
-        else if ( strncmp(tmp[0], "reneighbor", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->reneighbor = ival;
-        }
-        else if ( strncmp(tmp[0], "vlist_buffer", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->vlist_cut = val;
-        }
-        else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->nbr_cut = val;
-        }
-        else if ( strncmp(tmp[0], "thb_cutoff", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->thb_cut = val;
-        }
-        else if ( strncmp(tmp[0], "hbond_cutoff", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->hb_cut = val;
-        }
-        else if ( strncmp(tmp[0], "charge_method", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->charge_method = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_q_net", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->cm_q_net = val;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_type", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_type = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_max_iters", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_max_iters = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_restart", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_restart = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_q_err", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->cm_solver_q_err = val;
-        }
-        else if ( strncmp(tmp[0], "cm_domain_sparsity", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->cm_domain_sparsity = val;
-            if ( val < 1.0 )
+            if ( strncmp(tmp[0], "simulation_name", MAX_LINE) == 0 )
             {
-                control->cm_domain_sparsify_enabled = TRUE;
+                strncpy( control->sim_name, tmp[1], MAX_STR );
             }
-        }
-        else if ( strncmp(tmp[0], "cm_init_guess_extrap1", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_init_guess_extrap1 = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_init_guess_extrap2", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_init_guess_extrap2 = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_comp_type", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_pre_comp_type = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_comp_refactor", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_pre_comp_refactor = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_comp_droptol", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->cm_solver_pre_comp_droptol = val;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_comp_sweeps", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_pre_comp_sweeps = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_comp_sai_thres", MAX_LINE) == 0 )
-        {
-            val = atof( tmp[1] );
-            control->cm_solver_pre_comp_sai_thres = val;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_app_type", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_pre_app_type = ival;
-        }
-        else if ( strncmp(tmp[0], "cm_solver_pre_app_jacobi_iters", MAX_LINE) == 0 )
-        {
-            ival = atoi( tmp[1] );
-            control->cm_solver_pre_app_jacobi_iters = ival;
-        }
-        else if ( strncmp(tmp[0], "temp_init", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->T_init = val;
-
-            if ( control->T_init < 0.001 )
+            //else if( strncmp(tmp[0], "restart", MAX_LINE) == 0 ) {
+            //  ival = atoi(tmp[1]);
+            //  control->restart = ival;
+            //}
+            else if ( strncmp(tmp[0], "restart_format", MAX_LINE) == 0 )
             {
-                control->T_init = 0.001;
+                ival = atoi(tmp[1]);
+                out_control->restart_format = ival;
             }
-        }
-        else if ( strncmp(tmp[0], "temp_final", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->T_final = val;
-
-            if ( control->T_final < 0.1 )
+            else if ( strncmp(tmp[0], "restart_freq", MAX_LINE) == 0 )
             {
-                control->T_final = 0.1;
+                ival = atoi(tmp[1]);
+                out_control->restart_freq = ival;
             }
-        }
-        else if ( strncmp(tmp[0], "t_mass", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            /* convert t_mass from fs to ps */
-            control->Tau_T = val * 1.e-3;
-        }
-        else if ( strncmp(tmp[0], "t_mode", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->T_mode = ival;
-        }
-        else if ( strncmp(tmp[0], "t_rate", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->T_rate = val;
-        }
-        else if ( strncmp(tmp[0], "t_freq", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->T_freq = val;
-        }
-        else if ( strncmp(tmp[0], "pressure", MAX_LINE) == 0 )
-        {
-            if ( control->ensemble == iNPT )
+            else if ( strncmp(tmp[0], "random_vel", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->random_vel = ival;
+            }
+            else if ( strncmp(tmp[0], "reposition_atoms", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->reposition_atoms = ival;
+            }
+            else if ( strncmp(tmp[0], "ensemble_type", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->ensemble = ival;
+            }
+            else if ( strncmp(tmp[0], "nsteps", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->nsteps = ival;
+            }
+            else if ( strncmp(tmp[0], "dt", MAX_LINE) == 0 )
             {
                 val = atof(tmp[1]);
-                control->P[0] = control->P[1] = control->P[2] = val;
+                control->dt = val * 1.e-3;  // convert dt from fs to ps!
+            }
+            else if ( strncmp(tmp[0], "periodic_boundaries", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->periodic_boundaries = ival;
+            }
+            else if ( strncmp(tmp[0], "geo_format", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->geo_format = ival;
+            }
+            else if ( strncmp(tmp[0], "restrict_bonds", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->restrict_bonds = ival;
+            }
+            else if ( strncmp(tmp[0], "tabulate_long_range", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->tabulate = ival;
+            }
+            else if ( strncmp(tmp[0], "reneighbor", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->reneighbor = ival;
             }
-            else if ( control->ensemble == sNPT )
+            else if ( strncmp(tmp[0], "vlist_buffer", MAX_LINE) == 0 )
             {
                 val = atof(tmp[1]);
-                control->P[0] = val;
-
-                val = atof(tmp[2]);
-                control->P[1] = val;
-
-                val = atof(tmp[3]);
-                control->P[2] = val;
+                control->vlist_cut = val;
             }
-        }
-        else if ( strncmp(tmp[0], "p_mass", MAX_LINE) == 0 )
-        {
-            if ( control->ensemble == iNPT )
+            else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
             {
                 val = atof(tmp[1]);
-                control->Tau_P[0] = val * 1.e-3;   // convert p_mass from fs to ps
+                control->nbr_cut = val;
             }
-            else if ( control->ensemble == sNPT )
+            else if ( strncmp(tmp[0], "thb_cutoff", MAX_LINE) == 0 )
             {
                 val = atof(tmp[1]);
-                control->Tau_P[0] = val * 1.e-3;   // convert p_mass from fs to ps
-
-                val = atof(tmp[2]);
-                control->Tau_P[1] = val * 1.e-3;   // convert p_mass from fs to ps
+                control->thb_cut = val;
+            }
+            else if ( strncmp(tmp[0], "hbond_cutoff", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->hb_cut = val;
+            }
+            else if ( strncmp(tmp[0], "charge_method", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->charge_method = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_q_net", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_q_net = val;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_type", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_type = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_max_iters", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_max_iters = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_restart", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_restart = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_q_err", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_solver_q_err = val;
+            }
+            else if ( strncmp(tmp[0], "cm_domain_sparsity", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_domain_sparsity = val;
+                if ( val < 1.0 )
+                {
+                    control->cm_domain_sparsify_enabled = TRUE;
+                }
+            }
+            else if ( strncmp(tmp[0], "cm_init_guess_extrap1", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_init_guess_extrap1 = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_init_guess_extrap2", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_init_guess_extrap2 = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_type", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_pre_comp_type = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_refactor", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_pre_comp_refactor = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_droptol", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_solver_pre_comp_droptol = val;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_sweeps", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_pre_comp_sweeps = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_comp_sai_thres", MAX_LINE) == 0 )
+            {
+                val = atof( tmp[1] );
+                control->cm_solver_pre_comp_sai_thres = val;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_app_type", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_pre_app_type = ival;
+            }
+            else if ( strncmp(tmp[0], "cm_solver_pre_app_jacobi_iters", MAX_LINE) == 0 )
+            {
+                ival = atoi( tmp[1] );
+                control->cm_solver_pre_app_jacobi_iters = ival;
+            }
+            else if ( strncmp(tmp[0], "temp_init", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->T_init = val;
 
-                val = atof(tmp[3]);
-                control->Tau_P[2] = val * 1.e-3;   // convert p_mass from fs to ps
+                if ( control->T_init < 0.001 )
+                {
+                    control->T_init = 0.001;
+                }
             }
-        }
-        else if ( strncmp(tmp[0], "pt_mass", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->Tau_PT = val * 1.e-3;  // convert pt_mass from fs to ps
-        }
-        else if ( strncmp(tmp[0], "compress", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->compressibility = val;
-        }
-        else if ( strncmp(tmp[0], "press_mode", MAX_LINE) == 0 )
-        {
-            val = atoi(tmp[1]);
-            control->press_mode = val;
-        }
-        else if ( strncmp(tmp[0], "remove_CoM_vel", MAX_LINE) == 0 )
-        {
-            val = atoi(tmp[1]);
-            control->remove_CoM_vel = val;
-        }
-        else if ( strncmp(tmp[0], "debug_level", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->debug_level = ival;
-        }
-        else if ( strncmp(tmp[0], "energy_update_freq", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->energy_update_freq = ival;
-        }
-        else if ( strncmp(tmp[0], "write_freq", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->write_steps = ival;
-        }
-        else if ( strncmp(tmp[0], "traj_compress", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->traj_compress = ival;
+            else if ( strncmp(tmp[0], "temp_final", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->T_final = val;
 
-            if ( out_control->traj_compress )
-                out_control->write = (int (*)(FILE *, const char *, ...)) gzprintf;
-            else out_control->write = fprintf;
-        }
-        else if ( strncmp(tmp[0], "traj_format", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->traj_format = ival;
+                if ( control->T_final < 0.1 )
+                {
+                    control->T_final = 0.1;
+                }
+            }
+            else if ( strncmp(tmp[0], "t_mass", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                /* convert t_mass from fs to ps */
+                control->Tau_T = val * 1.e-3;
+            }
+            else if ( strncmp(tmp[0], "t_mode", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->T_mode = ival;
+            }
+            else if ( strncmp(tmp[0], "t_rate", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->T_rate = val;
+            }
+            else if ( strncmp(tmp[0], "t_freq", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->T_freq = val;
+            }
+            else if ( strncmp(tmp[0], "pressure", MAX_LINE) == 0 )
+            {
+                if ( control->ensemble == iNPT )
+                {
+                    val = atof(tmp[1]);
+                    control->P[0] = control->P[1] = control->P[2] = val;
+                }
+                else if ( control->ensemble == sNPT )
+                {
+                    val = atof(tmp[1]);
+                    control->P[0] = val;
+
+                    val = atof(tmp[2]);
+                    control->P[1] = val;
+
+                    val = atof(tmp[3]);
+                    control->P[2] = val;
+                }
+            }
+            else if ( strncmp(tmp[0], "p_mass", MAX_LINE) == 0 )
+            {
+                if ( control->ensemble == iNPT )
+                {
+                    val = atof(tmp[1]);
+                    control->Tau_P[0] = val * 1.e-3;   // convert p_mass from fs to ps
+                }
+                else if ( control->ensemble == sNPT )
+                {
+                    val = atof(tmp[1]);
+                    control->Tau_P[0] = val * 1.e-3;   // convert p_mass from fs to ps
+
+                    val = atof(tmp[2]);
+                    control->Tau_P[1] = val * 1.e-3;   // convert p_mass from fs to ps
+
+                    val = atof(tmp[3]);
+                    control->Tau_P[2] = val * 1.e-3;   // convert p_mass from fs to ps
+                }
+            }
+            else if ( strncmp(tmp[0], "pt_mass", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->Tau_PT = val * 1.e-3;  // convert pt_mass from fs to ps
+            }
+            else if ( strncmp(tmp[0], "compress", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->compressibility = val;
+            }
+            else if ( strncmp(tmp[0], "press_mode", MAX_LINE) == 0 )
+            {
+                val = atoi(tmp[1]);
+                control->press_mode = val;
+            }
+            else if ( strncmp(tmp[0], "remove_CoM_vel", MAX_LINE) == 0 )
+            {
+                val = atoi(tmp[1]);
+                control->remove_CoM_vel = val;
+            }
+            else if ( strncmp(tmp[0], "debug_level", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->debug_level = ival;
+            }
+            else if ( strncmp(tmp[0], "energy_update_freq", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->energy_update_freq = ival;
+            }
+            else if ( strncmp(tmp[0], "write_freq", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->write_steps = ival;
+            }
+            else if ( strncmp(tmp[0], "traj_compress", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->traj_compress = ival;
 
-            if ( out_control->traj_format == 0 )
+                if ( out_control->traj_compress )
+                    out_control->write = (int (*)(FILE *, const char *, ...)) gzprintf;
+                else out_control->write = fprintf;
+            }
+            else if ( strncmp(tmp[0], "traj_format", MAX_LINE) == 0 )
             {
-                out_control->write_header =
-                    (int (*)( reax_system*, control_params*,
-                              static_storage*, void* )) Write_Custom_Header;
-                out_control->append_traj_frame =
-                    (int (*)(reax_system*, control_params*, simulation_data*,
-                             static_storage*, reax_list **, void*)) Append_Custom_Frame;
+                ival = atoi(tmp[1]);
+                out_control->traj_format = ival;
+
+                if ( out_control->traj_format == 0 )
+                {
+                    out_control->write_header =
+                        (int (*)( reax_system*, control_params*,
+                                  static_storage*, void* )) Write_Custom_Header;
+                    out_control->append_traj_frame =
+                        (int (*)(reax_system*, control_params*, simulation_data*,
+                                 static_storage*, reax_list **, void*)) Append_Custom_Frame;
+                }
+                else if ( out_control->traj_format == 1 )
+                {
+                    out_control->write_header =
+                        (int (*)( reax_system*, control_params*,
+                                  static_storage*, void* )) Write_xyz_Header;
+                    out_control->append_traj_frame =
+                        (int (*)( reax_system*,  control_params*, simulation_data*,
+                                  static_storage*, reax_list **, void* )) Append_xyz_Frame;
+                }
             }
-            else if ( out_control->traj_format == 1 )
+            else if ( strncmp(tmp[0], "traj_title", MAX_LINE) == 0 )
             {
-                out_control->write_header =
-                    (int (*)( reax_system*, control_params*,
-                              static_storage*, void* )) Write_xyz_Header;
-                out_control->append_traj_frame =
-                    (int (*)( reax_system*,  control_params*, simulation_data*,
-                              static_storage*, reax_list **, void* )) Append_xyz_Frame;
+                strncpy( out_control->traj_title, tmp[1], 81 );
+            }
+            else if ( strncmp(tmp[0], "atom_info", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->atom_format += ival * 4;
+            }
+            else if ( strncmp(tmp[0], "atom_velocities", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->atom_format += ival * 2;
+            }
+            else if ( strncmp(tmp[0], "atom_forces", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->atom_format += ival * 1;
+            }
+            else if ( strncmp(tmp[0], "bond_info", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->bond_info = ival;
+            }
+            else if ( strncmp(tmp[0], "angle_info", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                out_control->angle_info = ival;
+            }
+            else if ( strncmp(tmp[0], "test_forces", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+            }
+            else if ( strncmp(tmp[0], "molec_anal", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->molec_anal = ival;
+            }
+            else if ( strncmp(tmp[0], "freq_molec_anal", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->freq_molec_anal = ival;
+            }
+            else if ( strncmp(tmp[0], "bond_graph_cutoff", MAX_LINE) == 0 )
+            {
+                val = atof(tmp[1]);
+                control->bg_cut = val;
+            }
+            else if ( strncmp(tmp[0], "ignore", MAX_LINE) == 0 )
+            {
+                control->num_ignored = atoi(tmp[1]);
+                for ( i = 0; i < control->num_ignored; ++i )
+                    control->ignore[atoi(tmp[i + 2])] = 1;
+            }
+            else if ( strncmp(tmp[0], "dipole_anal", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->dipole_anal = ival;
+            }
+            else if ( strncmp(tmp[0], "freq_dipole_anal", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->freq_dipole_anal = ival;
+            }
+            else if ( strncmp(tmp[0], "diffusion_coef", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->diffusion_coef = ival;
+            }
+            else if ( strncmp(tmp[0], "freq_diffusion_coef", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->freq_diffusion_coef = ival;
+            }
+            else if ( strncmp(tmp[0], "restrict_type", MAX_LINE) == 0 )
+            {
+                ival = atoi(tmp[1]);
+                control->restrict_type = ival;
+            }
+            else
+            {
+                fprintf( stderr, "WARNING: unknown parameter %s\n", tmp[0] );
+                exit( UNKNOWN_OPTION );
             }
-        }
-        else if ( strncmp(tmp[0], "traj_title", MAX_LINE) == 0 )
-        {
-            strncpy( out_control->traj_title, tmp[1], 81 );
-        }
-        else if ( strncmp(tmp[0], "atom_info", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->atom_format += ival * 4;
-        }
-        else if ( strncmp(tmp[0], "atom_velocities", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->atom_format += ival * 2;
-        }
-        else if ( strncmp(tmp[0], "atom_forces", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->atom_format += ival * 1;
-        }
-        else if ( strncmp(tmp[0], "bond_info", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->bond_info = ival;
-        }
-        else if ( strncmp(tmp[0], "angle_info", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            out_control->angle_info = ival;
-        }
-        else if ( strncmp(tmp[0], "test_forces", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-        }
-        else if ( strncmp(tmp[0], "molec_anal", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->molec_anal = ival;
-        }
-        else if ( strncmp(tmp[0], "freq_molec_anal", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->freq_molec_anal = ival;
-        }
-        else if ( strncmp(tmp[0], "bond_graph_cutoff", MAX_LINE) == 0 )
-        {
-            val = atof(tmp[1]);
-            control->bg_cut = val;
-        }
-        else if ( strncmp(tmp[0], "ignore", MAX_LINE) == 0 )
-        {
-            control->num_ignored = atoi(tmp[1]);
-            for ( i = 0; i < control->num_ignored; ++i )
-                control->ignore[atoi(tmp[i + 2])] = 1;
-        }
-        else if ( strncmp(tmp[0], "dipole_anal", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->dipole_anal = ival;
-        }
-        else if ( strncmp(tmp[0], "freq_dipole_anal", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->freq_dipole_anal = ival;
-        }
-        else if ( strncmp(tmp[0], "diffusion_coef", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->diffusion_coef = ival;
-        }
-        else if ( strncmp(tmp[0], "freq_diffusion_coef", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->freq_diffusion_coef = ival;
-        }
-        else if ( strncmp(tmp[0], "restrict_type", MAX_LINE) == 0 )
-        {
-            ival = atoi(tmp[1]);
-            control->restrict_type = ival;
-        }
-        else
-        {
-            fprintf( stderr, "WARNING: unknown parameter %s\n", tmp[0] );
-            exit( UNKNOWN_OPTION );
         }
     }
 
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 787837988b059dd97043fda691404c8bb77ae4a9..5b2f0a7d2a7ca29d85a74b86b1fda3db248a15d5 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -136,7 +136,7 @@ static void Compute_NonBonded_Forces( reax_system *system, control_params *contr
 #endif
 
     t_start = Get_Time( );
-    Compute_Charges( system, control, data, workspace, lists[FAR_NBRS], out_control );
+    Compute_Charges( system, control, data, workspace, out_control );
     t_elapsed = Get_Timing_Info( t_start );
     data->timing.cm += t_elapsed;
 
@@ -185,7 +185,8 @@ static void Compute_Total_Force( reax_system *system, control_params *control,
             {
                 if ( i < bonds->select.bond_list[pj].nbr )
                 {
-                    if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT)
+                    if ( control->ensemble == NVE || control->ensemble == nhNVT
+                            || control->ensemble == bNVT )
                     {
                         Add_dBond_to_Forces( i, pj, system, data, workspace, lists );
                     }
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index c46e9bc6dfc3938b606cfe84f29cb84165161fd8..2444c43f2ee0bbafba924c34a75d9af6a1e00f3e 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -494,7 +494,8 @@ void Four_Body_Interactions( reax_system *system, control_params *control,
 #endif
                                         bo_kl->Cdbo += (CEtors6 + CEconj3);
 
-                                        if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT )
+                                        if ( control->ensemble == NVE || control->ensemble == nhNVT
+                                                || control->ensemble == bNVT )
                                         {
                                             /* dcos_theta_ijk */
                                             rvec_ScaledAdd( *f_i,
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index cc0ce302ecb46c17b71cfb87810ecbae5ff2d6c3..6670dfc06f23b20b3f2bbf718d26d288dcef611c 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -149,7 +149,7 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
         break;
 
 
-    case NVT:
+    case nhNVT:
         data->N_f = 3 * system->N + 1;
         //control->Tau_T = 100 * data->N_f * K_B * control->T_final;
         if ( !control->restart || (control->restart && control->random_vel) )
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 4e948406e770ec744970955ca1e6d37f4e4306f4..af2d87bd2676d8dd705c5db5be254e846e5b43e3 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -58,11 +58,11 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
         rvec_ScaledSum( dx, dt, system->atoms[i].v,
-                0.5 * dt_sqr * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt_sqr * F_CONV * inv_m, system->atoms[i].f );
         Inc_on_T3( system->atoms[i].x, dx, &( system->box ) );
 
         rvec_ScaledAdd( system->atoms[i].v,
-                0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
 
 #if defined(DEBUG_FOCUS)
@@ -83,7 +83,7 @@ void Velocity_Verlet_NVE(reax_system *system, control_params *control,
     {
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
         rvec_ScaledAdd( system->atoms[i].v,
-                0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
 
 #if defined(DEBUG_FOCUS)
@@ -118,7 +118,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
         rvec_ScaledSum( dx, dt - 0.5 * dt_sqr * therm->v_xi, system->atoms[i].v,
-                0.5 * dt_sqr * inv_m * -F_CONV, system->atoms[i].f );
+                -0.5 * dt_sqr * inv_m * F_CONV, system->atoms[i].f );
 
         Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
 
@@ -150,14 +150,14 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
         rvec_Scale( workspace->v_const[i],
                     1.0 - 0.5 * dt * therm->v_xi, system->atoms[i].v );
         rvec_ScaledAdd( workspace->v_const[i],
-                        0.5 * dt * inv_m * -F_CONV, workspace->f_old[i] );
+                -0.5 * dt * inv_m * F_CONV, workspace->f_old[i] );
         rvec_ScaledAdd( workspace->v_const[i],
-                        0.5 * dt * inv_m * -F_CONV, system->atoms[i].f );
+                -0.5 * dt * inv_m * F_CONV, system->atoms[i].f );
 #if defined(DEBUG)
         fprintf( stderr, "atom%d: inv_m=%f, C1=%f, C2=%f, v_const=%f %f %f\n",
-                 i, inv_m, 1.0 - 0.5 * dt * therm->v_xi,
-                 0.5 * dt * inv_m * -F_CONV, workspace->v_const[i][0],
-                 workspace->v_const[i][1], workspace->v_const[i][2] );
+                i, inv_m, 1.0 - 0.5 * dt * therm->v_xi,
+                -0.5 * dt * inv_m * F_CONV, workspace->v_const[i][0],
+                workspace->v_const[i][1], workspace->v_const[i][2] );
 #endif
     }
 
@@ -219,30 +219,17 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
     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",
-    //control->Tau_T, control->Tau_P, dt / control->Tau_T, dt / control->Tau_P );
-    fprintf( stderr, "step %d: ", data->step );
-#endif
-
     /* velocity verlet, 1st part */
     for ( i = 0; i < system->N; i++ )
     {
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
         /* Compute x(t + dt) */
         rvec_ScaledSum( dx, dt, system->atoms[i].v,
-                        0.5 * -F_CONV * inv_m * SQR(dt), system->atoms[i].f );
+                -0.5 * F_CONV * inv_m * SQR(dt), system->atoms[i].f );
         Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
         /* Compute v(t + dt/2) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * -F_CONV * inv_m * dt, system->atoms[i].f );
-        /*fprintf( stderr, "%6d   %15.8f %15.8f %15.8f   %15.8f %15.8f %15.8f\n",
-          workspace->orig_id[i],
-          system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[0],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[1],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[2] ); */
+                -0.5 * F_CONV * inv_m * dt, system->atoms[i].f );
     }
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "verlet1 - " );
@@ -265,13 +252,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
         /* Compute v(t + dt) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
-        /* fprintf( stderr, "%6d   %15f %15f %15f   %15.8f %15.8f %15.8f\n",
-           workspace->orig_id[i],
-           system->atoms[i].v[0], system->atoms[i].v[1], system->atoms[i].v[2],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[0],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[1],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[2] );*/
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
     Compute_Kinetic_Energy( system, data );
     Compute_Pressure_Isotropic( system, control, data, out_control );
@@ -350,17 +331,11 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
         /* Compute x(t + dt) */
         rvec_ScaledSum( dx, dt, system->atoms[i].v,
-                        0.5 * -F_CONV * inv_m * SQR(dt), system->atoms[i].f );
+                -0.5 * F_CONV * inv_m * SQR(dt), system->atoms[i].f );
         Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
         /* Compute v(t + dt/2) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * -F_CONV * inv_m * dt, system->atoms[i].f );
-        /*fprintf( stderr, "%6d   %15.8f %15.8f %15.8f   %15.8f %15.8f %15.8f\n",
-          workspace->orig_id[i],
-          system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[0],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[1],
-          0.5 * SQR(dt) * -F_CONV * inv_m * system->atoms[i].f[2] ); */
+                -0.5 * F_CONV * inv_m * dt, system->atoms[i].f );
     }
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "verlet1 - " );
@@ -383,19 +358,12 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
         /* Compute v(t + dt) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
-        /* fprintf( stderr, "%6d   %15f %15f %15f   %15.8f %15.8f %15.8f\n",
-           workspace->orig_id[i],
-           system->atoms[i].v[0], system->atoms[i].v[1], system->atoms[i].v[2],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[0],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[1],
-           0.5 * dt * -F_CONV * inv_m * system->atoms[i].f[2] );*/
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
+
     Compute_Kinetic_Energy( system, data );
+
     Compute_Pressure_Isotropic( system, control, data, out_control );
-#if defined(DEBUG_FOCUS)
-    fprintf( stderr, "verlet2 - " );
-#endif
 
     /* pressure scaler */
     for ( d = 0; d < 3; ++d )
@@ -460,26 +428,26 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
     real dt_sqr = SQR(dt);
     rvec dx;
 
-    for (i = 0; i < system->N; i++)
+    for ( i = 0; i < system->N; i++ )
     {
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
-        // Compute x(t + dt)
+        /* Compute x(t + dt) */
         rvec_ScaledSum( dx, dt, system->atoms[i].v,
-                        0.5 * dt_sqr * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt_sqr * F_CONV * inv_m, system->atoms[i].f );
         Inc_on_T3_Gen( system->atoms[i].x, dx, &(system->box) );
 
-        // Compute v(t + dt/2)
+        /* Compute v(t + dt/2) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        -0.5 * dt * data->therm.xi, system->atoms[i].v );
+                -0.5 * dt * data->therm.xi, system->atoms[i].v );
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
 
     // Compute zeta(t + dt/2), E_Kininetic(t + dt/2)
     // IMPORTANT: What will be the initial value of zeta? and what is g?
     data->therm.xi += 0.5 * dt * control->Tau_T  *
-                      ( 2.0 * data->E_Kin - data->N_f * K_B * control->T );
+        ( 2.0 * data->E_Kin - data->N_f * K_B * control->T );
 
     Reset( system, control, data, workspace );
     fprintf(out_control->log, "reset-");
@@ -490,7 +458,7 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
     fprintf(out_control->log, "nbrs-");
     fflush( out_control->log );
 
-    /* Compute_Charges( system, control, workspace, lists[FAR_NBRS], out_control );
+    /* Compute_Charges( system, control, workspace, out_control );
        fprintf(out_control->log,"qeq-"); fflush( out_control->log ); */
 
     Compute_Forces( system, control, data, workspace, lists, out_control,
@@ -504,14 +472,14 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
     {
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
-        // compute v(t + dt)
+        /* compute v(t + dt) */
         rvec_ScaledAdd( system->atoms[i].v,
-                        -0.5 * dt * data->therm.xi, system->atoms[i].v );
+                -0.5 * dt * data->therm.xi, system->atoms[i].v );
         rvec_ScaledAdd( system->atoms[i].v,
-                        0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
     }
 
-    // Compute zeta(t + dt)
+    /* Compute zeta(t + dt) */
     data->therm.xi += 0.5 * dt * control->Tau_T  * ( 2.0 * data->E_Kin -
                       data->N_f * K_B * control->T );
 
@@ -567,12 +535,12 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
     {
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
-        // Compute x(t + dt)
-        rvec_ScaledSum( workspace->a[i], -F_CONV * inv_m, system->atoms[i].f,
-                        -( (2.0 + 3.0 / data->N_f) * iso_bar->v_eps + therm->v_xi ),
-                        system->atoms[i].v );
+        /* Compute x(t + dt) */
+        rvec_ScaledSum( workspace->a[i], -1.0 * F_CONV * inv_m, system->atoms[i].f,
+                -1.0 * ( (2.0 + 3.0 / data->N_f) * iso_bar->v_eps + therm->v_xi ),
+                system->atoms[i].v );
         rvec_ScaledSum( dx, dt, system->atoms[i].v,
-                        0.5 * dt_sqr, workspace->a[i] );
+                0.5 * dt_sqr, workspace->a[i] );
         Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
         rvec_Scale( system->atoms[i].x, exp_deps, system->atoms[i].x );
     }
@@ -593,7 +561,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
     fprintf(out_control->log, "nbrs-");
     fflush( out_control->log );
 
-    /* Compute_Charges( system, control, workspace, lists[FAR_NBRS], out_control );
+    /* Compute_Charges( system, control, workspace, out_control );
        fprintf(out_control->log,"qeq-"); fflush( out_control->log ); */
 
     Compute_Forces( system, control, data, workspace, lists,
@@ -609,15 +577,14 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
         inv_m = 1.0 / system->reaxprm.sbp[system->atoms[i].type].mass;
 
         rvec_ScaledSum( dv, 0.5 * dt, workspace->a[i],
-                        0.5 * dt * -F_CONV * inv_m, system->atoms[i].f );
+                -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
         rvec_Add( dv, system->atoms[i].v );
         rvec_Scale( workspace->v_const[i], exp_deps, dv );
 
-        P_int_const += ( -F_CONV *
-                         rvec_Dot( system->atoms[i].f, system->atoms[i].x ) );
+        P_int_const += ( -1.0 * F_CONV * rvec_Dot( system->atoms[i].f, system->atoms[i].x ) );
 
-        E_kin += (0.5 * system->reaxprm.sbp[system->atoms[i].type].mass *
-                  rvec_Dot( system->atoms[i].v, system->atoms[i].v ) );
+        E_kin += (0.5 * system->reaxprm.sbp[system->atoms[i].type].mass
+                * rvec_Dot( system->atoms[i].v, system->atoms[i].v ) );
     }
 
     // Compute initial p_int
@@ -716,15 +683,15 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
         atom = &(system->atoms[i]);
         inv_m = 1.0 / system->reaxprm.sbp[atom->type].mass;
         /* Compute x(t + dt) */
-        rvec_ScaledSum( dx, dt, atom->v, 0.5 * -F_CONV * inv_m * SQR(dt), atom->f );
+        rvec_ScaledSum( dx, dt, atom->v, -0.5 * F_CONV * inv_m * SQR(dt), atom->f );
 
         /* bNVT fix - Metin's suggestion */
         /* ORIGINAL CHANGE -- CHECK THE branch serial-bnvt for the fix */
         //rvec_Add( atom->x, dx );
-        Inc_on_T3( atom->x, dx, &( system->box ) );
+        Inc_on_T3( atom->x, dx, &system->box );
 
         /* Compute v(t + dt/2) */
-        rvec_ScaledAdd( atom->v, 0.5 * -F_CONV * inv_m * dt, atom->f );
+        rvec_ScaledAdd( atom->v, -0.5 * F_CONV * inv_m * dt, atom->f );
     }
 
 #if defined(DEBUG_FOCUS)
@@ -748,7 +715,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
         atom = &(system->atoms[i]);
         inv_m = 1.0 / system->reaxprm.sbp[atom->type].mass;
         /* Compute v(t + dt) */
-        rvec_ScaledAdd( atom->v, 0.5 * dt * -F_CONV * inv_m, atom->f );
+        rvec_ScaledAdd( atom->v, -0.5 * dt * F_CONV * inv_m, atom->f );
     }
 #if defined(DEBUG_FOCUS)
     fprintf(stderr, "step%d: verlet2 done\n", data->step);
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 26870dee539af412a5529c04e530fb5fdd97ff04..754e0c93b030ef976f41b350f08dc3e9e34c9de6 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -129,18 +129,19 @@
 #define LOOSE_ZONE     (0.75)
 
 
-/* config params */
+/* ensemble type */
 enum ensemble
 {
     NVE = 0,
-    NVT = 1,
-    NPT = 2,
+    bNVT = 1,
+    nhNVT = 2,
     sNPT = 3,
     iNPT = 4,
-    ensNR = 5,
-    bNVT = 6,
+    NPT = 5,
+    ens_N = 6,
 };
 
+/* interaction list type */
 enum interaction_list_offets
 {
     FAR_NBRS = 0,
@@ -154,6 +155,7 @@ enum interaction_list_offets
     LIST_N = 8,
 };
 
+/* interaction type */
 enum interaction_type
 {
     TYP_VOID = 0,
@@ -167,6 +169,7 @@ enum interaction_type
     TYP_N = 8,
 };
 
+/* error codes for simulation termination */
 enum errors
 {
     FILE_NOT_FOUND = -10,
@@ -190,13 +193,15 @@ enum molecular_analysis_type
     NUM_ANALYSIS = 3,
 };
 
-enum restart_format
+/* restart file format */
+enum restart_formats
 {
     WRITE_ASCII = 0,
     WRITE_BINARY = 1,
     RF_N = 2,
 };
 
+/* geometry file format */
 enum geo_formats
 {
     CUSTOM = 0,
@@ -243,6 +248,15 @@ enum pre_app
 };
 
 
+/* atom types as pertains to hydrogen bonding */
+enum hydrogen_bonding_atom_types
+{
+    NON_H_BONDING_ATOM = 0,
+    H_ATOM = 1,
+    H_BONDING_ATOM = 2,
+};
+
+
 typedef double real;
 typedef real rvec[3];
 typedef int ivec[3];
@@ -603,22 +617,52 @@ typedef struct
     int freq_diffusion_coef;
     int restrict_type;
 
+    /* method for computing atomic charges */
     unsigned int charge_method;
+    /* frequency (in terms of simulation time steps) at which to
+     * re-compute atomic charge distribution */
+    int charge_freq;
+    /* iterative linear solver type */
     unsigned int cm_solver_type;
+    /* system net charge */
     real cm_q_net;
+    /* max. iterations for linear solver */
     unsigned int cm_solver_max_iters;
+    /* max. iterations before restarting in specific solvers, e.g., GMRES(k) */
     unsigned int cm_solver_restart;
+    /* error tolerance of solution produced by charge distribution
+     * sparse iterative linear solver */
     real cm_solver_q_err;
+    /* ratio used in computing sparser charge matrix,
+     * between 0.0 and 1.0 */
     real cm_domain_sparsity;
+    /* TRUE if enabled, FALSE otherwise */
     unsigned int cm_domain_sparsify_enabled;
+    /* order of spline extrapolation used for computing initial guess
+     * to linear solver */
     unsigned int cm_init_guess_extrap1;
+    /* order of spline extrapolation used for computing initial guess
+     * to linear solver */
     unsigned int cm_init_guess_extrap2;
+    /* preconditioner type for linear solver */
     unsigned int cm_solver_pre_comp_type;
+    /* frequency (in terms of simulation time steps) at which to recompute
+     * incomplete factorizations */
     unsigned int cm_solver_pre_comp_refactor;
+    /* drop tolerance of incomplete factorization schemes (ILUT, ICHOLT, etc.)
+     * used for preconditioning the iterative linear solver used in charge distribution */
     real cm_solver_pre_comp_droptol;
+    /* num. of sweeps for computing preconditioner factors
+     * in fine-grained iterative methods (FG-ICHOL, FG-ILU) */
     unsigned int cm_solver_pre_comp_sweeps;
+    /* relative num. of non-zeros to charge matrix used to
+     * compute the sparse approximate inverse preconditioner,
+     * between 0.0 and 1.0 */
     real cm_solver_pre_comp_sai_thres;
+    /* preconditioner application type */
     unsigned int cm_solver_pre_app_type;
+    /* num. of iterations used to apply preconditioner via
+     * Jacobi relaxation scheme (truncated Neumann series) */
     unsigned int cm_solver_pre_app_jacobi_iters;
 
     int molec_anal;
@@ -671,24 +715,42 @@ typedef struct
 
 typedef struct
 {
+    /* start time of event */
     real start;
+    /* end time of event */
     real end;
+    /* total elapsed time of event */
     real elapsed;
-
+    /* total simulation time */
     real total;
+    /* neighbor list generation time */
     real nbrs;
+    /* force initialization time */
     real init_forces;
+    /* bonded force calculation time */
     real bonded;
+    /* non-bonded force calculation time */
     real nonb;
+    /* atomic charge distribution calculation time */
     real cm;
+    /**/
     real cm_sort_mat_rows;
+    /**/
     real cm_solver_pre_comp;
+    /**/
     real cm_solver_pre_app;
+    /* num. of steps in iterative linear solver for charge distribution */
     int cm_solver_iters;
+    /**/
     real cm_solver_spmv;
+    /**/
     real cm_solver_vector_ops;
+    /**/
     real cm_solver_orthog;
+    /**/
     real cm_solver_tri_solve;
+    /* num. of retries in main sim. loop */
+    int num_retries;
 } reax_timing;
 
 
@@ -953,17 +1015,30 @@ typedef struct
     rvec *dDeltap_self;
 
     /* charge method storage */
+    /* charge matrix */
     sparse_matrix *H;
+    /* charge matrix (full) */
     sparse_matrix *H_full;
+    /* sparser charge matrix */
     sparse_matrix *H_sp;
+    /* permuted charge matrix (graph coloring) */
     sparse_matrix *H_p;
+    /* sparsity pattern of charge matrix, used in
+     * computing a sparse approximate inverse preconditioner */
     sparse_matrix *H_spar_patt;
+    /* sparsity pattern of charge matrix (full), used in
+     * computing a sparse approximate inverse preconditioner */
     sparse_matrix *H_spar_patt_full;
+    /* sparse approximate inverse preconditioner */
     sparse_matrix *H_app_inv;
+    /* incomplete Cholesky or LU preconditioner */
     sparse_matrix *L;
+    /* incomplete Cholesky or LU preconditioner */
     sparse_matrix *U;
-    real *droptol;
+    /* Jacobi preconditioner */
     real *Hdia_inv;
+    /* row drop tolerences for incomplete Cholesky preconditioner */
+    real *droptol;
     real *b;
     real *b_s;
     real *b_t;
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 5837e3ea5e36fbd2dd7d20e791fe468f5800bd87..8e8fc825f09202e6a179275bf358f34efe618f96 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -477,7 +477,8 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
                                     }
 
 
-                                    if ( control->ensemble == NVE || control->ensemble == NVT  || control->ensemble == bNVT)
+                                    if ( control->ensemble == NVE || control->ensemble == nhNVT
+                                            || control->ensemble == bNVT )
                                     {
                                         rvec_ScaledAdd( *f_i, CEval8, p_ijk->dcos_di );
                                         rvec_ScaledAdd( *f_j, CEval8, p_ijk->dcos_dj );
@@ -820,7 +821,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
 #endif
                             bo_ij->Cdbo += CEhb1;
 
-                            if ( control->ensemble == NVE || control->ensemble == NVT  || control->ensemble == bNVT)
+                            if ( control->ensemble == NVE || control->ensemble == nhNVT
+                                    || control->ensemble == bNVT )
                             {
                                 /* dcos terms */
                                 rvec_ScaledAdd( *f_i, +CEhb2, dcos_theta_di );
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index 9a04bb71651d64aeafb43e2a6afa1ecafdea6b1c..8fb702c3eef0e19de3402e921883c6232f92c0db 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -220,11 +220,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
             {
                 if ( far_nbrs->select.far_nbr_list[pj].d <= control->r_cut )
                 {
-                    nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
+                    nbr_pj = &far_nbrs->select.far_nbr_list[pj];
                     j = nbr_pj->nbr;
                     r_ij = nbr_pj->d;
-                    twbp = &(system->reaxprm.tbp[ system->atoms[i].type ]
-                             [ system->atoms[j].type ]);
+                    twbp = &system->reaxprm.tbp[ system->atoms[i].type ]
+                             [ system->atoms[j].type ];
                     self_coef = (i == j) ? 0.5 : 1.0; // for supporting small boxes!
 
                     /* Calculate Taper and its derivative */
@@ -301,7 +301,8 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
                     CEclmb = self_coef * C_ele * system->atoms[i].q * system->atoms[j].q *
                              ( dTap -  Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
 
-                    if ( control->ensemble == NVE || control->ensemble == NVT || control->ensemble == bNVT )
+                    if ( control->ensemble == NVE || control->ensemble == nhNVT
+                            || control->ensemble == bNVT )
                     {
 #ifndef _OPENMP
                         rvec_ScaledAdd( system->atoms[i].f,
@@ -606,7 +607,8 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
                              t->CEclmb[r].a;
                     CEclmb *= self_coef * system->atoms[i].q * system->atoms[j].q;
 
-                    if ( control->ensemble == NVE || control->ensemble == NVT  || control->ensemble == bNVT)
+                    if ( control->ensemble == NVE || control->ensemble == nhNVT
+                            || control->ensemble == bNVT )
                     {
 #ifndef _OPENMP
                         rvec_ScaledAdd( system->atoms[i].f, -(CEvd + CEclmb), nbr_pj->dvec );