In [ ]:
import os
from collections import OrderedDict
import cantera as ct
import numpy as np
import pandas as pd

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

from rmgpy.species import Species

In [ ]:
def simulate(model, T, P, mol_frac, time):

    model.TPX = T, P, mol_frac

    reactor = ct.IdealGasReactor(model)
    
    simulation = ct.ReactorNet([reactor])

    simulation.advance(time)

    conditions = pd.Series()
    conditions['time'] = simulation.time
    conditions['temperature0'] = T
    conditions['temperature'] = model.T
    conditions['pressure0'] = P
    conditions['pressure'] = model.P
        
    species = pd.Series()
    for key, value in model.mole_fraction_dict().iteritems():
        species[key] = value
        
    return conditions, species

In [ ]:
def loadSpeciesDictionary(path):
    speciesDict = OrderedDict()
    with open(path, 'r') as f:
        adjlist = ''
        for line in f:
            if line.strip() == '' and adjlist.strip() != '':
                # Finish this adjacency list
                species = Species().fromAdjacencyList(adjlist)
                #species.generateResonanceIsomers()
                label = species.label
                speciesDict[label] = species
                adjlist = ''
            else:
                if "InChI" in line:
                    line = line.split()[0] + '\n'
                if '//' in line:
                    index = line.index('//')
                    line = line[0:index]
                adjlist += line
        else: #reach end of file
            if adjlist.strip() != '':
                species = Species().fromAdjacencyList(adjlist)
                #species.generateResonanceIsomers()
                label = species.label
                speciesDict[label] = species

    return speciesDict

In [ ]:
directory = '/Users/mjliu/Dropbox (MIT)/Research/Models/Naphthalene/RMG_sim/acetylene/run10'
filename = 'chem_mod'
chemkin_file = os.path.join(directory, filename + '.inp')
cantera_file = os.path.join(directory, filename + '.cti')

In [ ]:
species_dict = loadSpeciesDictionary(os.path.join(directory, 'species_dictionary.txt'))

In [ ]:
# Convert chemkin to cantera
from cantera import ck2cti
if os.path.exists(cantera_file):
    raise Exception('File already exists')
ck2cti.Parser().convertMech(chemkin_file, outName=cantera_file)

In [ ]:
model = ct.Solution(cantera_file)

In [ ]:
conditions_df = pd.DataFrame()
species_df = pd.DataFrame()

T = 1200
P = 0.5 * 101325
mol_frac = {'C2H2(2)':0.8, 'A2I1':0.2}

conditions, species = simulate(model, T, P, mol_frac, 0.0001)

In [ ]:
weights = []
for label, _ in species.iteritems():
    try:
        weights.append(1000*species_dict[label].molecule[0].getMolecularWeight())
    except KeyError:
        weights.append(0)
        print label

In [ ]:
species.sort_values(ascending=False)[0:10]

In [ ]:
plt.style.use('seaborn-poster')
#plt.rcParams['axes.labelsize'] = 24
#plt.rcParams['xtick.labelsize'] = 20
#plt.rcParams['ytick.labelsize'] = 20

fig = plt.figure()

plt.stem(weights, species, markerfmt=' ')

plt.xlabel('Mass')
plt.ylabel('Abundance')

In [ ]: