diff --git a/tools/driver.py b/tools/driver.py
index 0144e62cb45ad5b107deb565d6c5b1877a32725f..12a4a0882981ef73d8b252dbb360505e830d35f4 100644
--- a/tools/driver.py
+++ b/tools/driver.py
@@ -3,6 +3,123 @@
 from ctypes import *
 
 
+class BondOrderData(Structure):
+    _fields_ = [
+            ("BO", c_double),
+            ("BO_s", c_double),
+            ("BO_pi", c_double),
+            ("BO_pi2", c_double),
+            ("Cdbo", c_double),
+            ("Cdbopi", c_double),
+            ("Cdbopi2", c_double),
+            ("C1dbo", c_double),
+            ("C2dbo", c_double),
+            ("C3dbo", c_double),
+            ("C1dbopi", c_double),
+            ("C2dbopi", c_double),
+            ("C3dbopi", c_double),
+            ("C4dbopi", c_double),
+            ("C1dbopi2", c_double),
+            ("C2dbopi2", c_double),
+            ("C3dbopi2", c_double),
+            ("C4dbopi2", c_double),
+            ("dBOp", c_double * 3),
+            ("dln_BOp_s", c_double * 3),
+            ("dln_BOp_pi", c_double * 3),
+            ("dln_BOp_pi2", c_double * 3),
+            ]
+
+
+class ThreeBodyData(Structure):
+    _fields_ = [
+            ("thb", c_int),
+            ("pthb", c_int),
+            ("theta", c_double),
+            ("cos_theta", c_double),
+            ("dcos_di", c_double * 3),
+            ("dcos_dj", c_double * 3),
+            ("dcos_dk", c_double * 3),
+            ]
+
+
+class BondData(Structure):
+    _fields_ = [
+            ("nbr", c_int),
+            ("sym_index", c_int),
+            ("dbond_index", c_int),
+            ("rel_box", c_int * 3),
+            ("d", c_double),
+            ("dvec", c_double * 3),
+            ("bo_data", BondOrderData),
+            ]
+
+
+class DBondData(Structure):
+    _fields_ = [
+            ("wrt", c_int),
+            ("dBO", c_double * 3),
+            ("dBOpi", c_double * 3),
+            ("dBOpi2", c_double * 3),
+            ]
+
+
+class DDeltaData(Structure):
+    _fields_ = [
+            ("wrt", c_int),
+            ("dVal", c_double * 3),
+            ]
+
+
+class FarNbrData(Structure):
+    _fields_ = [
+            ("nbr", c_int),
+            ("rel_box", c_int * 3),
+            ("d", c_double),
+            ("dvec", c_double * 3),
+            ]
+
+
+class NearNbrData(Structure):
+    _fields_ = [
+            ("nbr", c_int),
+            ("rel_box", c_int * 3),
+            ("d", c_double),
+            ("dvec", c_double * 3),
+            ]
+
+
+class HBondData(Structure):
+    _fields_ = [
+            ("nbr", c_int),
+            ("scl", c_int),
+            ("ptr", POINTER(FarNbrData)),
+            ]
+
+
+class ReaxListSelector(Union):
+    _fields_ = [
+            ("v", c_void_p),
+            ("three_body_interaction_data", POINTER(ThreeBodyData)),
+            ("bond_data", POINTER(BondData)),
+            ("dbond_data", POINTER(DBondData)),
+            ("dDelta_data", POINTER(DDeltaData)),
+            ("far_neighbor_data", POINTER(FarNbrData)),
+            ("near_neighbor_data", POINTER(NearNbrData)),
+            ("hbond_data", POINTER(HBondData)),
+            ]
+
+
+class ReaxList(Structure):
+    _fields_ = [
+            ("n", c_int),
+            ("total_intrs", c_int),
+            ("index", POINTER(c_int)),
+            ("end_index", POINTER(c_int)),
+            ("max_intrs", POINTER(c_int)),
+            ("select", ReaxListSelector),
+            ]
+
+
 class Thermostat(Structure):
     _fields_ = [
             ("T", c_double),
@@ -133,10 +250,35 @@ if __name__ == '__main__':
     cleanup.argtypes = [c_void_p]
     cleanup.restype = c_int
 
+    get_atoms = lib.get_atoms
+    get_atoms.argtypes = [c_void_p]
+    get_atoms.restype = POINTER(ReaxAtom)
+
+    CALLBACKFUNC = CFUNCTYPE(None, POINTER(ReaxAtom),
+            POINTER(SimulationData), POINTER(ReaxList))
+
+    def get_simulation_step_results(atoms, data, lists):
+        print("{0:24.15f} {1:24.15f} {2:24.15f}".format(
+            data[0].E_Tot, data[0].E_Kin, data[0].E_Pot))
+
+    setup_callback = lib.setup_callback
+    setup_callback.restype = c_int
+
     handle = setup(b"data/benchmarks/water/water_6540.pdb",
             b"data/benchmarks/water/ffield.water",
             b"environ/param.gpu.water")
 
-    simulate(handle)
+    ret = setup_callback(handle, CALLBACKFUNC(get_simulation_step_results))
+
+    print("{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
+    ret = simulate(handle)
+
+    atoms = get_atoms(handle)
+
+    print()
+    print("{0:9}|{1:24}|{2:24}|{3:24}|{4:24}".format("Atom Num", "x-Position", "y-Position", "z-Position", "Charge"))
+    for i in range(10):
+        print("{0:9d} {1:24.15f} {2:24.15f} {3:24.15f} {4:24.15f}".format(
+            i + 1, atoms[i].x[0], atoms[i].x[1], atoms[i].x[2], atoms[i].q))
 
     cleanup(handle)
diff --git a/tools/run_sim.py b/tools/run_sim.py
index e8c56f0d61a59f327651cbb96a968e56746eda3f..3f3af62e8840486ef17a47a5000f9819a0bcb9b1 100644
--- a/tools/run_sim.py
+++ b/tools/run_sim.py
@@ -79,9 +79,11 @@ class TestCase():
         fp_temp.write(lines)
         fp_temp.close()
 
-    def run(self, bin_file='sPuReMD/bin/spuremd', process_results=False):
-        base_dir = getcwd()
-        bin_path = path.join(base_dir, bin_file)
+    def run(self, binary, process_results=False):
+        args = binary.split()
+        args.append(self.__geo_file)
+        args.append( self.__ffield_file)
+        args.append("")
         env = dict(environ)
 
         write_header = True
@@ -111,14 +113,13 @@ class TestCase():
                 + '_paji' + param_dict['cm_solver_pre_app_jacobi_iters'] \
 		+ '_t' + param_dict['threads']
 
-
             if not process_results:
                 self._setup(param_dict, temp_file)
     
                 env['OMP_NUM_THREADS'] = param_dict['threads']
                 start = time()
-                proc_handle = Popen([bin_path, self.__geo_file, self.__ffield_file, temp_file], 
-                    stdout=PIPE, stderr=PIPE, env=env, universal_newlines=True)
+                args[-1] = temp_file
+                proc_handle = Popen(args, stdout=PIPE, stderr=PIPE, env=env, universal_newlines=True)
                 stdout, stderr = proc_handle.communicate()
                 stop = time()
                 if proc_handle.returncode < 0:
@@ -202,6 +203,8 @@ if __name__ == '__main__':
             ]
 
     parser = argparse.ArgumentParser(description='Run molecular dynamics simulations on specified data sets.')
+    parser.add_argument('-b', '--binary', metavar='binary', default=None, nargs=1,
+            help='Binary file to run.')
     parser.add_argument('-f', '--out_file', metavar='out_file', default=None, nargs=1,
             help='Output file to write results.')
     parser.add_argument('-r', '--process_results', default=False, action='store_true',
@@ -250,6 +253,11 @@ if __name__ == '__main__':
     else:
         result_file = 'result.txt'
 
+    if args.binary:
+        binary = args.binary[0]
+    else:
+        binary  = path.join(base_dir, 'sPuReMD/bin/spuremd')
+
     # overwrite default params, if supplied via command line args
     if args.params:
         for param in args.params:
@@ -342,4 +350,4 @@ if __name__ == '__main__':
                 geo_format=['1'], result_file=result_file))
 
     for test in test_cases:
-        test.run(bin_file='sPuReMD/bin/spuremd', process_results=args.process_results)
+        test.run(binary, process_results=args.process_results)