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 IPython.display import display


    

In [ ]:

    

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


    

In [ ]:

    

families = ['Intra_R_Add_Endocyclic_New', 'Intra_R_Add_Endocyclic']
# 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 [ ]:

    

combos = [
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[CH2]C')),
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[CH2]C=C')),
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[CH]=C')),
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[OH]')),
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[c]1ccccc1')),
    (Species().fromSMILES('c1ccccc1'), Species().fromSMILES('[CH2]c1ccccc1')),
]

for c in combos:
    for r in c:
        submit(r)


    

In [ ]:

    

combos = [
    (Species().fromSMILES('C1=CCCC1C[CH2]'),),
    (Species().fromSMILES('C1=CCCCC1CC[CH2]'),),
    (Species().fromSMILES('c1ccccc1CC[CH2]'),),
]

for c in combos:
    for r in c:
        submit(r)


    

In [ ]:

    

for c in combos:
    for r in c:
        for m in r.molecule:
            display(m)


    

In [ ]:

    




    

In [ ]:

    

cerm = CoreEdgeReactionModel()

for reactants in combos:
    result = reactSpecies(reactants)

    for rxn0 in result:
        rxn1 = cerm.makeNewReaction(rxn0)
    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 rxn0 in cerm.newReactionList:
    display(rxn0)
    print rxn0.template
    print rxn0.kinetics


    

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, 3, figsize=(12,3), squeeze=False)

# selected = [
#     [0, 1],
#     [2, 3],
#     [4, 5],
#     [6, 7],
#     [8, 11],
#     [12, 18],
# ]

selected = [
    [1, 2],
    [5, 6],
    [11, 15],
]

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

labels = [
    'Ring',
    'No Ring',
]

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):
        rate = cerm.newReactionList[index].kinetics
        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('aromatic_bonds.png', dpi=300)
plt.savefig('ring_attribute.png', dpi=300)


    

In [ ]: