In [ ]:

    

import rmgpy
from rmgpy.rmg.model import CoreEdgeReactionModel
from rmgpy.data.rmg import RMGDatabase
from rmgpy.rmg.react import *
from rmgpy.reaction import Reaction
from rmgpy.molecule.molecule import Molecule
from rmgpy.molecule.resonance import *
from rmgpy.species import Species
from rmgpy.thermo.thermoengine import submit
from rmgpy.kinetics.kineticsdata import KineticsData
from rmgpy.data.kinetics.family import TemplateReaction
from rmgpy.data.kinetics.depository import DepositoryReaction
from rmgpy.kinetics.arrhenius import Arrhenius

from IPython.display import display


    

In [ ]:

    

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


    

In [ ]:

    

families = ['R_Addition_MultipleBond']


    

In [ ]:

    

databasePath = rmgpy.settings['database.directory']

database = RMGDatabase()
database.load(
    path = databasePath,
    thermoLibraries = ['primaryThermoLibrary', 'C10H11'],
    reactionLibraries = [],
    seedMechanisms = [],
    kineticsFamilies = families,
)

for family in database.kinetics.families.itervalues():
    family.addKineticsRulesFromTrainingSet(thermoDatabase=database.thermo)
    family.fillKineticsRulesByAveragingUp(verbose=True)


    

In [ ]:

    

benz = Species().fromSMILES('c12ccccc1cccc2')
benz.generate_resonance_structures()
arom = benz.molecule[1]
display(arom)
keku = benz.molecule[0]
display(keku)

combos = [
    (Species(molecule=[keku]), Species().fromSMILES('[H]')),
    (Species(molecule=[arom]), Species().fromSMILES('[H]')),
    (Species(molecule=[keku]), Species().fromSMILES('[CH3]')),
    (Species(molecule=[arom]), Species().fromSMILES('[CH3]')),
    (Species(molecule=[keku]), Species().fromSMILES('[OH]')),
    (Species(molecule=[arom]), Species().fromSMILES('[OH]')),
    (Species(molecule=[keku]), Species().fromSMILES('[c]1ccccc1')),
    (Species(molecule=[arom]), Species().fromSMILES('[c]1ccccc1')),
]

for c in combos:
    submit(c[1])


    

In [ ]:

    

cerm = CoreEdgeReactionModel()

for reactants in combos:
    result = reactSpecies(reactants)
    print len(result)

    for rxn0 in result:
        rxn1 = cerm.makeNewReaction(rxn0, checkExisting=False)
    for rxn0 in cerm.newReactionList:
        cerm.applyKineticsToReaction(rxn0)
        if isinstance(rxn0.kinetics, KineticsData):
            rxn0.kinetics = reaction.kinetics.toArrhenius()
        if isinstance(rxn0,TemplateReaction) or isinstance(rxn0,DepositoryReaction):
            rxn0.fixBarrierHeight()


    

In [ ]:

    

for i, rxn0 in enumerate(cerm.newReactionList):
    print i
    display(rxn0)
    print rxn0.template
    print rxn0.kinetics


    

In [ ]:

    

# Aromatics
selected = [
    [8, 12],
#    [15, 19],
    [22, 26],
]

literature = [
    [Arrhenius(A=(2.195e-3,'m^3/(mol*s)'), n=2.88, Ea=(45655.81,'J/mol'), T0=(1,'K'))],
#    [Arrhenius(A=(2.29e12,'cm^3/(mol*s)'), n=0, Ea=(0.68,'kcal/mol'), T0=(1,'K'))],
    [Arrhenius(A=(9.5499e11,'cm^3/(mol*s)'), n=0, Ea=(4.308,'kcal/mol'), T0=(1,'K'))],
]

labels = [
    'Double',
    'Aromatic',
    'Literature',
]

filename = 'aromatic_bonds.png'


    

In [ ]:

    

for rxns in selected:
    for rxn in rxns:
        rxn0 = cerm.newReactionList[rxn]
        display(rxn0)
        print rxn0.template
        print rxn0.kinetics


    

In [ ]:

    

# plt.rcParams['font.sans-serif'] = ['Source Sans Pro']

fs = 16  # font size

pressure = 1e5 # Pa
temperature = np.linspace(298, 2000, 20)

plt.style.use('seaborn-talk')

fig, axarr = plt.subplots(1, 2, figsize=(8,3), squeeze=False)

colormap = mpl.cm.Set1

for i, indices in enumerate(selected):
    if i < 3:
        ax = axarr[0, i]
    else:
        ax = axarr[1, i - 3]
    
    for j, index in enumerate(indices + literature[i]):
        if isinstance(index, int):
            rate = cerm.newReactionList[index].kinetics
        else:
            rate = index
        kunits = rate.A.units
        print kunits
        # Evaluate kinetics
        k = []
        for t in temperature:
            if 'm^3' in kunits:
                k.append(1e6 * rate.getRateCoefficient(t, pressure))
            else:
                k.append(rate.getRateCoefficient(t, pressure))

        x = 1000 / temperature

        ax.semilogy(x, k, color=colormap(j))
    
    ax.set_xlabel('1000/T (K)', fontsize=fs)
    ax.set_ylabel('k (' + kunits + ')', fontsize=fs)
    ax.legend(labels, loc=3)

fig.tight_layout()
plt.savefig(filename, dpi=300)


    

In [ ]: