"""
Created on 15 Feb 2021
@author: Christian Rab
"""
from abc import ABC, abstractmethod
import collections
import decimal
import math
import os
import re
from timeit import default_timer as timer
import warnings
import prodimopy.utils as putils
from prodimopy.chemistry.elements_data import elements as elemdata
import csv
import numpy as np
from numpy.typing import NDArray
from tqdm import tqdm
import copy
import logging
logger = logging.getLogger(__name__)
_non_species = [
"PE",
"PHOT",
"PHOTON",
"Photon",
"XPHOT",
"M",
"ELECTR",
"PHOTON",
"GRAIN",
"FORM",
"CRPHOT",
"CRP",
"CR",
"SPONT",
"DESCR",
"DESPH",
"DESTH",
"e-",
"dust",
"RADASS",
"PUMP",
"XH", # some strange species in KIDA2024
]
[docs]
class Reaction(object):
"""
General representation of one chemical reaction.
"""
def __init__(self):
self.id: int = 0
""" : The Reaction Id as is. """
self.idU: int = 0
""" : The Umist reaction id ... it is unclear what this actually is. """
self.type: str = ""
""" : The Reaction type identifier as used in |prodimo| """
self.gasphase: bool = None
""" : Gas phase reaction or not. """
self.reactants: list[str] = list()
""" : The list of reaction reactants (Species names) """
self.products: list[str] = list()
""" : The list of reaction products (Species names) """
self.coeffs: list[NDArray[np.float64]] = list()
""" : The list of `([alpha,beta,gamma])` coefficients per temperature range. """
self.temps: list[NDArray[np.float64]] = list()
""" : The list of temperature ranges given by `(min,max)`. """
self.rate: float | None = None
""" : The real rate for a given set of physical conditions. Only exists if the network is loaded from an actual model. """
self.method: str | None = None
""" : And indicator on how the rated was derived. Usually UMIST method identifier (M: measured; C: calculated; E: estimated; L: literature). """
self.accuracy: str | None = None
""" : A string describing the quality of the reactions coefficients. """
self.comment: str = ""
""" : Any relevant comment for this particular reaction. """
# this is just for internal use to have coeff also as a decimal, to take into account the strange
# format of Reactions.in ... will only be filled in that case
self._coeff2 = None
[docs]
def arrhenius_rate(self) -> tuple[list[NDArray[np.float64]], list[NDArray[np.float64]]]:
"""
Calculate the Arrhenius rate for the given (from the data) temperature range of the reaction.
Returns
-------
tuple[list[NDArray[np.float64]], list[NDArray[np.float64]]] :
A tuple containing two lists of numpy arrays. The first list contains the temperature points for each temperature range,
the second list contains the corresponding Arrhenius rates for each temperature range.
"""
rates = list()
temps = list()
for coeff, temp in zip(self.coeffs, self.temps):
temppoints = np.logspace(np.log10(temp[0]), np.log10(temp[1]), 30, endpoint=True)
rate = coeff[0] * (temppoints / 300.0) ** coeff[1] * np.exp(-coeff[2] / temppoints)
temps.append(temppoints)
rates.append(rate)
return temps, rates
[docs]
def compare(self, reaction: "Reaction", specsonly: bool = False) -> dict[str, bool]:
"""
Compares the Reaction to the passed reaction.
The routines evaluates for a selection of the properties from `class:`Reaction` individually if it
is equal or not. See below for the return value.^
Parameters
----------
reaction
Another reaction to compare with.
specsonly
If True, only the reactants and products are compared. Default: `False`
In that case the returned dictionary only contains the keys `reactants` and `products`.
Returns
-------
dict[str, bool] :
A dictionary with boolean values with the following keys. The keys have the
same names as the properties of the :class:`Reaction` class.
* `type` is the reaction type the sames
* `gasphase` is the gasphase flag the same
* `reactants` does have the reaction the same reactants (order does not matter)
* `products` same as reactants but for the products
* `coeffs` are the coefficients (alpha,beta,gamma) the same
* `temps` are the temperature ranges the sames
..todo::
- the routine is quite slow if speconly=False
- deal with comparison of the network if actual rates only exist in one network
- make the tolerances for comparison of the rate a parameter somehow
"""
equal = dict()
equal["reactants"] = sorted(self.reactants) == sorted(reaction.reactants)
equal["products"] = sorted(self.products) == sorted(reaction.products)
if specsonly:
return equal
equal["type"] = self.type == reaction.type
equal["gasphase"] = self.gasphase == reaction.gasphase
equal["coeffs"] = True
for coff, cofNew in zip(self.coeffs, reaction.coeffs):
equal["coeffs"] = equal["coeffs"] and (coff == cofNew).all()
equal["temps"] = True
for temp, tempNew in zip(self.temps, reaction.temps):
equal["temps"] = equal["temps"] and (temp == tempNew).all()
# TODO: if not loaded in both reactions, the comparison does not make sense
# Maybe write an error.
if self.rate is not None and reaction.rate is not None:
# print(self.rate,reaction.rate,math.isclose(self.rate,reaction.rate,rel_tol=1e-20,abs_tol=0.0))
equal["rate"] = math.isclose(self.rate, reaction.rate, rel_tol=1e-4, abs_tol=0.0)
else:
equal["rate"] = True
return equal
def __str__(self):
out = (
f"{str(self.id):>4s} {str(self.idU):>5s} {self.type:3s} "
+ "{:32s}".format("{}".format(self.reactants))
+ " --> "
+ "{:32s}".format("{}".format(self.products))
)
for i, (coeff, temp) in enumerate(zip(self.coeffs, self.temps)):
coeffstr = f"{np.array2string(coeff, formatter={'float_kind': '{:+7.2e}'.format})} "
tempstr = f"{np.array2string(temp, formatter={'float_kind': '{:7.2e}'.format})} "
if i == 0:
out += coeffstr + tempstr
out += " " + self.comment
else:
out += "\n" + " " * 84 + coeffstr + tempstr
# Don't do the rate thing, is confusing as it is for one particular point
# if self.rate!=None:
# out=out+"{:13.5e}".format(self.rate)
return out
[docs]
class ReactionNetwork(ABC):
"""
General representation of a chemical reaction network.
.. todo::
- make it an abstract class that requires the child classes to implement the load_reactions routine.
- properly define what is an ice species, e.g. add routines such as `species_is_ice` which can be overwritten.
- properly define non-species, also allow to overwrite that in subclass e.g. routine species_is_nonspecies
"""
def __init__(
self, name: str, filename: str | None = None, filename_dictionary: str | None = None
):
"""
Initialize a new chemical reaction network.
Parameters
----------
name :
A name for the reaction network.
filename :
The name/path of the file containing the Reaction network. Can be `None` if the network is created from scratch.
Or if the data is read later on. If it is not `None` the network is loaded from the file in the constructor.
Default: `None`
filename_dictionary :
The name/path of the dictionary file containing the mapping of the names used in this reaction network, to the
new names. Can be `None` if no mapping is used. If it is not `None` the mapping is applied to the network in the constructor.
after reading the network. Default: `None`
Attributes
----------
"""
self.name: str = name
""" : A name for the reaction network """
self.filename: str | None = filename
""" : The name/path of the file containing the reaction data. """
self.filename_dictionary: str | None = filename_dictionary
""" : The name/path of the dictionary file containing the mapping of the names used in this reaction network. """
self.reactions: list[Reaction] = []
""" : The list of reactions included in the network. """
self.non_species: list[str] = _non_species
""" : A list of Species names that are considered as non species (not really a chemical species) but might be included in the network. """
self._species: list[str] = []
""" : Caches the list of species derived from the list of Reactions. Is filled on demand. """
if filename is not None:
print(f"Loading reactions for network {self.name} from file {filename} ...")
self.load_reactions(filename)
if filename_dictionary is not None:
print(
f"Mapping species names in network {self.name} using dictionary {filename_dictionary} ..."
)
self.map_species_names(filename_dictionary)
@property
def species(self) -> list[str]:
"""
The list of individual species in the network. The list is construction from the reactions data.
"""
# only load if not yet loaded.
if len(self._species) == 0:
# Create a species list
for reac in self.reactions:
for spec in reac.reactants + reac.products:
if spec not in self._species and spec not in self.non_species:
self._species.append(spec)
self._species.sort()
return self._species
@property
def elements(self) -> list[str]:
"""
The list of unique elements in the network, sorted by atomic number
.. todo::
- use :func:`parse_stoichiometry` to be consistent in the parsing of the species names.
"""
elset = set()
for spec in self.species:
# strip ice/ion suffixes, then extract capitalised tokens
formula = spec.rstrip("+-#")
for token in re.findall(r"[A-Z][a-z]*", formula):
if token in elemdata._by_symbol:
elset.add(token)
return sorted(elset, key=lambda s: elemdata[s].atomic_number)
def __str__(self):
return (
"Name: "
+ self.name
+ "\n"
+ " File : "
+ str(self.filename)
+ "\n"
+ " Number of Reactions: "
+ str(len(self.reactions))
+ "\n Number of Species: "
+ str(len(self.species))
+ "\n Number of Elements: "
+ str(len(self.elements))
)
[docs]
def duplicates(self) -> list[tuple[Reaction, Reaction]]:
"""
Returns possible duplicate reactions in the network.
Currently a two reactions are a potential duplicate if they have the same reactants and products.
Returns
-------
list[tuple[:class:`Reaction`, :class:`Reaction`]] :
Returns a list of tuples containing the duplicate Reactions or an empty list if there are
no duplicates
"""
seen = set()
dups = list()
with tqdm(total=len(self.reactions), desc="Checking for duplicates ...") as pbar:
for reac in self.reactions:
# first we check if the reaction is already in the duplicate list
indup = False
for i, dup in enumerate(dups):
if self._duplicatereaction(dup[0].compare(reac, specsonly=True)):
dup.append(reac)
dups[i] = dup
indup = True
break
if not indup:
for reacseen in seen:
if self._duplicatereaction(reac.compare(reacseen, specsonly=True)):
dups.append([reac, reacseen])
break
seen.add(reac)
pbar.update(1)
for dup in dups:
dup.sort(key=lambda r: r.id)
return dups
[docs]
def rename_species(self, old: str, new: str, variants: tuple[str] | None = ("")) -> int:
"""
Renames species with name `old` with species name `new` in the whole network.
Parameters
----------
old :
The species name of the species that should be renamed.
new :
The new name of the species.
variants :
Postfix to the species names that should also be considered for renaming.
For example if `variants=["","+","++","#"]` and the old name is `SP` and the new name is `SPNEW` the routine would rename
`SP,SP+,SP#` with `SPNEW,SPNEW+,SPNEW#` , respectively.
If `None` -> `variants=("")`, which means only identical names are considered,. Default: `("")`
log : boolean
log detailed changes on the screen. Default: False
Returns
-------
int :
The total number of changed Reactions.
"""
if variants is None or len(variants) == 0:
variants = ("",)
# number of changes done
nchanged = 0
for reac in self.reactions:
for variant in variants:
spec = old + variant
nspecex = new + variant
# indices of the products that match the species
ip = [i for i, e in enumerate(reac.products) if e == spec]
ir = [i for i, e in enumerate(reac.reactants) if e == spec]
if len(ip) > 0 or len(ir) > 0:
logger.debug(f"Replace {spec} with {nspecex} in reaction {reac.id}")
logger.debug(reac)
nchanged += 1
for i in ip:
reac.products[i] = nspecex
for i in ir:
reac.reactants[i] = nspecex
# reset internal species so that is is reloaded again in method property `species``
self._species = list()
outstr = f"Renamed {old} with {new} including variants {variants} in {nchanged} reactions."
# print(outstr)
logger.info(outstr)
return nchanged
[docs]
def map_species_names(self, filename_dictionary: str):
"""
Maps the species names in the reaction network to new names based on a provided dictionary file.
Note that only the variants ("","#") are considered (see :func:`~prodimopy.chemistry.network.rename_species`).
Parameters
----------
filename_dictionary :
The path to the dictionary file containing the mapping of the names used in this reaction network, to the
new names. E.g. Here KIDA to UMIST.
see also :func:`~prodimopy.chemistry.network.read_species_dictionary`.
"""
logger.debug(f"map_species_names: {filename_dictionary}")
self.filename_dictionary = filename_dictionary
mapping = read_species_dictionary(filename_dictionary)
for spec in mapping:
self.rename_species(spec, mapping[spec], variants=("", "#"))
[docs]
def compareSpecies(self, reactionNetwork):
"""
Compares the species list of two networks.
Parameters
----------
reactionNetwork : :class:`ReactionNetwork`
Another reaction network
Returns
-------
(tuple) :
* `boolean` ... equal or not
* list ... list of species only in old network
* list ... list of species only in new network
"""
specOld = self.species
specNew = reactionNetwork.species
# probably not very fast but fast enough I guess
onlyOld = list()
for spec in specOld:
if spec not in specNew:
onlyOld.append(spec)
onlyNew = list()
for spec in specNew:
if spec not in specOld:
onlyNew.append(spec)
# if both empty the species are the same,
return (not onlyOld and not onlyNew, onlyOld, onlyNew)
[docs]
def find_species(
self,
searchstr: str | list[str],
checkelements: bool = False,
ignore_non_species: bool = False,
) -> list[str]:
"""
Find all species of the network that contain the given search string.
Parameters
----------
searchstr :
The string or list of strings to search for in the species names.
checkelements :
If `True` the search string is interpreted as an element symbol and the routine checks
additionally the stoichiometry for that element in the species. Default: `False`
ignore_non_species :
If `True` the routine ignores species that are in the `non_species` list of the network. Default: `False`
Returns
-------
list[str]
A list of species names that contain the search string.
"""
if isinstance(searchstr, str):
searchstr = [searchstr]
if checkelements:
found_species = [
spec
for spec in self.species
if any(s in spec and s in parse_stoichiometry(spec) for s in searchstr)
]
else:
found_species = [spec for spec in self.species if any(s in spec for s in searchstr)]
if ignore_non_species:
found_species = [spec for spec in found_species if spec not in self.non_species]
return found_species
[docs]
def find_isomers(self, species: str | None = None) -> list[tuple[str]] | tuple[str]:
"""
Find all species that have isomers (i.e. the same stoichiometry) in the network.
If a species name is passed, the routine only returns the isomer group for this
particular species. The passed species does not have to be in the network.
If the species is the only isomer, an empty tuple is returned. If the species is not in the network, it is not included in the returned isomer group.
By default the routine searches for all species in the network for isomers, and returns
a list of tuples containing those isomer groups.
Ice species '#' are not considered as isomers, and are ignored in the search.
Parameters
----------
network :
The reaction network to analyze.
species :
The specific species name to find isomers for, by default `None` -> search for all isomer groups in the network
Returns
-------
list[tuple[str]] or tuple[str]
A list of tuples where each tuple contains isomer species names (sorted) or just a tuple (if only one species is passed).
"""
isomers = list()
done = list()
if species is not None:
speclist = [species]
else:
speclist = self.species
for spec in speclist:
stoich = parse_stoichiometry(spec)
isomergroup = None
if spec in done or spec.endswith("#"):
continue
for spec2 in self.species:
# FIXME: ice species are not isomers
if spec2 in done or spec2.endswith("#"):
continue
stoich2 = parse_stoichiometry(spec2)
if stoich == stoich2 and spec != spec2:
if isomergroup is None:
if spec in self.species:
isomergroup = [spec, spec2]
else:
isomergroup = [
spec2
] # species is not in the network so we don't add it to the isomer group
done.append(spec2)
done.append(spec)
else:
isomergroup.append(spec2)
done.append(spec2)
# append the group sorted
if species is not None:
if isomergroup is not None:
isomers = tuple(sorted(isomergroup))
else:
isomers = tuple()
elif isomergroup is not None:
isomers.append(tuple(sorted(isomergroup)))
return isomers
[docs]
def find(
self,
species: str | list[str],
reac_type: str | None = None,
formation: bool = True,
destruction: bool = True,
) -> list[Reaction]:
"""
Find all reactions that involve any of the given species.
Parameters
----------
species :
The species name(s) to search for in the reactions.
reac_type :
If provided, only reactions of this type will be returned.
formation :
If True, include formation reactions.
destruction :
If True, include destruction reactions.
Returns
-------
list[Reaction]
A list of reactions that involve the given species. Please note those are references to the original Reactions in the network,
so if you update those they are also updated in the network. If you want to avoid that you need to copy the reactions.
"""
if isinstance(species, str):
species = [species]
found_reactions = []
for reac in self.reactions:
if formation:
if any(spec in reac.products for spec in species):
if reac_type is None or reac.type == reac_type:
found_reactions.append(reac)
if destruction:
if any(spec in reac.reactants for spec in species):
if reac_type is None or reac.type == reac_type:
found_reactions.append(reac)
return found_reactions
[docs]
def remove(self, id: int | list[int] | Reaction | list[Reaction]) -> list[Reaction]:
"""
Remove a reactions or a list of reactions from the network.
Parameters
----------
id :
The reaction(s) to delete.
Returns
-------
list[:class:`Reaction`] :
The deleted reactions.
"""
if isinstance(id, int):
id = [id]
elif isinstance(id, Reaction):
id = [id.id]
elif isinstance(id, list):
if len(id) == 0: # get an empty list
return
if len(id) > 0 and isinstance(id[0], Reaction):
id = [reac.id for reac in id]
deleted = [reac for reac in self.reactions if reac.id in id]
self.reactions = [reac for reac in self.reactions if reac.id not in id]
# reset the internal species list so that it is reloaded again in method property `species``
self._species = list()
logger.info(f"Deleted {len(deleted)} reactions from the network {self.name}.")
print(f"Deleted {len(deleted)} reactions from the network {self.name}.")
return deleted
[docs]
def compare(
self,
reactionNetwork: "ReactionNetwork",
printresults: bool = False,
eqfunc: collections.abc.Callable[[dict[str, bool]],bool] | None = None,
chgfunc: collections.abc.Callable[[dict[str, bool]],bool] | None = None,
) -> tuple[
bool, list[tuple[Reaction, Reaction, dict[str, bool]]], list[Reaction], list[Reaction]
]:
"""
Compares the network to the given reaction network.
Currently simply prints out the results to stdout.
Parameters
----------
reactionNetwork :
The reaction network to compare with.
printresults :
Print results to stdout.
eqfunc :
Pass a function that decides if two reactions are equal.
The function must take as an argument the outcome of :func:`Reaction.compare`.
chgfunc :
Pass a function that decides if two reactions are the same reaction, but some other quantities differ (coefficients, type etc).
The function must take as an argument the outcome of :func:`Reaction.compare`.
Returns
-------
tuple[bool, list[tuple[Reaction, Reaction, dict[str, bool]]], list[Reaction], list[Reaction]] :
tuple containing:
- bool: are the networks equal or not
- list of changed reactions (both, and the equal dictionary)
- list of reactions only in this network
- list of reactions only in the passed network
"""
starttime = timer()
if printresults:
print("Comparing ", self.name, " to ", reactionNetwork.name, " ...")
# Copy the reaction list as we will change it
reactionsNew = reactionNetwork.reactions.copy()
eqf = self._equalreaction
if eqfunc is not None:
eqf = eqfunc
chgf = self._changedreaction
if chgfunc is not None:
chgf = chgfunc
# Compare the two databases
changed = list()
onlyold = list()
for reac in self.reactions:
found = False
for reacNew in reactionsNew:
equal = reac.compare(reacNew)
if eqf(equal): # found the equal reaction ... done
found = True
# now we can remove it is already found
reactionsNew.remove(reacNew)
break
elif chgf(equal): # found a changed reactions ... done
changed.append((reac, reacNew, equal))
found = True
# now we can remove it is already found
reactionsNew.remove(reacNew)
break
# collect the ones only in the old
if not found:
onlyold.append(reac)
# the rest must be the onlynew ones
onlynew = reactionsNew
equal = len(changed) == 0 and len(onlyold) == 0 and len(onlynew) == 0
logger.info("time: ", "{:4.2f}".format(timer() - starttime) + " s")
if printresults:
print(" ")
print("FOUND " + str(len(changed)) + " CHANGED REACTIONS: ")
for cr in changed:
print("OLD: ", cr[0])
print("NEW: ", cr[1])
print("EQUAL: ", cr[2])
print(" ")
print(" ")
print(" ")
print(
"FOUND "
+ str(len(onlyold))
+ " REACTIONS THAT ONLY EXIST IN NETWORK "
+ self.name
+ " (OLD): "
)
for r in onlyold:
print(r)
print(" ")
print(" ")
print(
"FOUND "
+ str(len(onlynew))
+ " REACTIONS THAT ONLY EXIST IN NETWORK "
+ reactionNetwork.name
+ " (NEW): "
)
for r in onlynew:
print(r)
print(" ")
return equal, changed, onlyold, onlynew
def _duplicatereaction(self, eq):
"""
Defines what is a duplicate reaction depending on the eq dictionary.
Parameters
----------
eq : dictionary
The equal dictionary see
"""
return eq["reactants"] and eq["products"]
def _equalreaction(self, eq):
"""
Are the reactions equal?
"""
return eq["reactants"] and eq["products"] and eq["coeffs"] and eq["rate"]
def _changedreaction(self, eq):
"""
Has the reaction changed?
"""
return (eq["reactants"] and eq["products"]) and (
not eq["coeffs"] or not eq["type"] or not eq["gasphase"] or not eq["rate"]
)
[docs]
@abstractmethod
def load_reactions(self, filename: str = None):
"""
In this routine the reading of the Reaction network needs to be implemented.
Fills in the `reactions` field in this class.
Parameters
----------
filename : str
The name/path of the file to read
"""
self._species = list() # reset the internal species list so that it is reloaded again in method property `species``
self.filename = filename
print(f"Loaded {len(self.reactions):<d} reactions from {self.filename}.")
# FIXME: maybe store multitemp in class
multitemp = 0
for reac in self.reactions:
if len(reac.temps) > 1:
multitemp += 1
print(f"Found {multitemp:<d} reactions with multiple temperature fits.")
print(
f"Network includes {len(self.species):<d} species and {len(self.elements):<d} elements."
)
[docs]
@abstractmethod
def load_rates(self, filename=None):
"""
In this routine the reading of rate coefficients need to be implemented.
Parameters
----------
filename : str
The name/path of the file to read
"""
pass
[docs]
def write_reactions(
self,
filename: str,
fmt: str = "csv",
coeffsfmts=("{:8.5}", "{:8.6}", "{:12.4f}"),
header: str = None,
):
"""
Writes the network to a file. The default format is the csv format as used in |prodimo| (i.e. Reactions.in.csv).
Multiple temperature fits are only supported in the csv format.
.. warning::
Not well tested ... so be careful.
Parameters
----------
filename :
The name/path of the output file. And existing one will be overwritten.
fmt :
- `csv` ... write the network in a csv format, as used by |prodimo|.
- `oldprodimo` ... write the network in the old |prodimo| format, as used in Reactions.in. This format is not recommended any more.
header :
A string that is written as a header to the file. Only used for `fmt='csv'`. The string will be written as is, but using `header.strip()`, it is recommended to
use triple-quoted strings. Default: `None`
"""
fw = open(filename, "w+")
spfmt = "{:8s}"
# TODO: use a csv writer for this.
if fmt == "csv":
if header:
fw.write(header.strip() + "\n")
for reac in self.reactions:
# Deal with multiple temperature fits, simply duplicate the reaction, with the same id, but different coefficients etc.
for coeff, temp in zip(reac.coeffs, reac.temps):
fw.write("{:5d}".format(reac.id).strip())
fw.write(",")
fw.write(reac.type.strip())
fw.write(",")
for i in range(3):
if i < len(reac.reactants):
fw.write(reac.reactants[i].strip())
else:
fw.write("")
fw.write(",")
for i in range(4):
if i < len(reac.products):
fw.write(reac.products[i].strip())
else:
fw.write("")
fw.write(",")
# assume that it is always positive
fw.write(coeffsfmts[0].format(coeff[0]).strip())
fw.write(",")
# can also be negative, but can also be an exponential
fw.write(coeffsfmts[1].format(coeff[1]).strip())
fw.write(",")
# assume that it is always positive
fw.write(coeffsfmts[2].format(coeff[2]).strip())
fw.write(",")
fw.write(reac.method.strip())
for i in range(2):
fw.write(",")
fw.write("{:6.0f}".format(temp[i]).strip())
fw.write(",")
fw.write(reac.accuracy.strip())
fw.write(",")
# first re
# properly escape the comment string
reac.comment.replace('"', '""')
reac.comment = '"' + reac.comment + '"'
fw.write(reac.comment.strip())
fw.write("\n")
elif fmt == "oldprodimo": # this is the old one, should not be used anymore
print(
"WARN: You ares requesting the old ProDiMo format, that one should not be used anymore. Use fmt='csv' instead."
)
# abcformat="{:+8.2E}"
for reac in self.reactions:
fw.write("{:5d}".format(reac.id))
fw.write(" ")
# Three reactants
for i in range(3):
if i < len(reac.reactants):
fw.write(spfmt.format(reac.reactants[i]))
else:
fw.write(spfmt.format(" "))
for i in range(4):
if i < len(reac.products):
fw.write(spfmt.format(reac.products[i]))
else:
fw.write(spfmt.format(" "))
# fw.write(" ")
# assume that it is always positive
fw.write(coeffsfmts[0].format(reac.coeffs[0]))
fw.write(" ")
# the B coefficient
# Check for numbers that don't comply to the format
dig = str(reac._coeff2).split(".")
ndig = 2
if len(dig) > 1:
ndig = len(dig[1])
# can also be negative
if reac.coeffs[1] < 0.0:
fmt = "{:+5.2F}"
else:
fmt = "{:5.2F}"
if ndig < 3:
fw.write(fmt.format(reac.coeffs[1]))
else:
# print(reac,reac._coeff2)
fw.write(str(reac._coeff2))
fw.write(" ")
# if the B coefficient has more digits than 5 we have a problem
# simple make the C coeff less long, like one would need to do if it is written by hand
lenstr2 = len(str(reac._coeff2))
lenstr3 = 11
if lenstr2 > 5:
lenstr3 = lenstr3 - (lenstr2 - 5)
# print(lenstr3)
fmt = "{:" + str(lenstr3).strip() + ".4f}"
# assume that it is always positive
fw.write(fmt.format(reac.coeffs[2]))
fw.write(" ")
fw.write(reac.comment)
fw.write("\n")
fw.close()
else:
raise ValueError(f"Unknown format {fmt} for writing reaction network.")
print("Chemical reaction network written to ", filename)
[docs]
class ReactionNetworkPout(ReactionNetwork):
"""
Implementation of ReactionNetwork for the Reactions.out of a |prodimo| model.
"""
def __init__(self, name: str = "ReactionNetworkPout", modeldir: str = None):
"""
Parameters
----------
modeldir : str
If set the init routine tries to load the Reactions.out and the
rates log file (at the moment) from the given `modeldir` directory. But continues if it does not work.
"""
if modeldir is not None:
if name == "ReactionNetworkPout":
name = putils.guess_modelname_from_path(modeldir)
filename = os.path.join(modeldir, "Reactions.out")
self.modeldir = modeldir
else:
filename = None
self.modeldir = None
super().__init__(name, filename=filename, filename_dictionary=None)
if self.modeldir is not None:
try:
# let the routine look for a possible rate coefficents file
self.load_rates(None)
except Exception as e:
logger.warning(f"Could not load rate coefficients {e}")
print()
[docs]
def load_reactions(self, filename: str = "Reactions.out"):
"""
Reads a reaction network in the format of the |prodimo| Reactions.out.
Fills in the `reactions` field in this class.
Parameters
----------
filename :
The name/path of the file to read.
"""
fr = open(filename)
self.filename = filename
lines = fr.readlines()
# stupid workaround for old format
oldformat = False
try:
for line in lines:
if not oldformat:
fields = line.split()
try:
tidx = int(fields[-6])
except ValueError:
print("WARN: Try to read an older format ... let's hope it works!")
oldformat = True
if oldformat:
# now insert some spaces at the right places for the coefficients
line = line[:23] + " " + line[23:]
line = line[:53] + " " + line[53:]
line = line[: -43 + 8] + " " + line[-43 + 8 :]
line = line[: -43 + 17] + " " + line[-43 + 17 :]
fields = line.split()
tidx = int(fields[-6])
# new Reaction
if tidx == 1:
reac = Reaction()
# new reactions, append it and fill it
self.reactions.append(reac)
reac.id = int(fields[0])
reac.idU = int(fields[1])
reac.type = fields[2]
if oldformat:
# FIXME: there seems to be different old formats, can be that old format is detected, but
# it still has the gas phase flag. So check for it
if fields[3].strip().endswith(":"):
reac.gasphase = fields[3].strip() == "T:"
reacprod = fields[4:-6]
else:
reac.gasphase = True
reacprod = fields[3:-6]
else:
reac.gasphase = fields[3].strip() == "T:"
reacprod = fields[4:-6]
reac.coeffs.append(
np.array([float(fields[-5]), float(fields[-4]), float(fields[-3])])
)
reac.temps.append(np.array([float(fields[-2]), float(fields[-1])]))
# print(reacprod)
nextisprod = False
for entry in reacprod:
if entry == "-->":
nextisprod = True
elif entry == "+":
continue
else:
if nextisprod:
reac.products.append(entry)
else:
reac.reactants.append(entry)
else:
reac = self.reactions[-1]
reac.coeffs.append(
np.array([float(fields[-5]), float(fields[-4]), float(fields[-3])])
)
reac.temps.append(np.array([float(fields[-2]), float(fields[-1])]))
except Exception as e:
print(line)
print(e)
fr.close()
raise e
# print(reac)
fr.close()
print("Loaded ", len(self.reactions), " Reactions from ", self.filename)
def _equalreaction(self, eq):
"""
Are the reactions equal?
"""
# here we also have the type
return super()._equalreaction(eq) and eq["type"]
[docs]
def load_rates(self, filename: str | None = None):
"""
Load the real rates for a given set of physical conditions.
Can be used for comparison.
Has to be run after load_reactions.
"""
if filename is None:
filename = os.path.join(self.modeldir, "rate_coeffs_1NZZ.log")
if not os.path.isfile(filename):
filename = os.path.join(self.modeldir, "MC_rate_coefficients.txt")
if not os.path.isfile(filename):
filename = None
logger.info("Could not find a rate coefficient file in ", self.modeldir)
return # Is okay does not have to be there
print("Loading rate coefficients from ", filename)
# By coincident that works for both rate_coeffs_1NZZ.log and MC_rate_coefficients.txt
rates = np.loadtxt(filename, usecols=[1])
for i, reaction in enumerate(self.reactions):
reaction.rate = rates[i]
[docs]
class ReactionNetworkPin(ReactionNetwork):
"""
Implementation of ReactionNetwork for the Reactions.in and Reactions.in.csv formats of |prodimo|.
As the csv format is the same as for the UMIST data files in ProDiMo it can/should also be used for those.
"""
def __init__(
self,
name: str,
filename: str | None = None,
filename_dictionary: str | None = None,
modeldir: str | None = None,
):
"""
Parameters
----------
name :
a name for the network.
filename :
The name/path of the file to read.
filename_dictionary :
The name/path of the dictionary file to read.
modeldir :
If set, the init routine tries to load the `Reactions.in` from within ``modeldir``.
"""
self.modeldir = modeldir
if modeldir is not None:
filename = None
for filename in ["Reactions.in", "Reactions.in.csv"]:
fullpath = os.path.join(modeldir, filename)
if os.path.isfile(fullpath):
filename = fullpath
break
super().__init__(name=name, filename=filename, filename_dictionary=filename_dictionary)
def _equalreaction(self, eq):
"""
Are the reactions equal?
"""
return eq["reactants"] and eq["products"] and eq["coeffs"]
def _changedreaction(self, eq):
"""
Has the reaction changed?
"""
return (eq["reactants"] and eq["products"]) and (not eq["coeffs"])
[docs]
def load_reactions(
self, filename: str = "Reactions.in", fmt: str | None = None, threeReactants: bool = True
):
"""
Reads a reaction network in the format of the |prodimo| Reactions.in or
Reactions.in.csv format.
Fills ``reactions`` field in this class.
.. todo::
- Include reading of T-dependent rates for the csv format.
- more sophisticated guessing of the file format.
Parameters
----------
filename : str
The name/path of the file to read.
fmt : str
Format of the file. Currently either `in` for old Reactions.in or `csv` for
the UMIST csv format. If `None` the routine tries to guess the format (primitive at the moment).
"""
# guess the format old .in or csv format
if fmt is None:
csv = filename.strip().endswith(".csv")
with open(filename) as f:
for line in f:
if line.strip() == "":
continue
if line.strip().startswith("#"):
continue
sep = re.sub(r"^\d+", "", line)[0]
if sep in [",", ":"]:
csv = True
break
if csv:
logger.info(f"Guessing format of {filename} as csv format.")
elif fmt == "csv":
csv = True
else:
csv = False
fr = open(filename)
self.filename = filename
lines = fr.readlines()
if csv:
sep = None
for line in lines:
try:
if line.strip() == "":
continue
if line.strip().startswith("#"):
continue
# Detect the separator, it is the first non-numeric character
if sep is None:
sep = re.sub(r"^\d+", "", line)[0]
# by using maxplit, everyting after the last relevant field, will be in one field and that is the comment
fields = line.split(sep, maxsplit=16)
id = int(fields[0])
temps = np.array([float(fields[13].strip()), float(fields[14].strip())])
coeffs = np.array([float(x.replace("D", "E")) for x in fields[9:12]])
# not sure if this is general, but in that case we assume it is an additional temperature range
if len(self.reactions) > 0 and id == self.reactions[-1].id:
self.reactions[-1].temps.append(temps)
self.reactions[-1].coeffs.append(coeffs)
# self.reactions[-1].coeffs=sorted(self.reactions[-1].coeffs, key=lambda x: x[2])
else: # new reaction,fill it and append it
reac = Reaction()
reac.id = id
reac.idU = None
reac.type = fields[1]
reac.method = fields[12].strip()
reac.temps.append(temps)
reac.accuracy = fields[15].strip()
# whatever comes after 16 we take as comment, but add the separator again
reac.comment = fields[-1].strip()
# strip the leading and trailing '"' if it is there, because that is what a csv reader would do
if reac.comment.startswith('"'):
reac.comment = reac.comment[1:]
if reac.comment.endswith('"'):
reac.comment = reac.comment[:-1]
# three reactants and four products
for i, sp in enumerate(fields[2:9]):
if sp == "":
continue
if i < 3:
reac.reactants.append(sp)
else:
reac.products.append(sp)
reac.coeffs.append(coeffs)
# still need to check again for multiple t-fits.
# need to go through the already loaded reactions and check
# Don't do this, this could be double arrenhius fits reactions
# found=False
# for oldreac in self.reactions:
# if sorted(oldreac.reactants) == sorted(reac.reactants) and sorted(oldreac.products) == sorted(reac.products):
# oldreac.temps.append(temps)
# oldreac.coeffs.append(coeffs)
# oldreac.comment += "//"+reac.comment
# found=True
# break
# if not found:
# self.reactions.append(reac)
self.reactions.append(reac)
except Exception as e:
print(line)
print(e)
fr.close()
raise e
else:
for line in lines:
if line.strip() == "":
continue
# do it similar to ProDiMo
# print(line)
idnum = line[0:5]
specs = line[6:61]
ABC = line[61:89].split()
comm = line[89:]
reac = Reaction()
reac.id = int(idnum.strip())
reac.idU = None
reac.type = None
reac.temps = None
# lensp=8
# for i in range(7):
# sp=specs[i*lensp:8+i*lensp].strip()
# if i<3 and sp!="":
# reac.reactants.append(sp)
# elif i>=3 and sp!="":
# reac.products.append(sp)
# very tricky, each species should take 8 characters, but the reading routine of PRoDiMo works
# also if one is only 7 characters (don't know how), but it seems in the end the species can have
# max 7 characters.
lensp = 8
idxnextsp = 0
for i in range(7):
sp = specs[idxnextsp : (idxnextsp + lensp)].strip()
# if (sp.strip()=="1"): return
# print(i,sp,idxnextsp,specs[idxnextsp:(idxnextsp+lensp)],"+",specs[idxnextsp+lensp],"+")
if i < 3 and sp != "":
reac.reactants.append(sp)
elif i >= 3 and sp != "":
reac.products.append(sp)
if i == 6: # last species
idxnextsp += lensp - 1
else:
idxnextsp += lensp
# if (idxnextsp+lensp)<61 and specs[idxnextsp+lensp]!=" ": idxnextsp-=1
reac.coeffs = np.array([[float(x.replace("D", "E")) for x in ABC]])
reac._coeff2 = decimal.Decimal(
ABC[1].strip()
) # This is just for output and only for Reactions.in
reac.comment = comm.strip()
self.reactions.append(reac)
fr.close()
super().load_reactions(filename=filename)
[docs]
def updateAdsorptionDesorptionReactions(
self, species: list[tuple[str, float]] | dict[str, float], comment=""
) -> list[Reaction]:
"""
Adds or updates existing adsorption and desorption reactions for the given species.
Only add and update operations are supported.
For e.g. CO# it would create the following reactions:
```
1281,IC,CO,dust,,CO#,dust,,,0.00E+00,+0.00000,0.0000,L,10,41000,B,
1282,DT,CO#,dust,,CO,dust,,,0.00E+00,+0.00000,899.0000,L,10,41000,B,Martin 899/1150GH06
1283,DC,CO#,CRP,,CO,,,,0.00E+00,+0.00000,899.0000,L,10,41000,B,GH06
1284,DP,CO#,PHOTON,,CO,,,,0.00E+00,+0.00000,899.0000,L,10,41000,B,GH06
```
Parameters
----------
species : list[str,float] or dict[str,float]
A list of species names or a dictionary with species names as keys and their corresponding binding energies (in K) as values.
Returns
-------
list[Reaction]
The list of the created and added reactions.
"""
if isinstance(species, dict):
species_list = list(species.items())
else:
species_list = species
created_reactions = []
for spec, binding_energy in species_list:
if "#" in spec:
gasspec = spec.replace("#", "")
icespec = spec
else:
gasspec = spec
icespec = spec + "#"
checkreac = self.find([icespec], reac_type="IC", formation=True, destruction=False)
if len(checkreac) > 0:
print(
f"INFO: Adsorption reaction id {checkreac[0].id} for {icespec} already exists, skipping it."
)
else:
adsorption_reac = Reaction()
adsorption_reac.id = len(self.reactions) + 1
adsorption_reac.type = "IC"
adsorption_reac.reactants = [gasspec, "dust"]
adsorption_reac.products = [icespec, "dust"]
adsorption_reac.coeffs = np.array([[0.0, 0.0, 0.0]]) # Placeholder coefficients
adsorption_reac.temps = np.array([[10.0, 41000.0]]) # Placeholder temperatures
adsorption_reac.method = "L"
adsorption_reac.accuracy = "B"
adsorption_reac.comment = comment
self.reactions.append(adsorption_reac)
created_reactions.append(adsorption_reac)
for type, reactant, product in zip(
["DT", "DC", "DP"], ["dust", "CRP", "PHOTON"], ["dust", None, None]
):
checkreac = self.find([icespec], reac_type=type, formation=False, destruction=True)
if len(checkreac) > 0:
print(
f"INFO: Desorption reaction id {checkreac[0].id} for {icespec} already exists, updating binding energy. "
)
checkreac[0].coeffs = np.array(
[[0.0, 0.0, binding_energy]]
) # Update coefficients
else:
# Create desorption reaction
desorption_reac = Reaction()
desorption_reac.id = len(self.reactions) + 1
desorption_reac.type = type
desorption_reac.reactants = [icespec, reactant]
if product is None:
desorption_reac.products = [gasspec]
else:
desorption_reac.products = [gasspec, product]
desorption_reac.coeffs = np.array(
[[0.0, 0.0, binding_energy]]
) # Placeholder coefficients
desorption_reac.temps = np.array([[10.0, 41000.0]]) # Placeholder temperatures
desorption_reac.method = "L"
desorption_reac.accuracy = "B"
desorption_reac.comment = comment
self.reactions.append(desorption_reac)
created_reactions.append(desorption_reac)
return created_reactions
[docs]
def load_rates(self, filename=None):
print("Not useful for Reactions.in")
[docs]
class ReactionNetworkUMIST2022(ReactionNetwork):
"""
Implementation of ReactionNetwork for the UMIST 2022 database format.
"""
def _equalreaction(self, eq):
"""
Are the reactions equal?
"""
return eq["reactants"] and eq["products"] and eq["coeffs"]
def _changedreaction(self, eq):
"""
Has the reaction changed?
"""
return (eq["reactants"] and eq["products"]) and (not eq["coeffs"])
[docs]
def load_reactions(self, filename: str = "rate22_final.rates.csv"):
"""
Reads a UMIST2022 reaction network in the format of the UMIST 2022 csv format
Fills ``reactions`` field in this class.
Parameters
----------
filename : str
The name/path of the file to read.
"""
self.filename = filename
multitemp = 0
with open(filename) as csvfile:
reader = csv.reader(csvfile, delimiter=":", quotechar='"', strict=False)
for fields in reader:
if "".join(fields) == "":
print("empty line")
continue
# allow for a comment, must be first non-whitespace character
if fields[0].strip().startswith("#"):
print("comment Line")
continue
ilast = 17 # ist the index of the last field, for 1 temp fits
reac = Reaction()
reac.id = int(fields[0])
reac.idU = int(fields[0])
reac.type = fields[1]
reac.method = fields[14].strip()
reac.accuracy = fields[15].strip()
reac.comment = fields[ilast - 1].strip() + " " + fields[ilast].strip()
# two reactants and 4 products
for i, sp in enumerate(fields[2:8]):
if sp.strip() == "":
continue
if i < 2:
reac.reactants.append(sp)
else:
reac.products.append(sp)
itemp = int(fields[8]) # number of temperature fits
reac.coeffs.append(np.array([float(x.replace("D", "E")) for x in fields[9:12]]))
reac.temps.append(np.array([float(fields[12].strip()), float(fields[13].strip())]))
# more than one temperature fit
for it in range(1, itemp):
istart = 18 + 9 * (it - 1) # index of the start of the next fit
reac.coeffs.append(
np.array([float(x.replace("D", "E")) for x in fields[istart : istart + 3]])
)
reac.temps.append(
np.array(
[float(fields[istart + 3].strip()), float(fields[istart + 4].strip())]
)
)
reac.comment += (
" TEMP: " + fields[istart + 8].strip() + fields[istart + 9].strip()
)
if itemp > 1:
multitemp += 1
self.reactions.append(reac)
# check if we got all reactions, in UMIST 2022
ids = np.array([reac.id for reac in self.reactions])
ids = ids[1:] - ids[0:-1]
if np.max(ids) > 1:
print(
"It seems reaction: "
+ str(np.argmax(ids) + 2)
+ " is missing, check the file format and the reading routine!"
)
super(ReactionNetworkUMIST2022, self).load_reactions(filename=filename)
[docs]
def load_rates(self, filename=None):
print("Not useful for UMIST 2022")
[docs]
class ReactionNetworkKIDA(ReactionNetwork):
"""
Implementation of :class:`~prodimopy.chemistry.network.ReactionNetwork` for the KIDA database format.
"""
def __init__(self, name, filename: str | None = None, filename_dictionary: str | None = None):
"""
Parameters
----------
name :
a name for the network.
filename :
The name/path of the file to read.
filename_dictionary :
The name/path of the dictionary file to read.
"""
self.field_dict = {
"R1": "",
"R2": "",
"R3": "",
"P1": "",
"P2": "",
"P3": "",
"P4": "",
"P5": "",
"alpha": 0.0,
"beta": 0.0,
"gamma": 0.0,
"uncertainty": 0.0,
"uncertainty_temp": 0.0,
"uncertainty_type": "",
"itype": 0, # Reaction type
"Tmin": 0,
"Tmax": 0,
"formula": 0,
"id": 0,
"numbersof_abg": 0,
"recommendation": 0,
"comment": "", # not part of kida but usefull if one changes/adds reactions
}
self._dict_reactions = list()
super(ReactionNetworkKIDA, self).__init__(
name, filename=filename, filename_dictionary=filename_dictionary
)
def _equalreaction(self, eq):
"""
Are the reactions equal?
"""
return eq["reactants"] and eq["products"] and eq["coeffs"]
def _changedreaction(self, eq):
"""
Has the reaction changed?
"""
return (eq["reactants"] and eq["products"]) and (not eq["coeffs"])
def _dict_from_reaction(self, reac: Reaction) -> dict:
"""
Converts a Reaction object to a dictionary with the same structure as used in the KIDA database.
Parameters
----------
reac : Reaction
The Reaction object to convert.
Returns
-------
dict
A dictionary representation of the reaction.
"""
reac_dict = dict(self.field_dict) # Create a new dictionary for each reaction
reac_dict["R1"] = reac.reactants[0] if len(reac.reactants) > 0 else ""
reac_dict["R2"] = reac.reactants[1] if len(reac.reactants) > 1 else ""
reac_dict["R3"] = reac.reactants[2] if len(reac.reactants) > 2 else ""
reac_dict["P1"] = reac.products[0] if len(reac.products) > 0 else ""
reac_dict["P2"] = reac.products[1] if len(reac.products) > 1 else ""
reac_dict["P3"] = reac.products[2] if len(reac.products) > 2 else ""
reac_dict["P4"] = reac.products[3] if len(reac.products) > 3 else ""
reac_dict["P5"] = reac.products[4] if len(reac.products) > 4 else ""
reac_dict["alpha"] = reac.coeffs[0][0] if len(reac.coeffs) > 0 else 0.0
reac_dict["beta"] = reac.coeffs[0][1] if len(reac.coeffs) > 0 else 0.0
reac_dict["gamma"] = reac.coeffs[0][2] if len(reac.coeffs) > 0 else 0.0
reac_dict["Tmin"] = int(reac.temps[0][0]) if len(reac.temps) > 0 else 0
reac_dict["Tmax"] = int(reac.temps[0][1]) if len(reac.temps) > 0 else 0
reac_dict["itype"] = int(reac.type) if reac.type is not None else -1
reac_dict["id"] = int(reac.id) if reac.id is not None else -1
reac_dict["recommendation"] = reac.method
reac_dict["comment"] = reac.comment
# Additional fields can be filled as needed
return reac_dict
[docs]
def load_reactions(self, filename: str = "KIDA_uva_2024.dat"):
"""
Reads a KIDA reaction network in the format of the KIDA database.
Fills `reactions` list in this object.
Parameters
----------
filename :
The name/path of the file to read.
"""
# Format: 3(a11),1x,5(a11),1x,3(e10.3,1x),e8.2,1x,e8.2,1x,a4,1x,i2,1x,i6,1x,i6,1x,i2,1x,i5,i2,i3
pattern = r"^(.{11})(.{11})(.{11}) (.{11})(.{11})(.{11})(.{11})(.{11}) (.{10}) (.{10}) (.{10}) (.{8}) (.{8}) (.{4}) (.{2}) (.{6}) (.{6}) (.{2}) (.{5})(.{2})(.{3})(.*)$"
with open(filename) as csvfile:
for line in csvfile:
if line.startswith("#") or line.startswith("!") or line.strip() == "":
continue
match = re.match(pattern, line)
if match is not None:
if len(match.groups()) != 22:
raise ValueError("Unexpected number of fields in line: " + line)
reac_dict = dict(self.field_dict) # Create a new dictionary for each reaction
for field, key in zip(match.groups(), self.field_dict.keys()):
dtype = type(reac_dict[key])
(value,) = map(dtype, [field.strip()])
reac_dict[key] = value
# if there is only one product, add Photon as second product, to be compatible with ProDiMo
if reac_dict["P2"].strip() == "":
print(
"Fix single product reaction: ",
reac_dict["id"],
reac_dict["R1"],
reac_dict["R2"],
reac_dict["P1"],
)
reac_dict["P2"] = "Photon"
self._dict_reactions.append(reac_dict)
# check for multiple t fits
foundmultitemp = False
# check backwards a usually the previous reaction is what we want
for reac in self.reactions[::-1]:
if reac_dict["id"] == reac.id:
reac.coeffs.append(
np.array(
[reac_dict["alpha"], reac_dict["beta"], reac_dict["gamma"]]
)
)
reac.temps.append(np.array([reac_dict["Tmin"], reac_dict["Tmax"]]))
# reac.method=reac_dict["recommentation"]
foundmultitemp = True
break
if not foundmultitemp:
reac = Reaction()
reac.id = reac_dict["id"]
reac.idU = reac_dict["id"]
for key in ["R1", "R2", "R3"]:
if len(reac_dict[key].strip()) > 0:
reac.reactants.append(reac_dict[key])
for key in ["P1", "P2", "P3", "P4", "P5"]:
if len(reac_dict[key].strip()) > 0:
reac.products.append(reac_dict[key])
reac.coeffs.append(
np.array([reac_dict["alpha"], reac_dict["beta"], reac_dict["gamma"]])
)
reac.temps.append(np.array([reac_dict["Tmin"], reac_dict["Tmax"]]))
reac.method = reac_dict["recommendation"]
reac.type = str(reac_dict["itype"])
self.reactions.append(reac)
else: # no match
print("Could not read reaction:")
print(line)
super(ReactionNetworkKIDA, self).load_reactions(filename=filename)
[docs]
def write_reactions(self, filename="Reactions.in.new", fmt="KIDA"):
"""
Writes a KIDA reaction network in the format of the KIDA database.
Parameters
----------
filename : str
The name/path of the file to write.
fmt : str
Format of the file. Currently only `KIDA` is supported.
"""
if fmt != "KIDA":
super().write_reactions(filename=filename, fmt=fmt)
else:
# we use here the dictionary, but apply skip reactions that have been deleted
# we need to sort the dictionary list, to make sure multipel T-fits are next to each other
self._dict_reactions.sort(key=lambda x: x["id"])
with open(filename, "w") as fw:
fw.write(
"! Reactants Products A B C xxxxxxxxxxxxxxxxxxxxx ITYPE Tmin Tmax formula ID xxxxx\n"
)
for reac_dict in self._dict_reactions:
found = False
for reac in self.reactions:
if reac_dict["id"] == reac.id:
found = True
break
if not found:
continue
line = ""
for r in ["R1", "R2", "R3"]:
line += f"{reac_dict[r]:11s}"
line += " "
for p in ["P1", "P2", "P3", "P4", "P5"]:
line += f"{reac_dict[p]:11s}"
line += " "
for c in ["alpha", "beta", "gamma"]:
line += f"{reac_dict[c]:10.3e} "
line += f"{reac_dict['uncertainty']:8.2e} {reac_dict['uncertainty_temp']:8.2e} {reac_dict['uncertainty_type']:4s} "
line += f"{reac_dict['itype']:2d} {reac_dict['Tmin']:6d} {reac_dict['Tmax']:6d} {reac_dict['formula']:2d} {reac_dict['id']:5d}{reac_dict['numbersof_abg']:2d}{reac_dict['recommendation']:3d}"
if reac_dict["comment"].strip() != "":
line += f" {reac_dict['comment']}"
fw.write(line + "\n")
[docs]
def load_rates(self, filename=None):
print("Not useful for KIDA")
# Some utility functions
[docs]
def fit_double_arrhenius_ranges(reactions: list[Reaction], min_points: int = 3) -> Reaction | None:
"""
Fit two single Arrhenius laws with non-overlapping, contiguous temperature
ranges to the total (summed) rate of a reaction with a double/triple ... Arrhenius fit.
The split point between the two ranges is chosen to minimise the total
log-space least-squares residual across both segments.
Returns a new Reaction with the coefficents for the two temperature ranges, The returned Reaction is a new
object identical to the input (first Reaction in `reactions` list) except that its `coeffs` and `temps` are replaced
by the two fitted ranges (one entry each), and its comment is extended with "Newly fitted with 2 temperature fit".
Routine was partly written by Claude.
Parameters
----------
reactions :
A list of Reaction objects that have the same reactants and products, but multie Arrhenius fits for the
same temperature range (only one per Reaction). It is assumend that the actual total rate of the Reaction
is the sum of those Rreactions rates. And this is what will be fitted.
min_points :
Minimum number of positive-rate points required in each segment to perform the fit. If the total number of positive-rate points is less than `2 * min_points`, the function returns `None` and does not perform the fit.
Returns
-------
A new Reaction object with the fitted coefficients and temperature ranges, or `None` if the fit cannot be performed (too few positive-rate points).
"""
rate_total = None
try:
# check if temperature ranges are equal, otherwise this does not work
t, r = reactions[0].arrhenius_rate()
for reac in reactions:
if len(reac.temps) > 1:
raise ValueError(
"Input reactions must have single Arrhenius fits (one temperature range each)."
)
t_i, r_i = reac.arrhenius_rate()
if (len(t_i[0]) != len(t[0])) or not np.allclose(t_i[0], t[0]):
raise ValueError("Temperature ranges of input reactions must be the same.")
for reac in reactions:
t, r = reac.arrhenius_rate()
if rate_total is None:
rate_total = r[0]
else:
rate_total += r[0]
# per definition we assume here the the t ranges are the same, which it should be otherwise no double arrenhius
t_all = t[0]
pos = rate_total > 0
if pos.sum() < 2 * min_points:
return None
T_pos = t_all[pos]
r_pos = rate_total[pos]
n = len(T_pos)
def _fit_segment(T, r):
"""Fit Arrhenius to a segment. Returns (coeffs, residual) or (None, inf).
coeffs = [log(alpha), beta, gamma] in the log-space design matrix."""
if len(T) < 3:
return None, np.inf
y = np.log(r)
A = np.column_stack([np.ones(len(T)), np.log(T / 300.0), -1.0 / T])
c, res, _, _ = np.linalg.lstsq(A, y, rcond=None)
res_val = float(res[0]) if len(res) > 0 else float(np.sum((y - A @ c) ** 2))
return c, res_val
# Search for the split index that minimises the combined residual
best_res = np.inf
best_i = None
best_c1 = best_c2 = None
for i in range(min_points, n - min_points + 1):
c1, r1 = _fit_segment(T_pos[:i], r_pos[:i])
c2, r2 = _fit_segment(T_pos[i:], r_pos[i:])
if c1 is None or c2 is None:
continue
if r1 + r2 < best_res:
best_res = r1 + r2
best_i = i
best_c1, best_c2 = c1, c2
if best_i is None:
return None
# Shared boundary: geometric mean of the two adjacent grid points so
# that the two ranges are contiguous (t1_hi == t2_lo) with no gap.
T_split = np.sqrt(T_pos[best_i - 1] * T_pos[best_i])
# Build and return a new Reaction object with the fitted coeffs/temps
new_reac = copy.copy(reactions[0]) # copy the first reaction as a template
new_reac.coeffs = [
np.array([float(np.exp(best_c1[0])), float(best_c1[1]), float(best_c1[2])]),
np.array([float(np.exp(best_c2[0])), float(best_c2[1]), float(best_c2[2])]),
]
new_reac.temps = [
np.array([T_pos[0], T_split]),
np.array([T_split, T_pos[-1]]),
]
new_reac.comment = (reactions[0].comment or "") + " Newly fitted with 2 temperature fit"
except Exception as e:
print(f"Fit did not work: {e}")
return None
return new_reac
# FIXME: these are not general, but keep it here for the moment
[docs]
def parse_stoichiometry(species: str) -> dict[str, int]:
"""
Parse the species names and calculates the soichiometry and charge of the species.
Several special cases are handled:
- ice suffix '#' is stripped before parsing.
- `exc` (e.g. H2exc) is stripped before parsing
- l-, c-, t- isomer prefixes are stripped before parsing
- `Z` and `PAH` are treated as elements
Emits a UserWarning for any characters that cannot be matched to a known element.
Parameters
----------
species :
The species name to parse.
Returns
-------
dict[str, int]
A dictionary containing the stoichiometry and charge of the species.
Examples
--------
>>> parse_stoichiometry('CH3ON')
({'H': 3, 'C': 1, 'N': 1, 'O': 1,'charge': 0})
>>> parse_stoichiometry('H3O+')
({'H': 3, 'O': 1, 'charge': 1})
>>> parse_stoichiometry('H-')
({'H': 1, 'charge': -1})
>>> parse_stoichiometry('H2O#')
({'H': 2, 'O': 1, 'charge': 0})
"""
# Sort symbols longest-first so multi-char symbols (Fe, Mg, Na, Cl, ...)
# are matched before single-char ones (C, N, O, ...) in the regex.
_symbols_sorted = sorted(list(elemdata._by_symbol.keys()), key=len, reverse=True)
_stoich_pattern = re.compile("(" + "|".join(re.escape(s) for s in _symbols_sorted) + r")(\d*)")
isomerprefix = ["l", "c", "t"]
# Determine charge from +/- characters before stripping them
# to count the charges, we need to remove the l- etc.
formula = species
for prefix in isomerprefix:
formula = formula.replace(prefix + "-", "")
# calculate the charge
charge = formula.count("+") - formula.count("-")
# str
for iprefix in isomerprefix:
if formula.startswith(iprefix + "-"):
formula = formula[len(iprefix) + 1 :]
formula = formula.replace("exc", "")
# Strip ice suffix and charge markers
clean = re.sub(r"[#+\-]", "", formula)
stoich = {}
unknown = []
pos = 0
while pos < len(clean):
m = _stoich_pattern.match(clean, pos)
if m:
sym = m.group(1)
count = int(m.group(2)) if m.group(2) else 1
stoich[sym] = stoich.get(sym, 0) + count
pos = m.end()
else:
unknown.append(clean[pos])
pos += 1
if unknown:
warnings.warn(
f"parse_stoichiometry: unrecognised character(s) {''.join(unknown)!r} "
f"in formula {species!r}",
UserWarning,
stacklevel=2,
)
# Sort by atomic mass
stoich = dict(sorted(stoich.items(), key=lambda kv: elemdata[kv[0]].mass))
# now add the charge to the dictionary as a special key
stoich["charge"] = charge
return stoich
[docs]
def readSpeciesin(filename="Species.in"):
"""
Reads the list of species and their abundances from a ProDiMo Species.in file.
Parameters
----------
filename : str, optional
The input file name, by default "Species.in"
Returns
-------
list[tuple[str, float]]
A list of species names and their abundances.
"""
species = []
with open(filename, "r") as f:
lines = f.readlines()
nspec = int(lines[0].strip())
for line in lines[1 : nspec + 1]:
sp = line.split()
species.append((sp[0].strip(), float(sp[1].strip())))
return species
[docs]
def updateSpeciesin(filename: str, newSpecies: list[tuple[str, float]] | dict[str, float]):
"""
Updates the list of species and their abundances in a ProDiMo Species.in file.
Supports only add an overwrite operations.
Parameters
----------
filename : str
The input file name.
newSpecies : list[tuple[str, float]] | dict[str, float]
A list or dictionary of new species names and their abundances.
"""
# Read existing species
existingSpecies = readSpeciesin(filename)
# Update or add new species
if isinstance(newSpecies, dict):
newSpecies = list(newSpecies.items())
for spec, abundance in newSpecies:
found = False
for i, (existingSpec, existingAbundance) in enumerate(existingSpecies):
if existingSpec == spec:
existingSpecies[i] = (spec, abundance)
found = True
print(f"Updating already existing species {spec}.")
break
if not found:
existingSpecies.append((spec, abundance))
print(f"Adding new species {spec}.")
with open(filename, "w") as f:
f.write(str(len(existingSpecies)) + "\n")
for spec, abundance in existingSpecies:
f.write(f"{spec:<14}{abundance:15.7} \n")
[docs]
def writesSpeciesin(
network: ReactionNetwork, filename="Species.in", excludeElems=[], excludeSpecies=[]
):
"""
Writes the list of species in the chemical network to a file in the ProDiMo format.
All species have zero (very low) initial concentrations.
Parameters
----------
network :
A Chemical network
filename :
The output file name, by default "Species.in"
"""
# The initial abundances are from McElroy+ 2012.
initial_abund = {"H2": 0.99995, "H": 0.00005, "O": 1.0, "N": 1.0, "He": 1.0}
elems = network.elements
# first exclude all elements from the elems
elemsToWrite = []
for el in elems:
if el in excludeElems:
continue
elemsToWrite.append(el)
elemsToWritePlus = [el + "+" for el in elemsToWrite if el not in initial_abund]
elemsToWrite = sorted(elemsToWrite, reverse=True)
# first filter out the species that contain the elements to exclude
# FIXME: things like Fe and F might not be treated properly
speciesToWrite = []
for sp in network.species:
sp_copy = sp
sp_copy = sp_copy.replace("+", "")
sp_copy = sp_copy.replace("-", "")
sp_copy = sp_copy.replace("++", "")
for el in elemsToWrite:
if el in sp_copy:
sp_copy = sp_copy.replace(el, "")
if (sp_copy.strip() == "" or sp_copy.isdigit()) and sp not in excludeSpecies:
speciesToWrite.append(sp)
f = open(filename, "w")
f.write(str(len(speciesToWrite)) + "\n")
for sp in speciesToWrite:
if sp in initial_abund:
f.write(f"{sp:<14}{initial_abund[sp]:7.6f} \n")
elif sp in elemsToWritePlus:
f.write(f"{sp:<14}{1.0:2.1f} \n")
else:
f.write(f"{sp:<14}{1.0e-50:5.1E}\n")
f.close()
[docs]
def readElementsin(filename="Elements.in"):
"""
Reads the list of elements and their abundances and masses from a ProDiMo Elements.in file.
Parameters
----------
filename : str, optional
The input file name, by default "Elements.in"
Returns
-------
list[tuple[str, float, float]]
A list of element names and their abundances and masses.
"""
elements = []
with open(filename, "r") as f:
lines = f.readlines()
nelem = int(lines[0].strip())
for line in lines[1 : nelem + 1]:
el = line.split()
elements.append((el[0].strip(), float(el[1].strip()), float(el[2].strip())))
return elements
[docs]
def updateElementsin(filename: str, newElements: list[tuple[str, float, float]]):
"""
Updates the list of elements in a ProDiMo Elements.in file with new elements and their abundances and masses.
Parameters
----------
filename : str
The input file name, by default "Elements.in"
newElements : list[tuple[str, float, float]], optional
A list of new elements to add to the file, by default []
"""
elements = readElementsin(filename)
elements_dict = {el[0]: (el[1], el[2]) for el in elements}
for new_el in newElements:
if new_el[0] in elements_dict:
print(f"Updating already existing element {new_el[0]}.")
else:
print(f"Adding new element {new_el[0]}.")
elements_dict[new_el[0]] = (new_el[1], new_el[2])
# updated_elements = sorted(elements_dict.items(), key=lambda x: x[0])
# only replace the part in the file that is for the elements, to keep any comment.
lines = open(filename, "r").readlines()
linesout = []
nelmentsold = int(lines[0].strip())
linesout.append(str(len(elements_dict)) + "\n")
for el, (abundance, mass) in elements_dict.items():
linesout.append(f"{el:<3}{abundance:9.5f}{mass:8.3f} \n")
for line in lines[nelmentsold + 1 :]:
linesout.append(line)
with open(filename, "w") as f:
f.writelines(linesout)
[docs]
def writesElementsin(network: ReactionNetwork, filename="Elements.in", excludeElems=[]):
"""
Writes the list of elements in the chemical network to a file in the ProDiMo format.
Parameters
----------
network :
A Chemical network
filename : str, optional
The output file name, by default "Elements.in"
"""
definit = {
"H": 1.0,
"He": 0.9e-01,
"C": 1.4e-04,
"N": 7.5e-05,
"O": 3.2e-04,
"F": 2.0e-08,
"Na": 2.0e-09,
"Mg": 7.0e-09,
"Si": 8.0e-09,
"P": 3.0e-09,
"S": 8e-08,
"Fe": 3.0e-09,
"Cl": 4.0e-09, # up to here from McElroy+ 2013 Table 3
"Al": 5.0e-08, # Millar+ 2024 Table 10
"Ar": 1.0e-08,
"Ca": 5.0e-10,
"Ti": 5.0e-10,
}
elemsToWrite = []
for el in network.elements:
if el in excludeElems:
continue
elemsToWrite.append(el)
elemmass = {el: elemdata[el].mass for el in elemsToWrite}
elemsToWrite = sorted(elemsToWrite, key=lambda x: elemmass[x])
f = open(filename, "w")
f.write(str(len(elemsToWrite)) + "\n")
for el in elemsToWrite:
f.write(f"{el:<3}{np.log10(definit[el]) + 12:9.5f}{elemdata[el].mass:8.3f} \n")
f.close()
[docs]
def read_AdsorptionEnergies(filename="AdsorptionEnergies.in") -> dict[str, float | None]:
"""
Reads the adsorption energies data from a file and returns a dictionary mapping species to their adsorption energy values.
If the value given for a species in the file is not a valid float, the value will be set to `None`, the species is still
included in the dictionary.
Parameters
----------
filename : str, optional
The input file name, by default "AdsorptionEnergies.in"
Returns
-------
dict[str, float | None]
A dictionary mapping species names to their adsorption energy values (in K). If a value is not a valid float, it will be set to `None`.
"""
adsorption_energies = {}
with open(filename, "r") as f:
for line in f:
if line.strip().startswith("#"):
continue # Skip comment lines
fields = line.strip().split()
adsorption_energies[fields[0].strip()] = float(fields[1])
return adsorption_energies
[docs]
def read_species_dictionary(filename: str) -> dict[str, str]:
"""
Reads dictionary data from a file and returns a dictionary, mapping species to their mapped names.
The format has to be a two column file, with the first column being the species name and the second column being the mapped name, separated by whitespace characters.
After this two columns, separated by whitespace, any comment is also allowed, but will not be read.
The first line of the file will always be ignored, as it is assumed to be a header line.
A `#` as the first character of a line indicates a comment line, which will also be ignored.
"""
logger.info(f"Reading species dictionary from file: {filename}")
dictionary = {}
with open(filename, "r") as f:
for i, line in enumerate(f):
if i == 0 or line.strip().startswith("#"):
continue # Skip the header line and comment lines
name, mapped_name = line.split()[0:2]
dictionary[name.strip()] = mapped_name.strip()
return dictionary