From 4c3eb4a843982cdefa65003e36b1d5c47ce9bd7a Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Fri, 2 Feb 2018 11:29:22 -0500
Subject: [PATCH] sPuReMD: cleanup memory allocation routines.

---
 sPuReMD/src/charges.c   |  31 ++-----
 sPuReMD/src/ffield.c    |  63 +++++++++-----
 sPuReMD/src/forces.c    |   7 +-
 sPuReMD/src/geo_tools.c |  37 +++++----
 sPuReMD/src/grid.c      | 101 +++++++++++++++--------
 sPuReMD/src/init_md.c   | 107 ++++++++++++++++--------
 sPuReMD/src/lin_alg.c   | 177 ++++++++++++++++------------------------
 sPuReMD/src/lookup.c    |  57 ++++++++-----
 sPuReMD/src/spuremd.c   |   3 +-
 sPuReMD/src/tool_box.c  |  30 ++-----
 sPuReMD/src/tool_box.h  |   2 +-
 11 files changed, 333 insertions(+), 282 deletions(-)

diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index e75c2dc9..7a3900c0 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -259,13 +259,9 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
 //
 //    if ( ones == NULL )
 //    {
-//        if ( ( ones = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
-//            ( x = (real*) malloc( system->N * sizeof(real)) ) == NULL ||
-//            ( y = (real*) malloc( system->N * sizeof(real)) ) == NULL )
-//        {
-//            fprintf( stderr, "Not enough space for preconditioner computation. Terminating...\n" );
-//            exit( INSUFFICIENT_MEMORY );
-//        }
+//        ones = (real*) smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::ones" );
+//        x = (real*) smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::x" );
+//        y = (real*) smalloc( system->N * sizeof(real), "Compute_Preconditioner_EE::y" );
 //
 //        for ( i = 0; i < system->N; ++i )
 //        {
@@ -764,11 +760,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
         case DIAG_PC:
             if ( workspace->Hdia_inv == NULL )
             {
-                if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
-                {
-                    fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                workspace->Hdia_inv = (real *) scalloc( Hptr->n, sizeof( real ),
+                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
             }
             break;
 
@@ -916,11 +909,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
         case DIAG_PC:
             if ( workspace->Hdia_inv == NULL )
             {
-                if ( ( workspace->Hdia_inv = (real *) calloc( system->N_cm, sizeof( real ) ) ) == NULL )
-                {
-                    fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                workspace->Hdia_inv = (real *) scalloc( system->N_cm, sizeof( real ),
+                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
             }
             break;
 
@@ -1071,11 +1061,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
         case DIAG_PC:
             if ( workspace->Hdia_inv == NULL )
             {
-                if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
-                {
-                    fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                workspace->Hdia_inv = (real *) scalloc( Hptr->n, sizeof( real ),
+                        "Setup_Preconditioner_QEq::workspace->Hdiv_inv" );
             }
             break;
 
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 56d8a310..7aab5c52 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -33,11 +33,14 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
     int i, j, k, l, m, n, o, p, cnt;
     real val;
 
-    s = (char*) malloc( sizeof(char) * MAX_LINE );
-    tmp = (char**) malloc( sizeof(char*) * MAX_TOKENS );
+    s = (char*) smalloc( sizeof(char) * MAX_LINE,
+           "Read_Force_Field::s" );
+    tmp = (char**) smalloc( sizeof(char*) * MAX_TOKENS,
+           "Read_Force_Field::tmp" );
     for (i = 0; i < MAX_TOKENS; i++)
     {
-        tmp[i] = (char*) malloc( sizeof(char) * MAX_TOKEN_LEN );
+        tmp[i] = (char*) smalloc( sizeof(char) * MAX_TOKEN_LEN,
+               "Read_Force_Field::tmp[i]" );
     }
 
     /* reading first header comment */
@@ -57,7 +60,8 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
     }
 
     reax->gp.n_global = n;
-    reax->gp.l = (real*) malloc( sizeof(real) * n );
+    reax->gp.l = (real*) smalloc( sizeof(real) * n,
+           "Read_Force_Field::reax->gp-l" );
 
     /* see mytypes.h for mapping between l[i] and the lambdas used in ff */
     for (i = 0; i < n; i++)
@@ -85,48 +89,65 @@ char Read_Force_Field( FILE* fp, reax_interaction* reax )
 
     /* Allocating structures in reax_interaction */
     reax->sbp = (single_body_parameters*)
-                calloc( reax->num_atom_types, sizeof(single_body_parameters) );
+                scalloc( reax->num_atom_types, sizeof(single_body_parameters),
+                        "Read_Force_Field::reax->sbp" );
     reax->tbp = (two_body_parameters**)
-                calloc( reax->num_atom_types, sizeof(two_body_parameters*) );
+                scalloc( reax->num_atom_types, sizeof(two_body_parameters*),
+                        "Read_Force_Field::reax->tbp" );
     reax->thbp = (three_body_header***)
-                 calloc( reax->num_atom_types, sizeof(three_body_header**) );
+                 scalloc( reax->num_atom_types, sizeof(three_body_header**),
+                         "Read_Force_Field::reax->thbp" );
     reax->hbp = (hbond_parameters***)
-                calloc( reax->num_atom_types, sizeof(hbond_parameters**) );
+                scalloc( reax->num_atom_types, sizeof(hbond_parameters**),
+                        "Read_Force_Field::reax->hbp" );
     reax->fbp = (four_body_header****)
-                calloc( reax->num_atom_types, sizeof(four_body_header***) );
+                scalloc( reax->num_atom_types, sizeof(four_body_header***),
+                        "Read_Force_Field::reax->fbp" );
     tor_flag  = (char****)
-                calloc( reax->num_atom_types, sizeof(char***) );
+                scalloc( reax->num_atom_types, sizeof(char***),
+                        "Read_Force_Field::tor_flag" );
 
     for ( i = 0; i < reax->num_atom_types; i++ )
     {
         reax->tbp[i] = (two_body_parameters*)
-                       calloc( reax->num_atom_types, sizeof(two_body_parameters) );
+                       scalloc( reax->num_atom_types, sizeof(two_body_parameters),
+                               "Read_Force_Field::reax->tbp[i]" );
         reax->thbp[i] = (three_body_header**)
-                        calloc( reax->num_atom_types, sizeof(three_body_header*) );
+                        scalloc( reax->num_atom_types, sizeof(three_body_header*),
+                                "Read_Force_Field::reax->thbp[i]" );
         reax->hbp[i] = (hbond_parameters**)
-                       calloc( reax->num_atom_types, sizeof(hbond_parameters*) );
+                       scalloc( reax->num_atom_types, sizeof(hbond_parameters*),
+                               "Read_Force_Field::reax->hbp[i]" );
         reax->fbp[i] = (four_body_header***)
-                       calloc( reax->num_atom_types, sizeof(four_body_header**) );
+                       scalloc( reax->num_atom_types, sizeof(four_body_header**),
+                               "Read_Force_Field::reax->fbp[i]" );
         tor_flag[i]  = (char***)
-                       calloc( reax->num_atom_types, sizeof(char**) );
+                       scalloc( reax->num_atom_types, sizeof(char**),
+                               "Read_Force_Field::tor_flag[i]" );
 
         for ( j = 0; j < reax->num_atom_types; j++ )
         {
             reax->thbp[i][j] = (three_body_header*)
-                               calloc( reax->num_atom_types, sizeof(three_body_header) );
+                               scalloc( reax->num_atom_types, sizeof(three_body_header),
+                                       "Read_Force_Field::reax->thbp[i][j]" );
             reax->hbp[i][j] = (hbond_parameters*)
-                              calloc( reax->num_atom_types, sizeof(hbond_parameters) );
+                              scalloc( reax->num_atom_types, sizeof(hbond_parameters),
+                                      "Read_Force_Field::reax->hbp[i][j]" );
             reax->fbp[i][j] = (four_body_header**)
-                              calloc( reax->num_atom_types, sizeof(four_body_header*) );
+                              scalloc( reax->num_atom_types, sizeof(four_body_header*),
+                                      "Read_Force_Field::reax->fbp[i][j]" );
             tor_flag[i][j]  = (char**)
-                              calloc( reax->num_atom_types, sizeof(char*) );
+                              scalloc( reax->num_atom_types, sizeof(char*),
+                                      "Read_Force_Field::tor_flag[i][j]" );
 
             for (k = 0; k < reax->num_atom_types; k++)
             {
                 reax->fbp[i][j][k] = (four_body_header*)
-                                     calloc( reax->num_atom_types, sizeof(four_body_header) );
+                                     scalloc( reax->num_atom_types, sizeof(four_body_header),
+                                             "Read_Force_Field::reax->fbp[i][j][k]" );
                 tor_flag[i][j][k]  = (char*)
-                                     calloc( reax->num_atom_types, sizeof(char) );
+                                     scalloc( reax->num_atom_types, sizeof(char),
+                                             "Read_Force_Field::tor_flag[i][j][k]" );
             }
         }
     }
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 10e31fdc..4ad3069f 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -476,11 +476,8 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
             break;
 
         case ACKS2_CM:
-            if ( (X_diag = (real*) calloc(system->N, sizeof(real))) == NULL )
-            {
-                fprintf( stderr, "not enough memory for charge matrix. terminating.\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            X_diag = (real*) scalloc( system->N, sizeof(real),
+                    "Init_Charge_Matrix_Remaining_Entries::X_diag" );
 
             H->start[system->N] = *Htop;
             H_sp->start[system->N] = *H_sp_top;
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index fcde21be..69675645 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -247,12 +247,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
     }
 
     /* allocate memory for tokenizing pdb lines */
-    if ( Allocate_Tokenizer_Space( &s, &s1, &tmp ) == FAILURE )
-    {
-        fprintf( stderr, "Allocate_Tokenizer_Space: not enough memory!" );
-        fprintf( stderr, "terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Allocate_Tokenizer_Space( &s, &s1, &tmp );
 
     /* read box information */
     if ( Read_Box_Info( system, pdb, PDB ) == FAILURE )
@@ -567,12 +562,16 @@ char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
     }
 
     /* allocate memory for tokenizing biograf file lines */
-    line   = (char*)  malloc( sizeof(char)  * MAX_LINE );
-    backup = (char*)  malloc( sizeof(char)  * MAX_LINE );
-    tokens = (char**) malloc( sizeof(char*) * MAX_TOKENS );
+    line = (char*) smalloc( sizeof(char)  * MAX_LINE,
+           "Read_BGF::line" );
+    backup = (char*) smalloc( sizeof(char)  * MAX_LINE,
+           "Read_BGF::backup" );
+    tokens = (char**) smalloc( sizeof(char*) * MAX_TOKENS,
+           "Read_BGF::tokens" );
     for ( i = 0; i < MAX_TOKENS; i++ )
     {
-        tokens[i] = (char*) malloc( sizeof(char) * MAX_TOKEN_LEN );
+        tokens[i] = (char*) smalloc( sizeof(char) * MAX_TOKEN_LEN,
+               "Read_BGF::tokens[i]" );
     }
 
     /* count number of atoms in the pdb file */
@@ -599,20 +598,26 @@ char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
     fclose( bgf );
 
     /* memory allocations for atoms, atom maps, bond restrictions */
-    system->atoms = (reax_atom*) calloc( system->N, sizeof(reax_atom) );
+    system->atoms = (reax_atom*) scalloc( system->N, sizeof(reax_atom),
+            "Read_BGF::system->atoms" );
 
-    workspace->map_serials = (int*) calloc( MAX_ATOM_ID, sizeof(int) );
+    workspace->map_serials = (int*) scalloc( MAX_ATOM_ID, sizeof(int),
+            "Read_BGF::workspace->map_serials" );
     for ( i = 0; i < MAX_ATOM_ID; ++i )
     {
         workspace->map_serials[i] = -1;
     }
 
-    workspace->orig_id = (int*) calloc( system->N, sizeof(int) );
-    workspace->restricted  = (int*) calloc( system->N, sizeof(int) );
-    workspace->restricted_list = (int**) calloc( system->N, sizeof(int*) );
+    workspace->orig_id = (int*) scalloc( system->N, sizeof(int),
+            "Read_BGF::workspace->orig_id" );
+    workspace->restricted  = (int*) scalloc( system->N, sizeof(int),
+            "Read_BGF::workspace->restricted" );
+    workspace->restricted_list = (int**) scalloc( system->N, sizeof(int*),
+            "Read_BGF::workspace->restricted_list" );
     for ( i = 0; i < system->N; ++i )
     {
-        workspace->restricted_list[i] = (int*) calloc( MAX_RESTRICT, sizeof(int) );
+        workspace->restricted_list[i] = (int*) scalloc( MAX_RESTRICT, sizeof(int),
+                "Read_BGF::workspace->restricted_list[i]" );
     }
 
     /* start reading and processing bgf file */
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 3b7a45ce..9c04c3f2 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -76,33 +76,54 @@ void Allocate_Space_for_Grid( reax_system *system )
         * (2 * g->spread[1] + 1) * (2 * g->spread[2] + 1) + 3;
 
     /* allocate space for the new grid */
-    g->atoms = (int****) calloc( g->ncell[0], sizeof( int*** ));
-    g->top = (int***) calloc( g->ncell[0], sizeof( int** ));
-    g->mark = (int***) calloc( g->ncell[0], sizeof( int** ));
-    g->start = (int***) calloc( g->ncell[0], sizeof( int** ));
-    g->end = (int***) calloc( g->ncell[0], sizeof( int** ));
-    g->nbrs = (ivec****) calloc( g->ncell[0], sizeof( ivec*** ));
-    g->nbrs_cp = (rvec****) calloc( g->ncell[0], sizeof( rvec*** ));
+    g->atoms = (int****) scalloc( g->ncell[0], sizeof( int*** ),
+            "Allocate_Space_for_Grid::g->atoms" );
+    g->top = (int***) scalloc( g->ncell[0], sizeof( int** ),
+            "Allocate_Space_for_Grid::g->top" );
+    g->mark = (int***) scalloc( g->ncell[0], sizeof( int** ),
+            "Allocate_Space_for_Grid::g->mark" );
+    g->start = (int***) scalloc( g->ncell[0], sizeof( int** ),
+            "Allocate_Space_for_Grid::g->start" );
+    g->end = (int***) scalloc( g->ncell[0], sizeof( int** ),
+            "Allocate_Space_for_Grid::g->end" );
+    g->nbrs = (ivec****) scalloc( g->ncell[0], sizeof( ivec*** ),
+            "Allocate_Space_for_Grid::g->nbrs" );
+    g->nbrs_cp = (rvec****) scalloc( g->ncell[0], sizeof( rvec*** ),
+            "Allocate_Space_for_Grid::g->nbrs_cp" );
 
     for ( i = 0; i < g->ncell[0]; i++ )
     {
-        g->atoms[i] = (int***) calloc( g->ncell[1], sizeof( int** ));
-        g->top [i] = (int**) calloc( g->ncell[1], sizeof( int* ));
-        g->mark[i] = (int**) calloc( g->ncell[1], sizeof( int* ));
-        g->start[i] = (int**) calloc( g->ncell[1], sizeof( int* ));
-        g->end[i] = (int**) calloc( g->ncell[1], sizeof( int* ));
-        g->nbrs[i] = (ivec***) calloc( g->ncell[1], sizeof( ivec** ));
-        g->nbrs_cp[i] = (rvec***) calloc( g->ncell[1], sizeof( rvec** ));
+        g->atoms[i] = (int***) scalloc( g->ncell[1], sizeof( int** ),
+                "Allocate_Space_for_Grid::g->atoms[i]" );
+        g->top [i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
+                "Allocate_Space_for_Grid::g->top[i]" );
+        g->mark[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
+                "Allocate_Space_for_Grid::g->mark[i]" );
+        g->start[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
+                "Allocate_Space_for_Grid::g->start[i]" );
+        g->end[i] = (int**) scalloc( g->ncell[1], sizeof( int* ),
+                "Allocate_Space_for_Grid::g->end[i]" );
+        g->nbrs[i] = (ivec***) scalloc( g->ncell[1], sizeof( ivec** ),
+                "Allocate_Space_for_Grid::g->nbrs[i]" );
+        g->nbrs_cp[i] = (rvec***) scalloc( g->ncell[1], sizeof( rvec** ),
+                "Allocate_Space_for_Grid::g->nbrs_cp[i]" );
 
         for ( j = 0; j < g->ncell[1]; j++ )
         {
-            g->atoms[i][j] = (int**) calloc( g->ncell[2], sizeof( int* ));
-            g->top[i][j] = (int*) calloc( g->ncell[2], sizeof( int ));
-            g->mark[i][j] = (int*) calloc( g->ncell[2], sizeof( int ));
-            g->start[i][j] = (int*) calloc( g->ncell[2], sizeof( int ));
-            g->end[i][j] = (int*) calloc( g->ncell[2], sizeof( int ));
-            g->nbrs[i][j] = (ivec**) calloc( g->ncell[2], sizeof( ivec* ));
-            g->nbrs_cp[i][j] = (rvec**) calloc( g->ncell[2], sizeof( rvec* ));
+            g->atoms[i][j] = (int**) scalloc( g->ncell[2], sizeof( int* ),
+                    "Allocate_Space_for_Grid::g->atoms[i][j]" );
+            g->top[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
+                    "Allocate_Space_for_Grid::g->top[i][j]" );
+            g->mark[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
+                    "Allocate_Space_for_Grid::g->mark[i][j]" );
+            g->start[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
+                    "Allocate_Space_for_Grid::g->start[i][j]" );
+            g->end[i][j] = (int*) scalloc( g->ncell[2], sizeof( int ),
+                    "Allocate_Space_for_Grid::g->end[i][j]" );
+            g->nbrs[i][j] = (ivec**) scalloc( g->ncell[2], sizeof( ivec* ),
+                    "Allocate_Space_for_Grid::g->nbrs[i][j]" );
+            g->nbrs_cp[i][j] = (rvec**) scalloc( g->ncell[2], sizeof( rvec* ),
+                    "Allocate_Space_for_Grid::g->nbrs_cp[i][j]" );
 
             for ( k = 0; k < g->ncell[2]; k++ )
             {
@@ -110,8 +131,10 @@ void Allocate_Space_for_Grid( reax_system *system )
                 g->mark[i][j][k] = 0;
                 g->start[i][j][k] = 0;
                 g->end[i][j][k] = 0;
-                g->nbrs[i][j][k] = (ivec*) malloc( g->max_nbrs * sizeof( ivec ) );
-                g->nbrs_cp[i][j][k] = (rvec*) malloc( g->max_nbrs * sizeof( rvec ) );
+                g->nbrs[i][j][k] = (ivec*) smalloc( g->max_nbrs * sizeof( ivec ),
+                       "Allocate_Space_for_Grid::g->nbrs[i][j]][k]" );
+                g->nbrs_cp[i][j][k] = (rvec*) smalloc( g->max_nbrs * sizeof( rvec ),
+                       "Allocate_Space_for_Grid::g->nbrs_cp[i][j]][k]" );
 
                 for ( l = 0; l < g->max_nbrs; ++l )
                 {
@@ -135,7 +158,8 @@ void Allocate_Space_for_Grid( reax_system *system )
         {
             for ( k = 0; k < g->ncell[2]; k++ )
             {
-                g->atoms[i][j][k] = (int*) malloc( g->max_atoms * sizeof( int ) );
+                g->atoms[i][j][k] = (int*) smalloc( g->max_atoms * sizeof( int ),
+                       "Allocate_Space_for_Grid::g->atoms[i][j]][k]" );
 
                 for ( l = 0; l < g->max_atoms; ++l )
                 {
@@ -601,22 +625,31 @@ void Cluster_Atoms( reax_system *system, static_storage *workspace,
 
     num_H = 0;
     g = &( system->g );
-    new_atoms = (reax_atom*) calloc( system->N, sizeof(reax_atom) );
-    orig_id = (int  *) calloc( system->N, sizeof( int ) );
-    f_old = (rvec*) calloc( system->N, sizeof(rvec) );
-
-    s = (real**) calloc( 3, sizeof( real* ) );
-    t = (real**) calloc( 3, sizeof( real* ) );
+    new_atoms = (reax_atom*) scalloc( system->N, sizeof(reax_atom),
+            "Cluster_Atoms::new_atoms" );
+    orig_id = (int  *) scalloc( system->N, sizeof( int ),
+            "Cluster_Atoms::orig_id" );
+    f_old = (rvec*) scalloc( system->N, sizeof(rvec),
+            "Cluster_Atoms::f_old" );
+
+    s = (real**) scalloc( 3, sizeof( real* ),
+            "Cluster_Atoms::s" );
+    t = (real**) scalloc( 3, sizeof( real* ),
+            "Cluster_Atoms::t" );
     for ( i = 0; i < 3; ++i )
     {
-        s[i] = (real *) calloc( system->N, sizeof( real ) );
-        t[i] = (real *) calloc( system->N, sizeof( real ) );
+        s[i] = (real *) scalloc( system->N, sizeof( real ),
+                "Cluster_Atoms::s[i]" );
+        t[i] = (real *) scalloc( system->N, sizeof( real ),
+                "Cluster_Atoms::t[i]" );
     }
 
-    v = (real**) calloc( control->cm_solver_restart + 1, sizeof( real* ) );
+    v = (real**) scalloc( control->cm_solver_restart + 1, sizeof( real* ),
+            "Cluster_Atoms::v" );
     for ( i = 0; i < control->cm_solver_restart + 1; ++i )
     {
-        v[i] = (real *) calloc( system->N, sizeof( real ) );
+        v[i] = (real *) scalloc( system->N, sizeof( real ),
+                "Cluster_Atoms::v[i]" );
     }
 
     top = 0;
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index c5728825..96243efe 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -279,26 +279,43 @@ void Init_Workspace( reax_system *system, control_params *control,
     int i;
 
     /* Allocate space for hydrogen bond list */
-    workspace->hbond_index = (int *) malloc( system->N * sizeof( int ) );
+    workspace->hbond_index = (int *) smalloc( system->N * sizeof( int ),
+           "Init_Workspace::workspace->hbond_index" );
 
     /* bond order related storage  */
-    workspace->total_bond_order = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Deltap           = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Deltap_boc       = (real *) malloc( system->N * sizeof( real ) );
-    workspace->dDeltap_self     = (rvec *) malloc( system->N * sizeof( rvec ) );
-
-    workspace->Delta        = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Delta_lp         = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Delta_lp_temp    = (real *) malloc( system->N * sizeof( real ) );
-    workspace->dDelta_lp        = (real *) malloc( system->N * sizeof( real ) );
-    workspace->dDelta_lp_temp   = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Delta_e          = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Delta_boc        = (real *) malloc( system->N * sizeof( real ) );
-    workspace->nlp          = (real *) malloc( system->N * sizeof( real ) );
-    workspace->nlp_temp         = (real *) malloc( system->N * sizeof( real ) );
-    workspace->Clp          = (real *) malloc( system->N * sizeof( real ) );
-    workspace->CdDelta          = (real *) malloc( system->N * sizeof( real ) );
-    workspace->vlpex        = (real *) malloc( system->N * sizeof( real ) );
+    workspace->total_bond_order = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->bond_order" );
+    workspace->Deltap = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Deltap" );
+    workspace->Deltap_boc = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Deltap_boc" );
+    workspace->dDeltap_self = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->dDeltap_self" );
+
+    workspace->Delta = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Delta" );
+    workspace->Delta_lp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Delta_lp" );
+    workspace->Delta_lp_temp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Delta_lp_temp" );
+    workspace->dDelta_lp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->dDelta_lp" );
+    workspace->dDelta_lp_temp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->dDelta_lp_temp" );
+    workspace->Delta_e = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Delta_e" );
+    workspace->Delta_boc = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Delta_boc" );
+    workspace->nlp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->nlp" );
+    workspace->nlp_temp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->nlp_temp" );
+    workspace->Clp = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->Clp" );
+    workspace->CdDelta = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->CdDelta" );
+    workspace->vlpex = (real *) smalloc( system->N * sizeof( real ),
+           "Init_Workspace::workspace->vlpex" );
 
     /* charge method storage */
     switch ( control->charge_method )
@@ -448,12 +465,16 @@ void Init_Workspace( reax_system *system, control_params *control,
     }
 
     /* integrator storage */
-    workspace->a = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_old = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->v_const = (rvec *) malloc( system->N * sizeof( rvec ) );
+    workspace->a = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->a" );
+    workspace->f_old = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_old" );
+    workspace->v_const = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->v_const" );
 
 #ifdef _OPENMP
-    workspace->f_local = (rvec *) malloc( control->num_threads * system->N * sizeof( rvec ) );
+    workspace->f_local = (rvec *) smalloc( control->num_threads * system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_local" );
 #endif
 
     /* storage for analysis */
@@ -477,20 +498,34 @@ void Init_Workspace( reax_system *system, control_params *control,
     }
 
 #ifdef TEST_FORCES
-    workspace->dDelta = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_ele = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_vdw = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_bo = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_be = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_lp = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_ov = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_un = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_ang = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_coa = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_pen = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_hb = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_tor = (rvec *) malloc( system->N * sizeof( rvec ) );
-    workspace->f_con = (rvec *) malloc( system->N * sizeof( rvec ) );
+    workspace->dDelta = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->dDelta" );
+    workspace->f_ele = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_ele" );
+    workspace->f_vdw = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_vdw" );
+    workspace->f_bo = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_bo" );
+    workspace->f_be = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_be" );
+    workspace->f_lp = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_lp" );
+    workspace->f_ov = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_ov" );
+    workspace->f_un = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_un" );
+    workspace->f_ang = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_ang" );
+    workspace->f_coa = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_coa" );
+    workspace->f_pen = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_pen" );
+    workspace->f_hb = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_hb" );
+    workspace->f_tor = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_tor" );
+    workspace->f_con = (rvec *) smalloc( system->N * sizeof( rvec ),
+           "Init_Workspace::workspace->f_con" );
 #endif
 
     workspace->realloc.num_far = -1;
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 1794cd05..01cb1798 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -148,11 +148,8 @@ void Sort_Matrix_Rows( sparse_matrix * const A )
     //    #pragma omp parallel default(none) private(i, j, si, ei, temp) shared(stderr)
 #endif
     {
-        if ( ( temp = (sparse_matrix_entry *) malloc( A->n * sizeof(sparse_matrix_entry)) ) == NULL )
-        {
-            fprintf( stderr, "Not enough space for matrix row sort. Terminating...\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        temp = (sparse_matrix_entry *) smalloc( A->n * sizeof(sparse_matrix_entry),
+               "Sort_Matrix_Rows::temp" );
 
         /* sort each row of A using column indices */
 #ifdef _OPENMP
@@ -371,11 +368,8 @@ void Calculate_Droptol( const sparse_matrix * const A,
              * overhead per Sparse_MatVec call*/
             if ( droptol_local == NULL )
             {
-                if ( (droptol_local = (real*) malloc( omp_get_num_threads() * A->n * sizeof(real))) == NULL )
-                {
-                    fprintf( stderr, "Not enough space for droptol. Terminating...\n" );
-                    exit( INSUFFICIENT_MEMORY );
-                }
+                droptol_local = (real*) smalloc( omp_get_num_threads() * A->n * sizeof(real),
+                        "Calculate_Droptol::droptol_local" );
             }
         }
 
@@ -569,15 +563,16 @@ real SuperLU_Factorize( const sparse_matrix * const A,
     {
         SUPERLU_ABORT("Malloc fails for part_super__h[].");
     }
-    if ( ( (a = (real*) malloc( (2 * A->start[A->n] - A->n) * sizeof(real))) == NULL )
-            || ( (asub = (int_t*) malloc( (2 * A->start[A->n] - A->n) * sizeof(int_t))) == NULL )
-            || ( (xa = (int_t*) malloc( (A->n + 1) * sizeof(int_t))) == NULL )
-            || ( (Ltop = (unsigned int*) malloc( (A->n + 1) * sizeof(unsigned int))) == NULL )
-            || ( (Utop = (unsigned int*) malloc( (A->n + 1) * sizeof(unsigned int))) == NULL ) )
-    {
-        fprintf( stderr, "Not enough space for SuperLU factorization. Terminating...\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    a = (real*) smalloc( (2 * A->start[A->n] - A->n) * sizeof(real),
+            "SuperLU_Factorize::a" );
+    asub = (int_t*) smalloc( (2 * A->start[A->n] - A->n) * sizeof(int_t),
+            "SuperLU_Factorize::asub" );
+    xa = (int_t*) smalloc( (A->n + 1) * sizeof(int_t),
+            "SuperLU_Factorize::xa" );
+    Ltop = (unsigned int*) smalloc( (A->n + 1) * sizeof(unsigned int),
+            "SuperLU_Factorize::Ltop" );
+    Utop = (unsigned int*) smalloc( (A->n + 1) * sizeof(unsigned int),
+            "SuperLU_Factorize::Utop" );
     if ( Allocate_Matrix( &A_t, A->n, A->m ) == FAILURE )
     {
         fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
@@ -824,13 +819,12 @@ real ICHOLT( const sparse_matrix * const A, const real * const droptol,
 
     start = Get_Time( );
 
-    if ( ( Utop = (unsigned int*) malloc((A->n + 1) * sizeof(unsigned int)) ) == NULL ||
-            ( tmp_j = (int*) malloc(A->n * sizeof(int)) ) == NULL ||
-            ( tmp_val = (real*) malloc(A->n * sizeof(real)) ) == NULL )
-    {
-        fprintf( stderr, "[ICHOLT] Not enough memory for preconditioning matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    Utop = (unsigned int*) smalloc( (A->n + 1) * sizeof(unsigned int),
+            "ICHOLT::Utop" );
+    tmp_j = (int*) smalloc( A->n * sizeof(int),
+            "ICHOLT::Utop" );
+    tmp_val = (real*) smalloc( A->n * sizeof(real),
+            "ICHOLT::Utop" );
 
     // clear variables
     Ltop = 0;
@@ -976,10 +970,10 @@ real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
 
     start = Get_Time( );
 
-    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE ||
-            ( D = (real*) malloc(A->n * sizeof(real)) ) == NULL ||
-            ( D_inv = (real*) malloc(A->n * sizeof(real)) ) == NULL ||
-            ( Utop = (int*) malloc((A->n + 1) * sizeof(int)) ) == NULL )
+    D = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D" );
+    D_inv = (real*) smalloc( A->n * sizeof(real), "ICHOL_PAR::D_inv" );
+    Utop = (int*) smalloc( (A->n + 1) * sizeof(int), "ICHOL_PAR::Utop" );
+    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
     {
         fprintf( stderr, "not enough memory for ICHOL_PAR preconditioning matrices. terminating.\n" );
         exit( INSUFFICIENT_MEMORY );
@@ -1154,9 +1148,9 @@ real ILU_PAR( const sparse_matrix * const A, const unsigned int sweeps,
 
     start = Get_Time( );
 
-    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE ||
-            ( D = (real*) malloc(A->n * sizeof(real)) ) == NULL ||
-            ( D_inv = (real*) malloc(A->n * sizeof(real)) ) == NULL )
+    D = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D" );
+    D_inv = (real*) smalloc( A->n * sizeof(real), "ILU_PAR::D_inv" );
+    if ( Allocate_Matrix( &DAD, A->n, A->m ) == FAILURE )
     {
         fprintf( stderr, "[ILU_PAR] Not enough memory for preconditioning matrices. Terminating.\n" );
         exit( INSUFFICIENT_MEMORY );
@@ -1377,12 +1371,8 @@ real ILUT_PAR( const sparse_matrix * const A, const real * droptol,
         exit( INSUFFICIENT_MEMORY );
     }
 
-    if ( ( D = (real*) malloc(A->n * sizeof(real)) ) == NULL ||
-            ( D_inv = (real*) malloc(A->n * sizeof(real)) ) == NULL )
-    {
-        fprintf( stderr, "not enough memory for ILUT_PAR preconditioning matrices. terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    D = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D" );
+    D_inv = (real*) smalloc( A->n * sizeof(real), "ILUT_PAR::D_inv" );
 
 #ifdef _OPENMP
     #pragma omp parallel for schedule(static) \
@@ -1822,10 +1812,8 @@ static void Sparse_MatVec( const sparse_matrix * const A,
          * overhead per Sparse_MatVec call*/
         if ( b_local == NULL )
         {
-            if ( (b_local = (real*) malloc( omp_get_num_threads() * n * sizeof(real))) == NULL )
-            {
-                exit( INSUFFICIENT_MEMORY );
-            }
+            b_local = (real*) smalloc( omp_get_num_threads() * n * sizeof(real),
+                    "Sparse_MatVec::b_local" );
         }
     }
 
@@ -2087,22 +2075,18 @@ void tri_solve_level_sched( const sparse_matrix * const LU,
 
         if ( row_levels == NULL || level_rows == NULL || level_rows_cnt == NULL )
         {
-            if ( (row_levels = (unsigned int*) malloc((size_t)N * sizeof(unsigned int))) == NULL
-                    || (level_rows = (unsigned int*) malloc((size_t)N * sizeof(unsigned int))) == NULL
-                    || (level_rows_cnt = (unsigned int*) malloc((size_t)(N + 1) * sizeof(unsigned int))) == NULL )
-            {
-                fprintf( stderr, "Not enough space for triangular solve via level scheduling. Terminating...\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            row_levels = (unsigned int*) smalloc( (size_t)N * sizeof(unsigned int),
+                    "tri_solve_level_sched::row_levels" );
+            level_rows = (unsigned int*) smalloc( (size_t)N * sizeof(unsigned int),
+                    "tri_solve_level_sched::level_rows" );
+            level_rows_cnt = (unsigned int*) smalloc( (size_t)(N + 1) * sizeof(unsigned int),
+                    "tri_solve_level_sched::level_rows_cnt" );
         }
 
         if ( top == NULL )
         {
-            if ( (top = (unsigned int*) malloc((size_t)(N + 1) * sizeof(unsigned int))) == NULL )
-            {
-                fprintf( stderr, "Not enough space for triangular solve via level scheduling. Terminating...\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            top = (unsigned int*) smalloc( (size_t)(N + 1) * sizeof(unsigned int),
+                    "tri_solve_level_sched::top" );
         }
 
         /* find levels (row dependencies in substitutions) */
@@ -2338,12 +2322,10 @@ void graph_coloring( const sparse_matrix * const A, const TRIANGULARITY tri )
             }
         }
 
-        if ( (fb_color = (int*) malloc(sizeof(int) * MAX_COLOR)) == NULL ||
-                (conflict_local = (unsigned int*) malloc(sizeof(unsigned int) * A->n)) == NULL )
-        {
-            fprintf( stderr, "not enough memory for graph coloring. terminating.\n" );
-            exit( INSUFFICIENT_MEMORY );
-        }
+        fb_color = (int*) smalloc( sizeof(int) * MAX_COLOR,
+                 "graph_coloring::fb_color" );
+        conflict_local = (unsigned int*) smalloc( sizeof(unsigned int) * A->n,
+                "graph_coloring::fb_color" );
 
 #ifdef _OPENMP
         #pragma omp barrier
@@ -2526,11 +2508,7 @@ static void permute_vector( real * const x, const unsigned int n, const int inve
     {
         if ( x_p == NULL )
         {
-            if ( (x_p = (real*) malloc(sizeof(real) * n)) == NULL )
-            {
-                fprintf( stderr, "not enough memory for permuting vector. terminating.\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            x_p = (real*) smalloc( sizeof(real) * n, "permute_vector::x_p" );
         }
 
         if ( invert_map == TRUE )
@@ -2717,16 +2695,25 @@ sparse_matrix * setup_graph_coloring( sparse_matrix * const H )
 #endif
 
         /* internal storage for graph coloring (global to facilitate simultaneous access to OpenMP threads) */
-        if ( (color = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (to_color = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (conflict = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (conflict_cnt = (unsigned int*) malloc(sizeof(unsigned int) * (num_thread + 1))) == NULL ||
-                (recolor = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (color_top = (unsigned int*) malloc(sizeof(unsigned int) * (H->n + 1))) == NULL ||
-                (permuted_row_col = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (permuted_row_col_inv = (unsigned int*) malloc(sizeof(unsigned int) * H->n)) == NULL ||
-                (y_p = (real*) malloc(sizeof(real) * H->n)) == NULL ||
-                (Allocate_Matrix( &H_p, H->n, H->m ) == FAILURE ) ||
+        color = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::color" );
+        to_color = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::to_color" );
+        conflict = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::conflict" );
+        conflict_cnt = (unsigned int*) smalloc( sizeof(unsigned int) * (num_thread + 1),
+                "setup_graph_coloring::conflict_cnt" );
+        recolor = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::recolor" );
+        color_top = (unsigned int*) smalloc( sizeof(unsigned int) * (H->n + 1),
+                "setup_graph_coloring::color_top" );
+        permuted_row_col = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::premuted_row_col" );
+        permuted_row_col_inv = (unsigned int*) smalloc( sizeof(unsigned int) * H->n,
+                "setup_graph_coloring::premuted_row_col_inv" );
+        y_p = (real*) smalloc( sizeof(real) * H->n,
+                "setup_graph_coloring::y_p" );
+        if ( (Allocate_Matrix( &H_p, H->n, H->m ) == FAILURE ) ||
                 (Allocate_Matrix( &H_full, H->n, 2 * H->m - H->n ) == FAILURE ) )
         {
             fprintf( stderr, "not enough memory for graph coloring. terminating.\n" );
@@ -2773,27 +2760,15 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
     {
         if ( Dinv_b == NULL )
         {
-            if ( (Dinv_b = (real*) malloc(sizeof(real) * R->n)) == NULL )
-            {
-                fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            Dinv_b = (real*) smalloc( sizeof(real) * R->n, "jacobi_iter::Dinv_b" );
         }
         if ( rp == NULL )
         {
-            if ( (rp = (real*) malloc(sizeof(real) * R->n)) == NULL )
-            {
-                fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            rp = (real*) smalloc( sizeof(real) * R->n, "jacobi_iter::rp" );
         }
         if ( rp2 == NULL )
         {
-            if ( (rp2 = (real*) malloc(sizeof(real) * R->n)) == NULL )
-            {
-                fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
-                exit( INSUFFICIENT_MEMORY );
-            }
+            rp2 = (real*) smalloc( sizeof(real) * R->n, "jacobi_iter::rp2" );
         }
     }
 
@@ -2969,11 +2944,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
             {
                 if ( Dinv_L == NULL )
                 {
-                    if ( (Dinv_L = (real*) malloc(sizeof(real) * workspace->L->n)) == NULL )
-                    {
-                        fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
-                        exit( INSUFFICIENT_MEMORY );
-                    }
+                    Dinv_L = (real*) smalloc( sizeof(real) * workspace->L->n,
+                            "apply_preconditioner::Dinv_L" );
                 }
             }
 
@@ -2998,11 +2970,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
             {
                 if ( Dinv_U == NULL )
                 {
-                    if ( (Dinv_U = (real*) malloc(sizeof(real) * workspace->U->n)) == NULL )
-                    {
-                        fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
-                        exit( INSUFFICIENT_MEMORY );
-                    }
+                    Dinv_U = (real*) smalloc( sizeof(real) * workspace->U->n,
+                            "apply_preconditioner::Dinv_U" );
                 }
             }
 
@@ -3729,11 +3698,7 @@ real condest( const sparse_matrix * const L, const sparse_matrix * const U )
 
     N = L->n;
 
-    if ( (e = (real*) malloc(sizeof(real) * N)) == NULL )
-    {
-        fprintf( stderr, "Not enough memory for condest. Terminating.\n" );
-        exit( INSUFFICIENT_MEMORY );
-    }
+    e = (real*) smalloc( sizeof(real) * N, "condest::e" );
 
     memset( e, 1., N * sizeof(real) );
 
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 7ae8ab2e..7990169c 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -36,7 +36,8 @@ void Make_Lookup_Table( real xmin, real xmax, int n,
     t->dx = (xmax - xmin) / (n - 1);
     t->inv_dx = 1.0 / t->dx;
     t->a = (n - 1) / (xmax - xmin);
-    t->y = (real*) malloc(n * sizeof(real));
+    t->y = (real*) smalloc( n * sizeof(real),
+            "Make_Lookup_Table::t->y" );
 
     for (i = 0; i < n; i++)
     {
@@ -83,11 +84,16 @@ void Natural_Cubic_Spline( const real *h, const real *f,
     real *a, *b, *c, *d, *v;
 
     /* allocate space for the linear system */
-    a = (real*) malloc( n * sizeof(real) );
-    b = (real*) malloc( n * sizeof(real) );
-    c = (real*) malloc( n * sizeof(real) );
-    d = (real*) malloc( n * sizeof(real) );
-    v = (real*) malloc( n * sizeof(real) );
+    a = (real*) smalloc( n * sizeof(real),
+           "Natural_Cubic_Spline::a" );
+    b = (real*) smalloc( n * sizeof(real),
+           "Natural_Cubic_Spline::b" );
+    c = (real*) smalloc( n * sizeof(real),
+           "Natural_Cubic_Spline::c" );
+    d = (real*) smalloc( n * sizeof(real),
+           "Natural_Cubic_Spline::d" );
+    v = (real*) smalloc( n * sizeof(real),
+           "Natural_Cubic_Spline::v" );
 
     /* build the linear system */
     a[0] = 0.0;
@@ -153,11 +159,16 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast,
     real *a, *b, *c, *d, *v;
 
     /* allocate space for the linear system */
-    a = (real*) malloc( n * sizeof(real) );
-    b = (real*) malloc( n * sizeof(real) );
-    c = (real*) malloc( n * sizeof(real) );
-    d = (real*) malloc( n * sizeof(real) );
-    v = (real*) malloc( n * sizeof(real) );
+    a = (real*) smalloc( n * sizeof(real),
+           "Complete_Cubic_Spline::a" );
+    b = (real*) smalloc( n * sizeof(real),
+           "Complete_Cubic_Spline::b" );
+    c = (real*) smalloc( n * sizeof(real),
+           "Complete_Cubic_Spline::c" );
+    d = (real*) smalloc( n * sizeof(real),
+           "Complete_Cubic_Spline::d" );
+    v = (real*) smalloc( n * sizeof(real),
+           "Complete_Cubic_Spline::v" );
 
     /* build the linear system */
     a[0] = 0.0;
@@ -269,10 +280,12 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
 
     /* allocate Long-Range LookUp Table space based on
        number of atom types in the ffield file */
-    LR = (LR_lookup_table**) malloc( num_atom_types * sizeof(LR_lookup_table*) );
+    LR = (LR_lookup_table**) smalloc( num_atom_types * sizeof(LR_lookup_table*),
+           "Make_LR_Lookup_Table::LR" );
     for ( i = 0; i < num_atom_types; ++i )
     {
-        LR[i] = (LR_lookup_table*) malloc(num_atom_types * sizeof(LR_lookup_table));
+        LR[i] = (LR_lookup_table*) smalloc( num_atom_types * sizeof(LR_lookup_table),
+                "Make_LR_Lookup_Table::LR[i]");
     }
 
     /* most atom types in ffield file will not exist in the current
@@ -303,17 +316,23 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control )
                     LR[i][j].dx = dr;
                     LR[i][j].inv_dx = control->tabulate / control->r_cut;
                     LR[i][j].y = (LR_data*)
-                        malloc( LR[i][j].n * sizeof(LR_data) );
+                        smalloc( LR[i][j].n * sizeof(LR_data),
+                              "Make_LR_Lookup_Table::LR[i][j].y" );
                     LR[i][j].H = (cubic_spline_coef*)
-                        malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
+                        smalloc( LR[i][j].n * sizeof(cubic_spline_coef),
+                              "Make_LR_Lookup_Table::LR[i][j].H" );
                     LR[i][j].vdW = (cubic_spline_coef*)
-                        malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
+                        smalloc( LR[i][j].n * sizeof(cubic_spline_coef),
+                              "Make_LR_Lookup_Table::LR[i][j].vdW" );
                     LR[i][j].CEvd = (cubic_spline_coef*)
-                        malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
+                        smalloc( LR[i][j].n * sizeof(cubic_spline_coef),
+                              "Make_LR_Lookup_Table::LR[i][j].CEvd" );
                     LR[i][j].ele = (cubic_spline_coef*)
-                        malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
+                        smalloc( LR[i][j].n * sizeof(cubic_spline_coef),
+                              "Make_LR_Lookup_Table::LR[i][j].ele" );
                     LR[i][j].CEclmb = (cubic_spline_coef*)
-                        malloc( LR[i][j].n * sizeof(cubic_spline_coef) );
+                        smalloc( LR[i][j].n * sizeof(cubic_spline_coef),
+                              "Make_LR_Lookup_Table::LR[i][j].CEclmb" );
 
                     for ( r = 1; r <= control->tabulate; ++r )
                     {
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 2e791a30..82958abe 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -154,7 +154,8 @@ void static Read_System( char * const geo_file,
 int Setup( char ** args, reax_system * const system, control_params * const control,
         simulation_data * const data )
 {
-    lists = (reax_list*) malloc( sizeof(reax_list) * LIST_N );
+    lists = (reax_list*) smalloc( sizeof(reax_list) * LIST_N,
+           "Setup::lists" );
 
     Read_System( args[0], args[1], args[2], system, control,
             data, &workspace, &out_control );
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index c37ccc30..6b59d194 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -352,34 +352,22 @@ char *Get_Atom_Name( reax_system *system, int i )
 }
 
 
-int Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
+void Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
 {
     int i;
 
-    if ( (*line = (char*) malloc( sizeof(char) * MAX_LINE )) == NULL )
-    {
-        return FAILURE;
-    }
-
-    if ( (*backup = (char*) malloc( sizeof(char) * MAX_LINE )) == NULL )
-    {
-        return FAILURE;
-    }
-
-    if ( (*tokens = (char**) malloc( sizeof(char*) * MAX_TOKENS )) == NULL )
-    {
-        return FAILURE;
-    }
+    *line = (char*) smalloc( sizeof(char) * MAX_LINE,
+            "Allocate_Tokenizer_Space::*line" );
+    *backup = (char*) smalloc( sizeof(char) * MAX_LINE,
+            "Allocate_Tokenizer_Space::*backup" );
+    *tokens = (char**) smalloc( sizeof(char*) * MAX_TOKENS,
+            "Allocate_Tokenizer_Space::*tokens" );
 
     for ( i = 0; i < MAX_TOKENS; i++ )
     {
-        if ( ((*tokens)[i] = (char*) malloc(sizeof(char) * MAX_TOKEN_LEN)) == NULL )
-        {
-            return FAILURE;
-        }
+        (*tokens)[i] = (char*) smalloc(sizeof(char) * MAX_TOKEN_LEN,
+                "Allocate_Tokenizer_Space::(*tokens)[i]" );
     }
-
-    return SUCCESS;
 }
 
 
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index ba103f9a..d1b15f64 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -72,7 +72,7 @@ char *Get_Element( reax_system*, int );
 
 char *Get_Atom_Name( reax_system*, int );
 
-int Allocate_Tokenizer_Space( char**, char**, char*** );
+void Allocate_Tokenizer_Space( char**, char**, char*** );
 
 void Deallocate_Tokenizer_Space( char **, char **, char *** );
 
-- 
GitLab