From 72f6f666e9840bab6f325e793422e028eb328834 Mon Sep 17 00:00:00 2001 From: "Felipe S. S. Schneider" Date: Mon, 2 Oct 2023 13:19:00 -0300 Subject: [PATCH] refactor(overreact): clean up public api code (2) --- overreact/api.py | 110 +++++++++++++++---------------- overreact/coords.py | 153 ++++++++++++++++++++++---------------------- overreact/core.py | 50 +++++++-------- 3 files changed, 153 insertions(+), 160 deletions(-) diff --git a/overreact/api.py b/overreact/api.py index a105822d..c8d74180 100644 --- a/overreact/api.py +++ b/overreact/api.py @@ -1,11 +1,9 @@ -#!/usr/bin/env python3 # noqa: EXE001 - """ -This module contains the high-level application programming interface. +Module containing the high-level application programming interface. If you intend to use **overreact** as a library in a project, you should probably start here. -""" # noqa: D404 +""" from __future__ import annotations @@ -22,7 +20,7 @@ import logging import warnings -from typing import Optional, Union +from typing import TYPE_CHECKING import numpy as np from scipy.misc import derivative @@ -30,15 +28,17 @@ import overreact as rx from overreact import _constants as constants from overreact import coords, rates, tunnel -from overreact.core import Scheme # noqa: TCH001 + +if TYPE_CHECKING: + from overreact.core import Scheme logger = logging.getLogger(__name__) def get_internal_energies( compounds: dict, - qrrho: bool = True, # noqa: FBT001, FBT002 - temperature: float = 298.15, # noqa: RUF100 + qrrho: bool = True, + temperature: float = 298.15, ): """Obtain internal energies for compounds at a given temperature. @@ -70,10 +70,10 @@ def get_internal_energies( array([0. , 2.20053981]) """ - compounds = rx.io._check_compounds(compounds) # noqa: SLF001 + compounds = rx.io._check_compounds(compounds) internal_energies = [] for name in compounds: - logger.info(f"calculate internal energy: {name}") # noqa: G004 + logger.info(f"calculate internal energy: {name}") # TODO(schneiderfelipe): inertia might benefit from caching moments, _, _ = coords.inertia( @@ -95,9 +95,9 @@ def get_internal_energies( def get_enthalpies( compounds: dict, - qrrho: bool = True, # noqa: FBT001, FBT002 - temperature: float = 298.15, # noqa: RUF100 -): # noqa: RUF100 + qrrho: bool = True, + temperature: float = 298.15, +): """Obtain enthalpies for compounds at a given temperature. Parameters @@ -136,10 +136,10 @@ def get_enthalpies( >>> (enthalpies - zero_enthalpies) / constants.kcal array([2.78, 2.50]) """ - compounds = rx.io._check_compounds(compounds) # noqa: SLF001 + compounds = rx.io._check_compounds(compounds) enthalpies = [] for name in compounds: - logger.info(f"calculate enthalpy: {name}") # noqa: G004 + logger.info(f"calculate enthalpy: {name}") # TODO(schneiderfelipe): inertia might benefit from caching moments, _, _ = coords.inertia( @@ -159,11 +159,11 @@ def get_enthalpies( return np.array(enthalpies) -def get_entropies( # noqa: PLR0913 +def get_entropies( compounds: dict, - environment: Optional[str] = None, # noqa: UP007 + environment: str | None = None, method: str = "standard", - qrrho: bool = True, # noqa: FBT001, FBT002 + qrrho: bool = True, temperature: float = 298.15, pressure: float = constants.atm, ): @@ -220,10 +220,10 @@ def get_entropies( # noqa: PLR0913 >>> (sol_entropies - entropies) / constants.calorie array([-6.35360874, -6.35360874]) """ - compounds = rx.io._check_compounds(compounds) # noqa: SLF001 + compounds = rx.io._check_compounds(compounds) entropies = [] for name in compounds: - logger.info(f"calculate entropy: {name}") # noqa: G004 + logger.info(f"calculate entropy: {name}") if "point_group" in compounds[name]: point_group = compounds[name].point_group @@ -241,7 +241,7 @@ def get_entropies( # noqa: PLR0913 ) if environment is None: - environment = rx.core._get_environment(name) # noqa: SLF001 + environment = rx.core._get_environment(name) entropy = rx.thermo.calc_entropy( atommasses=compounds[name].atommasses, atomcoords=compounds[name].atomcoords, @@ -272,8 +272,8 @@ def get_entropies( # noqa: PLR0913 def _check_qrrho( - qrrho: Union[bool, tuple[bool, bool]], # noqa: UP007 -) -> tuple[bool, bool]: # noqa: RUF100 + qrrho: bool | tuple[bool, bool], +) -> tuple[bool, bool]: """Get options for QRRHO for both enthalpy and entropy. Parameters @@ -308,22 +308,21 @@ def _check_qrrho( """ if qrrho is True: return True, True - elif qrrho is False: # noqa: RET505 + elif qrrho is False: return False, False elif isinstance(qrrho, tuple): return qrrho else: - raise ValueError( # noqa: TRY003 - f"unrecognized QRRHO specification: {qrrho}", # noqa: EM102 - ) + msg = f"unrecognized QRRHO specification: {qrrho}" + raise ValueError(msg) -def get_freeenergies( # noqa: PLR0913 +def get_freeenergies( compounds: dict, bias: float = 0.0, - environment: Optional[str] = None, # noqa: UP007 + environment: str | None = None, method: str = "standard", - qrrho: Union[bool | tuple[bool, bool]] = True, # noqa: FBT002, UP007 + qrrho: bool | tuple[bool, bool] = True, temperature: float = 298.15, pressure: float = constants.atm, ): @@ -413,18 +412,18 @@ def get_freeenergies( # noqa: PLR0913 return enthalpies - temperature * entropies + np.asarray(bias) -def get_k( # noqa: PLR0913 +def get_k( scheme: Scheme, - compounds: Optional[dict] = None, # noqa: UP007 + compounds: dict | None = None, bias: float = 0.0, tunneling: str = "eckart", - qrrho: Union[bool | tuple[bool, bool]] = True, # noqa: FBT002, UP007 + qrrho: bool | tuple[bool, bool] = True, scale: str = "l mol-1 s-1", temperature: float = 298.15, pressure: float = constants.atm, - delta_freeenergies: Optional[float] = None, # noqa: UP007 - molecularity: Optional[float] = None, # noqa: UP007 - volume: Optional[float] = None, # noqa: UP007 + delta_freeenergies: float | None = None, + molecularity: float | None = None, + volume: float | None = None, ) -> float: r"""Obtain reaction rate constants for a given reaction scheme. @@ -540,11 +539,11 @@ def get_k( # noqa: PLR0913 ... bias=np.array([0.0, 0.0, -1.4, 0.0, 0.0]) * constants.kcal, ... temperature=300.0, scale="cm3 particle-1 s-1") array([1.1e-12]) - """ # noqa: E501 + """ qrrho_enthalpy, qrrho_entropy = _check_qrrho(qrrho) - scheme = rx.core._check_scheme(scheme) # noqa: SLF001 + scheme = rx.core._check_scheme(scheme) if compounds is not None: - compounds = rx.io._check_compounds(compounds) # noqa: SLF001 + compounds = rx.io._check_compounds(compounds) if delta_freeenergies is None: assert compounds is not None, "compounds could not be inferred" freeenergies = get_freeenergies( @@ -586,13 +585,12 @@ def get_k( # noqa: PLR0913 while i < len(scheme.is_half_equilibrium): if scheme.is_half_equilibrium[i]: pair = k[i : i + 2] - _K = pair[0] / pair[1] # noqa: N806 + _K = pair[0] / pair[1] denom = pair.min() - if denom == 0.0: # noqa: PLR2004 + if denom == 0.0: logger.warning( - "found half-equilibrium reaction with zero rate constant: " - "skipping equilibrium normalization", + "found half-equilibrium reaction with zero rate constant: skipping equilibrium normalization", ) denom = 1.0 @@ -608,7 +606,7 @@ def get_k( # noqa: PLR0913 i += 1 logger.info( - "(classical) reaction rate constants: " # noqa: G004 + "(classical) reaction rate constants: " f"{', '.join([f'{v:7.3g}' for v in k])} atm⁻ⁿ⁺¹·s⁻¹", ) if tunneling not in {"none", None}: @@ -625,11 +623,10 @@ def get_k( # noqa: PLR0913 # TODO(schneiderfelipe): when get_kappa accept deltas, this will # be probably unnecessary. logger.warning( - "assuming unitary tunneling coefficients due to incomplete " - "compound data", + "assuming unitary tunneling coefficients due to incomplete compound data", ) logger.info( - "(tunneling) reaction rate constants: " # noqa: G004 + "(tunneling) reaction rate constants: " f"{', '.join([f'{v:7.3g}' for v in k])} atm⁻ⁿ⁺¹·s⁻¹", ) @@ -652,7 +649,7 @@ def get_kappa( scheme: Scheme, compounds: dict, method: str = "eckart", - qrrho: bool = True, # noqa: FBT001, FBT002 + qrrho: bool = True, temperature: float = 298.15, ): r"""Obtain tunneling transmission coefficients at a given temperature. @@ -716,8 +713,8 @@ def get_kappa( >>> kappa * get_k(model.scheme, model.compounds, tunneling=None) array([8.e+10]) """ - scheme = rx.core._check_scheme(scheme) # noqa: SLF001 - compounds = rx.io._check_compounds(compounds) # noqa: SLF001 + scheme = rx.core._check_scheme(scheme) + compounds = rx.io._check_compounds(compounds) if method == "eckart": # NOTE(schneiderfelipe): We need electronic energies + ZPE here, so we @@ -751,7 +748,7 @@ def get_kappa( ) except RuntimeWarning as e: logger.warning( - f"using Wigner tunneling correction: {e}", # noqa: G004 + f"using Wigner tunneling correction: {e}", ) kappa = tunnel.wigner(vibfreq, temperature=temperature) elif method == "wigner": @@ -759,9 +756,8 @@ def get_kappa( elif method in {"none", None}: kappa = 1.0 else: - raise ValueError( # noqa: TRY003 - f"unavailable method: '{method}'", # noqa: EM102 - ) + msg = f"unavailable method: '{method}'" + raise ValueError(msg) kappas.append(kappa) @@ -769,19 +765,19 @@ def get_kappa( # somewhere else? vec_kappas = np.asarray(kappas).flatten() logger.info( - "(quantum) tunneling coefficients: " # noqa: G004 + "(quantum) tunneling coefficients: " f"{', '.join([f'{kappa:7.3g}' for kappa in vec_kappas])}", ) return vec_kappas -def get_drc( # noqa: PLR0913 +def get_drc( scheme, compounds, y0, t_span=None, method="RK23", - qrrho=True, # noqa: FBT002 + qrrho=True, scale="l mol-1 s-1", temperature=298.15, dx=1.5e3, # joules diff --git a/overreact/coords.py b/overreact/coords.py index ffa922b8..ceae7db2 100644 --- a/overreact/coords.py +++ b/overreact/coords.py @@ -1,9 +1,7 @@ -#!/usr/bin/env python3 # noqa: EXE001 - """Module dedicated to classifying molecules into point groups.""" -# TODO: add types to this module +# TODO(schneiderfelipe): add types to this module from __future__ import annotations __all__ = ["find_point_group", "symmetry_number"] @@ -25,10 +23,10 @@ # TODO(schneiderfelipe): alpha should depend on temperature? -def get_molecular_volume( # noqa: PLR0913 +def get_molecular_volume( atomnos, atomcoords, - full_output=False, # noqa: FBT002 + full_output=False, environment="water", method="garza", temperature=298.15, @@ -122,7 +120,7 @@ def get_molecular_volume( # noqa: PLR0913 >>> get_molecular_volume(data.atomnos, data.atomcoords, full_output=True, ... environment="benzene") (80., 593., 0.1) - """ # noqa: E501 + """ atomnos = np.atleast_1d(atomnos) _, _, atomcoords = inertia(np.ones_like(atomnos), atomcoords) vdw_radii = constants.vdw_radius(atomnos) @@ -136,7 +134,7 @@ def get_molecular_volume( # noqa: PLR0913 if full_output and method == "izato": cav_volumes = [] for _ in range(trials): - points = rx._misc.halton(n, 3) # noqa: SLF001 + points = rx._misc.halton(n, 3) points = v1 + points * (v2 - v1) tree = KDTree(points) @@ -156,16 +154,16 @@ def get_molecular_volume( # noqa: PLR0913 vdw_volume = np.mean(vdw_volumes) vdw_err = np.std(vdw_volumes) - logger.info(f"van der Waals volume = {vdw_volume} ± {vdw_err} ų") # noqa: G004 + logger.info(f"van der Waals volume = {vdw_volume} ± {vdw_err} ų") if full_output: if method == "izato": cav_volume = np.mean(cav_volumes) cav_err = np.std(cav_volumes) logger.debug( - f"Izato cavity volume = {cav_volume} ± {cav_err} ų", # noqa: G004 + f"Izato cavity volume = {cav_volume} ± {cav_err} ų", ) return (vdw_volume, cav_volume, max(vdw_err, cav_err)) - elif method == "garza": # noqa: RET505 + elif method == "garza": # TODO(schneiderfelipe): test for the following solvents: water, # pentane, hexane, heptane and octane. @@ -175,17 +173,18 @@ def get_molecular_volume( # noqa: PLR0913 temperature=temperature, pressure=pressure, ) - logger.debug(f"Garza cavity volume = {cav_volume} ų") # noqa: G004 + logger.debug(f"Garza cavity volume = {cav_volume} ų") return (vdw_volume, cav_volume, vdw_err) else: - raise ValueError(f"unavailable method: '{method}'") # noqa: EM102, TRY003 + msg = f"unavailable method: '{method}'" + raise ValueError(msg) return vdw_volume def _garza( vdw_volume, environment="water", - full_output=False, # noqa: FBT002 + full_output=False, temperature=298.15, pressure=constants.atm, ): @@ -239,7 +238,7 @@ def _garza( >>> _garza(100.0, full_output=True, environment="benzene") (665., 1.0, 1.07586575757374) """ - solvent = rx._misc._get_chemical(environment, temperature, pressure) # noqa: SLF001 + solvent = rx._misc._get_chemical(environment, temperature, pressure) # TODO(schneiderfelipe): things to do: # 1. check correctness of this function, @@ -252,22 +251,22 @@ def _garza( r_free = np.cbrt( solvent.Vm / (constants.angstrom**3 * constants.N_A) - solvent_volume, ) - r_M = np.cbrt(vdw_volume) # noqa: N806 + r_M = np.cbrt(vdw_volume) cav_volume = (r_M + r_free) ** 3 if not full_output: return cav_volume - r_S = np.cbrt(solvent_volume) # noqa: N806 + r_S = np.cbrt(solvent_volume) ratio = r_M / r_S area_free = r_free**2 - area_S_total = r_S**2 + area_free # noqa: N806 + area_S_total = r_S**2 + area_free x = max(area_free - r_M**2, 0.0) / area_S_total if np.isclose(x, 0.0): return cav_volume, 1.0, ratio - N_x = 4.0 * np.cbrt(cav_volume) ** 2 / area_S_total # noqa: N806 + N_x = 4.0 * np.cbrt(cav_volume) ** 2 / area_S_total return cav_volume, 1.0 + N_x * x / (1.0 - x), ratio @@ -374,11 +373,10 @@ def symmetry_number(point_group): elif pieces["letter"] == "s": symmetry_number = int(pieces["number"]) // 2 else: - raise ValueError( # noqa: TRY003 - f"unknown point group: '{point_group}'", # noqa: EM102 - ) + msg = f"unknown point group: '{point_group}'" + raise ValueError(msg) - logger.info(f"symmetry number = {symmetry_number}") # noqa: G004 + logger.info(f"symmetry number = {symmetry_number}") return symmetry_number @@ -424,7 +422,7 @@ def find_point_group(atommasses, atomcoords, proper_axes=None, rtol=0.0, atol=1. """ if len(atommasses) == 1: # atom point_group = "K" - elif len(atommasses) == 2: # diatomic molecule # noqa: PLR2004 + elif len(atommasses) == 2: # diatomic molecule point_group = "D∞h" if atommasses[0] == atommasses[1] else "C∞v" else: groups = _equivalent_atoms(atommasses, atomcoords) @@ -480,7 +478,7 @@ def find_point_group(atommasses, atomcoords, proper_axes=None, rtol=0.0, atol=1. atol=atol, ) - logger.info(f"point group = {point_group}") # noqa: G004 + logger.info(f"point group = {point_group}") return point_group @@ -493,11 +491,11 @@ def _find_point_group_linear(atomcoords, groups, rtol=0.0, atol=1.0e-2): """ if _has_inversion_center(atomcoords, groups, rtol=rtol, atol=atol): return "D∞h" - else: # noqa: RET505 + else: return "C∞v" -def _find_point_group_spheric( # noqa: PLR0913 +def _find_point_group_spheric( atomcoords, groups, axes, @@ -527,9 +525,9 @@ def _find_point_group_spheric( # noqa: PLR0913 ) for n, _ in proper_axes: - if n == 5: # noqa: PLR2004 + if n == 5: return "Ih" - elif n < 5: # noqa: PLR2004, RET505 + elif n < 5: break return "Oh" @@ -542,7 +540,7 @@ def _find_point_group_spheric( # noqa: PLR0913 # 1. doi:10.1016/0097-8485(76)80004-6 -def _find_point_group_asymmetric( # noqa: PLR0913 +def _find_point_group_asymmetric( atomcoords, groups, axes, @@ -578,7 +576,7 @@ def _find_point_group_asymmetric( # noqa: PLR0913 rtol=rtol, atol=atol, ) - elif rotor_class[1] in { # noqa: RET505 + elif rotor_class[1] in { "regular planar", "irregular planar", } or _get_mirror_planes( @@ -596,7 +594,7 @@ def _find_point_group_asymmetric( # noqa: PLR0913 return "C1" -def _find_point_group_symmetric( # noqa: PLR0913 +def _find_point_group_symmetric( atomcoords, groups, axes, @@ -624,7 +622,7 @@ def _find_point_group_symmetric( # noqa: PLR0913 count_twofold = 0 for n, _ in proper_axes: - if n == 2: # noqa: PLR2004 + if n == 2: count_twofold += 1 if n_principal == count_twofold: return _find_point_group_symmetric_dihedral( @@ -636,7 +634,7 @@ def _find_point_group_symmetric( # noqa: PLR0913 rtol=rtol, atol=atol, ) - if n < 2: # noqa: PLR2004 + if n < 2: break return _find_point_group_symmetric_nondihedral( atomcoords, @@ -652,7 +650,7 @@ def _find_point_group_symmetric( # noqa: PLR0913 # 1. doi:10.1016/0097-8485(76)80004-6 -def _find_point_group_symmetric_dihedral( # noqa: PLR0913 +def _find_point_group_symmetric_dihedral( atomcoords, groups, axes, @@ -689,15 +687,13 @@ def _find_point_group_symmetric_dihedral( # noqa: PLR0913 if mirror_axes: if mirror_axes[0][0] == "h": return f"D{proper_axes[0][0]}h" - elif ( # noqa: RET505 - len([v for c, v in mirror_axes if c == "v"]) == proper_axes[0][0] - ): # noqa: RET505, RUF100 + elif len([v for c, v in mirror_axes if c == "v"]) == proper_axes[0][0]: # all vertical mirror planes are dihedral for Dnd point groups return f"D{proper_axes[0][0]}d" return f"D{proper_axes[0][0]}" -def _find_point_group_symmetric_nondihedral( # noqa: PLR0913 +def _find_point_group_symmetric_nondihedral( atomcoords, groups, axes, @@ -734,9 +730,7 @@ def _find_point_group_symmetric_nondihedral( # noqa: PLR0913 if mirror_axes: if mirror_axes[0][0] == "h": return f"C{proper_axes[0][0]}h" - elif ( # noqa: RET505 - len([v for c, v in mirror_axes if c == "v"]) == proper_axes[0][0] - ): # noqa: RET505, RUF100 + elif len([v for c, v in mirror_axes if c == "v"]) == proper_axes[0][0]: return f"C{proper_axes[0][0]}v" improper_axes = _get_improper_axes( @@ -753,7 +747,7 @@ def _find_point_group_symmetric_nondihedral( # noqa: PLR0913 return f"C{proper_axes[0][0]}" -def _update_proper_axes( # noqa: PLR0913 +def _update_proper_axes( ax, axes, # found axes atomcoords, @@ -762,7 +756,7 @@ def _update_proper_axes( # noqa: PLR0913 rtol, atol, nondeg_axes=None, - normalize=False, # noqa: FBT002 + normalize=False, ): """Update axes with ax, and return it with added order (or None). @@ -798,7 +792,7 @@ def _update_proper_axes( # noqa: PLR0913 return axes, None -def _get_proper_axes( # noqa: C901, PLR0912, PLR0913 +def _get_proper_axes( atomcoords, groups, axes, @@ -917,7 +911,7 @@ def _get_proper_axes( # noqa: C901, PLR0912, PLR0913 nondeg_axes=nondeg_axes, normalize=True, ) - if rotor_class[0] == "spheric" and order == 5: # noqa: PLR2004 + if rotor_class[0] == "spheric" and order == 5: return sorted(found_axes, reverse=True) for b in group[:i]: @@ -933,11 +927,11 @@ def _get_proper_axes( # noqa: C901, PLR0912, PLR0913 nondeg_axes=nondeg_axes, normalize=True, ) - if rotor_class[0] == "spheric" and order == 5: # noqa: PLR2004 + if rotor_class[0] == "spheric" and order == 5: return sorted(found_axes, reverse=True) if rotor_class[0] == "spheric": - twofold_axes = [ax for o, ax in found_axes if o == 2] # noqa: PLR2004 + twofold_axes = [ax for o, ax in found_axes if o == 2] for i, ax_a in enumerate(twofold_axes): for ax_b in twofold_axes[:i]: ax = np.cross(ax_a, ax_b) @@ -952,7 +946,7 @@ def _get_proper_axes( # noqa: C901, PLR0912, PLR0913 nondeg_axes=nondeg_axes, normalize=True, ) - if order == 5: # noqa: PLR2004 + if order == 5: return sorted(found_axes, reverse=True) return sorted(found_axes, reverse=True) @@ -997,7 +991,7 @@ def _guess_orders(groups, rotor_class): return range(2, max_order + 1) -def _update_improper_axes( # noqa: PLR0913 +def _update_improper_axes( n, ax, axes, @@ -1005,7 +999,7 @@ def _update_improper_axes( # noqa: PLR0913 groups, rtol, atol, - normalize=False, # found axes # noqa: FBT002 + normalize=False, # found axes ): """Update axes with ax and return it. @@ -1033,7 +1027,7 @@ def _update_improper_axes( # noqa: PLR0913 return axes -def _get_improper_axes( # noqa: PLR0913 +def _get_improper_axes( atomcoords, groups, axes, @@ -1113,7 +1107,7 @@ def _get_improper_axes( # noqa: PLR0913 return sorted(found_axes, reverse=True) -def _update_mirror_axes( # noqa: PLR0913 +def _update_mirror_axes( ax, axes, # found axes atomcoords, @@ -1122,7 +1116,7 @@ def _update_mirror_axes( # noqa: PLR0913 atol, proper_axes, nondeg_axes=None, - normalize=False, # noqa: FBT002 + normalize=False, ): """Update axes with ax and return it. @@ -1169,7 +1163,7 @@ def _update_mirror_axes( # noqa: PLR0913 return axes -def _get_mirror_planes( # noqa: C901, PLR0913 +def _get_mirror_planes( atomcoords, groups, axes, @@ -1254,7 +1248,7 @@ def _kf(x): c, v = x if c: return -ord(c), v - else: # noqa: RET505 + else: return 0, v found_axes = [] @@ -1489,24 +1483,25 @@ def _operation(name, order=2, axis=None): if name == "e": return np.eye(3) - if name in {"c", "σ", "sigma", "s"}: # noqa: RUF001 + if name in {"c", "σ", "sigma", "s"}: # normalize axis axis = np.asarray(axis) axis = axis / np.linalg.norm(axis) if name in {"c", "s"}: rotation = Rotation.from_rotvec(2.0 * np.pi * axis / order).as_matrix() - if name in {"σ", "sigma", "s"}: # noqa: RUF001 + if name in {"σ", "sigma", "s"}: reflection = np.eye(3) - 2.0 * np.outer(axis, axis) if name == "c": return rotation - if name in {"σ", "sigma"}: # noqa: RUF001 + if name in {"σ", "sigma"}: return reflection if name == "s": return rotation @ reflection - raise ValueError(f"unknown operation: '{name}'") # noqa: EM102, TRY003 + msg = f"unknown operation: '{name}'" + raise ValueError(msg) def _classify_rotor(moments, rtol=0.0, atol=1.0e-2, slack=0.870): @@ -1683,13 +1678,14 @@ def gyradius(atommasses, atomcoords, method="iupac"): return np.sqrt( np.average(np.diag(atomcoords @ atomcoords.T), weights=atommasses), ) - elif method == "mean": # noqa: RET505 + elif method == "mean": return np.sqrt(np.mean(np.diag(atomcoords @ atomcoords.T))) else: - raise ValueError(f"unavailable method: '{method}'") # noqa: EM102, TRY003 + msg = f"unavailable method: '{method}'" + raise ValueError(msg) -def inertia(atommasses, atomcoords, align=True): # noqa: FBT002 +def inertia(atommasses, atomcoords, align=True): r"""Calculate primary moments and axes from the inertia tensor. Parameters @@ -1761,7 +1757,7 @@ def inertia(atommasses, atomcoords, align=True): # noqa: FBT002 moments, axes = np.linalg.eigh(inertia_tensor) if align: return inertia(atommasses, atomcoords @ axes, align=False) - logger.debug(f"moments = {moments} amu·Å²") # noqa: G004 + logger.debug(f"moments = {moments} amu·Å²") return moments, axes, atomcoords @@ -1824,15 +1820,15 @@ def calc_hessian(atommasses, atomcoords, vibfreqs, vibdisps): -0.01313383, -0.00678954, 0.25045664, 0.29443748]]) """ dof = 3 * len(atommasses) - L_cart = np.asarray(vibdisps).reshape((len(vibfreqs), dof)).T # noqa: N806 + L_cart = np.asarray(vibdisps).reshape((len(vibfreqs), dof)).T # this function is correct until here - L_cart = np.linalg.qr(L_cart, mode="complete")[0] # noqa: N806 + L_cart = np.linalg.qr(L_cart, mode="complete")[0] atommasses_sqrt = np.sqrt([mass for mass in atommasses for _ in range(3)]) - D = eckart_transform(atommasses, atomcoords) # noqa: N806 - M = np.diag(1.0 / atommasses_sqrt) # noqa: N806 - L = np.linalg.solve(M @ D, L_cart) # noqa: N806 + D = eckart_transform(atommasses, atomcoords) + M = np.diag(1.0 / atommasses_sqrt) + L = np.linalg.solve(M @ D, L_cart) assert np.allclose(M @ D @ L, L_cart), "L_cart is not orthogonal" @@ -1986,28 +1982,28 @@ def eckart_transform(atommasses, atomcoords): y *= np.sqrt(atommasses[:, np.newaxis]) z *= np.sqrt(atommasses[:, np.newaxis]) - D_trans = np.block( # noqa: N806 + D_trans = np.block( [x.reshape(1, dof).T, y.reshape(1, dof).T, z.reshape(1, dof).T], - ) # noqa: RUF100 - D_rot = np.array( # noqa: N806 + ) + D_rot = np.array( [ np.cross((atomcoords @ axes)[i], axes[:, j]) / np.sqrt(atommasses[i]) for i in range(natom) for j in range(3) ], ) - D = np.block([D_trans, D_rot]) # noqa: N806 + D = np.block([D_trans, D_rot]) return np.linalg.qr(D, mode="complete")[0] # NOTE(schneiderfelipe): thresh was found to be reasonable # when greater than or equal to 0.106. -def _equivalent_atoms( # noqa: C901 +def _equivalent_atoms( atommasses, atomcoords, method="cluster", thresh=0.106, - plot=False, # noqa: FBT002 + plot=False, ): """Generate groups of symmetry equivalent atoms. @@ -2069,7 +2065,7 @@ def _equivalent_atoms( # noqa: C901 """ if len(atommasses) == 1: # atom return [[0]] - elif len(atommasses) == 2: # diatomic molecule # noqa: PLR2004, RET505 + elif len(atommasses) == 2: # diatomic molecule if atommasses[0] == atommasses[1]: return [[0, 1]] return [[0], [1]] @@ -2082,14 +2078,14 @@ def _update_groups_with_condition(condition, groups): return groups if method == "cluster": - D = squareform(pdist(atomcoords)) # noqa: N806 + D = squareform(pdist(atomcoords)) omega = np.mean(D, axis=0) sigma = np.std(D, axis=0) delta = np.sqrt(np.sum(D**2, axis=0)) criteria = np.block([[omega], [sigma], [delta]]).T - Z = linkage(pdist(criteria), method="single") # noqa: N806 + Z = linkage(pdist(criteria), method="single") clusters = fcluster(Z, thresh, criterion="distance") # TODO(schneiderfelipe): this was for debug and should eventually be removed. @@ -2119,6 +2115,7 @@ def _update_groups_with_condition(condition, groups): for mass in np.unique(atommasses): groups = _update_groups_with_condition(atommasses == mass, groups) else: - raise ValueError(f"unavailable method: '{method}'") # noqa: EM102, TRY003 + msg = f"unavailable method: '{method}'" + raise ValueError(msg) return sorted(groups, key=len) diff --git a/overreact/core.py b/overreact/core.py index 55aaf025..8a75fc4c 100644 --- a/overreact/core.py +++ b/overreact/core.py @@ -1,14 +1,14 @@ -#!/usr/bin/env python3 # noqa: EXE001 - """Module dedicated to parsing and modeling of chemical reaction networks.""" +from __future__ import annotations + __all__ = ["Scheme", "parse_reactions"] import itertools import re -from typing import NamedTuple, Sequence, Union +from typing import NamedTuple import numpy as np @@ -23,16 +23,16 @@ class Scheme(NamedTuple): See `overreact.io.parse_model`. """ - compounds: Sequence[str] + compounds: list[str] """A descriptor of compounds.""" - reactions: Sequence[str] + reactions: list[str] """A descriptor of reactions.""" - is_half_equilibrium: Sequence[bool] + is_half_equilibrium: list[bool] """An indicator of whether a reaction is half-equilibrium.""" A: np.ndarray """A matrix of stoichiometric coefficients between reactants and products.""" B: np.ndarray - """A matrix of stoichiometric coefficients between reactants and transition states.""" # noqa: E501 + """A matrix of stoichiometric coefficients between reactants and transition states.""" _abbr_environment = { @@ -50,7 +50,7 @@ class Scheme(NamedTuple): } -def _check_scheme(scheme_or_text: Union[Scheme, str]) -> Scheme: +def _check_scheme(scheme_or_text: Scheme | str) -> Scheme: """Interface transparently between strings and schemes. Parameters @@ -81,7 +81,7 @@ def _check_scheme(scheme_or_text: Union[Scheme, str]) -> Scheme: return parse_reactions(scheme_or_text) -def get_transition_states(A, B, is_half_equilibrium): # noqa: N803 +def get_transition_states(A, B, is_half_equilibrium): """Return the indices of transition states for each reaction. Parameters @@ -399,7 +399,7 @@ def is_transition_state(name): return any(marker in name for marker in {"‡", "#"}) -def parse_reactions(text: Union[str, Sequence[str]]) -> Scheme: # noqa: C901 +def parse_reactions(text: str | list[str]) -> Scheme: """ Parse a kinetic model as a chemical reaction scheme. @@ -598,8 +598,8 @@ def parse_reactions(text: Union[str, Sequence[str]]) -> Scheme: # noqa: C901 tuple[str, str, bool, str], tuple[tuple[tuple[int, str], ...], bool], ] = {} - A = [] # coefficients between reactants and products # noqa: N806 - B = [] # coefficients between reactants and transition states # noqa: N806 + A = [] # coefficients between reactants and products + B = [] # coefficients between reactants and transition states def _add_reaction(reactants, products, is_half_equilibrium, transition): """Local helper function with side-effects.""" @@ -614,13 +614,13 @@ def _add_reaction(reactants, products, is_half_equilibrium, transition): # found new reaction reactions[(reactants, products, is_half_equilibrium, transition)] = None - A_vector = np.zeros(len(compounds)) # noqa: N806 + A_vector = np.zeros(len(compounds)) for coefficient, reactant in reactants: A_vector[compounds[reactant]] = -coefficient - B_vector = A_vector # noqa: N806 + B_vector = A_vector if transition is not None: - B_vector = A_vector.copy() # noqa: N806 + B_vector = A_vector.copy() # it's assumed that # 1. there's a singe transition compound, and @@ -634,7 +634,7 @@ def _add_reaction(reactants, products, is_half_equilibrium, transition): is_half_equilibrium and (products, reactants, is_half_equilibrium, transition) in reactions ): - B_vector = np.zeros(len(compounds)) # noqa: N806 + B_vector = np.zeros(len(compounds)) A.append(A_vector) B.append(B_vector) @@ -676,7 +676,7 @@ def _add_reaction(reactants, products, is_half_equilibrium, transition): else: after_transitions[reactants] = [products] continue - elif is_transition_state(products[-1][-1]): # noqa: RET507 + elif is_transition_state(products[-1][-1]): for after_products in after_transitions.get(products, []): _add_reaction(reactants, after_products, is_half_equilibrium, products) @@ -691,15 +691,15 @@ def _add_reaction(reactants, products, is_half_equilibrium, transition): return Scheme( compounds=tuple(compounds), reactions=tuple(_unparse_reactions(reactions)), - is_half_equilibrium=rx._misc.totuple( # noqa: SLF001 + is_half_equilibrium=rx._misc.totuple( [reaction[2] for reaction in reactions], ), - A=rx._misc.totuple( # noqa: SLF001 + A=rx._misc.totuple( np.block( [[vector, np.zeros(len(compounds) - len(vector))] for vector in A], ).T, ), - B=rx._misc.totuple( # noqa: SLF001 + B=rx._misc.totuple( np.block( [[vector, np.zeros(len(compounds) - len(vector))] for vector in B], ).T, @@ -752,7 +752,7 @@ def _parse_reactions(text): except AttributeError: lines = text for line in lines: - line = line.split("//")[0].strip() # noqa: PLW2901 + line = line.split("//")[0].strip() if not line: continue @@ -763,9 +763,9 @@ def _parse_reactions(text): pieces[2::2], ): if arrow == "<-": - reactants, products, arrow = products, reactants, "->" # noqa: PLW2901 - reactants = tuple(_parse_side(reactants)) # noqa: PLW2901 - products = tuple(_parse_side(products)) # noqa: PLW2901 + reactants, products, arrow = products, reactants, "->" + reactants = tuple(_parse_side(reactants)) + products = tuple(_parse_side(products)) if arrow == "<=>": yield reactants, products, True @@ -833,7 +833,7 @@ def _parse_side(side): """ for token in re.split(r"\s+\+\s+", side): - token = re.match( # noqa: PLW2901 + token = re.match( r"\s*(?P\d+)?\s*(?P[^\s]+)\s*", token, ).groupdict(1)