Fixes to Chemkin file writing and reading

examples/scripts/diffModels/Models/DFT_chem.inp

@@ -57,7 +57,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar  1               G   200.000  6000.000 1000.00      1

examples/scripts/diffModels/Models/Primary_chem.inp

@@ -53,7 +53,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar  1               G   200.000  6000.000 1000.00      1

ipython/data/pathway_analysis/minimal/chem.inp

@@ -20,7 +20,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar1                 G200.000   6000.000  1000.00       1

ipython/data/regression/new/chem_annotated.inp

@@ -25,7 +25,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar1                 G200.000   6000.000  1000.00       1

ipython/data/regression/old/chem_annotated.inp

@@ -25,7 +25,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar1                 G200.000   6000.000  1000.00       1

rmgpy/chemkin.pyx

@@ -52,6 +52,7 @@ from rmgpy.reaction import Reaction
 from rmgpy.rmg.pdep import PDepNetwork, PDepReaction
 from rmgpy.species import Species
 from rmgpy.thermo import NASAPolynomial, NASA
+from rmgpy.thermo.nasa cimport NASA, NASAPolynomial
 from rmgpy.transport import TransportData
 from rmgpy.util import make_output_subdirectory
 
@@ -1250,9 +1251,9 @@ def read_thermo_block(f, species_dict):
         meaningfulline, comment = remove_comment_from_line(line)
     Tmin = Tint = Tmax = None
     try:
-        Tmin = float(meaningfulline[0:9].strip())
-        Tint = float(meaningfulline[10:19].strip())
-        Tmax = float(meaningfulline[20:29].strip())
+        Tmin = float(meaningfulline[0:10].strip())
+        Tint = float(meaningfulline[10:20].strip())
+        Tmax = float(meaningfulline[20:30].strip())
         if [Tmin, Tint, Tmax] != [float(i) for i in meaningfulline.split()[0:3]]:
             logging.warning("Default temperature range line {0!r} may be badly formatted.".format(line))
             logging.warning("It should have Tmin in columns 1-10, Tmid in columns 11-20, and Tmax in columns 21-30")
@@ -1568,37 +1569,55 @@ def get_species_identifier(species):
 ################################################################################
 
 
-def write_thermo_entry(species, element_counts=None, verbose=True):
+def write_thermo_entry(species, element_counts=None, bint verbose=True):
     """
     Return a string representation of the NASA model readable by Chemkin.
     To use this method you must have exactly two NASA polynomials in your
     model, and you must use the seven-coefficient forms for each.
     """
-
-    thermo = species.get_thermo_data()
-
-    if not isinstance(thermo, NASA):
+    cdef NASA thermo
+    cdef NASAPolynomial poly_low, poly_high
+    cdef dict counts
+    cdef list sorted_elements, elements, short_lines
+    cdef bint extended_syntax
+    cdef int count, isotope
+    cdef str string, line, short_line, chemkin_name, symbol, elem_1, elem_2
+    cdef object thermo_data
+
+    thermo_data = species.get_thermo_data()
+
+    if not isinstance(thermo_data, NASA):
         raise ChemkinError('Cannot generate Chemkin string for species "{0}": '
                            'Thermodynamics data must be a NASA object.'.format(species))
+    thermo = <NASA>thermo_data
 
     assert len(thermo.polynomials) == 2
-    assert thermo.polynomials[0].Tmin.value_si < thermo.polynomials[1].Tmin.value_si
-    assert thermo.polynomials[0].Tmax.value_si == thermo.polynomials[1].Tmin.value_si
-    assert thermo.polynomials[0].cm2 == 0 and thermo.polynomials[0].cm1 == 0
-    assert thermo.polynomials[1].cm2 == 0 and thermo.polynomials[1].cm1 == 0
+    poly_low = thermo.polynomials[0]
+    poly_high = thermo.polynomials[1]
+    assert poly_low.Tmin.value_si < poly_high.Tmin.value_si
+    assert poly_low.Tmax.value_si == poly_high.Tmin.value_si
+    assert poly_low.cm2 == 0 and poly_low.cm1 == 0
+    assert poly_high.cm2 == 0 and poly_high.cm1 == 0
 
     # Determine the number of each type of element in the molecule
+    # Need to use the element's chemkin name, not the element symbol, because of isotopes.
+    # so we can't just use molecule[0].get_element_count().
     if element_counts is None:
-        element_counts = get_element_count(species.molecule[0])
+        counts = {}
+        for atom in species.molecule[0].vertices:
+            chemkin_name = atom.element.chemkin_name
+            counts[chemkin_name] = counts.get(chemkin_name, 0) + 1
+    else:
+        counts = element_counts
 
     # Sort the element_counts dictionary so that it's C, H, Al, B, Cl, D, etc.
     # if there's any C, else Al, B, Cl, D, H, if not. This is the "Hill" system
     # done by Molecule.get_formula
-    if 'C' in element_counts:
-        sorted_elements = sorted(element_counts, key = lambda e: {'C':'0','H':'1'}.get(e, e))
+    if 'C' in counts:
+        sorted_elements = sorted(counts, key=lambda e: {'C': '0', 'H': '1'}.get(e, e))
     else:
-        sorted_elements = sorted(element_counts)
-    element_counts = {e: element_counts[e] for e in sorted_elements}
+        sorted_elements = sorted(counts)
+    counts = {e: counts[e] for e in sorted_elements}
 
     string = ''
     # Write thermo comments
@@ -1613,9 +1632,9 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
                     string += "! {0}\n".format(line)
 
     # Compile element count string
-    extended_syntax = len(element_counts) > 4  # If there are more than 4 elements, use extended syntax
+    extended_syntax = len(counts) > 4  # If there are more than 4 elements, use extended syntax
     elements = []
-    for key, count in element_counts.items():
+    for key, count in counts.items():
         if isinstance(key, tuple):
             symbol, isotope = key
             chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
@@ -1645,31 +1664,31 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
     string += '{ident:<16}        {elem_1:<20}G{Tmin:>10.3f}{Tint:>10.3f}{Tmax:>8.2f}      1{elem_2}\n'.format(
         ident=get_species_identifier(species),
         elem_1=elem_1,
-        Tmin=thermo.polynomials[0].Tmin.value_si,
-        Tint=thermo.polynomials[1].Tmax.value_si,
-        Tmax=thermo.polynomials[0].Tmax.value_si,
+        Tmin=poly_low.Tmin.value_si,
+        Tint=poly_high.Tmax.value_si,
+        Tmax=poly_low.Tmax.value_si,
         elem_2=elem_2,
     )
 
     # Line 2
-    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E}    2\n'.format(thermo.polynomials[1].c0,
-                                                                                      thermo.polynomials[1].c1,
-                                                                                      thermo.polynomials[1].c2,
-                                                                                      thermo.polynomials[1].c3,
-                                                                                      thermo.polynomials[1].c4)
+    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E}    2\n'.format(poly_high.c0,
+                                                                                      poly_high.c1,
+                                                                                      poly_high.c2,
+                                                                                      poly_high.c3,
+                                                                                      poly_high.c4)
 
     # Line 3
-    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E}    3\n'.format(thermo.polynomials[1].c5,
-                                                                                      thermo.polynomials[1].c6,
-                                                                                      thermo.polynomials[0].c0,
-                                                                                      thermo.polynomials[0].c1,
-                                                                                      thermo.polynomials[0].c2)
+    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E}    3\n'.format(poly_high.c5,
+                                                                                      poly_high.c6,
+                                                                                      poly_low.c0,
+                                                                                      poly_low.c1,
+                                                                                      poly_low.c2)
 
     # Line 4
-    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}                   4\n'.format(thermo.polynomials[0].c3,
-                                                                                          thermo.polynomials[0].c4,
-                                                                                          thermo.polynomials[0].c5,
-                                                                                          thermo.polynomials[0].c6)
+    string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}                   4\n'.format(poly_low.c3,
+                                                                                          poly_low.c4,
+                                                                                          poly_low.c5,
+                                                                                          poly_low.c6)
 
     return string
 
@@ -2198,7 +2217,7 @@ def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_
 
     # Thermodynamics section
     f.write('THERM ALL\n')
-    f.write('    300.000  1000.000  5000.000\n\n')
+    f.write('   300.000  1000.000  5000.000\n\n')
     for spec in sorted_species:
         f.write(write_thermo_entry(spec, verbose=verbose))
         f.write('\n')

rmgpy/molecule/element.py

@@ -77,7 +77,7 @@ def __init__(self, number, symbol, name, mass, isotope=-1, chemkin_name=None):
         self.name = name
         self.mass = mass
         self.isotope = isotope
-        self.chemkin_name = chemkin_name or self.name
+        self.chemkin_name = chemkin_name or self.symbol
         if symbol in {'X','L','R','e'}:
             self.cov_radius = 0
         else:

rmgpy/solver/files/collider_model/chem.inp

@@ -20,7 +20,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 Ar                      Ar1                 G200.000   6000.000  1000.00       1
  2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2

rmgpy/solver/files/listener/chemkin/chem.inp

@@ -15,7 +15,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 Ar                      Ar1                 G200.000   6000.000  1000.00       1
  2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2

rmgpy/tools/data/diffmodels/chem1.inp

@@ -19,7 +19,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ethane(1)               C 2  H 6            G100.000   5000.000  954.52        1
  4.58991157E+00 1.41506360E-02-4.75954111E-06 8.60275134E-10-6.21700682E-14    2

rmgpy/tools/data/diffmodels/chem2.inp

@@ -19,7 +19,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ethane(1)               C 2  H 6            G100.000   5000.000  954.52        1
  4.58991157E+00 1.41506360E-02-4.75954111E-06 8.60275134E-10-6.21700682E-14    2

rmgpy/tools/data/diffmodels/chem3.inp

@@ -12,7 +12,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 ! GRI-Mech3.0
 ethane                  H   6C   2          G   100.000  5000.000 1002.57      1
  2.56122991E+00 1.83185519E-02-7.41959133E-06 1.38553537E-09-9.75184837E-14    2

rmgpy/tools/data/diffmodels/surf_model/chem_surface1.inp

@@ -21,7 +21,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: surfaceThermoPt111 Binding energy corrected by LSR () from Pt111
 X(1)                    X   1               G   100.000  5000.000 1554.83      1

rmgpy/tools/data/diffmodels/surf_model/chem_surface2.inp

@@ -20,7 +20,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: surfaceThermoPt111 Binding energy corrected by LSR () from Pt111
 X(1)                    X   1               G   100.000  5000.000 1554.81      1

rmgpy/tools/data/flux/chemkin/chem.inp

@@ -17,7 +17,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 Ar                      Ar1                 G200.000   6000.000  1000.00       1
  2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2

rmgpy/tools/data/regression/benchmark/chem_annotated.inp

@@ -25,7 +25,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar1                 G200.000   6000.000  1000.00       1

rmgpy/tools/data/regression/tested/chem_annotated.inp

@@ -25,7 +25,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Thermo library: primaryThermoLibrary
 Ar                      Ar1                 G200.000   6000.000  1000.00       1

rmgpy/tools/data/sim/liquid/chem.inp

@@ -49,7 +49,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 Ar                      Ar  1               G   100.000  5000.000 3142.77      1
  2.50000000E+00 3.32026413E-12-1.37161245E-15 2.48094222E-19-1.65799619E-23    2

rmgpy/tools/data/sim/mbSampled/chem.inp

@@ -37,7 +37,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! Narayanaswamy
 C6H5                    H   5C   6          G   100.000  5000.000  845.04      1

rmgpy/tools/data/sim/simple/chem.inp

@@ -19,7 +19,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ethane(1)               C 2  H 6            G100.000   5000.000  954.52        1
  4.58991157E+00 1.41506360E-02-4.75954111E-06 8.60275134E-10-6.21700682E-14    2

rmgpy/tools/data/various_kinetics/chem_annotated.inp

@@ -19,7 +19,7 @@ END
 
 
 THERM ALL
-    300.000  1000.000  5000.000
+   300.000  1000.000  5000.000
 
 ! 
 ! Thermo group additivity estimation: group(Cs-CsHHH) + gauche(Cs(CsRRR)) + other(R) + group(Cs-CsHHH) + gauche(Cs(CsRRR)) + other(R)