diff --git a/src/autochem/rate/data.py b/src/autochem/rate/data.py index 1d954024..0b316ed5 100644 --- a/src/autochem/rate/data.py +++ b/src/autochem/rate/data.py @@ -648,7 +648,10 @@ def fit( if not validate or bad_fit == "fill": try: with np.errstate(all="raise"): - obj(T=T) + vals = obj(T=T) + if not np.all(np.isfinite(vals)): + msg = "Fitted rate gives non-finite values over the input temperature range." + raise FloatingPointError(msg) except FloatingPointError as e: if A_fill is not None: return cls(order=order, A=A_fill, b=0, E=0) @@ -917,22 +920,25 @@ def fit( # noqa: PLR0913 (Options: "fill", replace with the fill value, or numpy.seterr options) :return: Rate fit """ - bad_fits = [ - "fill" if np.any(np.isclose(p, bad_fit_fill_pressures)) else bad_fit - for p in Ps - ] - k_data_fits = [ - ArrheniusRateFit.fit( + k_fits = [] + for k_data, P in zip(np.transpose(k_data), Ps, strict=True): + if validate: + print(f"Fitting Arrhenius rate for {P = }") + + bad_fit = ( + "fill" if np.any(np.isclose(P, bad_fit_fill_pressures)) else bad_fit + ) + k_fit = ArrheniusRateFit.fit( Ts=Ts, - ks=ks, + ks=k_data, A_fill=A_fill, - bad_fit=f, # pyright: ignore[reportArgumentType] + bad_fit=bad_fit, validate=validate, order=order, units=units, ) - for ks, f in zip(np.transpose(k_data), bad_fits, strict=True) - ] + k_fits.append(k_fit) + k_high_fit = None if k_high is not None: k_high_fit = ArrheniusRateFit.fit( @@ -943,9 +949,9 @@ def fit( # noqa: PLR0913 return cls( order=order, - As=[f.A for f in k_data_fits], - bs=[f.b for f in k_data_fits], - Es=[f.E for f in k_data_fits], + As=[f.A for f in k_fits], + bs=[f.b for f in k_fits], + Es=[f.E for f in k_fits], Ps=Ps, # pyright: ignore[reportArgumentType] ) diff --git a/src/autochem/therm/_species.py b/src/autochem/therm/_species.py index 0264a2b6..513ee426 100644 --- a/src/autochem/therm/_species.py +++ b/src/autochem/therm/_species.py @@ -1,8 +1,9 @@ """Thermodynamic data.""" import datetime +import re from collections.abc import Sequence -from typing import ClassVar, Literal +from typing import ClassVar, Literal, Self import altair as alt import numpy as np @@ -23,6 +24,17 @@ class Species(Scalable): name: str therm: Therm_ + def racemize( + self, enant_suffixes: tuple[str, str] = ("0", "1"), rac_suffix: str = "R" + ) -> Self: + """Racemize the thermodynamic functions for a chiral species.""" + suffix0, suffix1 = enant_suffixes + pattern = re.compile(rf"(?:{re.escape(suffix0)}|{re.escape(suffix1)})$") + if pattern.search(self.name): + self.name = pattern.sub(rac_suffix, self.name) + self.therm = self.therm.racemize() + return self + # Private attributes _scalers: ClassVar[Scalers] = {"therm": (lambda c, x: c * x)} diff --git a/src/autochem/therm/data.py b/src/autochem/therm/data.py index e12f6710..48b5d6f4 100644 --- a/src/autochem/therm/data.py +++ b/src/autochem/therm/data.py @@ -45,6 +45,11 @@ class BaseTherm(ThermCalculator, UnitManager, Frozen, Scalable, SubclassTyped, a formula: Formula_ charge: int = 0 + @abc.abstractmethod + def racemize(self) -> Self: + """Racemize the thermodynamic functions for a chiral species.""" + return self + def display( # noqa: PLR0913 self, props: Sequence[Literal["Cv", "Cp", "S", "H", "dH"]] = ("Cp", "S", "H"), @@ -227,6 +232,11 @@ def sort_by_temperature(self) -> Self: self.model_config["frozen"] = frozen return self + def racemize(self) -> Self: + """Racemize the thermodynamic functions for a chiral species.""" + self.Z0 = np.add(self.Z0, np.log(2)).tolist() + return self + # Replace this with a model validator (before or after), allowing extra arguments def __init__( self, @@ -530,6 +540,12 @@ def piecewise_conditions( calc_high.in_bounds(T, include_max=True), ] + def racemize(self) -> Self: + """Racemize the thermodynamic functions for a chiral species.""" + self.coeffs_low[-1] += np.log(2).item() + self.coeffs_high[-1] += np.log(2).item() + return self + def heat_capacity( self, T: ArrayLike, # noqa: N803 diff --git a/src/autochem/util/plot.py b/src/autochem/util/plot.py index 93f94dc8..7a19b344 100644 --- a/src/autochem/util/plot.py +++ b/src/autochem/util/plot.py @@ -166,7 +166,7 @@ def regular_scale(val_range: tuple[float, float]) -> alt.Scale: :param val_range: Range :return: Scale """ - return alt.Scale(domain=val_range) + return alt.Scale(domain=val_range, domainMax=val_range[-1], nice=True) def regular_scale_axis(val_range: tuple[float, float]) -> alt.Axis: diff --git a/src/automol/geom/_1conv.py b/src/automol/geom/_1conv.py index 711ec6ab..3ae11c4b 100644 --- a/src/automol/geom/_1conv.py +++ b/src/automol/geom/_1conv.py @@ -9,6 +9,7 @@ import pyparsing as pp from numpy.typing import ArrayLike from pyparsing import pyparsing_common as ppc +from rdkit.Chem.inchi import MolBlockToInchiAndAuxInfo from phydat import phycon from .. import vmat @@ -340,38 +341,10 @@ def inchi_with_numbers(geo, stereo=True, gra=None): gra = graph_base.set_stereo_from_geometry(gra, geo) mlf, key_map_inv = molfile_with_atom_mapping(gra, geo=geo) - rdm = rdkit_.from_molfile(mlf) - ich, aux_info = rdkit_.to_inchi(rdm, with_aux_info=True) + ich, aux_info = MolBlockToInchiAndAuxInfo(mlf) nums_lst = _parse_sort_order_from_aux_info(aux_info) nums_lst = tuple(tuple(map(key_map_inv.__getitem__, nums)) for nums in nums_lst) - - # Assuming the MolFile InChI works, the above code is all we need. What - # follows is to correct cases where it fails. - # This only appears to work sometimes, so when it doesn't, we fall back on - # the original inchi output. - if geo is not None: - gra = graph_base.set_stereo_from_geometry(gra, geo) - gra = graph_base.implicit(gra) - sub_ichs = inchi_base.split(ich) - - failed = False - - new_sub_ichs = [] - for sub_ich, nums in zip(sub_ichs, nums_lst): - sub_gra = graph_base.subgraph(gra, nums, stereo=True) - sub_ich = _connected_inchi_with_graph_stereo(sub_ich, sub_gra, nums) - if sub_ich is None: - failed = True - break - - new_sub_ichs.append(sub_ich) - - # If it worked, replace the InChI with our forced-stereo InChI. - if not failed: - ich = inchi_base.join(new_sub_ichs) - ich = inchi_base.standard_form(ich) - return ich, nums_lst