Skip to content

Commit df514c7

Browse files
authored
Merge pull request Auto-Mech#319 from snelliott/dev
Dev
2 parents 4d494c9 + 9a55bfc commit df514c7

File tree

5 files changed

+194
-18
lines changed

5 files changed

+194
-18
lines changed

mechanalyzer/builder/_update.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ def update_spc_dct(spc_infos, spc_dct, rename=False, enant_label=True,
4848
# when making name dictionary
4949
has_inf = False
5050
if spc_infos:
51-
if not isinstance(spc_ichs[0], str):
51+
if not isinstance(spc_infos[0], str):
5252
has_inf = True
5353

5454
_info_name_dct = ich_name_dct(spc_dct, incl_mult=has_inf, incl_chg=has_inf)

mechanalyzer/builder/rxn.py

+11-13
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,6 @@ def _print_species_info(rcts_lst):
3737
# Initialize new_spc (needed for sequential steps
3838
# new_spc_lst = None
3939
print('---------------------------------------------------------\n')
40-
4140
# Loop over the reaction series consisting of reactants and reaction type
4241
for sidx, series in enumerate(rxn_series):
4342
rct1_set, rct2_set, allowed_prds_lst, rxn_typs = series
@@ -59,11 +58,9 @@ def _print_species_info(rcts_lst):
5958
# Generate Reactions
6059
# for ichs in rct_ichs:
6160
# ini_rxns += generate_reactions(ichs, allowed_prd_ichs, rtyp)
62-
ini_rxns_infos_tmp = execute_function_in_parallel(
61+
ini_rxns_infos += execute_function_in_parallel(
6362
generate_reactions, rcts_lst, (mech_spc_dct, allowed_prds_lst, ReacInfo),
6463
nprocs=nprocs)
65-
print(ini_rxns_infos_tmp)
66-
print('does it make it here')
6764
# Add stereo
6865
if stereo:
6966
rxns = ()
@@ -89,11 +86,15 @@ def _print_species_info(rcts_lst):
8986
else:
9087
rxns_infos = ini_rxns_infos
9188

89+
print('rxn_infos')
90+
print(rxns_infos)
9291
# Update the mechanism objects with unique spc and rxns
9392
mech_spc_dct = update_spc_dct_from_reactions(
9493
rxns_infos, mech_spc_dct)
94+
print(list(mech_spc_dct.keys()))
9595
mech_rxn_dct = update_rxn_dct(
9696
rxns_infos, mech_rxn_dct, mech_spc_dct)
97+
print(list(mech_rxn_dct.keys()))
9798

9899
print('\n---------------------------------------------------------\n')
99100

@@ -104,8 +105,9 @@ def generate_reactions(
104105
mech_spc_dct, allowed_prds_lst, ReacInfo, rcts_lst, output_queue=None):
105106
""" For a given reactants
106107
"""
107-
108+
rxns_info = ()
108109
for rct_names in rcts_lst:
110+
print('rct names', rct_names)
109111
rct_ichs = [mech_spc_dct[name]['inchi'] for name in rct_names]
110112
# rct_mults = [mech_spc_dct[name]['mult'] for name in rct_names]
111113
_rct_smis = tuple(map(automol.chi.smiles, rct_ichs))
@@ -123,13 +125,12 @@ def generate_reactions(
123125
mech_spc_dct[name]['mult'],) for name in allowed_prds_lst)
124126
else:
125127
allowed_prds_info_lst = ()
126-
rcts_info = ((
128+
rcts_info = tuple((
127129
mech_spc_dct[name]['inchi'],
128130
mech_spc_dct[name]['charge'],
129131
mech_spc_dct[name]['mult'],) for name in rct_names)
130132

131133
# Build list of generated reactions including reactants and products
132-
rxn_infos = ()
133134
for pidx, prds_info in enumerate(prds_info_lst):
134135
# Don't add if requested products not found
135136
if allowed_prds_lst:
@@ -143,15 +144,12 @@ def generate_reactions(
143144
_prd_smis = tuple(map(automol.chi.smiles, (prd[0] for prd in prds_info)))
144145
log += f'\nFound Product(s) {pidx+1}: {_prd_smis}'
145146

146-
rxn_infos += ((rcts_info, prds_info, (None,)),)
147+
rxns_info += ((rcts_info, prds_info, (None,)),)
147148

148-
if not prds_info:
149+
if not prds_info_lst:
149150
log += '\nNO Product(s) Found'
150151
print(log)
151-
print('does it make it here')
152-
output_queue.put(rxn_infos)
153-
print('or here')
154-
# return rxn_ichs
152+
output_queue.put(rxns_info)
155153

156154

157155
# Functions to set lists for mech building step

mechanalyzer/cli/__init__.py

+92-2
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
import click
22

3-
from mechanalyzer.cli import run_sort, ste_mech
3+
from mechanalyzer.cli import run_sort, ste_mech, compare_rates
44

55

66
@click.group()
@@ -79,7 +79,97 @@ def sort(
7979
outgroups=outgroups,
8080
)
8181

82-
82+
@main.command()
83+
@click.option(
84+
"-m",
85+
"--mechs_yaml",
86+
default="mechs.yaml",
87+
show_default=True,
88+
help="yaml file, set up like:\\mech1:\\\trate_file: 'rate.ckin'\\\ttherm_file: 'therm.ckin'\\\tspecies_csv: 'species.csv'\\mech2: ...",
89+
)
90+
@click.option(
91+
"-o",
92+
"--plot_filename",
93+
default="rate_plot.pdf",
94+
show_default=True,
95+
help="name of output pdf of plots"
96+
)
97+
@click.option(
98+
"-f",
99+
"--out_txt_filename",
100+
default="ordering.txt",
101+
show_default=True,
102+
help="name of output text filename"
103+
)
104+
@click.option(
105+
"-d",
106+
"--job_path",
107+
default=".",
108+
show_default=True,
109+
help="directory for input/output files"
110+
)
111+
@click.option(
112+
"-t",
113+
"--temps_lst",
114+
default=None,
115+
type=lambda s: [float(temp) for temp in s.split(',')],
116+
show_default=True,
117+
help="array of temperatures, None defaults to 500--1500"
118+
)
119+
@click.option(
120+
"-p",
121+
"--pressures",
122+
default=None,
123+
type=lambda s: [float(press) for press in s.split(',')],
124+
show_default=True,
125+
help="array of pressures to plot, None defaults to (1, 10, 100)"
126+
)
127+
@click.option(
128+
"-s",
129+
"--sort_method",
130+
default="ratios",
131+
show_default=True,
132+
help="Sort the pdf plots by the ratio of the differences with 'ratios', or not at all with None"
133+
)
134+
@click.option(
135+
"-r",
136+
"--rev_rates",
137+
default=True,
138+
show_default=True,
139+
help="reverse rates of remaining mechanisms to make them match the direction in the first"
140+
)
141+
@click.option(
142+
"-l",
143+
"--remove_loners",
144+
default=True,
145+
show_default=True,
146+
help="True only plots rates that are in ALL mechanisms, False plots ALL rates in all mechanism"
147+
)
148+
def compare_mechanisms(
149+
mechs_yaml: str = 'mechs.yaml',
150+
plot_filename: str = 'rate_plot.pdf',
151+
out_txt_filename: str = 'ordering.txt',
152+
job_path: str = '.',
153+
temps_lst: list = [],
154+
pressures: list = [],
155+
sort_method: str = None,
156+
rev_rates: bool = True,
157+
remove_loners: bool = True,
158+
write_file: bool = False
159+
):
160+
"""Sort the reactions in a mechanism"""
161+
compare_rates.main(
162+
mechs_yaml=mechs_yaml,
163+
plot_filename=plot_filename,
164+
out_txt_filename=out_txt_filename,
165+
job_path=job_path,
166+
temps_lst=temps_lst,
167+
pressures=pressures,
168+
sort_method=sort_method,
169+
rev_rates=rev_rates,
170+
remove_loners=remove_loners,
171+
write_file=write_file,
172+
)
83173
@main.command()
84174
def expand():
85175
"""Expand stereochemistry for a mechanism"""

mechanalyzer/cli/compare_rates.py

+90
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
""" Script for running a comparison of rate constants between mechanisms
2+
"""
3+
4+
import os
5+
import sys
6+
import numpy
7+
import yaml
8+
import mechanalyzer.calculator.compare as compare
9+
import mechanalyzer.plotter.rates as plot_rates
10+
import mechanalyzer.plotter._util as util
11+
import mechanalyzer.parser.new_spc as spc_parser
12+
import mechanalyzer.parser.ckin_ as ckin_parser
13+
from ioformat import pathtools
14+
15+
16+
def read_mechs_yaml(mechs_yaml_file):
17+
""" read the yaml file that specifies each mechanism file
18+
for each mechanism, makes sure all files are given
19+
"""
20+
def _mech_is_fully_specified(mech_data):
21+
""" make sure all files are given in yaml for mechanism
22+
"""
23+
keys = ['rate_file', 'therm_file', 'species_csv']
24+
_fully_specified = True
25+
if not all(key in mech_data for key in keys):
26+
_fully_specified = False
27+
print(
28+
'missing in yaml: ',
29+
','.join([key for key in keys if key not in mech_data]))
30+
return _fully_specified
31+
32+
with open(mechs_yaml_file, 'r') as file:
33+
mechs = yaml.safe_load(file)
34+
mechs = {
35+
key: mech_files for key, mech_files in mechs.items()
36+
if _mech_is_fully_specified(mech_files)}
37+
38+
labels = []
39+
mech_files = []
40+
therm_files = []
41+
csv_files = []
42+
for key, mech in mechs.items():
43+
labels.append(key)
44+
mech_files.append(mech['rate_file'])
45+
therm_files.append(mech['therm_file'])
46+
csv_files.append(mech['species_csv'])
47+
48+
return labels, mech_files, therm_files, csv_files
49+
50+
51+
def main(
52+
mechs_yaml='mechs.yaml',
53+
plot_filename='rate_plot.pdf',
54+
out_txt_filename='ordering.txt',
55+
job_path='.',
56+
temps_lst=None, pressures=None,
57+
sort_method=None, rev_rates=True,
58+
remove_loners=True, write_file=False):
59+
60+
labels, mech_files, therm_files, csv_files = read_mechs_yaml(mechs_yaml)
61+
if temps_lst is None or len(temps_lst) < 2:
62+
temps_lst = [numpy.linspace(500, 1000, 16)]
63+
else:
64+
temps_lst = [numpy.linspace(temps_lst[0], temps_lst[1], 16)]
65+
if pressures is None or len(pressures) < 1:
66+
pressures = [1, 10, 100]
67+
68+
rxn_ktp_dcts = ckin_parser.load_rxn_ktp_dcts(
69+
mech_files, job_path, temps_lst, pressures)
70+
spc_therm_dcts = ckin_parser.load_spc_therm_dcts(
71+
therm_files, job_path, temps_lst[0]) # NOTE: taking first entry
72+
spc_dcts = spc_parser.load_mech_spc_dcts(csv_files, job_path)
73+
74+
# Get the algn_rxn_ktp_dct
75+
temps = temps_lst[0] # function receives a single Numpy array of temps
76+
algn_rxn_ktp_dct = compare.get_algn_rxn_ktp_dct(
77+
rxn_ktp_dcts, spc_therm_dcts, spc_dcts, temps, rev_rates=rev_rates,
78+
remove_loners=remove_loners, write_file=write_file)
79+
80+
# Run the plotter
81+
figs = plot_rates.build_plots(
82+
algn_rxn_ktp_dct,
83+
mech_names=labels,
84+
ratio_sort=bool(sort_method == 'ratios'))
85+
util.build_pdf(figs, filename=plot_filename, path=job_path)
86+
87+
# Write the ordered text file
88+
FSTR = compare.write_ordered_str(algn_rxn_ktp_dct, dct_type='rxn')
89+
pathtools.write_file(FSTR, job_path, out_txt_filename)
90+

mechanalyzer_bin/build_mech.py

-2
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,6 @@
2121
remove_comments='#', remove_whitespace=True)
2222

2323
FILE_DCT, RSERIES = mechanalyzer.parser.build_input_file(BLD_STR)
24-
2524
# Read the spc_dct and rxn_param_dct
2625
INP_SPC_STR = ioformat.pathtools.read_file(
2726
CWD, FILE_DCT['inp_spc'],
@@ -43,7 +42,6 @@
4342
rxn_param_dct = mechanalyzer.parser.mech.parse_mechanism(
4443
INP_MECH_STR, 'chemkin')
4544
isolate_spc, sort_lst, _ = mechanalyzer.parser.mech.parse_sort(SORT_STR)
46-
4745
# Generate the requested reactions
4846
mech_spc_dct, rxn_param_dct = mechanalyzer.builder.rxn.build_mechanism(
4947
mech_spc_dct, rxn_param_dct,

0 commit comments

Comments
 (0)