From 7ae444ccac5f3b7e6c7c953219ca898ea9509956 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 6 Dec 2017 11:04:15 -0500
Subject: [PATCH] sPuReMD: split into library and driver for optimization
 integration (NSGA-II). Fix bug with interaction lists being allocated for
 hydrogen bonds when no hydorgen atoms present. Fix minimum allocation size of
 3-body lists (uncovered with small lithium geometries). Refactor naming of
 conflicting data structures with NSGA-II. Update .gitignore.

---
 .gitignore                             |   6 +
 ar-lib                                 | 270 +++++++++++++++++++++++++
 sPuReMD/Makefile.am                    |  20 +-
 sPuReMD/aclocal.m4                     |  65 ++++++
 sPuReMD/configure.ac                   |   5 +-
 sPuReMD/src/allocate.c                 |  14 +-
 sPuReMD/src/allocate.h                 |   6 +-
 sPuReMD/src/analyze.c                  |  46 ++---
 sPuReMD/src/analyze.h                  |  12 +-
 sPuReMD/src/bond_orders.c              |  66 +++---
 sPuReMD/src/bond_orders.h              |  22 +-
 sPuReMD/src/charges.c                  |  22 +-
 sPuReMD/src/charges.h                  |   2 +-
 sPuReMD/src/control.c                  |   6 +-
 sPuReMD/src/driver.c                   |  53 +++++
 sPuReMD/src/forces.c                   |  30 +--
 sPuReMD/src/forces.h                   |   4 +-
 sPuReMD/src/four_body_interactions.c   |   4 +-
 sPuReMD/src/four_body_interactions.h   |   2 +-
 sPuReMD/src/geo_tools.c                |   2 +-
 sPuReMD/src/geo_tools.h                |   2 +-
 sPuReMD/src/init_md.c                  |  44 ++--
 sPuReMD/src/init_md.h                  |   5 +-
 sPuReMD/src/integrate.c                |  24 +--
 sPuReMD/src/integrate.h                |  16 +-
 sPuReMD/src/lin_alg.c                  |   2 +-
 sPuReMD/src/lin_alg.h                  |   2 +-
 sPuReMD/src/list.c                     |  14 +-
 sPuReMD/src/list.h                     |  14 +-
 sPuReMD/src/mytypes.h                  |  12 +-
 sPuReMD/src/neighbors.c                |  18 +-
 sPuReMD/src/neighbors.h                |   4 +-
 sPuReMD/src/print_utils.c              |  48 ++---
 sPuReMD/src/print_utils.h              |  34 ++--
 sPuReMD/src/reset_utils.c              |   8 +-
 sPuReMD/src/reset_utils.h              |   6 +-
 sPuReMD/src/single_body_interactions.c |   4 +-
 sPuReMD/src/single_body_interactions.h |   2 +-
 sPuReMD/src/{testmd.c => spuremd.c}    | 118 ++++++-----
 sPuReMD/src/spuremd.h                  |  38 ++++
 sPuReMD/src/three_body_interactions.c  |   8 +-
 sPuReMD/src/three_body_interactions.h  |   4 +-
 sPuReMD/src/tool_box.c                 |   4 +-
 sPuReMD/src/traj.c                     |   8 +-
 sPuReMD/src/traj.h                     |  10 +-
 sPuReMD/src/two_body_interactions.c    |  12 +-
 sPuReMD/src/two_body_interactions.h    |   6 +-
 47 files changed, 793 insertions(+), 331 deletions(-)
 create mode 100755 ar-lib
 create mode 100644 sPuReMD/src/driver.c
 rename sPuReMD/src/{testmd.c => spuremd.c} (62%)
 create mode 100644 sPuReMD/src/spuremd.h

diff --git a/.gitignore b/.gitignore
index a7a8c79..0e82e4b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -50,9 +50,13 @@ Makefile
 Makefile.in
 stamp-h1
 test-driver
+*.lo
+*.la
 
 # Compiled languages (C,C++,Fortran,...)
 *.o
+*.so
+*.a
 
 # flex/bison files
 *.yy.c
@@ -68,4 +72,6 @@ test-driver
 *.tar.gz
 *.pdf
 */bin
+*/lib
 *.txt
+*.zip
diff --git a/ar-lib b/ar-lib
new file mode 100755
index 0000000..463b9ec
--- /dev/null
+++ b/ar-lib
@@ -0,0 +1,270 @@
+#! /bin/sh
+# Wrapper for Microsoft lib.exe
+
+me=ar-lib
+scriptversion=2012-03-01.08; # UTC
+
+# Copyright (C) 2010-2014 Free Software Foundation, Inc.
+# Written by Peter Rosin <peda@lysator.liu.se>.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2, or (at your option)
+# any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# As a special exception to the GNU General Public License, if you
+# distribute this file as part of a program that contains a
+# configuration script generated by Autoconf, you may include it under
+# the same distribution terms that you use for the rest of that program.
+
+# This file is maintained in Automake, please report
+# bugs to <bug-automake@gnu.org> or send patches to
+# <automake-patches@gnu.org>.
+
+
+# func_error message
+func_error ()
+{
+  echo "$me: $1" 1>&2
+  exit 1
+}
+
+file_conv=
+
+# func_file_conv build_file
+# Convert a $build file to $host form and store it in $file
+# Currently only supports Windows hosts.
+func_file_conv ()
+{
+  file=$1
+  case $file in
+    / | /[!/]*) # absolute file, and not a UNC file
+      if test -z "$file_conv"; then
+	# lazily determine how to convert abs files
+	case `uname -s` in
+	  MINGW*)
+	    file_conv=mingw
+	    ;;
+	  CYGWIN*)
+	    file_conv=cygwin
+	    ;;
+	  *)
+	    file_conv=wine
+	    ;;
+	esac
+      fi
+      case $file_conv in
+	mingw)
+	  file=`cmd //C echo "$file " | sed -e 's/"\(.*\) " *$/\1/'`
+	  ;;
+	cygwin)
+	  file=`cygpath -m "$file" || echo "$file"`
+	  ;;
+	wine)
+	  file=`winepath -w "$file" || echo "$file"`
+	  ;;
+      esac
+      ;;
+  esac
+}
+
+# func_at_file at_file operation archive
+# Iterate over all members in AT_FILE performing OPERATION on ARCHIVE
+# for each of them.
+# When interpreting the content of the @FILE, do NOT use func_file_conv,
+# since the user would need to supply preconverted file names to
+# binutils ar, at least for MinGW.
+func_at_file ()
+{
+  operation=$2
+  archive=$3
+  at_file_contents=`cat "$1"`
+  eval set x "$at_file_contents"
+  shift
+
+  for member
+  do
+    $AR -NOLOGO $operation:"$member" "$archive" || exit $?
+  done
+}
+
+case $1 in
+  '')
+     func_error "no command.  Try '$0 --help' for more information."
+     ;;
+  -h | --h*)
+    cat <<EOF
+Usage: $me [--help] [--version] PROGRAM ACTION ARCHIVE [MEMBER...]
+
+Members may be specified in a file named with @FILE.
+EOF
+    exit $?
+    ;;
+  -v | --v*)
+    echo "$me, version $scriptversion"
+    exit $?
+    ;;
+esac
+
+if test $# -lt 3; then
+  func_error "you must specify a program, an action and an archive"
+fi
+
+AR=$1
+shift
+while :
+do
+  if test $# -lt 2; then
+    func_error "you must specify a program, an action and an archive"
+  fi
+  case $1 in
+    -lib | -LIB \
+    | -ltcg | -LTCG \
+    | -machine* | -MACHINE* \
+    | -subsystem* | -SUBSYSTEM* \
+    | -verbose | -VERBOSE \
+    | -wx* | -WX* )
+      AR="$AR $1"
+      shift
+      ;;
+    *)
+      action=$1
+      shift
+      break
+      ;;
+  esac
+done
+orig_archive=$1
+shift
+func_file_conv "$orig_archive"
+archive=$file
+
+# strip leading dash in $action
+action=${action#-}
+
+delete=
+extract=
+list=
+quick=
+replace=
+index=
+create=
+
+while test -n "$action"
+do
+  case $action in
+    d*) delete=yes  ;;
+    x*) extract=yes ;;
+    t*) list=yes    ;;
+    q*) quick=yes   ;;
+    r*) replace=yes ;;
+    s*) index=yes   ;;
+    S*)             ;; # the index is always updated implicitly
+    c*) create=yes  ;;
+    u*)             ;; # TODO: don't ignore the update modifier
+    v*)             ;; # TODO: don't ignore the verbose modifier
+    *)
+      func_error "unknown action specified"
+      ;;
+  esac
+  action=${action#?}
+done
+
+case $delete$extract$list$quick$replace,$index in
+  yes,* | ,yes)
+    ;;
+  yesyes*)
+    func_error "more than one action specified"
+    ;;
+  *)
+    func_error "no action specified"
+    ;;
+esac
+
+if test -n "$delete"; then
+  if test ! -f "$orig_archive"; then
+    func_error "archive not found"
+  fi
+  for member
+  do
+    case $1 in
+      @*)
+        func_at_file "${1#@}" -REMOVE "$archive"
+        ;;
+      *)
+        func_file_conv "$1"
+        $AR -NOLOGO -REMOVE:"$file" "$archive" || exit $?
+        ;;
+    esac
+  done
+
+elif test -n "$extract"; then
+  if test ! -f "$orig_archive"; then
+    func_error "archive not found"
+  fi
+  if test $# -gt 0; then
+    for member
+    do
+      case $1 in
+        @*)
+          func_at_file "${1#@}" -EXTRACT "$archive"
+          ;;
+        *)
+          func_file_conv "$1"
+          $AR -NOLOGO -EXTRACT:"$file" "$archive" || exit $?
+          ;;
+      esac
+    done
+  else
+    $AR -NOLOGO -LIST "$archive" | sed -e 's/\\/\\\\/g' | while read member
+    do
+      $AR -NOLOGO -EXTRACT:"$member" "$archive" || exit $?
+    done
+  fi
+
+elif test -n "$quick$replace"; then
+  if test ! -f "$orig_archive"; then
+    if test -z "$create"; then
+      echo "$me: creating $orig_archive"
+    fi
+    orig_archive=
+  else
+    orig_archive=$archive
+  fi
+
+  for member
+  do
+    case $1 in
+    @*)
+      func_file_conv "${1#@}"
+      set x "$@" "@$file"
+      ;;
+    *)
+      func_file_conv "$1"
+      set x "$@" "$file"
+      ;;
+    esac
+    shift
+    shift
+  done
+
+  if test -n "$orig_archive"; then
+    $AR -NOLOGO -OUT:"$archive" "$orig_archive" "$@" || exit $?
+  else
+    $AR -NOLOGO -OUT:"$archive" "$@" || exit $?
+  fi
+
+elif test -n "$list"; then
+  if test ! -f "$orig_archive"; then
+    func_error "archive not found"
+  fi
+  $AR -NOLOGO -LIST "$archive" || exit $?
+fi
diff --git a/sPuReMD/Makefile.am b/sPuReMD/Makefile.am
index a31309c..2531185 100644
--- a/sPuReMD/Makefile.am
+++ b/sPuReMD/Makefile.am
@@ -1,22 +1,20 @@
 ACLOCAL_AMFLAGS = -I ../m4
 
-bin_PROGRAMS = bin/spuremd
-
-bin_spuremd_SOURCES = src/ffield.c src/grid.c src/list.c src/lookup.c src/print_utils.c \
+lib_LTLIBRARIES = lib/libspuremd.la
+lib_libspuremd_la_SOURCES = src/ffield.c src/grid.c src/list.c src/lookup.c src/print_utils.c \
 		  src/reset_utils.c src/restart.c src/random.c src/tool_box.c src/traj.c \
 		  src/vector.c src/allocate.c src/analyze.c src/box.c src/system_props.c src/control.c \
 		  src/geo_tools.c src/neighbors.c src/lin_alg.c src/charges.c src/bond_orders.c \
 		  src/single_body_interactions.c src/two_body_interactions.c \
 		  src/three_body_interactions.c src/four_body_interactions.c src/forces.c \
-		  src/integrate.c src/init_md.c src/testmd.c 
+		  src/integrate.c src/init_md.c src/spuremd.c
+lib_libspuremd_la_CFLAGS = -I src
+lib_libspuremd_la_LDFLAGS = -version-info 1:0:0
 
-include_HEADERS = src/mytypes.h src/ffield.h src/grid.h src/list.h src/lookup.h src/print_utils.h \
-		  src/reset_utils.h src/restart.h src/random.h src/tool_box.h src/traj.h \
-		  src/vector.h src/allocate.h src/analyze.h src/box.h src/system_props.h src/control.h \
-		  src/geo_tools.h src/neighbors.h src/lin_alg.h src/charges.h src/bond_orders.h \
-		  src/single_body_interactions.h src/two_body_interactions.h \
-		  src/three_body_interactions.h src/four_body_interactions.h src/forces.h \
-		  src/integrate.h src/init_md.h
+bin_PROGRAMS = bin/spuremd
+bin_spuremd_SOURCES = src/driver.c
+bin_spuremd_CPPFLAGS = -I src
+bin_spuremd_LDADD = lib/libspuremd.la
 
 check_PROGRAMS =
 TESTS =
diff --git a/sPuReMD/aclocal.m4 b/sPuReMD/aclocal.m4
index 48e9bea..df75c6a 100644
--- a/sPuReMD/aclocal.m4
+++ b/sPuReMD/aclocal.m4
@@ -56,6 +56,66 @@ m4_ifndef([AC_AUTOCONF_VERSION],
   [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
 _AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
 
+# Copyright (C) 2011-2014 Free Software Foundation, Inc.
+#
+# This file is free software; the Free Software Foundation
+# gives unlimited permission to copy and/or distribute it,
+# with or without modifications, as long as this notice is preserved.
+
+# AM_PROG_AR([ACT-IF-FAIL])
+# -------------------------
+# Try to determine the archiver interface, and trigger the ar-lib wrapper
+# if it is needed.  If the detection of archiver interface fails, run
+# ACT-IF-FAIL (default is to abort configure with a proper error message).
+AC_DEFUN([AM_PROG_AR],
+[AC_BEFORE([$0], [LT_INIT])dnl
+AC_BEFORE([$0], [AC_PROG_LIBTOOL])dnl
+AC_REQUIRE([AM_AUX_DIR_EXPAND])dnl
+AC_REQUIRE_AUX_FILE([ar-lib])dnl
+AC_CHECK_TOOLS([AR], [ar lib "link -lib"], [false])
+: ${AR=ar}
+
+AC_CACHE_CHECK([the archiver ($AR) interface], [am_cv_ar_interface],
+  [AC_LANG_PUSH([C])
+   am_cv_ar_interface=ar
+   AC_COMPILE_IFELSE([AC_LANG_SOURCE([[int some_variable = 0;]])],
+     [am_ar_try='$AR cru libconftest.a conftest.$ac_objext >&AS_MESSAGE_LOG_FD'
+      AC_TRY_EVAL([am_ar_try])
+      if test "$ac_status" -eq 0; then
+        am_cv_ar_interface=ar
+      else
+        am_ar_try='$AR -NOLOGO -OUT:conftest.lib conftest.$ac_objext >&AS_MESSAGE_LOG_FD'
+        AC_TRY_EVAL([am_ar_try])
+        if test "$ac_status" -eq 0; then
+          am_cv_ar_interface=lib
+        else
+          am_cv_ar_interface=unknown
+        fi
+      fi
+      rm -f conftest.lib libconftest.a
+     ])
+   AC_LANG_POP([C])])
+
+case $am_cv_ar_interface in
+ar)
+  ;;
+lib)
+  # Microsoft lib, so override with the ar-lib wrapper script.
+  # FIXME: It is wrong to rewrite AR.
+  # But if we don't then we get into trouble of one sort or another.
+  # A longer-term fix would be to have automake use am__AR in this case,
+  # and then we could set am__AR="$am_aux_dir/ar-lib \$(AR)" or something
+  # similar.
+  AR="$am_aux_dir/ar-lib $AR"
+  ;;
+unknown)
+  m4_default([$1],
+             [AC_MSG_ERROR([could not determine $AR interface])])
+  ;;
+esac
+AC_SUBST([AR])dnl
+])
+
 # AM_AUX_DIR_EXPAND                                         -*- Autoconf -*-
 
 # Copyright (C) 2001-2014 Free Software Foundation, Inc.
@@ -1152,3 +1212,8 @@ AC_SUBST([am__untar])
 
 m4_include([../m4/acx_pthread.m4])
 m4_include([../m4/ax_compiler_vendor.m4])
+m4_include([../m4/libtool.m4])
+m4_include([../m4/ltoptions.m4])
+m4_include([../m4/ltsugar.m4])
+m4_include([../m4/ltversion.m4])
+m4_include([../m4/lt~obsolete.m4])
diff --git a/sPuReMD/configure.ac b/sPuReMD/configure.ac
index 312fa80..1a59944 100644
--- a/sPuReMD/configure.ac
+++ b/sPuReMD/configure.ac
@@ -8,8 +8,9 @@ AC_INIT([sPuReMD], [1.0], [ohearnku@msu.edu hma@msu.edu])
 AM_INIT_AUTOMAKE([1.15 subdir-objects -Wall -Werror foreign])
 # Enable silent build rules by default.
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])], [AC_SUBST([AM_DEFAULT_VERBOSITY],[1])])
-#LT_PREREQ([2.2])
-#LT_INIT([dlopen])
+AM_PROG_AR
+LT_PREREQ([2.2])
+LT_INIT([dlopen])
 
 AC_LANG([C])
 
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 8035410..15f4256 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -55,7 +55,7 @@ int PreAllocate_Space( reax_system *system, control_params *control,
 }
 
 
-void Reallocate_Neighbor_List( list *far_nbrs, int n, int num_intrs )
+void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs )
 {
     Delete_List( TYP_FAR_NEIGHBOR, far_nbrs );
 
@@ -131,7 +131,7 @@ int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name )
 
 
 int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
-        list *hbonds )
+        reax_list *hbonds )
 {
     int i, num_hbonds;
 
@@ -168,7 +168,7 @@ int Allocate_HBond_List( int n, int num_h, int *h_index, int *hb_top,
 }
 
 
-int Reallocate_HBonds_List(  int n, int num_h, int *h_index, list *hbonds )
+int Reallocate_HBonds_List( int n, int num_h, int *h_index, reax_list *hbonds )
 {
     int i;
     int *hb_top;
@@ -196,7 +196,7 @@ int Reallocate_HBonds_List(  int n, int num_h, int *h_index, list *hbonds )
 }
 
 
-int Allocate_Bond_List( int n, int *bond_top, list *bonds )
+int Allocate_Bond_List( int n, int *bond_top, reax_list *bonds )
 {
     int i, num_bonds;
 
@@ -227,7 +227,7 @@ int Allocate_Bond_List( int n, int *bond_top, list *bonds )
 }
 
 
-int Reallocate_Bonds_List( int n, list *bonds, int *num_bonds, int *est_3body )
+int Reallocate_Bonds_List( int n, reax_list *bonds, int *num_bonds, int *est_3body )
 {
     int i;
     int *bond_top;
@@ -256,7 +256,7 @@ int Reallocate_Bonds_List( int n, list *bonds, int *num_bonds, int *est_3body )
 }
 
 
-void Reallocate( reax_system *system, static_storage *workspace, list **lists,
+void Reallocate( reax_system *system, control_params *control, static_storage *workspace, reax_list **lists,
         int nbr_flag )
 {
     int i, j, k;
@@ -285,7 +285,7 @@ void Reallocate( reax_system *system, static_storage *workspace, list **lists,
         workspace->U = NULL;
     }
 
-    if ( realloc->hbonds > 0 )
+    if ( control->hb_cut > 0.0 && realloc->hbonds > 0 )
     {
         Reallocate_HBonds_List(system->N, workspace->num_H, workspace->hbond_index,
                                (*lists) + HBONDS );
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index e8fdc4d..b961521 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -26,14 +26,14 @@
 
 int PreAllocate_Space( reax_system*, control_params*, static_storage* );
 
-void Reallocate( reax_system*, static_storage*, list**, int );
+void Reallocate( reax_system*, control_params*, static_storage*, reax_list**, int );
 
 int Allocate_Matrix( sparse_matrix**, int, int );
 
 void Deallocate_Matrix( sparse_matrix* );
 
-int Allocate_HBond_List( int, int, int*, int*, list* );
+int Allocate_HBond_List( int, int, int*, int*, reax_list* );
 
-int Allocate_Bond_List( int, int*, list* );
+int Allocate_Bond_List( int, int*, reax_list* );
 
 #endif
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index d3582fe..4e2cbe2 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -59,11 +59,11 @@ typedef struct
 
 // copy bond list into old bond list
 void Copy_Bond_List( reax_system *system, control_params *control,
-                     list **lists )
+                     reax_list **lists )
 {
     int i, j, top_old;
-    list *new_bonds = (*lists) + BONDS;
-    list *old_bonds = (*lists) + OLD_BONDS;
+    reax_list *new_bonds = (*lists) + BONDS;
+    reax_list *old_bonds = (*lists) + OLD_BONDS;
 
     for ( top_old = 0, i = 0; i < system->N; ++i )
     {
@@ -89,11 +89,11 @@ void Copy_Bond_List( reax_system *system, control_params *control,
 
 
 // ASSUMPTION: Bond lists are sorted
-int Compare_Bond_Lists( int atom, control_params *control, list **lists )
+int Compare_Bond_Lists( int atom, control_params *control, reax_list **lists )
 {
     int oldp, newp;
-    list *new_bonds = (*lists) + BONDS;
-    list *old_bonds = (*lists) + OLD_BONDS;
+    reax_list *new_bonds = (*lists) + BONDS;
+    reax_list *old_bonds = (*lists) + OLD_BONDS;
 
     /*fprintf( stdout, "\n%d\nold_bonds:", atom );
       for( oldp = Start_Index( atom, old_bonds );
@@ -168,7 +168,7 @@ int Compare_Bond_Lists( int atom, control_params *control, list **lists )
 
 
 void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
-                   control_params *control, list *bonds, int print,
+                   control_params *control, reax_list *bonds, int print,
                    FILE *fout )
 {
     int i, start, end;
@@ -219,15 +219,15 @@ void Print_Molecule( reax_system *system, molecule *m, int mode, char *s )
 
 void Analyze_Molecules( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, FILE *fout )
+        reax_list **lists, FILE *fout )
 {
     int atom, i, j, k, l, flag;
     int *mark = workspace->mark;
     int *old_mark = workspace->old_mark;
     int num_old, num_new;
     char s[MAX_MOLECULE_SIZE * 10];
-    list *new_bonds = (*lists) + BONDS;
-    list *old_bonds = (*lists) + OLD_BONDS;
+    reax_list *new_bonds = (*lists) + BONDS;
+    reax_list *old_bonds = (*lists) + OLD_BONDS;
     molecule old_molecules[20], new_molecules[20];
 
     fprintf( fout, "molecular analysis @ %d\n", data->step );
@@ -314,8 +314,8 @@ void Analyze_Molecules( reax_system *system, control_params *control,
 
 
 void Report_Bond_Change( reax_system *system, control_params *control,
-                         static_storage *workspace,  list *old_bonds,
-                         list *new_bonds, int a1, int a2, int flag,
+                         static_storage *workspace,  reax_list *old_bonds,
+                         reax_list *new_bonds, int a1, int a2, int flag,
                          FILE *fout )
 {
     int i;
@@ -385,8 +385,8 @@ void Report_Bond_Change( reax_system *system, control_params *control,
 
 /* ASSUMPTION: Bond lists are sorted */
 void Compare_Bonding( int atom, reax_system *system, control_params *control,
-                      static_storage *workspace, list *old_bonds,
-                      list *new_bonds, FILE *fout )
+                      static_storage *workspace, reax_list *old_bonds,
+                      reax_list *new_bonds, FILE *fout )
 {
     int oldp, newp;
 
@@ -503,7 +503,7 @@ void Compare_Bonding( int atom, reax_system *system, control_params *control,
 
 
 void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
-                  control_params *control, list *bonds, int ignore )
+                  control_params *control, reax_list *bonds, int ignore )
 {
     int i, t, start, end, nbr;
     real bo;
@@ -529,7 +529,7 @@ void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
 
 void Analyze_Fragments( reax_system *system, control_params *control,
                         simulation_data *data, static_storage *workspace,
-                        list **lists, FILE *fout, int ignore )
+                        reax_list **lists, FILE *fout, int ignore )
 {
     int  atom, i, flag;
     int  *mark = workspace->mark;
@@ -538,8 +538,8 @@ void Analyze_Fragments( reax_system *system, control_params *control,
     char fragments[MAX_FRAGMENT_TYPES][MAX_ATOM_TYPES];
     int  fragment_count[MAX_FRAGMENT_TYPES];
     molecule m;
-    list *new_bonds = (*lists) + BONDS;
-    //list *old_bonds = (*lists) + OLD_BONDS;
+    reax_list *new_bonds = (*lists) + BONDS;
+    //reax_list *old_bonds = (*lists) + OLD_BONDS;
 
     /* fragment analysis */
     fprintf( fout, "step%d fragments\n", data->step );
@@ -594,14 +594,14 @@ void Analyze_Fragments( reax_system *system, control_params *control,
 
 void Analyze_Silica( reax_system *system, control_params *control,
                      simulation_data *data, static_storage *workspace,
-                     list **lists, FILE *fout )
+                     reax_list **lists, FILE *fout )
 {
     int atom, i, j, k, pi, pk, pk_j, newp, coord;
     int O_SI_O_count, SI_O_SI_count;
     int si_coord[10], ox_coord[10];
     real O_SI_O, SI_O_SI;
-    list *new_bonds = (*lists) + BONDS;
-    list *thb_intrs =  (*lists) + THREE_BODIES;
+    reax_list *new_bonds = (*lists) + BONDS;
+    reax_list *thb_intrs =  (*lists) + THREE_BODIES;
 
     Analyze_Fragments( system, control, data, workspace, lists, fout, 0 );
 
@@ -721,7 +721,7 @@ int Get_Type_of_Molecule( molecule *m )
 
 void Calculate_Dipole_Moment( reax_system *system, control_params *control,
                               simulation_data *data, static_storage *workspace,
-                              list *bonds, FILE *fout )
+                              reax_list *bonds, FILE *fout )
 {
     int i, atom, count;
     molecule m;
@@ -936,7 +936,7 @@ void Calculate_Density_Slice( reax_system *system, simulation_data *data,
 
 void Analysis( reax_system *system, control_params *control,
                simulation_data *data, static_storage *workspace,
-               list **lists, output_controls *out_control )
+               reax_list **lists, output_controls *out_control )
 {
     int steps;
 
diff --git a/sPuReMD/src/analyze.h b/sPuReMD/src/analyze.h
index e296e57..af2b056 100644
--- a/sPuReMD/src/analyze.h
+++ b/sPuReMD/src/analyze.h
@@ -25,21 +25,21 @@
 #include "mytypes.h"
 
 void Analysis( reax_system*, control_params*, simulation_data*,
-               static_storage*, list**, output_controls* );
+               static_storage*, reax_list**, output_controls* );
 
-//void Copy_Bond_List( reax_system*, control_params*, list** );
+//void Copy_Bond_List( reax_system*, control_params*, reax_list** );
 
 //void Analyze_Molecules( reax_system*, control_params*, simulation_data*,
-//                        static_storage*, list**, FILE* );
+//                        static_storage*, reax_list**, FILE* );
 
 //void Analyze_Bonding( reax_system*, control_params*, simulation_data*,
-//                      static_storage*, list**, FILE* );
+//                      static_storage*, reax_list**, FILE* );
 
 //void Analyze_Silica( reax_system*, control_params*, simulation_data*,
-//                     static_storage*, list**, FILE* );
+//                     static_storage*, reax_list**, FILE* );
 
 //void Calculate_Dipole_Moment( reax_system*, control_params*, simulation_data*,
-//                              static_storage *, list*, FILE* );
+//                              static_storage *, reax_list*, FILE* );
 
 //void Copy_Positions( reax_system*, static_storage* );
 
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index 7a0358e..52d5ad6 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -34,11 +34,11 @@ static inline real Cf45( real p1, real p2 )
 
 
 #ifdef TEST_FORCES
-void Get_dBO( reax_system *system, list **lists,
+void Get_dBO( reax_system *system, reax_list **lists,
         int i, int pj, real C, rvec *v )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     int start_pj, end_pj, k;
 
     pj = bonds->select.bond_list[pj].dbond_index;
@@ -53,11 +53,11 @@ void Get_dBO( reax_system *system, list **lists,
 }
 
 
-void Get_dBOpinpi2( reax_system *system, list **lists,
+void Get_dBOpinpi2( reax_system *system, reax_list **lists,
         int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -74,11 +74,11 @@ void Get_dBOpinpi2( reax_system *system, list **lists,
 }
 
 
-void Add_dBO( reax_system *system, list **lists,
+void Add_dBO( reax_system *system, reax_list **lists,
         int i, int pj, real C, rvec *v )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     int start_pj, end_pj, k;
 
     pj = bonds->select.bond_list[pj].dbond_index;
@@ -95,11 +95,11 @@ void Add_dBO( reax_system *system, list **lists,
 }
 
 
-void Add_dBOpinpi2( reax_system *system, list **lists,
+void Add_dBOpinpi2( reax_system *system, reax_list **lists,
         int i, int pj, real Cpi, real Cpi2, rvec *vpi, rvec *vpi2 )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -116,11 +116,11 @@ void Add_dBOpinpi2( reax_system *system, list **lists,
 }
 
 
-void Add_dBO_to_Forces( reax_system *system, list **lists,
+void Add_dBO_to_Forces( reax_system *system, reax_list **lists,
                         int i, int pj, real C )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     int start_pj, end_pj, k;
 
     pj = bonds->select.bond_list[pj].dbond_index;
@@ -135,11 +135,11 @@ void Add_dBO_to_Forces( reax_system *system, list **lists,
 }
 
 
-void Add_dBOpinpi2_to_Forces( reax_system *system, list **lists,
+void Add_dBOpinpi2_to_Forces( reax_system *system, reax_list **lists,
         int i, int pj, real Cpi, real Cpi2 )
 {
-    list *bonds = (*lists) + BONDS;
-    list *dBOs = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs = (*lists) + DBO;
     dbond_data *dbo_k;
     int start_pj, end_pj, k;
 
@@ -156,9 +156,9 @@ void Add_dBOpinpi2_to_Forces( reax_system *system, list **lists,
 }
 
 
-void Add_dDelta( reax_system *system, list **lists, int i, real C, rvec *v )
+void Add_dDelta( reax_system *system, reax_list **lists, int i, real C, rvec *v )
 {
-    list *dDeltas;
+    reax_list *dDeltas;
     int start, end, k;
 
     dDeltas = &((*lists)[DDELTA]);
@@ -173,9 +173,9 @@ void Add_dDelta( reax_system *system, list **lists, int i, real C, rvec *v )
 }
 
 
-void Add_dDelta_to_Forces( reax_system *system, list **lists, int i, real C )
+void Add_dDelta_to_Forces( reax_system *system, reax_list **lists, int i, real C )
 {
-    list *dDeltas;
+    reax_list *dDeltas;
     int start, end, k;
 
     dDeltas = &((*lists)[DDELTA]);
@@ -190,12 +190,12 @@ void Add_dDelta_to_Forces( reax_system *system, list **lists, int i, real C )
 }
 
 
-void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
+void Calculate_dBO( int i, int pj, static_storage *workspace, reax_list **lists,
         int *top )
 {
     int j, k, l, start_i, end_i, end_j;
     rvec dDeltap_self, dBOp;
-    list *bonds, *dBOs;
+    reax_list *bonds, *dBOs;
     bond_data *nbr_l, *nbr_k;
     bond_order_data *bo_ij;
     dbond_data *top_dbo;
@@ -331,9 +331,9 @@ void Calculate_dBO( int i, int pj, static_storage *workspace, list **lists,
 
 
 void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
-        simulation_data *data, static_storage *workspace, list **lists )
+        simulation_data *data, static_storage *workspace, reax_list **lists )
 {
-    list *bonds = (*lists) + BONDS;
+    reax_list *bonds = (*lists) + BONDS;
     bond_data *nbr_j, *nbr_k;
     bond_order_data *bo_ij, *bo_ji;
     dbond_coefficients coef;
@@ -521,9 +521,9 @@ void Add_dBond_to_Forces_NPT( int i, int pj, reax_system *system,
 
 
 void Add_dBond_to_Forces( int i, int pj, reax_system *system,
-        simulation_data *data, static_storage *workspace, list **lists )
+        simulation_data *data, static_storage *workspace, reax_list **lists )
 {
-    list *bonds = (*lists) + BONDS;
+    reax_list *bonds = (*lists) + BONDS;
     bond_data *nbr_j, *nbr_k;
     bond_order_data *bo_ij, *bo_ji;
     dbond_coefficients coef;
@@ -657,7 +657,7 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
 /* Locate j on i's list.
    This function assumes that j is there for sure!
    And this is the case given our method of neighbor generation*/
-int Locate_Symmetric_Bond( list *bonds, int i, int j )
+int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
 {
     int start = Start_Index(i, bonds);
     int end = End_Index(i, bonds);
@@ -719,11 +719,11 @@ int compare_bonds( const void *p1, const void *p2 )
    This can either be done in the general coordinator function or here */
 void Calculate_Bond_Orders( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     real p_lp1;
     real p_boc1, p_boc2;
-    list *bonds;
+    reax_list *bonds;
 
     p_lp1 = system->reaxprm.gp.l[15];
     p_boc1 = system->reaxprm.gp.l[0];
@@ -753,8 +753,8 @@ void Calculate_Bond_Orders( reax_system *system, control_params *control,
         int top_dbo, top_dDelta;
         dbond_data *pdbo;
         dDelta_data *ptop_dDelta;
-        list *dDeltas;
-        list *dBOs;
+        reax_list *dDeltas;
+        reax_list *dBOs;
 
         top_dbo = 0;
         top_dDelta = 0;
diff --git a/sPuReMD/src/bond_orders.h b/sPuReMD/src/bond_orders.h
index b6ca199..67592ed 100644
--- a/sPuReMD/src/bond_orders.h
+++ b/sPuReMD/src/bond_orders.h
@@ -33,23 +33,23 @@ typedef struct
 } dbond_coefficients;
 
 #ifdef TEST_FORCES
-void Get_dBO( reax_system*, list**, int, int, real, rvec* );
-void Get_dBOpinpi2( reax_system*, list**, int, int, real, real, rvec*, rvec* );
+void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+void Get_dBOpinpi2( reax_system*, reax_list**, int, int, real, real, rvec*, rvec* );
 
-void Add_dBO( reax_system*, list**, int, int, real, rvec* );
-void Add_dBOpinpi2( reax_system*, list**, int, int, real, real, rvec*, rvec* );
+void Add_dBO( reax_system*, reax_list**, int, int, real, rvec* );
+void Add_dBOpinpi2( reax_system*, reax_list**, int, int, real, real, rvec*, rvec* );
 
-void Add_dBO_to_Forces( reax_system*, list**, int, int, real );
-void Add_dBOpinpi2_to_Forces( reax_system*, list**, int, int, real, real );
+void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, real );
+void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**, int, int, real, real );
 
-void Add_dDelta( reax_system*, list**, int, real, rvec* );
-void Add_dDelta_to_Forces( reax_system *, list**, int, real );
+void Add_dDelta( reax_system*, reax_list**, int, real, rvec* );
+void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real );
 #endif
 
 void Add_dBond_to_Forces( int, int, reax_system*, simulation_data*,
-                          static_storage*, list** );
+                          static_storage*, reax_list** );
 void Add_dBond_to_Forces_NPT( int, int, reax_system*, simulation_data*,
-                              static_storage*, list** );
+                              static_storage*, reax_list** );
 void Calculate_Bond_Orders( reax_system*, control_params*, simulation_data*,
-                            static_storage*, list**, output_controls* );
+                            static_storage*, reax_list**, output_controls* );
 #endif
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 63ed940..5003d28 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1514,7 +1514,7 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
 static void Compute_Preconditioner_QEq( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     real time;
     sparse_matrix *Hptr;
@@ -1612,7 +1612,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
 //static void Compute_Preconditioner_EE( const reax_system * const system,
 //        const control_params * const control,
 //        simulation_data * const data, static_storage * const workspace,
-//        const list * const far_nbrs )
+//        const reax_list * const far_nbrs )
 //{
 //    int i, top;
 //    static real * ones = NULL, * x = NULL, * y = NULL;
@@ -1856,7 +1856,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
 static void Compute_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     real time;
     sparse_matrix *Hptr;
@@ -1958,7 +1958,7 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
 static void Compute_Preconditioner_ACKS2( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     real time;
     sparse_matrix *Hptr;
@@ -2062,7 +2062,7 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
 static void Setup_Preconditioner_QEq( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     int fillin;
     real time;
@@ -2207,7 +2207,7 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
 static void Setup_Preconditioner_EE( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     int fillin;
     real time;
@@ -2356,7 +2356,7 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
 static void Setup_Preconditioner_ACKS2( const reax_system * const system,
         const control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs )
+        const reax_list * const far_nbrs )
 {
     int fillin;
     real time;
@@ -2561,7 +2561,7 @@ static void Calculate_Charges_EE( const reax_system * const system,
  */
 static void QEq( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs, const output_controls * const out_control )
+        const reax_list * const far_nbrs, const output_controls * const out_control )
 {
     int iters;
 
@@ -2650,7 +2650,7 @@ static void QEq( reax_system * const system, control_params * const control,
  */
 static void EE( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs, const output_controls * const out_control )
+        const reax_list * const far_nbrs, const output_controls * const out_control )
 {
     int iters;
 
@@ -2722,7 +2722,7 @@ static void EE( reax_system * const system, control_params * const control,
  */
 static void ACKS2( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs, const output_controls * const out_control )
+        const reax_list * const far_nbrs, const output_controls * const out_control )
 {
     int iters;
 
@@ -2784,7 +2784,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
 
 void Compute_Charges( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        const list * const far_nbrs, const output_controls * const out_control )
+        const reax_list * const far_nbrs, const output_controls * const out_control )
 {
 #if defined(DEBUG_FOCUS)
     char fname[200];
diff --git a/sPuReMD/src/charges.h b/sPuReMD/src/charges.h
index 50f563c..e0520fc 100644
--- a/sPuReMD/src/charges.h
+++ b/sPuReMD/src/charges.h
@@ -25,7 +25,7 @@
 #include "mytypes.h"
 
 void Compute_Charges( reax_system* const, control_params* const, simulation_data* const,
-          static_storage* const, const list* const,
+          static_storage* const, const reax_list* const,
           const output_controls* const );
 
 #endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 127543e..ead9035 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -115,7 +115,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
                   static_storage*, void* )) Write_Custom_Header;
     out_control->append_traj_frame =
         (int (*)( reax_system*, control_params*, simulation_data*,
-                  static_storage*, list **, void* )) Append_Custom_Frame;
+                  static_storage*, reax_list **, void* )) Append_Custom_Frame;
 
     strcpy( out_control->traj_title, "default_title" );
     out_control->atom_format = 0;
@@ -450,7 +450,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
                               static_storage*, void* )) Write_Custom_Header;
                 out_control->append_traj_frame =
                     (int (*)(reax_system*, control_params*, simulation_data*,
-                             static_storage*, list **, void*)) Append_Custom_Frame;
+                             static_storage*, reax_list **, void*)) Append_Custom_Frame;
             }
             else if ( out_control->traj_format == 1 )
             {
@@ -459,7 +459,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
                               static_storage*, void* )) Write_xyz_Header;
                 out_control->append_traj_frame =
                     (int (*)( reax_system*,  control_params*, simulation_data*,
-                              static_storage*, list **, void* )) Append_xyz_Frame;
+                              static_storage*, reax_list **, void* )) Append_xyz_Frame;
             }
         }
         else if ( strcmp(tmp[0], "traj_title") == 0 )
diff --git a/sPuReMD/src/driver.c b/sPuReMD/src/driver.c
new file mode 100644
index 0000000..6d2ff6a
--- /dev/null
+++ b/sPuReMD/src/driver.c
@@ -0,0 +1,53 @@
+
+/*----------------------------------------------------------------------
+  SerialReax - Reax Force Field Simulator
+
+  Copyright (2010) Purdue University
+  Hasan Metin Aktulga, haktulga@cs.purdue.edu
+  Joseph Fogarty, jcfogart@mail.usf.edu
+  Sagar Pandit, pandit@usf.edu
+  Ananth Y Grama, ayg@cs.purdue.edu
+
+  This program is free software; you can redistribute it and/or
+  modify it under the terms of the GNU General Public License as
+  published by the Free Software Foundation; either version 2 of
+  the License, or (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+  See the GNU General Public License for more details:
+  <http://www.gnu.org/licenses/>.
+  ----------------------------------------------------------------------*/
+
+#include "mytypes.h"
+
+#include "spuremd.h"
+
+
+static void usage( char * argv[] )
+{
+    fprintf( stderr, "usage: ./%s geometry ffield control\n", argv[0] );
+}
+
+
+int main( int argc, char* argv[] )
+{
+    reax_system system;
+    control_params control;
+    simulation_data data;
+
+    if ( argc != 4 )
+    {
+        usage( argv );
+        exit( INVALID_INPUT );
+    }
+
+    Setup( argv + 1, &system, &control, &data );
+
+    Run( &system, &control, &data, TRUE );
+
+    Cleanup( &system, &control, &data, TRUE );
+
+    return 0;
+}
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 4ae2acc..10e31fd 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -47,7 +47,7 @@ typedef enum
 
 void Dummy_Interaction( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
 }
 
@@ -77,7 +77,7 @@ void Init_Bonded_Force_Functions( control_params *control )
 
 void Compute_Bonded_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
 
     int i;
@@ -133,7 +133,7 @@ void Compute_Bonded_Forces( reax_system *system, control_params *control,
 
 void Compute_NonBonded_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list** lists, output_controls *out_control )
+        reax_list** lists, output_controls *out_control )
 {
     real t_start, t_elapsed;
 
@@ -177,10 +177,10 @@ void Compute_NonBonded_Forces( reax_system *system, control_params *control,
 /* This version of Compute_Total_Force computes forces from coefficients
    accumulated by all interaction functions. Saves enormous time & space! */
 void Compute_Total_Force( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists )
+        simulation_data *data, static_storage *workspace, reax_list **lists )
 {
     int i;
-    list *bonds;
+    reax_list *bonds;
 
     bonds = (*lists) + BONDS;
 
@@ -230,11 +230,11 @@ void Compute_Total_Force( reax_system *system, control_params *control,
 }
 
 
-void Validate_Lists( static_storage *workspace, list **lists, int step, int n,
+void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
         int Hmax, int Htop, int num_bonds, int num_hbonds )
 {
     int i, flag;
-    list *bonds, *hbonds;
+    reax_list *bonds, *hbonds;
 
     bonds = *lists + BONDS;
     hbonds = *lists + HBONDS;
@@ -439,7 +439,7 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
 
 
 static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
-        control_params *control, list *far_nbrs,
+        control_params *control, reax_list *far_nbrs,
         sparse_matrix * H, sparse_matrix * H_sp,
         int * Htop, int * H_sp_top )
 {
@@ -610,7 +610,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
 
 void Init_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, j, pj;
     int start_i, end_i;
@@ -623,7 +623,7 @@ void Init_Forces( reax_system *system, control_params *control,
     real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
     real BO, BO_s, BO_pi, BO_pi2;
     sparse_matrix *H, *H_sp;
-    list *far_nbrs, *bonds, *hbonds;
+    reax_list *far_nbrs, *bonds, *hbonds;
     single_body_parameters *sbp_i, *sbp_j;
     two_body_parameters *twbp;
     far_neighbor_data *nbr_pj;
@@ -922,7 +922,7 @@ void Init_Forces( reax_system *system, control_params *control,
 
 void Init_Forces_Tab( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, j, pj;
     int start_i, end_i;
@@ -935,7 +935,7 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
     real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
     real BO, BO_s, BO_pi, BO_pi2;
     sparse_matrix *H, *H_sp;
-    list *far_nbrs, *bonds, *hbonds;
+    reax_list *far_nbrs, *bonds, *hbonds;
     single_body_parameters *sbp_i, *sbp_j;
     two_body_parameters *twbp;
     far_neighbor_data *nbr_pj;
@@ -1201,7 +1201,7 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
 
 
 void Estimate_Storage_Sizes( reax_system *system, control_params *control,
-        list **lists, int *Htop, int *hb_top, int *bond_top, int *num_3body )
+        reax_list **lists, int *Htop, int *hb_top, int *bond_top, int *num_3body )
 {
     int i, j, pj;
     int start_i, end_i;
@@ -1210,7 +1210,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
     real r_ij;
     real C12, C34, C56;
     real BO, BO_s, BO_pi, BO_pi2;
-    list *far_nbrs;
+    reax_list *far_nbrs;
     single_body_parameters *sbp_i, *sbp_j;
     two_body_parameters *twbp;
     far_neighbor_data *nbr_pj;
@@ -1304,7 +1304,7 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
 
 void Compute_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list** lists, output_controls *out_control )
+        reax_list** lists, output_controls *out_control )
 {
     real t_start, t_elapsed;
 
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index c57aa15..4777325 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -27,8 +27,8 @@
 void Init_Bonded_Force_Functions( control_params* );
 
 void Compute_Forces( reax_system*, control_params*, simulation_data*,
-                     static_storage*, list**, output_controls* );
+                     static_storage*, reax_list**, output_controls* );
 
-void Estimate_Storage_Sizes( reax_system*, control_params*, list**,
+void Estimate_Storage_Sizes( reax_system*, control_params*, reax_list**,
                              int*, int*, int*, int* );
 #endif
diff --git a/sPuReMD/src/four_body_interactions.c b/sPuReMD/src/four_body_interactions.c
index e0010e8..26a15b7 100644
--- a/sPuReMD/src/four_body_interactions.c
+++ b/sPuReMD/src/four_body_interactions.c
@@ -154,10 +154,10 @@ real Calculate_Omega( rvec dvec_ij, real r_ij, rvec dvec_jk, real r_jk,
 
 void Four_Body_Interactions( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     real p_tor2, p_tor3, p_tor4, p_cot2;
-    list *bonds, *thb_intrs;
+    reax_list *bonds, *thb_intrs;
     real e_tor_total, e_con_total;
 
     p_tor2 = system->reaxprm.gp.l[23];
diff --git a/sPuReMD/src/four_body_interactions.h b/sPuReMD/src/four_body_interactions.h
index 37136cd..88261c5 100644
--- a/sPuReMD/src/four_body_interactions.h
+++ b/sPuReMD/src/four_body_interactions.h
@@ -25,6 +25,6 @@
 #include "mytypes.h"
 
 void Four_Body_Interactions( reax_system*, control_params*, simulation_data*,
-                             static_storage*, list**, output_controls* );
+                             static_storage*, reax_list**, output_controls* );
 
 #endif
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index dbcb9ab..48eb0db 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -444,7 +444,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
    cause trouble, if so we'll have to rethink this approach
    Also, we do not write connect lines yet.
 */
-char Write_PDB( reax_system* system, list* bonds, simulation_data *data,
+char Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
         control_params *control, static_storage *workspace, output_controls *out_control )
 {
     int i, buffer_req, buffer_len;
diff --git a/sPuReMD/src/geo_tools.h b/sPuReMD/src/geo_tools.h
index 4c44e30..0339a16 100644
--- a/sPuReMD/src/geo_tools.h
+++ b/sPuReMD/src/geo_tools.h
@@ -123,7 +123,7 @@ char Read_PDB( char*, reax_system*, control_params*,
 char Read_BGF( char*, reax_system*, control_params*,
         simulation_data*, static_storage* );
 
-char Write_PDB( reax_system*, list*, simulation_data*,
+char Write_PDB( reax_system*, reax_list*, simulation_data*,
         control_params*, static_storage*, output_controls* );
 
 #endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 850c044..d4175d4 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -507,7 +507,7 @@ void Init_Workspace( reax_system *system, control_params *control,
 
 void Init_Lists( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, num_nbrs, num_bonds, num_3body, Htop, max_nnz;
     int *hb_top, *bond_top;
@@ -531,6 +531,7 @@ void Init_Lists( reax_system *system, control_params *control,
     num_3body = 0;
     Estimate_Storage_Sizes( system, control, lists, &Htop,
             hb_top, bond_top, &num_3body );
+    num_3body = MAX( num_3body, MIN_BONDS );
 
     switch ( control->charge_method )
     {
@@ -570,7 +571,7 @@ void Init_Lists( reax_system *system, control_params *control,
 #endif
 
     workspace->num_H = 0;
-    if ( control->hb_cut > 0 )
+    if ( control->hb_cut > 0.0 )
     {
         /* init H indexes */
         for ( i = 0; i < system->N; ++i )
@@ -586,8 +587,15 @@ void Init_Lists( reax_system *system, control_params *control,
             }
         }
 
-        Allocate_HBond_List( system->N, workspace->num_H, workspace->hbond_index,
-                hb_top, (*lists) + HBONDS );
+        if ( workspace->num_H == 0 )
+        {
+            control->hb_cut = 0.0;
+        }
+        else
+        {
+            Allocate_HBond_List( system->N, workspace->num_H, workspace->hbond_index,
+                    hb_top, (*lists) + HBONDS );
+        }
 
 #if defined(DEBUG_FOCUS)
         num_hbonds = hb_top[system->N - 1];
@@ -854,8 +862,9 @@ void Init_Out_Controls( reax_system *system, control_params *control,
 
 
 void Initialize( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
-        output_controls *out_control, evolve_function *Evolve )
+        simulation_data *data, static_storage *workspace, reax_list **lists,
+        output_controls *out_control, evolve_function *Evolve,
+        const int output_enabled )
 {
 #if defined(DEBUG)
     real start, end;
@@ -879,7 +888,10 @@ void Initialize( reax_system *system, control_params *control,
 
     Init_Lists( system, control, data, workspace, lists, out_control );
 
-    Init_Out_Controls( system, control, workspace, out_control );
+    if ( output_enabled == TRUE )
+    {
+        Init_Out_Controls( system, control, workspace, out_control );
+    }
 
     /* These are done in forces.c, only forces.c can see all those functions */
     Init_Bonded_Force_Functions( control );
@@ -1118,10 +1130,13 @@ void Finalize_Workspace( reax_system *system, control_params *control,
 }
 
 
-void Finalize_Lists( list **lists )
+void Finalize_Lists( control_params *control, reax_list **lists )
 {
     Delete_List( TYP_FAR_NEIGHBOR, (*lists) + FAR_NBRS );
-    Delete_List( TYP_HBOND, (*lists) + HBONDS );
+    if ( control->hb_cut > 0.0 )
+    {
+        Delete_List( TYP_HBOND, (*lists) + HBONDS );
+    }
     Delete_List( TYP_BOND, (*lists) + BONDS );
     Delete_List( TYP_THREE_BODY, (*lists) + THREE_BODIES );
 
@@ -1341,17 +1356,20 @@ void Finalize_Out_Controls( reax_system *system, control_params *control,
 
 
 void Finalize( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
-        output_controls *out_control )
+        simulation_data *data, static_storage *workspace, reax_list **lists,
+        output_controls *out_control, const int output_enabled )
 {
     if ( control->tabulate )
     {
         Finalize_LR_Lookup_Table( system, control );
     }
 
-    Finalize_Out_Controls( system, control, workspace, out_control );
+    if ( output_enabled == TRUE )
+    {
+        Finalize_Out_Controls( system, control, workspace, out_control );
+    }
 
-    Finalize_Lists( lists );
+    Finalize_Lists( control, lists );
 
     Finalize_Workspace( system, control, workspace );
 
diff --git a/sPuReMD/src/init_md.h b/sPuReMD/src/init_md.h
index 3d31189..0643e0d 100644
--- a/sPuReMD/src/init_md.h
+++ b/sPuReMD/src/init_md.h
@@ -26,10 +26,11 @@
 
 
 void Initialize( reax_system*, control_params*, simulation_data*,
-        static_storage*, list**, output_controls*, evolve_function* );
+        static_storage*, reax_list**, output_controls*, evolve_function*,
+        const int );
 
 void Finalize( reax_system*, control_params*, simulation_data*,
-        static_storage*, list**, output_controls* );
+        static_storage*, reax_list**, output_controls*, const int );
 
 
 #endif
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 7f42a36..567e96a 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -37,7 +37,7 @@
 
 void Velocity_Verlet_NVE(reax_system* system, control_params* control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, steps, renbr;
     real inv_m, dt, dt_sqr;
@@ -68,7 +68,7 @@ void Velocity_Verlet_NVE(reax_system* system, control_params* control,
     fprintf( stderr, "verlet1 - ");
 #endif
 
-    Reallocate( system, workspace, lists, renbr );
+    Reallocate( system, control, workspace, lists, renbr );
     Reset( system, control, data, workspace, lists );
     if ( renbr )
     {
@@ -92,7 +92,7 @@ void Velocity_Verlet_NVE(reax_system* system, control_params* control,
 
 
 void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params* control,
-        simulation_data *data, static_storage *workspace, list **lists,
+        simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control )
 {
     int i, itr, steps, renbr;
@@ -130,7 +130,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein(reax_system* system, control_params*
     fprintf( stderr, "verlet1 - " );
 #endif
 
-    Reallocate( system, workspace, lists, renbr );
+    Reallocate( system, control, workspace, lists, renbr );
     Reset( system, control, data, workspace, lists );
     if ( renbr )
     {
@@ -209,7 +209,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
         control_params* control,
         simulation_data *data,
         static_storage *workspace,
-        list **lists,
+        reax_list **lists,
         output_controls *out_control )
 {
     int i, steps, renbr;
@@ -248,7 +248,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
     fprintf( stderr, "verlet1 - " );
 #endif
 
-    Reallocate( system, workspace, lists, renbr );
+    Reallocate( system, control, workspace, lists, renbr );
     Reset( system, control, data, workspace, lists );
     if ( renbr )
     {
@@ -324,7 +324,7 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
         control_params* control,
         simulation_data *data,
         static_storage *workspace,
-        list **lists,
+        reax_list **lists,
         output_controls *out_control )
 {
     int i, d, steps, renbr;
@@ -363,7 +363,7 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
     fprintf( stderr, "verlet1 - " );
 #endif
 
-    Reallocate( system, workspace, lists, renbr );
+    Reallocate( system, control, workspace, lists, renbr );
     Reset( system, control, data, workspace, lists );
     if ( renbr )
     {
@@ -450,7 +450,7 @@ void Velocity_Verlet_Nose_Hoover_NVT(reax_system* system,
                                      control_params* control,
                                      simulation_data *data,
                                      static_storage *workspace,
-                                     list **lists,
+                                     reax_list **lists,
                                      output_controls *out_control )
 {
     int i;
@@ -524,7 +524,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
                                     control_params* control,
                                     simulation_data *data,
                                     static_storage *workspace,
-                                    list **lists,
+                                    reax_list **lists,
                                     output_controls *out_control )
 {
     int i, itr;
@@ -701,7 +701,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
                                     control_params* control,
                                     simulation_data *data,
                                     static_storage *workspace,
-                                    list **lists,
+                                    reax_list **lists,
                                     output_controls *out_control
                                   )
 {
@@ -739,7 +739,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
     fprintf(stderr, "step%d: verlet1 done\n", data->step);
 #endif
 
-    Reallocate( system, workspace, lists, renbr );
+    Reallocate( system, control, workspace, lists, renbr );
     Reset( system, control, data, workspace, lists );
 
     if ( renbr )
diff --git a/sPuReMD/src/integrate.h b/sPuReMD/src/integrate.h
index 55d9878..1a821fe 100644
--- a/sPuReMD/src/integrate.h
+++ b/sPuReMD/src/integrate.h
@@ -25,25 +25,25 @@
 #include "mytypes.h"
 
 void Velocity_Verlet_NVE( reax_system*, control_params*, simulation_data*,
-                          static_storage*, list**, output_controls* );
+                          static_storage*, reax_list**, output_controls* );
 void Velocity_Verlet_Nose_Hoover_NVT( reax_system*, control_params*,
                                       simulation_data*, static_storage*,
-                                      list**, output_controls* );
+                                      reax_list**, output_controls* );
 void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system*, control_params*,
         simulation_data*, static_storage*,
-        list**, output_controls* );
+        reax_list**, output_controls* );
 void Velocity_Verlet_Flexible_NPT( reax_system*, control_params*,
                                    simulation_data*, static_storage*,
-                                   list**, output_controls* );
+                                   reax_list**, output_controls* );
 void Velocity_Verlet_Isotropic_NPT( reax_system*, control_params*,
                                     simulation_data*, static_storage*,
-                                    list**, output_controls* );
+                                    reax_list**, output_controls* );
 void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system*, control_params*,
         simulation_data*, static_storage*,
-        list**, output_controls* );
+        reax_list**, output_controls* );
 void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system*, control_params*,
         simulation_data*,
-        static_storage*, list**,
+        static_storage*, reax_list**,
         output_controls* );
 
 //upon Adri's request moved from parallel code to serial code
@@ -51,7 +51,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* ,
                                     control_params* ,
                                     simulation_data *,
                                     static_storage *,
-                                    list **,
+                                    reax_list **,
                                     output_controls *
                                   );
 #endif
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 576c01b..2094d4f 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -142,7 +142,7 @@ static void Sparse_MatVec( const sparse_matrix * const A,
  * A: stored in CSR
  * A_t: stored in CSR
  */
-void Transpose( const sparse_matrix * const A, sparse_matrix const *A_t )
+void Transpose( const sparse_matrix * const A, sparse_matrix * const A_t )
 {
     unsigned int i, j, pj, *A_t_top;
 
diff --git a/sPuReMD/src/lin_alg.h b/sPuReMD/src/lin_alg.h
index a314851..3946f60 100644
--- a/sPuReMD/src/lin_alg.h
+++ b/sPuReMD/src/lin_alg.h
@@ -31,7 +31,7 @@ typedef enum
 } TRIANGULARITY;
 
 
-void Transpose( const sparse_matrix * const, sparse_matrix const * );
+void Transpose( const sparse_matrix * const, sparse_matrix * const );
 
 void Transpose_I( sparse_matrix * const );
 
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index 351ee88..e6fbd2f 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -24,7 +24,7 @@
 #include "tool_box.h"
 
 
-void Make_List( int n, int total_intrs, int type, list* l )
+void Make_List( int n, int total_intrs, int type, reax_list* l )
 {
     l->n = n;
     l->total_intrs = total_intrs;
@@ -89,16 +89,10 @@ void Make_List( int n, int total_intrs, int type, list* l )
 }
 
 
-void Delete_List( int type, list* l )
+void Delete_List( int type, reax_list* l )
 {
-    if ( l->index != NULL )
-    {
-        sfree( l->index, "Delete_List::l->index" );
-    }
-    if ( l->end_index != NULL )
-    {
-        sfree( l->end_index, "Delete_List::l->end_index" );
-    }
+    sfree( l->index, "Delete_List::l->index" );
+    sfree( l->end_index, "Delete_List::l->end_index" );
 
     switch ( type )
     {
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index a270158..9cd53cb 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -25,36 +25,36 @@
 #include "mytypes.h"
 
 
-void Make_List( int, int, int, list* );
+void Make_List( int, int, int, reax_list* );
 
-void Delete_List( int, list* );
+void Delete_List( int, reax_list* );
 
 
-static inline int Num_Entries( int i, list* l )
+static inline int Num_Entries( int i, reax_list* l )
 {
     return l->end_index[i] - l->index[i];
 }
 
 
-static inline int Start_Index( int i, list *l )
+static inline int Start_Index( int i, reax_list *l )
 {
     return l->index[i];
 }
 
 
-static inline int End_Index( int i, list *l )
+static inline int End_Index( int i, reax_list *l )
 {
     return l->end_index[i];
 }
 
 
-static inline void Set_Start_Index( int i, int val, list *l )
+static inline void Set_Start_Index( int i, int val, reax_list *l )
 {
     l->index[i] = val;
 }
 
 
-static inline void Set_End_Index( int i, int val, list *l )
+static inline void Set_End_Index( int i, int val, reax_list *l )
 {
     l->end_index[i] = val;
 }
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index ccc54b8..f75ad16 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -80,7 +80,9 @@
   #define NAN_REAL(a) (0)
 #endif
 
-#define PI            3.14159265
+#if !defined(PI)
+  #define PI            3.14159265
+#endif
 #define C_ele          332.06371
 //#define K_B         503.398008   // kcal/mol/K
 #define K_B             0.831687   // amu A^2 / ps^2 / K
@@ -993,7 +995,7 @@ typedef struct
         /* hydrogen bond interaction list */
         hbond_data *hbond_list;
     } select;
-} list;
+} reax_list;
 
 
 typedef struct
@@ -1025,7 +1027,7 @@ typedef struct
     /* trajectory output function pointer definitions */
     int (* write_header)( reax_system*, control_params*, static_storage*, void* );
     int (* append_traj_frame)(reax_system*, control_params*,
-            simulation_data*, static_storage*, list **, void* );
+            simulation_data*, static_storage*, reax_list **, void* );
     int (* write)( FILE *, const char *, ... );
 
 #ifdef TEST_ENERGY
@@ -1116,11 +1118,11 @@ typedef struct
 
 /* Function pointer definitions */
 typedef void (*interaction_function)(reax_system*, control_params*,
-        simulation_data*, static_storage*, list**, output_controls*);
+        simulation_data*, static_storage*, reax_list**, output_controls*);
 
 typedef void (*evolve_function)(reax_system*, control_params*,
         simulation_data*, static_storage*,
-        list**, output_controls*);
+        reax_list**, output_controls*);
 
 
 /* Global variables */
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 5e23558..c0351c8 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -54,7 +54,7 @@ static inline real DistSqr_to_CP( rvec cp, rvec x )
 
 void Generate_Neighbor_Lists( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, j, k, l, m, itr;
     int x, y, z;
@@ -64,7 +64,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
     ivec *nbrs;
     rvec *nbrs_cp;
     grid *g;
-    list *far_nbrs;
+    reax_list *far_nbrs;
     far_neighbor_data *nbr_data;
     real t_start, t_elapsed;
 
@@ -174,7 +174,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
 
 
 int Estimate_NumNeighbors( reax_system *system, control_params *control,
-                           static_storage *workspace, list **lists )
+                           static_storage *workspace, reax_list **lists )
 {
     int  i, j, k, l, m, itr;
     int  x, y, z;
@@ -349,7 +349,7 @@ static inline int can_Bond( static_storage *workspace, int atom1, int atom2 )
 
 
 /* check if atom2 is on atom1's near neighbor list */
-static inline int is_Near_Neighbor( list *near_nbrs, int atom1, int atom2 )
+static inline int is_Near_Neighbor( reax_list *near_nbrs, int atom1, int atom2 )
 {
     int i;
 
@@ -367,7 +367,7 @@ static inline int is_Near_Neighbor( list *near_nbrs, int atom1, int atom2 )
 
 void Generate_Neighbor_Lists( reax_system *system, control_params *control,
                               simulation_data *data, static_storage *workspace,
-                              list **lists, output_controls *out_control )
+                              reax_list **lists, output_controls *out_control )
 {
     int  i, j, k;
     int  x, y, z;
@@ -377,9 +377,9 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
     int   c, count;
     int   grid_top;
     grid *g = &( system->g );
-    list *far_nbrs = (*lists) + FAR_NBRS;
+    reax_list *far_nbrs = (*lists) + FAR_NBRS;
     //int   hb_type1, hb_type2;
-    //list *hbonds = (*lists) + HBOND;
+    //reax_list *hbonds = (*lists) + HBOND;
     //int   top_hbond1, top_hbond2;
     get_far_neighbors_function Get_Far_Neighbors;
     far_neighbor_data new_nbrs[125];
@@ -617,7 +617,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
 
 void Generate_Neighbor_Lists( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int  i, j, k, l, m, itr;
     int  x, y, z;
@@ -627,7 +627,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
     ivec *nbrs;
     rvec *nbrs_cp;
     grid *g;
-    list *far_nbrs;
+    reax_list *far_nbrs;
     get_far_neighbors_function Get_Far_Neighbors;
     far_neighbor_data new_nbrs[125];
 
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index e849b1e..fde27eb 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -25,8 +25,8 @@
 #include "mytypes.h"
 
 void Generate_Neighbor_Lists( reax_system*, control_params*, simulation_data*,
-                              static_storage*, list**, output_controls* );
+                              static_storage*, reax_list**, output_controls* );
 
 int Estimate_NumNeighbors( reax_system*, control_params*,
-                           static_storage*, list** );
+                           static_storage*, reax_list** );
 #endif
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 684c2cb..1ef0a52 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -31,19 +31,19 @@
 #ifdef TEST_FORCES
 void Dummy_Printer( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
 }
 
 
 void Print_Bond_Orders( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int  i, pj, pk;
     bond_order_data *bo_ij;
-    list *bonds = (*lists) + BONDS;
-    list *dBOs  = (*lists) + DBO;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *dBOs  = (*lists) + DBO;
     dbond_data *dbo_k;
 
     /* bond orders */
@@ -110,7 +110,7 @@ void Print_Bond_Orders( reax_system *system, control_params *control,
 
 void Print_Bond_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
 
@@ -128,7 +128,7 @@ void Print_Bond_Forces( reax_system *system, control_params *control,
 
 void Print_LonePair_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
 
@@ -147,7 +147,7 @@ void Print_LonePair_Forces( reax_system *system, control_params *control,
 
 
 void Print_OverUnderCoor_Forces( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
+        simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control )
 {
     int i;
@@ -187,7 +187,7 @@ void Print_OverUnderCoor_Forces( reax_system *system, control_params *control,
 
 void Print_Three_Body_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int j;
 
@@ -245,7 +245,7 @@ void Print_Three_Body_Forces( reax_system *system, control_params *control,
 
 
 void Print_Hydrogen_Bond_Forces( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
+        simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control )
 {
     int j;
@@ -266,7 +266,7 @@ void Print_Hydrogen_Bond_Forces( reax_system *system, control_params *control,
 
 void Print_Four_Body_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int j;
 
@@ -304,7 +304,7 @@ void Print_Four_Body_Forces( reax_system *system, control_params *control,
 
 void Print_vdW_Coulomb_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int  i;
 
@@ -343,7 +343,7 @@ void Print_vdW_Coulomb_Forces( reax_system *system, control_params *control,
 
 void Compare_Total_Forces( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
 
@@ -399,12 +399,12 @@ void Init_Force_Test_Functions( )
 
 /* near nbrs contain both i-j, j-i nbrhood info */
 void Print_Near_Neighbors( reax_system *system, control_params *control,
-        static_storage *workspace, list **lists )
+        static_storage *workspace, reax_list **lists )
 {
     int i, j, id_i, id_j;
     char fname[MAX_STR];
     FILE *fout;
-    list *near_nbrs = &((*lists)[NEAR_NBRS]);
+    reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
 
     sprintf( fname, "%s.near_nbrs", control->sim_name );
     fout = fopen( fname, "w" );
@@ -433,12 +433,12 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
 
 /* near nbrs contain both i-j, j-i nbrhood info */
 void Print_Near_Neighbors2( reax_system *system, control_params *control,
-        static_storage *workspace, list **lists )
+        static_storage *workspace, reax_list **lists )
 {
     int i, j, id_i, id_j;
     char fname[MAX_STR];
     FILE *fout;
-    list *near_nbrs = &((*lists)[NEAR_NBRS]);
+    reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
 
     sprintf( fname, "%s.near_nbrs_lgj", control->sim_name );
     fout = fopen( fname, "w" );
@@ -468,12 +468,12 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
 
 /* far nbrs contain only i-j nbrhood info, no j-i. */
 void Print_Far_Neighbors( reax_system *system, control_params *control,
-        static_storage *workspace, list **lists )
+        static_storage *workspace, reax_list **lists )
 {
     int i, j, id_i, id_j;
     char fname[MAX_STR];
     FILE *fout;
-    list *far_nbrs = &((*lists)[FAR_NBRS]);
+    reax_list *far_nbrs = &((*lists)[FAR_NBRS]);
 
     sprintf( fname, "%s.far_nbrs", control->sim_name );
     fout = fopen( fname, "w" );
@@ -513,12 +513,12 @@ int fn_qsort_intcmp( const void *a, const void *b )
 
 
 void Print_Far_Neighbors2( reax_system *system, control_params *control,
-        static_storage *workspace, list **lists )
+        static_storage *workspace, reax_list **lists )
 {
     int i, j, id_i, id_j;
     char fname[MAX_STR];
     FILE *fout;
-    list *far_nbrs = &((*lists)[FAR_NBRS]);
+    reax_list *far_nbrs = &((*lists)[FAR_NBRS]);
 
     sprintf( fname, "%s.far_nbrs_lgj", control->sim_name );
     fout = fopen( fname, "w" );
@@ -548,7 +548,7 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
 
 void Print_Total_Force( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
 #if !defined(TEST_FORCES)
@@ -575,7 +575,7 @@ void Print_Total_Force( reax_system *system, control_params *control,
 
 void Output_Results( reax_system *system, control_params *control,
     simulation_data *data, static_storage *workspace,
-    list **lists, output_controls *out_control )
+    reax_list **lists, output_controls *out_control )
 {
     int i, type_i;
     real e_pol, q, f_update;
@@ -937,7 +937,7 @@ void Print_Sparse_Matrix_Binary( sparse_matrix *A, char *fname )
 }
 
 
-void Print_Bonds( reax_system *system, list *bonds, char *fname )
+void Print_Bonds( reax_system *system, reax_list *bonds, char *fname )
 {
     int i, pj;
     bond_data *pbond;
@@ -962,7 +962,7 @@ void Print_Bonds( reax_system *system, list *bonds, char *fname )
 }
 
 
-void Print_Bond_List2( reax_system *system, list *bonds, char *fname )
+void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
 {
     int i, j, id_i, id_j, nbr, pj;
     FILE *f = fopen( fname, "w" );
diff --git a/sPuReMD/src/print_utils.h b/sPuReMD/src/print_utils.h
index 8c15ebb..fe75d19 100644
--- a/sPuReMD/src/print_utils.h
+++ b/sPuReMD/src/print_utils.h
@@ -25,7 +25,7 @@
 #include "mytypes.h"
 
 typedef void (*print_interaction)(reax_system*, control_params*, simulation_data*,
-                                  static_storage*, list**, output_controls*);
+                                  static_storage*, reax_list**, output_controls*);
 print_interaction Print_Interactions[NO_OF_INTERACTIONS];
 
 char *Get_Element( reax_system*, int );
@@ -33,19 +33,19 @@ char *Get_Element( reax_system*, int );
 char *Get_Atom_Name( reax_system*, int );
 
 void Print_Near_Neighbors( reax_system*, control_params*, static_storage*,
-                           list** );
+                           reax_list** );
 
 void Print_Far_Neighbors( reax_system*, control_params*, static_storage*,
-                          list** );
+                          reax_list** );
 
 void Print_Total_Force( reax_system*, control_params*, simulation_data*,
-                        static_storage*, list**, output_controls* );
+                        static_storage*, reax_list**, output_controls* );
 
 void Output_Results( reax_system*, control_params*, simulation_data*,
-                     static_storage*, list**, output_controls* );
+                     static_storage*, reax_list**, output_controls* );
 
 void Print_Bond_Orders( reax_system*, control_params*, simulation_data*,
-                        static_storage*, list**, output_controls* );
+                        static_storage*, reax_list**, output_controls* );
 
 void Print_Linear_System( reax_system*, control_params*, static_storage*, int );
 
@@ -59,28 +59,28 @@ void Print_Sparse_Matrix2( sparse_matrix*, char*, char* );
 
 void Print_Sparse_Matrix_Binary( sparse_matrix*, char* );
 
-void Print_Bonds( reax_system*, list*, char* );
-void Print_Bond_List2( reax_system*, list*, char* );
+void Print_Bonds( reax_system*, reax_list*, char* );
+void Print_Bond_List2( reax_system*, reax_list*, char* );
 
 #ifdef TEST_FORCES
 void Dummy_Printer( reax_system*, control_params*, simulation_data*,
-                    static_storage*, list**, output_controls* );
+                    static_storage*, reax_list**, output_controls* );
 void Print_Bond_Forces( reax_system*, control_params*, simulation_data*,
-                        static_storage*, list**, output_controls* );
+                        static_storage*, reax_list**, output_controls* );
 void Print_LonePair_Forces( reax_system*, control_params*, simulation_data*,
-                            static_storage*, list**, output_controls* );
+                            static_storage*, reax_list**, output_controls* );
 void Print_OverUnderCoor_Forces(reax_system*, control_params*, simulation_data*,
-                                static_storage*, list**, output_controls*);
+                                static_storage*, reax_list**, output_controls*);
 void Print_Three_Body_Forces( reax_system*, control_params*, simulation_data*,
-                              static_storage*, list**, output_controls* );
+                              static_storage*, reax_list**, output_controls* );
 void Print_Hydrogen_Bond_Forces(reax_system*, control_params*, simulation_data*,
-                                static_storage*, list**, output_controls*);
+                                static_storage*, reax_list**, output_controls*);
 void Print_Four_Body_Forces( reax_system*, control_params*, simulation_data*,
-                             static_storage*, list**, output_controls* );
+                             static_storage*, reax_list**, output_controls* );
 void Print_vdW_Coulomb_Forces( reax_system*, control_params*, simulation_data*,
-                               static_storage*, list**, output_controls* );
+                               static_storage*, reax_list**, output_controls* );
 void Compare_Total_Forces( reax_system*, control_params*, simulation_data*,
-                           static_storage*, list**, output_controls* );
+                           static_storage*, reax_list**, output_controls* );
 void Init_Force_Test_Functions( );
 #endif
 
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index bf5330d..0aed2a9 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -115,11 +115,11 @@ void Reset_Workspace( reax_system *system, static_storage *workspace )
 
 
 void Reset_Neighbor_Lists( reax_system *system, control_params *control,
-                           static_storage *workspace, list **lists )
+                           static_storage *workspace, reax_list **lists )
 {
     int i, tmp;
-    list *bonds = (*lists) + BONDS;
-    list *hbonds = (*lists) + HBONDS;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *hbonds = (*lists) + HBONDS;
 
     for ( i = 0; i < system->N; ++i )
     {
@@ -145,7 +145,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
 
 
 void Reset( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists  )
+        simulation_data *data, static_storage *workspace, reax_list **lists  )
 {
     Reset_Atoms( system );
 
diff --git a/sPuReMD/src/reset_utils.h b/sPuReMD/src/reset_utils.h
index aadfe73..4355460 100644
--- a/sPuReMD/src/reset_utils.h
+++ b/sPuReMD/src/reset_utils.h
@@ -37,12 +37,12 @@ void Reset_Test_Forces( reax_system*, static_storage* );
 void Reset_Workspace( reax_system*, static_storage* );
 
 void Reset_Neighbor_Lists( reax_system*, control_params*,
-                           static_storage*, list** );
+                           static_storage*, reax_list** );
 
 void Reset( reax_system*, control_params*, simulation_data*,
-            static_storage*, list** );
+            static_storage*, reax_list** );
 
-//void Reset_Neighbor_Lists( reax_system*, static_storage*, list** );
+//void Reset_Neighbor_Lists( reax_system*, static_storage*, reax_list** );
 
 void Reset_Grid( grid* );
 
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index 9da3e21..b6edce1 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -27,7 +27,7 @@
 
 
 void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
+        simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control )
 {
     int i, j, pj, type_i, type_j;
@@ -45,7 +45,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
     two_body_parameters *twbp;
     bond_data *pbond;
     bond_order_data *bo_ij;
-    list *bonds = (*lists) + BONDS;
+    reax_list *bonds = (*lists) + BONDS;
 
     /* Initialize parameters */
     p_lp1 = system->reaxprm.gp.l[15];
diff --git a/sPuReMD/src/single_body_interactions.h b/sPuReMD/src/single_body_interactions.h
index 41cc649..8283c62 100644
--- a/sPuReMD/src/single_body_interactions.h
+++ b/sPuReMD/src/single_body_interactions.h
@@ -26,5 +26,5 @@
 
 void LonePair_OverUnder_Coordination_Energy( reax_system*, control_params*,
         simulation_data*, static_storage*,
-        list**, output_controls* );
+        reax_list**, output_controls* );
 #endif
diff --git a/sPuReMD/src/testmd.c b/sPuReMD/src/spuremd.c
similarity index 62%
rename from sPuReMD/src/testmd.c
rename to sPuReMD/src/spuremd.c
index 9ef2774..2e791a3 100644
--- a/sPuReMD/src/testmd.c
+++ b/sPuReMD/src/spuremd.c
@@ -36,9 +36,14 @@
 #include "vector.h"
 
 
+static static_storage workspace;
+static reax_list *lists;
+static output_controls out_control;
+
+
 static void Post_Evolve( reax_system * const system, control_params * const control,
         simulation_data * const data, static_storage * const workspace,
-        list ** const lists, output_controls * const out_control )
+        reax_list ** const lists, output_controls * const out_control )
 {
     int i;
     rvec diff, cross;
@@ -90,12 +95,14 @@ void static Read_System( char * const geo_file,
 
     if ( (ffield = fopen( ffield_file, "r" )) == NULL )
     {
-        fprintf( stderr, "Error opening the ffield file!\n" );
+        fprintf( stderr, "[ERROR] Error opening the ffield file!\n" );
+        fprintf( stderr, "    [INFO] (%s)\n", ffield_file );
         exit( FILE_NOT_FOUND );
     }
     if ( (ctrl = fopen( control_file, "r" )) == NULL )
     {
-        fprintf( stderr, "Error opening the ffield file!\n" );
+        fprintf( stderr, "[ERROR] Error opening the ffield file!\n" );
+        fprintf( stderr, "    [INFO] (%s)\n", control_file );
         exit( FILE_NOT_FOUND );
     }
 
@@ -130,7 +137,7 @@ void static Read_System( char * const geo_file,
     }
     else
     {
-        fprintf( stderr, "unknown geo file format. terminating!\n" );
+        fprintf( stderr, "[ERROR] unknown geo file format. terminating!\n" );
         exit( INVALID_GEO );
     }
 
@@ -144,84 +151,93 @@ void static Read_System( char * const geo_file,
 }
 
 
-static void usage( char* argv[] )
+int Setup( char ** args, reax_system * const system, control_params * const control,
+        simulation_data * const data )
 {
-    fprintf(stderr, "usage: ./%s geometry ffield control\n", argv[0]);
+    lists = (reax_list*) malloc( sizeof(reax_list) * LIST_N );
+
+    Read_System( args[0], args[1], args[2], system, control,
+            data, &workspace, &out_control );
+
+    return SUCCESS;
 }
 
 
-int main( int argc, char* argv[] )
+int Run( reax_system * const system, control_params * const control, simulation_data * const data,
+      const int output_enabled )
 {
-    reax_system system;
-    control_params control;
-    simulation_data data;
-    static_storage workspace;
-    list *lists;
-    output_controls out_control;
-    evolve_function Evolve;
     int steps;
+    evolve_function Evolve;
 
-    if ( argc != 4 )
-    {
-        usage(argv);
-        exit( INVALID_INPUT );
-    }
-
-    lists = (list*) malloc( sizeof(list) * LIST_N );
-
-    Read_System( argv[1], argv[2], argv[3], &system, &control,
-            &data, &workspace, &out_control );
-
-    Initialize( &system, &control, &data, &workspace, &lists,
-            &out_control, &Evolve );
+    Initialize( system, control, data, &workspace, &lists,
+            &out_control, &Evolve, output_enabled );
 
     /* compute f_0 */
     //if( control.restart == 0 ) {
-    Reset( &system, &control, &data, &workspace, &lists );
-    Generate_Neighbor_Lists( &system, &control, &data, &workspace, 
+    Reset( system, control, data, &workspace, &lists );
+    Generate_Neighbor_Lists( system, control, data, &workspace, 
             &lists, &out_control );
 
     //fprintf( stderr, "total: %.2f secs\n", data.timing.nbrs);
-    Compute_Forces(&system, &control, &data, &workspace, &lists, &out_control);
-    Compute_Kinetic_Energy( &system, &data );
-    Output_Results(&system, &control, &data, &workspace, &lists, &out_control);
-    ++data.step;
+    Compute_Forces( system, control, data, &workspace, &lists, &out_control );
+    Compute_Kinetic_Energy( system, data );
+    if ( output_enabled == TRUE )
+    {
+        Output_Results( system, control, data, &workspace, &lists, &out_control );
+    }
+    ++data->step;
     //}
     //
     
-    for ( ; data.step <= control.nsteps; data.step++ )
+    for ( ; data->step <= control->nsteps; data->step++ )
     {
-        if ( control.T_mode )
+        if ( control->T_mode )
+        {
+            Temperature_Control( control, data, &out_control );
+        }
+
+        Evolve( system, control, data, &workspace, &lists, &out_control );
+        Post_Evolve( system, control, data, &workspace, &lists, &out_control );
+
+        if ( output_enabled == TRUE )
         {
-            Temperature_Control( &control, &data, &out_control );
+            Output_Results( system, control, data, &workspace, &lists, &out_control );
+            Analysis( system, control, data, &workspace, &lists, &out_control );
         }
-        Evolve( &system, &control, &data, &workspace, &lists, &out_control );
-        Post_Evolve( &system, &control, &data, &workspace, &lists, &out_control );
-        Output_Results( &system, &control, &data, &workspace, &lists, &out_control );
-        Analysis( &system, &control, &data, &workspace, &lists, &out_control );
 
-        steps = data.step - data.prev_steps;
+        steps = data->step - data->prev_steps;
         if ( steps && out_control.restart_freq &&
-                steps % out_control.restart_freq == 0 )
+                steps % out_control.restart_freq == 0 &&
+                output_enabled == TRUE )
         {
-            Write_Restart( &system, &control, &data, &workspace, &out_control );
+            Write_Restart( system, control, data, &workspace, &out_control );
         }
     }
 
-    if ( out_control.write_steps > 0 )
+    if ( out_control.write_steps > 0 && output_enabled == TRUE )
     {
         fclose( out_control.trj );
-        Write_PDB( &system, &(lists[BONDS]), &data, &control, &workspace, &out_control );
+        Write_PDB( system, &(lists[BONDS]), data, control, &workspace, &out_control );
     }
 
-    data.timing.end = Get_Time( );
-    data.timing.elapsed = Get_Timing_Info( data.timing.start );
-    fprintf( out_control.log, "total: %.2f secs\n", data.timing.elapsed );
+    data->timing.end = Get_Time( );
+    data->timing.elapsed = Get_Timing_Info( data->timing.start );
+    if ( output_enabled == TRUE )
+    {
+        fprintf( out_control.log, "total: %.2f secs\n", data->timing.elapsed );
+    }
+
+    return SUCCESS;
+}
+
 
-    Finalize( &system, &control, &data, &workspace, &lists,
-            &out_control );
+int Cleanup( reax_system * const system, control_params * const control,
+        simulation_data * const data, const int output_enabled )
+{
+    Finalize( system, control, data, &workspace, &lists, &out_control,
+           output_enabled );
 
     sfree( lists, "main::lists" );
 
-    return 0;
+    return SUCCESS;
 }
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
new file mode 100644
index 0000000..7cad1a6
--- /dev/null
+++ b/sPuReMD/src/spuremd.h
@@ -0,0 +1,38 @@
+/*----------------------------------------------------------------------
+  SerialReax - Reax Force Field Simulator
+
+  Copyright (2010) Purdue University
+  Hasan Metin Aktulga, haktulga@cs.purdue.edu
+  Joseph Fogarty, jcfogart@mail.usf.edu
+  Sagar Pandit, pandit@usf.edu
+  Ananth Y Grama, ayg@cs.purdue.edu
+
+  This program is free software; you can redistribute it and/or
+  modify it under the terms of the GNU General Public License as
+  published by the Free Software Foundation; either version 2 of
+  the License, or (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+  See the GNU General Public License for more details:
+  <http://www.gnu.org/licenses/>.
+  ----------------------------------------------------------------------*/
+
+#ifndef __SPUREMD_H_
+#define __SPUREMD_H_
+
+#include "mytypes.h"
+
+
+int Setup( char **, reax_system * const, control_params * const,
+        simulation_data * const );
+
+int Run( reax_system * const, control_params * const,
+        simulation_data * const, const int );
+
+int Cleanup( reax_system * const, control_params * const,
+        simulation_data * const, const int );
+
+
+#endif
diff --git a/sPuReMD/src/three_body_interactions.c b/sPuReMD/src/three_body_interactions.c
index 533c652..159e0cc 100644
--- a/sPuReMD/src/three_body_interactions.c
+++ b/sPuReMD/src/three_body_interactions.c
@@ -80,10 +80,10 @@ void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk,
    played by j which sits in the middle of the other two. */
 void Three_Body_Interactions( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     real *total_bo;
-    list *bonds, *thb_intrs;
+    reax_list *bonds, *thb_intrs;
     bond_data *bond_list;
     three_body_interaction_data *thb_list;
     real p_pen2, p_pen3, p_pen4;
@@ -667,7 +667,7 @@ void Three_Body_Interactions( reax_system *system, control_params *control,
 
 void Hydrogen_Bonds( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     real e_hb_total;
 
@@ -694,7 +694,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
         bond_order_data *bo_ij;
         bond_data *pbond_ij;
         far_neighbor_data *nbr_jk;
-        list *bonds, *hbonds;
+        reax_list *bonds, *hbonds;
         bond_data *bond_list;
         hbond_data *hbond_list;
         rvec *f_i, *f_j, *f_k;
diff --git a/sPuReMD/src/three_body_interactions.h b/sPuReMD/src/three_body_interactions.h
index 55a1188..c654751 100644
--- a/sPuReMD/src/three_body_interactions.h
+++ b/sPuReMD/src/three_body_interactions.h
@@ -25,10 +25,10 @@
 #include "mytypes.h"
 
 void Three_Body_Interactions( reax_system*, control_params*, simulation_data*,
-                              static_storage*, list**, output_controls* );
+                              static_storage*, reax_list**, output_controls* );
 
 void Hydrogen_Bonds( reax_system*, control_params*, simulation_data*,
-                     static_storage*, list**, output_controls* );
+                     static_storage*, reax_list**, output_controls* );
 
 void Calculate_Theta( rvec, real, rvec, real, real*, real* );
 
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 61ab5e7..c37ccc3 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -402,12 +402,12 @@ int Tokenize( char* s, char*** tok )
 {
     char test[MAX_LINE];
     char *sep = "\t \n!=";
-    char *word;
+    char *word, *saveptr;
     int count = 0;
 
     strncpy( test, s, MAX_LINE );
 
-    for ( word = strtok(test, sep); word; word = strtok(NULL, sep) )
+    for ( word = strtok_r(test, sep, &saveptr); word; word = strtok_r(NULL, sep, &saveptr) )
     {
         strncpy( (*tok)[count], word, MAX_LINE );
         count++;
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index 017d657..bfdd2fc 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -164,7 +164,7 @@ int Write_Custom_Header( reax_system *system, control_params *control,
 
 int Append_Custom_Frame( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i, j, pi, pk, pk_j;
     int write_atoms, write_bonds, write_angles;
@@ -172,8 +172,8 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
     int frame_globals_len, num_bonds, num_thb_intrs;
     real P;
     char buffer[2048];
-    list *bonds = (*lists) + BONDS;
-    list *thb_intrs =  (*lists) + THREE_BODIES;
+    reax_list *bonds = (*lists) + BONDS;
+    reax_list *thb_intrs =  (*lists) + THREE_BODIES;
     bond_data *bo_ij;
 
     /* IMPORTANT: This whole part will go to init_trj after finalized! */
@@ -534,7 +534,7 @@ int Write_xyz_Header( reax_system *system, control_params *control,
 
 int Append_xyz_Frame( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
 
diff --git a/sPuReMD/src/traj.h b/sPuReMD/src/traj.h
index 0e75905..e45921d 100644
--- a/sPuReMD/src/traj.h
+++ b/sPuReMD/src/traj.h
@@ -164,11 +164,11 @@ char Write_Traj_Header( FILE*, int, char**, char**, control_params* );
           control_params* control,
           simulation_data* data,
           static_storage* workspace,
-          list** lists,
+          reax_list** lists,
           char** various flags);
 */
 int Push_Traj_Frame( /*gzfile*/ FILE*, reax_system*, control_params*,
-                                simulation_data*, static_storage*, list**, char** );
+                                simulation_data*, static_storage*, reax_list**, char** );
 
 /*
   Append_Traj_Frame( gzfile file,
@@ -176,13 +176,13 @@ int Push_Traj_Frame( /*gzfile*/ FILE*, reax_system*, control_params*,
                         control_params* control,
                 simulation_data* data,
                 static_storage* workspace,
-                list** lists,
+                reax_list** lists,
                 char** various flags);
 */
 int Append_Custom_Frame( reax_system*, control_params*, simulation_data*,
-                         static_storage*, list**, output_controls* );
+                         static_storage*, reax_list**, output_controls* );
 int Append_xyz_Frame   ( reax_system*, control_params*, simulation_data*,
-                         static_storage*, list**, output_controls* );
+                         static_storage*, reax_list**, output_controls* );
 
 
 void Read_Traj( output_controls*, char * );
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index 49fa32b..e533d18 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -29,11 +29,11 @@
 
 void Bond_Energy( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
     real gp3, gp4, gp7, gp10, gp37, ebond_total;
-    list *bonds;
+    reax_list *bonds;
 
     bonds = (*lists) + BONDS;
     gp3 = system->reaxprm.gp.l[3];
@@ -168,11 +168,11 @@ void Bond_Energy( reax_system *system, control_params *control,
 
 void vdW_Coulomb_Energy( reax_system *system, control_params *control,
         simulation_data *data, static_storage *workspace,
-        list **lists, output_controls *out_control )
+        reax_list **lists, output_controls *out_control )
 {
     int i;
     real p_vdW1, p_vdW1i;
-    list *far_nbrs;
+    reax_list *far_nbrs;
     real e_vdW_total, e_ele_total;
 
     p_vdW1 = system->reaxprm.gp.l[28];
@@ -520,11 +520,11 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
 
 
 void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
-        simulation_data *data, static_storage *workspace, list **lists,
+        simulation_data *data, static_storage *workspace, reax_list **lists,
         output_controls *out_control )
 {
     int steps, update_freq, update_energies;
-    list *far_nbrs;
+    reax_list *far_nbrs;
     real e_vdW_total, e_ele_total;
 
     far_nbrs = (*lists) + FAR_NBRS;
diff --git a/sPuReMD/src/two_body_interactions.h b/sPuReMD/src/two_body_interactions.h
index bee9779..b176b5e 100644
--- a/sPuReMD/src/two_body_interactions.h
+++ b/sPuReMD/src/two_body_interactions.h
@@ -25,10 +25,10 @@
 #include <mytypes.h>
 
 void Bond_Energy( reax_system*, control_params*, simulation_data*,
-                  static_storage*, list**, output_controls* );
+                  static_storage*, reax_list**, output_controls* );
 void vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
-                         static_storage*, list**, output_controls* );
+                         static_storage*, reax_list**, output_controls* );
 void LR_vdW_Coulomb( reax_system*, control_params*, int, int, real, LR_data* );
 void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*, simulation_data*,
-                                   static_storage*, list**, output_controls* );
+                                   static_storage*, reax_list**, output_controls* );
 #endif
-- 
GitLab