In [ ]:
import os.path
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
%matplotlib inline
import rmgpy
from rmgpy.data.rmg import RMGDatabase
from rmgpy.species import Species
from rmgpy.chemkin import readThermoEntry
In [ ]:
databasePath = rmgpy.settings['database.directory']
database = RMGDatabase()
database.load(
path = databasePath,
thermoLibraries = [],
reactionLibraries = [],
seedMechanisms = [],
kineticsFamilies = 'none'
)
In [ ]:
spec = Species().fromSMILES("C#C")
spec.generateResonanceIsomers()
spec.thermo = database.thermo.getThermoData(spec)
display(spec)
print spec.thermo.comment
In [ ]:
# NASA polynomials
data = [
"""C2HBr T04/04C 2.H 1.BR 1. 0.G 200.000 6000.000 1000. 1
6.55399311E+00 3.37962726E-03-1.18362410E-06 1.87797808E-10-1.11059116E-14 2
3.17495713E+04-8.20269727E+00 1.10795098E+00 3.21065018E-02-6.02244383E-05 3
5.45400888E-08-1.86034151E-11 3.26428366E+04 1.67414085E+01 3.39671249E+04 4
""",
"""C2HCl T05/08C 2.H 1.CL 1. 0.G 200.000 6000.000 1000. 1
6.52865585E+00 3.32425623E-03-1.14637403E-06 1.79972218E-10-1.05639468E-14 2
2.51378884E+04-9.16499932E+00 1.25077097E+00 3.10939695E-02-5.78728028E-05 3
5.20651866E-08-1.76611780E-11 2.59985454E+04 1.50044210E+01 2.73367422E+04 4
""",
"""C2HI H 1C 2I 1 G 100.000 5000.000 876.60 1
6.34079734E+00 3.80618529E-03-1.55706912E-06 2.74740006E-10-1.80899125E-14 2
3.12245678E+04-6.32368449E+00 2.57050130E+00 2.10098737E-02-3.09944354E-05 3
2.26616327E-08-6.40248787E-12 3.18855950E+04 1.13693573E+01 4
"""
]
data2 = [
"""CL J 6/82CL 1 0 0 0G 200.000 6000.000 1000. 1
2.94658358E+00-3.85985408E-04 1.36139388E-07-2.17032923E-11 1.28751025E-15 2
1.36970327E+04 3.11330136E+00 2.26062480E+00 1.54154399E-03-6.80283622E-07 3
-1.59972975E-09 1.15416636E-12 1.38552986E+04 6.57020799E+00 1.45891941E+04 4
""",
"""BR J 6/82BR 1 0 0 0G 200.000 6000.000 1000. 1
0.20866945E+01 0.71459733E-03-0.27080691E-06 0.41519029E-10-0.23016335E-14 2
0.12857696E+05 0.90837335E+01 0.24820782E+01 0.18570465E-03-0.64313029E-06 3
0.84642045E-09-0.30137068E-12 0.12709455E+05 0.68740409E+01 0.13453589E+05 4
""",
"""I J 6/82I 1 0 0 0G 200.000 6000.000 1000. 1
2.61667712E+00-2.66010320E-04 1.86060150E-07-3.81927472E-11 2.52036053E-15 2
1.20582790E+04 6.87896653E+00 2.50041683E+00-4.48046831E-06 1.69962536E-08 3
-2.67708030E-11 1.48927452E-14 1.20947990E+04 7.49816581E+00 1.28402035E+04 4
"""
]
In [ ]:
labels = [
'H',
'Br',
'Cl',
'I',
]
labels = [
'Cl',
'Br',
'I',
]
In [ ]:
thermo = [spec.thermo]
for entry in data:
thermo.append(readThermoEntry(entry)[1])
thermo2 = [readThermoEntry(entry)[1] for entry in data2]
In [ ]:
tlist = np.linspace(300, 3000, 100)
Cpall = []
Hall = []
Sall = []
Gall = []
for entry in thermo2:
Cplist = np.zeros_like(tlist)
Hlist = np.zeros_like(tlist)
Slist = np.zeros_like(tlist)
Glist = np.zeros_like(tlist)
for i, t in enumerate(tlist):
Cplist[i] = entry.getHeatCapacity(t) / 4.184
Hlist[i] = entry.getEnthalpy(t) / 4184
Slist[i] = entry.getEntropy(t) / 4.184
Glist[i] = entry.getFreeEnergy(t) / 4184
Cpall.append(Cplist)
Hall.append(Hlist)
Sall.append(Slist)
Gall.append(Glist)
In [ ]:
for item, label in zip(Cpall, labels):
plt.plot(tlist, item, label=label)
plt.xlabel('Temperature (K)')
plt.ylabel('Heat Capacity (cal/mol-K)')
plt.legend(loc=4)
plt.show()
In [ ]:
for item, label in zip(Hall, labels):
plt.plot(tlist, item, label=label)
plt.xlabel('Temperature (K)')
plt.ylabel('Enthalpy (kcal/mol-K)')
plt.legend(loc=2)
plt.show()
In [ ]:
for item, label in zip(Sall, labels):
plt.plot(tlist, item, label=label)
plt.xlabel('Temperature (K)')
plt.ylabel('Entropy (cal/mol-K)')
plt.legend(loc=2)
plt.show()
In [ ]:
for item, label in zip(Gall, labels):
plt.plot(tlist, item, label=label)
plt.xlabel('Temperature (K)')
plt.ylabel('Free Energy (kcal/mol-K)')
plt.legend(loc=3)
plt.show()
In [ ]: