Chemkin was not writing isotopes correctly.

In this patch, we skip using the Molecule.get_element_count()
because that returns just the element.symbol, and deuterium is H.

Author: rwest

rmgpy/chemkin.pyx

@@ -1588,8 +1588,13 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
     assert thermo.polynomials[1].cm2 == 0 and thermo.polynomials[1].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])
+        element_counts = {}
+        for atom in species.molecule[0].atoms:
+            element = atom.element.chemkin_name
+            element_counts[element] = element_counts.get(element, 0) + 1
 
     # 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

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:

test/rmgpy/chemkinTest.py

@@ -791,6 +791,20 @@ def test_read_thermo_block(self):
         assert formula == {"H": 6, "C": 2}
         assert self.nasa.is_identical_to(thermo)
 
+    def test_write_thermo_block_for_isotope_uses_chemkin_name(self):
+        """Isotopic atoms should use their Chemkin names in thermo composition."""
+        deuterium = Species().from_adjacency_list(
+            "1 H u0 p0 c0 i2 {2,S}\n"
+            "2 H u0 p0 c0 {1,S}"
+        )
+        deuterium.thermo = self.nasa
+
+        result = write_thermo_entry(deuterium, verbose=False)
+
+        first_line = result.splitlines()[0]
+        assert "D   1" in first_line
+        assert "H   1" in first_line
+
     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(