diff --git a/INSTALL.rst b/INSTALL.rst index 95ee8023..d182179e 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -14,11 +14,13 @@ The package, together with the above dependencies, can be installed from pip install overreact -Optionally, extra functionality is provided by -`thermo `_:: +Optionally, extra functionality is provided such as a command-line interface +and solvent properties:: - pip install 'overreact[thermo]' + pip install 'overreact[cli,solvents]' -This last line installs thermo as well. thermo is used to calculate the -dynamic viscosity of solvents in the context of the +This last line installs `Rich `_ +and `thermo `_ as well. +Rich is used in the command-line interface, and thermo is used +to calculate the dynamic viscosity of solvents in the context of the :doc:`tutorials/collins-kimball` for diffusion-limited reactions. diff --git a/overreact/cli.py b/overreact/cli.py index cfccc607..03a52f0d 100644 --- a/overreact/cli.py +++ b/overreact/cli.py @@ -2,6 +2,7 @@ """Command-line interface.""" +import shutil import argparse import logging import os @@ -9,21 +10,22 @@ import sys import numpy as np +from rich.table import Table +from rich.table import Column +from rich.console import Console +from rich.markdown import Markdown +from rich import box +from overreact import __version__ from overreact import api from overreact import constants from overreact import core from overreact import io -levels = [logging.WARNING, logging.INFO, logging.DEBUG] - -def summarize_model( - model, quantities=None, savepath=None, plot=False, qrrho=True, temperature=298.15 -): - """Produce a string describing a model. - - The returned string contains a final line break. +# TODO(schneiderfelipe): test this class +class Report: + """Produce a report object based on a model. Parameters ---------- @@ -36,232 +38,253 @@ def summarize_model( temperature : array-like, optional Absolute temperature in Kelvin. - Returns - ------- - str - Examples -------- >>> model = api.parse_model("data/ethane/B97-3c/model.jk") - >>> print(summarize_model(model)) - - (READ) REACTIONS - ------------------------------------------------------------------------------ - S -> E‡ -> S - - (PARSED) REACTIONS - -------------------------------------------------------------------------------------------------------------------------------------- - no reactants via‡ products half equilib.? - -- --------------------------------------------------- ------------ --------------------------------------------------- -------------- - 0 S E‡ S - - COMPOUNDS - ------------------------------------------------------------------------------------------------------------------------- - no compound elec. energy spin mult. smallest vibfreqs original logfile - [Eh] [cm⁻¹] - -- ----------------- ----------------- ---------- ------------------------- --------------------------------------------- - 0 S -79.788170457691 1 307.6, 825.4, 826.1 data/ethane/B97-3c/staggered.out - 1 E‡ -79.783894160233 1 -298.9, 902.2, 902.5 data/ethane/B97-3c/eclipsed.out - - CALCULATED THERMOCHEMISTRY (COMPOUNDS) - ----------------------------------------------------------------------------------- - no compound mass Gcorr(298.15K) Ucorr(298.15K) Hcorr(298.15K) S(298.15K) - [amu] [kcal/mol] [kcal/mol] [kcal/mol] [cal/mol·K] - -- ----------------- ------ -------------- -------------- -------------- ---------- - 0 S 30.07 33.01 48.63 49.22 54.40 - 1 E‡ 30.07 32.95 48.15 48.74 52.96 - - CALCULATED THERMOCHEMISTRY (REACTIONS) - ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - no reaction Δmass‡ ΔG‡ ΔE‡ ΔU‡ ΔH‡ ΔS‡ Δmass° ΔG° ΔE° ΔU° ΔH° ΔS° - [amu] [kcal/mol] [kcal/mol] [kcal/mol] [kcal/mol] [cal/mol·K] [amu] [kcal/mol] [kcal/mol] [kcal/mol] [kcal/mol] [cal/mol·K] - -- ------------------------------------------------------------------- ------ ---------- ---------- ---------- ---------- ----------- ------ ---------- ---------- ---------- ---------- ----------- - 0 S -> S 0.00 2.63 2.68 2.20 2.20 -1.44 0.00 0.00 0.00 0.00 0.00 0.00 - - CALCULATED KINETICS - --------------------------------------------------------------------------------------------------------------------------------------- - no reaction half equilib.? k k k - [M⁻ⁿ⁺¹·s⁻¹] [(cm³/particle)ⁿ⁻¹·s⁻¹] [atm⁻ⁿ⁺¹·s⁻¹] - -- ------------------------------------------------------------------- -------------- ----------- ----------------------- ------------- - 0 S -> S 8.2e+10 8.2e+10 8.2e+10 + >>> Report(model) # doctest: +SKIP + ╔══════════════════════════════════════════════════════════════════════════════╗ + ║ overreact 1.0 ║ + ╚══════════════════════════════════════════════════════════════════════════════╝ + + Construct precise chemical microkinetic models from first principles + + + Input description + + (read) reactions + + S -> E‡ -> S + + (parsed) reactions + ╷ ╷ ╷ ╷ + no │ reactant │ via‡ │ product │ half equilib.? + ╶────┼──────────┼──────┼─────────┼────────────────╴ + 0 │ S │ E‡ │ S │ + ╵ ╵ ╵ ╵ + logfiles + ╷ ╷ + no │ compound │ path + ╶────┼──────────┼──────────────────────────────────╴ + 0 │ S │ data/ethane/B97-3c/staggered.out + 1 │ E‡ │ data/ethane/B97-3c/eclipsed.out + ╵ ╵ + compounds + ╷ ╷ ╷ ╷ + no │ compound │ elec. energy │ spin mult. │ smallest vibfreqs + │ │ [Eₕ] │ │ [cm⁻¹] + ╶────┼──────────┼───────────────────┼────────────┼────────────────────────╴ + 0 │ S │ -79.788170457691 │ 1 │ +307.6, +825.4, +826.1 + 1 │ E‡ │ -79.783894160233 │ 1 │ -298.9, +902.2, +902.5 + ╵ ╵ ╵ ╵ + Temperature = 298.15 K + + Output section + + calculated thermochemistry (compounds) + ╷ ╷ ╷ ╷ ╷ ╷ + no │ compound │ mass │ Gᶜᵒʳʳ │ Uᶜᵒʳʳ │ Hᶜᵒʳʳ │ S + │ │ [amu] │ [kcal/mol] │ [kcal/mol] │ [kcal/mol] │ [cal/mol·… + ╶────┼──────────┼────────┼─────────────┼────────────┼─────────────┼────────────╴ + 0 │ S │ 30.07 │ 3… │ … │ 4… │ 54.40 + 1 │ E‡ │ 30.07 │ 3… │ … │ 4… │ 52.96 + ╵ ╵ ╵ ╵ ╵ ╵ + calculated thermochemistry (reactions°) + ╷ ╷ ╷ ╷ ╷ ╷ ╷ + no │ reaction │ Δmass° │ ΔG° │ ΔE° │ ΔU° │ ΔH° │ ΔS° + │ │ [amu] │ [kcal/m… │ [kcal/m… │ [kcal/m… │ [kcal/m… │ [cal/m… + ╶────┼──────────┼────────┼──────────┼──────────┼──────────┼──────────┼─────────╴ + 0 │ S -> S │ 0.00 │ 0… │ 0… │ 0… │ 0… │ … + ╵ ╵ ╵ ╵ ╵ ╵ ╵ + calculated thermochemistry (reactions‡) + ╷ ╷ ╷ ╷ ╷ ╷ ╷ + no │ reaction │ Δmass‡ │ ΔG‡ │ ΔE‡ │ ΔU‡ │ ΔH‡ │ ΔS‡ + │ │ [amu] │ [kcal/m… │ [kcal/m… │ [kcal/m… │ [kcal/m… │ [cal/m… + ╶────┼──────────┼────────┼──────────┼──────────┼──────────┼──────────┼─────────╴ + 0 │ S -> S │ 0.00 │ 2… │ 2… │ 2… │ 2… │ … + ╵ ╵ ╵ ╵ ╵ ╵ ╵ + calculated kinetics + ╷ ╷ ╷ ╷ ╷ + no │ reaction │ half equilib.? │ k │ k │ k + │ │ │ [M⁻ⁿ⁺¹·s⁻¹] │ [(cm³/partic… │ [atm⁻ⁿ⁺¹·s⁻¹] + ╶────┼──────────┼────────────────┼─────────────┼───────────────┼───────────────╴ + 0 │ S -> S │ │ 8.2e+10 │ 8.2e+10 │ 8.2e+10 + ╵ ╵ ╵ ╵ ╵ """ - sections = [] - sections.append(_summarize_scheme(model.scheme)) - sections.append(_summarize_compounds(model.compounds)) - sections.append( - _summarize_thermochemistry( - model.scheme, model.compounds, qrrho=qrrho, temperature=temperature - ) - ) - sections.append( - _summarize_kinetics( - model.scheme, - model.compounds, - quantities=quantities, - savepath=savepath, - plot=plot, - qrrho=qrrho, - temperature=temperature, - ) - ) - return "".join(sections) - - -def _summarize_scheme(scheme): - """Produce a string describing a reaction scheme. - This is meant to be used from within `summarize_model`. The returned string - contains a final line break. - - Parameters - ---------- - scheme : Scheme - - Returns - ------- - str - """ - scheme = core._check_scheme(scheme) - reactions = _format_table( - [[row] for row in core.unparse_reactions(scheme).split("\n")], - title="(read) reactions", - length=[78], - ) - - transition_states = core.get_transition_states( - scheme.A, scheme.B, scheme.is_half_equilibrium - ) - - reaction_rows = [ - [ - "no", - "reactants".center(51), - "via‡".center(12), - "products".center(51), - "half equilib.?", - ] - ] - reaction_rows.append(["-" * len(field) for field in reaction_rows[0]]) - for i, reaction in enumerate(scheme.reactions): - reactants, _, products = re.split(r"\s*(->|<=>|<-)\s*", reaction) - row = [f"{i:2d}", reactants, None, products, None] - if transition_states[i] is not None: - row[2] = scheme.compounds[transition_states[i]] - elif scheme.is_half_equilibrium[i]: - row[4] = True - reaction_rows.append(row) - parsed_reactions = _format_table(reaction_rows, title="(parsed) reactions") - - return reactions + parsed_reactions - - -def _summarize_compounds(compounds): - """Produce a string describing compounds. - - This is meant to be used from within `summarize_model`. The returned string - contains a final line break. - - Parameters - ---------- - compounds : dict-like + def __init__( + self, + model, + quantities=None, + savepath=None, + plot=False, # TODO(schneiderfelipe): change to do_plot + qrrho=True, # TODO(schneiderfelipe): change to use_qrrho + temperature=298.15, + box_style=box.SIMPLE, + ): + self.model = model + self.quantities = quantities + self.savepath = savepath + self.plot = plot + self.qrrho = qrrho + self.temperature = temperature + self.box_style = box_style + + def __rich_console__(self, console, options): + """ + Implement Rich Console protocol. + + This works by yielding from generators. + + Yields + ------ + renderable + """ + yield from self._yield_scheme() + yield from self._yield_compounds() + yield from self._yield_thermochemistry() + yield from self._yield_kinetics() + + def _yield_scheme(self): + """Produce a renderables describing the reaction scheme. + + This is meant to be used from within `__rich_console__`. + + Yields + ------ + renderable + """ + scheme = core._check_scheme(self.model.scheme) + + raw_table = Table( + title="(read) reactions", box=self.box_style, show_header=False + ) + raw_table.add_column(justify="center") + for r in core.unparse_reactions(scheme).split("\n"): + raw_table.add_row(r) + yield raw_table - Returns - ------- - str + transition_states = core.get_transition_states( + scheme.A, scheme.B, scheme.is_half_equilibrium + ) - Raises - ------ - ValueError - If at least one compound has no data defined. - """ - undefined_compounds = [] - for name in compounds: - if not compounds[name]: - undefined_compounds.append(name) - if undefined_compounds: - raise ValueError(f"undefined compounds: {', '.join(undefined_compounds)}") - - compound_rows = [ - [ - "no", - "compound".center(17), - "elec. energy".center(17), - "spin mult.", - "smallest vibfreqs".center(25), - "original logfile".center(45), - ], - [None, None, "[Eh]", None, "[cm⁻¹]", None], - ] - compound_rows.append(["-" * len(field) for field in compound_rows[0]]) - for i, (name, data) in enumerate(compounds.items()): - compound_rows.append( - [ + parsed_table = Table( + Column("no", justify="center"), + Column("reactant", justify="center"), + Column("via‡", justify="center"), + Column("product", justify="center"), + Column("half equilib.?", justify="center"), + title="(parsed) reactions", + box=self.box_style, + ) + for i, reaction in enumerate(scheme.reactions): + reactants, _, products = re.split(r"\s*(->|<=>|<-)\s*", reaction) + # TODO(schneiderfelipe): should we use "No" instead of None for + # "half-equilib.?"? + row = [f"{i:2d}", reactants, None, products, None] + if transition_states[i] is not None: + row[2] = scheme.compounds[transition_states[i]] + elif scheme.is_half_equilibrium[i]: + row[4] = True + parsed_table.add_row(*row) + yield parsed_table + + def _yield_compounds(self): + """Produce a renderables describing the compounds. + + This is meant to be used from within `__rich_console__`. + + Yields + ------ + renderable + + Raises + ------ + ValueError + If at least one compound has no data defined. + """ + undefined_compounds = [] + for name in self.model.compounds: + if not self.model.compounds[name]: + undefined_compounds.append(name) + if undefined_compounds: + raise ValueError(f"undefined compounds: {', '.join(undefined_compounds)}") + + logfiles_table = Table( + Column("no", justify="center"), + Column("compound", justify="center"), + Column("path", justify="center"), + title="logfiles", + box=self.box_style, + ) + compounds_table = Table( + Column("no", justify="center"), + Column("compound", justify="center"), + Column("elec. energy\n\[Eₕ]", justify="center"), + Column("spin mult.", justify="center"), + Column("smallest vibfreqs\n\[cm⁻¹]", justify="center"), + title="compounds", + box=self.box_style, + ) + for i, (name, data) in enumerate(self.model.compounds.items()): + logfiles_table.add_row( f"{i:2d}", name, - f"{data.energy / (constants.hartree * constants.N_A):17.12f}", - data.mult, - ", ".join([f"{vibfreq:7.1f}" for vibfreq in data.vibfreqs[:3]]), # TODO(schneiderfelipe): show only the file name and inform # the absolute path to folder (as a bash variable) somewhere # else. data.logfile, - ] - ) - - return _format_table(compound_rows, title="compounds") - - -def _summarize_thermochemistry(scheme, compounds, qrrho=True, temperature=298.15): - """Produce a string describing the thermochemistry of a reaction scheme. + ) + compounds_table.add_row( + f"{i:2d}", + name, + f"{data.energy / (constants.hartree * constants.N_A):17.12f}", + f"{data.mult}", + ", ".join([f"{vibfreq:+.1f}" for vibfreq in data.vibfreqs[:3]]), + ) + yield logfiles_table + yield compounds_table - This is meant to be used from within `summarize_model`. The returned string - contains a final line break. + def _yield_thermochemistry(self): + """Produce a renderables describing the thermochemistry of the reaction scheme. - Parameters - ---------- - scheme : Scheme - compounds : dict-like - qrrho : bool, optional - Apply both the quasi-rigid rotor harmonic oscilator (QRRHO) - approximations of M. Head-Gordon (enthalpy correction, see - doi:10.1021/jp509921r) and S. Grimme (entropy correction, see - doi:10.1002/chem.201200497) on top of the classical RRHO. - temperature : array-like, optional - Absolute temperature in Kelvin. + This is meant to be used from within `__rich_console__`. - Returns - ------- - str - """ - scheme = core._check_scheme(scheme) + Yields + ------ + renderable + """ + scheme = core._check_scheme(self.model.scheme) - molecular_masses = np.array( - [np.sum(data.atommasses) for name, data in compounds.items()] - ) - energies = np.array([data.energy for name, data in compounds.items()]) - internal_energies = api.get_internal_energies( - compounds, qrrho=qrrho, temperature=temperature - ) - enthalpies = api.get_enthalpies(compounds, qrrho=qrrho, temperature=temperature) - entropies = api.get_entropies(compounds, qrrho=qrrho, temperature=temperature) - freeenergies = enthalpies - temperature * entropies - - compound_rows = [ - [ - "no", - "compound".center(17), - "mass".center(6), - f"Gcorr({temperature}K)", - f"Ucorr({temperature}K)", - f"Hcorr({temperature}K)", - f"S({temperature}K)", - ], - [None, None, "[amu]", "[kcal/mol]", "[kcal/mol]", "[kcal/mol]", "[cal/mol·K]"], - ] - compound_rows.append(["-" * len(field) for field in compound_rows[0]]) - for i, (name, data) in enumerate(compounds.items()): - compound_rows.append( - [ + molecular_masses = np.array( + [np.sum(data.atommasses) for name, data in self.model.compounds.items()] + ) + energies = np.array( + [data.energy for name, data in self.model.compounds.items()] + ) + internal_energies = api.get_internal_energies( + self.model.compounds, qrrho=self.qrrho, temperature=self.temperature + ) + enthalpies = api.get_enthalpies( + self.model.compounds, qrrho=self.qrrho, temperature=self.temperature + ) + entropies = api.get_entropies( + self.model.compounds, qrrho=self.qrrho, temperature=self.temperature + ) + freeenergies = enthalpies - self.temperature * entropies + + compounds_table = Table( + Column("no", justify="center"), + Column("compound", justify="center"), + Column("mass\n\[amu]", justify="center"), + Column("Gᶜᵒʳʳ\n\[kcal/mol]", justify="center", style="green"), + Column("Uᶜᵒʳʳ\n\[kcal/mol]", justify="center"), + Column("Hᶜᵒʳʳ\n\[kcal/mol]", justify="center"), + Column("S\n\[cal/mol·K]", justify="center"), + title="calculated thermochemistry (compounds)", + box=self.box_style, + ) + for i, (name, data) in enumerate(self.model.compounds.items()): + compounds_table.add_row( f"{i:2d}", name, f"{molecular_masses[i]:6.2f}", @@ -269,248 +292,210 @@ def _summarize_thermochemistry(scheme, compounds, qrrho=True, temperature=298.15 f"{(internal_energies[i] - data.energy) / constants.kcal:14.2f}", f"{(enthalpies[i] - data.energy) / constants.kcal:14.2f}", f"{entropies[i] / constants.calorie:10.2f}", - ] + ) + yield compounds_table + + delta_mass = api.get_delta(scheme.A, molecular_masses) + delta_energies = api.get_delta(scheme.A, energies) + delta_internal_energies = api.get_delta(scheme.A, internal_energies) + delta_enthalpies = api.get_delta(scheme.A, enthalpies) + delta_entropies = api.get_delta(scheme.A, entropies) + delta_freeenergies = api.get_delta(scheme.A, freeenergies) + + delta_activation_mass = api.get_delta(scheme.B, molecular_masses) + delta_activation_energies = api.get_delta(scheme.B, energies) + delta_activation_internal_energies = api.get_delta(scheme.B, internal_energies) + delta_activation_enthalpies = api.get_delta(scheme.B, enthalpies) + delta_activation_entropies = api.get_delta(scheme.B, entropies) + delta_activation_freeenergies = api.get_delta(scheme.B, freeenergies) + + circ_table = Table( + Column("no", justify="center"), + Column("reaction", justify="center"), + Column("Δmass°\n\[amu]", justify="center"), + Column("ΔG°\n\[kcal/mol]", justify="center", style="green"), + Column("ΔE°\n\[kcal/mol]", justify="center"), + Column("ΔU°\n\[kcal/mol]", justify="center"), + Column("ΔH°\n\[kcal/mol]", justify="center"), + Column("ΔS°\n\[cal/mol·K]", justify="center"), + title="calculated thermochemistry (reactions°)", + box=self.box_style, ) - compounds_table = _format_table( - compound_rows, title="calculated thermochemistry (compounds)" - ) - - delta_mass = api.get_delta(scheme.A, molecular_masses) - delta_energies = api.get_delta(scheme.A, energies) - delta_internal_energies = api.get_delta(scheme.A, internal_energies) - delta_enthalpies = api.get_delta(scheme.A, enthalpies) - delta_entropies = api.get_delta(scheme.A, entropies) - delta_freeenergies = api.get_delta(scheme.A, freeenergies) - - delta_activation_mass = api.get_delta(scheme.B, molecular_masses) - delta_activation_energies = api.get_delta(scheme.B, energies) - delta_activation_internal_energies = api.get_delta(scheme.B, internal_energies) - delta_activation_enthalpies = api.get_delta(scheme.B, enthalpies) - delta_activation_entropies = api.get_delta(scheme.B, entropies) - delta_activation_freeenergies = api.get_delta(scheme.B, freeenergies) - - reaction_rows = [ - [ - "no", - "reaction".center(67), - "Δmass‡", - "ΔG‡".center(10), - "ΔE‡".center(10), - "ΔU‡".center(10), - "ΔH‡".center(10), - "ΔS‡".center(11), - "Δmass°", - "ΔG°".center(10), - "ΔE°".center(10), - "ΔU°".center(10), - "ΔH°".center(10), - "ΔS°".center(11), - ], - [ - None, - None, - "[amu]", - "[kcal/mol]", - "[kcal/mol]", - "[kcal/mol]", - "[kcal/mol]", - "[cal/mol·K]", - "[amu]", - "[kcal/mol]", - "[kcal/mol]", - "[kcal/mol]", - "[kcal/mol]", - "[cal/mol·K]", - ], - ] - reaction_rows.append(["-" * len(field) for field in reaction_rows[0]]) - for i, reaction in enumerate(scheme.reactions): - if scheme.is_half_equilibrium[i]: - row = [ - f"{i:2d}", - reaction, - None, - None, - None, - None, - None, - None, - f"{delta_mass[i]:6.2f}", - f"{delta_freeenergies[i] / constants.kcal:10.2f}", - f"{delta_energies[i] / constants.kcal:10.2f}", - f"{delta_internal_energies[i] / constants.kcal:10.2f}", - f"{delta_enthalpies[i] / constants.kcal:10.2f}", - f"{delta_entropies[i] / constants.calorie:11.2f}", - ] - else: - row = [ - f"{i:2d}", - reaction, - f"{delta_activation_mass[i]:6.2f}", - f"{delta_activation_freeenergies[i] / constants.kcal:10.2f}", - f"{delta_activation_energies[i] / constants.kcal:10.2f}", - f"{delta_activation_internal_energies[i] / constants.kcal:10.2f}", - f"{delta_activation_enthalpies[i] / constants.kcal:10.2f}", - f"{delta_activation_entropies[i] / constants.calorie:11.2f}", - f"{delta_mass[i]:6.2f}", - f"{delta_freeenergies[i] / constants.kcal:10.2f}", - f"{delta_energies[i] / constants.kcal:10.2f}", - f"{delta_internal_energies[i] / constants.kcal:10.2f}", - f"{delta_enthalpies[i] / constants.kcal:10.2f}", - f"{delta_entropies[i] / constants.calorie:11.2f}", - ] - - reaction_rows.append(row) - reactions_table = _format_table( - reaction_rows, title="calculated thermochemistry (REACTIONS)" - ) - - return compounds_table + reactions_table - - -def _summarize_kinetics( - scheme, - compounds, - quantities=None, - savepath=None, - plot=False, - qrrho=True, - temperature=298.15, -): - # TODO(schneiderfelipe): apply other corrections to k - k = { - "M⁻ⁿ⁺¹·s⁻¹": api.get_k( - scheme, compounds, qrrho=qrrho, temperature=temperature, scale="l mol-1 s-1" - ), - "(cm³/particle)ⁿ⁻¹·s⁻¹": api.get_k( - scheme, - compounds, - qrrho=qrrho, - temperature=temperature, - scale="cm3 particle-1 s-1", - ), - "atm⁻ⁿ⁺¹·s⁻¹": api.get_k( - scheme, compounds, qrrho=qrrho, temperature=temperature, scale="atm-1 s-1" - ), - } - - reaction_rows = [ - ["no", "reaction".center(67), "half equilib.?"] - + ["k".center(len(scale) + 2) for scale in k], - [None, None, None] + [f"[{scale}]" for scale in k], - ] - reaction_rows.append(["-" * len(field) for field in reaction_rows[0]]) - for i, reaction in enumerate(scheme.reactions): - row = [f"{i:2d}", reaction, None] + [f"{k[scale][i]:7.2g}" for scale in k] - if scheme.is_half_equilibrium[i]: - row[2] = True - - reaction_rows.append(row) - reactions_table = _format_table(reaction_rows, title="calculated kinetics") - - if quantities is not None and quantities: - # TODO(schneiderfelipe): apply post-processing to scheme, k (with functions - # that receive a scheme, k and return a scheme, k). One that solves the pH - # problem is welcome: get a scheme, k and, for each reaction in it, remove - # the H+ and multiplies the reaction rate constants by the proper - # concentration if there is H+ in the reactants. - # TODO(schneiderfelipe): encapsulate everything in a function that depends - # on the freeenergies as first parameter - dydt = api.get_dydt(scheme, k) - - y0 = np.zeros(len(scheme.compounds)) - for spec in quantities: - fields = spec.split(":", 1) - name = fields[0] - try: - quantity = float(fields[1]) - except (IndexError, ValueError): - raise ValueError( - f"badly formatted quantities: '{' '.join(quantities)}'" + dagger_table = Table( + Column("no", justify="center"), + Column("reaction", justify="center"), + Column("Δmass‡\n\[amu]", justify="center"), + Column("ΔG‡\n\[kcal/mol]", justify="center", style="green"), + Column("ΔE‡\n\[kcal/mol]", justify="center"), + Column("ΔU‡\n\[kcal/mol]", justify="center"), + Column("ΔH‡\n\[kcal/mol]", justify="center"), + Column("ΔS‡\n\[cal/mol·K]", justify="center"), + title="calculated thermochemistry (reactions‡)", + box=self.box_style, + ) + for i, reaction in enumerate(scheme.reactions): + if scheme.is_half_equilibrium[i]: + circ_row = [ + f"{i:2d}", + reaction, + f"{delta_mass[i]:6.2f}", + f"{delta_freeenergies[i] / constants.kcal:10.2f}", + f"{delta_energies[i] / constants.kcal:10.2f}", + f"{delta_internal_energies[i] / constants.kcal:10.2f}", + f"{delta_enthalpies[i] / constants.kcal:10.2f}", + f"{delta_entropies[i] / constants.calorie:11.2f}", + ] + dagger_row = [ + f"{i:2d}", + reaction, + None, + None, + None, + None, + None, + None, + ] + else: + circ_row = [ + f"{i:2d}", + reaction, + f"{delta_mass[i]:6.2f}", + f"{delta_freeenergies[i] / constants.kcal:10.2f}", + f"{delta_energies[i] / constants.kcal:10.2f}", + f"{delta_internal_energies[i] / constants.kcal:10.2f}", + f"{delta_enthalpies[i] / constants.kcal:10.2f}", + f"{delta_entropies[i] / constants.calorie:11.2f}", + ] + dagger_row = [ + f"{i:2d}", + reaction, + f"{delta_activation_mass[i]:6.2f}", + f"{delta_activation_freeenergies[i] / constants.kcal:10.2f}", + f"{delta_activation_energies[i] / constants.kcal:10.2f}", + f"{delta_activation_internal_energies[i] / constants.kcal:10.2f}", + f"{delta_activation_enthalpies[i] / constants.kcal:10.2f}", + f"{delta_activation_entropies[i] / constants.calorie:11.2f}", + ] + + circ_table.add_row(*circ_row) + dagger_table.add_row(*dagger_row) + yield circ_table + yield dagger_table + + def _yield_kinetics(self): + """Produce a renderables describing the kinetics of the system. + + This is meant to be used from within `__rich_console__`. + + Yields + ------ + renderable + """ + # TODO(schneiderfelipe): apply other corrections to k (such as + # diffusion control) + k = { + "M⁻ⁿ⁺¹·s⁻¹": api.get_k( + self.model.scheme, + self.model.compounds, + qrrho=self.qrrho, + temperature=self.temperature, + scale="l mol-1 s-1", + ), + "(cm³/particle)ⁿ⁻¹·s⁻¹": api.get_k( + self.model.scheme, + self.model.compounds, + qrrho=self.qrrho, + temperature=self.temperature, + scale="cm3 particle-1 s-1", + ), + "atm⁻ⁿ⁺¹·s⁻¹": api.get_k( + self.model.scheme, + self.model.compounds, + qrrho=self.qrrho, + temperature=self.temperature, + scale="atm-1 s-1", + ), + } + + kinetics_table = Table( + *( + [ + Column("no", justify="center"), + Column("reaction", justify="center"), + Column("half equilib.?", justify="center"), + ] + + [Column(f"k\n\[{scale}]", justify="center") for scale in k] + ), + title="calculated kinetics", + box=self.box_style, + ) + for i, reaction in enumerate(self.model.scheme.reactions): + row = [f"{i:2d}", reaction, None] + [f"{k[scale][i]:7.2g}" for scale in k] + if self.model.scheme.is_half_equilibrium[i]: + row[2] = True + + kinetics_table.add_row(*row) + yield kinetics_table + + if self.quantities is not None and self.quantities: + scale = "M⁻ⁿ⁺¹·s⁻¹" + + # TODO(schneiderfelipe): apply post-processing to scheme, k (with functions + # that receive a scheme, k and return a scheme, k). One that solves the pH + # problem is welcome: get a scheme, k and, for each reaction in it, remove + # the H+ and multiplies the reaction rate constants by the proper + # concentration if there is H+ in the reactants. + # TODO(schneiderfelipe): encapsulate everything in a function that depends + # on the freeenergies as first parameter + dydt = api.get_dydt(self.model.scheme, k[scale]) + + y0 = np.zeros(len(self.model.scheme.compounds)) + for spec in self.quantities: + fields = spec.split(":", 1) + name = fields[0] + try: + quantity = float(fields[1]) + except (IndexError, ValueError): + raise ValueError( + f"badly formatted quantities: '{' '.join(self.quantities)}'" + ) + + # TODO(schneiderfelipe): the following is inefficient but probably OK + y0[self.model.scheme.compounds.index(name)] = quantity + + # TODO(schneiderfelipe): be clever and automatically select a + # suitable max time for simulation + t_span = [0, 1e5 / np.min(k[scale])] + + # TODO(schneiderfelipe): we can use rates as well (third returned value) + t, y, _ = api.get_y(dydt, y0=y0, t_span=t_span, method="Radau") + if self.savepath is not None: + np.savetxt( + self.savepath, + np.block([t[:, np.newaxis], y.T]), + header=f"t,{','.join(self.model.scheme.compounds)}", + delimiter=",", ) + yield f"CSV file saved to {self.savepath}" - # TODO(schneiderfelipe): the following is inefficient but probably OK - y0[scheme.compounds.index(name)] = quantity - - t, y, r = api.get_y(dydt, y0=y0, method="Radau") - if plot: - import matplotlib.pyplot as plt - - for i, name in enumerate(scheme.compounds): - if not core.is_transition_state(name): - plt.plot(t, y[i], label=name) + if self.plot: + import matplotlib.pyplot as plt - plt.legend() - plt.xlabel("Time (s)") - plt.ylabel("Concentration (M)") - plt.show() - - if savepath is not None: - np.savetxt( - savepath, - np.block([t[:, np.newaxis], y.T]), - header=f"t,{','.join(scheme.compounds)}", - delimiter=",", - ) - - # TODO(schneiderfelipe): implement the degree of rate control - - return reactions_table - - -def _create_banner(text, title_char="-", width=78): - banner = f"\n\n{text.center(width)}" - return banner + f"\n{title_char * width}\n" + for i, name in enumerate(self.model.scheme.compounds): + if not core.is_transition_state(name): + plt.plot(t, y[i], label=name) - -def _format_table(rows, length=None, sep=" ", title=None): - """Format a table for printing. - - The width of each column is taken from the length of the first row. Nones - are omitted. - - Parameters - ---------- - rows : sequence of sequence of str - length : sequence of int, optional - sep : str, optional - - Returns - ------- - str - - Examples - -------- - >>> print(_format_table([["one ", "two"], [1, 2], ["hello", "world"]])) - one two - 1 2 - hello world - """ - default_length = [len(field) for field in rows[0]] - if length is None: - length = default_length - else: - length = [a if a is not None else b for a, b in zip(length, default_length)] - - lines = [] - for row in rows: - line = [ - str(field).center(length[i]) if field is not None else " " * length[i] - for i, field in enumerate(row) - ] - lines.append(sep.join(line)) - - table = "\n".join(lines) - if title is not None: - return ( - _create_banner(title.upper(), width=np.sum(length) + len(length) - 1) - + table - ) - return table + "\n" + plt.legend() + plt.xlabel("Time (s)") + plt.ylabel("Concentration (M)") + plt.show() def main(): """Command-line interface.""" + console = Console(width=max(105, shutil.get_terminal_size()[0])) + levels = [logging.WARNING, logging.INFO, logging.DEBUG] + # TODO(schneiderfelipe): test and docs parser = argparse.ArgumentParser( description="Interface for building and modifying models." @@ -560,6 +545,29 @@ def main(): # files. This might be useful for some of the more complex operations I # want to be able to do in the future. args = parser.parse_args() + + console.print( + Markdown( + f""" +# overreact {__version__} +Construct precise chemical microkinetic models from first principles + +Inputs: +- Path = {args.path} +- Quantities = {args.quantities} +- Verbose level = {args.verbose} +- Compile? = {args.compile} +- Plot? = {args.plot} +- QRRHO? = {args.qrrho} +- Temperature = {args.temperature} K +- Pressure = {args.pressure} Pa + +Parsing and calculating… + """ + ), + justify="center", + ) + logging.basicConfig( level=levels[min(len(levels) - 1, args.verbose)], stream=sys.stdout ) @@ -567,16 +575,15 @@ def main(): handler.setFormatter(io.InterfaceFormatter("%(message)s")) model = io.parse_model(args.path, force_compile=args.compile) - print( - summarize_model( - model, - quantities=args.quantities, - savepath=os.path.splitext(args.path)[0] + ".csv", - plot=args.plot, - qrrho=args.qrrho, - temperature=args.temperature, - ) + report = Report( + model, + quantities=args.quantities, + savepath=os.path.splitext(args.path)[0] + ".csv", + plot=args.plot, + qrrho=args.qrrho, + temperature=args.temperature, ) + console.print(report, justify="center") if __name__ == "__main__": diff --git a/overreact/misc.py b/overreact/misc.py index 078851f7..1a13fbfb 100644 --- a/overreact/misc.py +++ b/overreact/misc.py @@ -448,6 +448,8 @@ def _get_chemical( ): """Wrap `thermo.Chemical`. + This function is used for obtaining property values and requires the `thermo` package. + All parameters are passed to `thermo.Chemical` and the returned object is returned. diff --git a/overreact/rates.py b/overreact/rates.py index 5623bb80..c3e2edd3 100644 --- a/overreact/rates.py +++ b/overreact/rates.py @@ -13,7 +13,7 @@ def liquid_viscosity(id, temperature=298.15, pressure=constants.atm): """Dynamic viscosity of a solvent. - This function uses the `thermo` package. + This function requires the `thermo` package for obtaining property values. Parameters ---------- diff --git a/requirements.txt b/requirements.txt index 6db59adc..f7405bfc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ -cclib==1.6.3 +cclib>=1.6.3 scipy>=1.4.0 matplotlib>=2.1.1 pytest>=5.2.1 +rich>=9.2.0 thermo>=0.1.39 diff --git a/setup.py b/setup.py index 813d9bba..ad2838aa 100644 --- a/setup.py +++ b/setup.py @@ -8,13 +8,13 @@ setup( name="overreact", version=open("VERSION").read().strip(), - description="A Python package for constructing microkinetic models", + description="Construct precise chemical microkinetic models from first principles", author="Felipe Silveira de Souza Schneider", author_email="schneider.felipe@posgrad.ufsc.br", url="https://github.com/schneiderfelipe/overreact", packages=find_packages(), entry_points={"console_scripts": ["overreact = overreact.cli:main"]}, install_requires=["cclib>=1.6.3", "scipy>=1.4.0"], - extras_require={"thermo": ["thermo>=0.1.39"]}, - tests_require=["pytest>=5.2.1"], + extras_require={"cli": ["rich>=9.2.0"], "solvents": ["thermo>=0.1.39"],}, + tests_require=["matplotlib>=2.1.1", "pytest>=5.2.1"], ) diff --git a/tests/test_thermo_gas.py b/tests/test_thermo_gas.py index a86ae18c..e32221ba 100644 --- a/tests/test_thermo_gas.py +++ b/tests/test_thermo_gas.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 -"""Tests for gas phase using the thermo module.""" +"""Tests for gas phase using the `thermo` module.""" import numpy as np import pytest diff --git a/tests/test_thermo_solv.py b/tests/test_thermo_solv.py index 99d24a50..b5fa9f25 100644 --- a/tests/test_thermo_solv.py +++ b/tests/test_thermo_solv.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 -"""Tests for solvation using the thermo module.""" +"""Tests for solvation using the `thermo` module.""" import pytest