diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 2531612a65ccfe2019e97bb20b743dc84cabc384..92c1d80dba39f85a214db63c377b210d741f2843 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -56,6 +56,14 @@ void PreAllocate_Space( reax_system * const system,
                         "PreAllocate_Space::workspace->restricted_list[i]" );
             }
         }
+
+        if ( control->geo_format == BGF
+                || control->geo_format == ASCII_RESTART
+                || control->geo_format == BINARY_RESTART )
+        {
+            workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
+                    "Read_BGF::workspace->map_serials" );
+        }
     }
     else
     {
@@ -156,20 +164,16 @@ static void Reallocate_Matrix( sparse_matrix *H, int n, int n_max, int m )
 }
 
 
-void Allocate_HBond_List( int n, int num_h, int num_h_max,
-        int const * const h_index, int * const hb_top, reax_list * const hbond_list )
+void Initialize_HBond_List( int n, int const * const h_index,
+        int * const hb_top, reax_list * const hbond_list )
 {
-    int i, num_hbonds;
+    int i;
 
-    num_hbonds = 0;
     /* find starting indexes for each H and the total number of hbonds */
     for ( i = 1; i < n; ++i )
     {
         hb_top[i] += hb_top[i - 1];
     }
-    num_hbonds = hb_top[n - 1];
-
-    Make_List( num_h, num_h_max, num_hbonds, TYP_HBOND, hbond_list );
 
     for ( i = 0; i < n; ++i )
     {
@@ -187,18 +191,22 @@ void Allocate_HBond_List( int n, int num_h, int num_h_max,
 }
 
 
-static void Reallocate_HBonds_List( int n, int num_h, int num_h_max,
+static void Reallocate_Initialize_HBond_List( int n, int num_h, int num_h_max,
         int *h_index, reax_list *hbond_list )
 {
-    int i;
-    int *hb_top;
+    int i, num_hbonds, *hb_top;
+
+    hb_top = scalloc( n, sizeof(int),
+            "Reallocate_Initialize_HBond_List::hb_top" );
+    num_hbonds = 0;
 
-    hb_top = scalloc( n, sizeof(int), "Reallocate_HBonds_List::hb_top" );
     for ( i = 0; i < n; ++i )
     {
         if ( h_index[i] >= 0 )
         {
-            hb_top[i] = MAX(Num_Entries(h_index[i], hbond_list) * SAFE_HBONDS, MIN_HBONDS);
+            hb_top[i] = MAX( Num_Entries( h_index[i], hbond_list ) * SAFE_HBONDS,
+                    MIN_HBONDS );
+            num_hbonds += hb_top[i];
         }
     }
 
@@ -206,31 +214,28 @@ static void Reallocate_HBonds_List( int n, int num_h, int num_h_max,
     {
         Delete_List( TYP_HBOND, hbond_list );
     }
+    Make_List( num_h, num_h_max, num_hbonds, TYP_HBOND, hbond_list );
 
-    Allocate_HBond_List( n, num_h, num_h_max, h_index, hb_top, hbond_list );
+    Initialize_HBond_List( n, h_index, hb_top, hbond_list );
 
-    sfree( hb_top, "Reallocate_HBonds_List::hb_top" );
+    sfree( hb_top, "Reallocate_Initialize_HBond_List::hb_top" );
 }
 
 
-void Allocate_Bond_List( int n, int n_max, int * const bond_top,
+void Initialize_Bond_List( int * const bond_top,
         reax_list * const bond_list )
 {
-    int i, num_bonds;
+    int i;
 
-    num_bonds = 0;
     /* find starting indexes for each atom and the total number of bonds */
-    for ( i = 1; i < n; ++i )
+    for ( i = 1; i < bond_list->n; ++i )
     {
         bond_top[i] += bond_top[i - 1];
     }
-    num_bonds = bond_top[n - 1];
-
-    Make_List( n, n_max, num_bonds, TYP_BOND, bond_list );
 
     Set_Start_Index( 0, 0, bond_list );
     Set_End_Index( 0, 0, bond_list );
-    for ( i = 1; i < n; ++i )
+    for ( i = 1; i < bond_list->n; ++i )
     {
         Set_Start_Index( i, bond_top[i - 1], bond_list );
         Set_End_Index( i, bond_top[i - 1], bond_list );
@@ -238,30 +243,34 @@ void Allocate_Bond_List( int n, int n_max, int * const bond_top,
 }
 
 
-static void Reallocate_Bonds_List( int n, int n_max, reax_list *bond_list,
-        int *num_bonds, int *est_3body )
+static void Reallocate_Initialize_Bond_List( int n, int n_max,
+        reax_list *bond_list, int *num_bonds, int *est_3body )
 {
     int i;
     int *bond_top;
 
-    bond_top = (int *) scalloc( n, sizeof(int), "Reallocate_Bonds_List::hb_top" );
+    bond_top = (int *) scalloc( n, sizeof(int),
+            "Reallocate_Initialize_Bond_List::hb_top" );
+    *num_bonds = 0;
     *est_3body = 0;
 
     for ( i = 0; i < n; ++i )
     {
         *est_3body += SQR( Num_Entries( i, bond_list ) );
         bond_top[i] = MAX( Num_Entries( i, bond_list ) * 2, MIN_BONDS );
+        *num_bonds += bond_top[i];
     }
 
     if ( bond_list->allocated == TRUE )
     {
         Delete_List( TYP_BOND, bond_list );
     }
+    Make_List( n, n_max, (int) CEIL( *num_bonds * SAFE_ZONE ),
+            TYP_BOND, bond_list );
 
-    Allocate_Bond_List( n, n_max, bond_top, bond_list );
-    *num_bonds = bond_top[n - 1];
+    Initialize_Bond_List( bond_top, bond_list );
 
-    sfree( bond_top, "Reallocate_Bonds_List::bond_top" );
+    sfree( bond_top, "Reallocate_Initialize_Bond_List::bond_top" );
 }
 
 
@@ -296,15 +305,15 @@ void Reallocate( reax_system * const system, control_params const * const contro
 
     if ( control->hbond_cut > 0.0 && realloc->hbonds > 0 )
     {
-        Reallocate_HBonds_List( system->N, workspace->num_H, workspace->num_H_max,
-                workspace->hbond_index, lists[HBONDS] );
+        Reallocate_Initialize_HBond_List( system->N, workspace->num_H,
+                workspace->num_H_max, workspace->hbond_index, lists[HBONDS] );
         realloc->hbonds = -1;
     }
 
     num_bonds = est_3body = -1;
     if ( realloc->bonds > 0 )
     {
-        Reallocate_Bonds_List( system->N, system->N_max, lists[BONDS], &num_bonds, &est_3body );
+        Reallocate_Initialize_Bond_List( system->N, system->N_max, lists[BONDS], &num_bonds, &est_3body );
         realloc->bonds = -1;
         realloc->num_3body = MAX( realloc->num_3body, est_3body );
     }
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index 1716bf42ecb307686a1f6b935f09c7298b1e605e..05ee7fdc4802cae80cc2e94313543b487d5c5cfd 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -35,10 +35,10 @@ void Allocate_Matrix( sparse_matrix * const , int, int, int );
 
 void Deallocate_Matrix( sparse_matrix * const  );
 
-void Allocate_HBond_List( int, int, int, int const * const, int * const,
+void Initialize_HBond_List( int, int const * const, int * const,
         reax_list * const );
 
-void Allocate_Bond_List( int, int, int * const, reax_list * const );
+void Initialize_Bond_List( int * const, reax_list * const );
 
 
 #endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index db62f2a51507cfa26636188df75240e53899d278..8dfc1af0687f9f058201676859a89f66903dd7c1 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -526,7 +526,7 @@ void Set_Control_Defaults( reax_system * const system,
         control->ignore[i] = 0;
     }
 
-    out_control->log_update_freq = 0;
+    out_control->log_update_freq = 1;
     out_control->write_steps = 0;
     out_control->traj_compress = 0;
     out_control->write = &fprintf;
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 4f929e9abf6762830697e166a1f265ec71837baf..f01e1e4a1beb03aeb93a86dfae2b4d8eb0158ccb 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -302,9 +302,9 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
 
     /* count atoms and allocate storage */
     n = Count_Atoms( system, geo, CUSTOM );
-    if ( system->prealloc_allocated == FALSE || n > system->N )
+    if ( system->prealloc_allocated == FALSE || n > system->N_max )
     {
-        PreAllocate_Space( system, control, workspace, n );
+        PreAllocate_Space( system, control, workspace, (int) CEIL( SAFE_ZONE * n ) );
     }
     system->N = n;
 
@@ -367,9 +367,9 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
     }
 
     n = Count_Atoms( system, pdb, PDB );
-    if ( system->prealloc_allocated == FALSE || n > system->N )
+    if ( system->prealloc_allocated == FALSE || n > system->N_max )
     {
-        PreAllocate_Space( system, control, workspace, n );
+        PreAllocate_Space( system, control, workspace, (int) CEIL( SAFE_ZONE * n ) );
     }
     system->N = n;
 
@@ -506,13 +506,15 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
 //                        "CONECT line exceeds max num restrictions allowed.\n" );
 
                 /* read bond restrictions */
-                // if( is_Valid_Serial( workspace, pdb_serial = sstrtol( tmp[1]), __FILE__, __LINE__ ) )
+                // pdb_serial = sstrtol( tmp[1], __FILE__, __LINE__ );
+                // if ( is_Valid_Serial( workspace->map_serials[ pdb_serial ] ) == TRUE )
                 //   ratom = workspace->map_serials[ pdb_serial ];
 
                 // workspace->restricted[ ratom ] = c1 - 2;
-                // for( i = 2; i < c1; ++i )
+                // for ( i = 2; i < c1; ++i )
                 //  {
-                //    if( is_Valid_Serial(workspace, pdb_serial = sstrtol( tmp[i], __FILE__, __LINE__ )) )
+                //    pdb_serial = sstrtol( tmp[i], __FILE__, __LINE__ );
+                //    if ( is_Valid_Serial( workspace->map_serials[ pdb_serial ] ) == TRUE )
                 //        workspace->restricted_list[ ratom ][ i-2 ] =
                 //          workspace->map_serials[ pdb_serial ];
                 //  }
@@ -659,14 +661,12 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
 
     sfclose( bgf, "Read_BGF::bgf" );
 
-    if ( system->prealloc_allocated == FALSE || n > system->N )
+    if ( system->prealloc_allocated == FALSE || n > system->N_max )
     {
-        PreAllocate_Space( system, control, workspace, n );
+        PreAllocate_Space( system, control, workspace, (int) CEIL( SAFE_ZONE * n ) );
     }
     system->N = n;
 
-    workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
-            "Read_BGF::workspace->map_serials" );
     for ( i = 0; i < MAX_ATOM_ID; ++i )
     {
         workspace->map_serials[i] = -1;
@@ -786,7 +786,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
 
                 /* read bond restrictions */
                 bgf_serial = sstrtol( tokens[1], __FILE__, __LINE__ );
-                if ( is_Valid_Serial( workspace, bgf_serial ) )
+                if ( is_Valid_Serial( workspace->map_serials[ bgf_serial ] ) == TRUE )
                 {
                     ratom = workspace->map_serials[ bgf_serial ];
                 }
@@ -794,7 +794,8 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
                 workspace->restricted[ ratom ] = token_cnt - 2;
                 for ( i = 2; i < token_cnt; ++i )
                 {
-                    if ( is_Valid_Serial( workspace, bgf_serial = sstrtol( tokens[i], __FILE__, __LINE__ )) )
+                    bgf_serial = sstrtol( tokens[i], __FILE__, __LINE__ );
+                    if ( is_Valid_Serial( workspace->map_serials[ bgf_serial ] ) == TRUE )
                     {
                         workspace->restricted_list[ratom][i - 2] =
                             workspace->map_serials[bgf_serial];
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 16eb339829d948b619bc6163aae14d7f04d415dc..309ee0b186481dc19afbd7b21262003fe9ad0076 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -787,7 +787,9 @@ static void Init_Lists( reax_system *system, control_params *control,
         {
             Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
         }
-        Make_List( system->N, system->N_max, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+        Make_List( system->N, system->N_max, 
+                MAX( num_nbrs, lists[FAR_NBRS]->total_intrs ),
+                TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
     }
     else
     {
@@ -895,8 +897,9 @@ static void Init_Lists( reax_system *system, control_params *control,
             {
                 workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
 
-                Allocate_HBond_List( system->N, workspace->num_H, workspace->num_H_max,
-                        workspace->hbond_index, hb_top, lists[HBONDS] );
+                Make_List( workspace->num_H, workspace->num_H_max,
+                        (int) CEIL( SAFE_ZONE * num_hbonds ),
+                        TYP_HBOND, lists[HBONDS] );
             }
             else if ( workspace->num_H_max < workspace->num_H
                     || lists[HBONDS]->total_intrs < num_hbonds )
@@ -910,21 +913,24 @@ static void Init_Lists( reax_system *system, control_params *control,
                 {
                     Delete_List( TYP_HBOND, lists[HBONDS] );
                 }
-                Allocate_HBond_List( system->N, workspace->num_H, workspace->num_H_max,
-                        workspace->hbond_index, hb_top, lists[HBONDS] );
+                Make_List( workspace->num_H, workspace->num_H_max,
+                        MAX( num_hbonds, lists[HBONDS]->total_intrs ),
+                        TYP_HBOND, lists[HBONDS] );
             }
             else
             {
                 lists[HBONDS]->n = workspace->num_H;
             }
+
+            Initialize_HBond_List( system->N, workspace->hbond_index, hb_top, lists[HBONDS] );
         }
+    }
 
 #if defined(DEBUG_FOCUS)
-        fprintf( stderr, "estimated storage - num_hbonds: %d\n", num_hbonds );
-        fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
-                 num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
+    fprintf( stderr, "estimated storage - num_hbonds: %d\n", num_hbonds );
+    fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
+             num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
 #endif
-    }
 
     num_bonds = 0;
     for ( i = 0; i < system->N; ++i )
@@ -935,7 +941,8 @@ static void Init_Lists( reax_system *system, control_params *control,
     /* bonds list */
     if ( lists[BONDS]->allocated == FALSE )
     {
-        Allocate_Bond_List( system->N, system->N_max, bond_top, lists[BONDS] );
+        Make_List( system->N, system->N_max, (int) CEIL( num_bonds * SAFE_ZONE ),
+                TYP_BOND, lists[BONDS] );
     }
     else if ( realloc == TRUE || lists[BONDS]->total_intrs < num_bonds )
     {
@@ -943,13 +950,17 @@ static void Init_Lists( reax_system *system, control_params *control,
         {
             Delete_List( TYP_BOND, lists[BONDS] );
         }
-        Allocate_Bond_List( system->N, system->N_max, bond_top, lists[BONDS] );
+        Make_List( system->N, system->N_max,
+                MAX( num_bonds, lists[BONDS]->total_intrs ),
+                TYP_BOND, lists[BONDS] );
     }
     else
     {
         lists[BONDS]->n = system->N;
     }
 
+    Initialize_Bond_List( bond_top, lists[BONDS] );
+
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "estimated storage - num_bonds: %d\n", num_bonds );
     fprintf( stderr, "memory allocated: bonds = %ldMB\n",
@@ -968,17 +979,10 @@ static void Init_Lists( reax_system *system, control_params *control,
         {
             Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
         }
-
-        if ( lists[THREE_BODIES]->n_max < num_bonds )
-        {
-            Make_List( num_bonds, num_bonds, num_3body,
-                    TYP_THREE_BODY, lists[THREE_BODIES] );
-        }
-        else
-        {
-            Make_List( num_bonds, lists[THREE_BODIES]->n_max, num_3body,
-                    TYP_THREE_BODY, lists[THREE_BODIES] );
-        }
+        Make_List( MAX( num_bonds, lists[THREE_BODIES]->n_max),
+                MAX( num_bonds, lists[THREE_BODIES]->n_max),
+                MAX( num_3body, lists[THREE_BODIES]->total_intrs ),
+                TYP_THREE_BODY, lists[THREE_BODIES] );
     }
     else
     {
@@ -1409,23 +1413,41 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
     }
 
     if ( workspace->H.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H );
+    }
     if ( workspace->H_full.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_full );
+    }
     if ( workspace->H_sp.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_sp );
+    }
     if ( workspace->H_p.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_p );
+    }
     if ( workspace->H_spar_patt.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_spar_patt );
+    }
     if ( workspace->H_spar_patt_full.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_spar_patt_full );
+    }
     if ( workspace->H_app_inv.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->H_app_inv );
+    }
     if ( workspace->L.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->L );
+    }
     if ( workspace->U.allocated == TRUE )
+    {
         Deallocate_Matrix( &workspace->U );
+    }
 
     for ( i = 0; i < 5; ++i )
     {
diff --git a/sPuReMD/src/restart.c b/sPuReMD/src/restart.c
index 34ad174c5fa00b503ce9f012facac8e6692f738e..43ef92027138ace8de2199b232732c068fdd23a8 100644
--- a/sPuReMD/src/restart.c
+++ b/sPuReMD/src/restart.c
@@ -108,18 +108,13 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
              system->box.box[2][0], system->box.box[2][1], system->box.box[2][2] );
 #endif
 
-    if ( system->prealloc_allocated == FALSE || res_header.N > system->N )
+    if ( system->prealloc_allocated == FALSE || res_header.N > system->N_max )
     {
-        PreAllocate_Space( system, control, workspace, res_header.N );
+        PreAllocate_Space( system, control, workspace,
+                (int) CEIL( SAFE_ZONE * res_header.N ) );
     }
     system->N = res_header.N;
 
-    if ( system->prealloc_allocated == FALSE )
-    {
-        workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
-                "Read_Binary_Restart::workspace->map_serials" );
-    }
-
     for ( i = 0; i < MAX_ATOM_ID; ++i )
     {
         workspace->map_serials[i] = -1;
@@ -222,18 +217,12 @@ void Read_ASCII_Restart( const char * const fname, reax_system *system,
              system->box.box[2][0], system->box.box[2][1], system->box.box[2][2] );
 #endif
 
-    if ( system->prealloc_allocated == FALSE || n > system->N )
+    if ( system->prealloc_allocated == FALSE || n > system->N_max )
     {
-        PreAllocate_Space( system, control, workspace, n );
+        PreAllocate_Space( system, control, workspace, (int) CEIL( SAFE_ZONE * n ) );
     }
     system->N = n;
 
-    if ( system->prealloc_allocated == FALSE )
-    {
-        workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
-                "Read_ASCII_Restart::workspace->map_serials" );
-    }
-
     for ( i = 0; i < MAX_ATOM_ID; ++i )
     {
         workspace->map_serials[i] = -1;
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index d9bfb001d943e5a85cbf9f61915661c197aac75f..8f38bc46a743607e181b05b1833fff1fef75c331 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -72,13 +72,7 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
         Compute_Kinetic_Energy( system, data );
     }
 
-    if ( (out_control->log_update_freq > 0
-                && data->step % out_control->log_update_freq == 0)
-            || (out_control->write_steps > 0
-                && data->step % out_control->write_steps == 0) )
-    {
-        Compute_Total_Energy( data );
-    }
+    Compute_Total_Energy( data );
 
     if ( control->compute_pressure == TRUE && control->ensemble != sNPT
             && control->ensemble != iNPT && control->ensemble != aNPT )
@@ -98,21 +92,27 @@ static void Read_Input_Files( const char * const geo_file,
         const char * const ffield_file, const char * const control_file,
         reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        output_controls * const out_control )
+        output_controls * const out_control, int reset )
 {
     if ( ffield_file != NULL )
     {
         Read_Force_Field( ffield_file, system, &system->reax_param );
     }
 
-    Set_Control_Defaults( system, control, out_control );
+    if ( reset == FALSE || control_file != NULL )
+    {
+        Set_Control_Defaults( system, control, out_control );
+    }
 
     if ( control_file != NULL )
     {
         Read_Control_File( control_file, system, control, out_control );
     }
 
-    Set_Control_Derived_Values( system, control );
+    if ( reset == FALSE || control_file != NULL )
+    {
+        Set_Control_Derived_Values( system, control );
+    }
 
     if ( geo_file != NULL )
     {
@@ -151,6 +151,70 @@ static void Read_Input_Files( const char * const geo_file,
 }
 
 
+static void Allocate_Top_Level_Structs( spuremd_handle ** handle )
+{
+    int i;
+
+    /* top-level allocation */
+    *handle = smalloc( sizeof(spuremd_handle), "Allocate_Top_Level_Structs::handle" );
+
+    /* second-level allocations */
+    (*handle)->system = smalloc( sizeof(reax_system),
+           "Allocate_Top_Level_Structs::handle->system" );
+
+    (*handle)->control = smalloc( sizeof(control_params),
+           "Allocate_Top_Level_Structs::handle->control" );
+
+    (*handle)->data = smalloc( sizeof(simulation_data),
+           "Allocate_Top_Level_Structs::handle->data" );
+
+    (*handle)->workspace = smalloc( sizeof(static_storage),
+           "Allocate_Top_Level_Structs::handle->workspace" );
+
+    (*handle)->lists = smalloc( sizeof(reax_list *) * LIST_N,
+           "Allocate_Top_Level_Structs::handle->lists" );
+    for ( i = 0; i < LIST_N; ++i )
+    {
+        (*handle)->lists[i] = smalloc( sizeof(reax_list),
+                "Allocate_Top_Level_Structs::handle->lists[i]" );
+    }
+    (*handle)->out_control = smalloc( sizeof(output_controls),
+           "Allocate_Top_Level_Structs::handle->out_control" );
+}
+
+
+static void Initialize_Top_Level_Structs( spuremd_handle * handle )
+{
+    int i;
+
+    /* top-level initializations */
+    handle->output_enabled = TRUE;
+    handle->realloc = TRUE;
+    handle->callback = NULL;
+    handle->data->sim_id = 0;
+
+    /* second-level initializations */
+    handle->system->prealloc_allocated = FALSE;
+    handle->system->ffield_params_allocated = FALSE;
+    handle->system->g.allocated = FALSE;
+
+    handle->workspace->H.allocated = FALSE;
+    handle->workspace->H_full.allocated = FALSE;
+    handle->workspace->H_sp.allocated = FALSE;
+    handle->workspace->H_p.allocated = FALSE;
+    handle->workspace->H_spar_patt.allocated = FALSE;
+    handle->workspace->H_spar_patt_full.allocated = FALSE;
+    handle->workspace->H_app_inv.allocated = FALSE;
+    handle->workspace->L.allocated = FALSE;
+    handle->workspace->U.allocated = FALSE;
+
+    for ( i = 0; i < LIST_N; ++i )
+    {
+        handle->lists[i]->allocated = FALSE;
+    }
+}
+
+
 /* Allocate top-level data structures and parse input files
  * for the first simulation
  *
@@ -161,58 +225,91 @@ static void Read_Input_Files( const char * const geo_file,
 void * setup( const char * const geo_file, const char * const ffield_file,
         const char * const control_file )
 {
-    int i;
     spuremd_handle *spmd_handle;
 
-    /* top-level allocation */
-    spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
-            "setup::spmd_handle" );
-
-    /* second-level allocations */
-    spmd_handle->system = smalloc( sizeof(reax_system),
-           "Setup::spmd_handle->system" );
-    spmd_handle->system->prealloc_allocated = FALSE;
-    spmd_handle->system->ffield_params_allocated = FALSE;
-    spmd_handle->system->g.allocated = FALSE;
-
-    spmd_handle->control = smalloc( sizeof(control_params),
-           "Setup::spmd_handle->control" );
-
-    spmd_handle->data = smalloc( sizeof(simulation_data),
-           "Setup::spmd_handle->data" );
-
-    spmd_handle->workspace = smalloc( sizeof(static_storage),
-           "Setup::spmd_handle->workspace" );
-    spmd_handle->workspace->H.allocated = FALSE;
-    spmd_handle->workspace->H_full.allocated = FALSE;
-    spmd_handle->workspace->H_sp.allocated = FALSE;
-    spmd_handle->workspace->H_p.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
-    spmd_handle->workspace->H_app_inv.allocated = FALSE;
-    spmd_handle->workspace->L.allocated = FALSE;
-    spmd_handle->workspace->U.allocated = FALSE;
-
-    spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
-           "Setup::spmd_handle->lists" );
-    for ( i = 0; i < LIST_N; ++i )
-    {
-        spmd_handle->lists[i] = smalloc( sizeof(reax_list),
-                "Setup::spmd_handle->lists[i]" );
-        spmd_handle->lists[i]->allocated = FALSE;
-    }
-    spmd_handle->out_control = smalloc( sizeof(output_controls),
-           "Setup::spmd_handle->out_control" );
+    Allocate_Top_Level_Structs( &spmd_handle );
+    Initialize_Top_Level_Structs( spmd_handle );
 
-    spmd_handle->output_enabled = TRUE;
-    spmd_handle->realloc = TRUE;
-    spmd_handle->callback = NULL;
-    spmd_handle->data->sim_id = 0;
+    /* note: assign here to avoid compiler warning
+     * of uninitialized usage in PreAllocate_Space */
+    spmd_handle->system->N_max = 0;
 
     Read_Input_Files( geo_file, ffield_file, control_file,
             spmd_handle->system, spmd_handle->control,
             spmd_handle->data, spmd_handle->workspace,
-            spmd_handle->out_control );
+            spmd_handle->out_control, FALSE );
+
+    spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+
+    return (void *) spmd_handle;
+}
+
+
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * num_atoms: num. atoms in this simulation
+ * types: integer representation of atom element (type)
+ *  NOTE: must match the 0-based index from section 2 in the ReaxFF parameter file
+ * sim_box_info: simulation box information, where the entries are
+ *  - box length per dimension (3 entries)
+ *  - angles per dimension (3 entries)
+ * pos: coordinates of atom positions (consecutively arranged), in Angstroms
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+void * setup2( int num_atoms, const int * const atom_type,
+        const double * const pos, const double * const sim_box_info,
+        const char * const ffield_file, const char * const control_file )
+{
+    int i;
+//    char atom_name[9];
+    rvec x;
+    spuremd_handle *spmd_handle;
+
+    Allocate_Top_Level_Structs( &spmd_handle );
+    Initialize_Top_Level_Structs( spmd_handle );
+
+    /* override default */
+    spmd_handle->output_enabled = FALSE;
+
+    Read_Input_Files( NULL, ffield_file, control_file,
+            spmd_handle->system, spmd_handle->control,
+            spmd_handle->data, spmd_handle->workspace,
+            spmd_handle->out_control, FALSE );
+
+    spmd_handle->system->N = num_atoms;
+    /* note: assign here to avoid compiler warning
+     * of uninitialized usage in PreAllocate_Space */
+    spmd_handle->system->N_max = 0;
+
+    PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+            spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
+
+    Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
+            sim_box_info[3], sim_box_info[4], sim_box_info[5],
+            &spmd_handle->system->box );
+
+    for ( i = 0; i < spmd_handle->system->N; ++i )
+    {
+        x[0] = pos[3 * i];
+        x[1] = pos[3 * i + 1];
+        x[2] = pos[3 * i + 2];
+
+        Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+        spmd_handle->workspace->orig_id[i] = i + 1;
+//        spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+//                element, sizeof(element) );
+        spmd_handle->system->atoms[i].type = atom_type[i];
+//        strncpy( spmd_handle->system->atoms[i].name, atom_name,
+//                sizeof(spmd_handle->system->atoms[i].name) - 1 );
+//        spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+        rvec_Copy( spmd_handle->system->atoms[i].x, x );
+        rvec_MakeZero( spmd_handle->system->atoms[i].v );
+        rvec_MakeZero( spmd_handle->system->atoms[i].f );
+        spmd_handle->system->atoms[i].q = 0.0;
+    }
 
     spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
 
@@ -293,17 +390,7 @@ int simulate( const void * const handle )
                     spmd_handle->data, spmd_handle->out_control );
         }
 
-        if ( spmd_handle->output_enabled == TRUE || spmd_handle->callback != NULL )
-        {
-            if ( ((spmd_handle->out_control->log_update_freq > 0
-                        && spmd_handle->data->step % spmd_handle->out_control->log_update_freq == 0)
-                    || (spmd_handle->out_control->write_steps > 0
-                        && spmd_handle->data->step % spmd_handle->out_control->write_steps == 0))
-                || spmd_handle->callback != NULL )
-            {
-                Compute_Total_Energy( spmd_handle->data );
-            }
-        }
+        Compute_Total_Energy( spmd_handle->data );
 
         if ( spmd_handle->output_enabled == TRUE )
         {
@@ -473,7 +560,7 @@ int reset( const void * const handle, const char * const geo_file,
         Read_Input_Files( geo_file, ffield_file, control_file,
                 spmd_handle->system, spmd_handle->control,
                 spmd_handle->data, spmd_handle->workspace,
-                spmd_handle->out_control );
+                spmd_handle->out_control, TRUE );
 
         if ( spmd_handle->system->N > spmd_handle->system->N_max )
         {
@@ -495,6 +582,104 @@ int reset( const void * const handle, const char * const geo_file,
 }
 
 
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * num_atoms: num. atoms in this simulation
+ * types: integer representation of atom element (type)
+ *  NOTE: must match the 0-based index from section 2 in the ReaxFF parameter file
+ * sim_box_info: simulation box information, where the entries are
+ *  - box length per dimension (3 entries)
+ *  - angles per dimension (3 entries)
+ * pos: coordinates of atom positions (consecutively arranged), in Angstroms
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+int reset2( const void * const handle, int num_atoms,
+        const int * const atom_type, const double * const pos,
+        const double * const sim_box_info, const char * const ffield_file,
+        const char * const control_file )
+{
+    int i, ret;
+    rvec x;
+    spuremd_handle *spmd_handle;
+
+    ret = SPUREMD_FAILURE;
+
+    if ( handle != NULL )
+    {
+        spmd_handle = (spuremd_handle*) handle;
+
+        /* close files used in previous simulation */
+        if ( spmd_handle->output_enabled == TRUE )
+        {
+            Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
+                    spmd_handle->workspace, spmd_handle->out_control );
+        }
+
+        spmd_handle->realloc = FALSE;
+        spmd_handle->data->sim_id++;
+
+        Read_Input_Files( NULL, ffield_file, control_file,
+                spmd_handle->system, spmd_handle->control,
+                spmd_handle->data, spmd_handle->workspace,
+                spmd_handle->out_control, TRUE );
+
+        spmd_handle->system->N = num_atoms;
+
+        if ( spmd_handle->system->prealloc_allocated == FALSE
+                || spmd_handle->system->N > spmd_handle->system->N_max )
+        {
+            PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+                    spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
+        }
+
+        Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
+                sim_box_info[3], sim_box_info[4], sim_box_info[5],
+                &spmd_handle->system->box );
+
+        for ( i = 0; i < spmd_handle->system->N; ++i )
+        {
+            x[0] = pos[3 * i];
+            x[1] = pos[3 * i + 1];
+            x[2] = pos[3 * i + 2];
+
+            Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+            spmd_handle->workspace->orig_id[i] = i + 1;
+//            spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+//                    element, sizeof(element) );
+            spmd_handle->system->atoms[i].type = atom_type[i];
+//            strncpy( spmd_handle->system->atoms[i].name, atom_name,
+//                    sizeof(spmd_handle->system->atoms[i].name) - 1 );
+//            spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+            rvec_Copy( spmd_handle->system->atoms[i].x, x );
+            rvec_MakeZero( spmd_handle->system->atoms[i].v );
+            rvec_MakeZero( spmd_handle->system->atoms[i].f );
+            spmd_handle->system->atoms[i].q = 0.0;
+        }
+
+        if ( spmd_handle->system->N > spmd_handle->system->N_max )
+        {
+            /* deallocate everything which needs more space
+             * (i.e., structures whose space is a function of the number of atoms),
+             * except for data structures allocated while parsing input files */
+            Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+                    spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+                    spmd_handle->output_enabled, TRUE );
+
+            spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+            spmd_handle->realloc = TRUE;
+        }
+
+        ret = SPUREMD_SUCCESS;
+    }
+
+    return ret;    
+}
+
+
 /* Getter for atom positions
  *
  * handle: pointer to wrapper struct with top-level data structures
@@ -684,6 +869,28 @@ int get_system_info( const void * const handle, double * const e_pot,
 }
 
 
+/* Getter for total energy
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * e_tot: system total energy, in kcal / mol (reference from caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_total_energy( const void * const handle, double * const e_tot )
+{
+    int ret;
+
+    ret = get_system_info( handle, e_tot, NULL, NULL, NULL, NULL, NULL );
+
+    if ( ret == SUCCESS )
+    {
+        ret = SPUREMD_SUCCESS;
+    }
+
+    return ret;
+}
+
+
 /* Setter for writing output to files
  *
  * handle: pointer to wrapper struct with top-level data structures
@@ -768,57 +975,26 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
     rvec x;
     spuremd_handle *spmd_handle;
 
-    /* top-level allocation */
-    spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
-            "setup::spmd_handle" );
-
-    /* second-level allocations */
-    spmd_handle->system = smalloc( sizeof(reax_system),
-           "Setup::spmd_handle->system" );
-    spmd_handle->system->prealloc_allocated = FALSE;
-    spmd_handle->system->ffield_params_allocated = FALSE;
-    spmd_handle->system->g.allocated = FALSE;
-
-    spmd_handle->control = smalloc( sizeof(control_params),
-           "Setup::spmd_handle->control" );
-
-    spmd_handle->data = smalloc( sizeof(simulation_data),
-           "Setup::spmd_handle->data" );
-
-    spmd_handle->workspace = smalloc( sizeof(static_storage),
-           "Setup::spmd_handle->workspace" );
-    spmd_handle->workspace->H.allocated = FALSE;
-    spmd_handle->workspace->H_full.allocated = FALSE;
-    spmd_handle->workspace->H_sp.allocated = FALSE;
-    spmd_handle->workspace->H_p.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
-    spmd_handle->workspace->H_app_inv.allocated = FALSE;
-    spmd_handle->workspace->L.allocated = FALSE;
-    spmd_handle->workspace->U.allocated = FALSE;
-
-    spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
-           "Setup::spmd_handle->lists" );
-    for ( i = 0; i < LIST_N; ++i )
-    {
-        spmd_handle->lists[i] = smalloc( sizeof(reax_list),
-                "Setup::spmd_handle->lists[i]" );
-        spmd_handle->lists[i]->allocated = FALSE;
-    }
-    spmd_handle->out_control = smalloc( sizeof(output_controls),
-           "Setup::spmd_handle->out_control" );
+    Allocate_Top_Level_Structs( &spmd_handle );
+    Initialize_Top_Level_Structs( spmd_handle );
 
+    /* override default */
     spmd_handle->output_enabled = FALSE;
-    spmd_handle->realloc = TRUE;
-    spmd_handle->callback = NULL;
-    spmd_handle->data->sim_id = 0;
+
+    Read_Input_Files( NULL, ffield_file, control_file,
+            spmd_handle->system, spmd_handle->control,
+            spmd_handle->data, spmd_handle->workspace,
+            spmd_handle->out_control, FALSE );
 
     spmd_handle->system->N_qm = qm_num_atoms;
     spmd_handle->system->N_mm = mm_num_atoms;
     spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
+    /* note: assign here to avoid compiler warning
+     * of uninitialized usage in PreAllocate_Space */
+    spmd_handle->system->N_max = 0;
 
     PreAllocate_Space( spmd_handle->system, spmd_handle->control,
-            spmd_handle->workspace, spmd_handle->system->N );
+            spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
 
     Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
             sim_box_info[3], sim_box_info[4], sim_box_info[5],
@@ -872,11 +1048,6 @@ void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
         spmd_handle->system->atoms[i].qmmm_mask = FALSE;
     }
 
-    Read_Input_Files( NULL, ffield_file, control_file,
-            spmd_handle->system, spmd_handle->control,
-            spmd_handle->data, spmd_handle->workspace,
-            spmd_handle->out_control );
-
     spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
 
     return (void *) spmd_handle;
@@ -930,12 +1101,21 @@ int reset_qmmm( const void * const handle,
         spmd_handle->realloc = FALSE;
         spmd_handle->data->sim_id++;
 
+        Read_Input_Files( NULL, ffield_file, control_file,
+                spmd_handle->system, spmd_handle->control,
+                spmd_handle->data, spmd_handle->workspace,
+                spmd_handle->out_control, TRUE );
+
         spmd_handle->system->N_qm = qm_num_atoms;
         spmd_handle->system->N_mm = mm_num_atoms;
         spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
 
-        PreAllocate_Space( spmd_handle->system, spmd_handle->control,
-                spmd_handle->workspace, spmd_handle->system->N );
+        if ( spmd_handle->system->prealloc_allocated == FALSE
+                || spmd_handle->system->N > spmd_handle->system->N_max )
+        {
+            PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+                    spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
+        }
 
         Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
                 sim_box_info[3], sim_box_info[4], sim_box_info[5],
@@ -987,11 +1167,6 @@ int reset_qmmm( const void * const handle,
             spmd_handle->system->atoms[i].qmmm_mask = FALSE;
         }
 
-        Read_Input_Files( NULL, ffield_file, control_file,
-                spmd_handle->system, spmd_handle->control,
-                spmd_handle->data, spmd_handle->workspace,
-                spmd_handle->out_control );
-
         if ( spmd_handle->system->N > spmd_handle->system->N_max )
         {
             /* deallocate everything which needs more space
@@ -1225,57 +1400,26 @@ void setup_qmmm_( void ** handle, const int * const qm_num_atoms,
     rvec x;
     spuremd_handle *spmd_handle;
 
-    /* top-level allocation */
-    spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
-            "setup::spmd_handle" );
-
-    /* second-level allocations */
-    spmd_handle->system = smalloc( sizeof(reax_system),
-           "Setup::spmd_handle->system" );
-    spmd_handle->system->prealloc_allocated = FALSE;
-    spmd_handle->system->ffield_params_allocated = FALSE;
-    spmd_handle->system->g.allocated = FALSE;
-
-    spmd_handle->control = smalloc( sizeof(control_params),
-           "Setup::spmd_handle->control" );
-
-    spmd_handle->data = smalloc( sizeof(simulation_data),
-           "Setup::spmd_handle->data" );
-
-    spmd_handle->workspace = smalloc( sizeof(static_storage),
-           "Setup::spmd_handle->workspace" );
-    spmd_handle->workspace->H.allocated = FALSE;
-    spmd_handle->workspace->H_full.allocated = FALSE;
-    spmd_handle->workspace->H_sp.allocated = FALSE;
-    spmd_handle->workspace->H_p.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt.allocated = FALSE;
-    spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
-    spmd_handle->workspace->H_app_inv.allocated = FALSE;
-    spmd_handle->workspace->L.allocated = FALSE;
-    spmd_handle->workspace->U.allocated = FALSE;
-
-    spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
-           "Setup::spmd_handle->lists" );
-    for ( i = 0; i < LIST_N; ++i )
-    {
-        spmd_handle->lists[i] = smalloc( sizeof(reax_list),
-                "Setup::spmd_handle->lists[i]" );
-        spmd_handle->lists[i]->allocated = FALSE;
-    }
-    spmd_handle->out_control = smalloc( sizeof(output_controls),
-           "Setup::spmd_handle->out_control" );
+    Allocate_Top_Level_Structs( &spmd_handle );
+    Initialize_Top_Level_Structs( spmd_handle );
 
+    /* override default */
     spmd_handle->output_enabled = FALSE;
-    spmd_handle->realloc = TRUE;
-    spmd_handle->callback = NULL;
-    spmd_handle->data->sim_id = 0;
+
+    Read_Input_Files( NULL, ffield_file, control_file,
+            spmd_handle->system, spmd_handle->control,
+            spmd_handle->data, spmd_handle->workspace,
+            spmd_handle->out_control, FALSE );
 
     spmd_handle->system->N_qm = *qm_num_atoms;
     spmd_handle->system->N_mm = *mm_num_atoms;
     spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
+    /* note: assign here to avoid compiler warning
+     * of uninitialized usage in PreAllocate_Space */
+    spmd_handle->system->N_max = 0;
 
     PreAllocate_Space( spmd_handle->system, spmd_handle->control,
-            spmd_handle->workspace, spmd_handle->system->N );
+            spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
 
     Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
             sim_box_info[3], sim_box_info[4], sim_box_info[5],
@@ -1329,11 +1473,6 @@ void setup_qmmm_( void ** handle, const int * const qm_num_atoms,
         spmd_handle->system->atoms[i].qmmm_mask = FALSE;
     }
 
-    Read_Input_Files( NULL, ffield_file, control_file,
-            spmd_handle->system, spmd_handle->control,
-            spmd_handle->data, spmd_handle->workspace,
-            spmd_handle->out_control );
-
     spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
 
     *handle = (void *) spmd_handle;
@@ -1383,12 +1522,21 @@ void reset_qmmm_( const void * const handle,
         spmd_handle->realloc = FALSE;
         spmd_handle->data->sim_id++;
 
+        Read_Input_Files( NULL, ffield_file, control_file,
+                spmd_handle->system, spmd_handle->control,
+                spmd_handle->data, spmd_handle->workspace,
+                spmd_handle->out_control, TRUE );
+
         spmd_handle->system->N_qm = *qm_num_atoms;
         spmd_handle->system->N_mm = *mm_num_atoms;
         spmd_handle->system->N = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
 
-        PreAllocate_Space( spmd_handle->system, spmd_handle->control,
-                spmd_handle->workspace, spmd_handle->system->N );
+        if ( spmd_handle->system->prealloc_allocated == FALSE
+                || spmd_handle->system->N > spmd_handle->system->N_max )
+        {
+            PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+                    spmd_handle->workspace, (int) CEIL( SAFE_ZONE * spmd_handle->system->N ) );
+        }
 
         Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
                 sim_box_info[3], sim_box_info[4], sim_box_info[5],
@@ -1440,11 +1588,6 @@ void reset_qmmm_( const void * const handle,
             spmd_handle->system->atoms[i].qmmm_mask = FALSE;
         }
 
-        Read_Input_Files( NULL, ffield_file, control_file,
-                spmd_handle->system, spmd_handle->control,
-                spmd_handle->data, spmd_handle->workspace,
-                spmd_handle->out_control );
-
         if ( spmd_handle->system->N > spmd_handle->system->N_max )
         {
             /* deallocate everything which needs more space
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index b0f53a8d2591ea646f023a123a277cacc68ee7ad..63d0ff5a5b49b40aab4a815b3a2f985e103e1519 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -36,6 +36,10 @@ extern "C"  {
 void * setup( const char * const, const char * const,
         const char * const );
 
+void * setup2( int, const int * const,
+        const double * const, const double * const,
+        const char * const, const char * const );
+
 int setup_callback( const void * const, const callback_function );
 
 int simulate( const void * const );
@@ -45,6 +49,10 @@ int cleanup( const void * const );
 int reset( const void * const, const char * const,
         const char * const, const char * const );
 
+int reset2( const void * const, int, const int * const,
+        const double * const, const double * const,
+        const char * const, const char * const );
+
 int get_atom_positions( const void * const, double * const );
 
 int get_atom_velocities( const void * const, double * const );
@@ -57,6 +65,8 @@ int get_system_info( const void * const, double * const,
         double * const, double * const, double * const,
         double * const, double * const );
 
+int get_total_energy( const void * const, double * const );
+
 int set_output_enabled( const void * const, const int );
 
 int set_control_parameter( const void * const, const char * const,
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 1766879c3e1bf454df4b602f82b0eb5af17192a0..45411f1bcb3a1582698196a00ebcbf32fcb940eb 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -154,9 +154,9 @@ void Make_Point( real x, real y, real z, rvec* p )
 }
 
 
-int is_Valid_Serial( static_storage *workspace, int serial )
+int is_Valid_Serial( int serial )
 {
-    if( workspace->map_serials[ serial ] < 0 )
+    if ( serial < 0 )
     {
         fprintf( stderr, "[ERROR] CONECT line includes invalid serial number %d.\n", serial );
         fprintf( stderr, "[ERROR] Please correct the input file. Terminating...\n" );
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index 2c232b1b729b00ac92af63f515590f8e94d46872..88e333c9a9a6b8d0a0a7a69240e0523177883c67 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -37,7 +37,7 @@ int is_Inside_Box( simulation_box *, rvec );
 /* from geo_tools.h */
 void Make_Point( real, real, real, rvec * );
 
-int is_Valid_Serial( static_storage *, int );
+int is_Valid_Serial( int );
 
 int Check_Input_Range( int, int, int, char * );