Refactor writing of elements section of Chemkin thermo entries

For species with 4 elements or less, put them in the first line
For species with 5 elements or more, put them on the next line
using the space delimited syntax.

Although the CHEMKIN manual specifies that one may put
the 5th element atom count in columns 74-78 of row 1,
in practice nobody does this, and Cantera cannot read it
if we do, so this commit stops us from doing this.
Instead if there are 5 or more elements present, we use
the extended syntax (that we were using for 6+ elements).

For more details see discussion at
Cantera/cantera#656

rmgpy/chemkin.pyx

@@ -1518,35 +1518,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,95 @@ def test_mark_duplicate_reactions(self):
         self.assertEqual(duplicate_flags, expected_flags)
 
 
+class TestThermoWrite(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,
+        )
+
+    def test_write_thermo_block(self):
+        """Test that we can write a normal thermo block"""
+        species = Species(smiles='CC')
+        species.thermo = self.nasa
+
+        expected = """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
+"""
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(expected, result)
+
+    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
+
+        expected = """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
+"""
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(expected, result)
+
+    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
+
+        expected = """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
+"""
+
+        result = write_thermo_entry(species, verbose=False)
+
+        self.assertEqual(expected, result)
+
+
 class TestReadReactionComments(unittest.TestCase):
     @classmethod
     def setUpClass(cls):