Fix bugs in Chemkin thermo entry element reading/writing

rmgpy/chemkin.pyx

@@ -80,7 +80,7 @@ def read_thermo_entry(entry, Tmin=0, Tint=0, Tmax=0):
 
     comment = lines[0][len(species):24].strip()
     formula = {}
-    for i in [24, 29, 34, 39, 74]:
+    for i in [24, 29, 34, 39, 73]:
         element, count = lines[0][i:i + 2].strip(), lines[0][i + 2:i + 5].strip()
         if element:
             try:
@@ -103,6 +103,21 @@ def read_thermo_entry(entry, Tmin=0, Tint=0, Tmax=0):
                         raise
             if count != 0:  # Some people put garbage elements in, with zero count. Ignore these. Allow negative counts though (eg. negative one electron)
                 formula[element] = count
+
+    # Parsing for extended elemental composition syntax, adapted from Cantera ck2cti.py
+    if lines[0].rstrip().endswith('&'):
+        complines = []
+        for i in range(len(lines)-1):
+            if lines[i].rstrip().endswith('&'):
+                complines.append(lines[i+1])
+            else:
+                break
+        lines = [lines[0]] + lines[i+1:]
+        elements = ' '.join(line.rstrip('&\n') for line in complines).split()
+        formula = {}
+        for i in range(0, len(elements), 2):
+            formula[elements[i].capitalize()] = int(elements[i+1])
+
     phase = lines[0][44]
     if phase.upper() != 'G':
         logging.warning("Was expecting gas phase thermo data for {0}. Skipping thermo data.".format(species))
@@ -1497,6 +1512,15 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
     if element_counts is None:
         element_counts = get_element_count(species.molecule[0])
 
+    # 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))
+    else:
+        sorted_elements = sorted(element_counts)
+    element_counts = {e: element_counts[e] for e in sorted_elements}
+
     string = ''
     # Write thermo comments
     if verbose:
@@ -1509,35 +1533,44 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
                 else:
                     string += "! {0}\n".format(line)
 
-    # Line 1
-    string += '{0:<16}        '.format(get_species_identifier(species))
-    if len(element_counts) <= 4:
-        # Use the original Chemkin syntax for the element counts
-        for key, count in element_counts.items():
-            if isinstance(key, tuple):
-                symbol, isotope = key
-                chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
-            else:
-                chemkin_name = key
-            string += '{0!s:<2}{1:>3d}'.format(chemkin_name, count)
-        string += '     ' * (4 - len(element_counts))
-    else:
-        string += '     ' * 4
-    string += 'G{0:>10.3f}{1:>10.3f}{2:>8.2f}      1'.format(thermo.polynomials[0].Tmin.value_si,
-                                                             thermo.polynomials[1].Tmax.value_si,
-                                                             thermo.polynomials[0].Tmax.value_si)
-    if len(element_counts) > 4:
-        string += '&\n'
+    # Compile element count string
+    extended_syntax = len(element_counts) > 4  # If there are more than 4 elements, use extended syntax
+    elements = []
+    for key, count in element_counts.items():
+        if isinstance(key, tuple):
+            symbol, isotope = key
+            chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
+        else:
+            chemkin_name = key
+        if extended_syntax:
+            # Create a list of alternating elements and counts
+            elements.extend([chemkin_name, str(count)])
+        else:
+            # Create a list of 5-column wide formatted element counts, e.g. 'C  10'
+            elements.append('{0!s:<2}{1:>3d}'.format(chemkin_name, count))
+    if extended_syntax:
         # Use the new-style Chemkin syntax for the element counts
+        # Place all elements in space delimited format on new line
         # This will only be recognized by Chemkin 4 or later
-        for key, count in element_counts.items():
-            if isinstance(key, tuple):
-                symbol, isotope = key
-                chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
-            else:
-                chemkin_name = key
-            string += '{0!s:<2}{1:>3d}'.format(chemkin_name, count)
-    string += '\n'
+        elem_1 = ' ' * 20
+        elem_2 = '&\n' + ' '.join(elements)
+    else:
+        # Use the original Chemkin syntax for the element counts
+        # Place up to 4 elements in columns 24-43 of the first line
+        # (don't use the space in columns 74-78 for the 5th element
+        #  because nobody else does and Cantera can't read it)
+        elem_1 = ''.join(elements)
+        elem_2 = ''
+
+    # Line 1
+    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,
+        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,

rmgpy/chemkinTest.py

@@ -33,15 +33,16 @@
 
 import rmgpy
 from rmgpy.chemkin import get_species_identifier, load_chemkin_file, load_transport_file, mark_duplicate_reactions, \
-    read_kinetics_entry, read_reaction_comments, read_thermo_entry, save_chemkin_file, save_species_dictionary, save_transport_file
+    read_kinetics_entry, read_reaction_comments, read_thermo_entry, save_chemkin_file, save_species_dictionary, \
+    save_transport_file, write_thermo_entry
 from rmgpy.chemkin import _remove_line_breaks, _process_duplicate_reactions
 from rmgpy.data.kinetics import LibraryReaction
 from rmgpy.exceptions import ChemkinError
 from rmgpy.kinetics.arrhenius import Arrhenius, MultiArrhenius
 from rmgpy.kinetics.chebyshev import Chebyshev
 from rmgpy.reaction import Reaction
 from rmgpy.species import Species
-from rmgpy.thermo import NASA
+from rmgpy.thermo import NASA, NASAPolynomial
 from rmgpy.transport import TransportData
 
 
@@ -428,6 +429,120 @@ def test_mark_duplicate_reactions(self):
         self.assertEqual(duplicate_flags, expected_flags)
 
 
+class TestThermoReadWrite(unittest.TestCase):
+
+    def setUp(self):
+        """This method is run once before each test."""
+        coeffs_low = [4.03055, -0.00214171, 4.90611e-05, -5.99027e-08, 2.38945e-11, -11257.6, 3.5613]
+        coeffs_high = [-0.307954, 0.0245269, -1.2413e-05, 3.07724e-09, -3.01467e-13, -10693, 22.628]
+        Tmin = 300.
+        Tmax = 3000.
+        Tint = 650.73
+        E0 = -782292.  # J/mol.
+        comment = "C2H6"
+        self.nasa = NASA(
+            polynomials=[
+                NASAPolynomial(coeffs=coeffs_low, Tmin=(Tmin, "K"), Tmax=(Tint, "K")),
+                NASAPolynomial(coeffs=coeffs_high, Tmin=(Tint, "K"), Tmax=(Tmax, "K")),
+            ],
+            Tmin=(Tmin, "K"),
+            Tmax=(Tmax, "K"),
+            E0=(E0, "J/mol"),
+            comment=comment,
+        )
+
+        # Chemkin entries for testing - note that the values are all the same
+        self.entry1 = """C2H6                    C   2H   6          G   300.000  3000.000  650.73      1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+        self.entry2 = """CH3NO2X                                     G   300.000  3000.000  650.73      1&
+C 1 H 3 N 1 O 2 X 1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+        self.entry3 = """CH3NO2SX                                    G   300.000  3000.000  650.73      1&
+C 1 H 3 N 1 O 2 S 1 X 1
+-3.07954000E-01 2.45269000E-02-1.24130000E-05 3.07724000E-09-3.01467000E-13    2
+-1.06930000E+04 2.26280000E+01 4.03055000E+00-2.14171000E-03 4.90611000E-05    3
+-5.99027000E-08 2.38945000E-11-1.12576000E+04 3.56130000E+00                   4
+"""
+
+    def test_write_thermo_block(self):
+        """Test that we can write a normal thermo block"""
+        species = Species(smiles='CC')
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry1)
+
+    def test_read_thermo_block(self):
+        """Test that we can read a normal thermo block"""
+        species, thermo, formula = read_thermo_entry(self.entry1)
+
+        self.assertEqual(species, 'C2H6')
+        self.assertEqual(formula, {'H': 6, 'C': 2})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+    def test_write_thermo_block_5_elem(self):
+        """Test that we can write a thermo block for a species with 5 elements"""
+        species = Species().from_adjacency_list("""
+1 O u0 p3 c-1 {3,S}
+2 O u0 p2 c0 {3,D}
+3 N u0 p0 c+1 {1,S} {2,D} {4,S}
+4 C u0 p0 c0 {3,S} {5,S} {6,S} {7,S}
+5 H u0 p0 c0 {4,S}
+6 H u0 p0 c0 {4,S}
+7 H u0 p0 c0 {4,S}
+8 X u0 p0 c0
+""")
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry2)
+
+    def test_read_thermo_block_5_elem(self):
+        """Test that we can read a thermo block with 5 elements"""
+        species, thermo, formula = read_thermo_entry(self.entry2)
+
+        self.assertEqual(species, 'CH3NO2X')
+        self.assertEqual(formula, {'X': 1, 'C': 1, 'O': 2, 'H': 3, 'N': 1})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+    def test_write_thermo_block_6_elem(self):
+        """Test that we can write a thermo block for a species with 6 elements"""
+        species = Species().from_adjacency_list("""
+1 O u0 p3 c-1 {2,S}
+2 N u0 p0 c+1 {1,S} {3,D} {4,S}
+3 O u0 p2 c0 {2,D}
+4 C u0 p0 c0 {2,S} {5,S} {6,S} {7,S}
+5 S u0 p2 c0 {4,S} {8,S}
+6 H u0 p0 c0 {4,S}
+7 H u0 p0 c0 {4,S}
+8 H u0 p0 c0 {5,S}
+9 X u0 p0 c0
+""")
+        species.thermo = self.nasa
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(result, self.entry3)
+
+    def test_read_thermo_block_6_elem(self):
+        """Test that we can read a thermo block with 6 elements"""
+        species, thermo, formula = read_thermo_entry(self.entry3)
+
+        self.assertEqual(species, 'CH3NO2SX')
+        self.assertEqual(formula, {'X': 1, 'C': 1, 'O': 2, 'H': 3, 'N': 1, 'S': 1})
+        self.assertTrue(self.nasa.is_identical_to(thermo))
+
+
 class TestReadReactionComments(unittest.TestCase):
     @classmethod
     def setUpClass(cls):

rmgpy/molecule/molecule.py

@@ -1144,7 +1144,9 @@ def get_formula(self):
             symbol = atom.element.symbol
             element_dict[symbol] = element_dict.get(symbol, 0) + 1
 
-        # Use the Hill system to generate the formula
+        # Use the Hill system to generate the formula.
+        # If you change this algorithm consider also updating 
+        # the chemkin.write_thermo_entry method
         formula = ''
 
         # Carbon and hydrogen always come first if carbon is present