diff --git a/environ/param.gpu.water b/environ/param.gpu.water
index 90288d3da7b9a108722c74a061ce47dd1d76c5e8..edd2784aa4f0a7a25df38e415d070f9cc392083c 100644
--- a/environ/param.gpu.water
+++ b/environ/param.gpu.water
@@ -35,7 +35,7 @@ compress                0.008134                ! in ps^2 * A / amu ( 4.5X10^(-5
 press_mode              0                       ! 0: internal + external pressure, 1: ext only, 2: int only
 
 geo_format              1                       ! 0: xyz, 1: pdb, 2: bgf
-write_freq              1000                    ! write trajectory after so many steps
+write_freq              0                       ! write trajectory after so many steps
 traj_compress           0                       ! 0: no compression  1: uses zlib to compress trajectory output
 traj_format             0                       ! 0: our own format (below options apply to this only), 1: xyz, 2: bgf, 3: pdb
 traj_title              WATER_NVT               ! (no white spaces)
diff --git a/tools/run_sim.py b/tools/run_sim.py
index de1c0372609af6c4eccce3c8a2b8dd95a531a03e..39342fdb03a3f71591453199b561fbc9ca9ba7c8 100644
--- a/tools/run_sim.py
+++ b/tools/run_sim.py
@@ -1,91 +1,197 @@
 #!/usr/bin/env python
 
 from fileinput import input
+from itertools import product
 from re import sub
 from subprocess import Popen, PIPE
-from os import getcwd, environ, path, rename
+from os import getcwd, environ, path, remove, rename, rmdir
 from sys import exit
+from tempfile import mkdtemp
 from time import time
 
-base_dir = getcwd()
-control_dir = path.join(base_dir, 'environ')
-data_dir = path.join(base_dir, 'data/benchmarks')
-puremd_args = [ \
-#  (path.join(data_dir, 'water/water_6540.pdb'), path.join(data_dir, 'water/ffield.water'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'water/water_78480.pdb'), path.join(data_dir, 'water/ffield.water'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'water/water_327000.geo'), path.join(data_dir, 'water/ffield.water'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'bilayer/bilayer_56800.pdb'), path.join(data_dir, 'bilayer/ffield-bio'), path.join(control_dir, 'param.gpu.water')), \
-  (path.join(data_dir, 'dna/dna_19733.pdb'), path.join(data_dir, 'dna/ffield-dna'), path.join(control_dir, 'param.gpu.water')), \
-  (path.join(data_dir, 'silica/silica_6000.pdb'), path.join(data_dir, 'silica/ffield-bio'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'silica/silica_72000.geo'), path.join(data_dir, 'silica/ffield-bio'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'silica/silica_300000.geo'), path.join(data_dir, 'silica/ffield-bio'), path.join(control_dir, 'param.gpu.water')), \
-#  (path.join(data_dir, 'petn/petn_48256.pdb'), path.join(data_dir, 'petn/ffield.petn'), path.join(control_dir, 'param.gpu.water')), \
-]
-
-header_fmt_str = '{:20}|{:5}|{:5}|{:5}|{:10}|{:10}|{:10}|{:10}|{:10}|{:10}'
-body_fmt_str = '{:20} {:5} {:5} {:5} {:10.6f} {:10.6f} {:10.6f} {:10.6f} {:10} {:10.6f}'
-
-steps = ['20'] #['100']
-qeq_solver_tol = ['1e-6'] #['1e-6', '1e-8', '1e-10', '1e-14']
-pre_comp_type = ['2'] #['0', '1', '2']
-threads = ['1', '2', '4', '12', '24']
-patterns = [ \
-  lambda l, x: sub(r'(?P<key>simulation_name\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
-  lambda l, x: sub(r'(?P<key>nsteps\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
-  lambda l, x: sub(r'(?P<key>qeq_solver_q_err\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
-  lambda l, x: sub(r'(?P<key>pre_comp_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
-  lambda l, x: sub(r'(?P<key>geo_format\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
-]
-d = dict(environ)
-
-print header_fmt_str.format('Data Set', 'Steps', 'Q Tol', 'Pre T', 'Pre Comp',
-    'Pre App', 'Iters', 'SpMV', 'Threads', 'Time (s)')
-
-for i in xrange(len(puremd_args)):
-  for s in xrange(len(steps)):
-    for j in xrange(len(qeq_solver_tol)):
-      for t in xrange(len(pre_comp_type)):
-        for k in xrange(len(threads)):
-          for line in input(puremd_args[i][2], inplace=1):
-            line = line.rstrip() 
-            line = patterns[0](line, path.basename(puremd_args[i][0]).split('.')[0]
-                + '_' + 'tol' + qeq_solver_tol[j] + '_precomp' + pre_comp_type[t] + '_thread' + threads[k])
-            line = patterns[1](line, steps[s])
-            line = patterns[2](line, qeq_solver_tol[j])
-            line = patterns[3](line, pre_comp_type[t])
-            print line
-  
-          d['OMP_NUM_THREADS'] = threads[k]
-          start = time()
-          p = Popen([path.join(base_dir, 'sPuReMD/bin/spuremd')] + list(puremd_args[i]), 
-              stdout=PIPE, stderr=PIPE, env=d)
-          stdout, stderr = p.communicate()
-          stop = time()
-  
-          iters = 0.
-          pre_comp = 0.
-          pre_app = 0.
-          spmv = 0.
-          cnt = 0
-          with open(path.basename(puremd_args[i][0]).split('.')[0]
-              + '_' + 'tol' + qeq_solver_tol[j] + '_precomp' + pre_comp_type[t] + '_thread'
-  	      + threads[k] + '.log', 'r') as fp:
-            for line in fp:
-              line = line.split()
-              try:
-                iters = iters + float(line[7])
-                pre_comp = pre_comp + float(line[8])
-                pre_app = pre_app + float(line[9])
-                spmv = spmv + float(line[10])
-                cnt = cnt + 1
-              except Exception:
+
+class TestCase():
+    def __init__(self, geo_file, ffield_file, control_file, params={}, result_header_fmt='',
+            result_header='', result_body_fmt='', result_file='results.txt'):
+        self.__geo_file = geo_file
+        self.__ffield_file = ffield_file
+        self.__control_file = control_file
+        self.__param_names = sorted(params.keys())
+        self.__params = params
+        self.__result_header_fmt = result_header_fmt
+        self.__result_header = result_header
+        self.__result_body_fmt = result_body_fmt
+        self.__result_file = result_file
+        self.__control_res = { \
+                'name': lambda l, x: sub(
+                    r'(?P<key>simulation_name\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'ensemble_type': lambda l, x: sub(
+                    r'(?P<key>ensemble_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'nsteps': lambda l, x: sub(
+                    r'(?P<key>nsteps\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'qeq_solver_type': lambda l, x: sub(
+                    r'(?P<key>qeq_solver_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'qeq_solver_q_err': lambda l, x: sub(
+                    r'(?P<key>qeq_solver_q_err\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'pre_comp_type': lambda l, x: sub(
+                    r'(?P<key>pre_comp_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'pre_comp_refactor': lambda l, x: sub(
+                    r'(?P<key>pre_comp_refactor\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'pre_comp_sweeps': lambda l, x: sub(
+                    r'(?P<key>pre_comp_sweeps\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'pre_app_type': lambda l, x: sub(
+                    r'(?P<key>pre_app_type\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'pre_app_jacobi_iters': lambda l, x: sub(
+                    r'(?P<key>pre_app_jacobi_iters\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+                'geo_format': lambda l, x: sub(
+                    r'(?P<key>geo_format\s+)\S+(?P<comment>.*)', r'\g<key>%s\g<comment>' % x, l), \
+        }
+
+    def _setup(self, param, temp_file):
+        fp = open(self.__control_file, 'r')
+	lines = fp.read()
+	fp.close()
+        for k in self.__control_res.keys():
+            try:
+                lines = self.__control_res[k](lines, param[k])
+            except KeyError:
                 pass
-          cnt = cnt - 1
-          iters = iters / cnt
-          pre_comp = pre_comp / cnt
-          pre_app = pre_app / cnt
-          spmv = spmv / cnt
-  
-          print body_fmt_str.format(path.basename(puremd_args[i][0]).split('.')[0], 
-              steps[s], qeq_solver_tol[j], pre_comp_type[t], pre_comp, pre_app, iters, spmv,
-  	      threads[k], stop - start)
+        fp_temp = open(temp_file, 'w', 0)
+	fp_temp.write(lines)
+	fp_temp.close()
+
+    def run(self, bin_file='sPuReMD/bin/spuremd'):
+        base_dir = getcwd()
+	bin_path = path.join(base_dir, bin_file)
+        env = dict(environ)
+
+        write_header = True
+        if path.exists(self.__result_file):
+            write_header = False
+        fout = open(self.__result_file, 'a', 0)
+        if write_header:
+            fout.write(self.__result_header_fmt.format(*self.__result_header))
+
+        temp_dir = mkdtemp()
+	temp_file = path.join(temp_dir, path.basename(self.__control_file))
+
+        for p in product(*[self.__params[k] for k in self.__param_names]):
+            param_dict = dict((k, v) for (k, v) in zip(self.__param_names, p))
+            param_dict['name'] = path.basename(self.__geo_file).split('.')[0] \
+                + '_step' + param_dict['nsteps'] + '_tol' + param_dict['qeq_solver_q_err'] \
+                + '_precomp' + param_dict['pre_comp_type'] + '_thread' + param_dict['threads']
+
+            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)
+            #TODO: handle outputs?
+            stdout, stderr = proc_handle.communicate()
+            stop = time()
+
+            self._process_result(fout, stop - start, param_dict)
+
+        fout.close()
+        remove(temp_file)
+	rmdir(temp_dir)
+
+    def _process_result(self, fout, time, param):
+        qeq = 0.
+        iters = 0.
+        pre_comp = 0.
+        pre_app = 0.
+        spmv = 0.
+        cnt = 0
+        with open(param['name'] + '.log', 'r') as fp:
+            for line in fp:
+                line = line.split()
+                try:
+                    qeq = qeq + float(line[6])
+                    iters = iters + float(line[7])
+                    pre_comp = pre_comp + float(line[8])
+                    pre_app = pre_app + float(line[9])
+                    spmv = spmv + float(line[10])
+                    cnt = cnt + 1
+                except Exception:
+                    pass
+            cnt = cnt - 1
+            qeq = qeq / cnt
+            iters = iters / cnt
+            pre_comp = pre_comp / cnt
+            pre_app = pre_app / cnt
+            spmv = spmv / cnt
+
+        fout.write(self.__result_body_fmt.format(path.basename(self.__geo_file).split('.')[0], 
+            param['nsteps'], param['qeq_solver_q_err'], param['pre_comp_type'], pre_comp, pre_app, iters, spmv,
+            qeq, param['threads'], time))
+
+
+if __name__ == '__main__':
+    base_dir = getcwd()
+    control_dir = path.join(base_dir, 'environ')
+    data_dir = path.join(base_dir, 'data/benchmarks')
+
+    header_fmt_str = '{:20}|{:5}|{:5}|{:5}|{:10}|{:10}|{:10}|{:10}|{:10}|{:10}|{:10}\n'
+    header_str = ['Data Set', 'Steps', 'Q Tol', 'Pre T', 'Pre Comp',
+		    'Pre App', 'Iters', 'SpMV', 'QEq', 'Threads', 'Time (s)']
+    body_fmt_str = '{:20} {:5} {:5} {:5} {:10.6f} {:10.6f} {:10.6f} {:10.6f} {:10.6f} {:10} {:10.6f}\n'
+
+    params = {
+            'ensemble_type': ['0'],
+            'nsteps': ['20', '100', '500', '1000'],
+            'qeq_solver_type': ['0'],
+            'qeq_solver_q_err': ['1e-6', '1e-10'],
+#            'qeq_solver_q_err': ['1e-6', '1e-8', '1e-10', '1e-14'],
+            'pre_comp_type': ['0'],
+#            'pre_comp_type': ['0', '1', '2'],
+            'pre_comp_refactor': ['1'],
+#            'pre_comp_sweeps': ['3'],
+#            'pre_app_type': ['0'],
+#            'pre_app_jacobi_iters': ['50'],
+            'threads': ['12'],
+#            'threads': ['1', '2', '4', '12', '24'],
+            'geo_format': ['1'],
+    }
+
+    test_cases = [
+            TestCase(path.join(data_dir, 'water/water_6540.pdb'),
+                path.join(data_dir, 'water/ffield.water'),
+                path.join(control_dir, 'param.gpu.water'),
+                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+#            TestCase(path.join(data_dir, 'water/water_78480.pdb'),
+#                path.join(data_dir, 'water/ffield.water'),
+#                path.join(control_dir, 'param.gpu.water'),
+#                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+#            TestCase(path.join(data_dir, 'water/water_327000.geo'),
+#                path.join(data_dir, 'water/ffield.water'),
+#                path.join(control_dir, 'param.gpu.water'),
+#                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+            TestCase(path.join(data_dir, 'bilayer/bilayer_56800.pdb'),
+                path.join(data_dir, 'bilayer/ffield-bio'),
+                path.join(control_dir, 'param.gpu.water'),
+                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+            TestCase(path.join(data_dir, 'dna/dna_19733.pdb'),
+                path.join(data_dir, 'dna/ffield-dna'),
+                path.join(control_dir, 'param.gpu.water'),
+                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+            TestCase(path.join(data_dir, 'silica/silica_6000.pdb'),
+                path.join(data_dir, 'silica/ffield-bio'),
+                path.join(control_dir, 'param.gpu.water'),
+                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+#            TestCase(path.join(data_dir, 'silica/silica_72000.geo'),
+#                path.join(data_dir, 'silica/ffield-bio'),
+#                path.join(control_dir, 'param.gpu.water'),
+#                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+#            TestCase(path.join(data_dir, 'silica/silica_300000.geo'),
+#                path.join(data_dir, 'silica/ffield-bio'),
+#                path.join(control_dir, 'param.gpu.water'),
+#                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+            TestCase(path.join(data_dir, 'petn/petn_48256.pdb'),
+                path.join(data_dir, 'petn/ffield.petn'),
+                path.join(control_dir, 'param.gpu.water'),
+                params=params, result_header_fmt=header_fmt_str, result_header = header_str, result_body_fmt=body_fmt_str),
+                        ]
+    for test in test_cases:
+        test.run(bin_file='sPuReMD/bin/spuremd')