diff --git a/pychemkin/chemkin.py b/pychemkin/chemkin.py
index c10ad29..c159f4d 100644
--- a/pychemkin/chemkin.py
+++ b/pychemkin/chemkin.py
@@ -184,7 +184,8 @@ def postprocess(self, sens=False, rop=False, all=False, transpose=True):
     def writeInputHomogeneousBatch(self,problemType, reactants, temperature, pressure, endTime, 
                       Continuations=False, typeContinuation = None, Tlist = [], Plist = [],
                       variableVolume=False, variableVolumeProfile = None, 
-                      solverTimeStepProfile = None, sensitivity=[], rop=[]):
+                      solverTimeStepProfile = None, sensitivity=[], rop=[],
+                      atol=1e-20, rtol=1e-8):
         """
         Write input file for homogeneous batch reactor
         """
@@ -260,10 +261,10 @@ def writeInputHomogeneousBatch(self,problemType, reactants, temperature, pressur
 ADAP   ! Save Additional Adaptive Points
 ASTEPS 20   ! Use Solver Integration Steps
 ATLS 1.0E-6   ! Sensitivity Absolute Tolerance
-ATOL 1.0E-20   ! Absolute Tolerance
+ATOL {0}   ! Absolute Tolerance
 MAXIT 4   ! Maximum Number of Iterations
 RTLS 0.0001   ! Sensitivity Relative Tolerance
-RTOL 1.0E-8   ! Relative Tolerance""")
+RTOL {1}   ! Relative Tolerance""".format(atol, rtol))
         
         if solverTimeStepProfile:
             with open(solverTimeStepProfile, 'rb') as csvfile2:
@@ -321,7 +322,8 @@ def writeInputHomogeneousBatch(self,problemType, reactants, temperature, pressur
     def writeInputPlugFlow(self,problemType, reactants, 
                            startingAxialPosition, endingAxialPosition, diameter, momentum=True, massflowrate = None, sccmflowrate = None, 
                            temperature  = None, pressure  = None,  
-                           temperatureProfile = None, pressureProfile = None, sensitivity=[], rop=[]):
+                           temperatureProfile = None, pressureProfile = None, sensitivity=[], rop=[],
+                           atol=1e-9, rtol=1e-6):
     
         """
         Write input file for typical plug flow
@@ -427,14 +429,14 @@ def writeInputPlugFlow(self,problemType, reactants,
 ACHG 0.0   ! Maximum Absolute Change in Site Fractions
 ADAP   ! Save Additional Adaptive Points
 ATLS 1.0E-6   ! Sensitivity Absolute Tolerance
-ATOL 1.0E-9   ! Absolute Tolerance
+ATOL {0}  ! Absolute Tolerance
 MAXIT 4   ! Maximum Number of Iterations
 NNEG   ! Force Non-negative Solution
 PSV 1.0E-8   ! Scaling Factor for Relaxing Surface Equations (cm/sec)
 RCHG 1.0E-6   ! Maximum Relative Change in Site Fractions
 RTLS 0.0001   ! Sensitivity Relative Tolerance
-RTOL 1.0E-6   ! Relative Tolerance
-TSTP 1.0   ! Initial Integration Step (cm""")          
+RTOL {1}   ! Relative Tolerance
+TSTP 1.0   ! Initial Integration Step (cm""".format(atol, rtol))          
 
         input_stream+=("""
 ! 
@@ -638,21 +640,22 @@ def getMoleFraction(ckcsvFile, species=[]):
         for row in reader:
             label = row[0].strip()
             tokens = label.split('_')
-            if tokens[0] == 'Time':
+            if label.startswith('Time'):
                 tdata = numpy.array([float(r) for r in row[2:]], numpy.float)
                 tunits = row[1].strip()[1:-1].lower()
                 tdata *= {'sec': 1.0, 'min': 60., 'hr': 3600., 'msec': 1e-3, 'microsec': 1e-6}[tunits]
                 timeData.append(tdata)
             
-            if tokens[0] == 'Distance':      
+            if label.startswith('Distance'):      
                 ddata = numpy.array([float(r) for r in row[2:]], numpy.float)
                 dunits = row[1].strip()[1:-1].lower()
                 ddata *= {'cm': 1.0, 'mm': 0.1, 'm': 100.}[dunits]
                 distanceData.append(ddata)
 
-            if label.startswith('Mole_fraction'):
+            if label.startswith('Mole_fraction_'):
+                species_label = label[14:]
                 for spec in species:
-                    if tokens[2] == spec:
+                    if species_label == spec:
                         specData[spec].append(numpy.array([float(r) for r in row[2:]], numpy.float))
     if timeData:
         return timeData, specData
@@ -661,6 +664,41 @@ def getMoleFraction(ckcsvFile, species=[]):
     else:
         return
 
+def getLumpedMoleFraction(ckcsvFile, lumped_species_dict={}):
+    """
+    This method is trying to aggregate MoleFraction for lumped species
+    the input variable `lumped_species_dict`: 
+        key is name by user e.g., Undecene, 
+        value is a list of all the lumped species, 
+        e.g., ['C=CCCCCCCCCC','CC=CCCCCCCCC','CCC=CCCCCCCC','CCCC=CCCCCCC','CCCCC=CCCCCC']
+    The output is same as getMoleFraction() except that the keys of lumpedSpecData are 
+    the keys of lumped_species_dict
+    """
+    species_list = []
+    for lumped_key in lumped_species_dict:
+        lumped_species = lumped_species_dict[lumped_key] 
+        if not isinstance(lumped_species, list):
+            lumped_species = [lumped_species]
+        species_list.extend(lumped_species)
+
+    leadData, specData = getMoleFraction(ckcsvFile, species_list)
+
+    lumpedSpecData = {}
+
+    for lumped_key in lumped_species_dict:
+        lumped_species = lumped_species_dict[lumped_key]
+        if not isinstance(lumped_species, list):
+            lumpedSpecData[lumped_key] = specData[lumped_species]
+        else:
+            molfrac_sum = specData[lumped_species[0]]
+            for spc in lumped_species[1:]:
+                molfrac_spc = specData[spc]
+                sln_num = len(molfrac_spc)
+                for sln in range(sln_num):
+                    molfrac_sum[sln] = numpy.add(molfrac_sum[sln], molfrac_spc[sln])
+            lumpedSpecData[lumped_key] = molfrac_sum
+    return leadData, lumpedSpecData
+
 def getTotalMoles(ckcsvFile):
     """
     Return the time and total moles data from the given CKCSV file.