In [ ]:

    

import os.path
import rmgpy
from IPython.display import display
from rmgpy.data.thermo import ThermoDatabase
from rmgpy.chemkin import loadSpeciesDictionary
from rmgpy.thermo.thermoengine import submit
from rmgpy.cnn_framework.predictor import Predictor


    

In [ ]:

    

temperature = 298


    

In [ ]:

    

databasePath = os.path.join(rmgpy.settings['database.directory'], 'thermo')

database = ThermoDatabase()
database.load(
    path = databasePath,
    libraries = ['PAHLibrary'],
    )


    

In [ ]:

    

model = '/home/mjliu/Code/RMG-Py/examples/cnn/evaluate/test_model'

h298_predictor = Predictor()

predictor_input = os.path.join(model,
                              'predictor_input.py')

h298_predictor.load_input(predictor_input)

param_path = os.path.join(model,
                         'saved_model',
                         'full_train.h5')
h298_predictor.load_parameters(param_path)


    

In [ ]:

    

speciesList = loadSpeciesDictionary('/home/mjliu/Documents/PAHThermo/species_dictionary.txt')


    

In [ ]:

    

labels = []
gavThermo = {}
gavEnthalpy = []
gavEntropy = []
calcThermo = {}
calcEnthalpy = []
calcEntropy = []
cnnEnthalpy = []

for label in speciesList.iterkeys():
    labels.append(label)
    spec = speciesList[label]
    display(spec)
    
    gavThermo[label] = database.getThermoDataFromGroups(spec)
    
    gavEnthalpy.append(gavThermo[label].getEnthalpy(temperature) / 4184)
    gavEntropy.append(gavThermo[label].getEntropy(temperature) / 4.184)

    calcThermo[label] = database.getThermoData(spec)
    
    calcEnthalpy.append(calcThermo[label].getEnthalpy(temperature) / 4184)
    calcEntropy.append(calcThermo[label].getEntropy(temperature) / 4.184)
    
    cnnEnthalpy.append(h298_predictor.predict(spec.molecule[0]))


    

In [ ]:

    

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

from bokeh.plotting import figure, show
from bokeh.io import output_notebook
output_notebook()


    

In [ ]:

    

p = figure(plot_width=500, plot_height=400)

# add a square renderer with a size, color, and alpha
p.circle(gavEnthalpy, calcEnthalpy, size=10, color="green", alpha=0.5)

x = np.array([0, 200])
y = x
p.line(x=x, y=y, line_width=2, color='#636363')
p.line(x=x, y=y+10, line_width=2,line_dash="dashed", color='#bdbdbd')
p.line(x=x, y=y-10, line_width=2, line_dash="dashed", color='#bdbdbd')

p.xaxis.axis_label = "H298 GAV (kcal/mol/K)"
p.yaxis.axis_label = "H298 G3(MP2)//B3 (kcal/mol/K)"
p.xaxis.axis_label_text_font_style = "normal"
p.yaxis.axis_label_text_font_style = "normal"
p.xaxis.axis_label_text_font_size = "16pt"
p.yaxis.axis_label_text_font_size = "16pt"
p.xaxis.major_label_text_font_size = "12pt"
p.yaxis.major_label_text_font_size = "12pt"

show(p)


    

In [ ]:

    

p = figure(plot_width=500, plot_height=400)

# add a square renderer with a size, color, and alpha
p.circle(cnnEnthalpy[:15], calcEnthalpy[:15], size=10, color="green", alpha=0.5)

x = np.array([0, 200])
y = x
p.line(x=x, y=y, line_width=2, color='#636363')
p.line(x=x, y=y+10, line_width=2,line_dash="dashed", color='#bdbdbd')
p.line(x=x, y=y-10, line_width=2, line_dash="dashed", color='#bdbdbd')

p.xaxis.axis_label = "H298 CNN (kcal/mol/K)"
p.yaxis.axis_label = "H298 G3(MP2)//B3 (kcal/mol/K)"
p.xaxis.axis_label_text_font_style = "normal"
p.yaxis.axis_label_text_font_style = "normal"
p.xaxis.axis_label_text_font_size = "16pt"
p.yaxis.axis_label_text_font_size = "16pt"
p.xaxis.major_label_text_font_size = "12pt"
p.yaxis.major_label_text_font_size = "12pt"

show(p)


    

In [ ]:

    

p = figure(plot_width=500, plot_height=400)

# add a square renderer with a size, color, and alpha
p.circle(gavEntropy, calcEntropy, size=10, color="green", alpha=0.5)

x = np.array([0, 140])
y = x
p.line(x=x, y=y, line_width=2, color='#636363')
p.line(x=x, y=y+10, line_width=2,line_dash="dashed", color='#bdbdbd')
p.line(x=x, y=y-10, line_width=2, line_dash="dashed", color='#bdbdbd')

p.xaxis.axis_label = "S298 GAV (cal/mol/K)"
p.yaxis.axis_label = "S298 G3(MP2)//B3 (cal/mol/K)"
p.xaxis.axis_label_text_font_style = "normal"
p.yaxis.axis_label_text_font_style = "normal"
p.xaxis.axis_label_text_font_size = "16pt"
p.yaxis.axis_label_text_font_size = "16pt"
p.xaxis.major_label_text_font_size = "12pt"
p.yaxis.major_label_text_font_size = "12pt"

show(p)


    

In [ ]:

    

hDiff = np.subtract(calcEnthalpy,gavEnthalpy)
plt.plot(np.sort(hDiff),'bo')
plt.show()


    

In [ ]:

    

print np.mean(np.absolute(np.subtract(calcEnthalpy,gavEnthalpy)))
print np.mean(np.absolute(np.subtract(calcEntropy,gavEntropy)))


    

In [ ]:

    

sDiff = np.subtract(calcEntropy,gavEntropy)
plt.plot(np.sort(sDiff),'bo')
plt.show()


    

In [ ]:

    

plt.plot(hDiff, sDiff, 'bo')
plt.xlabel('Enthalpy Error')
plt.ylabel('Entropy Error')
plt.show()


    

In [ ]:

    

indices = np.where(sDiff>120)[0]
for index in indices:
    print labels[index]
    #print speciesList1[labels[index]].thermo.comment


    

In [ ]:

    

indices = np.where(np.logical_and(sDiff<40,sDiff>20))[0]
for index in indices:
    print labels[index]
    #print speciesList1[labels[index]].thermo.comment


    

In [ ]:

    

indices = np.where(np.logical_and(sDiff<20,sDiff>-20))[0]
for index in indices:
    print labels[index]
    #print speciesList1[labels[index]].thermo.comment


    

In [ ]:

    

indices = np.where(sDiff<-20)[0]
for index in indices:
    print labels[index]
    #print speciesList1[labels[index]].thermo.comment


    

In [ ]: