Skip to content

Commit 8a0b7fe

Browse files
rwestrwest
authored
Merge pull request #2955 to fix Chemkin I/O
Fixes and optimizations to Chemkin file writing and reading
2 parents 7d729a8 + fc1c6bb commit 8a0b7fe

29 files changed

Lines changed: 130 additions & 66 deletions

File tree

‎examples/scripts/diffModels/Models/DFT_chem.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ END
5757

5858

5959
THERM ALL
60-
300.000 1000.000 5000.000
60+
300.000 1000.000 5000.000
6161

6262
! Thermo library: primaryThermoLibrary
6363
Ar Ar 1 G 200.000 6000.000 1000.00 1

‎examples/scripts/diffModels/Models/Primary_chem.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,7 +53,7 @@ END
5353

5454

5555
THERM ALL
56-
300.000 1000.000 5000.000
56+
300.000 1000.000 5000.000
5757

5858
! Thermo library: primaryThermoLibrary
5959
Ar Ar 1 G 200.000 6000.000 1000.00 1

‎ipython/data/pathway_analysis/minimal/chem.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ END
2020

2121

2222
THERM ALL
23-
300.000 1000.000 5000.000
23+
300.000 1000.000 5000.000
2424

2525
! Thermo library: primaryThermoLibrary
2626
Ar Ar1 G200.000 6000.000 1000.00 1

‎ipython/data/regression/new/chem_annotated.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ END
2525

2626

2727
THERM ALL
28-
300.000 1000.000 5000.000
28+
300.000 1000.000 5000.000
2929

3030
! Thermo library: primaryThermoLibrary
3131
Ar Ar1 G200.000 6000.000 1000.00 1

‎ipython/data/regression/old/chem_annotated.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ END
2525

2626

2727
THERM ALL
28-
300.000 1000.000 5000.000
28+
300.000 1000.000 5000.000
2929

3030
! Thermo library: primaryThermoLibrary
3131
Ar Ar1 G200.000 6000.000 1000.00 1

‎rmgpy/chemkin.pyx‎

Lines changed: 56 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ from rmgpy.reaction import Reaction
5252
from rmgpy.rmg.pdep import PDepNetwork, PDepReaction
5353
from rmgpy.species import Species
5454
from rmgpy.thermo import NASAPolynomial, NASA
55+
from rmgpy.thermo.nasa cimport NASA, NASAPolynomial
5556
from rmgpy.transport import TransportData
5657
from rmgpy.util import make_output_subdirectory
5758

@@ -1250,9 +1251,9 @@ def read_thermo_block(f, species_dict):
12501251
meaningfulline, comment = remove_comment_from_line(line)
12511252
Tmin = Tint = Tmax = None
12521253
try:
1253-
Tmin = float(meaningfulline[0:9].strip())
1254-
Tint = float(meaningfulline[10:19].strip())
1255-
Tmax = float(meaningfulline[20:29].strip())
1254+
Tmin = float(meaningfulline[0:10].strip())
1255+
Tint = float(meaningfulline[10:20].strip())
1256+
Tmax = float(meaningfulline[20:30].strip())
12561257
if [Tmin, Tint, Tmax] != [float(i) for i in meaningfulline.split()[0:3]]:
12571258
logging.warning("Default temperature range line {0!r} may be badly formatted.".format(line))
12581259
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):
15681569
################################################################################
15691570

15701571

1571-
def write_thermo_entry(species, element_counts=None, verbose=True):
1572+
def write_thermo_entry(species, element_counts=None, bint verbose=True):
15721573
"""
15731574
Return a string representation of the NASA model readable by Chemkin.
15741575
To use this method you must have exactly two NASA polynomials in your
15751576
model, and you must use the seven-coefficient forms for each.
15761577
"""
1577-
1578-
thermo = species.get_thermo_data()
1579-
1580-
if not isinstance(thermo, NASA):
1578+
cdef NASA thermo
1579+
cdef NASAPolynomial poly_low, poly_high
1580+
cdef dict counts
1581+
cdef list sorted_elements, elements, short_lines
1582+
cdef bint extended_syntax
1583+
cdef int count, isotope
1584+
cdef str string, line, short_line, chemkin_name, symbol, elem_1, elem_2
1585+
cdef object thermo_data
1586+
1587+
thermo_data = species.get_thermo_data()
1588+
1589+
if not isinstance(thermo_data, NASA):
15811590
raise ChemkinError('Cannot generate Chemkin string for species "{0}": '
15821591
'Thermodynamics data must be a NASA object.'.format(species))
1592+
thermo = <NASA>thermo_data
15831593

15841594
assert len(thermo.polynomials) == 2
1585-
assert thermo.polynomials[0].Tmin.value_si < thermo.polynomials[1].Tmin.value_si
1586-
assert thermo.polynomials[0].Tmax.value_si == thermo.polynomials[1].Tmin.value_si
1587-
assert thermo.polynomials[0].cm2 == 0 and thermo.polynomials[0].cm1 == 0
1588-
assert thermo.polynomials[1].cm2 == 0 and thermo.polynomials[1].cm1 == 0
1595+
poly_low = thermo.polynomials[0]
1596+
poly_high = thermo.polynomials[1]
1597+
assert poly_low.Tmin.value_si < poly_high.Tmin.value_si
1598+
assert poly_low.Tmax.value_si == poly_high.Tmin.value_si
1599+
assert poly_low.cm2 == 0 and poly_low.cm1 == 0
1600+
assert poly_high.cm2 == 0 and poly_high.cm1 == 0
15891601

15901602
# Determine the number of each type of element in the molecule
1603+
# Need to use the element's chemkin name, not the element symbol, because of isotopes.
1604+
# so we can't just use molecule[0].get_element_count().
15911605
if element_counts is None:
1592-
element_counts = get_element_count(species.molecule[0])
1606+
counts = {}
1607+
for atom in species.molecule[0].vertices:
1608+
chemkin_name = atom.element.chemkin_name
1609+
counts[chemkin_name] = counts.get(chemkin_name, 0) + 1
1610+
else:
1611+
counts = element_counts
15931612

15941613
# Sort the element_counts dictionary so that it's C, H, Al, B, Cl, D, etc.
15951614
# if there's any C, else Al, B, Cl, D, H, if not. This is the "Hill" system
15961615
# done by Molecule.get_formula
1597-
if 'C' in element_counts:
1598-
sorted_elements = sorted(element_counts, key = lambda e: {'C':'0','H':'1'}.get(e, e))
1616+
if 'C' in counts:
1617+
sorted_elements = sorted(counts, key=lambda e: {'C': '0', 'H': '1'}.get(e, e))
15991618
else:
1600-
sorted_elements = sorted(element_counts)
1601-
element_counts = {e: element_counts[e] for e in sorted_elements}
1619+
sorted_elements = sorted(counts)
1620+
counts = {e: counts[e] for e in sorted_elements}
16021621

16031622
string = ''
16041623
# Write thermo comments
@@ -1613,9 +1632,9 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
16131632
string += "! {0}\n".format(line)
16141633

16151634
# Compile element count string
1616-
extended_syntax = len(element_counts) > 4 # If there are more than 4 elements, use extended syntax
1635+
extended_syntax = len(counts) > 4 # If there are more than 4 elements, use extended syntax
16171636
elements = []
1618-
for key, count in element_counts.items():
1637+
for key, count in counts.items():
16191638
if isinstance(key, tuple):
16201639
symbol, isotope = key
16211640
chemkin_name = get_element(symbol, isotope=isotope).chemkin_name
@@ -1645,31 +1664,31 @@ def write_thermo_entry(species, element_counts=None, verbose=True):
16451664
string += '{ident:<16} {elem_1:<20}G{Tmin:>10.3f}{Tint:>10.3f}{Tmax:>8.2f} 1{elem_2}\n'.format(
16461665
ident=get_species_identifier(species),
16471666
elem_1=elem_1,
1648-
Tmin=thermo.polynomials[0].Tmin.value_si,
1649-
Tint=thermo.polynomials[1].Tmax.value_si,
1650-
Tmax=thermo.polynomials[0].Tmax.value_si,
1667+
Tmin=poly_low.Tmin.value_si,
1668+
Tint=poly_high.Tmax.value_si,
1669+
Tmax=poly_low.Tmax.value_si,
16511670
elem_2=elem_2,
16521671
)
16531672

16541673
# Line 2
1655-
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E} 2\n'.format(thermo.polynomials[1].c0,
1656-
thermo.polynomials[1].c1,
1657-
thermo.polynomials[1].c2,
1658-
thermo.polynomials[1].c3,
1659-
thermo.polynomials[1].c4)
1674+
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E} 2\n'.format(poly_high.c0,
1675+
poly_high.c1,
1676+
poly_high.c2,
1677+
poly_high.c3,
1678+
poly_high.c4)
16601679

16611680
# Line 3
1662-
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E} 3\n'.format(thermo.polynomials[1].c5,
1663-
thermo.polynomials[1].c6,
1664-
thermo.polynomials[0].c0,
1665-
thermo.polynomials[0].c1,
1666-
thermo.polynomials[0].c2)
1681+
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E}{4:< 15.8E} 3\n'.format(poly_high.c5,
1682+
poly_high.c6,
1683+
poly_low.c0,
1684+
poly_low.c1,
1685+
poly_low.c2)
16671686

16681687
# Line 4
1669-
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E} 4\n'.format(thermo.polynomials[0].c3,
1670-
thermo.polynomials[0].c4,
1671-
thermo.polynomials[0].c5,
1672-
thermo.polynomials[0].c6)
1688+
string += '{0:< 15.8E}{1:< 15.8E}{2:< 15.8E}{3:< 15.8E} 4\n'.format(poly_low.c3,
1689+
poly_low.c4,
1690+
poly_low.c5,
1691+
poly_low.c6)
16731692

16741693
return string
16751694

@@ -2198,7 +2217,7 @@ def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_
21982217

21992218
# Thermodynamics section
22002219
f.write('THERM ALL\n')
2201-
f.write(' 300.000 1000.000 5000.000\n\n')
2220+
f.write(' 300.000 1000.000 5000.000\n\n')
22022221
for spec in sorted_species:
22032222
f.write(write_thermo_entry(spec, verbose=verbose))
22042223
f.write('\n')

‎rmgpy/molecule/element.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ def __init__(self, number, symbol, name, mass, isotope=-1, chemkin_name=None):
7777
self.name = name
7878
self.mass = mass
7979
self.isotope = isotope
80-
self.chemkin_name = chemkin_name or self.name
80+
self.chemkin_name = chemkin_name or self.symbol
8181
if symbol in {'X','L','R','e'}:
8282
self.cov_radius = 0
8383
else:

‎rmgpy/solver/files/collider_model/chem.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ END
2020

2121

2222
THERM ALL
23-
300.000 1000.000 5000.000
23+
300.000 1000.000 5000.000
2424

2525
Ar Ar1 G200.000 6000.000 1000.00 1
2626
2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ END
1515

1616

1717
THERM ALL
18-
300.000 1000.000 5000.000
18+
300.000 1000.000 5000.000
1919

2020
Ar Ar1 G200.000 6000.000 1000.00 1
2121
2.50000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2

‎rmgpy/tools/data/diffmodels/chem1.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ END
1919

2020

2121
THERM ALL
22-
300.000 1000.000 5000.000
22+
300.000 1000.000 5000.000
2323

2424
ethane(1) C 2 H 6 G100.000 5000.000 954.52 1
2525
4.58991157E+00 1.41506360E-02-4.75954111E-06 8.60275134E-10-6.21700682E-14 2

0 commit comments

Comments
 (0)