From 34be51d61e409c5910c531113209ffbc0bd1be41 Mon Sep 17 00:00:00 2001 From: "Felipe S. S. Schneider" <37125+schneiderfelipe@users.noreply.github.com> Date: Thu, 26 Nov 2020 13:27:18 -0300 Subject: [PATCH] Use bias in the CLI and make model objects immutable (#66) * Use bias in the CLI * Make model objects immutable and improve dotdict --- data | 2 +- docsrc/quickstart.rst | 6 +- docsrc/tutorials/1-rate_constant.rst | 10 +- docsrc/tutorials/2-equilibrium-constants.rst | 4 +- docsrc/tutorials/4-curtin-hammett.rst | 2 +- .../degree-of-rate-control.rst | 4 +- .../1-basic-reaction-simulation.rst | 16 +- overreact/api.py | 109 ++--- overreact/cli.py | 117 ++---- overreact/core.py | 373 +++++++++--------- overreact/io.py | 345 +++++++++------- overreact/misc.py | 38 +- tests/test_core.py | 46 +-- tests/test_regressions.py | 14 +- 14 files changed, 579 insertions(+), 507 deletions(-) diff --git a/data b/data index b1581f78..0b9562c8 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit b1581f78dbe3f20e169842605974dc07791d7491 +Subproject commit 0b9562c8077a92b0c66698822dcb11d9864e9f83 diff --git a/docsrc/quickstart.rst b/docsrc/quickstart.rst index bab8cc4d..f1220b4d 100644 --- a/docsrc/quickstart.rst +++ b/docsrc/quickstart.rst @@ -7,9 +7,9 @@ any thinkable reaction model: >>> from overreact import api >>> scheme = api.parse_reactions("S -> E‡ -> S") >>> scheme -Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], +Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), ...) The "‡" symbol is used to indicate transition states (but the "#" symbol is diff --git a/docsrc/tutorials/1-rate_constant.rst b/docsrc/tutorials/1-rate_constant.rst index f40a98e2..6aef44fb 100644 --- a/docsrc/tutorials/1-rate_constant.rst +++ b/docsrc/tutorials/1-rate_constant.rst @@ -58,11 +58,11 @@ is a transformation from absolute state energies to deltas of energy relative to the reactions in the model: >>> scheme.B -[[-1., -1.], - [1., 0.], - [0., 0.], - [0., 1.], - [0., 0.]] +((-1., -1.), + (1., 0.), + (0., 0.), + (0., 1.), + (0., 0.)) >>> api.get_delta(scheme.B, [1.92, 19.18, 2.15, 20.78, 3.40]) array([17.26, 18.86]) diff --git a/docsrc/tutorials/2-equilibrium-constants.rst b/docsrc/tutorials/2-equilibrium-constants.rst index cc0da03b..20515b01 100644 --- a/docsrc/tutorials/2-equilibrium-constants.rst +++ b/docsrc/tutorials/2-equilibrium-constants.rst @@ -18,8 +18,8 @@ So let's check it: ... Cd2p + 4 MeNH2 <=> [Cd(MeNH2)4]2p ... """) >>> scheme.compounds, scheme.reactions -(['Cd2p', 'MeNH2', '[Cd(MeNH2)4]2p'], - ['Cd2p + 4 MeNH2 -> [Cd(MeNH2)4]2p', '[Cd(MeNH2)4]2p -> Cd2p + 4 MeNH2']) +(('Cd2p', 'MeNH2', '[Cd(MeNH2)4]2p'), + ('Cd2p + 4 MeNH2 -> [Cd(MeNH2)4]2p', '[Cd(MeNH2)4]2p -> Cd2p + 4 MeNH2')) >>> dydt = simulate.get_dydt(scheme, [K, 1.]) >>> y, r = simulate.get_y(dydt, y0=[0., 0., 1.]) >>> y(y.t_max) diff --git a/docsrc/tutorials/4-curtin-hammett.rst b/docsrc/tutorials/4-curtin-hammett.rst index 8b828045..85dfca4b 100644 --- a/docsrc/tutorials/4-curtin-hammett.rst +++ b/docsrc/tutorials/4-curtin-hammett.rst @@ -18,7 +18,7 @@ The selectivity is defined as the ratio between both products: S = \frac{c_{\ce{P_1}}}{c_{\ce{P_2}}} >>> scheme.compounds -['I1', 'I2', 'T1‡', 'P1', 'T2‡', 'P2'] +('I1', 'I2', 'T1‡', 'P1', 'T2‡', 'P2') >>> freeenergy = [0.0, -0.5, 10.0, 1.0, 11.0, 1.0] >>> k = rates.eyring(constants.kcal * api.get_delta(scheme.B, freeenergy)) >>> k diff --git a/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst b/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst index 6b2072f4..29e41113 100644 --- a/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst +++ b/docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst @@ -12,9 +12,9 @@ Here we try to reproduce the analysis on the Langmuir-Hinshelwood mechanism of ... AB* <=> AB + * # fast to equilibrium ... """) >>> scheme.compounds -['A', '*', 'A*', 'B', 'B*', 'AB*‡', 'AB*', 'AB'] +('A', '*', 'A*', 'B', 'B*', 'AB*‡', 'AB*', 'AB') >>> scheme.reactions -['A + * -> A*', 'A* -> A + *', 'B + * -> B*', 'B* -> B + *', 'A* + B* -> AB* + *', 'AB* -> AB + *', 'AB + * -> AB*'] +('A + * -> A*', 'A* -> A + *', 'B + * -> B*', 'B* -> B + *', 'A* + B* -> AB* + *', 'AB* -> AB + *', 'AB + * -> AB*') where ``*`` denotes a free surface site. diff --git a/docsrc/tutorials/simulation/1-basic-reaction-simulation.rst b/docsrc/tutorials/simulation/1-basic-reaction-simulation.rst index fd4af588..fddd213f 100644 --- a/docsrc/tutorials/simulation/1-basic-reaction-simulation.rst +++ b/docsrc/tutorials/simulation/1-basic-reaction-simulation.rst @@ -18,11 +18,11 @@ The above model translates to the following in overreact: Among the return values we have compounds and reactions: >>> scheme.compounds -['E', 'S', 'ES', 'ES‡', 'P'] +('E', 'S', 'ES', 'ES‡', 'P') >>> scheme.reactions -['E + S -> ES', +('E + S -> ES', 'ES -> E + S', - 'ES -> E + P'] + 'ES -> E + P') All compounds and reactions are parsed in the order they appear in the model. Observe that the :math:`\ce{ES^{\ddagger}}` transition state was correctly @@ -90,11 +90,11 @@ matrix, whose entry :math:`A_{ij}` stores the coefficient of the i-th compound in the j-th reaction: >>> scheme.A -[[-1., 1., 1.], - [-1., 1., 0.], - [1., -1., -1.], - [0., 0., 0.], - [0., 0., 1.]] +((-1., 1., 1.), + (-1., 1., 0.), + (1., -1., -1.), + (0., 0., 0.), + (0., 0., 1.)) Observe that negative (positive) entries represent reactants (products) and a line full of zeros represents transition states (which is to say that they are diff --git a/overreact/api.py b/overreact/api.py index 25c6cb4a..68af547b 100644 --- a/overreact/api.py +++ b/overreact/api.py @@ -27,6 +27,7 @@ from overreact._thermo import get_delta from overreact._thermo import get_reaction_entropies from overreact._thermo import change_reference_state +from overreact.misc import cache logger = logging.getLogger(__name__) @@ -262,62 +263,11 @@ def get_freeenergies( """ enthalpies = get_enthalpies(compounds, qrrho=qrrho, temperature=temperature) entropies = get_entropies(compounds, qrrho=qrrho, temperature=temperature) + # TODO(schneiderfelipe): log the contribution of bias return enthalpies - temperature * entropies + _np.asanyarray(bias) -# TODO(schneiderfelipe): this should probably be deprecated but it is very -# good as a starting point for sensitivity analyses in general. -def get_drc( - scheme, - compounds, - y0, - t_span=None, - method="Radau", - qrrho=True, - scale="l mol-1 s-1", - temperature=298.15, - dx=1.5e3, # joules - order=3, -): - """Calculate the degree of rate control for a single compound. - - Examples - -------- - >>> model = parse_model("data/tanaka1996/UMP2/6-311G(2df,2pd)/model.jk") - """ - x0 = _np.zeros(len(scheme.compounds)) - - def func(t, x=0.0, i=-1): - bias = _np.copy(x0) - bias[i] += x - - k = get_k( - scheme, - compounds=compounds, - bias=bias, - # tunneling=tunneling, - qrrho=qrrho, - scale=scale, - temperature=temperature, - # pressure=pressure, - # delta_freeenergies=delta_freeenergies, - # molecularity=molecularity, - # volume=volume, - ) - _, r = get_y(get_dydt(scheme, k), y0=y0, t_span=t_span, method=method) - - return _np.log(r(t)) - - def drc(t, i=-1): - return ( - -constants.R - * temperature - * _derivative(lambda x: func(t, x, i), x0=0.0, dx=dx, n=1, order=order) - ) - - return drc - - +# @cache def get_k( scheme, compounds=None, @@ -544,3 +494,56 @@ def get_kappa(scheme, compounds, method="eckart", temperature=298.15): f"{', '.join([f'{kappa:7.3g}' for kappa in kappas])}" ) return kappas + + +# TODO(schneiderfelipe): this should probably be deprecated but it is very +# good as a starting point for sensitivity analyses in general. +def get_drc( + scheme, + compounds, + y0, + t_span=None, + method="Radau", + qrrho=True, + scale="l mol-1 s-1", + temperature=298.15, + dx=1.5e3, # joules + order=3, +): + """Calculate the degree of rate control for a single compound. + + Examples + -------- + >>> model = parse_model("data/tanaka1996/UMP2/6-311G(2df,2pd)/model.jk") + """ + x0 = _np.zeros(len(scheme.compounds)) + + def func(t, x=0.0, i=-1): + bias = _np.copy(x0) + bias[i] += x + + k = get_k( + scheme, + compounds=compounds, + bias=bias, + # tunneling=tunneling, + qrrho=qrrho, + scale=scale, + temperature=temperature, + # pressure=pressure, + # delta_freeenergies=delta_freeenergies, + # molecularity=molecularity, + # volume=volume, + ) + _, r = get_y(get_dydt(scheme, k), y0=y0, t_span=t_span, method=method) + + return _np.log(r(t)) + + def drc(t, i=-1): + return ( + -constants.R + * temperature + * _derivative(lambda x: func(t, x, i), x0=0.0, dx=dx, n=1, order=order) + ) + + return drc diff --git a/overreact/cli.py b/overreact/cli.py index c39a1e7a..f543f2d8 100644 --- a/overreact/cli.py +++ b/overreact/cli.py @@ -43,75 +43,25 @@ class Report: Examples -------- + >>> from rich import print >>> model = api.parse_model("data/ethane/B97-3c/model.jk") - >>> 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 - ╵ ╵ ╵ ╵ ╵ + >>> print(Report(model)) + ╭──────────────────╮ + │ (read) reactions │ + │ │ + │ S -> E‡ -> S │ + │ │ + ╰──────────────────╯ + (parsed) reactions + + no reactant(s) via‡ product(s) half equilib.? + ─────────────────────────────────────────────────────── + 0 S E‡ S No + + ... + + Only in the table above, all Gibbs free energies were biased by 0.0 J/mol. + For half-equilibria, only ratios make sense. """ def __init__( @@ -122,6 +72,7 @@ def __init__( plot=False, # TODO(schneiderfelipe): change to do_plot qrrho=True, # TODO(schneiderfelipe): change to use_qrrho temperature=298.15, + bias=0.0, box_style=box.SIMPLE, ): self.model = model @@ -130,6 +81,7 @@ def __init__( self.plot = plot self.qrrho = qrrho self.temperature = temperature + self.bias = bias self.box_style = box_style def __rich_console__(self, console, options): @@ -183,11 +135,11 @@ def _yield_scheme(self): reactants, _, products = re.split(r"\s*(->|<=>|<-)\s*", reaction) # TODO(schneiderfelipe): should we use "No" instead of None for # "half-equilib.?"? - row = [f"{i:d}", reactants, None, products, "NO"] + row = [f"{i:d}", reactants, None, products, "No"] if transition_states[i] is not None: row[2] = scheme.compounds[transition_states[i]] elif scheme.is_half_equilibrium[i]: - row[4] = "YES" + row[4] = "Yes" parsed_table.add_row(*row) yield parsed_table @@ -289,7 +241,6 @@ def _yield_thermochemistry(self): freeenergies, api.get_freeenergies( self.model.compounds, - # bias=bias, qrrho=self.qrrho, temperature=self.temperature, # pressure=pressure, @@ -434,11 +385,15 @@ def _yield_kinetics(self): renderable """ # TODO(schneiderfelipe): apply other corrections to k (such as - # diffusion control) + # diffusion control). + # TODO(schneiderfelipe): use pressure. + # TODO(schneiderfelipe): support changing the type of tunneling + # correction. k = { "M⁻ⁿ⁺¹·s⁻¹": api.get_k( self.model.scheme, self.model.compounds, + bias=self.bias, qrrho=self.qrrho, temperature=self.temperature, scale="l mol-1 s-1", @@ -446,6 +401,7 @@ def _yield_kinetics(self): "(cm³/particle)ⁿ⁻¹·s⁻¹": api.get_k( self.model.scheme, self.model.compounds, + bias=self.bias, qrrho=self.qrrho, temperature=self.temperature, scale="cm3 particle-1 s-1", @@ -453,6 +409,7 @@ def _yield_kinetics(self): "atm⁻ⁿ⁺¹·s⁻¹": api.get_k( self.model.scheme, self.model.compounds, + bias=self.bias, qrrho=self.qrrho, temperature=self.temperature, scale="atm-1 s-1", @@ -472,12 +429,17 @@ def _yield_kinetics(self): box=self.box_style, ) for i, reaction in enumerate(self.model.scheme.reactions): - row = [f"{i:d}", reaction, "NO"] + [f"{k[scale][i]:.3g}" for scale in k] + row = [f"{i:d}", reaction, "No"] + [f"{k[scale][i]:.3g}" for scale in k] if self.model.scheme.is_half_equilibrium[i]: - row[2] = "YES" + row[2] = "Yes" kinetics_table.add_row(*row) yield kinetics_table + yield Markdown( + "Only in the table above, all Gibbs free energies were biased by " + f"{self.bias} J/mol." + ) + yield Markdown("For **half-equilibria**, only ratios make sense.") if self.concentrations is not None and self.concentrations: scale = "M⁻ⁿ⁺¹·s⁻¹" @@ -499,7 +461,8 @@ def _yield_kinetics(self): quantity = float(fields[1]) except (IndexError, ValueError): raise ValueError( - f"badly formatted concentrations: '{' '.join(self.concentrations)}'" + "badly formatted concentrations: " + f"'{' '.join(self.concentrations)}'" ) # TODO(schneiderfelipe): the following is inefficient but probably OK @@ -552,6 +515,7 @@ def main(): levels = [logging.WARNING, logging.INFO, logging.DEBUG] # TODO(schneiderfelipe): test and docs + # TODO(schneiderfelipe): improve help page parser = argparse.ArgumentParser( description="Interface for building and modifying models." ) @@ -564,6 +528,7 @@ def main(): ), nargs="*", ) + parser.add_argument("-b", "--bias", type=float, default=0.0) parser.add_argument("-T", "--temperature", type=float, default=298.15) # TODO(schneiderfelipe): support pressure specification! parser.add_argument("-p", "--pressure", type=float, default=constants.atm) @@ -616,6 +581,7 @@ def main(): - QRRHO? = {args.qrrho} - Temperature = {args.temperature} K - Pressure = {args.pressure} Pa +- Bias = {args.bias / constants.kcal} kcal/mol Parsing and calculating… """ @@ -637,6 +603,7 @@ def main(): plot=args.plot, qrrho=args.qrrho, temperature=args.temperature, + bias=args.bias, ) console.print(report, justify="left") diff --git a/overreact/core.py b/overreact/core.py index 21a9da32..3228cde1 100644 --- a/overreact/core.py +++ b/overreact/core.py @@ -8,6 +8,8 @@ import numpy as np +from overreact import misc as _misc + Scheme = namedtuple("Scheme", "compounds reactions is_half_equilibrium A B") _abbr_environment = { @@ -39,17 +41,17 @@ def _check_scheme(scheme_or_text): Examples -------- >>> _check_scheme("A -> B") - Scheme(compounds=['A', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [1.]], - B=[[-1.], [1.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (1.,)), + B=((-1.,), (1.,))) >>> _check_scheme(_check_scheme("A -> B")) - Scheme(compounds=['A', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [1.]], - B=[[-1.], [1.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (1.,)), + B=((-1.,), (1.,))) """ if isinstance(scheme_or_text, Scheme): return scheme_or_text @@ -72,48 +74,48 @@ def get_transition_states(A, B, is_half_equilibrium): -------- >>> scheme = parse_reactions("A -> B") >>> print(scheme) - Scheme(compounds=['A', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [1.]], - B=[[-1.], [1.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (1.,)), + B=((-1.,), (1.,))) >>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium) - [None] + (None,) >>> scheme = parse_reactions("S -> E‡ -> S") >>> print(scheme) - Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], - A=[[0.], [0.]], - B=[[-1.], [1.]]) + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((0.,), (0.,)), + B=((-1.,), (1.,))) >>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium) - [1] + (1,) >>> scheme = parse_reactions("E + S <=> ES -> ES‡ -> E + P") >>> print(scheme) - Scheme(compounds=['E', 'S', 'ES', 'ES‡', 'P'], - reactions=['E + S -> ES', 'ES -> E + S', 'ES -> E + P'], - is_half_equilibrium=[True, True, False], - A=[[-1., 1., 1.], - [-1., 1., 0.], - [1., -1., -1.], - [0., 0., 0.], - [0., 0., 1.]], - B=[[-1., 0., 0.], - [-1., 0., 0.], - [1., 0., -1.], - [0., 0., 1.], - [0., 0., 0.]]) + Scheme(compounds=('E', 'S', 'ES', 'ES‡', 'P'), + reactions=('E + S -> ES', 'ES -> E + S', 'ES -> E + P'), + is_half_equilibrium=(True, True, False), + A=((-1., 1., 1.), + (-1., 1., 0.), + (1., -1., -1.), + (0., 0., 0.), + (0., 0., 1.)), + B=((-1., 0., 0.), + (-1., 0., 0.), + (1., 0., -1.), + (0., 0., 1.), + (0., 0., 0.))) >>> get_transition_states(scheme.A, scheme.B, scheme.is_half_equilibrium) - [None, None, 3] + (None, None, 3) """ tau = np.asanyarray(B) - np.asanyarray(A) > 0 # transition state matrix - return [ + return tuple( x if not is_half_equilibrium[i] and tau[:, i].any() else None for i, x in enumerate(np.argmax(tau, axis=0)) - ] + ) # TODO(schneiderfelipe): some of the more esoteric doctests should become real @@ -137,80 +139,87 @@ def unparse_reactions(scheme): Examples -------- - >>> unparse_reactions(Scheme(compounds=['A', 'B'], reactions=['A -> B'], - ... is_half_equilibrium=np.array([False]), - ... A=np.array([[-1.], - ... [ 1.]]), - ... B=np.array([[-1.], - ... [ 1.]]))) + >>> unparse_reactions(Scheme(compounds=('A', 'B'), reactions=('A -> B',), + ... is_half_equilibrium=(False,), + ... A=((-1.,), + ... ( 1.,)), + ... B=((-1.,), + ... ( 1.,)))) 'A -> B' - >>> unparse_reactions(Scheme(compounds=['A', 'B'], - ... reactions=['A -> B', 'B -> A'], - ... is_half_equilibrium=np.array([ True, True]), - ... A=np.array([[-1., 1.], - ... [ 1., -1.]]), - ... B=np.array([[-1., 0.], - ... [ 1., 0.]]))) + >>> unparse_reactions(Scheme(compounds=('A', 'B'), + ... reactions=('A -> B', 'B -> A'), + ... is_half_equilibrium=(True, True), + ... A=((-1., 1.), + ... ( 1., -1.)), + ... B=((-1., 0.), + ... ( 1., 0.)))) 'A <=> B' - >>> unparse_reactions(Scheme(compounds=['A', 'A‡', 'B'], - ... reactions=['A -> B'], - ... is_half_equilibrium=np.array([False]), - ... A=np.array([[-1.], - ... [ 0.], - ... [ 1.]]), - ... B=np.array([[-1.], - ... [ 1.], - ... [ 0.]]))) + >>> unparse_reactions(Scheme(compounds=('A', 'A‡', 'B'), + ... reactions=('A -> B',), + ... is_half_equilibrium=(False,), + ... A=((-1.,), + ... ( 0.,), + ... ( 1.,)), + ... B=((-1.,), + ... ( 1.,), + ... ( 0.,)))) 'A -> A‡ -> B' - >>> print(unparse_reactions(Scheme(compounds=['B', 'B‡', 'C', 'D', "B'‡", - ... 'E', 'A'], - ... reactions=['B -> C', 'B -> D', 'B -> E', - ... 'A -> C', 'A -> D'], - ... is_half_equilibrium=np.array([False, False, - ... False, False, - ... False]), - ... A=np.array([[-1., -1., -1., 0., 0.], - ... [ 0., 0., 0., 0., 0.], - ... [ 1., 0., 0., 1., 0.], - ... [ 0., 1., 0., 0., 1.], - ... [ 0., 0., 0., 0., 0.], - ... [ 0., 0., 1., 0., 0.], - ... [ 0., 0., 0., -1., -1.]]), - ... B=np.array([[-1., -1., -1., 0., 0.], - ... [ 1., 1., 0., 1., 1.], - ... [ 0., 0., 0., 0., 0.], - ... [ 0., 0., 0., 0., 0.], - ... [ 0., 0., 1., 0., 0.], - ... [ 0., 0., 0., 0., 0.], - ... [ 0., 0., 0., -1., -1.]])))) + >>> print(unparse_reactions(Scheme(compounds=('B', 'B‡', 'C', 'D', "B'‡", + ... 'E', 'A'), + ... reactions=('B -> C', 'B -> D', 'B -> E', + ... 'A -> C', 'A -> D'), + ... is_half_equilibrium=(False, False, False, + ... False, False), + ... A=((-1., -1., -1., 0., 0.), + ... ( 0., 0., 0., 0., 0.), + ... ( 1., 0., 0., 1., 0.), + ... ( 0., 1., 0., 0., 1.), + ... ( 0., 0., 0., 0., 0.), + ... ( 0., 0., 1., 0., 0.), + ... ( 0., 0., 0., -1., -1.)), + ... B=((-1., -1., -1., 0., 0.), + ... ( 1., 1., 0., 1., 1.), + ... ( 0., 0., 0., 0., 0.), + ... ( 0., 0., 0., 0., 0.), + ... ( 0., 0., 1., 0., 0.), + ... ( 0., 0., 0., 0., 0.), + ... ( 0., 0., 0., -1., -1.))))) B -> B‡ -> C B -> B‡ -> D B -> B'‡ -> E A -> B‡ -> C A -> B‡ -> D - >>> print(unparse_reactions(Scheme(compounds=['A', 'A‡', 'B'], - ... reactions=['A -> B', 'A -> B'], - ... is_half_equilibrium=np.array([False, False]), - ... A=np.array([[-1., -1.], - ... [ 0., 0.], - ... [ 1., 1.]]), - ... B=np.array([[-1., -1.], - ... [ 1., 0.], - ... [ 0., 1.]])))) + >>> print(unparse_reactions(Scheme(compounds=('A', 'A‡', 'B'), + ... reactions=('A -> B', 'A -> B'), + ... is_half_equilibrium=(False, False), + ... A=((-1., -1.), + ... ( 0., 0.), + ... ( 1., 1.)), + ... B=((-1., -1.), + ... ( 1., 0.), + ... ( 0., 1.))))) A -> A‡ -> B A -> B - >>> unparse_reactions(Scheme(compounds=['A', 'A‡', "A'‡", 'B'], - ... reactions=["A -> A'‡"], - ... is_half_equilibrium=np.array([False]), - ... A=np.array([[-1.], - ... [ 0.], - ... [ 1.], - ... [ 0.]]), - ... B=np.array([[-1.], - ... [ 1.], - ... [ 0.], - ... [ 0.]]))) + >>> unparse_reactions(Scheme(compounds=('A', 'A‡', "A'‡", 'B'), + ... reactions=("A -> A'‡",), + ... is_half_equilibrium=(False,), + ... A=((-1.,), + ... ( 0.,), + ... ( 1.,), + ... ( 0.,)), + ... B=((-1.,), + ... ( 1.,), + ... ( 0.,), + ... ( 0.,)))) "A -> A‡ -> A'‡" + >>> unparse_reactions(Scheme(compounds=('S', 'E‡'), + ... reactions=('S -> S',), + ... is_half_equilibrium=(False,), + ... A=((0.0,), + ... (0.0,)), + ... B=((-1.0,), + ... (1.0,)))) + 'S -> E‡ -> S' """ scheme = _check_scheme(scheme) transition_states = get_transition_states( @@ -400,36 +409,36 @@ def parse_reactions(text): returned object has the following attributes: >>> scheme.compounds - ['A', 'B'] + ('A', 'B') >>> scheme.reactions - ['A -> B'] + ('A -> B',) >>> scheme.is_half_equilibrium - [False] + (False,) >>> scheme.A - [[-1.], [1.]] + ((-1.,), (1.,)) >>> scheme.B - [[-1.], [1.]] + ((-1.,), (1.,)) The same reaction can be specified in reverse order: >>> parse_reactions("B <- A // reverse reaction of the above") - Scheme(compounds=['A', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [1.]], - B=[[-1.], [1.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (1.,)), + B=((-1.,), (1.,))) Equilibria produce twice as many direct reactions, while the B matrix defines an energy relationship for only one of each pair: >>> parse_reactions("A <=> B // an equilibrium") - Scheme(compounds=['A', 'B'], - reactions=['A -> B', 'B -> A'], - is_half_equilibrium=[True, True], - A=[[-1., 1.], - [1., -1.]], - B=[[-1., 0.], - [1., 0.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B', 'B -> A'), + is_half_equilibrium=(True, True), + A=((-1., 1.), + (1., -1.)), + B=((-1., 0.), + (1., 0.))) Adding twice the same reaction results in a single reaction being added. This of course also works with equilibria (extra whitespaces are ignored): @@ -440,13 +449,13 @@ def parse_reactions(text): ... A -> B <- A ... B <- A -> B ... ''') - Scheme(compounds=['A', 'B'], - reactions=['A -> B', 'B -> A'], - is_half_equilibrium=[True, True], - A=[[-1., 1.], - [1., -1.]], - B=[[-1., 0.], - [1., 0.]]) + Scheme(compounds=('A', 'B'), + reactions=('A -> B', 'B -> A'), + is_half_equilibrium=(True, True), + A=((-1., 1.), + (1., -1.)), + B=((-1., 0.), + (1., 0.))) Transition states are specified with an asterisk at the end. They are shown among compounds, but the matrix A ensures they'll never have a non-zero @@ -454,31 +463,31 @@ def parse_reactions(text): matrix: >>> parse_reactions("A -> A‡ -> B") - Scheme(compounds=['A', 'A‡', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [0.], [1.]], - B=[[-1.], [1.], [0.]]) + Scheme(compounds=('A', 'A‡', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (0.,), (1.,)), + B=((-1.,), (1.,), (0.,))) This gives the same result as above: >>> parse_reactions("A -> A‡ -> B <- A‡ <- A") - Scheme(compounds=['A', 'A‡', 'B'], - reactions=['A -> B'], - is_half_equilibrium=[False], - A=[[-1.], [0.], [1.]], - B=[[-1.], [1.], [0.]]) + Scheme(compounds=('A', 'A‡', 'B'), + reactions=('A -> B',), + is_half_equilibrium=(False,), + A=((-1.,), (0.,), (1.,)), + B=((-1.,), (1.,), (0.,))) It is possible to define a reaction whose product is the same as the reactant. This is found in phenomena such as ammonia inversion or the methyl rotation in ethane: >>> parse_reactions("S -> E‡ -> S") - Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], - A=[[0.], [0.]], - B=[[-1.], [1.]]) + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((0.,), (0.,)), + B=((-1.,), (1.,))) As such, a column full of zeros in the A matrix corresponds to a reaction with zero net change. As can be seen, overreact allows for very general @@ -491,23 +500,23 @@ def parse_reactions(text): ... B -> B'‡ -> E // this is a classical competitive reaction ... A -> B‡ ... ''') - Scheme(compounds=['B', 'B‡', 'C', 'D', "B'‡", 'E', 'A'], - reactions=['B -> C', 'B -> D', 'B -> E', 'A -> C', 'A -> D'], - is_half_equilibrium=[False, False, False, False, False], - A=[[-1., -1., -1., 0., 0.], - [0., 0., 0., 0., 0.], - [1., 0., 0., 1., 0.], - [0., 1., 0., 0., 1.], - [0., 0., 0., 0., 0.], - [0., 0., 1., 0., 0.], - [0., 0., 0., -1., -1.]], - B=[[-1., -1., -1., 0., 0.], - [1., 1., 0., 1., 1.], - [0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0.], - [0., 0., 1., 0., 0.], - [0., 0., 0., 0., 0.], - [0., 0., 0., -1., -1.]]) + Scheme(compounds=('B', 'B‡', 'C', 'D', "B'‡", 'E', 'A'), + reactions=('B -> C', 'B -> D', 'B -> E', 'A -> C', 'A -> D'), + is_half_equilibrium=(False, False, False, False, False), + A=((-1., -1., -1., 0., 0.), + (0., 0., 0., 0., 0.), + (1., 0., 0., 1., 0.), + (0., 1., 0., 0., 1.), + (0., 0., 0., 0., 0.), + (0., 0., 1., 0., 0.), + (0., 0., 0., -1., -1.)), + B=((-1., -1., -1., 0., 0.), + (1., 1., 0., 1., 1.), + (0., 0., 0., 0., 0.), + (0., 0., 0., 0., 0.), + (0., 0., 1., 0., 0.), + (0., 0., 0., 0., 0.), + (0., 0., 0., -1., -1.))) The following is a borderline case but both reactions should be considered different since they define different processes: @@ -516,15 +525,15 @@ def parse_reactions(text): ... A -> A‡ -> B ... A -> B ... ''') - Scheme(compounds=['A', 'A‡', 'B'], - reactions=['A -> B', 'A -> B'], - is_half_equilibrium=[False, False], - A=[[-1., -1.], - [0., 0.], - [1., 1.]], - B=[[-1., -1.], - [1., 0.], - [0., 1.]]) + Scheme(compounds=('A', 'A‡', 'B'), + reactions=('A -> B', 'A -> B'), + is_half_equilibrium=(False, False), + A=((-1., -1.), + (0., 0.), + (1., 1.)), + B=((-1., -1.), + (1., 0.), + (0., 1.))) The following is correct bevahior. In fact, the reactions are badly defined: if more than one transition state are chained, the following @@ -536,11 +545,11 @@ def parse_reactions(text): barrier be defined in such a weird case: >>> parse_reactions("A -> A‡ -> A'‡ -> B") - Scheme(compounds=['A', 'A‡', "A'‡", 'B'], - reactions=["A -> A'‡"], - is_half_equilibrium=[False], - A=[[-1.], [0.], [1.], [0.]], - B=[[-1.], [1.], [0.], [0.]]) + Scheme(compounds=('A', 'A‡', "A'‡", 'B'), + reactions=("A -> A'‡",), + is_half_equilibrium=(False,), + A=((-1.,), (0.,), (1.,), (0.,)), + B=((-1.,), (1.,), (0.,), (0.,))) """ # TODO(schneiderfelipe): we rely on the ordering of compounds and, as such, # we can only support Python 3.6 and onward. @@ -634,15 +643,19 @@ def _add_reaction(reactants, products, is_half_equilibrium, transition): _add_reaction(reactants, products, is_half_equilibrium, None) return Scheme( - compounds=list(compounds), - reactions=list(_unparse_reactions(reactions)), - is_half_equilibrium=np.array([reaction[2] for reaction in reactions]).tolist(), - A=np.block( - [[vector, np.zeros(len(compounds) - len(vector))] for vector in A] - ).T.tolist(), - B=np.block( - [[vector, np.zeros(len(compounds) - len(vector))] for vector in B] - ).T.tolist(), + compounds=tuple(compounds), + reactions=tuple(_unparse_reactions(reactions)), + is_half_equilibrium=_misc.totuple([reaction[2] for reaction in reactions]), + A=_misc.totuple( + np.block( + [[vector, np.zeros(len(compounds) - len(vector))] for vector in A] + ).T + ), + B=_misc.totuple( + np.block( + [[vector, np.zeros(len(compounds) - len(vector))] for vector in B] + ).T + ), ) diff --git a/overreact/io.py b/overreact/io.py index 142b73e7..fd9766c0 100644 --- a/overreact/io.py +++ b/overreact/io.py @@ -4,6 +4,7 @@ from collections.abc import MutableMapping from collections import defaultdict +from copy import deepcopy as _deepcopy import json import logging import os @@ -14,7 +15,7 @@ from overreact import constants from overreact import core as _core -from overreact import misc +from overreact import misc as _misc with warnings.catch_warnings(): warnings.simplefilter("ignore") @@ -44,7 +45,7 @@ def parse_model(path, force_compile=False): Returns ------- - model : dict-like + model : immutable dict-like Raises ------ @@ -55,26 +56,27 @@ def parse_model(path, force_compile=False): -------- >>> model = parse_model("data/ethane/B97-3c/model.jk") >>> model.scheme - Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], - A=[[1.], [0.]], - B=[[-1.], [1.]]) + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((0.0,), (0.0,)), + B=((-1.0,), (1.0,))) >>> model.compounds["S"] {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.633588, 2.520693, -4.8e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.633588, 2.520693, -4.8e-05), ..., - [-5.832852, 3.674431, 0.363239]], - 'vibfreqs': [307.57, 825.42, ..., 3071.11, 3071.45], - 'vibdisps': [[[-1.7e-05, 3.4e-05, 5.4e-05], + (-5.832852, 3.674431, 0.363239)), + 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45), + 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05), ..., - [-0.011061, -0.030431, -0.027036]]]} + (-0.011061, -0.030431, -0.027036)))} - >>> model_from_source = parse_model("data/ethane/B97-3c/model.k") + >>> model_from_source = parse_model("data/ethane/B97-3c/model.k", + ... force_compile=True) >>> model_from_source == model True @@ -85,7 +87,7 @@ def parse_model(path, force_compile=False): if not path.endswith((".k", ".jk")): path = path + ".jk" logger.warning(f"assuming jk-file in {path}") - name, ext = os.path.splitext(path) + name, _ = os.path.splitext(path) path_jk = name + ".jk" if not force_compile and os.path.isfile(path_jk): @@ -96,6 +98,7 @@ def parse_model(path, force_compile=False): logger.info(f"parsing k-file in {path_k}") if not os.path.isfile(path_k): raise FileNotFoundError + model = _parse_source(path_k) with open(path_jk, "w") as f: logger.info(f"writing jk-file in {path_jk}") @@ -117,43 +120,43 @@ def _parse_model(file_or_path): Returns ------- - model : dict-like + model : immutable dict-like Examples -------- >>> model = _parse_model("data/ethane/B97-3c/model.jk") >>> model.scheme - Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], - A=[[1.], [0.]], - B=[[-1.], [1.]]) + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((1.,), (0.,)), + B=((-1.,), (1.,))) >>> model.compounds["S"] {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.633588, 2.520693, -4.8e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.633588, 2.520693, -4.8e-05), ..., - [-5.832852, 3.674431, 0.363239]], - 'vibfreqs': [307.57, 825.42, ..., 3071.11, 3071.45], - 'vibdisps': [[[-1.7e-05, 3.4e-05, 5.4e-05], + (-5.832852, 3.674431, 0.363239)), + 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45), + 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05), ..., - [-0.011061, -0.030431, -0.027036]]]} + (-0.011061, -0.030431, -0.027036)))} >>> model.compounds["E‡"] {'logfile': 'data/ethane/B97-3c/eclipsed.out', 'energy': -209472585.3539883, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.640622, 2.51993, -1.6e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.640622, 2.51993, -1.6e-05), ..., - [-5.730333, 2.893778, 0.996894]], - 'vibfreqs': [-298.94, 902.19, ..., 3077.75, 3078.05], - 'vibdisps': [[[-6.7e-05, 2.3e-05, 3.4e-05], + (-5.730333, 2.893778, 0.996894)), + 'vibfreqs': (-298.94, 902.19, ..., 3077.75, 3078.05), + 'vibdisps': (((-6.7e-05, 2.3e-05, 3.4e-05), ..., - [0.133315, 0.086028, 0.35746]]]} + (0.133315, 0.086028, 0.35746)))} """ try: model = json.load(file_or_path) @@ -161,14 +164,10 @@ def _parse_model(file_or_path): with open(file_or_path, "r") as stream: model = json.load(stream) - model = dotdict(model) - if "scheme" in model: - model.scheme = _core.parse_reactions(model.scheme) - if "compounds" in model: - for key in model.compounds: - model.compounds[key] = dotdict(model.compounds[key]) - return model + model["scheme"] = _core.parse_reactions(model["scheme"]) + + return dotdict(model) def _parse_source(file_path_or_str): @@ -183,7 +182,7 @@ def _parse_source(file_path_or_str): Returns ------- - json : str + model : immutable dict-like Notes ----- @@ -192,6 +191,14 @@ def _parse_source(file_path_or_str): Examples -------- >>> model = _parse_source("data/ethane/B97-3c/model.k") + >>> model.scheme + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((0.0,), + (0.0,)), + B=((-1.0,), + (1.0,))) >>> print(_unparse_model(model)) { "scheme": [ @@ -255,16 +262,17 @@ def _parse_source(file_path_or_str): keys of compounds match. This implementation detail is crucial for the proper internal bevahiour of overreact: - >>> model = parse_model("data/perez-soto2020/NoRI/GFN2-xTB/model.k") + >>> model = parse_model("data/perez-soto2020/NoRI/GFN2-xTB/model.k", + ... force_compile=True) >>> model.scheme.compounds - ['Benzaldehyde(dcm)', 'NButylamine(dcm)', 'A_N(dcm)', 'A_N_N(dcm)', + ('Benzaldehyde(dcm)', 'NButylamine(dcm)', 'A_N(dcm)', 'A_N_N(dcm)', 'Water(dcm)', 'A_N_W(dcm)', 'A_N_N_W(dcm)', 'A_N_W_W(dcm)', 'TS1_#(dcm)', 'Hemiaminal(dcm)', 'TS2_#(dcm)', 'I_W(dcm)', 'TS1N_#(dcm)', 'Int_N(dcm)', 'TS2N_#(dcm)', 'I_N_W(dcm)', 'TS1W_#(dcm)', 'Int_W(dcm)', 'TS2W_#(dcm)', 'I_W_W(dcm)', 'TS1NW_#(dcm)', 'Int_N_W(dcm)', 'TS2NW_#(dcm)', 'I_N_W_W(dcm)', 'TS1WW_#(dcm)', 'Int_W_W(dcm)', 'TS2WW_#(dcm)', - 'I_W_W_W(dcm)', 'Imine(dcm)'] - >>> list(model.compounds) == model.scheme.compounds + 'I_W_W_W(dcm)', 'Imine(dcm)') + >>> tuple(model.compounds) == model.scheme.compounds True """ name = None @@ -312,7 +320,7 @@ def _unparse_source(model): Parameters ---------- - model : dict-like + model : immutable dict-like Returns ------- @@ -329,7 +337,7 @@ def _unparse_source(model): $end $compounds S: - logfile=data/ethane/B97-3c/staggered.out + logfile="data/ethane/B97-3c/staggered.out" energy=-209483812.77142256 mult=1 atomnos=[6, 6, 1, 1, 1, 1, 1, 1] @@ -338,7 +346,7 @@ def _unparse_source(model): vibfreqs=[307.57, 825.42, ..., 3071.11, 3071.45] vibdisps=[[[-1.7e-05, 3.4e-05, 5.4e-05], ..., [..., -0.027036]]] E‡: - logfile=data/ethane/B97-3c/eclipsed.out + logfile="data/ethane/B97-3c/eclipsed.out" energy=-209472585.3539883 mult=1 ... @@ -349,37 +357,37 @@ def _unparse_source(model): >>> model = _parse_source(source) >>> model.scheme - Scheme(compounds=['S', 'E‡'], - reactions=['S -> S'], - is_half_equilibrium=[False], - A=[[1.], [0.]], - B=[[-1.], [1.]]) + Scheme(compounds=('S', 'E‡'), + reactions=('S -> S',), + is_half_equilibrium=(False,), + A=((1.,), (0.,)), + B=((-1.,), (1.,))) >>> model.compounds["S"] {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.633588, 2.520693, -4.8e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.633588, 2.520693, -4.8e-05), ..., - [-5.832852, 3.674431, 0.363239]], - 'vibfreqs': [307.57, 825.42, ..., 3071.11, 3071.45], - 'vibdisps': [[[-1.7e-05, 3.4e-05, 5.4e-05], + (-5.832852, 3.674431, 0.363239)), + 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45), + 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05), ..., - [-0.011061, -0.030431, -0.027036]]]} + (-0.011061, -0.030431, -0.027036)))} >>> model.compounds["E‡"] {'logfile': 'data/ethane/B97-3c/eclipsed.out', 'energy': -209472585.3539883, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.640622, 2.51993, -1.6e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.640622, 2.51993, -1.6e-05), ..., - [-5.730333, 2.893778, 0.996894]], - 'vibfreqs': [-298.94, 902.19, ..., 3077.75, 3078.05], - 'vibdisps': [[[-6.7e-05, 2.3e-05, 3.4e-05], + (-5.730333, 2.893778, 0.996894)), + 'vibfreqs': (-298.94, 902.19, ..., 3077.75, 3078.05), + 'vibdisps': (((-6.7e-05, 2.3e-05, 3.4e-05), ..., - [0.133315, 0.086028, 0.35746]]]} + (0.133315, 0.086028, 0.35746)))} """ source = "// generated by overreact\n" if "scheme" in model: @@ -392,7 +400,10 @@ def _unparse_source(model): for compound in model.compounds: source += f" {compound}:\n" for key in model.compounds[compound]: - source += f" {key}={model.compounds[compound][key]}\n" + inline_json = json.dumps( + model.compounds[compound][key], ensure_ascii=False + ) + source += f" {key}={inline_json}\n" source += "$end\n" return source @@ -405,7 +416,7 @@ def _unparse_model(model): Parameters ---------- - model : dict-like + model : immutable dict-like Returns ------- @@ -473,9 +484,11 @@ def _unparse_model(model): } } """ - model = dotdict(model.copy()) + # create a new mutable object to avoid side effects + model = dict(model.copy()) + if "scheme" in model: - model.scheme = _core.unparse_reactions(model.scheme).split("\n") + model["scheme"] = _core.unparse_reactions(model["scheme"]).split("\n") return json.dumps(model, indent=2, ensure_ascii=False) @@ -496,28 +509,28 @@ def _check_compounds(compounds): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.633588, 2.520693, -4.8e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.633588, 2.520693, -4.8e-05), ..., - [-5.832852, 3.674431, 0.363239]], - 'vibfreqs': [307.57, 825.42, ..., 3071.11, 3071.45], - 'vibdisps': [[[-1.7e-05, 3.4e-05, 5.4e-05], + (-5.832852, 3.674431, 0.363239)), + 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45), + 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05), ..., - [-0.011061, -0.030431, -0.027036]]]}} + (-0.011061, -0.030431, -0.027036)))}} >>> _check_compounds(_check_compounds({"S": "data/ethane/B97-3c/staggered.out"})) {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008], - 'atomcoords': [[-7.633588, 2.520693, -4.8e-05], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, 12.011, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008), + 'atomcoords': ((-7.633588, 2.520693, -4.8e-05), ..., - [-5.832852, 3.674431, 0.363239]], - 'vibfreqs': [307.57, 825.42, ..., 3071.11, 3071.45], - 'vibdisps': [[[-1.7e-05, 3.4e-05, 5.4e-05], + (-5.832852, 3.674431, 0.363239)), + 'vibfreqs': (307.57, 825.42, ..., 3071.11, 3071.45), + 'vibdisps': (((-1.7e-05, 3.4e-05, 5.4e-05), ..., - [-0.011061, -0.030431, -0.027036]]]}} + (-0.011061, -0.030431, -0.027036)))}} """ for name in compounds: if isinstance(compounds[name], str): @@ -540,7 +553,7 @@ def parse_compounds(text, path=("",), select=None): Returns ------- - compounds : dict-like + compounds : immutable dict-like Examples -------- @@ -549,7 +562,7 @@ def parse_compounds(text, path=("",), select=None): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}} >>> compounds = parse_compounds({"S": "data/ethane/B97-3c/staggered.out"}) @@ -557,7 +570,7 @@ def parse_compounds(text, path=("",), select=None): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}} >>> compounds = parse_compounds(["S: data/ethane/B97-3c/staggered.out", @@ -566,12 +579,12 @@ def parse_compounds(text, path=("",), select=None): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}, 'E‡': {'logfile': 'data/ethane/B97-3c/eclipsed.out', 'energy': -209472585.3539883, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}} >>> compounds = parse_compounds('''S: data/ethane/B97-3c/staggered.out @@ -580,12 +593,12 @@ def parse_compounds(text, path=("",), select=None): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}, 'E‡': {'logfile': 'data/ethane/B97-3c/eclipsed.out', 'energy': -209472585.3539883, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}} >>> compounds = parse_compounds({"S": "data/ethane/B97-3c/staggered.out", @@ -594,12 +607,12 @@ def parse_compounds(text, path=("",), select=None): {'S': {'logfile': 'data/ethane/B97-3c/staggered.out', 'energy': -209483812.77142256, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}, 'E‡': {'logfile': 'data/ethane/B97-3c/eclipsed.out', 'energy': -209472585.3539883, 'mult': 1, - 'atomnos': [6, 6, 1, 1, 1, 1, 1, 1], + 'atomnos': (6, 6, 1, 1, 1, 1, 1, 1), ...}} """ if isinstance(text, dict): @@ -609,7 +622,7 @@ def parse_compounds(text, path=("",), select=None): except AttributeError: lines = text name = None - compounds = defaultdict(dotdict) + compounds = defaultdict(dict) for line in lines: if ":" in line: name, line = [x.strip() for x in line.split(":", 1)] @@ -626,6 +639,7 @@ def parse_compounds(text, path=("",), select=None): # TODO(schneiderfelipe): test "nested" logfiles, i.e., # DNPLO-CCSD(T)/def2-TZVP(-f)//revPBE-D4-gCP/def2-SVP, etc. success = False + value = value.strip('"') for p in path: try: logger.info(f"trying to read {os.path.join(p, value)}") @@ -645,7 +659,7 @@ def parse_compounds(text, path=("",), select=None): # TODO(schneiderfelipe): this workaround still allow unused compounds # to be parsed! This should change in the future. compounds = {name: compounds[name] for name in select} - return dict(compounds) + return dotdict(compounds) def read_logfile(path): @@ -671,41 +685,40 @@ def read_logfile(path): {'logfile': 'data/symmetries/benzene.out', 'energy': -609176691.0746485, 'mult': 1, - 'atomnos': [6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1], - 'atommasses': [12.011, ..., 1.008], - 'atomcoords': [[1.856873, 1.370646, -0.904202], + 'atomnos': (6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1), + 'atommasses': (12.011, ..., 1.008), + 'atomcoords': ((1.856873, 1.370646, -0.904202), ..., - [0.794845, 1.265216, -1.190139]], - 'vibfreqs': [397.68, 397.72, ..., 3099.78, 3109.35], - 'vibdisps': [[[-0.030923, -0.073676, 0.142038], + (0.794845, 1.265216, -1.190139)), + 'vibfreqs': (397.68, 397.72, ..., 3099.78, 3109.35), + 'vibdisps': (((-0.030923, -0.073676, 0.142038), ..., - [0.400329, 0.039742, 0.107782]]]} + (0.400329, 0.039742, 0.107782)))} >>> read_logfile("data/symmetries/water.out") {'logfile': 'data/symmetries/water.out', 'energy': -200411626.9934847, 'mult': 1, - 'atomnos': [8, 1, 1], - 'atommasses': [15.999, 1.008, 1.008], - 'atomcoords': [[0.178164, 0.697853, 0.297922], ..., [..., -0.407083]], - 'vibfreqs': [1619.07, 3671.68, 3769.05], - 'vibdisps': [[[0.037749, 0.058895, 0.00295], ..., [..., -0.509409]]], - 'hessian': [[0.22070725835, 0.14857971186, -0.20392672101, ...], + 'atomnos': (8, 1, 1), + 'atommasses': (15.999, 1.008, 1.008), + 'atomcoords': ((0.178164, 0.697853, 0.297922), ..., (..., -0.407083)), + 'vibfreqs': (1619.07, 3671.68, 3769.05), + 'vibdisps': (((0.037749, 0.058895, 0.00295), ..., (..., -0.509409))), + 'hessian': ((0.22070725835, 0.14857971186, -0.20392672101, ...), ..., - [..., -0.01263596771, 0.23474251331, 0.28651782768]]} + (..., -0.01263596771, 0.23474251331, 0.28651782768))} """ try: ccdata = ccopen(path).parse() - ccdata.listify() data = { "logfile": path, # This energy may lack dispersion, solvation, correlation, etc. "energy": ccdata.scfenergies[-1] * constants.eV * constants.N_A, "mult": ccdata.mult, - "atomnos": ccdata.atomnos, - "atommasses": ccdata.atommasses, - "atomcoords": ccdata.atomcoords[-1], - "vibfreqs": ccdata.vibfreqs, - "vibdisps": ccdata.vibdisps, + "atomnos": _misc.totuple(ccdata.atomnos), + "atommasses": _misc.totuple(ccdata.atommasses), + "atomcoords": _misc.totuple(ccdata.atomcoords[-1]), + "vibfreqs": _misc.totuple(ccdata.vibfreqs), + "vibdisps": _misc.totuple(ccdata.vibdisps), } # TODO(schneiderfelipe): only run the code below if we know it is an # ORCA logfile. I want the code in this final section to be very @@ -805,29 +818,29 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 {'energy': -609176691.0746485} >>> _read_orca_logfile("data/tanaka1996/UMP2/6-311G(2df,2pd)/Cl·.out") {'energy': ..., - 'hessian': [[...]]} + 'hessian': ((...))} >>> _read_orca_logfile("data/symmetries/benzene.out", minimal=False) {'energy': -609176691.0746485, 'logfile': 'data/symmetries/benzene.out', 'mult': 1, - 'atomnos': [6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1], - 'atomcoords': [[1.856873, 1.370646, -0.904202], + 'atomnos': (6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1), + 'atomcoords': ((1.856873, 1.370646, -0.904202), ..., - [0.794845, 1.265216, -1.190139]], - 'atommasses': [12.011, ..., 1.008], - 'vibfreqs': [397.68, 397.72, ..., 3099.78, 3109.35]} + (0.794845, 1.265216, -1.190139)), + 'atommasses': (12.011, ..., 1.008), + 'vibfreqs': (397.68, 397.72, ..., 3099.78, 3109.35)} >>> _read_orca_logfile("data/tanaka1996/UMP2/6-311G(2df,2pd)/Cl·.out", ... minimal=False) {'energy': -1206916696.9980993, - 'hessian': [[2.344970859e-11, -1.1953313121e-14, -6.1524162557e-17], - [-1.1953313121e-14, 2.1263734888e-11, 9.2241716473e-17], - [-6.1524162557e-17, 9.2241716473e-17, 2.1024502755e-11]], + 'hessian': ((2.344970859e-11, -1.1953313121e-14, -6.1524162557e-17), + (-1.1953313121e-14, 2.1263734888e-11, 9.2241716473e-17), + (-6.1524162557e-17, 9.2241716473e-17, 2.1024502755e-11)), 'logfile': 'data/tanaka1996/UMP2/6-311G(2df,2pd)/Cl·.out', 'mult': 2, - 'atomnos': [17], - 'atomcoords': [[0.0, 0.0, 0.0]], - 'atommasses': [35.453], - 'vibfreqs': array([], dtype=float64)} + 'atomnos': (17,), + 'atomcoords': ((0.0, 0.0, 0.0),), + 'atommasses': (35.453,), + 'vibfreqs': ()} """ # TODO(schneiderfelipe): return an empty dict-like on error. atomcoords = None @@ -859,7 +872,7 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 while len(line) > 1: atom, x, y, z = line.split() if atom[-1] != ">": - atomnos.append(misc.atomic_number[atom]) + atomnos.append(_misc.atomic_number[atom]) atomcoords.append([float(x), float(y), float(z)]) line = next(file) @@ -878,7 +891,7 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 == "> coreless ECP center with (optional) point charge" ): break - no, lb, za, frag, mass, x, y, z = line.split() + _, lb, _, _, mass, x, y, z = line.split() if lb[-1] != ">": atommasses.append(float(mass)) line = next(file) @@ -909,7 +922,7 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 if hessian is None: try: hessian = _read_orca_hess(path.replace(".out", ".hess")) - data.update({"hessian": hessian.tolist()}) + data.update({"hessian": _misc.totuple(hessian)}) except FileNotFoundError: pass @@ -929,20 +942,22 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 for _ in range(n): line = next(file) atom, x, y, z = line.split() - atomnos.append(misc.atomic_number[atom]) + atomnos.append(_misc.atomic_number[atom]) atomcoords.append([float(x), float(y), float(z)]) - data.update({"atomnos": atomnos, "atomcoords": atomcoords}) + data.update( + {"atomnos": _misc.totuple(atomnos), "atomcoords": _misc.totuple(atomcoords)} + ) if atommasses is None: # TODO(schneiderfelipe): show a warning about guessing masses from the # periodic table. atommasses = [] for n in atomnos: - atommasses.append(misc.atomic_mass[n]) - data.update({"atommasses": atommasses}) + atommasses.append(_misc.atomic_mass[n]) + data.update({"atommasses": _misc.totuple(atommasses)}) if vibfreqs is not None: - data.update({"vibfreqs": vibfreqs}) + data.update({"vibfreqs": _misc.totuple(vibfreqs)}) return data @@ -951,21 +966,57 @@ def _read_orca_logfile(path, minimal=True): # noqa: C901 class dotdict(dict): """Access dictionary attributes through dot.notation. + This object is meant to be immutable, so that it can be hashed. + + Raises + ------ + NotImplementedError + If one attempts to change a value. + Examples -------- - >>> mydict = {"val": "it works"} - >>> nested_dict = {"val": "nested works too"} - >>> mydict = dotdict(mydict) + >>> mydict = dotdict({ + ... "val": "it works like a dict", + ... "nested": { + ... "val": "nested works too" + ... }, + ... }) >>> mydict.val - 'it works' - >>> mydict.nested = dotdict(nested_dict) + 'it works like a dict' >>> mydict.nested.val 'nested works too' + + The constructor actually works recursively: + + >>> type(mydict.nested) + """ + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + for key, val in self.items(): + if isinstance(val, (list, np.ndarray)): + super().__setitem__(key, _misc.totuple(val)) + elif isinstance(val, dict): + super().__setitem__(key, dotdict(val)) + __getattr__ = dict.get - __setattr__ = dict.__setitem__ - __delattr__ = dict.__delitem__ + + def __setitem__(self, key, value): + raise NotImplementedError("dotdict objects are immutable") + + # https://stackoverflow.com/a/1151686/4039050 + # https://stackoverflow.com/a/1151705/4039050 + def __hash__(self): + return hash(self._key()) + + # https://stackoverflow.com/a/16162138/4039050 + def _key(self): + return (frozenset(self), frozenset(self.items())) + + def __eq__(self, other): + return self._key() == other._key() # https://stackoverflow.com/a/61144084/4039050 diff --git a/overreact/misc.py b/overreact/misc.py index 5c895321..f3d394d8 100644 --- a/overreact/misc.py +++ b/overreact/misc.py @@ -5,13 +5,14 @@ Ideally, the functions here will be transfered to other modules in the future. """ -from functools import lru_cache +from functools import lru_cache as cache import numpy as np from scipy.stats import cauchy as _cauchy from scipy.stats import norm as _norm from overreact import constants +from overreact import core as _core def _find_package(package): @@ -562,6 +563,39 @@ def broaden_spectrum( return s +# https://stackoverflow.com/a/10016613 +def totuple(a): + """Convert a numpy.array into nested tuples. + + Parameters + ---------- + a : array-like + + Returns + ------- + tuple + + Examples + -------- + >>> array = np.array(((2,2),(2,-2))) + >>> totuple(array) + ((2, 2), (2, -2)) + """ + # we don't touch some types, and this includes namedtuples + if isinstance(a, (int, float, str, _core.Scheme)): + return a + + try: + a = a.tolist() + except AttributeError: + pass + + try: + return tuple(totuple(i) for i in a) + except TypeError: + return a + + def halton(num, dim=None, jump=1, cranley_patterson=True): """Calculate Halton low-discrepancy sequences. @@ -677,7 +711,7 @@ def _is_prime(num): return primes -@lru_cache(maxsize=1000000) +@cache(maxsize=1000000) def _vdc(n, b=2): """Help haltonspace.""" res, denom = 0, 1 diff --git a/tests/test_core.py b/tests/test_core.py index c8d1455a..6424665b 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -10,22 +10,22 @@ def test_parse_works(): """Test parsing of reactions.""" scheme = core.parse_reactions("A -> B // a direct reaction") - assert ["A", "B"] == scheme[0] - assert ["A -> B"] == scheme[1] - assert [False] == scheme[2] + assert ("A", "B") == scheme[0] + assert ("A -> B",) == scheme[1] + assert (False,) == scheme[2] assert np.all(scheme[3] == scheme[4]) assert np.all(np.array([[-1], [1]]) == scheme[3]) scheme = core.parse_reactions("B <- A // reverse reaction of the above") - assert ["A", "B"] == scheme[0] - assert ["A -> B"] == scheme[1] - assert [False] == scheme[2] + assert ("A", "B") == scheme[0] + assert ("A -> B",) == scheme[1] + assert (False,) == scheme[2] assert np.all(scheme[3] == scheme[4]) assert np.all(np.array([[-1], [1]]) == scheme[3]) scheme = core.parse_reactions("A <=> B // an equilibrium") - assert ["A", "B"] == scheme[0] - assert ["A -> B", "B -> A"] == scheme[1] + assert ("A", "B") == scheme[0] + assert ("A -> B", "B -> A") == scheme[1] assert np.all(np.array([True, True]) == scheme[2]) assert np.all(np.array([[-1.0, 1.0], [1.0, -1.0]]) == scheme[3]) assert np.all(np.array([[-1.0, 0.0], [1.0, 0.0]]) == scheme[4]) @@ -36,25 +36,25 @@ def test_parse_works(): A -> B <- A // reactions B <- A -> B""" ) - assert ["A", "B"] == scheme[0] - assert ["A -> B", "B -> A"] == scheme[1] + assert ("A", "B") == scheme[0] + assert ("A -> B", "B -> A") == scheme[1] assert np.all(np.array([True, True]) == scheme[2]) assert np.all(np.array([[-1.0, 1.0], [1.0, -1.0]]) == scheme[3]) assert np.all(np.array([[-1.0, 0.0], [1.0, 0.0]]) == scheme[4]) scheme = core.parse_reactions("A -> A‡ -> B // a transition state") - assert ["A", "A‡", "B"] == scheme[0] - assert ["A -> B"] == scheme[1] - assert [False] == scheme[2] + assert ("A", "A‡", "B") == scheme[0] + assert ("A -> B",) == scheme[1] + assert (False,) == scheme[2] assert np.all(np.array([[-1.0], [0.0], [1.0]]) == scheme[3]) assert np.all(np.array([[-1.0], [1.0], [0.0]]) == scheme[4]) scheme = core.parse_reactions( "A -> A‡ -> B <- A‡ <- A // (should be) same as above" ) - assert ["A", "A‡", "B"] == scheme[0] - assert ["A -> B"] == scheme[1] - assert [False] == scheme[2] + assert ("A", "A‡", "B") == scheme[0] + assert ("A -> B",) == scheme[1] + assert (False,) == scheme[2] assert np.all(np.array([[-1.0], [0.0], [1.0]]) == scheme[3]) assert np.all(np.array([[-1.0], [1.0], [0.0]]) == scheme[4]) @@ -66,8 +66,8 @@ def test_parse_works(): A -> B‡ """ ) - assert ["B", "B‡", "C", "D", "B'‡", "E", "A"] == scheme[0] - assert ["B -> C", "B -> D", "B -> E", "A -> C", "A -> D"] == scheme[1] + assert ("B", "B‡", "C", "D", "B'‡", "E", "A") == scheme[0] + assert ("B -> C", "B -> D", "B -> E", "A -> C", "A -> D") == scheme[1] assert np.all(np.array([False, False, False, False, False]) == scheme[2]) assert np.all( np.array( @@ -103,8 +103,8 @@ def test_parse_works(): A -> A‡ -> B // this is a tricky example A -> B // but it's better to be explicit""" ) - assert ["A", "A‡", "B"] == scheme[0] - assert ["A -> B", "A -> B"] == scheme[1] + assert ("A", "A‡", "B") == scheme[0] + assert ("A -> B", "A -> B") == scheme[1] assert np.all(np.array([False, False]) == scheme[2]) assert np.all(np.array([[-1.0, -1.0], [0.0, 0.0], [1.0, 1.0]]) == scheme[3]) assert np.all(np.array([[-1.0, -1.0], [1.0, 0.0], [0.0, 1.0]]) == scheme[4]) @@ -114,9 +114,9 @@ def test_parse_works(): // the policy is to not chain transition states A -> A‡ -> A'‡ -> B // this is weird""" ) - assert ["A", "A‡", "A'‡", "B"] == scheme[0] - assert ["A -> A'‡"] == scheme[1] - assert [False] == scheme[2] + assert ("A", "A‡", "A'‡", "B") == scheme[0] + assert ("A -> A'‡",) == scheme[1] + assert (False,) == scheme[2] assert np.all(np.array([[-1.0], [0.0], [1.0], [0.0]]) == scheme[3]) assert np.all(np.array([[-1.0], [1.0], [0.0], [0.0]]) == scheme[4]) diff --git a/tests/test_regressions.py b/tests/test_regressions.py index 1c89237c..70f767a7 100644 --- a/tests/test_regressions.py +++ b/tests/test_regressions.py @@ -15,12 +15,12 @@ def test_solubility_of_acetic_acid(): """Reproduce literature data for AcOH(g) <=> AcOH(aq). - Data is as cited in doi:10.1021/jp810292n and is experimental except when - otherwise indicated in the comments. + Data is as cited in doi:10.1021/jp810292n and doi:10.1063/1.1416902, and + is experimental except when otherwise indicated in the comments. """ model = api.parse_model("data/acetate/model.k") temperature = 298.15 - pK = 4.76 # doi:10.1063/1.1416902 + pK = 4.756 # doi:10.1063/1.1416902 enthalpies = api.get_enthalpies(model.compounds, temperature=temperature) entropies = api.get_entropies(model.compounds, temperature=temperature) @@ -33,17 +33,21 @@ def test_solubility_of_acetic_acid(): delta_freeenergies = delta_enthalpies - temperature * delta_entropies assert delta_freeenergies / constants.kcal == pytest.approx( - [6.50, -6.50, -7.94, 7.94, -74.77, 74.77], 3e-4 + [6.49, -6.49, -7.94, 7.94, -74.77, 74.77], 5e-4 ) concentration_correction = temperature * api.change_reference_state( sign=-1.0, temperature=temperature, ) - # doi:10.1021/jp810292n + + # the following tests the solvation free energy from doi:10.1021/jp810292n assert delta_freeenergies[2] / constants.kcal == pytest.approx( -6.70 + concentration_correction / constants.kcal, 8e-2 ) + + # the following tests the reaction free energy from doi:10.1063/1.1416902 + assert delta_freeenergies[0] == pytest.approx(27.147 * constants.kilo, 8e-4) assert delta_freeenergies[0] == pytest.approx( -constants.R * temperature * np.log(10 ** -pK), 7e-4 )