Skip to content

Commit

Permalink
Merge pull request #1636 from ReactionMechanismGenerator/ck_thermo
Browse files Browse the repository at this point in the history
Fix bugs in Chemkin thermo entry element reading/writing
  • Loading branch information
mliu49 authored Dec 2, 2019
2 parents 7e62adf + 44d7554 commit a87d79f
Show file tree
Hide file tree
Showing 3 changed files with 181 additions and 31 deletions.
89 changes: 61 additions & 28 deletions rmgpy/chemkin.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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))
Expand Down Expand Up @@ -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:
Expand All @@ -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,
Expand Down
119 changes: 117 additions & 2 deletions rmgpy/chemkinTest.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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):
Expand Down
4 changes: 3 additions & 1 deletion rmgpy/molecule/molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit a87d79f

Please sign in to comment.