Skip to content

Commit 7d3edc1

Browse files
committed
polish console out, fn docstrings
1 parent 427aaf3 commit 7d3edc1

5 files changed

Lines changed: 111 additions & 83 deletions

File tree

navicat_mikimo/helper.py

Lines changed: 36 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99

1010

1111
def yesno(question):
12-
"""Simple Yes/No Function, Ruben's code"""
12+
"""Simple Yes/No Function, From Volcanic"""
1313
prompt = f"{question} ? (y/n): "
1414
ans = input(prompt).strip().lower()
1515
if ans not in ["y", "n"]:
@@ -21,7 +21,6 @@ def yesno(question):
2121

2222

2323
def check_existence(wdir, verb):
24-
2524
kinetic_mode = False
2625
k_exist = False
2726

@@ -35,10 +34,10 @@ def check_existence(wdir, verb):
3534
c0_exist = any([last_row_index.lower() in rows_to_search])
3635
k_exist = all([column in df.columns for column in columns_to_search])
3736
if not (c0_exist):
38-
logging.critical("Initial concentration not found in rxn_network.csv")
37+
logging.critical("Initial concentration not found in rxn_network.csv.")
3938

4039
else:
41-
logging.critical("rxn_network.csv not found")
40+
logging.critical("rxn_network.csv not found.")
4241

4342
filename = f"{wdir}reaction_data"
4443
extensions = [".csv", ".xls", ".xlsx"]
@@ -47,21 +46,21 @@ def check_existence(wdir, verb):
4746

4847
if energy_exist:
4948
if verb > 2:
50-
print("reaction_data file exists")
49+
print("Found reaction data file.")
5150
if k_exist:
52-
print("Both energy data and rate constants are provided")
51+
print("Both energy data and rate constants are provided.")
5352
else:
5453
if k_exist:
55-
print("reaction_data file not found, but rate constants are provided")
54+
print("reaction_data file not found, but rate constants are provided.")
5655
else:
5756
logging.critical(
58-
"reaction_data file not found and rate constants are not provided"
57+
"reaction_data file not found and rate constants are not provided."
5958
)
6059

6160
if os.path.exists(f"{wdir}kinetic_data.csv") or os.path.exists(
6261
f"{wdir}kinetic_data.xlsx"
6362
):
64-
kinetic_mode = yesno("kinetic_profile.csv exists, toggle to kinetic mode?")
63+
kinetic_mode = yesno("kinetic_profile.csv exists, toggle the kinetic mode?")
6564

6665
return kinetic_mode
6766

@@ -86,9 +85,9 @@ def check_km_inp(df, df_network, mode="energy"):
8685
if last_row_index.lower() in ["initial_conc", "c0", "initial conc"]:
8786
initial_conc = df_network.iloc[-1:].to_numpy()[0]
8887
df_network = df_network.drop(df_network.index[-1])
89-
logging.info("Initial conditions found")
88+
logging.info("Initial conditions found.")
9089
else:
91-
logging.critical("Initial conditions not found")
90+
logging.critical("Initial conditions not found.")
9291

9392
states_network = df_network.columns.to_numpy()[:]
9493
states_profile = df.columns.to_numpy()[1:]
@@ -116,19 +115,19 @@ def check_km_inp(df, df_network, mode="energy"):
116115
clear = False
117116
logging.warning(
118117
f"""\n{state} cannot be found in the reaction data, if it is in different name,
119-
change it to be the same in both reaction data and the network"""
118+
change it to be the same in both reaction data and the network."""
120119
)
121120
elif mode == "kinetic":
122121
if int((df.shape[1] - 1) / 2) != df_network.shape[0]:
123122
clear = False
124123
logging.critical(
125-
"Number of rate constants in the profile doesn't match with number of steps in the network input"
124+
"Number of rate constants in the profile doesn't match with number of steps in the network input."
126125
)
127126

128127
# initial conc
129128
if len(states_network) != len(initial_conc):
130129
clear = False
131-
logging.warning("\nYour initial conc seems wrong")
130+
logging.warning("\nYour initial conc seems wrong.")
132131

133132
# check network sanity
134133
mask = (~df_network.isin([-1, 1])).all(axis=1)
@@ -143,21 +142,20 @@ def check_km_inp(df, df_network, mode="energy"):
143142
if np.any(mask_R):
144143
clear = False
145144
logging.warning(
146-
f"\nThe reactant location: {states_network[r_indices[mask_R]]} appears wrong"
145+
f"\nThe reactant location: {states_network[r_indices[mask_R]]} appears wrong."
147146
)
148147

149148
mask_P = (~df_network.iloc[:, p_indices].isin([1])).all(axis=0).to_numpy()
150149
if np.any(mask_P):
151150
clear = False
152151
logging.warning(
153-
f"\nThe product location: {states_network[p_indices[mask_P]]} appears wrong"
152+
f"\nThe product location: {states_network[p_indices[mask_P]]} appears wrong."
154153
)
155154

156155
return clear
157156

158157

159158
def preprocess_data_mkm(arguments, mode):
160-
161159
parser = argparse.ArgumentParser(
162160
prog="mikimo",
163161
description="Perform microkinetic simulations and generate microkinetic volcano plots for homogeneous catalytic reactions.",
@@ -295,7 +293,7 @@ def preprocess_data_mkm(arguments, mode):
295293
type=int,
296294
default=1,
297295
help="""Run mode (default: 1)
298-
0: Run mkm for every profile in the input.
296+
0: Run mkm for every profile in the reaction data input.
299297
1: Construct microkinetic volcano plot.
300298
2: Construct microkinetic activity/selectivity map.
301299
""",
@@ -361,24 +359,24 @@ def preprocess_data_mkm(arguments, mode):
361359

362360
if t_finals is None:
363361
if verb > 0:
364-
print("No time input, use the default value of 54800 s (1d)")
362+
print("No time input, use the default value of 54800 s (1d).")
365363
t_finals = 54800
366364
elif len(t_finals) == 1:
367365
t_finals = t_finals[0]
368366
elif len(t_finals) > 1:
369367
if verb > 0:
370-
print("t_final is a range, use the first value")
368+
print("t_final is a range, use the first value.")
371369
t_finals = t_finals[0]
372370

373371
if temperatures is None:
374372
if verb > 0:
375-
print("No temperature input, use the default value of 298.15 K")
373+
print("No temperature input, use the default value of 298.15 K.")
376374
temperatures = 298.15
377375
elif len(temperatures) == 1:
378376
temperatures = temperatures[0]
379377
elif len(temperatures) > 1:
380378
if verb > 0:
381-
print("temperature is a range, use the first value")
379+
print("Temperature input is a range, use the first value.")
382380
temperatures = temperatures[0]
383381

384382
rnx = f"{wdir}/rxn_network.csv"
@@ -424,10 +422,10 @@ def preprocess_data_mkm(arguments, mode):
424422
clear = check_km_inp(df, df_network)
425423

426424
if not (clear):
427-
print("\nRecheck your reaction network and your reaction data\n")
425+
print("\nRecheck your reaction network and your reaction data.\n")
428426
else:
429427
if verb > 0:
430-
print("\nKM input is clear\n")
428+
print("\nMKM inputs are clean\n")
431429
return (
432430
dg,
433431
df_network,
@@ -486,16 +484,16 @@ def preprocess_data_mkm(arguments, mode):
486484
df = pd.read_csv(filename_csv)
487485
clear = check_km_inp(df, df_network)
488486
if not (clear):
489-
print("\nRecheck your reaction network and your reaction data\n")
487+
print("\nRecheck your reaction network and your reaction data.\n")
490488
else:
491489
if verb > 0:
492-
print("\nKM inputs are clear\n")
490+
print("\nMKM inputs are clean.\n")
493491

494492
tags = np.array([str(tag) for tag in df.columns[1:]], dtype=object)
495493
if ncore == -1:
496494
ncore = multiprocessing.cpu_count()
497495
if verb > 2:
498-
print(f"Use {ncore} cores for parallel computing")
496+
print(f"Use {ncore} cores for parallel computing.")
499497

500498
if plotmode == 0 and comp_ci:
501499
plotmode = 1
@@ -603,10 +601,10 @@ def preprocess_data_mkm(arguments, mode):
603601
clear = check_km_inp(df, df_network)
604602

605603
if not (clear):
606-
print("\nRecheck your reaction network and your reaction data\n")
604+
print("\nRecheck your reaction network and your reaction data.\n")
607605
else:
608606
if verb > 0:
609-
print("\nKM inputs are clear\n")
607+
print("\nMKM inputs are clean\n")
610608
return (
611609
dg,
612610
df_network,
@@ -629,7 +627,7 @@ def process_data_mkm(
629627
dg: np.ndarray, df_network: pd.DataFrame, tags: List[str], states: List[str]
630628
) -> Tuple[np.ndarray, List[np.ndarray], List[float], List[np.ndarray], np.ndarray]:
631629
"""
632-
Processes data for kinetic modeling.
630+
Processes data for micokinetic modeling.
633631
634632
Args:
635633
dg (numpy.array): free energy profile.
@@ -722,7 +720,7 @@ def process_data_mkm(
722720

723721

724722
def test_process_data_mkm():
725-
# Test data
723+
"""test data processor"""
726724
dg = np.array(
727725
[
728726
0.0,
@@ -828,9 +826,13 @@ def test_process_data_mkm():
828826
)
829827

830828
# Test the function
831-
initial_conc, energy_profile_all, dgr_all, coeff_TS_all, rxn_network_all = process_data_mkm(
832-
dg, df_network, tags, states
833-
)
829+
(
830+
initial_conc,
831+
energy_profile_all,
832+
dgr_all,
833+
coeff_TS_all,
834+
rxn_network_all,
835+
) = process_data_mkm(dg, df_network, tags, states)
834836

835837
# Compare the results
836838
assert np.array_equal(initial_conc, initial_conc_expected)

navicat_mikimo/kinetic_solver.py

Lines changed: 30 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -17,18 +17,18 @@
1717
warnings.filterwarnings("ignore")
1818

1919

20-
def erying(
20+
def eyring(
2121
dG_ddag: Union[float, np.ndarray], temperature: float
2222
) -> Union[float, np.ndarray]:
2323
"""
2424
Calculates the rate constant given the energy barrier and temperature based on Eyring equation.
2525
2626
Args:
27-
dG_ddag (float or array-like): Energy barrier or barriers.
27+
dG_ddag (float or array-like): Energy barrier(s).
2828
temperature (float): Temperature in Kelvin.
2929
3030
Returns:
31-
float or array-like: Erying rate constant or constants.
31+
float or array-like: Eyring rate constant(s).
3232
3333
"""
3434
R_: float = R * (1 / calorie) * (1 / kilo)
@@ -47,7 +47,7 @@ def get_k(
4747
4848
Parameters:
4949
energy_profile (array-like): The relative free energies profile (in kcal/mol).
50-
dgr (float): Free energy of the reaction.
50+
dgr (float): Free energy of the reaction (in kcal/mol).
5151
coeff_TS (one-hot array): One-hot encoding of the elements that are "TS" along the reaction coordinate.
5252
temperature (float): Temperature in Kelvin. Default is 298.15.
5353
@@ -57,7 +57,6 @@ def get_k(
5757
"""
5858

5959
def get_dG_ddag(energy_profile, dgr, coeff_TS):
60-
6160
# compute all dG_ddag in the profile
6261
n_S = energy_profile.size
6362
n_TS = np.count_nonzero(coeff_TS)
@@ -104,8 +103,8 @@ def get_dG_ddag(energy_profile, dgr, coeff_TS):
104103
energy_profile_reverse = np.insert(energy_profile_reverse, 0, 0)
105104
dG_ddag_reverse = get_dG_ddag(energy_profile_reverse, -dgr, coeff_TS_reverse)
106105

107-
k_forward = erying(dG_ddag_forward, temperature)
108-
k_reverse = erying(dG_ddag_reverse, temperature)
106+
k_forward = eyring(dG_ddag_forward, temperature)
107+
k_reverse = eyring(dG_ddag_reverse, temperature)
109108

110109
return k_forward, k_reverse[::-1]
111110

@@ -132,9 +131,9 @@ def calc_k(
132131
k_forward_all = []
133132
k_reverse_all = []
134133

135-
for i in range(len(energy_profile_all)):
134+
for energy_profile, dgr, coeff_TS in zip(energy_profile_all, dgr_all, coeff_TS_all):
136135
k_forward, k_reverse = get_k(
137-
energy_profile_all[i], dgr_all[i], coeff_TS_all[i], temperature=temperature
136+
energy_profile, dgr, coeff_TS, temperature=temperature
138137
)
139138
k_forward_all.extend(k_forward)
140139
k_reverse_all.extend(k_reverse)
@@ -266,7 +265,7 @@ def wrapper(t, y):
266265
# dy_dt[violate_up_idx] = dy_dt[violate_up_idx] + (boundary[violate_up_idx, 1] - y[violate_up_idx])/2
267266
dy_dt[violate_low_idx] = 0
268267
dy_dt[violate_up_idx] = 0
269-
except TypeError as e:
268+
except TypeError as err:
270269
y_ = np.array(y._value)
271270
dy_dt_ = np.array(dy_dt._value)
272271
# arraybox failure
@@ -312,7 +311,7 @@ def calc_km(
312311
ks=None,
313312
) -> Tuple[np.ndarray, Union[str, scipy.integrate._ivp.ivp.OdeResult]]:
314313
"""
315-
Calculate the kinetic model (KM) simulation.
314+
Perform MKM simulation.
316315
317316
Parameters:
318317
energy_profile_all (List): List of energy profiles for all steps.
@@ -527,7 +526,6 @@ def calc_km(
527526
upper = np.min(initial_conc[R_idx] * Rp_)
528527

529528
if report_as_yield:
530-
531529
c_target_yield = c_target_t / upper * 100
532530
c_target_yield[c_target_yield > 100] = 100
533531
c_target_yield[c_target_yield < 0] = 0
@@ -545,7 +543,6 @@ def calc_km(
545543

546544

547545
def test_get_k():
548-
549546
# Test case 1: CoL6 pincer co2
550547
energy_profile = np.array([0.0, 14.6, 0.5, 20.1, -1.7, 20.1])
551548
dgr = 2.2
@@ -740,10 +737,18 @@ def test_calc_dX_dt():
740737

741738

742739
def main():
743-
744-
dg, df_network, tags, states, t_final, temperature, x_scale, more_species_mkm, wdir, ks = preprocess_data_mkm(
745-
sys.argv[2:], mode="mkm_solo"
746-
)
740+
(
741+
dg,
742+
df_network,
743+
tags,
744+
states,
745+
t_final,
746+
temperature,
747+
x_scale,
748+
more_species_mkm,
749+
wdir,
750+
ks,
751+
) = preprocess_data_mkm(sys.argv[2:], mode="mkm_solo")
747752

748753
if ks is not None:
749754
t_span = (0, t_final)
@@ -769,9 +774,13 @@ def main():
769774
ks=ks,
770775
)
771776
else:
772-
initial_conc, energy_profile_all, dgr_all, coeff_TS_all, rxn_network_all = process_data_mkm(
773-
dg, df_network, tags, states
774-
)
777+
(
778+
initial_conc,
779+
energy_profile_all,
780+
dgr_all,
781+
coeff_TS_all,
782+
rxn_network_all,
783+
) = process_data_mkm(dg, df_network, tags, states)
775784
t_span = (0, t_final)
776785
_, result_solve_ivp = calc_km(
777786
energy_profile_all,
@@ -805,4 +814,4 @@ def main():
805814
for i in p_indices:
806815
print("--[{}]: {:.4f}--".format(states[i], result_solve_ivp.y[i][-1]))
807816

808-
print("\nWords that have faded to gray are colored like cappuccino\n")
817+
print("\nWords that have faded to gray are colored like cappuccino.\n")

0 commit comments

Comments
 (0)