diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 3b760d682c96a22f223f316f305041bea0f80ad7..87312167ca54ce2aa9258dcb0f1fbea60d2c6492 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -949,7 +949,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, fillin );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, fillin );
             }
-            else if ( workspace->L.m < fillin || realloc == TRUE )
+            else if ( workspace->L.m < fillin || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -970,7 +971,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -991,7 +993,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -1086,7 +1089,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -1113,7 +1117,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -1210,7 +1215,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
@@ -1239,7 +1245,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
                 Allocate_Matrix( &workspace->L, Hptr->n, Hptr->n_max, Hptr->m );
                 Allocate_Matrix( &workspace->U, Hptr->n, Hptr->n_max, Hptr->m );
             }
-            else if ( workspace->L.m < Hptr->m || realloc == TRUE )
+            else if ( workspace->L.m < Hptr->m || workspace->L.n_max < system->N_cm_max
+                    || realloc == TRUE )
             {
                 Deallocate_Matrix( &workspace->L );
                 Deallocate_Matrix( &workspace->U );
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index f57684d7381bcf835b6e8badd5d4408bb0ff7f96..736e543990a2469d2529a0019685c6fa84912e19 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -374,23 +374,35 @@ static void Init_Workspace( reax_system *system, control_params *control,
     {
         case QEQ_CM:
             system->N_cm = system->N;
-            system->N_cm_max = system->N_max;
+            if ( realloc == TRUE || system->N_cm > system->N_cm_max )
+            {
+                system->N_cm_max = system->N_max;
+            }
             break;
         case EE_CM:
             if ( system->num_molec_charge_constraints == 0 )
             {
                 system->N_cm = system->N + 1;
-                system->N_cm_max = system->N_max + 1;
+                if ( realloc == TRUE || system->N_cm > system->N_cm_max )
+                {
+                    system->N_cm_max = system->N_max + 1;
+                }
             }
             else
             {
                 system->N_cm = system->N + system->num_molec_charge_constraints;
-                system->N_cm_max = system->N_max + system->num_molec_charge_constraints;
+                if ( realloc == TRUE || system->N_cm > system->N_cm_max )
+                {
+                    system->N_cm_max = system->N_max + system->num_molec_charge_constraints;
+                }
             }
             break;
         case ACKS2_CM:
             system->N_cm = 2 * system->N + 2;
-            system->N_cm_max = 2 * system->N_max + 2;
+            if ( realloc == TRUE || system->N_cm > system->N_cm_max )
+            {
+                system->N_cm_max = 2 * system->N_max + 2;
+            }
             break;
         default:
             fprintf( stderr, "[ERROR] Unknown charge method type. Terminating...\n" );
@@ -845,7 +857,8 @@ static void Init_Lists( reax_system *system, control_params *control,
     {
         Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, max_nnz );
     }
-    else if ( realloc == TRUE || workspace->H.m < max_nnz )
+    else if ( realloc == TRUE || workspace->H.m < max_nnz
+            || workspace->H.n_max < system->N_cm_max )
     {
         if ( workspace->H.allocated == TRUE )
         {
@@ -866,7 +879,8 @@ static void Init_Lists( reax_system *system, control_params *control,
          *   (non-bonded, hydrogen, 3body, etc.) */
         Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max, max_nnz );
     }
-    else if ( realloc == TRUE || workspace->H_sp.m < max_nnz )
+    else if ( realloc == TRUE || workspace->H_sp.m < max_nnz
+            || workspace->H.n_max < system->N_cm_max )
     {
         if ( workspace->H_sp.allocated == TRUE )
         {
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index c6abc8347b38295c43c3f086c2de94310764ec18..98b71d9af175b31ebe97b988dffa479150bcd562 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -171,9 +171,13 @@ static void compute_full_sparse_matrix( const sparse_matrix * const A,
     {
         Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
     }
-    else if ( A_full->m < 2 * A->m - A->n || realloc == TRUE )
+    else if ( A_full->m < 2 * A->m - A->n || A_full->n_max < A->n_max
+            || realloc == TRUE )
     {
-        Deallocate_Matrix( A_full );
+        if ( A_full->allocated == TRUE )
+        {
+            Deallocate_Matrix( A_full );
+        }
         Allocate_Matrix( A_full, A->n, A->n_max, 2 * A->m - A->n );
     }
 
@@ -236,9 +240,13 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A,
     {
         Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
     }
-    else if ( A_spar_patt->m < A->m || realloc == TRUE )
+    else if ( A_spar_patt->m < A->m || A_spar_patt->n_max < A->n_max
+            || realloc == TRUE )
     {
-        Deallocate_Matrix( A_spar_patt );
+        if ( A_spar_patt->allocated == TRUE )
+        {
+            Deallocate_Matrix( A_spar_patt );
+        }
         Allocate_Matrix( A_spar_patt, A->n, A->n_max, A->m );
     }
 
@@ -355,9 +363,13 @@ void setup_sparse_approx_inverse( const sparse_matrix * const A,
         Allocate_Matrix( A_app_inv, A_spar_patt_full->n,
                 A_spar_patt_full->n_max, A_spar_patt_full->m );
     }
-    else if ( A_app_inv->m < A->m || realloc == TRUE )
+    else if ( A_app_inv->m < A->m || A_app_inv->n_max < A->n_max
+            || realloc == TRUE )
     {
-        Deallocate_Matrix( A_app_inv );
+        if ( A_app_inv->allocated == TRUE )
+        {
+            Deallocate_Matrix( A_app_inv );
+        }
 
         /* A_app_inv has the same sparsity pattern
          * as A_spar_patt_full (omit non-zero values) */
@@ -2634,9 +2646,12 @@ void setup_graph_coloring( const control_params * const control,
     {
         Allocate_Matrix( H_p, H->n, H->n_max, H->m );
     }
-    else if ( H_p->m < H->m || realloc == TRUE )
+    else if ( H_p->m < H->m || H_p->n_max < H->n_max || realloc == TRUE )
     {
-        Deallocate_Matrix( H_p );
+        if ( H_p->allocated == TRUE )
+        {
+            Deallocate_Matrix( H_p );
+        }
         Allocate_Matrix( H_p, H->n, H->n_max, H->m );
     }